hdf5 halo fix for mpi

igfuw · Feb 10, 2023 · 484eaf3 · 484eaf3
1 parent 8b21e92
commit 484eaf3
Showing 1 changed file with 18 additions and 13 deletions.
diff --git a/libmpdata++/output/hdf5.hpp b/libmpdata++/output/hdf5.hpp
@@ -55,10 +55,10 @@ namespace libmpdataxx
               H5::PredType::NATIVE_FLOAT,
         flttype_output = H5::PredType::NATIVE_FLOAT; // using floats not to waste disk space
 
-      blitz::TinyVector<hsize_t, parent_t::n_dims> cshape, shape, chunk, srfcshape, srfcchunk, offst, shape_h, chunk_h, offst_h;
+      blitz::TinyVector<hsize_t, parent_t::n_dims> cshape, shape, chunk, srfcshape, srfcchunk, offst, shape_h, chunk_h, offst_h, shape_mem_h, offst_mem_h;
       H5::DSetCreatPropList params;
 
-      H5::DataSpace sspace, cspace, srfcspace, sspace_h;
+      H5::DataSpace sspace, cspace, srfcspace, sspace_h, sspace_mem_h;
 #if defined(USE_MPI)
       hid_t fapl_id;
 #endif
@@ -97,11 +97,13 @@ namespace libmpdataxx
           // x,y,z
           offst   = 0;
           offst_h = 0;
+          offst_mem_h = 0;
 
           for (int d = 0; d < parent_t::n_dims; ++d)
           {
-            shape[d]   = this->mem->distmem.grid_size[d];
-            shape_h[d] = this->mem->distmem.grid_size[d] + 2 * this->halo;
+            shape[d]   = this->mem->distmem.grid_size[d];                        // shape of arrays stored in file
+            shape_h[d] = this->mem->distmem.grid_size[d] + 2 * this->halo;       // shape of arrays with halos stored in files
+            shape_mem_h[d] = this->mem->grid_size[d].length() + 2 * this->halo;  // shape of the array with halo stored in memory of given MPI rank
           }
 
           chunk   = shape;
@@ -113,10 +115,11 @@ namespace libmpdataxx
           srfcshape = shape;
           *(srfcshape.end()-1) = 1;
 
-          sspace    = H5::DataSpace(parent_t::n_dims, shape.data());
-          sspace_h  = H5::DataSpace(parent_t::n_dims, shape_h.data());
-          srfcspace = H5::DataSpace(parent_t::n_dims, srfcshape.data());
-          cspace    = H5::DataSpace(parent_t::n_dims, cshape.data());
+          sspace        = H5::DataSpace(parent_t::n_dims, shape.data());
+          sspace_h      = H5::DataSpace(parent_t::n_dims, shape_h.data());
+          sspace_mem_h  = H5::DataSpace(parent_t::n_dims, shape_mem_h.data());
+          srfcspace     = H5::DataSpace(parent_t::n_dims, srfcshape.data());
+          cspace        = H5::DataSpace(parent_t::n_dims, cshape.data());
 
 #if defined(USE_MPI)
           if (this->mem->distmem.size() > 1)
@@ -133,8 +136,10 @@ namespace libmpdataxx
             if (this->mem->distmem.rank() == this->mem->distmem.size() - 1)
               cshape[0] += 1;
 
-            offst[0]   = this->mem->grid_size[0].first();
-            offst_h[0] = this->mem->distmem.rank() == 0 ? 0 : this->mem->grid_size[0].first() + this->halo;
+            offst[0]       = this->mem->grid_size[0].first();
+            offst_h[0]     = this->mem->distmem.rank() == 0 ? 0 : this->mem->grid_size[0].first() + this->halo;
+            if (this->mem->distmem.rank() > 0)
+              offst_mem_h[0] = this->halo;
 
             // chunk size has to be common to all processes !
             // TODO: set to 1? Test performance...
@@ -402,11 +407,11 @@ namespace libmpdataxx
         // revert to default chunk
         params.setChunk(parent_t::n_dims, chunk.data());
 
-        record_dsc_helper(aux, arr);
-
         auto space = aux.getSpace();
         space.selectHyperslab(H5S_SELECT_SET, shape_h.data(), offst_h.data());
-        aux.write(arr.data(), flttype_solver, H5::DataSpace(parent_t::n_dims, shape_h.data()), space, dxpl_id);
+        sspace_mem_h.selectHyperslab(H5S_SELECT_SET, shape_h.data(), offst_mem_h.data());
+
+        aux.write(arr.data(), flttype_solver, sspace_mem_h, space, dxpl_id);
       }
 
       void record_scalar_hlpr(const std::string &name, const std::string &group_name, typename solver_t::real_t data, H5::H5File hdf)