diff options
author | tradke <tradke@6a38eb6e-646e-4a02-a296-d141613ad6c4> | 2001-06-12 10:42:07 +0000 |
---|---|---|
committer | tradke <tradke@6a38eb6e-646e-4a02-a296-d141613ad6c4> | 2001-06-12 10:42:07 +0000 |
commit | b572512daf1b2fea0c9251a4f47924f5463b8d0d (patch) | |
tree | c087c289e9aca928827979486eb5799684a84ca6 | |
parent | 30c6cb6dc7c54b0b9af2a1d18e2759af97fa70e0 (diff) |
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
-rw-r--r-- | src/RobinBoundary.c | 269 |
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); } } |