diff options
Diffstat (limited to 'src/ScalarBoundary.c')
-rw-r--r-- | src/ScalarBoundary.c | 238 |
1 files changed, 187 insertions, 51 deletions
diff --git a/src/ScalarBoundary.c b/src/ScalarBoundary.c index cebd85e..1a7cb03 100644 --- a/src/ScalarBoundary.c +++ b/src/ScalarBoundary.c @@ -37,6 +37,13 @@ static int ApplyBndScalar (const cGH *GH, CCTK_REAL scalar, int first_var, int num_vars); +static int OldApplyBndScalar (const cGH *GH, + int stencil_dir, + const int *stencil_alldirs, + int dir, + CCTK_REAL scalar, + int first_var, + int num_vars); /******************************************************************** ******************** External Routines ************************ @@ -51,6 +58,7 @@ static int ApplyBndScalar (const cGH *GH, the Scalar boundary condition @enddesc @calls ApplyBndScalar + @var GH @vdesc Pointer to CCTK grid hierarchy @vtype const cGH * @@ -72,6 +80,11 @@ static int ApplyBndScalar (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 +93,35 @@ static int ApplyBndScalar (const cGH *GH, @returntype int @returndesc return code of @seeroutine ApplyBndScalar + -21 error reading boundary width array from table + -22 wrong size boundary width array in table @endreturndesc @@*/ -int BndScalar(const cGH *GH, int num_vars, int *vars, int *faces, int *tables) +int BndScalar(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 ApplyBndScalar */ - /*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 stencil in all directions */ + int dir; /* direction in which to apply bc */ CCTK_REAL scalar; - 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,28 +132,12 @@ int BndScalar(const cGH *GH, int num_vars, int *vars, int *faces, int *tables) CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, "Faces specification %d for Scalar boundary conditions on " "%s is not implemented yet. " - "Applying radiative bcs to all (external) faces.", faces[i], + "Applying Scalar bcs to all (external) faces.", faces[i], CCTK_VarName(vars[i])); } dir = 0; /* apply bc to all faces */ /* Set up default arguments for ApplyBndScalar */ - /* 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 */ scalar = 0.; /* Look on table for possible non-default arguments @@ -148,20 +152,58 @@ int BndScalar(const cGH *GH, int num_vars, int *vars, int *faces, int *tables) "Invalid table handle passed for Scalar boundary " "conditions for %s. Using all default values.", CCTK_VarName(vars[i])); + } + + /* 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 { - /* Stencil width array */ - Util_TableGetIntArray(tables[i], gdim, stencil_alldirs, "STENCIL WIDTH"); + for (k=0; k<2*gdim; ++k) + { + width_alldirs[k] = widths[i]; + } } - if ((retval = ApplyBndScalar(GH, 0, stencil_alldirs, dir, scalar, vars[i], + /* Apply the boundary condition */ + if ((retval = ApplyBndScalar(GH, 0, width_alldirs, dir, scalar, vars[i], j)) < 0) { CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, "ApplyBndScalar() returned %d", retval); } } - free(stencil_alldirs); + free(width_alldirs); return retval; } @@ -347,7 +389,7 @@ int BndScalarVI (const cGH *GH, if (vi >= 0 && vi < CCTK_NumVars ()) { - retval = ApplyBndScalar (GH, -1, stencil, 0, scalar, vi, 1); + retval = OldApplyBndScalar (GH, -1, stencil, 0, scalar, vi, 1); } else { @@ -498,7 +540,7 @@ int BndScalarGI (const cGH *GH, first_vi = CCTK_FirstVarIndexI (gi); if (first_vi >= 0) { - retval = ApplyBndScalar (GH, -1, stencil, 0, scalar, first_vi, + retval = OldApplyBndScalar (GH, -1, stencil, 0, scalar, first_vi, CCTK_NumVarsInGroupI (gi)); } else @@ -918,28 +960,28 @@ void CCTK_FCALL CCTK_FNAME (BndScalarVN) if (gdim > 0) \ { \ /* lower x */ \ - SCALAR_BOUNDARY_TYPED (doBC[0], stencil[0], lssh[1], lssh[2], \ + SCALAR_BOUNDARY_TYPED (doBC[0], widths[0], lssh[1], lssh[2], \ i, j, k, left_cctk_type, right_cctk_type); \ /* upper x */ \ - SCALAR_BOUNDARY_TYPED (doBC[1], stencil[0], lssh[1], lssh[2], \ + SCALAR_BOUNDARY_TYPED (doBC[1], widths[1], lssh[1], lssh[2], \ lssh[0]-i-1, j, k, left_cctk_type, right_cctk_type);\ } \ if (gdim > 1) \ { \ /* lower y */ \ - SCALAR_BOUNDARY_TYPED (doBC[2], lssh[0], stencil[1], lssh[2], \ + SCALAR_BOUNDARY_TYPED (doBC[2], lssh[0], widths[2], lssh[2], \ i, j, k, left_cctk_type, right_cctk_type); \ /* upper y */ \ - SCALAR_BOUNDARY_TYPED (doBC[3], lssh[0], stencil[1], lssh[2], \ + SCALAR_BOUNDARY_TYPED (doBC[3], lssh[0], widths[3], lssh[2], \ i, lssh[1]-j-1, k, left_cctk_type, right_cctk_type);\ } \ if (gdim > 2) \ { \ /* lower z */ \ - SCALAR_BOUNDARY_TYPED (doBC[4], lssh[0], lssh[1], stencil[2], \ + SCALAR_BOUNDARY_TYPED (doBC[4], lssh[0], lssh[1], widths[4], \ i, j, k, left_cctk_type, right_cctk_type); \ /* upper z */ \ - SCALAR_BOUNDARY_TYPED (doBC[5], lssh[0], lssh[1], stencil[2], \ + SCALAR_BOUNDARY_TYPED (doBC[5], lssh[0], lssh[1], widths[5], \ i, j, lssh[2]-k-1, left_cctk_type, right_cctk_type);\ } \ } @@ -967,13 +1009,13 @@ void CCTK_FCALL CCTK_FNAME (BndScalarVN) @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 @@ -1010,13 +1052,13 @@ void CCTK_FCALL CCTK_FNAME (BndScalarVN) 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 @endreturndesc @@*/ static int ApplyBndScalar (const cGH *GH, - int stencil_dir, - const int *stencil_alldirs, + int width_dir, + const int *in_widths, int dir, CCTK_REAL scalar, int first_var, @@ -1025,7 +1067,8 @@ static int ApplyBndScalar (const cGH *GH, int i, j, k; int gindex, gdim; int var, 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; /* check the direction parameter */ @@ -1050,19 +1093,20 @@ static int ApplyBndScalar (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, "ApplyBndScalar: NULL pointer passed " - "for stencil width array"); + "for boundary width array"); return (-3); } @@ -1176,3 +1220,95 @@ static int ApplyBndScalar (const cGH *GH, return(0); } + + +/*@@ + @routine OldApplyBndScalar + @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 boundary width in direction dir + @vtype int + @vio in + @endvar + @var stencil_alldirs + @vdesc boundary 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 scalar + @vdesc scalar value to set the boundaries to + @vtype CCTK_REAL + @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 + ApplyBndScalar + @returntype int + @returndesc + returncode from @seeroutine ApplyBndScalar + @endreturndesc +@@*/ + +static int OldApplyBndScalar(const cGH *GH, + int stencil_dir, + const int *stencil_alldirs, + int dir, + CCTK_REAL scalar, + 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 = ApplyBndScalar(GH, stencil_dir, boundary_widths, dir, scalar, + first_var, num_vars); + + return retval; +} |