Miguel's new code for applying Robin BC, converted into C.

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

1 files changed, 168 insertions, 101 deletions

diff --git a/src/RobinBoundary.c b/src/RobinBoundary.c
index fe8427f..6d78c01 100644
--- a/src/RobinBoundary.c
+++ b/src/RobinBoundary.c
@@ -1,7 +1,7 @@
  /*@@
    @file      RobinBoundary.c
    @date      July 6th 2000
-   @author    Gabrielle Allen, Gerd Lanfermann
+   @author    Miguel Alcubierre, Gabrielle Allen, Gerd Lanfermann
    @desc
               Routines for Robin boundary conditions
    @enddesc
@@ -394,85 +394,84 @@ void CCTK_FCALL CCTK_FNAME (BndRobinVN)
 /* maximum dimension we can deal with */
 #define MAXDIM  3
 
-/* macro to compute the linear index of a 3D point */
-#define INDEX_3D(lsh, i, j, k)      ((i) + (lsh)[0]*((j) + (lsh)[1]*(k)))
+/* macro to compute x*x */
+#define SQR(x)  ((x) * (x))
+
 
 /*@@
-   @routine    ROBIN_BOUNDARY_3D_TYPED
-   @date       Tue 10 April 2001
+   @routine    SET_LINEAR_INDICES
+   @date       Thu 7 June 2001
    @author     Thomas Radke
    @desc
-               Macro to apply Robin boundary conditions to a 3D variable
-               of given datatype
+               Macro to set the linear indices for the source and destination
+               element of the current grid variable
    @enddesc
 
-   @var        doBC
-   @vdesc      flag telling whether to apply boundary conditions or not
-   @vtype      int
-   @vio        in
-   @endvar
-   @var        iend, jend
-   @vdesc      upper ranges for the loopers
-   @vtype      int
-   @vio        in
-   @endvar
-   @var        ii, jj, kk
-   @vdesc      indices of the current grid point
-   @vtype      int
-   @vio        in
-   @endvar
-   @var        offset
-   @vdesc      offset to the next neighbouring grid point
+   @var        i
+   @vdesc      x index to use
    @vtype      int
    @vio        in
    @endvar
+@@*/
+#define SET_LINEAR_INDICES(i)                                                 \
+{                                                                             \
+  dst = CCTK_GFINDEX3D (GH, i, j, k);                                         \
+  src = CCTK_GFINDEX3D (GH, i+dx, j+dy, k+dz);                                \
+  distance = dist[abs(dx) + 2*abs (dy) + 4*abs (dz)];                         \
+}
+
+
+/*@@
+   @routine    ROBIN_BOUNDARY_TYPED_3D
+   @date       Thu 7 June 2001
+   @author     Thomas Radke
+   @desc
+               Macro to apply Robin boundary conditions to a 3D variable
+               of given datatype
+   @enddesc
+
    @var        cctk_type
    @vdesc      CCTK datatype of the variable
    @vtype      <cctk_type>
    @vio        in
    @endvar
 @@*/
