CarpetIOHDF5: add option to not output the buffer zones

  This is implemented by re-computing the allactive bboxset that dh::regrid
  computes and outputting only the part of each component that intersects
  allactive. This changes the number of components since the intersection
  might be non-rectangular-box.

5 files changed, 509 insertions, 35 deletions

diff --git a/Carpet/CarpetIOHDF5/param.ccl b/Carpet/CarpetIOHDF5/param.ccl
index 2054d6918..11a7323e6 100644
--- a/Carpet/CarpetIOHDF5/param.ccl
+++ b/Carpet/CarpetIOHDF5/param.ccl
@@ -372,6 +372,10 @@ BOOLEAN output_boundary_points "Output outer boundary points (assuming that ther
 {
 } "yes"
 
+BOOLEAN output_buffer_points "Output refinement buffer points" STEERABLE = ALWAYS
+{
+} "yes"
+
 BOOLEAN out3D_ghosts "Output ghost zones (DEPRECATED)" STEERABLE = ALWAYS
 {
 } "yes"
diff --git a/Carpet/CarpetIOHDF5/src/CarpetIOHDF5.hh b/Carpet/CarpetIOHDF5/src/CarpetIOHDF5.hh
index e605a2ddc..f42c8f23b 100644
--- a/Carpet/CarpetIOHDF5/src/CarpetIOHDF5.hh
+++ b/Carpet/CarpetIOHDF5/src/CarpetIOHDF5.hh
@@ -147,6 +147,13 @@ namespace CarpetIOHDF5
 
   // Everything is a class template, so that it can easily be
   // instantiated for all output dimensions
