CarpetIOHDF5: output index files for sliced output if requested

3 files changed, 84 insertions, 13 deletions

diff --git a/Carpet/CarpetIOHDF5/src/CarpetIOHDF5.cc b/Carpet/CarpetIOHDF5/src/CarpetIOHDF5.cc
index 0ce40eee4..40a2a5d56 100644
--- a/Carpet/CarpetIOHDF5/src/CarpetIOHDF5.cc
+++ b/Carpet/CarpetIOHDF5/src/CarpetIOHDF5.cc
@@ -82,6 +82,10 @@ static int WriteAttribute (hid_t const group,
                            char const * const name,
                            char const * const * const svalues,
                            int const nvalues);
+static int WriteAttribute (hid_t const group,
+                           char const * const name,
+                           hsize_t const * const svalues,
+                           int const nvalues);
 static int WriteLargeAttribute (hid_t const group,
                                 char const * const name,
                                 char const * const svalue);
@@ -333,7 +337,9 @@ int AddSliceAttributes(const cGH* const cctkGH,
                        const vector<int>& ioffsetdenom,
                        const vector<int>& bbox,
                        const vector<int>& nghostzones,
-                       hid_t& dataset)
+                       hid_t& dataset,
+                       const vector<hsize_t>& shape,
+                       const bool is_index)
 {
   int error_count = 0;
 
@@ -360,6 +366,10 @@ int AddSliceAttributes(const cGH* const cctkGH,
     int const map_is_cartesian = MultiPatch_MapIsCartesian (map);
     error_count += WriteAttribute(dataset, "MapIsCartesian", map_is_cartesian);
   }
+  if (is_index) {
+    error_count += WriteAttribute(dataset, "h5shape", &shape[0], shape.size());
+  }
+
 
   return error_count;
 }
@@ -1473,6 +1483,27 @@ static int WriteAttribute (hid_t const group,
 }
 
 
+// Write an array of hsize_t attributes
+static int WriteAttribute (hid_t const group,
+                           char const * const name,
+                           hsize_t const * const svalues,
+                           int const nvalues)
+{
+  hid_t dataspace, attribute;
+  int error_count = 0;
+  
+  hsize_t const size = nvalues;
+  HDF5_ERROR (dataspace = H5Screate_simple (1, & size, NULL));
+  HDF5_ERROR (attribute = H5Acreate (group, name, H5T_NATIVE_HSIZE,
+                                     dataspace, H5P_DEFAULT));
+  HDF5_ERROR (H5Awrite (attribute, H5T_NATIVE_HSIZE, svalues));
+  HDF5_ERROR (H5Aclose (attribute));
+  HDF5_ERROR (H5Sclose (dataspace));
+  
+  return error_count;
+}
+
+
 // Write an array of string attributes
 static int WriteAttribute (hid_t const group,
                            char const * const name,
diff --git a/Carpet/CarpetIOHDF5/src/CarpetIOHDF5.hh b/Carpet/CarpetIOHDF5/src/CarpetIOHDF5.hh
index 42c7d6622..31d97bc9a 100644
--- a/Carpet/CarpetIOHDF5/src/CarpetIOHDF5.hh
+++ b/Carpet/CarpetIOHDF5/src/CarpetIOHDF5.hh
@@ -145,7 +145,9 @@ namespace CarpetIOHDF5
                          const vector<int>& ioffsetdenom,
                          const vector<int>& bbox,
                          const vector<int>& nghostzones,
-                         hid_t& dataset);
+                         hid_t& dataset,
+                         const vector<hsize_t>& shape,
+                         const bool is_index);
 
   // returns an HDF5 datatype corresponding to the given CCTK datatype
   hid_t CCTKtoHDF5_Datatype (const cGH* const cctkGH,
@@ -220,10 +222,12 @@ namespace CarpetIOHDF5
                          const vect<int,outdim>& dirs,
                          bool is_new_file,
                          bool truncate_file,
-                         hid_t& file);
+                         hid_t& file,
+                         hid_t& index_file);
 
      static int WriteHDF5 (const cGH* cctkGH,
                            hid_t& file,
+                           hid_t& index_file,
                            vector<gdata*> const gfdatas,
                            const bbox<int,dim>& gfext,
                            const int vi,
@@ -240,7 +244,8 @@ namespace CarpetIOHDF5
                            const vect<CCTK_REAL,dim>& coord_upper);
 
     static int CloseFile (const cGH* cctkGH,
-                          hid_t& file);
+                          hid_t& file,
+                          hid_t& index_file);
 
     static ivect GetOutputOffset (const cGH* cctkGH, int m,
                                   const vect<int,outdim>& dirs);