-#define ROBIN_BOUNDARY_3D_TYPED(doBC,                                         \
-                                iend, jend,                                   \
-                                ii, jj, kk,                                   \
-                                offset,                                       \
-                                cctk_type)                                    \
+#define ROBIN_BOUNDARY_TYPED_3D(cctk_type)                                    \
 {                                                                             \
-  if (doBC)                                                                   \
-  {                                                                           \
-    for (j = 0; j < jend; j++)                                                \
-    {                                                                         \
-      for (i = 0; i < iend; i++)                                              \
-      {                                                                       \
-        int _i = CCTK_GFINDEX3D (GH, ii, jj, kk);                             \
-        cctk_type *_var = (cctk_type *) (GH->data[var][timelvl]) + _i;        \
-        const CCTK_REAL *_r = radius + _i;                                    \
-        /* the native C compiler on DEC Alphas has problems with automatic    \
-           type-casting of integers to doubles for inlined math routines */   \
-        const double dpow = (double) npow;                                    \
+  cctk_type *data;                                                            \
+  double u_src, u_dst, aux;                                                   \
                                                                               \
                                                                               \
-        _var[0] = (cctk_type) (                                               \
-                   (  1 * pow (_r[3*offset], dpow) * (_var[3*offset] - finf)  \
-                    - 3 * pow (_r[2*offset], dpow) * (_var[2*offset] - finf)  \
-                    + 3 * pow (_r[1*offset], dpow) * (_var[1*offset] - finf)  \
-                   ) /    pow (_r[0], dpow)                                   \
-                  + finf);                                                    \
-      }                                                                       \
-    }                                                                         \
+  /* avoid the else branch with the expensive sqrt() operation if possible */ \
+  if (abs (dx) + abs (dy) + abs (dz) == 1)                                    \
+  {                                                                           \
+    u_dst = dx ? x[dst] : (dy ? y[dst] : z[dst]);                             \
+    u_src = dx ? x[src] : (dy ? y[src] : z[src]);                             \
+  }                                                                           \
+  else                                                                        \
+  {                                                                           \
+    u_dst = sqrt (SQR (dx * x[dst]) + SQR (dy * y[dst]) + SQR (dz * z[dst])); \
+    u_src = sqrt (SQR (dx * x[src]) + SQR (dy * y[src]) + SQR (dz * z[src])); \
   }                                                                           \
+                                                                              \
+  aux = decay * distance * (u_src / SQR (r[src]) + u_dst / SQR (r[dst]));     \
+                                                                              \
+  data = (cctk_type *) GH->data[var][0];                                      \
+  data[dst] = (cctk_type) ((2*aux*  data[src] + finf*(1 - aux)) / (1 + aux)); \
 }
 
 
 /*@@
    @routine    ROBIN_BOUNDARY
-   @date       Sat 20 Jan 2001
+   @date       Thu 7 June 2001
    @author     Thomas Radke
    @desc
                Macro to apply Robin boundary conditions to a variable
                of a given datatype in all directions
                Currently it is limited up to 3D variables only.
    @enddesc
+   @calls      SET_LINEAR_INDICES
+               ROBIN_BOUNDARY_TYPED_3D
 
    @var        cctk_type
    @vdesc      CCTK datatype of the variable
@@ -482,6 +481,11 @@ void CCTK_FCALL CCTK_FNAME (BndRobinVN)
 @@*/
 #define ROBIN_BOUNDARY(cctk_type)                                             \
 {                                                                             \
+  int i, j, k;                                                                \
+  int dx, dy, dz;                                                             \
+  int src, dst, distance;                                                     \
+                                                                              \
+                                                                              \
   /* check the dimensionality */                                              \
   if (gdim != 3)                                                              \
   {                                                                           \
@@ -491,25 +495,65 @@ void CCTK_FCALL CCTK_FNAME (BndRobinVN)
     return (-5);                                                              \
   }                                                                           \
                                                                               \
-  /* now set the boundaries face by face */                                   \
-  /* lower x */                                                               \
-  ROBIN_BOUNDARY_3D_TYPED (doBC[0], lsh[1], lsh[2], 0, i, j,                  \
-                           offset[0], cctk_type);                             \
-  /* upper x */                                                               \
-  ROBIN_BOUNDARY_3D_TYPED (doBC[1], lsh[1], lsh[2], lsh[0]-1, i, j,           \
-                           -offset[0], cctk_type);                            \
-  /* lower y */                                                               \
-  ROBIN_BOUNDARY_3D_TYPED (doBC[2], lsh[0], lsh[2], i, 0, j,                  \
-                           offset[1], cctk_type);                             \
-  /* upper y */                                                               \
-  ROBIN_BOUNDARY_3D_TYPED (doBC[3], lsh[0], lsh[2], i, lsh[1]-1, j,           \
-                           -offset[1], cctk_type);                            \
-  /* lower z */                                                               \
-  ROBIN_BOUNDARY_3D_TYPED (doBC[4], lsh[0], lsh[1], i, j, 0,                  \
-                           offset[2], cctk_type);                             \
-  /* upper z */                                                               \
-  ROBIN_BOUNDARY_3D_TYPED (doBC[5], lsh[0], lsh[1], i, j, lsh[2]-1,           \
-                           -offset[2], cctk_type);                            \
+  /* outermost loop over all z points */                                      \
+  for (k = 0; k < GH->cctk_lsh[2]; k++)                                       \
+  {                                                                           \
+    if ((k == 0 && ! doBC[4]) || (k == GH->cctk_lsh[2]-1 && ! doBC[5]))       \
+    {                                                                         \
+      continue;                                                               \
+    }                                                                         \
+                                                                              \
+    dz = 0;                                                                   \
+    if (k == 0 && doBC[4])                                                    \
+    {                                                                         \
+      dz = +1;                                                                \
+    }                                                                         \
+    else if (k == GH->cctk_lsh[2]-1 && doBC[5])                               \
+    {                                                                         \
+      dz = -1;                                                                \
+    }                                                                         \
+                                                                              \
+    /* middle loop over all y points */                                       \
+    for (j = 0; j < GH->cctk_lsh[1]; j++)                                     \
+    {                                                                         \
+      dy = 0;                                                                 \
+      if (j == 0 && doBC[2])                                                  \
+      {                                                                       \
+        dy = +1;                                                              \
+      }                                                                       \
+      else if (j == GH->cctk_lsh[1]-1 && doBC[3])                             \
+      {                                                                       \
+        dy = -1;                                                              \
+      }                                                                       \
+                                                                              \
+      /* lower x */                                                           \
+      if (doBC[0])                                                            \
+      {                                                                       \
+        dx = +1;                                                              \
+        SET_LINEAR_INDICES (0);                                               \
+        ROBIN_BOUNDARY_TYPED_3D (cctk_type);                                  \
+      }                                                                       \
+                                                                              \
+      /* lower/upper y and/or z */                                            \
+      if (dy || dz)                                                           \
+      {                                                                       \
+        dx = 0;                                                               \
+        SET_LINEAR_INDICES (1);                                               \
+        for (i = 1; i < GH->cctk_lsh[0]-1; i++, src++, dst++)                 \
+        {                                                                     \
+          ROBIN_BOUNDARY_TYPED_3D (cctk_type);                                \
+        }                                                                     \
+      }                                                                       \
+                                                                              \
+      /* upper x */                                                           \
+      if (doBC[1])                                                            \
+      {                                                                       \
+        dx = -1;                                                              \
+        SET_LINEAR_INDICES (GH->cctk_lsh[0]-1);                               \
+        ROBIN_BOUNDARY_TYPED_3D (cctk_type);                                  \
+      }                                                                       \
+    }                                                                         \
+  }                                                                           \
 }
 
 
@@ -522,8 +566,8 @@ void CCTK_FCALL CCTK_FNAME (BndRobinVN)
                given by their indices
                This routine is called by the various BndRobinXXX wrappers.
 
-               Although it is currently limited to handle 1D, 2D, or 3D
-               variables only it can easily be extended for higher dimensions
+               Although it is currently limited to handle 3D variables only
+               it can easily be extended for higher dimensions
                by adapting the appropriate macros.
    @enddesc
 
@@ -538,7 +582,7 @@ void CCTK_FCALL CCTK_FNAME (BndRobinVN)
    @vio        in
    @endvar
    @var        finf
-   @vdesc      value of f at infimum
+   @vdesc      value of f at infinity
    @vtype      CCTK_REAL
    @vio        in
    @endvar
@@ -576,6 +620,7 @@ void CCTK_FCALL CCTK_FNAME (BndRobinVN)
                -3 if stencil width is other than 1
                -4 if variable type is not supported
                -5 if variable dimension is other than 3D
+               -6 if no coordinate information is available
    @endreturndesc
 @@*/
 static int ApplyBndRobin (const cGH *GH,
@@ -585,19 +630,17 @@ static int ApplyBndRobin (const cGH *GH,
                           int first_var,
                           int num_vars)
 {
-  int i, j;
-  int var, vtype, gindex, gdim, timelvl;
-  int doBC[2*MAXDIM], lsh[MAXDIM], offset[MAXDIM];
+  int var, vtype, dim, gdim;
+  int doBC[2*MAXDIM];
   SymmetryGHex *sGHex;
-  const CCTK_REAL *radius;
   char coord_system_name[20];
+  double decay;
+  const CCTK_REAL *x, *y, *z, *r;
+  double dist[8];
 
 
-  /* get the group index of the variables */
-  gindex = CCTK_GroupIndexFromVarI (first_var);
-
   /* get the number of dimensions and the variables' type */
-  gdim  = CCTK_GroupDimI (gindex);
+  gdim  = CCTK_GroupDimI (CCTK_GroupIndexFromVarI (first_var));
   vtype = CCTK_VarTypeI (first_var);
 
   /* make sure we can deal with this number of dimensions */
@@ -615,30 +658,55 @@ static int ApplyBndRobin (const cGH *GH,
     return (-2);
   }
 
-  for (i = 0; i < gdim; i++)
+  for (dim = 0; dim < gdim; dim++)
   {
-    if (stencil[i] != 1)
+    if (stencil[dim] != 1)
     {
-      CCTK_WARN (1, "ApplyBndRobin: Stencil width must be 1 for Robin boundary conditions");
+      CCTK_WARN (1, "ApplyBndRobin: Stencil width must be 1 "
+                    "for Robin boundary conditions");
       return (-3);
     }
-    offset[i] = i == 0 ? 1 : offset[i-1] * GH->cctk_lsh[i-1];
   }
 
-  /* initialize array for variables with less dimensions than MAXDIM
-     so that we can use the INDEX_3D macro later on */
-  memset (lsh, 0, sizeof (lsh));
-
-  /* get the current timelevel */
-  timelvl = 0;
+  /* Robin boundaries need the underlying grid coordinates */
+  sprintf (coord_system_name, "cart%dd", gdim);
+  if (CCTK_CoordSystemHandle (coord_system_name) < 0)
+  {
+    CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                "ApplyBndRobin: Couldn't get coordinates from '%s'",
+                coord_system_name);
+    return (-6);
+  }
+  x = GH->data[CCTK_CoordIndex (-1, "x", coord_system_name)][0];
+  y = GH->data[CCTK_CoordIndex (-1, "y", coord_system_name)][0];
+  z = GH->data[CCTK_CoordIndex (-1, "z", coord_system_name)][0];
 
-  /* Robin boundaries need the underlying radius coordinate */
   sprintf (coord_system_name, "spher%dd", gdim);
-  radius = GH->data[CCTK_CoordIndex (-1, "r", coord_system_name)][0];
+  if (CCTK_CoordSystemHandle (coord_system_name) < 0)
+  {
+    CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                "ApplyBndRobin: Couldn't get coordinates from '%s'",
+                coord_system_name);
+    return (-6);
+  }
+  r = GH->data[CCTK_CoordIndex (-1, "r", coord_system_name)][0];
 
   /* see if we have a symmetry array */
   sGHex = (SymmetryGHex *) CCTK_GHExtension (GH, "Symmetry");
 
+  /* get the decay rate as a double */
+  decay = 0.25 * (double) npow;
+
+  /* precompute the distance to all 8 neighbors in a 3D grid */
+  dist[0] = 0;                                    /* not used */
+  dist[1] = GH->cctk_delta_space[0];
+  dist[2] = GH->cctk_delta_space[1];
+  dist[3] = sqrt (SQR (dist[1]) + SQR (dist[2]));
+  dist[4] = GH->cctk_delta_space[2];
+  dist[5] = sqrt (SQR (dist[1]) + SQR (dist[4]));
+  dist[6] = sqrt (SQR (dist[2]) + SQR (dist[4]));
+  dist[7] = sqrt (SQR (dist[1]) + SQR (dist[2]) + SQR (dist[4]));
+
   /* now loop over all variables */
   for (var = first_var; var < first_var + num_vars; var++)
   {
@@ -650,17 +718,16 @@ static int ApplyBndRobin (const cGH *GH,
     memset (doBC, 1, sizeof (doBC));
     if (sGHex)
     {
-      for (i = 0; i < 2 * gdim; i++)
+      for (dim = 0; dim < 2 * gdim; dim++)
       {
-        doBC[i] = sGHex->GFSym[var][i] == GFSYM_NOSYM ||
-                  sGHex->GFSym[var][i] == GFSYM_UNSET;
+        doBC[dim] = sGHex->GFSym[var][dim] == GFSYM_NOSYM ||
+                    sGHex->GFSym[var][dim] == GFSYM_UNSET;
       }
     }
-    for (i = 0; i < gdim; i++)
+    for (dim = 0; dim < gdim; dim++)
     {
-      lsh[i]       = GH->cctk_lsh[i];
-      doBC[i*2]   &= GH->cctk_lsh[i] > 1 && GH->cctk_bbox[i*2];
-      doBC[i*2+1] &= GH->cctk_lsh[i] > 1 && GH->cctk_bbox[i*2+1];
+      doBC[dim*2]   &= GH->cctk_lsh[dim] > 1 && GH->cctk_bbox[dim*2];
+      doBC[dim*2+1] &= GH->cctk_lsh[dim] > 1 && GH->cctk_bbox[dim*2+1];
     }
 
     switch (CCTK_VarTypeI (var))
@@ -706,8 +773,8 @@ static int ApplyBndRobin (const cGH *GH,
 
       default:
         CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
-                    "ApplyBndRobin: Unsupported variable type %d for variable '%s'",
-                    CCTK_VarTypeI (var), CCTK_VarName (var));
+                    "ApplyBndRobin: Unsupported variable type %d for "
+                    "variable '%s'", CCTK_VarTypeI (var), CCTK_VarName (var));
         return (-4);
     }
   }