Enhance support for multi-patch simulations

Carpet:

Ensure that at most one of GetDomainSpecificatio or
MultiPatch_GetDomainSpecification is defined.

Allow the boundary case of having zero regions on a refinement level.

CarpetLib:

Allow the boundary case of having zero components on a patch, if there
are still more than zero components overall.

CarpetIOASCII:

Loop over the components using explicit light-weight for loops instead of
Carpet's looping macros.  This is faster, since less state information
needs to be updated.

Correct an inconsistency in converting between integer indices and
coordinates in multi-patch simulations.

CarpetRegrid2:

Take multi-patch systems into account.

2 files changed, 159 insertions, 36 deletions

diff --git a/Carpet/CarpetRegrid2/interface.ccl b/Carpet/CarpetRegrid2/interface.ccl
index 7b55ed523..aa6f799da 100644
--- a/Carpet/CarpetRegrid2/interface.ccl
+++ b/Carpet/CarpetRegrid2/interface.ccl
@@ -47,6 +47,42 @@ CCTK_INT FUNCTION ConvertFromPhysicalBoundary \
      CCTK_REAL IN  ARRAY spacing)
 USES FUNCTION ConvertFromPhysicalBoundary
 
+# The location of the boundary points
+CCTK_INT FUNCTION MultiPatch_GetBoundarySpecification \
+    (CCTK_INT IN        map,                          \
+     CCTK_INT IN        size,                         \
+     CCTK_INT OUT ARRAY nboundaryzones,               \
+     CCTK_INT OUT ARRAY is_internal,                  \
+     CCTK_INT OUT ARRAY is_staggered,                 \
+     CCTK_INT OUT ARRAY shiftout)
+USES FUNCTION MultiPatch_GetBoundarySpecification
+
+# The overall size of the domain
+CCTK_INT FUNCTION MultiPatch_GetDomainSpecification \
+    (CCTK_INT  IN        map,                       \
+     CCTK_INT  IN        size,                      \
+     CCTK_REAL OUT ARRAY physical_min,              \
+     CCTK_REAL OUT ARRAY physical_max,              \
+     CCTK_REAL OUT ARRAY interior_min,              \
+     CCTK_REAL OUT ARRAY interior_max,              \
+     CCTK_REAL OUT ARRAY exterior_min,              \
+     CCTK_REAL OUT ARRAY exterior_max,              \
+     CCTK_REAL OUT ARRAY spacing)
+USES FUNCTION MultiPatch_GetDomainSpecification
+
+# Conversion between boundary types
+CCTK_INT FUNCTION MultiPatch_ConvertFromPhysicalBoundary \
+    (CCTK_INT  IN        map,                            \
+     CCTK_INT  IN        size,                           \
+     CCTK_REAL IN  ARRAY physical_min,                   \
+     CCTK_REAL IN  ARRAY physical_max,                   \
+     CCTK_REAL OUT ARRAY interior_min,                   \
+     CCTK_REAL OUT ARRAY interior_max,                   \
+     CCTK_REAL OUT ARRAY exterior_min,                   \
+     CCTK_REAL OUT ARRAY exterior_max,                   \
+     CCTK_REAL IN  ARRAY spacing)
+USES FUNCTION MultiPatch_ConvertFromPhysicalBoundary
+
 
 
 # The true prototype of the routine below:
