CarpetIOHDF5: Add parameter "use_grid_structure_from_checkpoint"

Add a parameter "use_grid_structure_from_checkpoint" that reads the
grid structure from the checkpoint file, and sets up the Carpet grid
hierarchy accordingly.

The Carpet grid hierarchy is written unconditionally to all checkpoint
files.

darcs-hash:20060413202124-dae7b-f97e6aac2267ebc5f5e3867cbf78ca52bbd33016.gz

6 files changed, 138 insertions, 40 deletions

diff --git a/Carpet/CarpetIOHDF5/interface.ccl b/Carpet/CarpetIOHDF5/interface.ccl
index 100f4dcf6..82104d1a5 100644
--- a/Carpet/CarpetIOHDF5/interface.ccl
+++ b/Carpet/CarpetIOHDF5/interface.ccl
@@ -6,6 +6,8 @@ includes header: CarpetIOHDF5.hh in CarpetIOHDF5.hh
 
 uses include header: carpet.hh
 
+uses include header: defs.hh
+
 uses include header: bbox.hh
 uses include header: vect.hh
 
@@ -13,6 +15,7 @@ uses include header: data.hh
 uses include header: gdata.hh
 
 uses include header: ggf.hh
+uses include header: gh.hh
 
 
 
diff --git a/Carpet/CarpetIOHDF5/param.ccl b/Carpet/CarpetIOHDF5/param.ccl
index a190469a2..a62f28d2c 100644
--- a/Carpet/CarpetIOHDF5/param.ccl
+++ b/Carpet/CarpetIOHDF5/param.ccl
@@ -86,7 +86,10 @@ BOOLEAN checkpoint_next "Checkpoint at next iteration ?" STEERABLE = ALWAYS
 {
 } "no"
 
+BOOLEAN use_reflevels_from_checkpoint "Use 'CarpetRegrid::refinement_levels' from the checkpoint file rather than from the parameter file ?" STEERABLE = RECOVER
+{
+} "no"
 
-BOOLEAN use_reflevels_from_checkpoint "Use 'CarpetRegrid::refinement_levels' from the checkpoint file rather than from the parameter file ?" STEERABLE = ALWAYS
+BOOLEAN use_grid_structure_from_checkpoint "Use the grid structure stored in the checkpoint file" STEERABLE = RECOVER
 {
 } "no"
diff --git a/Carpet/CarpetIOHDF5/schedule.ccl b/Carpet/CarpetIOHDF5/schedule.ccl
index 0ff1986a5..b586376ad 100644
--- a/Carpet/CarpetIOHDF5/schedule.ccl
+++ b/Carpet/CarpetIOHDF5/schedule.ccl
@@ -57,6 +57,15 @@ if (! CCTK_Equals (recover, "no") && *recover_file)
     } "Overwrite 'CarpetRegird::refinement_levels' with the number of levels found in the checkpoint file"
   }
 
+  if (use_grid_structure_from_checkpoint)
+  {
+    schedule CarpetIOHDF5_RecoverGridStructure at RECOVER_PARAMETERS
+    {
+      LANG:C
+      OPTIONS: meta
+    } "Grid structure recovery routine"
+  }
+
   schedule CarpetIOHDF5_CloseFiles at POST_RECOVER_VARIABLES
   {
     LANG: C
diff --git a/Carpet/CarpetIOHDF5/src/CarpetIOHDF5.cc b/Carpet/CarpetIOHDF5/src/CarpetIOHDF5.cc
index 4392fa4e3..fc96024d6 100644
--- a/Carpet/CarpetIOHDF5/src/CarpetIOHDF5.cc
+++ b/Carpet/CarpetIOHDF5/src/CarpetIOHDF5.cc
@@ -1,7 +1,8 @@
-#include <assert.h>
-
+#include <cassert>
+#include <cstring>
 #include <map>
 #include <sstream>
+#include <string>
 
 #include "util_Table.h"
 #include "cctk.h"
@@ -14,6 +15,9 @@
 
 #include "CarpetIOHDF5.hh"
 
+#include "defs.hh"
+#include "gh.hh"
+
 
 namespace CarpetIOHDF5
 {
@@ -65,19 +69,17 @@ int CarpetIOHDF5_Startup (void)
 }
 
 
-int CarpetIOHDF5_Init (const cGH* const cctkGH)
+void CarpetIOHDF5_Init (CCTK_ARGUMENTS)
 {
   DECLARE_CCTK_ARGUMENTS;
 
   *this_iteration = -1;
   *next_output_iteration = 0;
   *next_output_time = cctk_time;
-
-  return (0);
 }
 
 
-int CarpetIOHDF5_InitialDataCheckpoint (const cGH* const cctkGH)
+void CarpetIOHDF5_InitialDataCheckpoint (CCTK_ARGUMENTS)
 {
   DECLARE_CCTK_PARAMETERS;
 
@@ -88,24 +90,23 @@ int CarpetIOHDF5_InitialDataCheckpoint (const cGH* const cctkGH)
     CCTK_INFO ("---------------------------------------------------------");
   }
   int retval = Checkpoint (cctkGH, CP_INITIAL_DATA);
-
-  return (retval);
 }
 
 
