CarpetIOF5: Various improvements

13 files changed, 489 insertions, 121 deletions

diff --git a/CarpetDev/CarpetIOF5/interface.ccl b/CarpetDev/CarpetIOF5/interface.ccl
index de552a6b2..59cea71d5 100644
--- a/CarpetDev/CarpetIOF5/interface.ccl
+++ b/CarpetDev/CarpetIOF5/interface.ccl
@@ -4,6 +4,7 @@ IMPLEMENTS: CarpetIOF5
 
 INHERITS: grid
 
+USES INCLUDE: cacheinfo.hh
 USES INCLUDE: carpet.hh
 USES INCLUDE: CarpetTimers.hh
 
diff --git a/CarpetDev/CarpetIOF5/par/iof5-array.par b/CarpetDev/CarpetIOF5/par/iof5-array.par
index 8c9fd2a03..26de7054d 100644
--- a/CarpetDev/CarpetIOF5/par/iof5-array.par
+++ b/CarpetDev/CarpetIOF5/par/iof5-array.par
@@ -29,8 +29,6 @@ Cactus::cctk_itlast = 2
 
 IO::out_dir = $parfile
 
-Carpet::processor_topology = "balanced"
-
 Carpet::domain_from_coordbase = yes
 CartGrid3D::type              = "coordbase"
 CoordBase::domainsize         = "minmax"
@@ -40,9 +38,9 @@ CoordBase::zmin               = -1.0
 CoordBase::xmax               = +1.0
 CoordBase::ymax               = +1.0
 CoordBase::zmax               = +1.0
-CoordBase::dx                 =  0.25
-CoordBase::dy                 =  0.25
-CoordBase::dz                 =  0.25
+CoordBase::dx                 =  0.025
+CoordBase::dy                 =  0.025
+CoordBase::dz                 =  0.025
 
 InterpToArray::nscalars           = 1
 InterpToArray::scalar_vars[0]     = "grid::r"
diff --git a/CarpetDev/CarpetIOF5/par/iof5-multipatch-kerrschild-input.par b/CarpetDev/CarpetIOF5/par/iof5-multipatch-kerrschild-input.par
index c28ade8a0..d6c5cbe5c 100644
--- a/CarpetDev/CarpetIOF5/par/iof5-multipatch-kerrschild-input.par
+++ b/CarpetDev/CarpetIOF5/par/iof5-multipatch-kerrschild-input.par
@@ -145,6 +145,8 @@ CarpetIOASCII::out1D_vars         = "
 
 
 
+CarpetIOF5::verbose = yes
+
 CarpetIOF5::out_every = 1024
 CarpetIOF5::out_vars  = "
         ADMBase::metric
diff --git a/CarpetDev/CarpetIOF5/par/iof5-multipatch-kerrschild.par b/CarpetDev/CarpetIOF5/par/iof5-multipatch-kerrschild.par
index 171ef919e..8ddd56f94 100644
--- a/CarpetDev/CarpetIOF5/par/iof5-multipatch-kerrschild.par
+++ b/CarpetDev/CarpetIOF5/par/iof5-multipatch-kerrschild.par
@@ -150,6 +150,8 @@ CarpetIOASCII::out1D_vars         = "
 
 
 
+CarpetIOF5::verbose = yes
+
 CarpetIOF5::out_every = 1024
 CarpetIOF5::out_vars  = "
         ADMBase::metric
diff --git a/CarpetDev/CarpetIOF5/param.ccl b/CarpetDev/CarpetIOF5/param.ccl
index b83b74675..bd8f14332 100644
--- a/CarpetDev/CarpetIOF5/param.ccl
+++ b/CarpetDev/CarpetIOF5/param.ccl
@@ -154,10 +154,19 @@ STRING checkpoint_dir "Checkpoint directory (overrides IO::checkpoint_dir)" STEE
 
 
 
+BOOLEAN use_chunks "Always use chunked layout" STEERABLE=always
+{
+} "yes"
+
+INT max_chunksize "Maximum chunk size (in bytes)" STEERABLE=always
+{
+  1:* :: ""
+} 131072                        # 128 kB
+
 INT compression_level "Compression level" STEERABLE=always
 {
   -1  :: "no compression"
-  0:9 :: "higher numbers compress better"
+  0:9 :: "smaller numbers compress faster, larger numbers compress better"
 } 1
 
 BOOLEAN use_checksums "Include a checksum for the data" STEERABLE=always
diff --git a/CarpetDev/CarpetIOF5/schedule.ccl b/CarpetDev/CarpetIOF5/schedule.ccl
index da8745b56..e633aef18 100644
--- a/CarpetDev/CarpetIOF5/schedule.ccl
+++ b/CarpetDev/CarpetIOF5/schedule.ccl
@@ -46,31 +46,3 @@ schedule CarpetIOF5_TerminationCheckpoint at TERMINATE
   LANG: C
   OPTIONS: global
 } "Termination checkpoint routine"
-
-
-
-#SCHEDULE F5_Output AT initial
-#{
-#  LANG: C
-#  OPTIONS: global-late
-#} "Create an output file"
-#
-#
-#
-#SCHEDULE F5_Poison AT initial AFTER F5_Output BEFORE F5_Input
-#{
-#  LANG: C
-#  OPTIONS: global-late
-#} "Poison all variables"
-#
-#SCHEDULE F5_Input AT initial AFTER F5_Output
-#{
-#  LANG: C
-#  OPTIONS: global-late
-#} "Read from file"
-#
-#SCHEDULE F5_Check AT initial AFTER F5_Input
-#{
-#  LANG: C
-#  OPTIONS: global-late
-#} "Check all variables for poison"
diff --git a/CarpetDev/CarpetIOF5/src/distribute.cc b/CarpetDev/CarpetIOF5/src/distribute.cc
index 5fd8e5a29..ab7183df3 100644
--- a/CarpetDev/CarpetIOF5/src/distribute.cc
+++ b/CarpetDev/CarpetIOF5/src/distribute.cc
@@ -2,6 +2,7 @@
 
 #include <cctk_Parameters.h>
 
