diff options
Diffstat (limited to 'src/rotatingsymmetry180.c')
| -rw-r--r-- | src/rotatingsymmetry180.c | 146 |
1 files changed, 139 insertions, 7 deletions
diff --git a/src/rotatingsymmetry180.c b/src/rotatingsymmetry180.c index 7404770..fdb69e9 100644 --- a/src/rotatingsymmetry180.c +++ b/src/rotatingsymmetry180.c @@ -39,6 +39,7 @@ int BndRot180VI (cGH const * restrict const cctkGH, int fake_bbox[6]; CCTK_REAL x0[3], dx[3]; + CCTK_REAL symbnd[3]; /* location of symmetry boundary */ CCTK_REAL origin[3], dorigin[3]; int avoid_origin[3], iorigin[3]; int offset[3]; /* offset 0..1 due to avoid_origin */ @@ -274,11 +275,36 @@ int BndRot180VI (cGH const * restrict const cctkGH, alldirs[0] = 0; alldirs[1] = 1; + /* Find location of symmetry boundary */ + if (use_coordbase) { + CCTK_INT const size = 3; + CCTK_REAL physical_min[3]; + CCTK_REAL physical_max[3]; + CCTK_REAL interior_min[3]; + CCTK_REAL interior_max[3]; + CCTK_REAL exterior_min[3]; + CCTK_REAL exterior_max[3]; + CCTK_REAL spacing; + GetDomainSpecification + (size, + physical_min, physical_max, + interior_min, interior_max, + exterior_min, exterior_max, + & spacing); + symbnd[0] = physical_min[0]; + symbnd[1] = physical_min[1]; + symbnd[2] = physical_min[2]; + } else { + symbnd[0] = symmetry_boundary_x; + symbnd[1] = symmetry_boundary_y; + symbnd[2] = 0; /* unused */ + } + /* Find grid point that corresponds to the origin */ for (q=0; q<2; ++q) { d = alldirs[q]; - /* x0 + dx * origin == 0 */ - origin[d] = - x0[d] / dx[d]; + /* x0 + dx * origin == symbnd */ + origin[d] = (symbnd[d] - x0[d]) / dx[d]; dorigin[d] = origin[d] - floor(origin[d]); if (fabs(dorigin[d]) < 1.0e-6 || fabs(dorigin[d] - 1.0) < 1.0e-6) { avoid_origin[d] = 0; @@ -307,6 +333,15 @@ int BndRot180VI (cGH const * restrict const cctkGH, "xyz"[otherdir]); } + if (iorigin[dir] != data.nghostzones[dir]) { + CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, + "The group \"%s\" has in the %c-direction %d symmetry points (grid points outside of the symmetry boundary). " + "This is not equal to the number of ghost zones, which is %d. " + "The number of symmetry points must be equal to the number of ghost zones.", + CCTK_GroupNameFromVarI(vis[var]), "xyz"[dir], + iorigin[dir], data.nghostzones[dir]); + } + if (data.gsh[dir] < 2*data.nghostzones[dir] + offset[dir]) { assert (nvars > 0); CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, @@ -318,6 +353,47 @@ int BndRot180VI (cGH const * restrict const cctkGH, 2*data.nghostzones[dir] + offset[dir]); } + if (poison_boundaries) { + /* poison destination grid points */ + if (cctkGH->cctk_bbox[2*dir]) { + + int imin[3], imax[3]; + int i, j, k; + for (d=0; d<3; ++d) { + imin[d] = xferinfo[d].dst.off; + imax[d] = xferinfo[d].dst.off + xferinfo[d].dst.len; + imin[d] -= cctk_lbnd[d]; + imax[d] -= cctk_lbnd[d]; + if (imin[d] < 0) imin[d] = 0; + if (imax[d] >= cctk_lsh[d]) imax[d] = cctk_lsh[d]; + } + + var = 0; + + assert (group.dim == 3); + switch (group.vartype) { + case CCTK_VARIABLE_INT: + /* do nothing */ + break; + case CCTK_VARIABLE_REAL: + for (k=imin[0]; k<imax[2]; ++k) { + for (j=imin[1]; j<imax[1]; ++j) { + for (i=imin[2]; i<imax[0]; ++i) { + const int ind = CCTK_GFINDEX3D(cctkGH, i, j, k); + ((CCTK_REAL *) varptrs[var]) [ind] = poison_value; + } + } + } + break; + case CCTK_VARIABLE_COMPLEX: + /* do nothing */ + break; + default: + assert (0); + } /* switch grouptype */ + } /* if bbox */ + } /* if poison_boundaries */ + /* Allocate slab transfer description */ xferinfo = malloc(group.dim * sizeof *xferinfo); assert (xferinfo); @@ -361,6 +437,48 @@ int BndRot180VI (cGH const * restrict const cctkGH, nvars, vartypes, varptrs, vartypes, varptrs); assert (!ierr); + if (poison_boundaries) { + /* check destination grid points for poison */ + if (cctkGH->cctk_bbox[2*dir]) { + + int imin[3], imax[3]; + int i, j, k; + for (d=0; d<3; ++d) { + imin[d] = xferinfo[d].dst.off; + imax[d] = xferinfo[d].dst.off + xferinfo[d].dst.len; + imin[d] -= cctk_lbnd[d]; + imax[d] -= cctk_lbnd[d]; + if (imin[d] < 0) imin[d] = 0; + if (imax[d] >= cctk_lsh[d]) imax[d] = cctk_lsh[d]; + } + + var = 0; + + assert (group.dim == 3); + switch (group.vartype) { + case CCTK_VARIABLE_INT: + /* do nothing */ + break; + case CCTK_VARIABLE_REAL: + for (k=imin[0]; k<imax[2]; ++k) { + for (j=imin[1]; j<imax[1]; ++j) { + for (i=imin[2]; i<imax[0]; ++i) { + const int ind = CCTK_GFINDEX3D(cctkGH, i, j, k); + CCTK_REAL const val = ((CCTK_REAL const *) varptrs[var]) [ind]; + assert (val != poison_value); + } + } + } + break; + case CCTK_VARIABLE_COMPLEX: + /* do nothing */ + break; + default: + assert (0); + } /* switch grouptype */ + } /* if bbox */ + } /* if poison_boundaries */ + /* take parity into account */ if (cctkGH->cctk_bbox[2*dir]) { for (var=0; var<nvars; ++var) { @@ -401,17 +519,31 @@ int BndRot180VI (cGH const * restrict const cctkGH, } } break; + case CCTK_VARIABLE_COMPLEX: + { + CCTK_COMPLEX const czero = CCTK_Cmplx (0.0, 0.0); + for (k=imin[0]; k<imax[2]; ++k) { + for (j=imin[1]; j<imax[1]; ++j) { + for (i=imin[2]; i<imax[0]; ++i) { + const int ind = CCTK_GFINDEX3D(cctkGH, i, j, k); + CCTK_COMPLEX * restrict const ptr = + & ((CCTK_COMPLEX *) varptrs[var]) [ind]; + * ptr = CCTK_CmplxSub (czero, * ptr); + } + } + } + } + break; default: assert (0); } /* switch grouptype */ } } /* for var */ - - free (xferinfo); - - } /* if global_bbox */ + } /* if bbox */ + + free (xferinfo); - } /* if bbox */ + } /* if global_bbox */ free (gis); free (varptrs); |