Support for new boundary API, which places an integer boundary_width

argument on calls to Select*ForBC*.  If a negative value is given then
the boundary condition will look into the table for a 2d element
integer array specifying the width of each boundary face.

A number bug fixes, including the restriction to only execute a
boundary condition on a single group at a time.  (consider that the
next group may have different staggering)

Improvements of debugging code.

Wrapper functions OldApplyBnd<BC> added, which support the old
d-element 'stencil width' array API.


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

1 files changed, 213 insertions, 66 deletions

diff --git a/src/RadiationBoundary.c b/src/RadiationBoundary.c
index 4df443b..88d6b1e 100644
--- a/src/RadiationBoundary.c
+++ b/src/RadiationBoundary.c
@@ -95,6 +95,15 @@ static int ApplyBndRadiative (const cGH *GH,
                               int first_var_to,
                               int first_var_from,
                               int num_vars);
+static int OldApplyBndRadiative (const cGH *GH,
+				 int stencil_dir,
+				 const int *stencil_alldirs,
+				 int dir,
+				 CCTK_REAL var0,
+				 CCTK_REAL speed,
+				 int first_var_to,
+				 int first_var_from,
+				 int num_vars);
 
 /********************************************************************
  ********************    External Routines   ************************
@@ -106,9 +115,10 @@ static int ApplyBndRadiative (const cGH *GH,
    @author     David Rideout
    @desc
                Top level function which is registered as handling 
-	       this boundary condition
+	       the Radiative boundary condition
    @enddesc
    @calls      ApplyBndRadiative
+
    @var        GH
    @vdesc      Pointer to CCTK grid hierarchy
    @vtype      const cGH *
@@ -130,6 +140,11 @@ static int ApplyBndRadiative (const cGH *GH,
    @vtype      int
    @vio        in
    @endvar
+   @var        widths
+   @vdesc      array of boundary widths for each variable
+   @vtype      int
+   @vio        in
+   @endvar
    @var        table_handles
    @vdesc      array of table handles which hold extra arguments
    @vtype      int
@@ -137,23 +152,23 @@ static int ApplyBndRadiative (const cGH *GH,
    @endvar
    @returntype int
    @returndesc
-               return code of @seeroutine ApplyBndRadiative <BR>
+               return code of @seeroutine ApplyBndRadiative
+               -21 error reading boundary width array from table
+	       -22 wrong size boundary width array in table
    @endreturndesc
 @@*/
 
