Merge branch 'master' of carpetgit@carpetcode.dyndns.org:carpet

14 files changed, 294 insertions, 116 deletions

diff --git a/Carpet/Carpet/param.ccl b/Carpet/Carpet/param.ccl
index 94b4c7e4e..74b2be338 100644
--- a/Carpet/Carpet/param.ccl
+++ b/Carpet/Carpet/param.ccl
@@ -440,6 +440,12 @@ CCTK_REAL aspect_ratio_z "Desired aspect ratio for each processor's domain" STEE
   (0:* :: ""
 } 1.0
 
+CCTK_INT min_points_per_proc "Minimum number of grid points per processor" STEERABLE=always
+{
+  -1  :: "infinity"
+  1:* :: "that many"
+} -1
+
 
 
 CCTK_INT num_threads "Number of threads per process" STEERABLE=recover
diff --git a/Carpet/Carpet/src/Recompose.cc b/Carpet/Carpet/src/Recompose.cc
index ac4732523..2d6c7beac 100644
--- a/Carpet/Carpet/src/Recompose.cc
+++ b/Carpet/Carpet/src/Recompose.cc
@@ -105,8 +105,8 @@ namespace Carpet {
           // assert (all(regsss.at(rl).at(c).at(ml).extent.lower() <=
           //             regsss.at(rl).at(c).at(ml).extent.upper()));
           // Check strides
-          const ivect str
-            = (maxspacereflevelfact / spacereffacts.at(rl) * ipow(mgfact, ml));
+          ivect const str =
+            (maxspacereflevelfact / spacereffacts.at(rl) * ipow(mgfact, ml));
           assert (all(regsss.at(ml).at(rl).at(c).extent.stride() % str == 0));
           // Check alignments
           assert (all(regsss.at(ml).at(rl).at(c).extent.lower() % str == 0));
@@ -1314,6 +1314,8 @@ namespace Carpet {
                               vector<vector<region_t> > & superregss,
                               vector<vector<region_t> > & regss)
   {
+    DECLARE_CCTK_PARAMETERS;
+    
     if (DEBUG) cout << "SRMA enter" << endl;
     
     int const nmaps = superregss.size();
@@ -1358,7 +1360,21 @@ namespace Carpet {
       }
     }
     
-    const int nprocs = CCTK_nProcs (cctkGH);
+    int const real_nprocs = CCTK_nProcs (cctkGH);
+    if (DEBUG) cout << "SRMA real_nprocs " << real_nprocs << endl;
+    
+    // Deactivate some processors if there are too many
+    int nprocs;
+    if (min_points_per_proc < 0) {
+      nprocs = real_nprocs;
+    } else {
+      CCTK_REAL mycost = 0;
+      for (int r=0; r<nregs; ++r) {
+        mycost += prod (cost (superregs.at(r)));
+      }
+      int const goodnprocs = int (floor (mycost / min_points_per_proc));
+      nprocs = max (1, min (real_nprocs, goodnprocs));
+    }
     if (DEBUG) cout << "SRMA nprocs " << nprocs << endl;
     
     // ncomps: number of components per processor
@@ -1453,7 +1469,6 @@ namespace Carpet {
         for (ipfulltree::iterator fti (* regf); not fti.done(); ++ fti) {
           pseudoregion_t & preg = (* fti).payload();
           preg.component = tmpncomps.at(m)++;
-          // preg.processor /= ncomps;
         }
       }
       for (int m=0; m<nmaps; ++m) {
@@ -1481,6 +1496,11 @@ namespace Carpet {
       assert (m>=0 and m<nmaps);
       superregss.at(m).push_back (superregs.at(r));
     }
+    // Output regions
+    if (DEBUG) {
+      cout << "SRMA superregss " << superregss << endl;
+      cout << "SRMA regss " << regss << endl;
+    }
     // Check sizes
     for (int m=0; m<nmaps; ++m) {
       assert (int(regss.at(m).size()) == myncomps.at(m) + empty_comps.at(m));
diff --git a/Carpet/Carpet/src/SetupGH.cc b/Carpet/Carpet/src/SetupGH.cc
index 0c9c37958..ed3f01bfe 100644
--- a/Carpet/Carpet/src/SetupGH.cc
+++ b/Carpet/Carpet/src/SetupGH.cc
@@ -854,6 +854,45 @@ namespace Carpet {
       assert (0);
     }
     
+    // Add empty regions if there are fewer regions than processors
+    {
+      int const nprocs = CCTK_nProcs (cctkGH);
+      int const oldsize = regs.size();
+      if (oldsize < nprocs) {
+        // Ensure that there are at least nprocs components
+        
+        // Construct an empty bbox that is just to the right of the
+        // last bbox, and which has the correct stride
+        assert (not regs.empty());
+        ibbox const & lbox = (*regs.rbegin()).extent;
+        ibbox const ebox
+          (lbox.upper() + lbox.stride(), lbox.upper(), lbox.stride());
+        
+        regs.resize (nprocs);
+        for (int c=oldsize; c<nprocs; ++c) {
+          region_t empty;
+          empty.extent = ebox;
+          empty.processor = c;
+          regs.at(c) = empty;
+        }
+      }
+    }
+    
+    // Check processor distribution
+    {
+      int const nprocs = CCTK_nProcs (cctkGH);
+      vector<bool> used (nprocs, false);
+      for (int c=0; c<nprocs; ++c) {
+        int const p = regs.at(c).processor;
+        assert (p >= 0 and p < nprocs);
+        assert (not used.at(p));
+        used.at(p) = true;
+      }
+      for (int c=0; c<nprocs; ++c) {
+        assert (used.at(c));
+      }
+    }
+    
     // Only one refinement level
     vector<vector<region_t> > superregss(1);
     superregss.at(rl) = superregs;
diff --git a/Carpet/CarpetInterp2/interface.ccl b/Carpet/CarpetInterp2/interface.ccl
index fc822c9e1..1a8e1c177 100644
--- a/Carpet/CarpetInterp2/interface.ccl
+++ b/Carpet/CarpetInterp2/interface.ccl
@@ -41,7 +41,24 @@ CCTK_INT FUNCTION                                \
          CCTK_INT              IN  npoints,      \
          CCTK_POINTER_TO_CONST IN  globalcoords, \
          CCTK_INT        ARRAY OUT patch,        \
-         CCTK_POINTER          IN localcoords,   \
-         CCTK_POINTER          IN dadx,          \
-         CCTK_POINTER          IN ddadxdx)
+         CCTK_POINTER          IN  localcoords,  \
+         CCTK_POINTER          IN  dadx,         \
+         CCTK_POINTER          IN  ddadxdx)
 USES FUNCTION MultiPatch_GlobalToLocal
+
+
+
+CCTK_INT FUNCTION                                      \
+    InterpGridArrays                                   \
+        (CCTK_POINTER_TO_CONST IN cctkGH,              \
+         CCTK_INT              IN N_dims,              \
+         CCTK_INT              IN order,               \
+         CCTK_INT              IN N_interp_points,     \
+         CCTK_POINTER_TO_CONST IN interp_coords,       \
+         CCTK_INT              IN N_input_arrays,      \
+         CCTK_INT        ARRAY IN input_array_indices, \
+         CCTK_INT              IN N_output_arrays,     \
+         CCTK_POINTER          IN output_arrays)
+PROVIDES FUNCTION InterpGridArrays      \
+    WITH CarpetInterp2_InterpGridArrays \
+    LANGUAGE C
diff --git a/Carpet/CarpetInterp2/src/interp2.cc b/Carpet/CarpetInterp2/src/interp2.cc
new file mode 100644
index 000000000..f8a17df99
--- /dev/null
+++ b/Carpet/CarpetInterp2/src/interp2.cc
@@ -0,0 +1,78 @@
+#include <cctk.h>
+
+#include "fasterp.hh"
+
+
+
+namespace CarpetInterp2 {
+  
+  extern "C"
+  CCTK_INT
+  CarpetInterp2_InterpGridArrays (CCTK_POINTER_TO_CONST const cctkGH_,
+                                  CCTK_INT const N_dims,
+                                  CCTK_INT const order,
+                                  CCTK_INT const N_interp_points,
+                                  CCTK_POINTER_TO_CONST const interp_coords_,
+                                  CCTK_INT const N_input_arrays,
+                                  CCTK_INT const * const input_array_indices,
+                                  CCTK_INT const N_output_arrays,
+                                  CCTK_POINTER const output_arrays_)
+  {
+    // Check input values and convert types
+    cGH const * const cctkGH = static_cast<cGH const *> (cctkGH_);
+    assert (cctkGH);
+    
+    assert (N_dims == dim);
+    
+    assert (N_interp_points >= 0);
+    
+    CCTK_REAL const * const * const interp_coords =
+      static_cast<CCTK_REAL const * const *> (interp_coords_);
+    assert (interp_coords);
+    for (int d=0; d<dim; ++d) {
+      assert (N_interp_points==0 or interp_coords[d]);
+    }
+    
+    assert (N_input_arrays >= 0);
+    
+    assert (input_array_indices);
+    for (int n=0; n<N_input_arrays; ++n) {
+      assert (input_array_indices[n]>=0 and
+              input_array_indices[n]<CCTK_NumVars());
+    }
+    
+    assert (N_output_arrays >= 0);
+    assert (N_output_arrays == N_input_arrays);
+    
+    CCTK_REAL * const * const output_arrays =
+      static_cast<CCTK_REAL * const *> (output_arrays_);
+    assert (output_arrays);
+    for (int n=0; n<N_output_arrays; ++n) {
+      assert (N_output_arrays==0 or output_arrays[n]);
+    }
+    
+    // Set up interpolation
+    fasterp_glocs_t locations (N_interp_points);
+    for (int d=0; d<dim; ++d) {
+      for (int i=0; i<N_interp_points; ++i) {
+        locations.coords[d].AT(i) = interp_coords[d][i];
+      }
+    }
+    fasterp_setup_t const setup (cctkGH, locations, order);
+    
+    // Interpolate
+    vector<int> varinds (N_input_arrays);
+    for (int n=0; n<N_input_arrays; ++n) {
+      varinds.AT(n) = input_array_indices[n];
+    }
+    vector<CCTK_REAL *> values (N_input_arrays);
+    for (int n=0; n<N_input_arrays; ++n) {
+      values.AT(n) = output_arrays[n];
+    }
+    setup.interpolate (cctkGH, varinds, values);
+    
+    // Done
+    return 0;
+  }
+  
+} // namespace CarpetInterp2
diff --git a/Carpet/CarpetInterp2/src/make.code.defn b/Carpet/CarpetInterp2/src/make.code.defn
index d95df47d6..323b3652b 100644
--- a/Carpet/CarpetInterp2/src/make.code.defn
+++ b/Carpet/CarpetInterp2/src/make.code.defn
@@ -1,7 +1,7 @@
 # Main make.code.defn file for thorn CarpetInterp2
 
 # Source files in this directory
-SRCS = fasterp.cc
+SRCS = fasterp.cc interp2.cc
 
 # Subdirectories containing source files
 SUBDIRS = 
diff --git a/Carpet/CarpetLib/src/bbox.cc b/Carpet/CarpetLib/src/bbox.cc
index 6183f17e5..a14e1e9a3 100644
--- a/Carpet/CarpetLib/src/bbox.cc
+++ b/Carpet/CarpetLib/src/bbox.cc
@@ -160,7 +160,8 @@ bool bbox<T,D>::operator> (const bbox& b) const {
 template<class T, int D>
 bbox<T,D> bbox<T,D>::expand (const vect<T,D>& lo, const vect<T,D>& hi) const {
   // Allow expansion only into directions where the extent is not negative
-  assert (all(lower()<=upper() or (lo==T(0) and hi==T(0))));
+  // assert (all(lower()<=upper() or (lo==T(0) and hi==T(0))));
+  assert (all(shape()>=vect<T,D>(0) or (lo==T(0) and hi==T(0))));
   const vect<T,D> str = stride();
   const vect<T,D> lb = lower() - lo * str;
   const vect<T,D> ub = upper() + hi * str;
diff --git a/Carpet/CarpetMask/param.ccl b/Carpet/CarpetMask/param.ccl
index 8db65e41b..f63bceab3 100644
--- a/Carpet/CarpetMask/param.ccl
+++ b/Carpet/CarpetMask/param.ccl
@@ -36,6 +36,10 @@ CCTK_REAL excluded_radius[10] "radius of excluded region" STEERABLE=always
   0.0:* :: ""
 } -1.0
 
+BOOLEAN exclude_exterior[10] "exclude the exterior of this radius, otherwise the interior" STEERABLE=always
+{
+} no
+
 
 
 # Exclude spherical surfaces
diff --git a/Carpet/CarpetMask/src/mask_excluded.cc b/Carpet/CarpetMask/src/mask_excluded.cc
index 97268ef42..e842caf31 100644
--- a/Carpet/CarpetMask/src/mask_excluded.cc
+++ b/Carpet/CarpetMask/src/mask_excluded.cc
@@ -44,6 +44,8 @@ namespace CarpetMask {
         
         CCTK_REAL const r2 = pow (r0, 2);
         
+        bool const exterior = exclude_exterior[n];
+        
         for (int k = 0; k < cctk_lsh[2]; ++ k) {
           for (int j = 0; j < cctk_lsh[1]; ++ j) {
             for (int i = 0; i < cctk_lsh[0]; ++ i) {
@@ -53,7 +55,10 @@ namespace CarpetMask {
               CCTK_REAL const dy2 = pow (y[ind] - y0, 2);
               CCTK_REAL const dz2 = pow (z[ind] - z0, 2);
               
-              if (dx2 + dy2 + dz2 <= r2) {
+              if (exterior ?
+                  dx2 + dy2 + dz2 >= r2 :
+                  dx2 + dy2 + dz2 <= r2)
+              {
                 weight[ind] = 0.0;
               }
               
diff --git a/Carpet/CarpetReduce/configuration.ccl b/Carpet/CarpetReduce/configuration.ccl
index 50ee465ed..f81217393 100644
--- a/Carpet/CarpetReduce/configuration.ccl
+++ b/Carpet/CarpetReduce/configuration.ccl
@@ -1,5 +1,5 @@
 # Configuration definitions for thorn CarpetReduce
 
-REQUIRES Carpet CarpetLib
+REQUIRES Carpet CarpetLib LoopControl
 
-REQUIRES THORNS: Carpet CarpetLib
+REQUIRES THORNS: Carpet CarpetLib LoopControl
diff --git a/Carpet/CarpetReduce/interface.ccl b/Carpet/CarpetReduce/interface.ccl
index e19087aa9..467607cb5 100644
--- a/Carpet/CarpetReduce/interface.ccl
+++ b/Carpet/CarpetReduce/interface.ccl
@@ -10,6 +10,8 @@ uses include header: carpet.hh
 
 uses include header: carpet_typecase.hh
 
+uses include header: loopcontrol.h
+
 CCTK_INT FUNCTION \
     SymmetryTableHandleForGrid (CCTK_POINTER_TO_CONST IN cctkGH)
 REQUIRES FUNCTION SymmetryTableHandleForGrid
diff --git a/Carpet/CarpetReduce/src/mask_carpet.cc b/Carpet/CarpetReduce/src/mask_carpet.cc
index 5593dfbfa..679c85392 100644
--- a/Carpet/CarpetReduce/src/mask_carpet.cc
+++ b/Carpet/CarpetReduce/src/mask_carpet.cc
@@ -7,6 +7,8 @@
 
 #include <carpet.hh>
 
+#include <loopcontrol.h>
+
 #include "mask_carpet.hh"
 
 
@@ -133,15 +135,15 @@ namespace CarpetMask {
                 
                 // Set weight on the boundary to 1/2
                 assert (dim == 3);
-#pragma omp parallel for
-                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 ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
-                      weight[ind] *= 0.5;
-                    }
-                  }
-                }
+#pragma omp parallel
+                LC_LOOP3(CarpetMaskSetup_prolongation,
+                         i,j,k,
+                         imin[0],imin[1],imin[2], imax[0],imax[1],imax[2],
+                         cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
+                {
+                  int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
+                  weight[ind] *= 0.5;
+                } LC_ENDLOOP3(CarpetMaskSetup_prolongation);
                 
               } // if box not empty
               
@@ -196,15 +198,15 @@ namespace CarpetMask {
               
               // Set weight in the restricted region to 0
               assert (dim == 3);
-#pragma omp parallel for
-              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 ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
-                    weight[ind] = 0;
-                  }
-                }
-              }
+#pragma omp parallel
+              LC_LOOP3(CarpetMaskSetup_restriction,
+                       i,j,k,
+                       imin[0],imin[1],imin[2], imax[0],imax[1],imax[2],
+                       cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
+              {
+                int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
+                weight[ind] = 0;
+              } LC_ENDLOOP3(CarpetMaskSetup_restriction);
               
             } // if box not empty
               
@@ -214,15 +216,15 @@ namespace CarpetMask {
           vector<int> mask (cctk_lsh[0] * cctk_lsh[1] * cctk_lsh[2]);
           
           assert (dim == 3);
-#pragma omp parallel for
-          for (int k=0; k<cctk_lsh[2]; ++k) {
-            for (int j=0; j<cctk_lsh[1]; ++j) {
-              for (int i=0; i<cctk_lsh[0]; ++i) {
-                int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
-                mask[ind] = 0;
-              }
-            }
-          }
+#pragma omp parallel
+          LC_LOOP3(CarpetMaskSetup_restriction_boundary_init,
+                   i,j,k,
+                   0,0,0, cctk_lsh[0],cctk_lsh[1],cctk_lsh[2],
+                   cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
+          {
+            int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
+            mask[ind] = 0;
+          } LC_ENDLOOP3(CarpetMaskSetup_restriction_boundary_init);
           
           for (int d=0; d<dim; ++d) {
             for (ibset::const_iterator bi = boundaries[d].begin();
@@ -251,18 +253,18 @@ namespace CarpetMask {
                 
                 // Set weight on the boundary to 1/2
                 assert (dim == 3);
-#pragma omp parallel for
-                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 ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
-                      if (mask[ind] == 0) {
-                        mask[ind] = 1;
-                      }
-                      mask[ind] *= 2;
-                    }
+#pragma omp parallel
+                LC_LOOP3(CarpetMaskSetup_restriction_boundary_partial,
+                         i,j,k,
+                         imin[0],imin[1],imin[2], imax[0],imax[1],imax[2],
+                         cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
+                {
+                  int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
+                  if (mask[ind] == 0) {
+                    mask[ind] = 1;
                   }
-                }
+                  mask[ind] *= 2;
+                } LC_ENDLOOP3(CarpetMaskSetup_restriction_boundary_partial);
                 
               } // if box not empty
               
@@ -270,17 +272,17 @@ namespace CarpetMask {
           } // for d
           
           assert (dim == 3);
-#pragma omp parallel for
-          for (int k=0; k<cctk_lsh[2]; ++k) {
-            for (int j=0; j<cctk_lsh[1]; ++j) {
-              for (int i=0; i<cctk_lsh[0]; ++i) {
-                int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
-                if (mask[ind] > 0) {
-                  weight[ind] *= 1.0 - 1.0 / mask[ind];
-                }
-              }
+#pragma omp parallel
+          LC_LOOP3(CarpetMaskSetup_restriction_boundary_apply,
+                   i,j,k,
+                   0,0,0, cctk_lsh[0],cctk_lsh[1],cctk_lsh[2],
+                   cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
+          {
+            int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
+            if (mask[ind] > 0) {
+              weight[ind] *= 1.0 - 1.0 / mask[ind];
             }
-          }
+          } LC_ENDLOOP3(CarpetMaskSetup_restriction_boundary_apply);
           
         } END_LOCAL_COMPONENT_LOOP;
         
diff --git a/Carpet/CarpetReduce/src/mask_coords.c b/Carpet/CarpetReduce/src/mask_coords.c
index 654cfee40..a5ef4324d 100644
--- a/Carpet/CarpetReduce/src/mask_coords.c
+++ b/Carpet/CarpetReduce/src/mask_coords.c
@@ -4,6 +4,8 @@
 #include <cctk_Arguments.h>
 #include <cctk_Parameters.h>
 
+#include <loopcontrol.h>
+
 
 
 void
@@ -49,7 +51,7 @@ CoordBase_SetupMask (CCTK_ARGUMENTS)
       
       if (is_internal[2*d+f]) {
         /* The boundary extends inwards */
-        bnd_points[2*d+f] = shiftout[2*d+f];
+        bnd_points[2*d+f] = shiftout[2*d+f] - is_staggered[2*d+f];
       } else {
         /* The boundary extends outwards */
         bnd_points[2*d+f] =
@@ -113,17 +115,17 @@ CoordBase_SetupMask (CCTK_ARGUMENTS)
             CCTK_VInfo (CCTK_THORNSTRING,
                         "Setting boundary points in direction %d face %d to weight 0", d, f);
           }
-#pragma omp parallel for
-          for (int k=bmin[2]; k<bmax[2]; ++k) {
-            for (int j=bmin[1]; j<bmax[1]; ++j) {
-              for (int i=bmin[0]; i<bmax[0]; ++i) {
-                
-                int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
-                weight[ind] = 0.0;
-                
-              }
-            }
-          }
+#pragma omp parallel
+          LC_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);
+            weight[ind] = 0.0;
+            
+          } LC_ENDLOOP3(CoordBase_SetupMask);
           
           /* When the boundary is not staggered, then give the points
              on the boundary the weight 1/2 */
@@ -159,17 +161,17 @@ CoordBase_SetupMask (CCTK_ARGUMENTS)
                 CCTK_VInfo (CCTK_THORNSTRING,
                             "Setting non-staggered boundary points in direction %d face %d to weight 1/2", d, f);
               }
-#pragma omp parallel for
-              for (int k=bmin[2]; k<bmax[2]; ++k) {
-                for (int j=bmin[1]; j<bmax[1]; ++j) {
-                  for (int i=bmin[0]; i<bmax[0]; ++i) {
-                    
-                    int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
-                    weight[ind] *= 0.5;
-                    
-                  }
-                }
-              }
+#pragma omp parallel
+              LC_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);
+                weight[ind] *= 0.5;
+                
+              } LC_ENDLOOP3(CoordBase_SetupMask_boundary2);
               
             } /* if the domain is not empty */
             
diff --git a/Carpet/CarpetReduce/src/mask_init.c b/Carpet/CarpetReduce/src/mask_init.c
index 129f95dcc..c1545a4fe 100644
--- a/Carpet/CarpetReduce/src/mask_init.c
+++ b/Carpet/CarpetReduce/src/mask_init.c
@@ -4,6 +4,8 @@
 
 #include <util_Table.h>
 
+#include <loopcontrol.h>
+
 
 
 void
@@ -29,17 +31,17 @@ MaskBase_InitMask (CCTK_ARGUMENTS)
     CCTK_INFO ("Initialising to weight 1");
   }
   
-#pragma omp parallel for
-  for (int k=0; k<cctk_lsh[2]; ++k) {
-    for (int j=0; j<cctk_lsh[1]; ++j) {
-      for (int i=0; i<cctk_lsh[0]; ++i) {
-        
-        int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
-        weight[ind] = 1.0;
-        
-      }
-    }
-  }
+#pragma omp parallel
+  LC_LOOP3(MaskBase_InitMask_interior,
+           i,j,k,
+           0,0,0, cctk_lsh[0],cctk_lsh[1],cctk_lsh[2],
+           cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
+  {
+    
+    int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
+    weight[ind] = 1.0;
+    
+  } LC_ENDLOOP3(MaskBase_InitMask_interior);
   
   
   
@@ -69,17 +71,17 @@ MaskBase_InitMask (CCTK_ARGUMENTS)
         }
         
         /* Loop over the boundary slab */
-#pragma omp parallel for
-        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 ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
-              weight[ind] = 0.0;
-              
-            }
-          }
-        }
+#pragma omp parallel
+        LC_LOOP3(MaskBase_InitMask_ghosts,
+                 i,j,k,
+                 imin[0],imin[1],imin[2], imax[0],imax[1],imax[2],
+                 cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
+        {
+          
+          int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
+          weight[ind] = 0.0;
+          
+        } LC_ENDLOOP3(MaskBase_InitMask_ghosts);
         
       }
     }
@@ -138,17 +140,17 @@ MaskBase_InitMask (CCTK_ARGUMENTS)
           }
           
           /* Loop over the boundary slab */
-#pragma omp parallel for
-          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 ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
-                weight[ind] = 0.0;
-                
-              }
-            }
-          }
+#pragma omp parallel
+          LC_LOOP3(MaskBase_InitMask_boundary,
+                   i,j,k,
+                   imin[0],imin[1],imin[2], imax[0],imax[1],imax[2],
+                   cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
+          {
+            
+            int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
+            weight[ind] = 0.0;
+            
+          } LC_ENDLOOP3(MaskBase_InitMask_boundary);
           
         }
       }