diff --git a/Carpet/CarpetIOHDF5/src/OutputSlice.cc b/Carpet/CarpetIOHDF5/src/OutputSlice.cc
index 4f1620c0e..acc94493c 100644
--- a/Carpet/CarpetIOHDF5/src/OutputSlice.cc
+++ b/Carpet/CarpetIOHDF5/src/OutputSlice.cc
@@ -547,10 +547,10 @@ namespace CarpetIOHDF5 {
     const int m_max = groupdata.grouptype == CCTK_GF ? Carpet::maps : 1;
     for (int m = m_min; m < m_max; ++ m) {
 
-      hid_t file = -1;
+      hid_t file = -1, index_file = -1;
       int error_count = 0;
       error_count += OpenFile (cctkGH, m, vindex, numvars, alias, basefilename,
-                               dirs, is_new_file, truncate_file, file);
+                               dirs, is_new_file, truncate_file, file, index_file);
 
       // Find the output offset
       const ivect offset =
@@ -667,7 +667,7 @@ namespace CarpetIOHDF5 {
                    ext != exts.end(); 
                    ++ ext, ++ c_offset) {
                 error_count +=
-                  WriteHDF5 (cctkGH, file, tmpdatas, *ext, vindex,
+                  WriteHDF5 (cctkGH, file, index_file, tmpdatas, *ext, vindex,
                              offsets1[c_offset], dirs,
                              rl, ml, m, c, c_base + c_offset, tl,
                              coord_time, coord_lower[c_offset], 
@@ -688,7 +688,7 @@ namespace CarpetIOHDF5 {
         c_base += exts.setsize();
       } // for c
 
-      error_count += CloseFile (cctkGH, file);
+      error_count += CloseFile (cctkGH, file, index_file);
       if (error_count > 0 and abort_on_io_errors) {
         CCTK_WARN (0, "Aborting simulation due to previous I/O errors");
       }
@@ -769,7 +769,8 @@ namespace CarpetIOHDF5 {
                                 const vect<int,outdim>& dirs,
                                 const bool is_new_file,
                                 const bool truncate_file,
-                                hid_t& file)
+                                hid_t& file,
+                                hid_t& index_file)
   {
     DECLARE_CCTK_PARAMETERS;
 
@@ -800,7 +801,9 @@ namespace CarpetIOHDF5 {
       if (nioprocs > 1) {
         filenamebuf << ".file_" << dist::rank();
       }
-      filenamebuf << ".h5";
+      string index_filename(filenamebuf.str() + ".idx" + out_extension);
+      filenamebuf << out_extension;
+      
       // we need a persistent temporary here
       const string filenamestr = filenamebuf.str();
       const char* const filename = filenamestr.c_str();
@@ -816,11 +819,23 @@ namespace CarpetIOHDF5 {
       if (truncate_file or not file_exists) {
         HDF5_ERROR(file = H5Fcreate(filename, H5F_ACC_TRUNC, H5P_DEFAULT,
                                     H5P_DEFAULT));
+        if (output_index) {
+          HDF5_ERROR (index_file = H5Fcreate (index_filename.c_str(),
+                                              H5F_ACC_TRUNC, H5P_DEFAULT,
+                                              H5P_DEFAULT));
+        }
         // write metadata information
         error_count +=
           WriteMetadata(cctkGH, nioprocs, vindex, numvars, false, file);
+
+        if (output_index) {
+          error_count +=
+            WriteMetadata (cctkGH, nioprocs, vindex, numvars, false, index_file);
+        }
       } else {
         HDF5_ERROR (file = H5Fopen(filename, H5F_ACC_RDWR, H5P_DEFAULT));
+        if (output_index)
+          HDF5_ERROR (index_file = H5Fopen (index_filename.c_str(), H5F_ACC_RDWR, H5P_DEFAULT));
       }
       io_files += 1;
 
@@ -831,7 +846,7 @@ namespace CarpetIOHDF5 {
 
   template<int outdim>
   int IOHDF5<outdim>::CloseFile (const cGH* const cctkGH,
-                                 hid_t& file)
+                                 hid_t& file, hid_t& index_file)
   {
     DECLARE_CCTK_PARAMETERS;
 
@@ -841,6 +856,9 @@ namespace CarpetIOHDF5 {
       if (file >= 0) {
         HDF5_ERROR(H5Fclose(file));
       }
+      if (output_index and index_file >= 0) {
+        HDF5_ERROR (H5Fclose (index_file));
+      }
     }
     if (nioprocs > 1) {
       CCTK_REAL local[2], global[2];
@@ -1242,6 +1260,7 @@ namespace CarpetIOHDF5 {
   template<int outdim>
   int IOHDF5<outdim>::WriteHDF5 (const cGH* cctkGH,
                                  hid_t& file,
+                                 hid_t& indexfile,
                                  vector<gdata*> const gfdatas,
                                  const bbox<int,dim>& gfext,
                                  const int vi,
@@ -1457,6 +1476,8 @@ namespace CarpetIOHDF5 {
         if (request->check_exist) {
           H5E_BEGIN_TRY {
             H5Gunlink(file, datasetname.c_str());
+            if (indexfile != -1)
+              H5Gunlink (indexfile, datasetname.c_str());
           } H5E_END_TRY;
         }
         // free I/O request structure
@@ -1465,16 +1486,30 @@ namespace CarpetIOHDF5 {
         }
 
         // write the dataset
-        hid_t dataset;
+        hid_t dataset, index_dataset;
         HDF5_ERROR(dataset =
                    H5Dcreate (file, datasetname.c_str(),
                               slice_type, slice_space, plist));
+        if (indexfile != -1) {
+          HDF5_ERROR (index_dataset = H5Dcreate (indexfile, datasetname.c_str(),
+                                                 slice_type, slice_space, plist));
+        }
+
         HDF5_ERROR(H5Dwrite (dataset, mem_type, mem_space, H5S_ALL,
                              H5P_DEFAULT, gfdatas[n]->storage()));
         error_count +=
           AddSliceAttributes (cctkGH, fullname, rl, ml, m, tl, origin, delta, 
-                              iorigin, ioffset, ioffsetdenom, bbox, nghostzones, dataset);
+                              iorigin, ioffset, ioffsetdenom, bbox, nghostzones,
+                              dataset, slice_shape, false);
         HDF5_ERROR(H5Dclose (dataset));
+
+        if (indexfile != -1) {
+          error_count += 
+            AddSliceAttributes (cctkGH, fullname, rl, ml, m, tl, origin, delta, 
+                                iorigin, ioffset, ioffsetdenom, bbox, nghostzones,
+                                index_dataset, slice_shape, true);
+          HDF5_ERROR (H5Dclose (index_dataset));
+        }
         free (fullname);
 
         io_bytes +=