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, 169 insertions, 43 deletions

diff --git a/src/RobinBoundary.c b/src/RobinBoundary.c
index 3d2d67e..7c8020c 100644
--- a/src/RobinBoundary.c
+++ b/src/RobinBoundary.c
@@ -37,6 +37,12 @@ static int ApplyBndRobin (const cGH *GH,
                           int npow,
                           int first_var,
                           int num_vars);
+static int OldApplyBndRobin (const cGH *GH,
+			     const int *stencil,
+			     CCTK_REAL finf,
+			     int npow,
+			     int first_var,
+			     int num_vars);
 
 
 /********************************************************************
@@ -72,6 +78,11 @@ static int ApplyBndRobin (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
@@ -79,32 +90,40 @@ static int ApplyBndRobin (const cGH *GH,
    @endvar
    @returntype int
    @returndesc
-               return code of @seeroutine ApplyBndRobin <BR>
+               return code of @seeroutine ApplyBndRobin
+               -21 error reading boundary width array from table
+	       -22 wrong size boundary width array in table
    @endreturndesc
 @@*/
 
-int BndRobin(const cGH *GH, int num_vars, int *vars, int *faces, int *tables) 
+int BndRobin(const cGH *GH, int num_vars, int *vars, int *faces, int *widths, 
+	     int *tables) 
 {
-  int i, j, k, err, gdim, max_gdim, retval;
+  int i, j, k, gi, err, gdim, max_gdim, retval;
 
   /* variables to pass to ApplyBndRobin */
-  int *stencil; /* width of stencil in all directions */
-  CCTK_REAL finf; /* value of function at infinity */
-  int npow; /* decay rate */
+  int *width_alldirs; /* width of boundary in all directions */
+  CCTK_REAL finf;     /* value of function at infinity */
+  int npow;           /* decay rate */
 
 #ifdef DEBUG
   printf("BndRobin(): got passed GH=%p, num_vars=%d, var_indices[0]=%d, table_handles[0]=%d\n", (const void *) GH, num_vars, var_indices[0], table_handles[0]);
 #endif
 
-  retval = 0; stencil = 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;
     }
@@ -115,27 +134,11 @@ int BndRobin(const cGH *GH, int num_vars, int *vars, int *faces, int *tables)
       CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
 		 "Faces specification %d for Robin boundary conditions on "
 		 "%s is not implemented yet.  "
-		 "Applying radiative bcs to all (external) faces.", faces[i],
+		 "Applying Robin bcs to all (external) faces.", faces[i],
 		 CCTK_VarName(vars[i]));
     }
 
     /* Set up default arguments for ApplyBndRobin */
-    /* Allocate memory and populate stencil width array */
-    gdim = CCTK_GroupDimFromVarI(vars[i]);
-    if (!stencil)
-    {
-      stencil = (int *) malloc(gdim*sizeof(int));
-      max_gdim = gdim;
-    } else if (gdim > max_gdim)
-    {
-      realloc(stencil, gdim*sizeof(int));
-      max_gdim = gdim;
-    }
-    for (k=0; k<gdim; ++k)
-    {
-      stencil[k] = 1;
-    }
-    /* Defaults for remainder of arguments */
     finf = 1.;
     npow = 1;
 
@@ -153,14 +156,52 @@ int BndRobin(const cGH *GH, int num_vars, int *vars, int *faces, int *tables)
 		 CCTK_VarName(vars[i]));
     } else
     { 
-      /* Stencil width array */
-      Util_TableGetIntArray(tables[i], gdim, stencil, "STENCIL WIDTH");
-    
       /* Decay power */
       Util_TableGetInt(tables[i], &npow, "DECAY POWER");
     }
 