-int CarpetIOHDF5_EvolutionCheckpoint (const cGH* const cctkGH)
+void CarpetIOHDF5_EvolutionCheckpoint (CCTK_ARGUMENTS)
 {
   int retval = 0;
+  DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
 
 
   if (checkpoint and
-    ((checkpoint_every > 0 and cctkGH->cctk_iteration % checkpoint_every == 0) or
+    ((checkpoint_every > 0 and cctk_iteration % checkpoint_every == 0) or
      checkpoint_next)) {
     if (not CCTK_Equals (verbose, "none")) {
       CCTK_INFO ("---------------------------------------------------------");
       CCTK_VInfo (CCTK_THORNSTRING, "Dumping periodic checkpoint at "
-                  "iteration %d", cctkGH->cctk_iteration);
+                  "iteration %d", cctk_iteration);
       CCTK_INFO ("---------------------------------------------------------");
     }
 
@@ -115,27 +116,26 @@ int CarpetIOHDF5_EvolutionCheckpoint (const cGH* const cctkGH)
       CCTK_ParameterSet ("checkpoint_next", CCTK_THORNSTRING, "no");
     }
   }
-
-  return (retval);
 }
 
 
-int CarpetIOHDF5_TerminationCheckpoint (const cGH *const GH)
+void CarpetIOHDF5_TerminationCheckpoint (CCTK_ARGUMENTS)
 {
   int retval = 0;
+  DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
 
 
   if (checkpoint and checkpoint_on_terminate) {
-    if (last_checkpoint_iteration < GH->cctk_iteration) {
+    if (last_checkpoint_iteration < cctk_iteration) {
       if (not CCTK_Equals (verbose, "none")) {
         CCTK_INFO ("---------------------------------------------------------");
         CCTK_VInfo (CCTK_THORNSTRING, "Dumping termination checkpoint at "
-                    "iteration %d", GH->cctk_iteration);
+                    "iteration %d", cctk_iteration);
         CCTK_INFO ("---------------------------------------------------------");
       }
 
-      retval = Checkpoint (GH, CP_EVOLUTION_DATA);
+      retval = Checkpoint (cctkGH, CP_EVOLUTION_DATA);
     } else if (not CCTK_Equals (verbose, "none")) {
       CCTK_INFO ("---------------------------------------------------------");
       CCTK_VInfo (CCTK_THORNSTRING, "Termination checkpoint already dumped "
@@ -144,8 +144,6 @@ int CarpetIOHDF5_TerminationCheckpoint (const cGH *const GH)
       CCTK_INFO ("---------------------------------------------------------");
     }
   }
-
-  return (retval);
 }
 
 
@@ -957,6 +955,27 @@ static int WriteMetadata (const cGH *cctkGH, int nioprocs,
     free (parameters);
   }
 
+  // Save grid structure
+  if (called_from_checkpoint) {
+    vector <grid_structure_t> grid_structure (maps);
+    for (int m = 0; m < maps; ++ m) {
+      grid_structure.at(m).bbss = Carpet::vhh.at(m)->extents().at(0);
+      grid_structure.at(m).obss = Carpet::vhh.at(m)->outer_boundaries();
+    }
+    ostringstream gs_buf;
+    gs_buf << grid_structure;
+    string const gs_str = gs_buf.str();
+    size = gs_str.size() + 1;
+    char const * const gs_cstr = gs_str.c_str();
+    HDF5_ERROR (H5Sset_extent_simple (array_dataspace, 1, & size, NULL));
+    hid_t dataset;
+    HDF5_ERROR (dataset = H5Dcreate (group, GRID_STRUCTURE, H5T_NATIVE_UCHAR,
+                                     array_dataspace, H5P_DEFAULT));
+    HDF5_ERROR (H5Dwrite (dataset, H5T_NATIVE_CHAR, H5S_ALL, H5S_ALL,
+                          H5P_DEFAULT, gs_cstr));
+    HDF5_ERROR (H5Dclose (dataset));
+  }
+
   HDF5_ERROR (H5Tclose (datatype));
   HDF5_ERROR (H5Sclose (scalar_dataspace));
   HDF5_ERROR (H5Sclose (array_dataspace));
diff --git a/Carpet/CarpetIOHDF5/src/CarpetIOHDF5.hh b/Carpet/CarpetIOHDF5/src/CarpetIOHDF5.hh
index ab5624149..1101a2545 100644
--- a/Carpet/CarpetIOHDF5/src/CarpetIOHDF5.hh
+++ b/Carpet/CarpetIOHDF5/src/CarpetIOHDF5.hh
@@ -3,16 +3,18 @@
 
 
 #include <hdf5.h>
-#include "carpet.hh"
+#include "cctk_Arguments.h"
 #include "CactusBase/IOUtil/src/ioutil_Utils.h"
+#include "carpet.hh"
 
 
 // some macros for HDF5 group names
 #define METADATA_GROUP "Parameters and Global Attributes"
 #define ALL_PARAMETERS "All Parameters"
+#define GRID_STRUCTURE "Grid Structure"
 
 // atomic HDF5 datatypes for the generic CCTK datatypes
-// (the one for CCTK_COMPLEX is created at startup as a compound HDF5 datatype
+// (the one for CCTK_COMPLEX is created at startup as a compound HDF5 datatype)
 #define HDF5_CHAR   H5T_NATIVE_CHAR
 
 #ifdef  CCTK_REAL_PRECISION_16
@@ -110,13 +112,14 @@ namespace CarpetIOHDF5
   extern "C" {
 
     int CarpetIOHDF5_Startup (void);
-    int CarpetIOHDF5_RecoverParameters (void);
-    int CarpetIOHDF5_Init (const cGH* const);
-    int CarpetIOHDF5_SetNumRefinementLevels (void);
-    int CarpetIOHDF5_InitialDataCheckpoint (const cGH* const);
-    int CarpetIOHDF5_EvolutionCheckpoint (const cGH* const);
-    int CarpetIOHDF5_TerminationCheckpoint (const cGH* const);
-    void CarpetIOHDF5_CloseFiles (const cGH* const);
+    void CarpetIOHDF5_Init (CCTK_ARGUMENTS);
+    void CarpetIOHDF5_SetNumRefinementLevels (void);
+    void CarpetIOHDF5_CloseFiles (CCTK_ARGUMENTS);
+    void CarpetIOHDF5_InitialDataCheckpoint (CCTK_ARGUMENTS);
+    void CarpetIOHDF5_EvolutionCheckpoint (CCTK_ARGUMENTS);
+    void CarpetIOHDF5_TerminationCheckpoint (CCTK_ARGUMENTS);
+    void CarpetIOHDF5_RecoverParameters (void);
+    void CarpetIOHDF5_RecoverGridStructure (void);
 
   } // extern "C"
 
diff --git a/Carpet/CarpetIOHDF5/src/Input.cc b/Carpet/CarpetIOHDF5/src/Input.cc
index 098a872fa..225fc9b55 100644
--- a/Carpet/CarpetIOHDF5/src/Input.cc
+++ b/Carpet/CarpetIOHDF5/src/Input.cc
@@ -1,5 +1,7 @@
-#include <assert.h>
-#include <string.h>
+#include <cassert>
+#include <cstring>
+#include <sstream>
+#include <vector>
 
 #include "util_Table.h"
 #include "cctk.h"
@@ -9,6 +11,8 @@
 #include "CactusBase/IOUtil/src/ioGH.h"
 #include "CactusBase/IOUtil/src/ioutil_CheckpointRecovery.h"
 
+#include "defs.hh"
+
 
 namespace CarpetIOHDF5
 {
@@ -53,6 +57,8 @@ typedef struct {
   CCTK_REAL global_time;
   CCTK_REAL delta_time;
   vector<CCTK_REAL> mgleveltimes;  // [num_mglevels*num_reflevels]
+
+  vector<grid_structure_t> grid_structure; // [maps]
 } fileset_t;
 
 // list of checkpoint/filereader files
@@ -78,13 +84,55 @@ static int ReadVar (const cGH* const cctkGH,
 //////////////////////////////////////////////////////////////////////////////
 // Register with the Cactus Recovery Interface
 //////////////////////////////////////////////////////////////////////////////
-int CarpetIOHDF5_RecoverParameters (void)
+void CarpetIOHDF5_RecoverParameters ()
 {
-  int retval = IOUtil_RecoverParameters (Recover, ".h5", "HDF5");
-
-  return (retval);
+  IOUtil_RecoverParameters (Recover, ".h5", "HDF5");
 }
 
+//////////////////////////////////////////////////////////////////////////////
+// Recover the grid structure
+//////////////////////////////////////////////////////////////////////////////
+void CarpetIOHDF5_RecoverGridStructure ()
+{
+  cGH const * const cctkGH = 0; // fake it
+  
+  fileset_t & fileset = * filesets.begin();
+  
+  // Abort with an error if there is no grid structure in the
+  // checkpoint file
+  assert (fileset.grid_structure.size() == maps);
+  
+  for (int m = 0; m < maps; ++ m) {
+    grid_structure_t const & grid_structure = fileset.grid_structure.at(m);
+    
+    int const rls = grid_structure.bbss.size();
+    assert (grid_structure.obss.size() == rls);
+    
+    vector <vector <ibbox> >  bbss = grid_structure.bbss;
+    vector <vector <bbvect> > obss = grid_structure.obss;
+    vector <vector <int> >    pss (rls);
+    
+    for (int rl = 0; rl < rls; ++ rl) {
+      
+      vector <ibbox>  & bbs = bbss.at(rl);
+      vector <bbvect> & obs = obss.at(rl);
+      vector <int>    & ps  = pss.at(rl);
+      
+      // Make multiprocessor aware
+      Carpet::SplitRegions (cctkGH, bbs, obs, ps);
+      
+    } // for rl
+    
+    // Make multigrid aware
+    vector <vector <vector <ibbox> > > bbsss;
+    Carpet::MakeMultigridBoxes (cctkGH, bbss, obss, bbsss);
+    
+    // Regrid
+    Carpet::vhh.at(m)->recompose (bbsss, obss, pss, false);
+    Carpet::OutputGrids (cctkGH, m, * Carpet::vhh.at(m));
+    
+  } // for m
+}
 
 //////////////////////////////////////////////////////////////////////////////
 // Overwrite the "CarpetRegrid::refinement_levels"
@@ -93,7 +141,7 @@ int CarpetIOHDF5_RecoverParameters (void)
 // Note that this has to be done after parameter recovery in order to have
 // any effect of steering "CarpetRegrid::refinement_levels".
 //////////////////////////////////////////////////////////////////////////////
-int CarpetIOHDF5_SetNumRefinementLevels (void)
+void CarpetIOHDF5_SetNumRefinementLevels ()
 {
   DECLARE_CCTK_PARAMETERS;
 
@@ -114,18 +162,17 @@ int CarpetIOHDF5_SetNumRefinementLevels (void)
                                           buffer);
     assert (retval == 0);
   }
-
-  return (0);
 }
 
 
 //////////////////////////////////////////////////////////////////////////////
 // close all open checkpoint/filereader files after recovering grid variables
 //////////////////////////////////////////////////////////////////////////////
-void CarpetIOHDF5_CloseFiles (const cGH* const cctkGH)
+void CarpetIOHDF5_CloseFiles (CCTK_ARGUMENTS)
 {
-  int error_count = 0;
+  DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
+  int error_count = 0;
 
 
   for (list<fileset_t>::const_iterator set = filesets.begin();
@@ -586,6 +633,20 @@ static void ReadMetadata (fileset_t& fileset, hid_t file)
   HDF5_ERROR (H5Aread (attr, H5T_NATIVE_DOUBLE, &fileset.delta_time));
   HDF5_ERROR (H5Aclose (attr));
 
+  // Read grid structure if it is present
+  hid_t dataset;
+  H5E_BEGIN_TRY {
+    dataset = H5Dopen (metadata, METADATA_GROUP "/" GRID_STRUCTURE);
+  } H5E_END_TRY;
+  if (dataset >= 0) {
+    vector<char> gs_cstr (H5Dget_storage_size (dataset) + 1);
+    HDF5_ERROR (H5Dread (dataset, H5T_NATIVE_CHAR, H5S_ALL, H5S_ALL,
+                         H5P_DEFAULT, &gs_cstr.front()));
+    HDF5_ERROR (H5Dclose (dataset));
+    istringstream gs_buf (&gs_cstr.front());
+    gs_buf >> fileset.grid_structure;
+  }
+
   fileset.mgleveltimes.resize (fileset.num_mglevels * fileset.num_reflevels);
   for (int i = 0; i < fileset.num_mglevels; i++) {
     char buffer[32];