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, 176 insertions, 63 deletions

diff --git a/src/FlatBoundary.c b/src/FlatBoundary.c
index 0ac4727..53c2322 100644
--- a/src/FlatBoundary.c
+++ b/src/FlatBoundary.c
@@ -38,7 +38,12 @@ static int ApplyBndFlat (const cGH *GH,
                          int dir,
                          int first_var,
                          int num_vars);
-
+static int OldApplyBndFlat (const cGH *GH,
+			    int stencil_dir,
+			    const int *stencil_alldirs,
+			    int dir,
+			    int first_var,
+			    int num_vars);
 
 /********************************************************************
  ********************    External Routines   ************************
@@ -52,6 +57,7 @@ static int ApplyBndFlat (const cGH *GH,
 	       the Flat boundary condition
    @enddesc
    @calls      ApplyBndFlat
+
    @var        GH
    @vdesc      Pointer to CCTK grid hierarchy
    @vtype      const cGH *
@@ -73,6 +79,11 @@ static int ApplyBndFlat (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
@@ -80,28 +91,35 @@ static int ApplyBndFlat (const cGH *GH,
    @endvar
    @returntype int
    @returndesc
-               return code of @seeroutine ApplyBndScalar
+               return code of @seeroutine ApplyBndFlat
+               -21 error reading boundary width array from table
+	       -22 wrong size boundary width array in table
    @endreturndesc
 @@*/
 
-int BndFlat(const cGH *GH, int num_vars, int *vars, int *faces, int *tables) 
+int BndFlat(const cGH *GH, int num_vars, int *vars, int *faces, int *widths,
+	    int *tables)
 {
-  int i, j, k, gdim, max_gdim, err, retval;
+  int i, j, k, gi, gdim, max_gdim, err, retval;
 
   /* variables to pass to ApplyBndFlat */
-  /*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 */
 
-  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;
     }
@@ -112,50 +130,59 @@ int BndFlat(const cGH *GH, int num_vars, int *vars, int *faces, int *tables)
       CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
 		 "Faces specification %d for Flat boundary conditions on "
 		 "%s is not implemented yet.  "
-		 "Applying radiative bcs to all (external) faces.", faces[i],
+		 "Applying Flat bcs to all (external) faces.", faces[i],
 		 CCTK_VarName(vars[i]));
     }
     dir = 0; /* apply bc to all faces */
 
-    /* Set up default arguments for ApplyBndFlat */
-    /* Allocate memory and populate stencil width array */
-    gdim = CCTK_GroupDimFromVarI(vars[i]);
-    if (!stencil_alldirs)
+    /* Determine boundary width on all faces */
+    /* allocate memory for buffer */
+    gdim = CCTK_GroupDimI(gi);
+    if (!width_alldirs)
     {
-      stencil_alldirs = (int *) malloc(gdim*sizeof(int));
+      width_alldirs = (int *) malloc(2*gdim*sizeof(int));
       max_gdim = gdim;
     } else if (gdim > max_gdim)
     {
-      realloc(stencil_alldirs, gdim*sizeof(int));
+      width_alldirs = realloc(width_alldirs, gdim*sizeof(int));
       max_gdim = gdim;
     }
-    for (k=0; k<gdim; ++k)
-    {
-      stencil_alldirs[k] = 1;
-    }
 
-    /* Look on table for possible non-default arguments
-     * (If any of these table look-ups fail, the value will be unchanged 
-     * from its default value)
-     */
-    /* Stencil width array */
-    err = Util_TableGetIntArray(tables[i], gdim, stencil_alldirs, 
-				"STENCIL WIDTH");
-    if (err == UTIL_ERROR_BAD_HANDLE)
+    /* fill it with values, either from table or the boundary_width 
+       parameter */
+    if (widths[i]<0)
     {
-      CCTK_VWarn(5, __LINE__, __FILE__, CCTK_THORNSTRING,
-		 "Invalid table handle passed for Flat boundary "
-		 "conditions for %s.  Using all default values.", 
-		 CCTK_VarName(vars[i]));
+      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];
+      }
     }
 
-    if ((retval = ApplyBndFlat(GH, 0, stencil_alldirs, dir, vars[i], j)) < 0)
+    /* Apply the boundary condition */
+    if ((retval = ApplyBndFlat(GH, 0, width_alldirs, dir, vars[i], j)) < 0)
     {
       CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
 		 "ApplyBndFlat() returned %d", retval);
     }
   }
-  free(stencil_alldirs);
+  free(width_alldirs);
 
   return retval;
 }
