CarpetReduce: Use LoopControl

3 files changed, 110 insertions, 104 deletions

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