diff --git a/Carpet/CarpetRegrid2/src/regrid.cc b/Carpet/CarpetRegrid2/src/regrid.cc
index 983f7acba..c08705d29 100644
--- a/Carpet/CarpetRegrid2/src/regrid.cc
+++ b/Carpet/CarpetRegrid2/src/regrid.cc
@@ -174,6 +174,68 @@ namespace CarpetRegrid2 {
   
   
   
+  void
+  get_physical_boundary (rvect & physical_lower,
+                         rvect & physical_upper,
+                         rvect & spacing)
+  {
+    rvect interior_lower, interior_upper;
+    rvect exterior_lower, exterior_upper;
+    if (CCTK_IsFunctionAliased ("MultiPatch_GetDomainSpecification")) {
+      assert (Carpet::map >= 0);
+      CCTK_INT const ierr = MultiPatch_GetDomainSpecification
+        (Carpet::map, dim,
+         & physical_lower[0], & physical_upper[0],
+         & interior_lower[0], & interior_upper[0],
+         & exterior_lower[0], & exterior_upper[0],
+         & spacing[0]);
+      assert (not ierr);
+    } else if (CCTK_IsFunctionAliased ("GetDomainSpecification")) {
+      CCTK_INT const ierr = GetDomainSpecification
+        (dim,
+         & physical_lower[0], & physical_upper[0],
+         & interior_lower[0], & interior_upper[0],
+         & exterior_lower[0], & exterior_upper[0],
+         & spacing[0]);
+      assert (not ierr);
+    } else {
+      assert (0);
+    }
+  }
+  
+  void
+  calculate_exterior_boundary (rvect const & physical_lower,
+                               rvect const & physical_upper,
+                               rvect       & exterior_lower,
+                               rvect       & exterior_upper,
+                               rvect const & spacing)
+  {
+    rvect interior_lower, interior_upper;
+    if (CCTK_IsFunctionAliased ("MultiPatch_ConvertFromPhysicalBoundary")) {
+      assert (Carpet::map >= 0);
+      CCTK_INT const ierr = MultiPatch_ConvertFromPhysicalBoundary
+        (Carpet::map, dim,
+         & physical_lower[0], & physical_upper[0],
+         & interior_lower[0], & interior_upper[0],
+         & exterior_lower[0], & exterior_upper[0],
+         & spacing[0]);
+      assert (not ierr);
+    } else if (CCTK_IsFunctionAliased ("ConvertFromPhysicalBoundary")) {
+      CCTK_INT const ierr =
+        ConvertFromPhysicalBoundary
+        (dim,
+         & physical_lower[0], & physical_upper[0],
+         & interior_lower[0], & interior_upper[0],
+         & exterior_lower[0], & exterior_upper[0],
+         & spacing[0]);
+      assert (not ierr);
+    } else {
+      assert (0);
+    }
+  }
+  
+  
+  
   extern "C" {
     CCTK_INT
     CarpetRegrid2_Regrid (CCTK_POINTER_TO_CONST const cctkGH_,
@@ -208,6 +270,7 @@ namespace CarpetRegrid2 {
 #warning "TODO: check this (for Carpet, and maybe also for CartGrid3D)"
 #warning "TODO: (the check that these two are consistent should be in Carpet)"
     
+#if 0
     typedef vect<vect<CCTK_INT,2>,3> jjvect;
     jjvect nboundaryzones;
     jjvect is_internal;
@@ -222,7 +285,9 @@ namespace CarpetRegrid2 {
          & shiftout[0][0]);
       assert (not ierr);
     }
+#endif
     
+#if 0
     rvect physical_lower, physical_upper;
     rvect interior_lower, interior_upper;
     rvect exterior_lower, exterior_upper;
@@ -236,10 +301,16 @@ namespace CarpetRegrid2 {
          & spacing[0]);
       assert (not ierr);
     }
+#endif
+    rvect physical_lower, physical_upper;
+    rvect spacing;
+    get_physical_boundary (physical_lower, physical_upper,
+                           spacing);
     
     // Adapt spacing for convergence level
     spacing *= ipow ((CCTK_REAL) Carpet::mgfact, Carpet::basemglevel);
     
+#if 0
     {
       CCTK_INT const ierr =
         ConvertFromPhysicalBoundary
@@ -250,6 +321,11 @@ namespace CarpetRegrid2 {
          & spacing[0]);
       assert (not ierr);
     }
+#endif
+    rvect exterior_lower, exterior_upper;
+    calculate_exterior_boundary (physical_lower, physical_upper,
+                                 exterior_lower, exterior_upper,
+                                 spacing);
     
     //
     // Calculate the union of the bounding boxes for all levels
@@ -274,38 +350,43 @@ namespace CarpetRegrid2 {
     }
     regions.at(0).normalize();
     
-    // Loop over all centres
-    for (int n = 0; n < num_centres; ++ n) {
-      centre_description centre (cctkGH, n);
-      if (centre.active) {
-        
-        // Loop over all levels for this centre
-        for (int rl = 1; rl < centre.num_levels; ++ rl) {
-          
-          // Calculate a bbox for this region
-          rvect const rmin = centre.position - centre.radius.at(rl);
-          rvect const rmax = centre.position + centre.radius.at(rl);
-          
-          // Convert to an integer bbox
-          ivect const istride = hh.baseextents.at(0).at(rl).stride();
-          ivect const imin =
-            rpos2ipos (rmin, origin, scale, hh, rl)
-            - boundary_shiftout * istride;
-          ivect const imax =
-            rpos2ipos1 (rmax, origin, scale, hh, rl)
-            + boundary_shiftout * istride;
-          
-          ibbox const region (imin, imax, istride);
+    // Refine only patch 0
+    if (Carpet::map == 0) {
+      // Loop over all centres
+      for (int n = 0; n < num_centres; ++ n) {
+        centre_description centre (cctkGH, n);
+        if (centre.active) {
           
-          // Add this region to the list of regions
-          if (static_cast <int> (regions.size()) < rl+1) regions.resize (rl+1);
-          regions.at(rl) |= region;
-          regions.at(rl).normalize();
+          // Loop over all levels for this centre
+          for (int rl = 1; rl < centre.num_levels; ++ rl) {
+            
+            // Calculate a bbox for this region
+            rvect const rmin = centre.position - centre.radius.at(rl);
+            rvect const rmax = centre.position + centre.radius.at(rl);
+            
+            // Convert to an integer bbox
+            ivect const istride = hh.baseextents.at(0).at(rl).stride();
+            ivect const imin =
+              rpos2ipos (rmin, origin, scale, hh, rl)
+              - boundary_shiftout * istride;
+            ivect const imax =
+              rpos2ipos1 (rmax, origin, scale, hh, rl)
+              + boundary_shiftout * istride;
+            
+            ibbox const region (imin, imax, istride);
+            
+            // Add this region to the list of regions
+            if (static_cast <int> (regions.size()) < rl+1) {
+              regions.resize (rl+1);
+            }
+            regions.at(rl) |= region;
+            regions.at(rl).normalize();
+            
+          } // for rl
           
-        } // for rl
-        
-      } // if centre is active
-    }   // for n
+        } // if centre is active
+      }   // for n
+    }     // if map==0
     
     
       
@@ -448,6 +529,7 @@ namespace CarpetRegrid2 {
       rvect const level_physical_lower = physical_lower;
       rvect const level_physical_upper = physical_upper;
       rvect const level_spacing = spacing / rvect (hh.reffacts.at(rl));
+#if 0
       rvect level_interior_lower, level_interior_upper;
       rvect level_exterior_lower, level_exterior_upper;
       {
@@ -460,6 +542,11 @@ namespace CarpetRegrid2 {
            & level_spacing[0]);
         assert (not ierr);
       }
+#endif
+      rvect level_exterior_lower, level_exterior_upper;
+      calculate_exterior_boundary (level_physical_lower, level_physical_upper,
+                                   level_exterior_lower, level_exterior_upper,
+                                   level_spacing);
       
       ivect const level_exterior_ilower =
         rpos2ipos (level_exterior_lower, origin, scale, hh, rl);
@@ -975,20 +1062,20 @@ namespace CarpetRegrid2 {
       for (int m = 0; m < maps; ++ m) {
         maxrl = max (maxrl, rls.at(m));
       }
+      // All maps must have the same number of levels
+      for (int m = 0; m < maps; ++ m) {
+        regsss.at(m).resize (maxrl);
+      }
       
       // Make multiprocessor aware
       for (int rl = 0; rl < maxrl; ++ rl) {
         vector <vector <region_t> > regss (maps);
         for (int m = 0; m < maps; ++ m) {
-          if (rl < rls.at(m)) {
-            regss.at(m) = regsss.at(m).at(rl);
-          }
+          regss.at(m) = regsss.at(m).at(rl);
         }
         Carpet::SplitRegionsMaps (cctkGH, regss);
         for (int m = 0; m < maps; ++ m) {
-          if (rl < rls.at(m)) {
-            regsss.at(m).at(rl) = regss.at(m);
-          }
+          regsss.at(m).at(rl) = regss.at(m);
         }
       } // for rl