CarpetReduce: Correct handling staggered boundaries

1 files changed, 99 insertions, 104 deletions

diff --git a/Carpet/CarpetReduce/src/mask_coords.c b/Carpet/CarpetReduce/src/mask_coords.c
index b018dbce4..271e1b957 100644
--- a/Carpet/CarpetReduce/src/mask_coords.c
+++ b/Carpet/CarpetReduce/src/mask_coords.c
@@ -16,17 +16,17 @@ CoordBase_SetupMask (CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
-  CCTK_INT nboundaryzones[6];
-  CCTK_INT is_internal[6];
-  CCTK_INT is_staggered[6];
-  CCTK_INT shiftout[6];
+  CCTK_INT nboundaryzones[2*cctk_dim];
+  CCTK_INT is_internal[2*cctk_dim];
+  CCTK_INT is_staggered[2*cctk_dim];
+  CCTK_INT shiftout[2*cctk_dim];
   
-  int bnd_points[6];            /* points outside the domain */
-  int int_points[3];            /* global interior points */
+  int bnd_points[2*cctk_dim];   /* points outside the domain */
+  int int_points[cctk_dim];     /* global interior points */
   
-  int gmin[3], gmax[3];         /* global domain extent */
-  int imin[3], imax[3];         /* domain extent */
-  int bmin[3], bmax[3];         /* boundary extent */
+  int gmin[cctk_dim], gmax[cctk_dim]; /* global domain extent */
+  int imin[cctk_dim], imax[cctk_dim]; /* domain extent */
+  int bmin[cctk_dim], bmax[cctk_dim]; /* boundary extent */
   
   int ierr;
   
@@ -40,24 +40,25 @@ CoordBase_SetupMask (CCTK_ARGUMENTS)
     int const m = MultiPatch_GetMap (cctkGH);
     assert (m >= 0);
     ierr = MultiPatch_GetBoundarySpecification
-      (m, 6, nboundaryzones, is_internal, is_staggered, shiftout);
+      (m, 2*cctk_dim, nboundaryzones, is_internal, is_staggered, shiftout);
   } else {
     ierr = GetBoundarySpecification
-      (6, nboundaryzones, is_internal, is_staggered, shiftout);
+      (2*cctk_dim, nboundaryzones, is_internal, is_staggered, shiftout);
   }
   if (ierr != 0) {
-    CCTK_WARN (0, "Could not get boundary specification");
+    CCTK_WARN (CCTK_WARN_ABORT, "Could not get boundary specification");
   }
   
   
   
