diff options
| -rw-r--r-- | param.ccl | 11 | ||||
| -rw-r--r-- | src/periodic.c | 112 |
2 files changed, 117 insertions, 6 deletions
@@ -19,3 +19,14 @@ BOOLEAN periodic_y "Periodic boundary conditions in y-direction" BOOLEAN periodic_z "Periodic boundary conditions in z-direction" { } "no" + + + +BOOLEAN poison_boundaries "Fill the symmetry boundary with a poison value before the symmetry is applied, and check afterwards whether it has been overwritten" STEERABLE=always +{ +} "no" + +CCTK_INT poison_value "Integer value (0..255) used to poison new timelevels (with memset)" STEERABLE=always +{ + 0:255 :: "Must fit into a byte. Use 0 for zero, 255 for nan, and e.g. 113 for a large value." +} 254 diff --git a/src/periodic.c b/src/periodic.c index 9f672df..e347e86 100644 --- a/src/periodic.c +++ b/src/periodic.c @@ -1,6 +1,7 @@ #include <assert.h> #include <stdlib.h> #include <stdio.h> +#include <string.h> #include "cctk.h" #include "cctk_Parameters.h" #include "Slab.h" @@ -15,7 +16,8 @@ BndPeriodicVI (CCTK_POINTER_TO_CONST _GH, int const vi) { cGH const * restrict const cctkGH = _GH; - + DECLARE_CCTK_PARAMETERS; + cGroup group; cGroupDynamicData data; void * restrict varptr; @@ -24,8 +26,8 @@ BndPeriodicVI (CCTK_POINTER_TO_CONST _GH, int global_lbnd[3], global_ubnd[3]; int fake_bbox[6]; int gi; - int dir; - int d, f; + int dir, face; + int d; int ierr; /* Check arguments */ @@ -42,6 +44,8 @@ BndPeriodicVI (CCTK_POINTER_TO_CONST _GH, assert (group.grouptype == CCTK_GF); assert (group.disttype == CCTK_DISTRIB_DEFAULT); assert (group.stagtype == 0); + int const vartypesize = CCTK_VarTypeSize(group.vartype); + assert (vartypesize>0); if(group.dim != size) { CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, @@ -84,6 +88,46 @@ BndPeriodicVI (CCTK_POINTER_TO_CONST _GH, } } + if (poison_boundaries) { + /* poison destination grid points */ + + for (dir=0; dir<group.dim; ++dir) { + if (dir<3 && do_periodic[dir]) { + assert (stencil[dir] >= 0); + for (face=0; face<2; ++face) { + if (cctkGH->cctk_bbox[2*dir+face]) { + + int imin[3], imax[3]; + for (d=0; d<group.dim; ++d) { + imin[d] = 0; + imax[d] = data.lsh[d]; + } + for (d=group.dim; d<3; ++d) { + imin[d] = 0; + imax[d] = 1; + } + if (face==0) { + imax[dir] = imin[dir] + stencil[dir]; + } else { + imin[dir] = imax[dir] - stencil[dir]; + } + +#pragma omp parallel for + for (int k=imin[2]; k<imax[2]; ++k) { + for (int j=imin[1]; j<imax[1]; ++j) { + int const ind3d = imin[0] + data.lsh[0] * (j + data.lsh[1] * k); + memset ((char*)varptr + ind3d*vartypesize, + poison_value, (imax[0]-imin[0])*vartypesize); + } + } + + } + } + } + } + + } + /* Allocate slab transfer description */ xferinfo = malloc(group.dim * sizeof *xferinfo); assert (xferinfo); @@ -98,9 +142,9 @@ BndPeriodicVI (CCTK_POINTER_TO_CONST _GH, } /* Loop over faces */ - for (f=0; f<2; ++f) { + for (face=0; face<2; ++face) { - if (global_bbox[2*dir+f]) { + if (global_bbox[2*dir+face]) { /* Fill in slab transfer description */ for (d=0; d<group.dim; ++d) { @@ -146,7 +190,7 @@ BndPeriodicVI (CCTK_POINTER_TO_CONST _GH, for (step=0; step<num_steps; ++step) { - if (f==0) { + if (face==0) { /* Fill in lower face */ source = data.gsh[dir] - 2*stencil[dir] + step*step_size; while (source < stencil[dir]) { @@ -188,6 +232,62 @@ BndPeriodicVI (CCTK_POINTER_TO_CONST _GH, free (xferinfo); + if (poison_boundaries) { + /* check destination grid points for poison */ + + char poison[vartypesize]; + memset (poison, poison_value, vartypesize); + + for (dir=0; dir<group.dim; ++dir) { + if (dir<3 && do_periodic[dir]) { + assert (stencil[dir] >= 0); + for (face=0; face<2; ++face) { + if (cctkGH->cctk_bbox[2*dir+face]) { + + int imin[3], imax[3]; + for (d=0; d<group.dim; ++d) { + imin[d] = 0; + imax[d] = data.lsh[d]; + } + for (d=group.dim; d<3; ++d) { + imin[d] = 0; + imax[d] = 1; + } + if (face==0) { + imax[dir] = imin[dir] + stencil[dir]; + } else { + imin[dir] = imax[dir] - stencil[dir]; + } + + int nerrors = 0; +#pragma omp parallel for reduction(+: nerrors) + for (int k=imin[2]; k<imax[2]; ++k) { + for (int j=imin[1]; j<imax[1]; ++j) { + for (int i=imin[0]; i<imax[0]; ++i) { + int const ind3d = i + data.lsh[0] * (j + data.lsh[1] * k); + if (! memcmp((char*)varptr + ind3d*vartypesize, poison, + vartypesize)) + { + ++nerrors; + } + } + } + } + if (nerrors > 0) { + char *const fullname = CCTK_FullName(vi); + CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, + "Found poison on symmetry boundary of variable \"%s\" direction %d face %d after applying periodicity boundary condition", + fullname, dir, face); + free (fullname); + } + + } + } + } + } + + } + return 0; } |