-int BndRadiative(const cGH *GH, int num_vars, int *vars, int *faces,
-		 int *tables) 
+int BndRadiative(const cGH *GH, int num_vars, int *vars, int *faces, 
+		 int *widths, int *tables) 
 {
-  int i, j, k, value_type, value_size, err, 
-    gdim, max_gdim, retval;
+  int i, j, k, gi, gdim, max_gdim, value_type, value_size, err, retval;
   char *prev_time_level_name;
 
   /* variables to pass to ApplyBndRadiative */
-  /*int stencil_dir;*/ /* width of stencil in direction dir, not needed */
-  int *stencil_alldirs; /* width of stencil in all directions */
-  int dir; /* direction in which to apply bc */
+  int *width_alldirs;  /* width of boundary in all directions */
+  int dir;             /* direction in which to apply bc */
   CCTK_REAL limit, speed;
-  int prev_time_level; /* variable (index) which holds the previous time level */
+  int prev_time_level; /* variable index which holds the previous time level */
 
 #ifdef DEBUG
   printf("BndRadiative() got passed: GH=%p, num_vars=%d:\n", (const void *) GH,
@@ -168,15 +183,20 @@ int BndRadiative(const cGH *GH, int num_vars, int *vars, int *faces,
   /*  CCTK_WARN(0, "stopping code"); */
 #endif
 
-  retval = 0; stencil_alldirs = NULL; max_gdim = 0;
+  retval = 0; width_alldirs = NULL; max_gdim = 0;
 
   /* loop through variables, j at a time */
   for (i=0; i<num_vars; i+=j) {
-
     /* find other adjacent vars which are selected for identical bcs */
     j=1;
-    while (i+j<num_vars && vars[i+j]==vars[i]+j && tables[i+j]==tables[i] && 
-	   faces[i+j]==faces[i])
+    /* Since GFs are allowed to have different staggering, the best we
+       can do is find variables of the same group which are selected
+       for identical bcs.  If all GFs had the same staggering then we
+       could groups many GFs together. */
+    gi = CCTK_GroupIndexFromVarI (vars[i]);    
+    while (i+j<num_vars && CCTK_GroupIndexFromVarI(vars[i+j])==gi && 
+	   tables[i+j]==tables[i] && faces[i+j]==faces[i] && 
+	   widths[i+j]==widths[i])
     {
       ++j;
     }
@@ -187,27 +207,12 @@ int BndRadiative(const cGH *GH, int num_vars, int *vars, int *faces,
       CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
 		 "Faces specification %d for Radiative boundary conditions on "
 		 "%s is not implemented yet.  "
-		 "Applying radiative bcs to all (external) faces.", faces[i],
+		 "Applying Radiative bcs to all (external) faces.", faces[i],
 		 CCTK_VarName(vars[i]));
     }
     dir = 0; /* apply bc to all faces */
 
     /* Set up default arguments for ApplyBndRadiative */
-    /* Allocate memory and populate stencil width array */
-    gdim = CCTK_GroupDimFromVarI(vars[i]);
-    if (!stencil_alldirs)
-    {
-      stencil_alldirs = (int *) malloc(gdim*sizeof(int));
-      max_gdim = gdim;
-    } else if (gdim > max_gdim)
-    {
-      realloc(stencil_alldirs, gdim*sizeof(int));
-      max_gdim = gdim;
-    }
-    for (k=0; k<gdim; ++k)
-    {
-      stencil_alldirs[k] = 1;
-    }
     /* Defaults for remainder of arguments */
     limit = 0.;
     prev_time_level = vars[i];
@@ -253,14 +258,52 @@ int BndRadiative(const cGH *GH, int num_vars, int *vars, int *faces,
 	}
       }
 
-      /* Stencil width array */
-      Util_TableGetIntArray(tables[i], gdim, stencil_alldirs, "STENCIL WIDTH");
-    
       /* Wave speed */
       Util_TableGetReal(tables[i], &speed, "SPEED");
     }
 
-    if ((retval = ApplyBndRadiative(GH, 0, stencil_alldirs, dir, limit, speed,
+    /* Determine boundary width on all faces */
+    /* allocate memory for buffer */
+    gdim = CCTK_GroupDimI(gi);
+    if (!width_alldirs)
+    {
+      width_alldirs = (int *) malloc(2*gdim*sizeof(int));
+      max_gdim = gdim;
+    } else if (gdim > max_gdim)
+    {
+      width_alldirs = realloc(width_alldirs, gdim*sizeof(int));
+      max_gdim = gdim;
+    }
+
+    /* fill it with values, either from table or the boundary_width 
+       parameter */
+    if (widths[i]<0)
+    {
+      err = Util_TableGetIntArray(tables[i], gdim, width_alldirs, 
+				  "BOUNDARY WIDTH");
+      if (err<0)
+      {
+	CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
+		   "Error %d when reading boundary width array from table "
+		   "for %s", err, CCTK_VarName(vars[i]));
+	return -21;
+      } else if (err!=2*gdim)
+      {
+	CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
+		   "Boundary width array for %s has %d elements, but %d "
+		   "expected", CCTK_VarName(vars[i]), err, 2*gdim);
+	return -22;
+      }
+    } else
+    {
+      for (k=0; k<2*gdim; ++k)
+      {
+	width_alldirs[k] = widths[i];
+      }
+    }
+
+    /* Apply the boundary condition */
+    if ((retval = ApplyBndRadiative(GH, 0, width_alldirs, dir, limit, speed,
 				    vars[i], prev_time_level, j)) < 0)
     {
       CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
@@ -270,7 +313,7 @@ int BndRadiative(const cGH *GH, int num_vars, int *vars, int *faces,
 #ifdef DEBUG
   printf("BndRadiative(): returning %d\n",retval);
 #endif
-  free(stencil_alldirs);
+  free(width_alldirs);
 
   return retval;
 }
@@ -502,7 +545,7 @@ int BndRadiativeGI (const cGH *GH,
   first_vi_from = CCTK_FirstVarIndexI (gi_from);
   if (first_vi_to >= 0 && first_vi_from >= 0)
   {
-    retval = ApplyBndRadiative (GH, -1, stencil, 0, var0, speed,
+    retval = OldApplyBndRadiative (GH, -1, stencil, 0, var0, speed,
                                 first_vi_to, first_vi_from,
                                 CCTK_NumVarsInGroupI (gi_to));
   }
@@ -872,7 +915,7 @@ int BndRadiativeVI (const cGH *GH,
   num_vars = CCTK_NumVars ();
   if (vi_to >= 0 && vi_to < num_vars && vi_from >= 0 && vi_from < num_vars)
   {
-    retval = ApplyBndRadiative (GH, -1, stencil, 0, var0, speed,
+    retval = OldApplyBndRadiative (GH, -1, stencil, 0, var0, speed,
                                 vi_to, vi_from, 1);
   }
   else
@@ -1408,7 +1451,7 @@ void CCTK_FCALL CCTK_FNAME (BndRadiativeVN)
   /* Upper x-bound */                                                         \
   if (doBC[1])                                                                \
   {                                                                           \
-    UPPER_RADIATIVE_BOUNDARY_3D (lsh[0] - stencil[0], 0, 0,                   \
+    UPPER_RADIATIVE_BOUNDARY_3D (lsh[0] - stencil[1], 0, 0,                   \
                                  coords[0], coords[MAXDIM], offset[0],        \
                                  var_to, var_from, cctk_type);                \
   }                                                                           \
@@ -1416,7 +1459,7 @@ void CCTK_FCALL CCTK_FNAME (BndRadiativeVN)
   /* Lower y-bound */                                                         \
   if (doBC[2])                                                                \
   {                                                                           \
-    LOWER_RADIATIVE_BOUNDARY_3D (lsh[0], stencil[1], lsh[2],                  \
+    LOWER_RADIATIVE_BOUNDARY_3D (lsh[0], stencil[2], lsh[2],                  \
                                  coords[1], coords[MAXDIM], offset[1],        \
                                  var_to, var_from, cctk_type);                \
   }                                                                           \
@@ -1424,7 +1467,7 @@ void CCTK_FCALL CCTK_FNAME (BndRadiativeVN)
   /* Upper y-bound */                                                         \
   if (doBC[3])                                                                \
   {                                                                           \
-    UPPER_RADIATIVE_BOUNDARY_3D (0, lsh[1] - stencil[1], 0,                   \
+    UPPER_RADIATIVE_BOUNDARY_3D (0, lsh[1] - stencil[3], 0,                   \
                                  coords[1], coords[MAXDIM], offset[1],        \
                                  var_to, var_from, cctk_type);                \
   }                                                                           \
@@ -1432,7 +1475,7 @@ void CCTK_FCALL CCTK_FNAME (BndRadiativeVN)
   /* Lower z-bound */                                                         \
   if (doBC[4])                                                                \
   {                                                                           \
-    LOWER_RADIATIVE_BOUNDARY_3D (lsh[0], lsh[1], stencil[2],                  \
+    LOWER_RADIATIVE_BOUNDARY_3D (lsh[0], lsh[1], stencil[4],                  \
                                  coords[2], coords[MAXDIM], offset[2],        \
                                  var_to, var_from, cctk_type);                \
   }                                                                           \
@@ -1440,7 +1483,7 @@ void CCTK_FCALL CCTK_FNAME (BndRadiativeVN)
   /* Upper z-bound */                                                         \
   if (doBC[5])                                                                \
   {                                                                           \
-    UPPER_RADIATIVE_BOUNDARY_3D (0, 0, lsh[2] - stencil[2],                   \
+    UPPER_RADIATIVE_BOUNDARY_3D (0, 0, lsh[2] - stencil[5],                   \
                                  coords[2], coords[MAXDIM], offset[2],        \
                                  var_to, var_from, cctk_type);                \
   }                                                                           \
@@ -1466,13 +1509,13 @@ void CCTK_FCALL CCTK_FNAME (BndRadiativeVN)
    @vtype      const cGH *
    @vio        in
    @endvar
-   @var        stencil_dir
-   @vdesc      stencil width in direction dir
+   @var        width_dir
+   @vdesc      boundary width in direction dir
    @vtype      int
    @vio        in
    @endvar
-   @var        stencil_alldirs
-   @vdesc      stencil widths for all directions
+   @var        in_widths
+   @vdesc      boundary widths for all directions
    @vtype      int [ dimension of variable(s) ]
    @vio        in
    @endvar
@@ -1506,7 +1549,7 @@ void CCTK_FCALL CCTK_FNAME (BndRadiativeVN)
    @vtype      int
    @vio        in
    @endvar
-               CCTK_VarTypeI
+   @calls      CCTK_VarTypeI
                CCTK_GroupDimFromVarI
                RADIATIVE_BOUNDARY
    @history
@@ -1521,15 +1564,15 @@ void CCTK_FCALL CCTK_FNAME (BndRadiativeVN)
                 0 for success
                -1 if abs(direction) is greater than variables' dimension
                -2 if variable dimension is not supported
-               -3 if NULL pointer passed as stencil array
+               -3 if NULL pointer passed as boundary width array
                -4 if variable type is not supported
                -5 if variable dimension is other than 3D
                -6 if a coordinate is not found
    @endreturndesc
 @@*/
 static int ApplyBndRadiative (const cGH *GH,
-                              int stencil_dir,
-                              const int *stencil_alldirs,
+                              int width_dir,
+                              const int *in_widths,
                               int dir,
                               CCTK_REAL var0,
                               CCTK_REAL speed,
@@ -1544,10 +1587,9 @@ static int ApplyBndRadiative (const cGH *GH,
   SymmetryGHex *sGHex;
   char coord_system_name[10];
   CCTK_REAL dxyz[MAXDIM], rho[MAXDIM], *coords[MAXDIM+1];
-  int doBC[2*MAXDIM], stencil[MAXDIM], offset[MAXDIM];
+  int doBC[2*MAXDIM], widths[2*MAXDIM], offset[MAXDIM];
   CCTK_REAL dtv, dtvh, dtvvar0, dtvvar0H;
 
-
   /* check the direction parameter */
   if (abs (dir) > MAXDIM)
   {
@@ -1569,19 +1611,20 @@ static int ApplyBndRadiative (const cGH *GH,
     return (-2);
   }
 
-  /* set up stencil width array */
+  /* set up boundary width array */
   if (dir)
   {
-    stencil[abs (dir) - 1] = stencil_dir;
+    widths[2*(abs(dir)-1)] = width_dir;
+    widths[2*(abs(dir)-1)+1] = width_dir;
   }
-  else if (stencil_alldirs)
+  else if (in_widths)
   {
-    memcpy (stencil, stencil_alldirs, gdim * sizeof (int));
+    memcpy (widths, in_widths, 2* gdim * sizeof (int));
   }
   else
   {
     CCTK_WARN (1, "ApplyBndRadiative: NULL pointer passed "
-                  "for stencil width array");
+                  "for boundary width array");
     return (-3);
   }
 
@@ -1674,26 +1717,26 @@ static int ApplyBndRadiative (const cGH *GH,
     switch (CCTK_VarTypeI (var_to))
     {
       case CCTK_VARIABLE_CHAR:
-        RADIATIVE_BOUNDARY (GH->cctk_lsh, doBC, stencil, coords, offset,
+        RADIATIVE_BOUNDARY (GH->cctk_lsh, doBC, widths, coords, offset,
                             GH->data[var_to][timelvl_to],
                             GH->data[var_from][timelvl_from], gdim, CCTK_CHAR);
         break;
 
       case CCTK_VARIABLE_INT:
-        RADIATIVE_BOUNDARY (GH->cctk_lsh, doBC, stencil, coords, offset,
+        RADIATIVE_BOUNDARY (GH->cctk_lsh, doBC, widths, coords, offset,
                             GH->data[var_to][timelvl_to],
                             GH->data[var_from][timelvl_from], gdim, CCTK_INT);
         break;
 
       case CCTK_VARIABLE_REAL:
-        RADIATIVE_BOUNDARY (GH->cctk_lsh, doBC, stencil, coords, offset,
+        RADIATIVE_BOUNDARY (GH->cctk_lsh, doBC, widths, coords, offset,
                             GH->data[var_to][timelvl_to],
                             GH->data[var_from][timelvl_from], gdim, CCTK_REAL);
         break;
 
 #ifdef CCTK_INT2
       case CCTK_VARIABLE_INT2:
-        RADIATIVE_BOUNDARY (GH->cctk_lsh, doBC, stencil, coords, offset,
+        RADIATIVE_BOUNDARY (GH->cctk_lsh, doBC, widths, coords, offset,
                             GH->data[var_to][timelvl_to],
                             GH->data[var_from][timelvl_from], gdim, CCTK_INT2);
         break;
@@ -1701,7 +1744,7 @@ static int ApplyBndRadiative (const cGH *GH,
 
 #ifdef CCTK_INT4
       case CCTK_VARIABLE_INT4:
-        RADIATIVE_BOUNDARY (GH->cctk_lsh, doBC, stencil, coords, offset,
+        RADIATIVE_BOUNDARY (GH->cctk_lsh, doBC, widths, coords, offset,
                             GH->data[var_to][timelvl_to],
                             GH->data[var_from][timelvl_from], gdim, CCTK_INT4);
         break;
@@ -1709,7 +1752,7 @@ static int ApplyBndRadiative (const cGH *GH,
 
 #ifdef CCTK_INT8
       case CCTK_VARIABLE_INT8:
-        RADIATIVE_BOUNDARY (GH->cctk_lsh, doBC, stencil, coords, offset,
+        RADIATIVE_BOUNDARY (GH->cctk_lsh, doBC, widths, coords, offset,
                             GH->data[var_to][timelvl_to],
                             GH->data[var_from][timelvl_from], gdim, CCTK_INT8);
         break;
@@ -1717,7 +1760,7 @@ static int ApplyBndRadiative (const cGH *GH,
 
 #ifdef CCTK_REAL4
       case CCTK_VARIABLE_REAL4:
-        RADIATIVE_BOUNDARY (GH->cctk_lsh, doBC, stencil, coords, offset,
+        RADIATIVE_BOUNDARY (GH->cctk_lsh, doBC, widths, coords, offset,
                             GH->data[var_to][timelvl_to],
                             GH->data[var_from][timelvl_from], gdim, CCTK_REAL4);
         break;
@@ -1725,7 +1768,7 @@ static int ApplyBndRadiative (const cGH *GH,
 
 #ifdef CCTK_REAL8
       case CCTK_VARIABLE_REAL8:
-        RADIATIVE_BOUNDARY (GH->cctk_lsh, doBC, stencil, coords, offset,
+        RADIATIVE_BOUNDARY (GH->cctk_lsh, doBC, widths, coords, offset,
                             GH->data[var_to][timelvl_to],
                             GH->data[var_from][timelvl_from], gdim, CCTK_REAL8);
         break;
@@ -1733,7 +1776,7 @@ static int ApplyBndRadiative (const cGH *GH,
 
 #ifdef CCTK_REAL16
       case CCTK_VARIABLE_REAL16:
-        RADIATIVE_BOUNDARY (GH->cctk_lsh, doBC, stencil, coords, offset,
+        RADIATIVE_BOUNDARY (GH->cctk_lsh, doBC, widths, coords, offset,
                             GH->data[var_to][timelvl_to],
                             GH->data[var_from][timelvl_from], gdim,CCTK_REAL16);
         break;
@@ -1749,3 +1792,107 @@ static int ApplyBndRadiative (const cGH *GH,
 
   return (0);
 }
+
+
+/*@@
+   @routine    OldApplyBndRadiative
+   @date       5 May 2003
+   @author     David Rideout
+   @desc
+               The new boundary API expects a 2d-element array for the
+               boundary_widths (d=dimension of grid variable), while
+               the old API expects a d-element array.  This function
+               converts the old array to the new format.
+   @enddesc
+
+   @var        GH
+   @vdesc      Pointer to CCTK grid hierarchy
+   @vtype      const cGH *
+   @vio        in
+   @endvar
+   @var        width_dir
+   @vdesc      boundary width in direction dir
+   @vtype      int
+   @vio        in
+   @endvar
+   @var        in_widths
+   @vdesc      boundary widths for all directions
+   @vtype      int [ dimension of variable(s) ]
+   @vio        in
+   @endvar
+   @var        dir
+   @vdesc      direction to copy boundaries (0 for copying all directions)
+   @vtype      int
+   @vio        in
+   @endvar
+   @var        var0
+   @vdesc      asymptotic value of function at infinity
+   @vtype      CCTK_REAL
+   @vio        in
+   @endvar
+   @var        speed
+   @vdesc      wave speed
+   @vtype      CCTK_REAL
+   @vio        in
+   @endvar
+   @var        first_var_to
+   @vdesc      index of first variable to copy boundaries to
+   @vtype      int
+   @vio        in
+   @endvar
+   @var        first_var_from
+   @vdesc      index of first variable to copy boundaries from
+   @vtype      int
+   @vio        in
+   @endvar
+   @var        num_vars
+   @vdesc      number of variables
+   @vtype      int
+   @vio        in
+   @endvar
+
+   @calls      CCTK_GroupIndexFromVarI
+               ApplyBndRadiative
+   @returntype int
+   @returndesc
+               returncode from @seeroutine ApplyBndRadiative
+   @endreturndesc
+@@*/
+
+int OldApplyBndRadiative(const cGH *GH,
+			 int width_dir,
+			 const int *stencil_alldirs,
+			 int dir,
+			 CCTK_REAL var0,
+			 CCTK_REAL speed,
+			 int first_var_to,
+			 int first_var_from,
+			 int num_vars)
+{
+  int retval, *boundary_widths, dim, i;
+  static int warned;
+
+  /* Convert stencil_alldirs to new format */
+  dim = CCTK_GroupDimFromVarI(first_var_to);
+  boundary_widths = (int *) malloc(2*dim*sizeof(int));
+  for (i=0; i<2*dim; ++i)
+  {
+    boundary_widths[i] = stencil_alldirs[i/2];
+  }
+
+  /* Bug people for using the old interface */
+  if (!warned)
+  {
+    CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
+                "Copied older d-element array of boundary widths into the "
+		"newer 2d-element format.  Please use the new boundary "
+		"interface to avoid this.");
+    warned = 1;
+  }
+
+  /* Call ApplyBnd... with new boundary width array */
+  retval = ApplyBndRadiative(GH, width_dir, boundary_widths, dir, var0, speed, 
+			     first_var_to, first_var_from, num_vars);
+
+  return retval;
+}