CarpetIOF5: Limit number of I/O processes

3 files changed, 114 insertions, 61 deletions

diff --git a/CarpetDev/CarpetIOF5/param.ccl b/CarpetDev/CarpetIOF5/param.ccl
index a002cd016..b83b74675 100644
--- a/CarpetDev/CarpetIOF5/param.ccl
+++ b/CarpetDev/CarpetIOF5/param.ccl
@@ -91,6 +91,14 @@ BOOLEAN one_dir_per_file "Create one subdirectory per output file to reduce lock
 
 
 
+INT max_nioprocs "Maximum number of I/O processes" STEERABLE=always
+{
+  0   :: "unlimited"
+  1:* :: "that many processes"
+} 0
+
+
+
 BOOLEAN separate_single_component_tensors "Use a separated representation even for single-component tensors" STEERABLE=always
 {
 } "no"
diff --git a/CarpetDev/CarpetIOF5/src/distribute.cc b/CarpetDev/CarpetIOF5/src/distribute.cc
index 05658ee42..ae3139280 100644
--- a/CarpetDev/CarpetIOF5/src/distribute.cc
+++ b/CarpetDev/CarpetIOF5/src/distribute.cc
@@ -21,10 +21,7 @@ namespace CarpetIOF5 {
   
   int fragdesc_t::npoints() const
   {
-    int np = 1;
-    for (int d=0; d<dim; ++d) {
-      np *= imax[d] - imin[d] + 1;
-    }
+    int const np = prod(imax-imin+1);
     assert(np>0);
     return np;
   }
@@ -479,7 +476,6 @@ namespace CarpetIOF5 {
     
     gh const& hh = *Carpet::arrdata.at(groupindex).at(fd.map).hh;
     dh const& dd = *Carpet::arrdata.at(groupindex).at(fd.map).dd;
-    th const& tt = *Carpet::arrdata.at(groupindex).at(fd.map).tt;
     
     ibbox const& baseext =
       hh.baseextents.AT(fd.mglevel).AT(fd.reflevel);
@@ -512,11 +508,11 @@ namespace CarpetIOF5 {
         ptrdiff_t const ni = tm->fragdesc.imax[0] - tm->fragdesc.imin[0] + 1;
         ptrdiff_t const nj = tm->fragdesc.imax[1] - tm->fragdesc.imin[1] + 1;
         ptrdiff_t const nk = tm->fragdesc.imax[2] - tm->fragdesc.imin[2] + 1;
-        ptrdiff_t const di = 1;
-        ptrdiff_t const dj = ni;
-        ptrdiff_t const dk = ni * nj;
-        ptrdiff_t const np = ni * nj * nk;
-        tm->data.resize(np * vartypesize);
+        ptrdiff_t const di = vartypesize;
+        ptrdiff_t const dj = di * ni;
+        ptrdiff_t const dk = dj * nj;
+        ptrdiff_t const np = dk * nk;
+        tm->data.resize(np);
         bytes_allocated += tm->data.size();
         if (verbose) {
           CCTK_VInfo(CCTK_THORNSTRING,
@@ -527,21 +523,24 @@ namespace CarpetIOF5 {
         ptrdiff_t const nis = fd.imax[0] - fd.imin[0] + 1;
         ptrdiff_t const njs = fd.imax[1] - fd.imin[1] + 1;
         ptrdiff_t const nks = fd.imax[2] - fd.imin[2] + 1;
-        ptrdiff_t const dis = 1;
-        ptrdiff_t const djs = nis;
-        ptrdiff_t const dks = nis * njs;
-        ptrdiff_t const nps = nis * njs * nks;
+        ptrdiff_t const dis = vartypesize;
+        ptrdiff_t const djs = vartypesize * nis;
+        ptrdiff_t const dks = vartypesize * nis * njs;
+        ptrdiff_t const nps = vartypesize * nis * njs * nks;
         ptrdiff_t const i0s = tm->fragdesc.imin[0] - fd.imin[0];
         ptrdiff_t const j0s = tm->fragdesc.imin[1] - fd.imin[1];
         ptrdiff_t const k0s = tm->fragdesc.imin[2] - fd.imin[2];
         ptrdiff_t const ind0s = i0s * dis + j0s * djs + k0s * dks;
-        char const *const src = &((char const*)data)[vartypesize * ind0s];
+        assert(ind0s < nps);
+        char const *const src = &((char const*)data)[ind0s];
         
 #pragma omp parallel for //collapse(2)
         for (ptrdiff_t k=0; k<nk; ++k) {
           for (ptrdiff_t j=0; j<nj; ++j) {
             ptrdiff_t const ind = j*dj + k*dk;
             ptrdiff_t const inds = j*djs + k*dks;
+            assert(ind+ni*vartypesize <= np);
+            assert(inds+ni*vartypesize <= nps);
             memcpy(&dst[ind], &src[inds], ni*vartypesize);
           }
         }
@@ -568,7 +567,6 @@ namespace CarpetIOF5 {
     
     gh const& hh = *Carpet::arrdata.at(groupindex).at(fd.map).hh;
     dh const& dd = *Carpet::arrdata.at(groupindex).at(fd.map).dd;
-    th const& tt = *Carpet::arrdata.at(groupindex).at(fd.map).tt;
     ggf const& ff =
       *Carpet::arrdata.at(groupindex).at(fd.map).data.at(varoffset);
     int const lc = hh.get_local_component(fd.reflevel, fd.component);
@@ -592,10 +590,10 @@ namespace CarpetIOF5 {
     ptrdiff_t const ni = fd.imax[0] - fd.imin[0] + 1;
     ptrdiff_t const nj = fd.imax[1] - fd.imin[1] + 1;
     ptrdiff_t const nk = fd.imax[2] - fd.imin[2] + 1;
-    ptrdiff_t const di = 1;
-    ptrdiff_t const dj = ni;
-    ptrdiff_t const dk = ni * nj;
-    ptrdiff_t const np = ni * nj * nk;
+    ptrdiff_t const di = vartypesize;
+    ptrdiff_t const dj = di * ni;
+    ptrdiff_t const dk = dj * nj;
+    ptrdiff_t const np = dk * nk;
     char const *const src = (char const*)&tm->data[0];
     
     ivect const lbnd = (extr.lower() - baseext.lower()) / baseext.stride();
@@ -603,10 +601,10 @@ namespace CarpetIOF5 {
     ptrdiff_t const nid = lsh[0];
     ptrdiff_t const njd = lsh[1];
     ptrdiff_t const nkd = lsh[2];
-    ptrdiff_t const did = 1;
-    ptrdiff_t const djd = nid;
-    ptrdiff_t const dkd = nid * njd;
-    ptrdiff_t const npd = nid * njd * nkd;
+    ptrdiff_t const did = vartypesize;
+    ptrdiff_t const djd = vartypesize * nid;
+    ptrdiff_t const dkd = vartypesize * nid * njd;
+    ptrdiff_t const npd = vartypesize * nid * njd * nkd;
     ptrdiff_t const i0d = fd.imin[0] - lbnd[0];
     ptrdiff_t const j0d = fd.imin[1] - lbnd[1];
     ptrdiff_t const k0d = fd.imin[2] - lbnd[2];
@@ -619,6 +617,8 @@ namespace CarpetIOF5 {
       for (ptrdiff_t j=0; j<nj; ++j) {
         ptrdiff_t const indd = j*djd + k*dkd;
         ptrdiff_t const ind = j*dj + k*dk;
+        assert(indd+ni*vartypesize <= npd);
+        assert(ind+ni*vartypesize <= np);
         memcpy(&dst[indd], &src[ind], ni*vartypesize);
       }
     }
diff --git a/CarpetDev/CarpetIOF5/src/iof5.cc b/CarpetDev/CarpetIOF5/src/iof5.cc
index c9ebb028f..469a7c105 100644
--- a/CarpetDev/CarpetIOF5/src/iof5.cc
+++ b/CarpetDev/CarpetIOF5/src/iof5.cc
@@ -249,6 +249,29 @@ namespace CarpetIOF5 {
       // We don't know how to open multiple files yet
       assert(CCTK_EQUALS(file_content, "everything"));
       
+      // Determine number of I/O processes
+      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;
+      // We split processes into "I/O groups" which range from
+      // myioproc to myioproc+ioproc_every-1 (inclusive). Within each
+      // group, at most one process can perform I/O.
+      
+      // If I am not the first process in my I/O group, wait for a
+      // token from my predecessor
+      if (myproc > myioproc) {
+        MPI_Recv(NULL, 0, MPI_INT, myproc-1, 0, dist::comm(),
+                 MPI_STATUS_IGNORE);
+      }
+      
+      
+      
       // Open file
       static bool first_time = true;
       
@@ -257,17 +280,17 @@ namespace CarpetIOF5 {
       int const vindex = CCTK_VarIndex(varname);
       assert(vindex >= 0);
       string const basename = generate_basename(cctkGH, vindex);
-      int const myproc = CCTK_MyProc(cctkGH);
-      int const proc = myproc;
+      int const ioproc = myioproc / ioproc_every;
       string const name =
-        create_filename(cctkGH, basename, cctkGH->cctk_iteration, proc,
+        create_filename(cctkGH, basename, cctkGH->cctk_iteration, ioproc,
                         io_dir_output, first_time);
       
       indent_t indent;
-      cout << indent << "process=" << proc << "\n";
+      cout << indent << "I/O process=" << ioproc << "\n";
       
       enter_keep_file_open();
-      bool const truncate_file = first_time and IO_TruncateOutputFiles(cctkGH);
+      bool const truncate_file =
+        first_time and IO_TruncateOutputFiles(cctkGH) and myproc == myioproc;
       if (file < 0) {
         // Reuse file hid if file is already open
         file =
@@ -285,6 +308,16 @@ namespace CarpetIOF5 {
       // Close file
       leave_keep_file_open();
       
+      
+      
+      // If I am not the last process in my I/O group, send a token to
+      // my successor
+      if (myproc < max(myioproc + ioproc_every, nprocs) - 1) {
+        MPI_Send(NULL, 0, MPI_INT, myproc+1-1, 0, dist::comm());
+      }
+      
+      
+      
     } END_GLOBAL_MODE;
     
     return 0;                   // no error
@@ -478,38 +511,48 @@ namespace CarpetIOF5 {
       // expect exactly that many files
       int const myproc = CCTK_MyProc(cctkGH);
       int const nprocs = CCTK_nProcs(cctkGH);
-      // Loop over all (possible) files
-      for (int proc=myproc; ; proc+=nprocs) {
-        string const name =
-          create_filename(cctkGH, basename, cctkGH->cctk_iteration, proc,
-                          in_recovery ? io_dir_recover : io_dir_input, false);
-        
-        bool file_exists;
-        H5E_BEGIN_TRY {
-          file_exists = H5Fis_hdf5(name.c_str()) > 0;
-        } H5E_END_TRY;
-        if (not file_exists) {
-          notfoundproc = proc;
-          break;
+      
+      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) {
+        // Loop over all (possible) files
+        for (int ioproc = myioproc / ioproc_every; ; ioproc += nioprocs) {
+          string const name =
+            create_filename(cctkGH, basename, cctkGH->cctk_iteration, ioproc,
+                            in_recovery ? io_dir_recover : io_dir_input, false);
+          
+          bool file_exists;
+          H5E_BEGIN_TRY {
+            file_exists = H5Fis_hdf5(name.c_str()) > 0;
+          } H5E_END_TRY;
+          if (not file_exists) {
+            notfoundproc = ioproc;
+            break;
+          }
+          foundproc = ioproc;
+          
+          indent_t indent;
+          cout << indent << "I/O process=" << ioproc << "\n";
+          
+          hid_t const file = H5Fopen(name.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
+          assert(file >= 0);
+          
+          // Iterate over all time slices
+          bool const input_past_timelevels = in_recovery;
+          // TODO: read metadata when recoverying parameters
+          bool const input_metadata = false;
+          input(cctkGH, file, input_var, input_past_timelevels, input_metadata,
+                scatter);
+          
+          // Close file
+          herr = H5Fclose(file);
+          assert(not herr);
         }
-        foundproc = proc;
-        
-        indent_t indent;
-        cout << indent << "process=" << proc << "\n";
-        
-        hid_t const file = H5Fopen(name.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
-        assert(file >= 0);
-        
-        // Iterate over all time slices
-        bool const input_past_timelevels = in_recovery;
-        // TODO: read metadata when recoverying parameters
-        bool const input_metadata = false;
-        input(cctkGH, file, input_var, input_past_timelevels, input_metadata,
-              scatter);
-        
-        // Close file
-        herr = H5Fclose(file);
-        assert(not herr);
       }
       
       {
@@ -526,13 +569,15 @@ namespace CarpetIOF5 {
                      "Could not read input file \"%s\"", name.c_str());
           return 1;
         }
-        if (notfoundproc <= maxfoundproc) {
+        if (notfoundproc > -1 and notfoundproc <= maxfoundproc) {
           CCTK_VWarn(CCTK_WARN_ALERT, __LINE__, __FILE__, CCTK_THORNSTRING,
                      "Could not read file of process %d (but could read file of process %d)",
                      notfoundproc, maxfoundproc);
         }
       }
       
+      // The scatter object is implicitly destroyed here
+      
     } END_GLOBAL_MODE;
     
     return 0;                   // no error