+#include <cacheinfo.hh>
 #include <carpet.hh>
 
 #include <cassert>
@@ -52,10 +53,48 @@ namespace CarpetIOF5 {
   
   
   
+  /*** scatter_t::busy_tags_t *************************************************/
+  
+  scatter_t::busy_tags_t::busy_tags_t()
+  {
+    // do nothing
+  }
+  
+  scatter_t::busy_tags_t::~busy_tags_t()
+  {
+    // Ensure all tags are free
+    for (int tag=0; tag<(int)tags.size(); ++tag) {
+      assert(not tags.at(tag));
+    }
+  }
+  
+  int scatter_t::busy_tags_t::allocate_tag()
+  {
+    if (tags.empty()) {
+      tags.resize(max_concurrent_sends, false);
+    }
+    for (int tag=0; tag<(int)tags.size(); ++tag) {
+      if (not tags.at(tag)) {
+        tags.at(tag) = true;
+        return tag;
+      }
+    }
+    return -1;              // no free tag available
+  }
+  
+  void scatter_t::busy_tags_t::free_tag(int const tag)
+  {
+    assert(tags.at(tag));
+    tags.at(tag) = false;
+  }
+  
+  
+  
   /*** scatter_t **************************************************************/
   
   scatter_t::scatter_t(cGH const *const cctkGH_)
     : cctkGH(cctkGH_),
+      busy_tags(CCTK_nProcs(cctkGH)),
       num_received(0), num_sent(0), bytes_allocated(0),
       did_send_all(false),
       did_receive_sent_all(0), did_receive_all_sent_all(false)
@@ -136,7 +175,7 @@ namespace CarpetIOF5 {
     }
     
     if (verbose) {
-      CCTK_INFO("Destroying down global scatter object");
+      CCTK_INFO("Destroying global scatter object");
     }
     
     assert(public_recvs.empty());
@@ -152,7 +191,7 @@ namespace CarpetIOF5 {
   // consecutive)
   int child() { return dist::rank()*nsiblings() + 1; }
   // Get process id if my parent
-  int parent() { return (dist::rank()-1) / nsiblings(); }
+  int parent() { return dist::rank()==0 ? -1 : (dist::rank()-1) / nsiblings(); }
   
   // Check whether we should tell our parent that all our (ours and
   // our childrens') messages have been sent
@@ -164,7 +203,8 @@ namespace CarpetIOF5 {
       if (p < dist::size()) ++need_receive;
     }
     if (did_send_all and did_receive_sent_all == need_receive) {
-      if (dist::rank() == 0) {
+      int const p = parent();
+      if (p < 0) {
         // We are root; now broadcast this to all
         set_did_receive_all_sent_all();
       } else {
@@ -172,7 +212,6 @@ namespace CarpetIOF5 {
           CCTK_INFO("[Telling our parent that we and our children sent all our data]");
         }
         to_parent.state = fragdesc_t::state_sent_all;
-        int const p = parent();
         MPI_Request request;
         MPI_Isend(&to_parent, to_parent.num_ints(), MPI_INT, p, tag_desc,
                   dist::comm(), &request);
@@ -265,10 +304,25 @@ namespace CarpetIOF5 {
     while (not tosend.empty()) {
       list<transmission_t*>::iterator const tmi = tosend.begin();
       transmission_t *const tm = *tmi;
+      tosend.erase(tmi);
+      
+      // Choose tag
+      int const proc = tm->fragdesc.process;
+      int tag;
+      while (true) {
+        tag = busy_tags.at(proc).allocate_tag();
+        if (tag >= 0) break;
+        // No tag was available; wait for some progress before trying
+        // again (until the receiver has received some of our
+        // messages)
+        do_some_work(true);
+      }
+      tm->fragdesc.tag = tag_data_min + tag;
       
       if (verbose) {
         CCTK_VInfo(CCTK_THORNSTRING,
-                   "   Sending to process %d...", tm->fragdesc.process);
+                   "   Sending data to process %d with tag %d...",
+                   tm->fragdesc.process, tm->fragdesc.tag);
       }
       
       // Send descriptor and data
@@ -276,11 +330,10 @@ namespace CarpetIOF5 {
       MPI_Isend(&tm->fragdesc, tm->fragdesc.num_ints(), MPI_INT,
                 tm->fragdesc.process, tag_desc,
                 dist::comm(), &req);
+      assert(tm->request == MPI_REQUEST_NULL);
       MPI_Isend(&tm->data[0], tm->fragdesc.npoints(), tm->fragdesc.datatype(),
-                tm->fragdesc.process, tag_data,
+                tm->fragdesc.process, tm->fragdesc.tag,
                 dist::comm(), &tm->request);
-      
-      tosend.erase(tmi);
       sends.push_back(tm);
     }
     
@@ -352,8 +405,9 @@ namespace CarpetIOF5 {
       requests.push_back(tm->request);
       iterators.push_back(tmi);
     }
+    assert((int)requests.size() == nrequests);
     
-    // Wait for (or test for) some open requests
+    // Wait (or test) for some open requests
     int outcount;
     vector<int> indices(nrequests);
     vector<MPI_Status> statuses(nrequests);
@@ -377,32 +431,38 @@ namespace CarpetIOF5 {
     // Process all completed requests
     for (int n=0; n<outcount; ++n) {
       int const idx = indices.at(n);
+      MPI_Status const& status = statuses.at(n);
       
       if (idx < npublic_recvs) {
         
         // We received a new descriptor
-        int const source = statuses.at(idx).MPI_SOURCE;
-        
-        if (verbose) {
-          CCTK_VInfo(CCTK_THORNSTRING,
-                     "Receiving data from process %d", source);
-        }
+        int const source = status.MPI_SOURCE;
+        assert(source>=0 and source<dist::size());
         
         list<transmission_t*>::iterator const tmi = iterators.at(idx);
         transmission_t *const tm = *tmi;
+        assert(tm->request != MPI_REQUEST_NULL);
+        tm->request = MPI_REQUEST_NULL;
         public_recvs.erase(tmi);
         
         if (tm->fragdesc.state != fragdesc_t::state_normal) {
           handle_state_transition(tm->fragdesc);
         } else {
           
+          if (verbose) {
+            CCTK_VInfo(CCTK_THORNSTRING,
+                       "Receiving data from process %d with tag %d",
+                       source, tm->fragdesc.tag);
+          }
+          
           // Prepare receiving the data
           assert(tm->fragdesc.process == dist::rank());
           tm->data.resize(tm->fragdesc.npoints() * tm->fragdesc.vartypesize());
           bytes_allocated += tm->data.size();
+          assert(tm->request == MPI_REQUEST_NULL);
           MPI_Irecv(&tm->data[0],
                     tm->fragdesc.npoints(), tm->fragdesc.datatype(),
-                    source, tag_data,
+                    source, tm->fragdesc.tag,
                     dist::comm(), &tm->request);
           recvs.push_back(tm);
           if (verbose) {
@@ -418,6 +478,8 @@ namespace CarpetIOF5 {
         // We completed receiving a dataset; process it
         list<transmission_t*>::iterator const tmi = iterators.at(idx);
         transmission_t *const tm = *tmi;
+        assert(tm->request != MPI_REQUEST_NULL);
+        tm->request = MPI_REQUEST_NULL;
         
         if (verbose) {
           char *const fullname = CCTK_FullName(tm->fragdesc.varindex);
@@ -441,6 +503,8 @@ namespace CarpetIOF5 {
         // We completed sending a dataset; forget it
         list<transmission_t*>::iterator const tmi = iterators.at(idx);
         transmission_t *const tm = *tmi;
+        assert(tm->request != MPI_REQUEST_NULL);
+        tm->request = MPI_REQUEST_NULL;
         
         if (verbose) {
           char *const fullname = CCTK_FullName(tm->fragdesc.varindex);
@@ -449,6 +513,10 @@ namespace CarpetIOF5 {
           free(fullname);
         }
         
+        int const proc = tm->fragdesc.process;
+        int const tag = tm->fragdesc.tag - tag_data_min;
+        busy_tags.at(proc).free_tag(tag);
+        
         bytes_allocated -= tm->data.size();
         delete tm;
         if (verbose) {
@@ -484,11 +552,28 @@ namespace CarpetIOF5 {
                       baseext.lower() + fd.imax * baseext.stride(),
                       baseext.stride());
     
+    // Distribute load: don't start with sending to process 0, instead
+    // start with sending to (self+1)
+    int coffset;
+    {
+      int const nlc = hh.local_components(fd.reflevel);
+      if (nlc > 0) {
+        // find component of next processor
+        int c0 = hh.get_component(fd.reflevel, 0);
+        coffset = c0 + nlc;
+      } else {
+        // estimate component of next processor
+        int const nc = hh.components(fd.reflevel);
+        coffset = nc * dist::rank() / dist::size();
+      }
+    }
+    
     ibset done;
     list<transmission_t*> tosend;
     dh::light_cboxes const& light_cbox =
       dd.light_boxes.at(fd.mglevel).at(fd.reflevel);
-    for (int c=0; c<hh.components(fd.reflevel); ++c) {
+    for (int ci=0; ci<hh.components(fd.reflevel); ++ci) {
+      int const c = (ci + coffset) % hh.components(fd.reflevel);
       dh::light_dboxes const& light_box = light_cbox.at(c);
       ibbox const& intr = light_box.interior;
       ibbox const ovlp = mybox & intr;
@@ -567,11 +652,9 @@ namespace CarpetIOF5 {
     
     gh const& hh = *Carpet::arrdata.at(groupindex).at(fd.map).hh;
     dh const& dd = *Carpet::arrdata.at(groupindex).at(fd.map).dd;
-    ggf const& ff =
-      *Carpet::arrdata.at(groupindex).at(fd.map).data.at(varoffset);
+    ggf& ff = *Carpet::arrdata.at(groupindex).at(fd.map).data.at(varoffset);
     int const lc = hh.get_local_component(fd.reflevel, fd.component);
-    gdata const& data =
-      *ff.data_pointer(fd.timelevel, fd.reflevel, lc, fd.mglevel);
+    gdata& data = *ff.data_pointer(fd.timelevel, fd.reflevel, lc, fd.mglevel);
     
     ibbox const& baseext =
       hh.baseextents.AT(fd.mglevel).AT(fd.reflevel);
diff --git a/CarpetDev/CarpetIOF5/src/distribute.hh b/CarpetDev/CarpetIOF5/src/distribute.hh
index 420b584c4..96eaf6633 100644
--- a/CarpetDev/CarpetIOF5/src/distribute.hh
+++ b/CarpetDev/CarpetIOF5/src/distribute.hh
@@ -28,6 +28,7 @@ namespace CarpetIOF5 {
     
     // destination
     int component, process;
+    int tag;
     
     // meta-messages
     enum state_t { state_normal, state_sent_all, state_all_sent_all };
@@ -46,16 +47,33 @@ namespace CarpetIOF5 {
   // Scatter (distribute) Cactus variables
   class scatter_t {
     
-    enum tags { tag_desc, tag_data };
+    cGH const *const cctkGH;
+    
+    static int const max_concurrent_sends = 100;
+    enum tags_t {
+      tag_desc,
+      tag_data_min,
+      tag_data_max = tag_data_min + max_concurrent_sends
+    };
+    
+    class busy_tags_t {
+      vector<bool> tags;
+    public:
+      busy_tags_t();
+      ~busy_tags_t();
+      int allocate_tag();
+      void free_tag(int tag);
+    };
+    vector<busy_tags_t> busy_tags;
     
     struct transmission_t {
       fragdesc_t fragdesc;
       vector<char> data;
       MPI_Request request;
+      transmission_t(): request(MPI_REQUEST_NULL) {}
+      ~transmission_t() { assert(request==MPI_REQUEST_NULL); }
     };
     
-    cGH const *const cctkGH;
-    
     // Desired number of public receives to keep posted at all times
     static int const num_public_recvs = 10;
     list<transmission_t*> public_recvs, recvs, sends;
@@ -68,7 +86,7 @@ namespace CarpetIOF5 {
     fragdesc_t to_parent, to_children;
     
   public:
-    scatter_t(cGH const *const cctkGH_);
+    scatter_t(cGH const *cctkGH_);
     ~scatter_t();
     
     void send(fragdesc_t const& fd, void const *data);
diff --git a/CarpetDev/CarpetIOF5/src/input.cc b/CarpetDev/CarpetIOF5/src/input.cc
index b7f25b59f..eb2db82e5 100644
--- a/CarpetDev/CarpetIOF5/src/input.cc
+++ b/CarpetDev/CarpetIOF5/src/input.cc
@@ -119,15 +119,41 @@ namespace CarpetIOF5 {
       
       // TODO: synchronise all read grid functions
     }
-  
+    
+    
+    
     void read_grid(F5Path *const path)
     {
       indent_t indent;
       cout << indent << "grid=" << gridname << "\n";
+      
+      
+      
+      if (input_metadata) {
+        if (reflevel == 0) {
+          indent_t indent2;
+          cout << indent2 << "reading metadata\n";
+          // hid_t const metadata_group = path->Grid_hid;
+          ostringstream pathname;
+          pathname << FIBER_CONTENT_GRIDS << "/" << gridname;
+          hid_t group;
+          group = H5Gopen(path->ContentsGroup_hid, pathname.str().c_str(),
+                          H5P_DEFAULT);
+          assert(group >= 0);
+          read_metadata(cctkGH, group);
+          herr_t const herr = H5Gclose(group);
+          assert(not herr);
+        }
+      }
+      
+      
+      
       // F5iterate_vertex_fields(path, NULL, field_iterator, this, NULL, NULL);
       F5iterate_topologies(path, NULL, topology_iterator, this);
     }
     
+    
+    
     void read_topology(F5Path *const path)
     {
       indent_t indent;
@@ -151,9 +177,9 @@ namespace CarpetIOF5 {
         herr = H5Gget_linkval
           (path->Grid_hid, topologyname, sizeof linkval, linkval);
         assert(not herr);
-        indent_t indent;
-        cout << indent << "alias for topology \"" << linkval << "\"\n"
-             << indent << "ignoring this topology\n";
+        indent_t indent2;
+        cout << indent2 << "alias for topology \"" << linkval << "\"\n"
+             << indent2 << "ignoring this topology\n";
         return;
       }
       
@@ -185,6 +211,8 @@ namespace CarpetIOF5 {
       F5iterate_topology_fields(path, NULL, field_iterator, this, NULL, NULL); 
     }
     
+    
+    
     void read_field(F5Path *const path)
     {
       indent_t indent;
@@ -220,8 +248,8 @@ namespace CarpetIOF5 {
         F5Fclose(path);
         
       } else {
-        indent_t indent;
-        cout << indent << "ignoring this field\n";
+        indent_t indent2;
+        cout << indent2 << "ignoring this field\n";
       }
       // TODO: keep track of which fields have been read, and complain
       // about unread ones
@@ -229,6 +257,8 @@ namespace CarpetIOF5 {
       // complain about unread parts
     }
     
+    
+    
     void read_fragment(F5Path *const path)
     {
       indent_t indent;
@@ -293,6 +323,8 @@ namespace CarpetIOF5 {
       }
     }
     
+    
+    
     void read_variable(F5Path *const path, char const *const name,
                        int const var)
     {
@@ -340,6 +372,10 @@ namespace CarpetIOF5 {
       // fragment, and (b) there exists a group with this name. (We
       // probably could examine F5's internal state as well, but
       // looking for an HDF5 group is the most robust approach.)
+      
+      // Note: F5 has a function F5Fis_group(path) that returns
+      // whether the field is a group or a dataset.
+      
       bool fragment_is_group = false;
       if (fragmentname) {
         H5O_info_t info;
@@ -432,8 +468,8 @@ namespace CarpetIOF5 {
       
       ibbox const fbox(foff, foff+flen-1, 1);
       {
-        indent_t indent;
-        cout << indent
+        indent_t indent2;
+        cout << indent2
              << "dataset bbox is " << foff << ":" << foff+flen-1 << "\n";
       }
       
@@ -603,7 +639,7 @@ namespace CarpetIOF5 {
     // Grid structure
     string gs;
     ReadLargeAttribute(group, grid_structure, gs);
-    // TODO: set grid structure
+    deserialise_grid_structure(cctkGH, gs);
     
     herr = H5Gclose(group);
     assert(not herr);
diff --git a/CarpetDev/CarpetIOF5/src/iof5.cc b/CarpetDev/CarpetIOF5/src/iof5.cc
index 1141a7471..1642bfe0e 100644
--- a/CarpetDev/CarpetIOF5/src/iof5.cc
+++ b/CarpetDev/CarpetIOF5/src/iof5.cc
@@ -86,7 +86,7 @@ namespace CarpetIOF5 {
   
   // A mechanism to keep the HDF5 output file open across multiple
   // write operations:
-  hid_t file = H5I_INVALID_HID;
+  hid_t open_file = H5I_INVALID_HID;
   int keep_file_open = 0;
   void enter_keep_file_open()
   {
@@ -96,10 +96,10 @@ namespace CarpetIOF5 {
   {
     assert(keep_file_open > 0);
     --keep_file_open;
-    if (keep_file_open==0 and file>=0) {
-      herr_t const herr = H5Fclose(file);
+    if (keep_file_open==0 and open_file>=0) {
+      herr_t const herr = H5Fclose(open_file);
       assert(not herr);
-      file = H5I_INVALID_HID;
+      open_file = H5I_INVALID_HID;
       H5garbage_collect();
     }
   }
@@ -294,22 +294,22 @@ namespace CarpetIOF5 {
       enter_keep_file_open();
       bool const truncate_file =
         first_time and IO_TruncateOutputFiles(cctkGH) and myproc == myioproc;
-      if (file < 0) {
+      if (open_file < 0) {
         // Reuse file hid if file is already open
         hid_t fapl = H5Pcreate(H5P_FILE_ACCESS);
         H5Pset_fclose_degree(fapl, H5F_CLOSE_STRONG);
-        file =
+        open_file =
           truncate_file ?
           H5Fcreate(name.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, fapl) :
           H5Fopen  (name.c_str(), H5F_ACC_RDWR , fapl);
-        assert(file >= 0);
+        assert(open_file >= 0);
         H5Pclose(fapl);
       }
       first_time = false;
       
       vector<bool> output_var(CCTK_NumVars());
       output_var.at(vindex) = true;
-      output(cctkGH, file, output_var, false, true);
+      output(cctkGH, open_file, output_var, false, true);
       
       // Close file
       leave_keep_file_open();
@@ -471,16 +471,22 @@ namespace CarpetIOF5 {
     
     
     
+    assert(called_from == CP_RECOVER_PARAMETERS or
+           called_from == CP_RECOVER_DATA or
+           called_from == FILEREADER_DATA);
+    bool const in_recovery =
+      called_from == CP_RECOVER_PARAMETERS or
+      called_from == CP_RECOVER_DATA;
+    if (not in_recovery) {
+      assert(is_level_mode());
+      if (reflevel > 0) return 0; // no error
+    }
+    
+    
+    
     BEGIN_GLOBAL_MODE(cctkGH) {
       DECLARE_CCTK_ARGUMENTS;
       
-      assert(called_from == CP_RECOVER_PARAMETERS or
-             called_from == CP_RECOVER_DATA or
-             called_from == FILEREADER_DATA);
-      bool const in_recovery =
-        called_from == CP_RECOVER_PARAMETERS or
-        called_from == CP_RECOVER_DATA;
-      
       if (in_recovery) {
         CCTK_VInfo(CCTK_THORNSTRING, "F5::Input: recovering iteration %d",
                    recovery_iteration);
@@ -505,35 +511,28 @@ namespace CarpetIOF5 {
       
       
       
-      scatter_t scatter(cctkGH);
+      // Keep track of which files could be read, and which could not
+      int foundproc = -1, notfoundproc = -1;
       
-      // Open file
       // string const basename =
       //   in_recovery
       //   ? "checkpoint"
       //   : generate_basename(cctkGH, CCTK_VarIndex("grid::r"));
       string const basename = basefilename;
       
-      // Keep track of which files could be read, and which could not
-      int foundproc = -1, notfoundproc = -1;
-      
-      // TODO: Store how many processes contributed to the output, and
-      // expect exactly that many files
-      int const myproc = CCTK_MyProc(cctkGH);
-      int const nprocs = CCTK_nProcs(cctkGH);
-      
-      int const ioproc_every =
-        max_nioprocs == 0 ? 1 : (nprocs + max_nioprocs - 1) / max_nioprocs;
-      assert(ioproc_every > 0);
-      int const nioprocs = nprocs / ioproc_every;
-      assert(nioprocs > 0 and nioprocs <= max_nioprocs);
-      int const myioproc = myproc / ioproc_every * ioproc_every;
-      
-      if (myproc == myioproc) {
+      {
+        scatter_t scatter(cctkGH);
+        
+        // Iterate over files
+        
+        // TODO: Store how many processes contributed to the output,
+        // and expect exactly that many files
+        int const myproc = CCTK_MyProc(cctkGH);
+        int const nprocs = CCTK_nProcs(cctkGH);
         // Loop over all (possible) files
-        for (int ioproc = myioproc / ioproc_every; ; ioproc += nioprocs) {
+        for (int proc=myproc; ; proc+=nprocs) {
           string const name =
-            create_filename(cctkGH, basename, cctkGH->cctk_iteration, ioproc,
+            create_filename(cctkGH, basename, cctkGH->cctk_iteration, proc,
                             in_recovery ? io_dir_recover : io_dir_input, false);
           
           bool file_exists;
@@ -541,24 +540,21 @@ namespace CarpetIOF5 {
             file_exists = H5Fis_hdf5(name.c_str()) > 0;
           } H5E_END_TRY;
           if (not file_exists) {
-            notfoundproc = ioproc;
+            notfoundproc = proc;
             break;
           }
-          foundproc = ioproc;
+          foundproc = proc;
           
           indent_t indent;
-          cout << indent << "I/O process=" << ioproc << "\n";
+          cout << indent << "process=" << proc << "\n";
           
-          hid_t fapl = H5Pcreate(H5P_FILE_ACCESS);
-          H5Pset_fclose_degree(fapl, H5F_CLOSE_STRONG);
-          hid_t const file = H5Fopen(name.c_str(), H5F_ACC_RDONLY, fapl);
+          hid_t const file = H5Fopen(name.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
           assert(file >= 0);
-          H5Pclose(fapl);
           
           // Iterate over all time slices
           bool const input_past_timelevels = in_recovery;
-          // TODO: read metadata when recoverying parameters
-          bool const input_metadata = false;
+          // Read metadata when recoverying parameters
+          bool const input_metadata = called_from == CP_RECOVER_PARAMETERS;
           input(cctkGH, file, input_var, input_past_timelevels, input_metadata,
                 scatter);
           
@@ -567,6 +563,8 @@ namespace CarpetIOF5 {
           assert(not herr);
           H5garbage_collect();
         }
+        
+        // Destroy scatter object
       }
       
       {
@@ -644,10 +642,10 @@ namespace CarpetIOF5 {
       int const iter = strtol(p, &p, 10);
       if (!*p) continue;
       
-      // Read the process number
+      // Read (and ignore) the process number
       if (infix.compare(0, infix.length(), p, infix.length()) != 0) continue;
       p += infix.length();
-      int const proc = strtol(p, &p, 10);
+      /* int const proc = */ strtol(p, &p, 10);
       if (!*p) continue;
       
       // Finally check the file extension
diff --git a/CarpetDev/CarpetIOF5/src/iof5.hh b/CarpetDev/CarpetIOF5/src/iof5.hh
index c082206c3..a450ae812 100644
--- a/CarpetDev/CarpetIOF5/src/iof5.hh
+++ b/CarpetDev/CarpetIOF5/src/iof5.hh
@@ -264,13 +264,19 @@ namespace CarpetIOF5 {
   
   
   
+  // Write/read Cactus metadata to a particular location in an HDF5
+  // file
   void
-  write_metadata(cGH const *const cctkGH, hid_t const file);
+  write_metadata(cGH const *const cctkGH, hid_t const group);
+  void
+  read_metadata(cGH const *const cctkGH, hid_t const group);
   
   // Handle Carpet's grid structure (this should move to Carpet and/or
   // CarpetLib)
   string
   serialise_grid_structure(cGH const *const cctkGH);
+  void
+  deserialise_grid_structure(cGH const *const cctkGH, string const buf);
   
   
   
diff --git a/CarpetDev/CarpetIOF5/src/output.cc b/CarpetDev/CarpetIOF5/src/output.cc
index aad6174f2..d7c2f4649 100644
--- a/CarpetDev/CarpetIOF5/src/output.cc
+++ b/CarpetDev/CarpetIOF5/src/output.cc
@@ -104,8 +104,7 @@ namespace CarpetIOF5 {
     
     struct component_indices_t: map_indices_t {
       // elements >=dim remain undefined
-      ivect ash;
-      ivect lbnd, lsh, lghosts, ughosts;
+      ivect lbnd, lsh, ash, lghosts, ughosts;
       ivect imin, imax, ioff, ilen;
       rvect clower, cupper;
       
@@ -758,14 +757,37 @@ namespace CarpetIOF5 {
       // Dataset properties
       hid_t const prop = H5Pcreate(H5P_DATASET_CREATE);
       assert(prop >= 0);
+      // Use chunked I/O if requested or needed
+      bool const set_chunks =
+        use_chunks or compression_level >= 0 or use_checksums;
+      if (set_chunks) {
+        assert(ci.dim > 0);
+        int const chunksize = max_chunksize / H5Tget_size(type);
+        assert(chunksize > 0);
+        int cs = 1, csd = 1;
+        while (csd < chunksize) {
+          cs <<= 1; csd <<= ci.dim;
+          assert(ci.dim > 0);
+        }
+        ivect chunkshape = cs;
+        chunkshape[ci.dim-1] = cs / (csd / chunksize);
+        chunkshape = max(chunkshape, 1);
+        chunkshape = min(chunkshape, ci.ilen);
+        FAILWARN(H5Pset_chunk(prop, ci.dim,
+                              &v2h(chunkshape).reverse()[dim-ci.dim]));
+      }
       // Enable compression if requested
       if (compression_level >= 0) {
-        FAILWARN(H5Pset_chunk(prop, ci.dim, &v2h(ci.ilen).reverse()[0]));
+        assert(set_chunks);
+        indent_t indent2;
+        cout << indent2
+             << "using chunk size " << ci.ilen << " "
+             << "(" << prod(ci.ilen) << ")\n";
         FAILWARN(H5Pset_deflate(prop, compression_level));
       }
       // Enable checksums if requested
       if (use_checksums) {
-        FAILWARN(H5Pset_chunk(prop, ci.dim, &v2h(ci.ilen).reverse()[0]));
+        assert(set_chunks);
         FAILWARN(H5Pset_fletcher32(prop));
       }
       
@@ -773,7 +795,8 @@ namespace CarpetIOF5 {
         // Write single-component tensors into non-separated fractions
         // for convenience (could also use a separated compound
         // instead)
-        // TODO: Extent F5 API to allow writing non-fragmented datasets
+        // TODO: Extend F5 API to allow writing non-fragmented datasets
+        // TODO: Extend F5 API to accept hyperslab descriptors
         FAILWARN
           (F5Fwrite_fraction(path, name.c_str(),
                              ci.dim,
@@ -798,6 +821,7 @@ namespace CarpetIOF5 {
                     "This does not seem to be supported by F5ls "
                     "(or is implemented wrong in the F5 library).");
         }
+        // TODO: Extend F5 API to accept hyperslab descriptors
         FAILWARN
           (F5FSwrite_fraction(path, name.c_str(),
                               ci.dim,
diff --git a/CarpetDev/CarpetIOF5/src/util.cc b/CarpetDev/CarpetIOF5/src/util.cc
index 3d924e775..b2ef6fae2 100644
--- a/CarpetDev/CarpetIOF5/src/util.cc
+++ b/CarpetDev/CarpetIOF5/src/util.cc
@@ -381,7 +381,7 @@ namespace CarpetIOF5 {
   
   
   
-  char const *const grid_structure = "Grid Structure v5";
+  char const *const grid_structure = "Grid Structure v6";
   
   string
   serialise_grid_structure(cGH const *const cctkGH)
@@ -398,14 +398,233 @@ namespace CarpetIOF5 {
       buf << "regions:" << vhh.AT(m)->regions.AT(0) << ",";
       buf << "ghost_widths:" << vdd.AT(m)->ghost_widths << ",";
       buf << "buffer_widths:" << vdd.AT(m)->buffer_widths << ",";
+      buf << "overlap_widths:" << vdd.AT(m)->overlap_widths << ",";
       buf << "prolongation_orders_space:"
           << vdd.AT(m)->prolongation_orders_space << ",";
       buf << "},";
     }
     buf << "],";
+    buf << "global_time:" << global_time << ",";
+    buf << "delta_time:" << delta_time << ",";
     buf << "times:" << *tt << ",";
     buf << "}.";
     return buf.str();
   }
   
+  void
+  deserialise_grid_structure(cGH const *const cctkGH, string const str)
+  {
+    // Define variables for the metadata that will be read in
+    int const my_maps = maps;
+    vector<gh::rregs> my_superregions(my_maps);
+    vector<gh::rregs> my_regions(my_maps);
+    vector<vector<i2vect> > my_ghost_widths(my_maps);
+    vector<vector<i2vect> > my_buffer_widths(my_maps);
+    vector<vector<i2vect> > my_overlap_widths(my_maps);
+    vector<int> my_prolongation_orders_space(my_maps);
+    CCTK_REAL my_global_time, my_delta_time;
+    th my_tt(tt->h, tt->time_interpolation_during_regridding);
+    
+    // Read in the metadata
+    istringstream is(str);
+    skipws(is);
+    consume(is, grid_structure);
+    skipws(is);
+    consume(is, "{");
+    skipws(is);
+    consume(is, "maps");
+    skipws(is);
+    consume(is, '=');
+    skipws(is);
+    consume(is, '[');
+    for (int m=0; m<my_maps; ++m) {
+      skipws(is);
+      consume(is, "[");
+      int my_m;
+      is >> my_m;
+      assert(my_m == m);
+      skipws(is);
+      consume(is, ']');
+      skipws(is);
+      consume(is, '=');
+      skipws(is);
+      consume(is, '{');
+      skipws(is);
+      consume(is, "superregions:");
+      skipws(is);
+      is >> my_superregions.AT(m);
+      skipws(is);
+      consume(is, ',');
+      skipws(is);
+      consume(is, "regions:");
+      skipws(is);
+      // ignore region specification
+      // is >> my_regions.AT(m);
+      gh::rregs tmp_regions;
+      is >> tmp_regions;
+      skipws(is);
+      consume(is, ',');
+      skipws(is);
+      consume(is, "ghost_widths:");
+      skipws(is);
+      is >> my_ghost_widths.AT(m);
+      skipws(is);
+      consume(is, ',');
+      skipws(is);
+      consume(is, "buffer_widths:");
+      skipws(is);
+      is >> my_buffer_widths.AT(m);
+      skipws(is);
+      consume(is, ',');
+      skipws(is);
+      consume(is, "overlap_widths:");
+      skipws(is);
+      is >> my_overlap_widths.AT(m);
+      skipws(is);
+      consume(is, ',');
+      skipws(is);
+      consume(is, "prolongation_orders_space:");
+      skipws(is);
+      is >> my_prolongation_orders_space.AT(m);
+      skipws(is);
+      consume(is, ',');
+      skipws(is);
+      consume(is, '}');
+      skipws(is);
+      consume(is, ',');
+    }
+    skipws(is);
+    consume(is, ']');
+    skipws(is);
+    consume(is, ',');
+    skipws(is);
+    consume(is, "global_time:");
+    is >> my_global_time;
+    skipws(is);
+    consume(is, ',');
+    skipws(is);
+    consume(is, "delta_time:");
+    is >> my_delta_time;
+    skipws(is);
+    consume(is, ',');
+    skipws(is);
+    consume(is, "times:");
+    is >> my_tt;
+    skipws(is);
+    consume(is, ',');
+    skipws(is);
+    consume(is, '}');
+    skipws(is);
+    consume(is, '.');
+    
+    // Change Carpet's grid structure according to what has just been
+    // read
+    {
+      int type;
+      void const *const ptr =
+        CCTK_ParameterGet("regrid_in_level_mode", "Carpet", & type);
+      assert(ptr);
+      assert(type == PARAMETER_BOOLEAN);
+      CCTK_INT const regrid_in_level_mode = *static_cast<CCTK_INT const*>(ptr);
+      assert(regrid_in_level_mode);
+    }
+    
+    // Count refinement levels
+    vector<int> rls(maps);
+    for (int m=0; m<maps; ++m) {
+      rls.AT(m) = my_superregions.AT(m).size();
+    }
+    int maxrl = 0;
+    for (int m=0; m<maps; ++m) {
+      maxrl = max(maxrl, rls.AT(m));
+    }
+    // All maps must have the same number of refinement levels
+    for (int m=0; m<maps; ++m) {
+      my_superregions.AT(m).resize(maxrl);
+      my_regions.AT(m).resize(maxrl);
+    }
+    
+    // Make multiprocessor aware
+    for (int rl=0; rl<maxrl; ++rl) {
+      vector<gh::cregs> superregss(maps);
+      for (int m=0; m<maps; ++m) {
+        superregss.AT(m) = my_superregions.AT(m).AT(rl);
+      }
+      vector<gh::cregs> regss(maps);
+      Carpet::SplitRegionsMaps(cctkGH, superregss, regss);
+      for (int m=0; m<maps; ++m) {
+        my_superregions.AT(m).AT(rl) = superregss.AT(m);
+        my_regions.AT(m).AT(rl) = regss.AT(m);
+      }
+    } // for rl
+    
+    // Make multigrid aware
+    vector<gh::mregs> my_multigrid_regions(maps);
+    Carpet::MakeMultigridBoxesMaps(cctkGH, my_regions, my_multigrid_regions);
+    
+    // Regrid (ignoring the previous content of the grid hierarchy)
+    for (int m=0; m<maps; ++m) {
+      RegridMap
+        (cctkGH, m, my_superregions.AT(m), my_multigrid_regions.AT(m), false);
+    } // for m
+    
+    // Set up time hierarchy after RegridMap created it
+    for (int ml=0; ml<vhh.AT(0)->mglevels(); ++ ml) {
+      for (int rl=0; rl<vhh.AT(0)->reflevels(); ++ rl) {
+        for (int tl=0; tl<tt->timelevels; ++ tl) {
+          tt->set_time(ml, rl, tl, my_tt.get_time(ml, rl, tl));
+        }
+        tt->set_delta(ml, rl, my_tt.get_delta(ml, rl));
+      }
+    }
+    
+    PostRegrid(cctkGH);
+    
+    // Ensure we are in global mode
+    assert(is_global_mode());
+    global_time = my_global_time;
+    (const_cast<cGH*>(cctkGH))->cctk_time = global_time;
+    delta_time = my_delta_time;
+    timereflevelfact = timereffacts.AT(reflevels - 1);
+    spacereflevelfact = ivect(-get_deadbeef());
+    ivect::ref(cctkGH->cctk_levfac) = spacereflevelfact;
+    (const_cast<cGH*>(cctkGH))->cctk_timefac = timereflevelfact;
+    
+    // Recompose the grid structure (ignoring the previous content of
+    // the grid hierarchy)
+    for (int rl=0; rl<reflevels; ++rl) {
+      Recompose(cctkGH, rl, false);
+    }
+    
+    // Free old grid structure
+    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(vdd.AT(m)->ghost_widths.size() == my_ghost_widths.AT(m).size());
+      for (int rl=0; rl<int(my_ghost_widths.AT(m).size()); ++rl) {
+        assert(all(all(vdd.AT(m)->ghost_widths.AT(rl) ==
+                       my_ghost_widths.AT(m).AT(rl))));
+      }
+      assert(vdd.AT(m)->buffer_widths.size() == my_buffer_widths.AT(m).size());
+      for (int rl=0; rl<int(my_buffer_widths.AT(m).size()); ++rl) {
+        assert(all(all(vdd.AT(m)->buffer_widths.AT(rl) ==
+                       my_buffer_widths.AT(m).AT(rl))));
+      }
+      assert(vdd.AT(m)->overlap_widths.size() ==
+             my_overlap_widths.AT(m).size());
+      for (int rl=0; rl<int(my_overlap_widths.AT(m).size()); ++rl) {
+        assert(all(all(vdd.AT(m)->overlap_widths.AT(rl) ==
+                       my_overlap_widths.AT(m).AT(rl))));
+      }
+      for (int rl=0; rl<int(my_prolongation_orders_space.size()); ++rl) {
+        assert(vdd.AT(m)->prolongation_orders_space.AT(rl) ==
+               my_prolongation_orders_space.AT(rl));
+      }
+    }
+  }
+  
 } // end namespace CarpetIOF5