-  /* Calculate the number of boundary points */
-  for (int d=0; d<3; ++d) {
+  /* Calculate the number of boundary points. This excludes points
+     that are directly on the boundary. */
+  for (int d=0; d<cctk_dim; ++d) {
     for (int f=0; f<2; ++f) {
       
       if (is_internal[2*d+f]) {
         /* The boundary extends inwards */
-        bnd_points[2*d+f] = shiftout[2*d+f] - is_staggered[2*d+f];
+        bnd_points[2*d+f] = shiftout[2*d+f];
       } else {
         /* The boundary extends outwards */
         bnd_points[2*d+f] =
@@ -68,126 +69,120 @@ CoordBase_SetupMask (CCTK_ARGUMENTS)
   }
   
   /* Calculate the global extent of the domain */
-  for (int d=0; d<3; ++d) {
+  for (int d=0; d<cctk_dim; ++d) {
     gmin[d] = 0;
     gmax[d] = cctk_gsh[d];
   }
   
   /* Ensure that the domain specification is consistent */
-  for (int d=0; d<3; ++d) {
-    int_points[d] =
-      (gmax[d] - bnd_points[2*d]) - (gmin[d] + bnd_points[2*d+1]);
+  for (int d=0; d<cctk_dim; ++d) {
+    int_points[d] = (gmax[d] - bnd_points[2*d+1]) - (gmin[d] + bnd_points[2*d]);
     if (int_points[d] < 0) {
-      CCTK_WARN (0, "Number of internal grid points is negative");
+      CCTK_WARN (CCTK_WARN_ABORT, "Number of internal grid points is negative");
     }
   }
   
   
   
   /* Loop over all dimensions and faces */
-  for (int d=0; d<3; ++d) {
+  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]) {
         
-        /* If there are boundary that should be ignored */
-        if (bnd_points[2*d+f] >= 0) {
-          
-          if (bnd_points[2*d+f] < 0 || bnd_points[2*d+f] > cctk_lsh[d]) {
-            CCTK_WARN (0, "Illegal number of boundary points");
-          }
-          
-          /* Calculate the extent of the domain */
-          for (int dd=0; dd<3; ++dd) {
-            imin[dd] = 0;
-            imax[dd] = cctk_lsh[dd];
-          }
-          
-          /* Calculate the extent of the boundary */
-          for (int dd=0; dd<3; ++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 boundary points in direction %d face %d to weight 0 on level %d", d, f, reflevel);
-          }
+        if (bnd_points[2*d+f] < 0 || bnd_points[2*d+f] > cctk_lsh[d]) {
+          CCTK_WARN (CCTK_WARN_ABORT, "Illegal number of boundary points");
+        }
+        
+        /* 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 boundary points in direction %d face %d to weight 0 on level %d", d, f, reflevel);
+        }
 #pragma omp parallel
-          LC_LOOP3(CoordBase_SetupMask_boundary,
+        CCTK_LOOP3(CoordBase_SetupMask_boundary,
                    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_boundary);
+        
+        /* When the boundary is not staggered, then give the points
+           directly on the boundary the weight 1/2 */
+        if (! is_staggered[2*d+f]) {
+          
+          /* Check whether the domain is empty */
+          if (int_points[d] == 0) {
             
-            int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
-            iweight[ind] = 0;
+            /* The domain is empty. The correct thing to do would be
+               to set the weights to 0. But this is boring, because
+               then there are no interior points left, and all
+               reductions become trivial. Instead, leave the
+               non-staggered boundary points alone, and assume that
+               the user wants to perform a reduction in one fewer
+               dimensions. */
+            /* do nothing */
             
-          } LC_ENDLOOP3(CoordBase_SetupMask);
-          
-          /* When the boundary is not staggered, then give the points
-             on the boundary the weight 1/2 */
-          if (! is_staggered[2*d+f]) {
+          } else {
             
-            /* Check whether the domain is empty */
-            if (int_points[d] == 0) {
-              
-              /* The domain is empty.  The correct thing to do would
-                 be to set the weights to 0.  But this is boring,
-                 because then there are no interior points left, and
-                 all reductions become trivial.  Instead, leave the
-                 non-staggered boundary points alone, and assume that
-                 the user wants to perform a reduction in one fewer
-                 dimensions.  */
-              /* do nothing */
-              
+            /* Adapt the extent of the boundary */
+            if (f==0) {
+              /* lower face */
+              bmin[d] = bmax[d];
+              bmax[d] = bmin[d] + 1;
             } else {
-              
-              /* Adapt the extent of the boundary */
-              if (f==0) {
-                /* lower face */
-                bmin[d] = bmax[d];
-                bmax[d] = bmin[d] + 1;
-              } else {
-                /* upper face */
-                bmax[d] = bmin[d];
-                bmin[d] = bmax[d] - 1;
-              }
-              
-              /* Loop over the points next to boundary */
-              if (verbose) {
-                CCTK_VInfo (CCTK_THORNSTRING,
-                            "Setting non-staggered boundary points in direction %d face %d to weight 1/2 on level %d", d, f, reflevel);
-              }
-              int bmask = 0;
-              for (unsigned b=0; b<BMSK(cctk_dim); ++b) {
-                if (BGET(b,d)==f) {
-                  bmask = BSET(bmask, b);
-                }
+              /* upper face */
+              bmax[d] = bmin[d];
+              bmin[d] = bmax[d] - 1;
+            }
+            
+            /* Loop over the points next to boundary */
+            if (verbose) {
+              CCTK_VInfo (CCTK_THORNSTRING,
+                          "Setting non-staggered boundary points in direction %d face %d to weight 1/2 on level %d", d, f, reflevel);
+            }
+            unsigned const bits = BMSK(cctk_dim);
+            unsigned bmask = 0;
+            for (unsigned b=0; b<bits; ++b) {
+              if (BGET(b,d) == f) {
+                bmask = BSET(bmask, b);
               }
+            }
+            assert (BCNT(bmask) == bits/2);
 #pragma omp parallel
-              LC_LOOP3(CoordBase_SetupMask_boundary2,
+            CCTK_LOOP3(CoordBase_SetupMask_boundary2,
                        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] &= ~bmask;
-              } LC_ENDLOOP3(CoordBase_SetupMask_boundary2);
-              
-            } /* if the domain is not empty */
+            {
+              int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
+              iweight[ind] &= ~bmask;
+            } CCTK_ENDLOOP3(CoordBase_SetupMask_boundary2);
             
-          } /* if the boundary is not staggered */
+          } /* if the domain is not empty */
           
-        } /* if there are boundary points */
+        } /* if the boundary is not staggered */
         
       } /* if is outer boundary */
     } /* loop over faces */