Another fix for *_RADIATIVE_BOUNDARY_3D: dx was applied on all the faces.

git-svn-id: http://svn.cactuscode.org/arrangements/CactusBase/Boundary/trunk@246 6a38eb6e-646e-4a02-a296-d141613ad6c4

1 files changed, 78 insertions, 218 deletions

diff --git a/src/RadiationBoundary.c b/src/RadiationBoundary.c
index 7ed9d15..9033d68 100644
--- a/src/RadiationBoundary.c
+++ b/src/RadiationBoundary.c
@@ -1156,50 +1156,18 @@ void CCTK_FCALL CCTK_FNAME (BndRadiativeVN)
    @vtype      int
    @vio        in
    @endvar
-   @var        xyz
-   @vdesc      coordinate arrays for x, y, and z
-   @vtype      CCTK_REAL [ MAXDIM ]
-   @vio        in
-   @endvar
-   @var        radius
-   @vdesc      radius coordinate array
-   @vtype      CCTK_REAL []
-   @vio        in
-   @endvar
-   @var        offset
-   @vdesc      offset to the next element in this direction
+   @var        dim
+   @vdesc      dimension to apply BC
    @vtype      int
    @vio        in
    @endvar
-   @var        rho
-   @vdesc      wave speed multiplied by courant factor
-   @vtype      CCTK_REAL
-   @vio        in
-   @endvar
-   @var        var_to
-   @vdesc      target variable array
-   @vtype      cctk_type []
-   @vio        in
-   @endvar
-   @var        var_from
-   @vdesc      source variable array
-   @vtype      cctk_type []
-   @vio        in
-   @endvar
    @var        cctk_type
    @vdesc      CCTK datatypes of the source and target variable
    @vtype      <cctk_type>
    @vio        in
    @endvar
 @@*/
