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

1 files changed, 106 insertions, 6 deletions

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;
 }