Add parameter to poison periodic boundaries

Add a parameter to poison periodic boundaries before the periodicity
condition is applied. This allows detecting errors in the periodicity
boundary condition, which does not work if there are multiple local
components per process.


git-svn-id: http://svn.cactuscode.org/arrangements/CactusNumerical/Periodic/trunk@31 1bf05452-ddb3-4880-bfa1-00436340132b

2 files changed, 117 insertions, 6 deletions

diff --git a/param.ccl b/param.ccl
index fa1f5ca..d43e02b 100644
--- a/param.ccl
+++ b/param.ccl
@@ -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;
 }