-#define LOWER_RADIATIVE_BOUNDARY_3D(istart, jstart, kstart,                   \
-                                    xyz,                                      \
-                                    radius,                                   \
-                                    offset,                                   \
-                                    rho,                                      \
-                                    var_to,                                   \
-                                    var_from,                                 \
-                                    cctk_type)                                \
+#define LOWER_RADIATIVE_BOUNDARY_3D(istart, jstart, kstart, dim, cctk_type)   \
 {                                                                             \
   int _i, _j, _k;                                                             \
                                                                               \
@@ -1208,51 +1176,50 @@ void CCTK_FCALL CCTK_FNAME (BndRadiativeVN)
   {                                                                           \
     for (_j = jstart-1; _j >= 0; _j--)                                        \
     {                                                                         \
-      int _0 = 0*offset,                                                      \
-          _1 = 1*offset,                                                      \
-          _2 = 2*offset;                                                      \
-      int _index = CCTK_GFINDEX3D (GH, istart-1, _j, _k);                     \
-      CCTK_REAL *_radius = radius   + _index,                                 \
-                *_xyz    = xyz      + _index;                                 \
-      cctk_type *_to     = (cctk_type *) (var_to)   + _index,                 \
-                *_from   = (cctk_type *) (var_from) + _index;                 \
+      int _0 = 0*offset[dim],                                                 \
+          _1 = 1*offset[dim],                                                 \
+          _2 = 2*offset[dim];                                                 \
+      int _idx = CCTK_GFINDEX3D (GH, istart-1, _j, _k);                       \
+      CCTK_REAL *_r    = coords[MAXDIM]   + _idx,                             \
+                *_xyz  = coords[dim]      + _idx;                             \
+      cctk_type *_to   = (cctk_type *) GH->data[var_to][timelvl_to]   + _idx, \
+                *_from = (cctk_type *) GH->data[var_from][timelvl_from] + _idx;\
                                                                               \
                                                                               \
       for (_i = istart-1; _i >= 0; _i--)                                      \
       {                                                                       \
-        CCTK_REAL _radius0_inv = 1/_radius[_0],                               \
-                  _radius1_inv = 1/_radius[_1];                               \
+        CCTK_REAL _r0_inv = 1/_r[_0], _r1_inv = 1/_r[_1];                     \
                                                                               \
         if (radpower >= 0)                                                    \
         {                                                                     \
           CCTK_REAL H;                                                        \
                                                                               \
-          H  = 0.25 * radpower * dxyz[0]                                      \
-               * (_xyz[_0]*SQR (_radius0_inv) + _xyz[_1]*SQR (_radius1_inv)); \
+          H  = 0.25 * radpower * dxyz[dim]                                    \
+               * (_xyz[_0]*SQR (_r0_inv) + _xyz[_1]*SQR (_r1_inv));           \
           H  = (1 + H) / (1 - H);                                             \
           H *= dtv * (0.25*(_to[_1] + _to[_2] + _from[_1] + _from[_2]) - var0)\
-               + 0.5*(  _radius[_1] * (_to[_1] - _from[_1])                   \
-                      + _radius[_2] * (_to[_2] - _from[_2]))                  \
-               + 0.25*(_to[_2] - _to[_1] + _from[_2] - _from[_1]) * rho       \
-                     *(  SQR (_radius[_1]) / _xyz[_1]                         \
-                       + SQR (_radius[_2]) / _xyz[_2]);                       \
+               + 0.5*(  _r[_1] * (_to[_1] - _from[_1])                        \
+                      + _r[_2] * (_to[_2] - _from[_2]))                       \
+               + 0.25*(_to[_2] - _to[_1] + _from[_2] - _from[_1]) * rho[dim]  \
+                     *(  SQR (_r[_1]) / _xyz[_1]                              \
+                       + SQR (_r[_2]) / _xyz[_2]);                            \
           dtvvar0H = dtvvar0 + H;                                             \
         }                                                                     \
                                                                               \
         _to[_0] = (cctk_type) (                                               \
-                   (dtvvar0H * (  _xyz[_0] * SQR (_radius0_inv)               \
-                                + _xyz[_1] * SQR (_radius1_inv))              \
-                    - _to[_1]   * ( rho + _xyz[_1]*_radius1_inv               \
-                                                     *(1 + dtvh*_radius1_inv))\
-                    + _from[_0] * ( rho + _xyz[_0]*_radius0_inv               \
-                                                     *(1 - dtvh*_radius0_inv))\
-                    - _from[_1] * ( rho - _xyz[_1]*_radius1_inv               \
-                                                     *(1 - dtvh*_radius1_inv))\
+                   (dtvvar0H * (  _xyz[_0] * SQR (_r0_inv)                    \
+                                + _xyz[_1] * SQR (_r1_inv))                   \
+                    - _to[_1]   * ( rho[dim] + _xyz[_1]*_r1_inv               \
+                                                     *(1 + dtvh*_r1_inv))     \
+                    + _from[_0] * ( rho[dim] + _xyz[_0]*_r0_inv               \
+                                                     *(1 - dtvh*_r0_inv))     \
+                    - _from[_1] * ( rho[dim] - _xyz[_1]*_r1_inv               \
+                                                     *(1 - dtvh*_r1_inv))     \
                    )                                                          \
-                   /              (-rho + _xyz[_0]*_radius0_inv               \
-                                                     *(1 + dtvh*_radius0_inv))\
+                   /              (-rho[dim] + _xyz[_0]*_r0_inv               \
+                                                     *(1 + dtvh*_r0_inv))     \
                   );                                                          \
-        _radius--; _xyz--; _to--; _from--;                                    \
+        _r--; _xyz--; _to--; _from--;                                         \
       }                                                                       \
     }                                                                         \
   }                                                                           \
@@ -1272,50 +1239,18 @@ void CCTK_FCALL CCTK_FNAME (BndRadiativeVN)
    @vtype      int
    @vio        in
    @endvar
-   @var        xyz
-   @vdesc      coordinate arrays for x, y, and z
-   @vtype      CCTK_REAL [ MAXDIM ]
-   @vio        in
-   @endvar
-   @var        radius
-   @vdesc      radius coordinate array
-   @vtype      CCTK_REAL []
-   @vio        in
-   @endvar
-   @var        offset
-   @vdesc      offset to the next element in this direction
+   @var        dim
+   @vdesc      dimension to apply BC
    @vtype      int
    @vio        in
    @endvar
-   @var        rho
-   @vdesc      wave speed multiplied by courant factor
-   @vtype      CCTK_REAL
-   @vio        in
-   @endvar
-   @var        var_to
-   @vdesc      target variable array
-   @vtype      cctk_type []
-   @vio        in
-   @endvar
-   @var        var_from
-   @vdesc      source variable array
-   @vtype      cctk_type []
-   @vio        in
-   @endvar
    @var        cctk_type
    @vdesc      CCTK datatypes of the source and target variable
    @vtype      <cctk_type>
    @vio        in
    @endvar
 @@*/
-#define UPPER_RADIATIVE_BOUNDARY_3D(istart, jstart, kstart,                   \
-                                    xyz,                                      \
-                                    radius,                                   \
-                                    offset,                                   \
-                                    rho,                                      \
-                                    var_to,                                   \
-                                    var_from,                                 \
-                                    cctk_type)                                \
+#define UPPER_RADIATIVE_BOUNDARY_3D(istart, jstart, kstart, dim, cctk_type)   \
 {                                                                             \
   int _i, _j, _k;                                                             \
                                                                               \
@@ -1324,51 +1259,50 @@ void CCTK_FCALL CCTK_FNAME (BndRadiativeVN)
   {                                                                           \
     for (_j = jstart; _j < GH->cctk_lsh[1]; _j++)                             \
     {                                                                         \
-      int _0 = -0*offset,                                                     \
-          _1 = -1*offset,                                                     \
-          _2 = -2*offset;                                                     \
-      int _index = CCTK_GFINDEX3D (GH, istart, _j, _k);                       \
-      CCTK_REAL *_radius = radius   + _index,                                 \
-                *_xyz    = xyz      + _index;                                 \
-      cctk_type *_to     = (cctk_type *) (var_to)   + _index,                 \
-                *_from   = (cctk_type *) (var_from) + _index;                 \
+      int _0 = -0*offset[dim],                                                \
+          _1 = -1*offset[dim],                                                \
+          _2 = -2*offset[dim];                                                \
+      int _idx = CCTK_GFINDEX3D (GH, istart, _j, _k);                         \
+      CCTK_REAL *_r    = coords[MAXDIM]   + _idx,                             \
+                *_xyz  = coords[dim]      + _idx;                             \
+      cctk_type *_to   = (cctk_type *) GH->data[var_to][timelvl_to]   + _idx, \
+                *_from = (cctk_type *) GH->data[var_from][timelvl_from] + _idx;\
                                                                               \
                                                                               \
       for (_i = istart; _i < GH->cctk_lsh[0]; _i++)                           \
       {                                                                       \
-        CCTK_REAL _radius0_inv = 1/_radius[_0],                               \
-                  _radius1_inv = 1/_radius[_1];                               \
+        CCTK_REAL _r0_inv = 1/_r[_0], _r1_inv = 1/_r[_1];                     \
                                                                               \
         if (radpower >= 0)                                                    \
         {                                                                     \
           CCTK_REAL H;                                                        \
                                                                               \
-          H  = 0.25 * radpower * dxyz[0]                                      \
-               * (_xyz[_0]*SQR (_radius0_inv) + _xyz[_1]*SQR (_radius1_inv)); \
+          H  = 0.25 * radpower * dxyz[dim]                                    \
+               * (_xyz[_0]*SQR (_r0_inv) + _xyz[_1]*SQR (_r1_inv));           \
           H  = (1 - H) / (1 + H);                                             \
           H *= dtv * (0.25*(_to[_1] + _to[_2] + _from[_1] + _from[_2]) - var0)\
-               + 0.5 *(  _radius[_1] * (_to[_1] - _from[_1])                  \
-                       + _radius[_2] * (_to[_2] - _from[_2]))                 \
-               + 0.25*(_to[_1] - _to[_2] + _from[_1] - _from[_2]) * rho       \
-                     * (  SQR (_radius[_1]) / _xyz[_1]                        \
-                        + SQR (_radius[_2]) / _xyz[_2]);                      \
+               + 0.5 *(  _r[_1] * (_to[_1] - _from[_1])                       \
+                       + _r[_2] * (_to[_2] - _from[_2]))                      \
+               + 0.25*(_to[_1] - _to[_2] + _from[_1] - _from[_2]) * rho[dim]  \
+                     * (  SQR (_r[_1]) / _xyz[_1]                             \
+                        + SQR (_r[_2]) / _xyz[_2]);                           \
           dtvvar0H = dtvvar0 + H;                                             \
         }                                                                     \
                                                                               \
         _to[_0] = (cctk_type) (                                               \
-                   (dtvvar0H * (  _xyz[_0] * (SQR (_radius0_inv))             \
-                                + _xyz[_1] * (SQR (_radius1_inv)))            \
-                    + _to[_1]   * ( rho - _xyz[_1]*_radius1_inv               \
-                                                     *(1 + dtvh*_radius1_inv))\
-                    + _from[_0] * (-rho + _xyz[_0]*_radius0_inv               \
-                                                     *(1 - dtvh*_radius0_inv))\
-                    + _from[_1] * ( rho + _xyz[_1]*_radius1_inv               \
-                                                     *(1 - dtvh*_radius1_inv))\
+                   (dtvvar0H * (  _xyz[_0] * (SQR (_r0_inv))                  \
+                                + _xyz[_1] * (SQR (_r1_inv)))                 \
+                    + _to[_1]   * ( rho[dim] - _xyz[_1]*_r1_inv               \
+                                                     *(1 + dtvh*_r1_inv))     \
+                    + _from[_0] * (-rho[dim] + _xyz[_0]*_r0_inv               \
+                                                     *(1 - dtvh*_r0_inv))     \
+                    + _from[_1] * ( rho[dim] + _xyz[_1]*_r1_inv               \
+                                                     *(1 - dtvh*_r1_inv))     \
                    )                                                          \
-                   /              ( rho + _xyz[_0]*_radius0_inv               \
-                                                     *(1 + dtvh*_radius0_inv))\
+                   /              ( rho[dim] + _xyz[_0]*_r0_inv               \
+                                                     *(1 + dtvh*_r0_inv))     \
                   );                                                          \
-        _radius++; _xyz++; _to++; _from++;                                    \
+        _r++; _xyz++; _to++; _from++;                                         \
       }                                                                       \
     }                                                                         \
   }                                                                           \
@@ -1391,118 +1325,62 @@ void CCTK_FCALL CCTK_FNAME (BndRadiativeVN)
    @vtype      int [ dim ]
    @vio        in
    @endvar
-   @var        doBC
-   @vdesc      flags telling whether to apply BC in a given direction or not
-   @vtype      int [ 2*dim ]
-   @vio        in
-   @endvar
    @var        stencil
    @vdesc      stencils in every direction
    @vtype      int [ 2*dim ]
    @vio        in
    @endvar
-   @var        coords
-   @vdesc      coordinate arrays for x, y, z, and the radius
-   @vtype      CCTK_REAL [ MAXDIM+1 ]
-   @vio        in
-   @endvar
-   @var        offset
-   @vdesc      offsets to the next element in each dimension
-   @vtype      int [ dim ]
-   @vio        in
-   @endvar
-   @var        rho
-   @vdesc      wave speed multiplied by courant factor
-   @vtype      CCTK_REAL []
-   @vio        in
-   @endvar
-   @var        var_to
-   @vdesc      target variable array
-   @vtype      cctk_type []
-   @vio        in
-   @endvar
-   @var        var_from
-   @vdesc      source variable array
-   @vtype      cctk_type []
-   @vio        in
-   @endvar
-   @var        dim
-   @vdesc      dimension of the source and target variable
-   @vtype      int
-   @vio        in
-   @endvar
    @var        cctk_type
    @vdesc      CCTK datatypes of the source and target variable
    @vtype      <cctk_type>
    @vio        in
    @endvar
 @@*/
-#define RADIATIVE_BOUNDARY(lsh,                                               \
-                           doBC,                                              \
-                           stencil,                                           \
-                           coords,                                            \
-                           offset,                                            \
-                           rho,                                               \
-                           var_to,                                            \
-                           var_from,                                          \
-                           dim,                                               \
-                           cctk_type)                                         \
+#define RADIATIVE_BOUNDARY(lsh, stencil, cctk_type)                           \
 {                                                                             \
   /* check the dimensionality */                                              \
-  if (dim != 3)                                                               \
+  if (gdim != 3)                                                              \
   {                                                                           \
     CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,                      \
                 "ApplyBndRadiative: variable dimension of %d not supported",  \
-                dim);                                                         \
+                gdim);                                                        \
     return (-5);                                                              \
   }                                                                           \
                                                                               \
   /* Lower x-bound */                                                         \
   if (doBC[0])                                                                \
   {                                                                           \
-    LOWER_RADIATIVE_BOUNDARY_3D (stencil[0], lsh[1], lsh[2],                  \
-                                 coords[0], coords[MAXDIM], offset[0], rho[0],\
-                                 var_to, var_from, cctk_type);                \
+    LOWER_RADIATIVE_BOUNDARY_3D (stencil[0], lsh[1], lsh[2], 0, cctk_type);   \
   }                                                                           \
                                                                               \
   /* Upper x-bound */                                                         \
   if (doBC[1])                                                                \
   {                                                                           \
-    UPPER_RADIATIVE_BOUNDARY_3D (lsh[0] - stencil[1], 0, 0,                   \
-                                 coords[0], coords[MAXDIM], offset[0], rho[0],\
-                                 var_to, var_from, cctk_type);                \
+    UPPER_RADIATIVE_BOUNDARY_3D (lsh[0] - stencil[1], 0, 0, 0, cctk_type);    \
   }                                                                           \
                                                                               \
   /* Lower y-bound */                                                         \
   if (doBC[2])                                                                \
   {                                                                           \
-    LOWER_RADIATIVE_BOUNDARY_3D (lsh[0], stencil[2], lsh[2],                  \
-                                 coords[1], coords[MAXDIM], offset[1], rho[1],\
-                                 var_to, var_from, cctk_type);                \
+    LOWER_RADIATIVE_BOUNDARY_3D (lsh[0], stencil[2], lsh[2], 1, cctk_type);   \
   }                                                                           \
                                                                               \
   /* Upper y-bound */                                                         \
   if (doBC[3])                                                                \
   {                                                                           \
-    UPPER_RADIATIVE_BOUNDARY_3D (0, lsh[1] - stencil[3], 0,                   \
-                                 coords[1], coords[MAXDIM], offset[1], rho[1],\
-                                 var_to, var_from, cctk_type);                \
+    UPPER_RADIATIVE_BOUNDARY_3D (0, lsh[1] - stencil[3], 0, 1, cctk_type);    \
   }                                                                           \
                                                                               \
   /* Lower z-bound */                                                         \
   if (doBC[4])                                                                \
   {                                                                           \
-    LOWER_RADIATIVE_BOUNDARY_3D (lsh[0], lsh[1], stencil[4],                  \
-                                 coords[2], coords[MAXDIM], offset[2], rho[2],\
-                                 var_to, var_from, cctk_type);                \
+    LOWER_RADIATIVE_BOUNDARY_3D (lsh[0], lsh[1], stencil[4], 2, cctk_type);   \
   }                                                                           \
                                                                               \
   /* Upper z-bound */                                                         \
   if (doBC[5])                                                                \
   {                                                                           \
-    UPPER_RADIATIVE_BOUNDARY_3D (0, 0, lsh[2] - stencil[5],                   \
-                                 coords[2], coords[MAXDIM], offset[2], rho[2],\
-                                 var_to, var_from, cctk_type);                \
+    UPPER_RADIATIVE_BOUNDARY_3D (0, 0, lsh[2] - stencil[5], 2, cctk_type);    \
   }                                                                           \
 }
 
