Handle empty groups.

Reduce cctk_bbox to find out whether there are outer boundaries.

Fixes / changes to handling domains smaller than the stencil size.


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

1 files changed, 81 insertions, 73 deletions

diff --git a/src/periodic.c b/src/periodic.c
index 0f65ff0..b114919 100644
--- a/src/periodic.c
+++ b/src/periodic.c
@@ -26,6 +26,7 @@ BndPeriodicVI (cGH const * restrict const cctkGH,
   void * restrict varptr;
   struct xferinfo * restrict xferinfo;
   int do_periodic[3];
+  int global_bbox[6];
   int gi;
   int dir;
   int d, f;
@@ -66,6 +67,17 @@ BndPeriodicVI (cGH const * restrict const cctkGH,
   do_periodic[1] = periodic || periodic_y;
   do_periodic[2] = periodic || periodic_z;
   
+  {
+    int max_handle;
+    CCTK_REAL local[6], global[6];
+    max_handle = CCTK_ReductionArrayHandle ("maximum");
+    if (max_handle<0) CCTK_WARN (0, "Could not obtain reduction handle");
+    for (d=0; d<6; ++d) local[d] = cctkGH->cctk_bbox[d];
+    ierr = CCTK_ReduceLocArrayToArray1D
+      (cctkGH, -1, max_handle, local, global, 6, CCTK_VARIABLE_REAL);
+    for (d=0; d<6; ++d) global_bbox[d] = (int)global[d];
+  }
+  
   /* Allocate slab transfer description */
   xferinfo = malloc(group.dim * sizeof *xferinfo);
   assert (xferinfo);
@@ -88,70 +100,67 @@ BndPeriodicVI (cGH const * restrict const cctkGH,
       /* Loop over faces */
       for (f=0; f<2; ++f) {
         
-        /* Fill in slab transfer description */
-        for (d=0; d<group.dim; ++d) {
-          xferinfo[d].src.gsh         = data.gsh        [d];
-          xferinfo[d].src.lbnd        = data.lbnd       [d];
-          xferinfo[d].src.lsh         = data.lsh        [d];
-          xferinfo[d].src.lbbox       = data.bbox       [2*d  ];
-          xferinfo[d].src.ubbox       = data.bbox       [2*d+1];
-          xferinfo[d].src.nghostzones = data.nghostzones[d];
-          xferinfo[d].src.off         = 0;
-          xferinfo[d].src.str         = 1;
-          xferinfo[d].src.len         = data.gsh        [d];
-          
-          xferinfo[d].dst.gsh         = data.gsh        [d];
-          xferinfo[d].dst.lbnd        = data.lbnd       [d];
-          xferinfo[d].dst.lsh         = data.lsh        [d];
-          xferinfo[d].dst.lbbox       = data.bbox       [2*d  ];
-          xferinfo[d].dst.ubbox       = data.bbox       [2*d+1];
-          xferinfo[d].dst.nghostzones = data.nghostzones[d];
-          xferinfo[d].dst.off         = 0;
-          xferinfo[d].dst.str         = 1;
-          xferinfo[d].dst.len         = data.gsh        [d];
-          
-          xferinfo[d].xpose = d;
-          xferinfo[d].flip  = 0;
-        }
-        
-        {
-          int const domain_size = data.gsh[dir] - 2*stencil[dir];
-          int step_size = stencil[dir];
-          int num_steps = 1;
-          int step;
-          int source;
+        if (global_bbox[2*dir+f]) {
           
-          assert (domain_size > 0);
-          if (f==0 && domain_size < step_size) {
-            /* TODO: this could be made more efficient by taking
-               larger steps */
-            step_size = 1;
-            num_steps = stencil[dir];
+          /* Fill in slab transfer description */
+          for (d=0; d<group.dim; ++d) {
+            xferinfo[d].src.gsh         = data.gsh        [d];
+            xferinfo[d].src.lbnd        = data.lbnd       [d];
+            xferinfo[d].src.lsh         = data.lsh        [d];
+            xferinfo[d].src.lbbox       = data.bbox       [2*d  ];
+            xferinfo[d].src.ubbox       = data.bbox       [2*d+1];
+            xferinfo[d].src.nghostzones = data.nghostzones[d];
+            xferinfo[d].src.off         = 0;
+            xferinfo[d].src.str         = 1;
+            xferinfo[d].src.len         = data.gsh        [d];
+            
+            xferinfo[d].dst.gsh         = data.gsh        [d];
+            xferinfo[d].dst.lbnd        = data.lbnd       [d];
+            xferinfo[d].dst.lsh         = data.lsh        [d];
+            xferinfo[d].dst.lbbox       = data.bbox       [2*d  ];
+            xferinfo[d].dst.ubbox       = data.bbox       [2*d+1];
+            xferinfo[d].dst.nghostzones = data.nghostzones[d];
+            xferinfo[d].dst.off         = 0;
+            xferinfo[d].dst.str         = 1;
+            xferinfo[d].dst.len         = data.gsh        [d];
+            
+            xferinfo[d].xpose = d;
+            xferinfo[d].flip  = 0;
           }
           
-          for (step=0; step<num_steps; ++step) {
+          {
+            int const domain_size = data.gsh[dir] - 2*stencil[dir];
+            int step_size = stencil[dir];
+            int num_steps = 1;
+            int step;
+            int source;
             
-            if (cctkGH->cctk_bbox[2*dir+f]) {
+            assert (domain_size > 0);
+            if (f==0 && domain_size < step_size) {
+              /* TODO: this could be made more efficient by taking
+                 larger steps */
+              step_size = 1;
+              num_steps = stencil[dir];
+            }
+            assert (num_steps*step_size == stencil[dir]);
+            
+            for (step=0; step<num_steps; ++step) {
+              
               if (f==0) {
                 /* Fill in lower face */
 #if 0
-                xferinfo[dir].src.off = data.gsh[dir] - 2 * stencil[dir];
+                xferinfo[dir].src.off = data.gsh[dir] - 2*stencil[dir];
                 xferinfo[dir].src.len = stencil[dir];
                 xferinfo[dir].dst.len = stencil[dir];
-#endif
-                source = data.gsh[dir] - 2 * stencil[dir] + step;
-                while (source < stencil[dir]) {
-                  source += domain_size;
-                }
-                while (source >= data.gsh[dir]-stencil[dir]) {
-                  source -= domain_size;
-                }
-                assert (source >= stencil[dir]);
-                assert (source < data.gsh[dir]-stencil[dir]);
+#else
+                source = data.gsh[dir] - 2*stencil[dir] + step*step_size;
+                while (source < stencil[dir]) source += domain_size;
+                assert (source+step_size <= data.gsh[dir]-stencil[dir]);
                 xferinfo[dir].src.off = source;
                 xferinfo[dir].src.len = step_size;
-                xferinfo[dir].dst.off = step;
+                xferinfo[dir].dst.off = step*step_size;
                 xferinfo[dir].dst.len = step_size;
+#endif
               } else {
                 /* Fill in upper face */
                 xferinfo[dir].src.off = stencil[dir];
@@ -159,18 +168,16 @@ BndPeriodicVI (cGH const * restrict const cctkGH,
                 xferinfo[dir].dst.off = data.gsh[dir] - stencil[dir];
                 xferinfo[dir].dst.len = stencil[dir];
               }
-            } else {
-              xferinfo[dir].src.len = 0;
-              xferinfo[dir].dst.len = 0;
-            }
-            
-            ierr = Slab_Transfer
-              (cctkGH, group.dim, xferinfo, -1,
-               group.vartype, varptr, group.vartype, varptr);
-            assert (!ierr);
-            
-          } /* for step */
-        }
+              
+              ierr = Slab_Transfer
+                (cctkGH, group.dim, xferinfo, -1,
+                 group.vartype, varptr, group.vartype, varptr);
+              assert (!ierr);
+              
+            } /* for step */
+          }
+          
+        } /* if bbox */
         
       } /* for f */
       
@@ -201,18 +208,19 @@ BndPeriodicGI (cGH const * restrict const cctkGH,
                int const * restrict const stencil,
                int const gi)
 {
-  int vi1, vi2;
+  int v1, nv;
   int vi;
   int ierr;
   assert (gi>=0 && gi<CCTK_NumGroups());
-  vi1 = CCTK_FirstVarIndexI(gi);
-  assert (vi1>=0 && vi1<CCTK_NumVars());
-  vi2 = vi1 + CCTK_NumVarsInGroupI(gi);
-  assert (vi2>=vi1 && vi2<CCTK_NumVars());
-  assert (vi>=0 && vi<CCTK_NumVars());
-  for (vi=vi1; vi<vi2; ++vi) {
-    ierr = BndPeriodicVI (cctkGH, stencil, vi);
-    if (ierr) return ierr;
+  nv = CCTK_NumVarsInGroup(gi);
+  assert (nv>=0);
+  if (nv>0) {
+    v1 = CCTK_FirstVarIndexI(gi);
+    assert (v1>=0 && v1<CCTK_NumVars());
+    for (vi=v1; vi<v1+nv; ++vi) {
+      ierr = BndPeriodicVI (cctkGH, stencil, vi);
+      if (ierr) return ierr;
+    }
   }
   return 0;
 }