-    if ((retval = ApplyBndRobin(GH, stencil, finf, npow, vars[i], j)) < 0)
+    /* 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 = ApplyBndRobin(GH, width_alldirs, finf, npow, vars[i], j)) < 0)
     {
       CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
 		 "ApplyBndRobin() returned %d", retval);
@@ -169,7 +210,7 @@ int BndRobin(const cGH *GH, int num_vars, int *vars, int *faces, int *tables)
 #ifdef DEBUG
   printf("BndRobin(): returning %d\n",retval);
 #endif
-  free(stencil);
+  free(width_alldirs);
 
   return retval;
 }
@@ -263,7 +304,7 @@ int BndRobinGI (const cGH *GH,
   first_vi = CCTK_FirstVarIndexI (gi);
   if (first_vi >= 0)
   {
-    retval = ApplyBndRobin (GH, stencil, finf, npow, first_vi,
+    retval = OldApplyBndRobin (GH, stencil, finf, npow, first_vi,
                             CCTK_NumVarsInGroupI (gi));
   }
   else
@@ -419,7 +460,7 @@ int BndRobinVI (const cGH *GH,
 
   if (vi >= 0 && vi < CCTK_NumVars ())
   {
-    retval = ApplyBndRobin (GH, stencil, finf, npow, vi, 1);
+    retval = OldApplyBndRobin (GH, stencil, finf, npow, vi, 1);
   }
   else
   {
@@ -712,8 +753,8 @@ void CCTK_FCALL CCTK_FNAME (BndRobinVN)
    @vtype      const cGH *
    @vio        in
    @endvar
-   @var        stencil
-   @vdesc      stencil width array
+   @var        in_widths
+   @vdesc      boundary width array
    @vtype      int [ dimension of variable ]
    @vio        in
    @endvar
@@ -752,7 +793,7 @@ void CCTK_FCALL CCTK_FNAME (BndRobinVN)
    @returndesc
                 0 for success
                -1 if variable dimension is not supported
-               -2 if NULL pointer passed as stencil width array
+               -2 if NULL pointer passed as boundary width array
                -3 if stencil width is other than 1
                -4 if variable type is not supported
                -5 if variable dimension is other than 3D
@@ -760,7 +801,7 @@ void CCTK_FCALL CCTK_FNAME (BndRobinVN)
    @endreturndesc
 @@*/
 static int ApplyBndRobin (const cGH *GH,
-                          const int *stencil,
+                          const int *in_widths,
                           CCTK_REAL finf,
                           int npow,
                           int first_var,
@@ -787,16 +828,17 @@ static int ApplyBndRobin (const cGH *GH,
     return (-1);
   }
 
-  /* check the stencil width */
-  if (! stencil)
+  /* check the boundary width */
+  if (! in_widths)
   {
-    CCTK_WARN (1, "ApplyBndRobin: NULL pointer passed for stencil width array");
+    CCTK_WARN (1, "ApplyBndRobin: NULL pointer passed for boundary width "
+	       "array");
     return (-2);
   }
 
-  for (dim = 0; dim < gdim; dim++)
+  for (dim = 0; dim < 2*gdim; dim++)
   {
-    if (stencil[dim] != 1)
+    if (in_widths[dim] != 1)
     {
       CCTK_WARN (1, "ApplyBndRobin: Stencil width must be 1 "
                     "for Robin boundary conditions");
@@ -917,3 +959,87 @@ static int ApplyBndRobin (const cGH *GH,
 
   return(0);
 }
+
+
+/*@@
+   @routine    OldApplyBndRobin
+   @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        in_widths
+   @vdesc      boundary width array
+   @vtype      int [ dimension of variable ]
+   @vio        in
+   @endvar
+   @var        finf
+   @vdesc      value of f at infinity
+   @vtype      CCTK_REAL
+   @vio        in
+   @endvar
+   @var        npow
+   @vdesc      power of decay rate
+   @vtype      int
+   @vio        in
+   @endvar
+   @var        first_var
+   @vdesc      index of first variable to apply boundary conditions to
+   @vtype      int
+   @vio        in
+   @endvar
+   @var        num_vars
+   @vdesc      number of variables
+   @vtype      int
+   @vio        in
+   @endvar
+
+   @calls      CCTK_GroupIndexFromVarI
+               ApplyBndRobin
+   @returntype int
+   @returndesc
+               returncode from @seeroutine ApplyBndRobin
+   @endreturndesc
+@@*/
+static int OldApplyBndRobin (const cGH *GH,
+			     const int *in_widths,
+			     CCTK_REAL finf,
+			     int npow,
+			     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] = in_widths[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 = ApplyBndRobin(GH, boundary_widths, finf, npow, first_var, num_vars);
+
+  return retval;
+}