CarpetIOHDF5: Allow different numbers of ghost zones on different levels

Allow different numbers of ghost zones and different spatial
prolongation orders on different refinement levels.

This causes incompatible changes to the checkpoint file format.

4 files changed, 94 insertions, 10 deletions

diff --git a/Carpet/CarpetIOHDF5/src/CarpetIOHDF5.cc b/Carpet/CarpetIOHDF5/src/CarpetIOHDF5.cc
index a139cea7e..01888001a 100644
--- a/Carpet/CarpetIOHDF5/src/CarpetIOHDF5.cc
+++ b/Carpet/CarpetIOHDF5/src/CarpetIOHDF5.cc
@@ -1170,10 +1170,16 @@ int WriteMetadata (const cGH * const cctkGH, int const nioprocs,
     vector <vector <vector <region_t> > > grid_superstructure (maps);
     vector <vector <vector <region_t> > > grid_structure (maps);
     vector <vector <vector <CCTK_REAL> > > grid_times (maps);
+    vector <vector <i2vect> > grid_ghosts (maps);
+    vector <vector <i2vect> > grid_buffers (maps);
+    vector <vector <int> > grid_prolongation_orders (maps);
     for (int m = 0; m < maps; ++ m) {
       grid_superstructure.at(m) = vhh.at(m)->superregions;
       grid_structure.at(m) = vhh.at(m)->regions.at(0);
       grid_times.at(m).resize(mglevels);
+      grid_ghosts.at(m) = vdd.at(m)->ghost_widths;
+      grid_buffers.at(m) = vdd.at(m)->buffer_widths;
+      grid_prolongation_orders.at(m) = vdd.at(m)->prolongation_orders_space;
       for  (int ml = 0; ml < mglevels; ++ ml) {
         grid_times.at(m).at(ml).resize(vhh.at(m)->reflevels());
         for (int rl = 0; rl < vhh.at(m)->reflevels(); ++ rl) {
@@ -1186,12 +1192,15 @@ int WriteMetadata (const cGH * const cctkGH, int const nioprocs,
     // We could write this information only into one of the checkpoint
     // files (to save space), or write it into a separate metadata
     // file
-    gs_buf << grid_superstructure;
+    gs_buf << "grid_superstructure:" << grid_superstructure << ",";
     // We could omit the grid structure (to save space), or write it
     // only into one of the checkpoint files
-    gs_buf << grid_structure;
-    gs_buf << grid_times;
-    gs_buf << leveltimes;
+    gs_buf << "grid_structure:" << grid_structure << ",";
+    gs_buf << "grid_times:" << grid_times << ",";
+    gs_buf << "leveltimes:" << leveltimes << ",";
+    gs_buf << "grid_ghosts" << grid_ghosts << ",";
+    gs_buf << "grid_buffers" << grid_buffers << ",";
+    gs_buf << "grid_prolongation_orders" << grid_prolongation_orders << ".";
     string const gs_str = gs_buf.str();
     error_count += WriteLargeAttribute (group, GRID_STRUCTURE, gs_str.c_str());
   }
diff --git a/Carpet/CarpetIOHDF5/src/CarpetIOHDF5.hh b/Carpet/CarpetIOHDF5/src/CarpetIOHDF5.hh
index f31be614c..81ae3d943 100644
--- a/Carpet/CarpetIOHDF5/src/CarpetIOHDF5.hh
+++ b/Carpet/CarpetIOHDF5/src/CarpetIOHDF5.hh
@@ -17,7 +17,7 @@
 // some macros for HDF5 group names
 #define METADATA_GROUP "Parameters and Global Attributes"
 #define ALL_PARAMETERS "All Parameters"
-#define GRID_STRUCTURE "Grid Structure v2"
+#define GRID_STRUCTURE "Grid Structure v3"
 
 // atomic HDF5 datatypes for the generic CCTK datatypes
 // (the one for CCTK_COMPLEX is created at startup as a compound HDF5 datatype)
diff --git a/Carpet/CarpetIOHDF5/src/Input.cc b/Carpet/CarpetIOHDF5/src/Input.cc
index cc3358840..2414bc3ca 100644
--- a/Carpet/CarpetIOHDF5/src/Input.cc
+++ b/Carpet/CarpetIOHDF5/src/Input.cc
@@ -65,6 +65,9 @@ typedef struct {
   vector<vector<vector<region_t> > > grid_structure; // [map][reflevel][component]
   vector<vector<vector<CCTK_REAL> > > grid_times;    // [map][mglevel][reflevel]
   vector<vector<CCTK_REAL> > leveltimes;             // [mglevel][reflevel]
+  vector <vector <i2vect> > grid_ghosts; // [map]
+  vector <vector <i2vect> > grid_buffers; // [map]
+  vector <vector <int> > grid_prolongation_orders; // [map]
 
 } fileset_t;
 
@@ -212,6 +215,33 @@ void CarpetIOHDF5_RecoverGridStructure (CCTK_ARGUMENTS)
   }
   
   RegridFree (cctkGH, false);
+
+  // Check ghost and buffer widths and prolongation orders
+  // TODO: Instead of only checking them, set them during regridding.
+  // (This requires setting them during regridding instead of during setup.)
+  for (int m = 0; m < maps; ++ m) {
+    assert (fileset.grid_ghosts.at(m).size()
+            == vdd.at(m)->ghost_widths.size());
+    for (size_t rl = 0; rl < fileset.grid_ghosts.at(m).size(); ++ rl) {
+      assert (all (all (fileset.grid_ghosts.at(m).at(rl)
+                        == vdd.at(m)->ghost_widths.at(rl))));
+    }
+    assert (fileset.grid_buffers.at(m).size()
+            == vdd.at(m)->buffer_widths.size());
+    for (size_t rl = 0; rl < fileset.grid_buffers.at(m).size(); ++ rl) {
+      assert (all (all (fileset.grid_buffers.at(m).at(rl)
+                        == vdd.at(m)->buffer_widths.at(rl))));
+    }
+    assert (fileset.grid_prolongation_orders.at(m).size()
+            == vdd.at(m)->prolongation_orders_space.size());
+    for (size_t rl = 0;
+         rl < fileset.grid_prolongation_orders.at(m).size();
+         ++ rl)
+    {
+      assert (fileset.grid_prolongation_orders.at(m).at(rl)
+              == vdd.at(m)->prolongation_orders_space.at(rl));
+    }
+  }
 }
 
 //////////////////////////////////////////////////////////////////////////////
@@ -882,10 +912,55 @@ static void ReadMetadata (fileset_t& fileset, hid_t file)
                          H5P_DEFAULT, &gs_cstr.front()));
     HDF5_ERROR (H5Dclose (dataset));
     istringstream gs_buf (&gs_cstr.front());