@@ -1734,68 +1612,50 @@ static int ApplyBndRadiative (const cGH *GH,
     switch (CCTK_VarTypeI (var_to))
     {
       case CCTK_VARIABLE_CHAR:
-        RADIATIVE_BOUNDARY (GH->cctk_lsh, doBC, widths, coords, offset, rho,
-                            GH->data[var_to][timelvl_to],
-                            GH->data[var_from][timelvl_from], gdim, CCTK_CHAR);
+        RADIATIVE_BOUNDARY (GH->cctk_lsh, widths, CCTK_CHAR);
         break;
 
       case CCTK_VARIABLE_INT:
-        RADIATIVE_BOUNDARY (GH->cctk_lsh, doBC, widths, coords, offset, rho,
-                            GH->data[var_to][timelvl_to],
-                            GH->data[var_from][timelvl_from], gdim, CCTK_INT);
+        RADIATIVE_BOUNDARY (GH->cctk_lsh, widths, CCTK_INT);
         break;
 
       case CCTK_VARIABLE_REAL:
-        RADIATIVE_BOUNDARY (GH->cctk_lsh, doBC, widths, coords, offset, rho,
-                            GH->data[var_to][timelvl_to],
-                            GH->data[var_from][timelvl_from], gdim, CCTK_REAL);
+        RADIATIVE_BOUNDARY (GH->cctk_lsh, widths, CCTK_REAL);
         break;
 
 #ifdef CCTK_INT2
       case CCTK_VARIABLE_INT2:
-        RADIATIVE_BOUNDARY (GH->cctk_lsh, doBC, widths, coords, offset, rho,
-                            GH->data[var_to][timelvl_to],
-                            GH->data[var_from][timelvl_from], gdim, CCTK_INT2);
+        RADIATIVE_BOUNDARY (GH->cctk_lsh, widths, CCTK_INT2);
         break;
 #endif
 
 #ifdef CCTK_INT4
       case CCTK_VARIABLE_INT4:
-        RADIATIVE_BOUNDARY (GH->cctk_lsh, doBC, widths, coords, offset, rho,
-                            GH->data[var_to][timelvl_to],
-                            GH->data[var_from][timelvl_from], gdim, CCTK_INT4);
+        RADIATIVE_BOUNDARY (GH->cctk_lsh, widths, CCTK_INT4);
         break;
 #endif
 
 #ifdef CCTK_INT8
       case CCTK_VARIABLE_INT8:
-        RADIATIVE_BOUNDARY (GH->cctk_lsh, doBC, widths, coords, offset, rho,
-                            GH->data[var_to][timelvl_to],
-                            GH->data[var_from][timelvl_from], gdim, CCTK_INT8);
+        RADIATIVE_BOUNDARY (GH->cctk_lsh, widths, CCTK_INT8);
         break;
 #endif
 
 #ifdef CCTK_REAL4
       case CCTK_VARIABLE_REAL4:
-        RADIATIVE_BOUNDARY (GH->cctk_lsh, doBC, widths, coords, offset, rho,
-                            GH->data[var_to][timelvl_to],
-                            GH->data[var_from][timelvl_from], gdim, CCTK_REAL4);
+        RADIATIVE_BOUNDARY (GH->cctk_lsh, widths, CCTK_REAL4);
         break;
 #endif
 
 #ifdef CCTK_REAL8
       case CCTK_VARIABLE_REAL8:
-        RADIATIVE_BOUNDARY (GH->cctk_lsh, doBC, widths, coords, offset, rho,
-                            GH->data[var_to][timelvl_to],
-                            GH->data[var_from][timelvl_from], gdim, CCTK_REAL8);
+        RADIATIVE_BOUNDARY (GH->cctk_lsh, widths, CCTK_REAL8);
         break;
 #endif
 
 #ifdef CCTK_REAL16
       case CCTK_VARIABLE_REAL16:
-        RADIATIVE_BOUNDARY (GH->cctk_lsh, doBC, widths, coords, offset, rho,
-                            GH->data[var_to][timelvl_to],
-                            GH->data[var_from][timelvl_from], gdim,CCTK_REAL16);
+        RADIATIVE_BOUNDARY (GH->cctk_lsh, widths, CCTK_REAL16);
         break;
 #endif