@@ -324,7 +351,7 @@ int BndFlatGI (const cGH *GH,
   first_vi = CCTK_FirstVarIndexI (gi);
   if (first_vi >= 0)
   {
-    retval = ApplyBndFlat (GH, -1, stencil, 0, first_vi,
+    retval = OldApplyBndFlat (GH, -1, stencil, 0, first_vi,
                            CCTK_NumVarsInGroupI (gi));
   }
   else
@@ -598,7 +625,7 @@ int BndFlatVI (const cGH *GH,
 
   if (vi >= 0 && vi < CCTK_NumVars ())
   {
-    retval = ApplyBndFlat (GH, -1, stencil, 0, vi, 1);
+    retval = OldApplyBndFlat (GH, -1, stencil, 0, vi, 1);
   }
   else
   {
@@ -846,13 +873,13 @@ void CCTK_FCALL CCTK_FNAME (BndFlatVN)
    @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
@@ -889,19 +916,20 @@ void CCTK_FCALL CCTK_FNAME (BndFlatVN)
                 0 for success
                -1 if dimension is not supported
                -2 if direction parameter is invalid
-               -3 if stencil width array parameter is NULL
+               -3 if boundary width array parameter is NULL
    @endreturndesc
 @@*/
 static int ApplyBndFlat (const cGH *GH,
-                         int stencil_dir,
-                         const int *stencil_alldirs,
+                         int width_dir,
+                         const int *in_widths,
                          int dir,
                          int first_var,
                          int num_vars)
 {
   int i, j, k;
   int var, vtypesize, gindex, gdim, timelvl;
-  int doBC[2*MAXDIM], dstag[MAXDIM], lsh[MAXDIM], lssh[MAXDIM], stencil[MAXDIM];
+  int doBC[2*MAXDIM], dstag[MAXDIM], lsh[MAXDIM], lssh[MAXDIM], 
+    widths[2*MAXDIM];
   SymmetryGHex *sGHex;
 
 
@@ -929,18 +957,20 @@ static int ApplyBndFlat (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, "ApplyBndFlat: NULL pointer passed for stencil width array");
+    CCTK_WARN (1, "ApplyBndFlat: NULL pointer passed for boundary width "
+	       "array");
     return (-3);
   }
 
@@ -1005,11 +1035,11 @@ static int ApplyBndFlat (const cGH *GH,
       }
 #endif /* DEBUG_BOUNDARY */
       /* lower x */
-      FLAT_BOUNDARY (doBC[0], stencil[0], lssh[1], lssh[2],
-                     i, j, k, stencil[0], j, k);
+      FLAT_BOUNDARY (doBC[0], widths[0], lssh[1], lssh[2],
+                     i, j, k, widths[0], j, k);
       /* upper x */
-      FLAT_BOUNDARY (doBC[1], stencil[0], lssh[1], lssh[2],
-                     lssh[0]-i-1, j, k, lssh[0]-stencil[0]-1, j, k);
+      FLAT_BOUNDARY (doBC[1], widths[1], lssh[1], lssh[2],
+                     lssh[0]-i-1, j, k, lssh[0]-widths[1]-1, j, k);
 
     }
     if (gdim > 1)
@@ -1025,11 +1055,11 @@ static int ApplyBndFlat (const cGH *GH,
       }
 #endif /* DEBUG_BOUNDARY */
       /* lower y */
-      FLAT_BOUNDARY (doBC[2], lssh[0], stencil[1], lssh[2],
-                     i, j, k, i, stencil[1], k);
+      FLAT_BOUNDARY (doBC[2], lssh[0], widths[2], lssh[2],
+                     i, j, k, i, widths[2], k);
       /* upper y */
-      FLAT_BOUNDARY (doBC[3], lssh[0], stencil[1], lssh[2],
-                     i, lssh[1]-j-1, k, i, lssh[1]-stencil[1]-1, k);
+      FLAT_BOUNDARY (doBC[3], lssh[0], widths[3], lssh[2],
+                     i, lssh[1]-j-1, k, i, lssh[1]-widths[3]-1, k);
     }
     if (gdim > 2)
     {
@@ -1044,13 +1074,96 @@ static int ApplyBndFlat (const cGH *GH,
       }
 #endif /* DEBUG_BOUNDARY */
       /* lower z */
-      FLAT_BOUNDARY (doBC[4], lssh[0], lssh[1], stencil[2],
-                     i, j, k, i, j, stencil[2]);
+      FLAT_BOUNDARY (doBC[4], lssh[0], lssh[1], widths[4],
+                     i, j, k, i, j, widths[4]);
       /* upper z */
-      FLAT_BOUNDARY (doBC[5], lssh[0], lssh[1], stencil[2],
-                     i, j, lssh[2]-k-1, i, j, lssh[2]-stencil[2]-1);
+      FLAT_BOUNDARY (doBC[5], lssh[0], lssh[1], widths[5],
+                     i, j, lssh[2]-k-1, i, j, lssh[2]-widths[5]-1);
     }
   }
 
   return(0);
 }
+
+/*@@
+   @routine    OldApplyBndFlat
+   @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        stencil_dir
+   @vdesc      stencil width in direction dir
+   @vtype      int
+   @vio        in
+   @endvar
+   @var        stencil_alldirs
+   @vdesc      stencil widths for all directions
+   @vtype      int [ dimension of variable(s) ]
+   @vio        in
+   @endvar
+   @var        dir
+   @vdesc      direction to set boundaries (0 for setting all directions)
+   @vtype      int
+   @vio        in
+   @endvar
+   @var        first_var
+   @vdesc      index of first variable to apply boundaries to
+   @vtype      int
+   @vio        in
+   @endvar
+   @var        num_vars
+   @vdesc      number of variables
+   @vtype      int
+   @vio        in
+   @endvar
+
+   @calls      CCTK_GroupIndexFromVarI
+               CCTK_GroupDimI
+               ApplyBndFlat
+   @history
+   @returntype int
+   @returndesc
+               returncode from @seeroutine ApplyBndFlat
+   @endreturndesc
+@@*/
+
+int OldApplyBndFlat(const cGH *GH, int stencil_dir, const int *stencil_alldirs,
+		    int dir, int first_var, int num_vars)
+{
+  int retval, *boundary_widths, dim, i;
+  static int warned;
+
+  /* Convert stencil_alldirs to new format */
+  dim = CCTK_GroupDimFromVarI(first_var);
+  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 = ApplyBndFlat(GH, stencil_dir, boundary_widths, dir, first_var, 
+			num_vars);
+
+  return retval;
+}