CarpetReduce: Correct weight calculation for cell-centred grid hierarchies

2 files changed, 111 insertions, 33 deletions

diff --git a/Carpet/CarpetReduce/src/mask_carpet.cc b/Carpet/CarpetReduce/src/mask_carpet.cc
index 81c53e7b7..160471cb4 100644
--- a/Carpet/CarpetReduce/src/mask_carpet.cc
+++ b/Carpet/CarpetReduce/src/mask_carpet.cc
@@ -16,6 +16,7 @@
 
 
 
+#if 0
 #define LOOP_OVER_NEIGHBOURS(dir)               \
 {                                               \
   ivect dir_(-1);                               \
@@ -33,6 +34,7 @@
     }                                           \
   } while (not all (dir_ == -1));               \
 }
+#endif
 
 
 
@@ -112,7 +114,15 @@ namespace CarpetMask {
       ibset test_boxes;
       ibset test_cfboxes;
       
-      LOOP_OVER_NEIGHBOURS (shift) {
+      for (int neighbour=0; neighbour<ipow(3,dim); ++neighbour) {
+        ivect shift;
+        int itmp=neighbour;
+        for (int d=0; d<dim; ++d) {
+          shift[d] = itmp % 3 - 1; // [-1 ... +1]
+          itmp /= 3;
+        }
+        assert (itmp==0);
+        
         // In this loop, shift [1,1,1] denotes a convex corner of the
         // region which should be masked out, i.e. a region where only
         // a small bit (1/8) of the region should be masked out.
@@ -122,34 +132,51 @@ namespace CarpetMask {
         
         ibset boxes  = not_active;
         ibset fboxes = fine_active;
-        for (int d=0; d<dim; ++d) {
-          ivect const dir = ivect::dir(d);
-          fboxes = fboxes.shift(-dir) & fboxes & fboxes.shift(+dir);
-        }
-        for (int d=0; d<dim; ++d) {
-          // Calculate the boundary in direction d
-          ivect const dir = ivect::dir(d);
-          switch (shift[d]) {
-          case -1: {
-            // left boundary
-            boxes  = boxes.shift (-dir) - boxes;
-            fboxes = fboxes.shift(-dir) - fboxes;
-            break;
+        switch (hh.refcent) {
+        case vertex_centered: {
+          for (int d=0; d<dim; ++d) {
+            ivect const dir = ivect::dir(d);
+            fboxes = fboxes.shift(-dir) & fboxes & fboxes.shift(+dir);
           }
-          case 0: {
-            // interior
-            // do nothing
-            break;
-          }
-          case +1: {
-            // right boundary
-            boxes  = boxes.shift (+dir) - boxes;
-            fboxes = fboxes.shift(+dir) - fboxes;
-            break;
+          for (int d=0; d<dim; ++d) {
+            // Calculate the boundary in direction d
+            ivect const dir = ivect::dir(d);
+            switch (shift[d]) {
+            case -1: {
+              // left boundary
+              boxes  = boxes.shift (-dir) - boxes;
+              fboxes = fboxes.shift(-dir) - fboxes;
+              break;
+            }
+            case 0: {
+              // interior
+              // do nothing
+              break;
+            }
+            case +1: {
+              // right boundary
+              boxes  = boxes.shift (+dir) - boxes;
+              fboxes = fboxes.shift(+dir) - fboxes;
+              break;
+            }
+            default:
+              assert (0);
+            }
           }
-          default:
-            assert (0);
+          break;
+        }
+        case cell_centered: {
+          // Assume that all cell boundaries are aligned
+          if (all(shift == 0)) {
+            // do nothing
+          } else {
+            boxes  = ibset();
+            fboxes = ibset();
           }
+          break;
+        }
+        default:
+          assert (0);
         }
         boxes &= ext;
         ibset const cfboxes = fboxes.contracted_for(base) & ext;
@@ -213,16 +240,29 @@ namespace CarpetMask {
           } CCTK_ENDLOOP3(CarpetMaskSetup_restriction);
         } END_LOOP_OVER_BSET;
         
-      } END_LOOP_OVER_NEIGHBOURS;
+      } // for neighbours
       
       {
-        ibset const boxes = not_active.expand(ivect(1), ivect(1)) & ext;
+        ibset boxes  = not_active;
         ibset fboxes = fine_active;
-        for (int d=0; d<dim; ++d) {
-          ivect const dir = ivect::dir(d);
-          fboxes = fboxes.shift(-dir) & fboxes & fboxes.shift(+dir);
+        switch (hh.refcent) {
+        case vertex_centered: {
+          boxes = boxes.expand(ivect(1), ivect(1));
+          for (int d=0; d<dim; ++d) {
+            ivect const dir = ivect::dir(d);
+            fboxes = fboxes.shift(-dir) & fboxes & fboxes.shift(+dir);
+          }
+          fboxes = fboxes.expand(ivect(1), ivect(1));
+          break;
+        }
+        case cell_centered: {
+          // do nothing
+          break;
         }
-        fboxes = fboxes.expand(ivect(1), ivect(1));
+        default:
+          assert (0);
+        }
+        boxes &= ext;
         ibset const cfboxes = fboxes.contracted_for(base) & ext;
         if (not (test_boxes   == boxes  ) or
             not (test_cfboxes == cfboxes))
diff --git a/Carpet/CarpetReduce/src/mask_coords.c b/Carpet/CarpetReduce/src/mask_coords.c
index 271e1b957..d861a6d2b 100644
--- a/Carpet/CarpetReduce/src/mask_coords.c
+++ b/Carpet/CarpetReduce/src/mask_coords.c
@@ -87,6 +87,7 @@ CoordBase_SetupMask (CCTK_ARGUMENTS)
   /* Loop over all dimensions and faces */
   for (int d=0; d<cctk_dim; ++d) {
     for (int f=0; f<2; ++f) {
+      
       /* If this processor has the outer boundary */
       if (cctk_bbox[2*d+f]) {
         
@@ -184,7 +185,44 @@ CoordBase_SetupMask (CCTK_ARGUMENTS)
           
         } /* if the boundary is not staggered */
         
-      } /* if is outer boundary */
+      } else { /* if this is a ghost boundary */
+        
+        /* Calculate the extent of the local part of the domain */
+        for (int dd=0; dd<cctk_dim; ++dd) {
+          imin[dd] = 0;
+          imax[dd] = cctk_lsh[dd];
+        }
+        
+        /* Calculate the extent of the boundary */
+        for (int dd=0; dd<cctk_dim; ++dd) {
+          bmin[dd] = imin[dd];
+          bmax[dd] = imax[dd];
+        }
+        if (f==0) {
+          /* lower face */
+          bmax[d] = imin[d] + bnd_points[2*d+f];
+        } else {
+          /* upper face */
+          bmin[d] = imax[d] - bnd_points[2*d+f];
+        }
+        
+        /* Loop over the boundary */
+        if (verbose) {
+          CCTK_VInfo (CCTK_THORNSTRING,
+                      "Setting ghost points in direction %d face %d to weight 0 on level %d", d, f, reflevel);
+        }
+#pragma omp parallel
+        CCTK_LOOP3(CoordBase_SetupMask_ghost,
+                   i,j,k,
+                   bmin[0],bmin[1],bmin[2], bmax[0],bmax[1],bmax[2],
+                   cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
+        {
+          int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
+          iweight[ind] = 0;
+        } CCTK_ENDLOOP3(CoordBase_SetupMask_ghost);
+        
+      } /* if this is a ghost boundary */
+      
     } /* loop over faces */
   } /* loop over directions */