CarpetLib: Add routine to self-check regions and superregions

Add routine to self-check the datastructure representing regions and
superregion.

Call this routine after recomposing.

3 files changed, 134 insertions, 0 deletions

diff --git a/Carpet/Carpet/src/Recompose.cc b/Carpet/Carpet/src/Recompose.cc
index 8c9758039..b8bba7d96 100644
--- a/Carpet/Carpet/src/Recompose.cc
+++ b/Carpet/Carpet/src/Recompose.cc
@@ -1308,6 +1308,35 @@ namespace Carpet {
     } else {
       assert (0);
     }
+    
+    for (size_t m=0; m<superregss.size(); ++m) {
+      for (size_t r=0; r<superregss.AT(m).size(); ++r) {
+        bool const good_superregs = superregss.AT(m).AT(r).check_region(true);
+        if (not good_superregs) {
+          cout << "superregs[" << m << "][" << r << "]:\n"
+               << superregss.AT(m).AT(r) << "\n";
+          cout << "all superregions:\n" << superregss << "\n";
+          CCTK_VWarn(CCTK_WARN_ALERT, __LINE__, __FILE__, CCTK_THORNSTRING,
+                     "Superregions failed self-check for map %d superregion %d",
+                     int(m), int(r));
+        }
+      }
+    }
+    assert(regss.size() == superregss.size());
+    for (size_t m=0; m<regss.size(); ++m) {
+      for (size_t r=0; r<regss.AT(m).size(); ++r) {
+        bool const good_regs = regss.AT(m).AT(r).check_region(false);
+        if (not good_regs) {
+          cout << "regs[" << m << "][" << r << "]:\n"
+               << regss.AT(m).AT(r) << "\n";
+          cout << "all superregions:\n" << superregss << "\n";
+          cout << "all regions:\n" << regss << "\n";
+          CCTK_VWarn(CCTK_WARN_ALERT, __LINE__, __FILE__, CCTK_THORNSTRING,
+                     "Regions failed self-check for map %d superregion %d",
+                     int(m), int(r));
+        }
+      }
+    }
   }
   
   
diff --git a/Carpet/CarpetLib/src/region.cc b/Carpet/CarpetLib/src/region.cc
index a709364dd..09fb7d5f0 100644
--- a/Carpet/CarpetLib/src/region.cc
+++ b/Carpet/CarpetLib/src/region.cc
@@ -4,6 +4,7 @@
 
 #include "bboxset.hh"
 #include "defs.hh"
+#include "dist.hh"
 #include "mpi_string.hh"
 #include "region.hh"
 
@@ -280,6 +281,103 @@ operator>> (istream & is, region_t & reg)
 
 
 
+bool region_t::check_region(bool const is_superregion) const
+{
+  assert(invariant());
+  // TODO: empty regions need to be allowed (e.g. for empty grid
+  // arrays, or when the domain decomposition assigns no points to a
+  // particular process)
+  // if (extent.empty()) {
+  //   CCTK_WARN(CCTK_WARN_PICKY, "extent is empty");
+  //   return false;
+  // }
+  assert(map >= 0);
+  
+  if (is_superregion) {
+    // Checking a superregion: a tree without processor assignments
+    
+    if (processor != -1) {
+      CCTK_WARN(CCTK_WARN_PICKY, "processor number is defined");
+      return false;
+    }
+    if (not processors) {
+      CCTK_WARN(CCTK_WARN_PICKY, "processor tree not defined");
+      return false;
+    }
+    
+    // TODO: check outer_boundaries as well
+    ibset child_extents;
+    check_children(*processors, map, 0, child_extents);
+    if (child_extents != extent) {
+      CCTK_WARN(CCTK_WARN_PICKY, "child extents not equal to total extent");
+      return false;
+    }
+    
+  } else {
+    // Checking a regular region: no tree structure, but has a
+    // processor assignment
+    
+    if (processor < 0 or processor >= dist::size()) {
+      CCTK_WARN(CCTK_WARN_PICKY, "processor number not defined");
+      return false;
+    }
+    if (processors) {
+      CCTK_WARN(CCTK_WARN_PICKY, "processor tree is defined");
+      return false;
+    }
+    
+  }
+  
+  return true;
+}
+
+bool region_t::check_children(ipfulltree const& tree,
+                              int const parent_map,
+                              int const level,
+                              ibset& child_extents)
+  const
+{
+  if (tree.empty()) {
+    CCTK_VWarn(CCTK_WARN_PICKY, __LINE__, __FILE__, CCTK_THORNSTRING,
+               "level %d: tree is empty", level);
+    return false;
+  } else if (tree.is_branch()) {
+    for (ipfulltree::const_iterator i(tree); not i.done(); ++i) {
+      bool const bres = check_children(*i, parent_map, level+1, child_extents);
+      if (not bres) {
+        CCTK_VWarn(CCTK_WARN_PICKY, __LINE__, __FILE__, CCTK_THORNSTRING,
+                   "level %d: error in children", level);
+        return false;
+      }
+    }
+  } else if (tree.is_leaf()) {
+    ibbox const& ext = tree.payload().extent;
+    if (ext.empty()) {
+      CCTK_VWarn(CCTK_WARN_PICKY, __LINE__, __FILE__, CCTK_THORNSTRING,
+                 "level %d: tree leaf extent is empty", level);
+      return false;
+    }
+    if (map != parent_map) {
+      CCTK_VWarn(CCTK_WARN_PICKY, __LINE__, __FILE__, CCTK_THORNSTRING,
+                 "level %d: tree leaf map=%d, differs from parent map=%d",
+                 level, map, parent_map);
+      return false;
+    }
+    if (child_extents.intersects(ext)) {
+      CCTK_VWarn(CCTK_WARN_PICKY, __LINE__, __FILE__, CCTK_THORNSTRING,
+                 "level %d: tree leaf extent overlaps other tree extents",
+                 level);
+      return false;
+    }
+    child_extents += ext;
+  } else {
+    assert(0);
+  }
+  return true;
+}
+
+
+
 bool region_t::full_output = false;
 
 ostream &
diff --git a/Carpet/CarpetLib/src/region.hh b/Carpet/CarpetLib/src/region.hh
index 4c8bf25ac..e5f36011c 100644
--- a/Carpet/CarpetLib/src/region.hh
+++ b/Carpet/CarpetLib/src/region.hh
@@ -31,6 +31,13 @@ struct region_t {
   CCTK_REAL load () const CCTK_MEMBER_ATTRIBUTE_PURE;
   region_t split (CCTK_REAL ratio_new_over_old);
   
+  // Check whether a region is defined consistently
+  bool check_region(bool is_superregion) const;
+private:
+  bool check_children(ipfulltree const& tree, int parent_map, int level,
+                      ibset& child_extents) const;
+public:
+  
   // Output processor decomposition? (Off by default.)
   static bool full_output;
 };