+    
+    skipws (gs_buf);
+    consume (gs_buf, "grid_superstructure:");
+    skipws (gs_buf);
     gs_buf >> fileset.grid_superstructure;
+    skipws (gs_buf);
+    consume (gs_buf, ",");
+    
+    skipws (gs_buf);
+    consume (gs_buf, "grid_structure:");
+    skipws (gs_buf);
     gs_buf >> fileset.grid_structure;
+    skipws (gs_buf);
+    consume (gs_buf, ",");
+    
+    skipws (gs_buf);
+    consume (gs_buf, "grid_times:");
+    skipws (gs_buf);
     gs_buf >> fileset.grid_times;
+    skipws (gs_buf);
+    consume (gs_buf, ",");
+    
+    skipws (gs_buf);
+    consume (gs_buf, "grid_leveltimes:");
+    skipws (gs_buf);
     gs_buf >> fileset.leveltimes;
+    skipws (gs_buf);
+    consume (gs_buf, ",");
+    
+    skipws (gs_buf);
+    consume (gs_buf, "grid_ghosts:");
+    skipws (gs_buf);
+    gs_buf >> fileset.grid_ghosts;
+    skipws (gs_buf);
+    consume (gs_buf, ",");
+    
+    skipws (gs_buf);
+    consume (gs_buf, "grid_buffers:");
+    skipws (gs_buf);
+    gs_buf >> fileset.grid_buffers;
+    skipws (gs_buf);
+    consume (gs_buf, ",");
+    
+    skipws (gs_buf);
+    consume (gs_buf, "grid_prolongation_orders:");
+    skipws (gs_buf);
+    gs_buf >> fileset.grid_prolongation_orders;
+    skipws (gs_buf);
+    consume (gs_buf, ".");
   }
 
   fileset.mgleveltimes.resize (fileset.num_mglevels * fileset.num_reflevels);
diff --git a/Carpet/CarpetIOHDF5/src/OutputSlice.cc b/Carpet/CarpetIOHDF5/src/OutputSlice.cc
index b1698b2a1..712844b3e 100644
--- a/Carpet/CarpetIOHDF5/src/OutputSlice.cc
+++ b/Carpet/CarpetIOHDF5/src/OutputSlice.cc
@@ -39,7 +39,7 @@ namespace CarpetIOHDF5 {
   // routines which are independent of the output dimension
   static ibbox GetOutputBBox (const cGH* cctkGH,
                               int group,
-                              int m, int c,
+                              int rl, int m, int c,
                               const ibbox& ext);
 
   static void GetCoordinates (const cGH* cctkGH, int m,
@@ -548,7 +548,7 @@ namespace CarpetIOHDF5 {
         if (dist::rank() == proc or dist::rank() == ioproc) {
           
           const ibbox& data_ext = dd->boxes.at(ml).at(rl).at(c).exterior;
-          const ibbox ext = GetOutputBBox (cctkGH, group, m, c, data_ext);
+          const ibbox ext = GetOutputBBox (cctkGH, group, rl, m, c, data_ext);
           
           CCTK_REAL coord_time;
           rvect coord_lower, coord_upper;
@@ -950,7 +950,7 @@ namespace CarpetIOHDF5 {
   // Omit symmetry and ghost zones if requested
   ibbox GetOutputBBox (const cGH* const cctkGH,
                        const int group,
-                       const int m, const int c,
+                       const int rl, const int m, const int c,
                        const ibbox& ext)
   {
     DECLARE_CCTK_PARAMETERS;
@@ -984,8 +984,8 @@ namespace CarpetIOHDF5 {
     ivect hi = ext.upper();
     const ivect str = ext.stride();
 
-    const b2vect obnds       = vhh.at(m)->outer_boundaries(reflevel,c);
-    const i2vect ghost_width = arrdata.at(group).at(m).dd->ghost_width;
+    const b2vect obnds       = vhh.at(m)->outer_boundaries(rl,c);
+    const i2vect ghost_width = arrdata.at(group).at(m).dd->ghost_widths.AT(rl);
 
     for (int d=0; d<groupdim; ++d) {
       bool const output_lower_ghosts =