+  
+  // computes the active region the way dh::regrid does
+  void GetAllActive (const dh *dd, 
+                     const gh *hh, 
+                     int ml, 
+                     int rl, 
+                     ibset &allactive);
 
   template<int outdim>
   struct IOHDF5 {
@@ -218,6 +225,7 @@ namespace CarpetIOHDF5
                            const int ml,
                            const int m,
                            const int c,
+                           const int output_component,
                            const int tl,
                            const CCTK_REAL coord_time,
                            const vect<CCTK_REAL,dim>& coord_lower,
diff --git a/Carpet/CarpetIOHDF5/src/GetAllActive.cc b/Carpet/CarpetIOHDF5/src/GetAllActive.cc
new file mode 100644
index 000000000..6dcca4602
--- /dev/null
+++ b/Carpet/CarpetIOHDF5/src/GetAllActive.cc
@@ -0,0 +1,427 @@
+ /*@@
+   @file      GetAllActive.cc
+   @date      Sun Sep 19 20:34:27 EDT 2010
+   @author    Erik Schnetter, Roland Haas
+   @desc
+              a piece of dh::regrid to compute the active region
+   @enddesc
+ @@*/
+
+#include <cassert>
+
+#include <vector>
+
+#include "cctk.h"
+
+#include "CarpetIOHDF5.hh"
+
+namespace CarpetIOHDF5 {
+
+  /********************************************************************
+   *********************     External Routines   **********************
+   ********************************************************************/
+
+  void
+  GetAllActive(const dh *dd, const gh *hh, int ml, int rl, ibset &allactive);
+
+
+  /********************************************************************
+   *********************     Internal Routines     ********************
+   ********************************************************************/
+
+  static
+  void
+  assert_error (char const * restrict const checkstring,
+                char const * restrict const file, int const line,
+                int const ml, int const rl,
+                char const * restrict const message);
+  static
+  void
+  assert_error (char const * restrict const checkstring,
+                char const * restrict const file, int const line,
+                int const ml, int const rl, int const c,
+                char const * restrict const message);
+  static
+  void
+  assert_error (char const * restrict const checkstring,
+                char const * restrict const file, int const line,
+                int const ml, int const rl, int const c, int const cc,
+                char const * restrict const message);
+
+  /********************************************************************
+   *********************     Local Data   *****************************
+   ********************************************************************/
+
+  // accumulates any consistency errors
+  static bool there_was_an_error = false;
+
+
+   /*@@
+     @routine    
+     @date       Sun Sep 19 20:46:35 EDT 2010
+     @author     Erik Schnetter, Roland Haas
+     @desc
+                 Output warning for consistency checks for a refinement level.
+                
+     @enddesc
+     @history
+     @endhistory
+     @var        checkstring
+     @vdesc      the check that failed
+     @vtype      char const * restrict const
+     @vio        in
+     @endvar
+     @var        file
+     @vdesc      name of source file where error occured
+     @vtype      char const * restrict const
+     @vio        in
+     @endvar
+     @var        line
+     @vdesc      source line in file where error occured
+     @vtype      int const
+     @vio        in
+     @endvar
+     @var        ml
+     @vdesc      meta level where error occured
+     @vtype      int const
+     @vio        in
+     @endvar
+     @var        rl
+     @vdesc      refinement level where error occured
+     @vtype      int const
+     @vio        in
+     @endvar
+     @var        message
+     @vdesc      error message to display
+     @vtype      char const * restrict const
+     @vio        in
+     @endvar
+     @returntype void
+  @@*/
+  static
+  void
+  assert_error (char const * restrict const checkstring,
+                char const * restrict const file, int const line,
+                int const ml, int const rl,
+                char const * restrict const message)
+  {
+    CCTK_VWarn (CCTK_WARN_ALERT, __LINE__, __FILE__, CCTK_THORNSTRING,
+                "\n%s:%d:\n   [ml=%d rl=%d] The following grid structure consistency check failed:\n   %s\n   %s",
+                file, line, ml, rl, message, checkstring);
+    there_was_an_error = true;
+  }
+
+   /*@@
+     @routine    
+     @date       Sun Sep 19 20:46:35 EDT 2010
+     @author     Erik Schnetter, Roland Haas
+     @desc
+                 Output warning for consistency checks for a single component.
+                
+     @enddesc
+     @history
+     @endhistory
+     @var        checkstring
+     @vdesc      the check that failed
+     @vtype      char const * restrict const
+     @vio        in
+     @endvar
+     @var        file
+     @vdesc      name of source file where error occured
+     @vtype      char const * restrict const
+     @vio        in
+     @endvar
+     @var        line
+     @vdesc      source line in file where error occured
+     @vtype      int const
+     @vio        in
+     @endvar
+     @var        ml
+     @vdesc      meta level where error occured
+     @vtype      int const
+     @vio        in
+     @endvar
+     @var        rl
+     @vdesc      refinement level where error occured
+     @vtype      int const
+     @vio        in
+     @endvar
+     @var        c
+     @vdesc      component where error occured
+     @vtype      int const
+     @vio        in
+     @endvar
+     @var        message
+     @vdesc      error message to display
+     @vtype      char const * restrict const
+     @vio        in
+     @endvar
+     @returntype void
+  @@*/
+  static
+  void
+  assert_error (char const * restrict const checkstring,
+                char const * restrict const file, int const line,
+                int const ml, int const rl, int const c,
+                char const * restrict const message)
+  {
+    CCTK_VWarn (CCTK_WARN_ALERT, __LINE__, __FILE__, CCTK_THORNSTRING,
+                "\n%s:%d:\n   [ml=%d rl=%d c=%d] The following grid structure consistency check failed:\n   %s\n   %s",
+                file, line, ml, rl, c, message, checkstring);
+    there_was_an_error = true;
+  }
+
+   /*@@
+     @routine    
+     @date       Sun Sep 19 20:46:35 EDT 2010
+     @author     Erik Schnetter, Roland Haas
+     @desc
+                 Output warning for consistency checks between two components.
+                
+     @enddesc
+     @history
+     @endhistory
+     @var        checkstring
+     @vdesc      the check that failed
+     @vtype      char const * restrict const
+     @vio        in
+     @endvar
+     @var        file
+     @vdesc      name of source file where error occured
+     @vtype      char const * restrict const
+     @vio        in
+     @endvar
+     @var        line
+     @vdesc      source line in file where error occured
+     @vtype      int const
+     @vio        in
+     @endvar
+     @var        ml
+     @vdesc      meta level where error occured
+     @vtype      int const
+     @vio        in
+     @endvar
+     @var        rl
+     @vdesc      refinement level where error occured
+     @vtype      int const
+     @vio        in
+     @endvar
+     @var        c
+     @vdesc      first component where error occured
+     @vtype      int const
+     @vio        in
+     @endvar
+     @var        cc
+     @vdesc      second component where error occured
+     @vtype      int const
+     @vio        in
+     @endvar
+     @var        message
+     @vdesc      error message to display
+     @vtype      char const * restrict const
+     @vio        in
+     @endvar
+     @returntype void
+  @@*/
+  static
+  void
+  assert_error (char const * restrict const checkstring,
+                char const * restrict const file, int const line,
+                int const ml, int const rl, int const c, int const cc,
+                char const * restrict const message)
+  {
+    CCTK_VWarn (CCTK_WARN_ALERT, __LINE__, __FILE__, CCTK_THORNSTRING,
+                "\n%s:%d:\n   [ml=%d rl=%d c=%d cc=%d] The following grid structure consistency check failed:\n   %s\n   %s",
+                file, line, ml, rl, c, cc, message, checkstring);
+    there_was_an_error = true;
+  }
+
+#ifdef CARPET_OPTIMISE
+
+  // For highest efficiency, omit all self-checks
+#define ASSERT_rl(check, message)
+#define ASSERT_c(check, message)
+#define ASSERT_cc(check, message)
+
+#else
+
+#define ASSERT_rl(check, message)                                       \
+    do {                                                                  \
+      if (not (check)) {                                                  \
+        assert_error (#check, __FILE__, __LINE__, ml, rl, message);       \
+      }                                                                   \
+    } while (false)
+
+#define ASSERT_c(check, message)                                        \
+    do {                                                                  \
+      if (not (check)) {                                                  \
+        assert_error (#check, __FILE__, __LINE__, ml, rl, c, message);    \
+      }                                                                   \
+    } while (false)
+
+#define ASSERT_cc(check, message)                                       \
+    do {                                                                  \
+      if (not (check)) {                                                  \
+        assert_error (#check, __FILE__, __LINE__, ml, rl, c, cc, message); \
+      }                                                                   \
+    } while (false)
+
+#endif
+
+   /*@@
+     @routine    
+     @date       Sun Sep 19 20:46:35 EDT 2010
+     @author     Erik Schnetter, Roland Haas
+     @desc
+                 Recompute the union of all active regions. This routine is a
+                 stripped down copy of dh::regrid.
+     @enddesc
+     @history
+     @endhistory
+     @var        dh
+     @vdesc      data hirarchy for this group
+     @vtype      const dh *
+     @vio        in
+     @endvar
+     @var        hh
+     @vdesc      grid hirarchy
+     @vtype      const gh *
+     @vio        in
+     @endvar
+     @var        ml
+     @vdesc      meta level to compute active region on
+     @vtype      int
+     @vio        in
+     @endvar
+     @var        rl
+     @vdesc      refinement level to compute active region on
+     @vtype      int
+     @vio        in
+     @endvar
+     @var        allactive
+     @vdesc      receives the list of active boxes
+     @vtype      ibset &
+     @vio        out
+     @endvar
+     @returntype void
+  @@*/
+  void
+  GetAllActive(const dh *dd, const gh *hh, int ml, int rl, ibset &allactive)
+  {
+    DECLARE_CCTK_PARAMETERS;
+      
+    // All owned regions
+    ibset allowned;
+    
+    // Domain:
+    
+    ibbox const & domain_exterior = hh->baseextent(ml,rl);
+    // Variables may have size zero
+    // ASSERT_rl (not domain_exterior.empty(),
+    //            "The exterior of the domain must not be empty");
+    
+    i2vect const & buffer_width = dd->buffer_widths.AT(rl);
+    i2vect const & boundary_width = hh->boundary_width;
+    ASSERT_rl (all (all (boundary_width >= 0)),
+               "The gh boundary widths must not be negative");
+    
+    ibbox const domain_active = domain_exterior.expand (- boundary_width);
+    // Variables may have size zero
+    // ASSERT_rl (not domain_active.empty(),
+    //            "The active part of the domain must not be empty");
+    ASSERT_rl (domain_active <= domain_exterior,
+               "The active part of the domain must be contained in the exterior part of the domain");
+    
+    for (int c = 0; c < hh->components(rl); ++ c) {
+      
+      // Interior:
+      
+      ibbox intr;
+      intr = ibbox::poison();
+      
+      // The interior of the grid has the extent as specified by the
+      // regridding thorn
+      intr = hh->extent (ml,rl,c);
+      
+      // (The interior must not be empty)
+      // Variables may have size zero
+      // ASSERT_c (not intr.empty(),
+      //           "The interior must not be empty");
+      
+      // The interior must be contained in the domain
+      ASSERT_c (intr <= hh->baseextent(ml,rl),
+                "The interior must be contained in the domain");
+      
+      // All interiors must be disjunct
+#ifdef CARPET_DEBUG
+      for (int cc = 0; cc < c; ++ cc) {
+        ASSERT_cc (not intr.intersects (full_level.AT(cc).interior),
+                   "All interiors must be disjunct");
+      }
+#endif
+      
+      
+      
+      // Outer boundary faces:
+      
+      b2vect is_outer_boundary;
+      
+      // The outer boundary faces are where the interior extends up
+      // to the outer boundary of the domain.  It is not possible to
+      // check whether it extends past the active part of the
+      // domain, since this would be wrong when the outer boundary
+      // width is zero.
+      is_outer_boundary[0] = intr.lower() == domain_exterior.lower(); 
+      is_outer_boundary[1] = intr.upper() == domain_exterior.upper(); 
+      
+      
+      
+      // Owned region:
+      
+      ibbox owned;
+      owned = ibbox::poison();
+      
+      owned = intr.expand (i2vect (is_outer_boundary) * (- boundary_width));
+      
+      // (The owned region must not be empty)
+      // Variables may have size zero
+      // ASSERT_c (not owned.empty(),
+      //           "The owned region must not be empty");
+      
+      // The owned region must be contained in the active part of
+      // the domain
+      ASSERT_c (owned <= domain_active,
+                "The owned region must be contained in the active part of the domain");
+      
+      // All owned regions must be disjunct
+#ifdef CARPET_DEBUG
+      for (int cc = 0; cc < c; ++ cc) {
+        ASSERT_cc (not owned.intersects (full_level.AT(cc).owned),
+                   "All owned regions must be disjunct");
+      }
+#endif
+
+       allowned |= owned;
+       ASSERT_rl (allowned <= domain_active,
+                "The owned regions must be contained in the active part of the domain");
+       
+     } // for c
+     
+     // Enlarge active part of domain
+     i2vect const safedist = i2vect (0);
+     ibbox const domain_enlarged = domain_active.expand (safedist);
+     
+     // All not-owned regions
+     ibset const notowned = domain_enlarged - allowned;
+     
+     // All not-active points
+     ibset const notactive = notowned.expand (buffer_width);
+     
+     // All active points
+     allactive = allowned - notactive;
+
+     // Fail if any of the consistency checks failed
+     assert(there_was_an_error == false);
+  }
+
+} // namespace CarpetIOHDF5
diff --git a/Carpet/CarpetIOHDF5/src/OutputSlice.cc b/Carpet/CarpetIOHDF5/src/OutputSlice.cc
index 3011696d2..5225e8aee 100644
--- a/Carpet/CarpetIOHDF5/src/OutputSlice.cc
+++ b/Carpet/CarpetIOHDF5/src/OutputSlice.cc
@@ -40,16 +40,17 @@ namespace CarpetIOHDF5 {
 
 
   // routines which are independent of the output dimension
-  static ibbox GetOutputBBox (const cGH* cctkGH,
+  static ibset GetOutputBBoxes (const cGH* cctkGH,
                               int group,
                               int rl, int m, int c,
-                              const ibbox& ext);
+                              const ibbox& ext, const ibset& allactive);
 
   static void GetCoordinates (const cGH* cctkGH, int m,
                               const cGroup& groupdata,
-                              const ibbox& ext,
+                              const ibset& exts,
                               CCTK_REAL& coord_time,
-                              rvect& coord_lower, rvect& coord_upper);
+                              vector<rvect>& coord_lower, 
+                              vector<rvect>& coord_upper);
 
   static int GetGridOffset (const cGH* cctkGH, int m, int dir,
                             const char* itempl, const char* iglobal,
@@ -60,6 +61,7 @@ namespace CarpetIOHDF5 {
 
 
 
+
   // IO processor
   template<int outdim> int IOHDF5<outdim>::ioproc;
   template<int outdim> int IOHDF5<outdim>::nioprocs;
@@ -557,6 +559,10 @@ namespace CarpetIOHDF5 {
 
       const gh* const hh = arrdata.at(group).at(m).hh;
       const dh* const dd = arrdata.at(group).at(m).dd;
+      
+      // re-compute the active (non-buffere) region
+      ibset allactive;
+      GetAllActive (dd, hh, m, reflevel, allactive);
 
       // Traverse all components on this multigrid level, refinement
       // level, and map
@@ -567,35 +573,42 @@ namespace CarpetIOHDF5 {
         groupdata.disttype != CCTK_DISTRIB_CONSTANT ?
         CCTK_nProcs(cctkGH) :
         1;
-      for (int c = c_min; c < c_max; ++ c) {
+      for (int c = c_min, c_base = c_min; c < c_max; ++ c) {
         int const lc = hh->get_local_component(rl,c);
         int const proc = hh->processor(rl,c);
+        const ibbox& data_ext = dd->light_boxes.at(ml).at(rl).at(c).exterior;
+        const ibset exts = GetOutputBBoxes (cctkGH, group, rl, m, c, 
+                                            data_ext, allactive);
         // we have to take part in this output if we either own data to be
         // output or are the ioproc in the same group of processors as the
         // dataholder
         if (dist::rank() == proc or 
             dist::rank() == IOProcForProc(proc)) {
           
-          const ibbox& data_ext = dd->light_boxes.at(ml).at(rl).at(c).exterior;
-          const ibbox ext = GetOutputBBox (cctkGH, group, rl, m, c, data_ext);
-          
           CCTK_REAL coord_time;
-          rvect coord_lower, coord_upper;
-          GetCoordinates (cctkGH, m, groupdata, ext,
+          vector<rvect> coord_lower, coord_upper;
+          GetCoordinates (cctkGH, m, groupdata, exts,
                           coord_time, coord_lower, coord_upper);
 
           // Apply offset
-          ivect offset1;
-          if (groupdata.grouptype == CCTK_GF) {
-            const ibbox& baseext = hh->baseextents.at(ml).at(rl);
-            offset1 = baseext.lower() + offset * ext.stride();
-          } else {
-            offset1 = offset * ext.stride();
-          }
-          for (int d=0; d<outdim; ++d) {
-            if (dirs[d] < 3) {
-              offset1[dirs[d]] = ext.lower()[dirs[d]];
+          vector<ivect> offsets1;
+          offsets1.reserve(exts.setsize());
+          for (ibset::const_iterator ext = exts.begin(); 
+               ext != exts.end(); 
+               ++ ext) {
+            ivect offset1;
+            if (groupdata.grouptype == CCTK_GF) {
+              const ibbox& baseext = hh->baseextents.at(ml).at(rl);
+              offset1 = baseext.lower() + offset * ext->stride();
+            } else {
+              offset1 = offset * ext->stride();
+            }
+            for (int d=0; d<outdim; ++d) {
+              if (dirs[d] < 3) {
+                offset1[dirs[d]] = ext->lower()[dirs[d]];
+              }
             }
+            offsets1.push_back(offset1);
           }
 
           const int tl_min = 0;
@@ -648,12 +661,18 @@ namespace CarpetIOHDF5 {
 
             }
 
-              error_count +=
-                WriteHDF5 (cctkGH, file, tmpdatas, ext, vindex,
-                           offset1, dirs,
-                           rl, ml, m, c, tl,
-                           coord_time, coord_lower, coord_upper);
             if (dist::rank() == IOProcForProc(proc)) {
+              int c_offset = 0;
+              for (ibset::const_iterator ext = exts.begin(); 
+                   ext != exts.end(); 
+                   ++ ext, ++ c_offset) {
+                error_count +=
+                  WriteHDF5 (cctkGH, file, tmpdatas, *ext, vindex,
+                             offsets1[c_offset], dirs,
+                             rl, ml, m, c, c_base + c_offset, tl,
+                             coord_time, coord_lower[c_offset], 
+                             coord_upper[c_offset]);
+              }
             }
 
             if (proc != ioproc) {
@@ -665,6 +684,8 @@ namespace CarpetIOHDF5 {
           } // for tl
 
         }
+
+        c_base += exts.setsize();
       } // for c
 
       error_count += CloseFile (cctkGH, file);
@@ -982,11 +1003,11 @@ namespace CarpetIOHDF5 {
 
 
 
-  // Omit symmetry and ghost zones if requested
-  ibbox GetOutputBBox (const cGH* const cctkGH,
+  // Omit symmetry, ghost and buffer zones if requested
+  ibset GetOutputBBoxes (const cGH* const cctkGH,
                        const int group,
                        const int rl, const int m, const int c,
-                       const ibbox& ext)
+                       const ibbox& ext, const ibset& allactive)
   {
     DECLARE_CCTK_PARAMETERS;
 
@@ -1044,7 +1065,13 @@ namespace CarpetIOHDF5 {
       }
     }
 
-    return ibbox(lo,hi,str);
+    ibset exts(ibbox(lo,hi,str));
+    // do grid arrays have buffer zones?
+    if (not output_buffer_points) {
+      exts &= allactive;
+    }
+
+    return exts;
   }
 
 
@@ -1052,9 +1079,9 @@ namespace CarpetIOHDF5 {
   // Determine coordinates
   void GetCoordinates (const cGH* const cctkGH, const int m,
                        const cGroup& groupdata,
-                       const ibbox& ext,
+                       const ibset& exts,
                        CCTK_REAL& coord_time,
-                       rvect& coord_lower, rvect& coord_upper)
+                       vector<rvect>& coord_lower, vector<rvect>& coord_upper)
   {
     coord_time = cctkGH->cctk_time;
 
@@ -1077,8 +1104,15 @@ namespace CarpetIOHDF5 {
       }
     }
 
-    coord_lower = global_lower + coord_delta * rvect(ext.lower());
-    coord_upper = global_lower + coord_delta * rvect(ext.upper());
+    coord_lower.reserve(exts.setsize());
+    coord_upper.reserve(exts.setsize());
+
+    for (ibset::const_iterator ext = exts.begin(); 
+         ext != exts.end(); 
+         ++ ext) {
+      coord_lower.push_back(global_lower + coord_delta * rvect(ext->lower()));
+      coord_upper.push_back(global_lower + coord_delta * rvect(ext->upper()));
+    }
   }
 
 
@@ -1192,6 +1226,7 @@ namespace CarpetIOHDF5 {
                                  const int ml,
                                  const int m,
                                  const int c,
+                                 const int output_component,
                                  const int tl,
                                  const CCTK_REAL coord_time,
                                  const vect<CCTK_REAL,dim>& coord_lower,
@@ -1263,7 +1298,7 @@ namespace CarpetIOHDF5 {
     }
     if (groupdata.grouptype == CCTK_GF or
         groupdata.disttype  != CCTK_DISTRIB_CONSTANT) {
-      datasetname_suffix << " c=" << c;
+      datasetname_suffix << " c=" << output_component;
     }
 
     // enable compression and checksums if requested
diff --git a/Carpet/CarpetIOHDF5/src/make.code.defn b/Carpet/CarpetIOHDF5/src/make.code.defn
index 223c95158..e8babc042 100644
--- a/Carpet/CarpetIOHDF5/src/make.code.defn
+++ b/Carpet/CarpetIOHDF5/src/make.code.defn
@@ -1,7 +1,7 @@
 # Main make.code.defn file for thorn CarpetIOHDF5
 
 # Source files in this directory
-SRCS = CarpetIOHDF5.cc Input.cc Output.cc OutputSlice.cc
+SRCS = CarpetIOHDF5.cc Input.cc Output.cc OutputSlice.cc GetAllActive.cc
 
 # Extend CXXFLAGS if HDF5 library was built with LFS support
 ifneq ($(strip $(HDF5_LFS_FLAGS)),)