Optimize HDF5 I/O performance by using collective I/O operations

CHANGELOG.md

@@ -1,5 +1,9 @@
 # FANS Changelog
 
+## latest
+
+- Optimize HDF5 I/O performance by using collective I/O operations [#105](https://github.com/DataAnalyticsEngineering/FANS/pull/105)
+
 ## v0.5.0
 
 - Add explicit FE types (HEX8, HEX8R, BBAR) to JSON input [#96](https://github.com/DataAnalyticsEngineering/FANS/pull/96)

CMakeLists.txt

@@ -88,8 +88,8 @@ if ("$ENV{SETVARS_COMPLETED}" STREQUAL "1")
     )
 endif ()
 
-set(HDF5_ENABLE_PARALLEL ON)
-find_package(HDF5 REQUIRED COMPONENTS CXX)
+set(HDF5_PREFER_PARALLEL TRUE)
+find_package(HDF5 REQUIRED COMPONENTS C CXX)
 
 find_package(Eigen3 REQUIRED)
 
@@ -198,6 +198,7 @@ target_link_libraries(FANS_FANS PRIVATE m)
 # call to target_link_libraries(FANS_FANS PUBLIC HDF5::HDF5). But CMake 3.16
 # does not yet support this.
 target_include_directories(FANS_FANS PRIVATE ${HDF5_INCLUDE_DIRS})
+target_link_libraries(FANS_FANS PUBLIC ${HDF5_C_LIBRARIES})
 target_link_libraries(FANS_FANS PUBLIC ${HDF5_CXX_LIBRARIES})
 target_compile_definitions(FANS_FANS PUBLIC ${HDF5_DEFINITIONS})
 if (HDF5_IS_PARALLEL)

README.md

@@ -33,6 +33,7 @@ FANS has the following dependencies:
   ```bash
   apt-get install \
       libhdf5-dev \
+      libhdf5-openmpi-dev \
       libopenmpi-dev \
       libeigen3-dev \
       libfftw3-dev \

docker/Dockerfile

@@ -36,6 +36,7 @@ RUN apt-get update -qq && apt-get install -y --no-install-recommends \
     file \
 # FANS dependencies \
     libhdf5-dev \
+    libhdf5-openmpi-dev \
     libopenmpi-dev \
     libeigen3-dev \
     libfftw3-dev \

include/material_models/GBDiffusion.h

@@ -174,14 +174,7 @@ class GBDiffusion : public ThermalModel, public LinearModel<1, 3> {
                     GBnormals_field[element_idx * 3 + 2] = GBnormals[3 * mat_index + 2];
                 }
             }
-            for (int i = 0; i < solver.world_size; ++i) {
-                if (i == solver.world_rank) {
-                    char name[5096];
-                    sprintf(name, "%s/load%i/time_step%i/GBnormals", solver.dataset_name, load_idx, time_idx);
-                    reader.WriteSlab<double>(GBnormals_field, 3, resultsFileName, name);
-                }
-                MPI_Barrier(MPI_COMM_WORLD);
-            }
+            reader.writeSlab("GBnormals", resultsFileName, solver.dataset_name, load_idx, time_idx, GBnormals_field, 3);
             FANS_free(GBnormals_field);
         }
     }

include/material_models/J2Plasticity.h

@@ -251,38 +251,9 @@ inline void J2Plasticity::postprocess(Solver<3, 6> &solver, Reader &reader, cons
         mean_kinematic_hardening_variable.segment(n_str * elem_idx, n_str) = psi_bar_t[elem_idx].rowwise().mean();
     }
 
-    if (find(reader.resultsToWrite.begin(), reader.resultsToWrite.end(), "plastic_strain") != reader.resultsToWrite.end()) {
-        for (int i = 0; i < solver.world_size; ++i) {
-            if (i == solver.world_rank) {
-                char name[5096];
-                sprintf(name, "%s/load%i/time_step%i/plastic_strain", solver.dataset_name, load_idx, time_idx);
-                reader.WriteSlab<double>(mean_plastic_strain.data(), n_str, resultsFileName, name);
-            }
-            MPI_Barrier(MPI_COMM_WORLD);
-        }
-    }
-
-    if (find(reader.resultsToWrite.begin(), reader.resultsToWrite.end(), "isotropic_hardening_variable") != reader.resultsToWrite.end()) {
-        for (int i = 0; i < solver.world_size; ++i) {
-            if (i == solver.world_rank) {
-                char name[5096];
-                sprintf(name, "%s/load%i/time_step%i/isotropic_hardening_variable", solver.dataset_name, load_idx, time_idx);
-                reader.WriteSlab<double>(mean_isotropic_hardening_variable.data(), 1, resultsFileName, name);
-            }
-            MPI_Barrier(MPI_COMM_WORLD);
-        }
-    }
-
-    if (find(reader.resultsToWrite.begin(), reader.resultsToWrite.end(), "kinematic_hardening_variable") != reader.resultsToWrite.end()) {
-        for (int i = 0; i < solver.world_size; ++i) {
-            if (i == solver.world_rank) {
-                char name[5096];
-                sprintf(name, "%s/load%i/time_step%i/kinematic_hardening_variable", solver.dataset_name, load_idx, time_idx);
-                reader.WriteSlab<double>(mean_kinematic_hardening_variable.data(), n_str, resultsFileName, name);
-            }
-            MPI_Barrier(MPI_COMM_WORLD);
-        }
-    }
+    reader.writeSlab("plastic_strain", resultsFileName, solver.dataset_name, load_idx, time_idx, mean_plastic_strain.data(), n_str);
+    reader.writeSlab("isotropic_hardening_variable", resultsFileName, solver.dataset_name, load_idx, time_idx, mean_isotropic_hardening_variable.data(), 1);
+    reader.writeSlab("kinematic_hardening_variable", resultsFileName, solver.dataset_name, load_idx, time_idx, mean_kinematic_hardening_variable.data(), n_str);
 }
 
 #endif // J2PLASTICITY_H

include/material_models/PseudoPlastic.h

@@ -54,16 +54,7 @@ class PseudoPlastic : public SmallStrainMechModel {
             element_plastic_flag(elem_idx) = plastic_flag[elem_idx].cast<float>().mean();
         }
 
-        if (find(reader.resultsToWrite.begin(), reader.resultsToWrite.end(), "plastic_flag") != reader.resultsToWrite.end()) {
-            for (int i = 0; i < solver.world_size; ++i) {
-                if (i == solver.world_rank) {
-                    char name[5096];
-                    sprintf(name, "%s/load%i/time_step%i/plastic_flag", solver.dataset_name, load_idx, time_idx);
-                    reader.WriteSlab<float>(element_plastic_flag.data(), 1, resultsFileName, name);
-                }
-                MPI_Barrier(MPI_COMM_WORLD);
-            }
-        }
+        reader.writeSlab("plastic_flag", resultsFileName, solver.dataset_name, load_idx, time_idx, element_plastic_flag.data(), 1);
     }
 
   protected:

include/reader.h

@@ -5,7 +5,11 @@
 #include <map>
 #include <string>
 #include <vector>
+#include <unsupported/Eigen/CXX11/Tensor>
 #include "mixedBCs.h"
+#include "hdf5.h"
+#include "H5FDmpio.h"
+#include "mpi.h"
 
 using namespace std;
 
@@ -58,6 +62,15 @@ class Reader {
     // void ReadHDF5(char file_name[], char dset_name[]);
     void safe_create_group(hid_t file, const char *const name);
 
+    // Convenience methods to check if a result should be written and write it
+    template <typename T>
+    void writeData(const char *fieldName, const char *file_name, const char *dataset_name,
+                   int load_idx, int time_idx, T *data, hsize_t *dims, int ndims);
+
+    template <typename T>
+    void writeSlab(const char *fieldName, const char *file_name, const char *dataset_name,
+                   int load_idx, int time_idx, T *data, int size);
+
     template <typename T>
     void WriteSlab(T *data, int _howmany, const char *file_name, const char *dset_name);
 
@@ -176,7 +189,7 @@ void Reader::WriteSlab(
     /* 1. open or create the HDF5 file                                  */
     /*------------------------------------------------------------------*/
     hid_t plist_id = H5Pcreate(H5P_FILE_ACCESS);
-    // H5Pset_fapl_mpio(plist_id, comm, info);   /* if you need MPI-IO */
+    H5Pset_fapl_mpio(plist_id, MPI_COMM_WORLD, MPI_INFO_NULL);
 
     /* temporarily silence “file not found” during H5Fopen */
     herr_t (*old_func)(hid_t, void *);
@@ -243,43 +256,73 @@ void Reader::WriteSlab(
     /*------------------------------------------------------------------*/
     /* 4. transpose local slab  [X][Y][Z][k]  ->  [Z][Y][X][k]          */
     /*------------------------------------------------------------------*/
-    const size_t NxLoc = static_cast<size_t>(local_n0);
-    const size_t NyLoc = static_cast<size_t>(Ny);
-    const size_t NzLoc = static_cast<size_t>(Nz);
-    const size_t kLoc  = static_cast<size_t>(kDim);
-
-    size_t         slabElems = NxLoc * NyLoc * NzLoc * kLoc;
-    std::vector<T> tmp(slabElems); /* automatic RAII buffer */
-
-    for (size_t x = 0; x < NxLoc; ++x)
-        for (size_t y = 0; y < NyLoc; ++y)
-            for (size_t z = 0; z < NzLoc; ++z) {
-                size_t srcBase = (((x * NyLoc) + y) * NzLoc + z) * kLoc; /* X-major */
-                size_t dstBase = (((z * NyLoc) + y) * NxLoc + x) * kLoc; /* Z-major */
-                std::memcpy(&tmp[dstBase], &data[srcBase], kLoc * sizeof(T));
-            }
+    const Eigen::Index NxLoc     = static_cast<Eigen::Index>(local_n0);
+    const Eigen::Index NyLoc     = static_cast<Eigen::Index>(Ny);
+    const Eigen::Index NzLoc     = static_cast<Eigen::Index>(Nz);
+    const Eigen::Index kLoc      = static_cast<Eigen::Index>(kDim);
+    const size_t       slabElems = static_cast<size_t>(NxLoc * NyLoc * NzLoc * kLoc);
+
+    T *tmp = static_cast<T *>(fftw_malloc(slabElems * sizeof(T)));
+    if (!tmp) {
+        throw std::bad_alloc();
+    }
+
+    Eigen::TensorMap<Eigen::Tensor<const T, 4, Eigen::RowMajor>>
+        input_tensor(data, NxLoc, NyLoc, NzLoc, kLoc);
+
+    Eigen::TensorMap<Eigen::Tensor<T, 4, Eigen::RowMajor>>
+        output_tensor(tmp, NzLoc, NyLoc, NxLoc, kLoc);
+
+    // Permute dimensions: (0,1,2,3) -> (2,1,0,3) swaps X and Z
+    Eigen::array<Eigen::Index, 4> shuffle_dims = {2, 1, 0, 3};
+    output_tensor                              = input_tensor.shuffle(shuffle_dims);
 
     /*------------------------------------------------------------------*/
     /* 5. MEMORY dataspace matches FILE slab exactly                    */
     /*------------------------------------------------------------------*/
     hid_t memspace = H5Screate_simple(4, fcount, nullptr);
 
     plist_id = H5Pcreate(H5P_DATASET_XFER);
-    // H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE);
+    H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE);
 
     herr_t status = H5Dwrite(dset_id, data_type,
-                             memspace, filespace, plist_id, tmp.data());
+                             memspace, filespace, plist_id, tmp);
     if (status < 0)
         throw std::runtime_error("WriteSlab: H5Dwrite failed");
 
     /*------------------------------------------------------------------*/
     /* 6. tidy up                                                       */
     /*------------------------------------------------------------------*/
+    fftw_free(tmp);
     H5Sclose(memspace);
     H5Sclose(filespace);
     H5Pclose(plist_id);
     H5Dclose(dset_id);
     H5Fclose(file_id);
 }
 
+template <typename T>
+void Reader::writeData(const char *fieldName, const char *file_name, const char *dataset_name,
+                       int load_idx, int time_idx, T *data, hsize_t *dims, int ndims)
+{
+    if (world_rank != 0 || std::find(resultsToWrite.begin(), resultsToWrite.end(), fieldName) == resultsToWrite.end()) {
+        return;
+    }
+    char name[5096];
+    snprintf(name, sizeof(name), "%s/load%i/time_step%i/%s", dataset_name, load_idx, time_idx, fieldName);
+    WriteData(data, file_name, name, dims, ndims);
+}
+
+template <typename T>
+void Reader::writeSlab(const char *fieldName, const char *file_name, const char *dataset_name,
+                       int load_idx, int time_idx, T *data, int size)
+{
+    if (std::find(resultsToWrite.begin(), resultsToWrite.end(), fieldName) == resultsToWrite.end()) {
+        return;
+    }
+    char name[5096];
+    snprintf(name, sizeof(name), "%s/load%i/time_step%i/%s", dataset_name, load_idx, time_idx, fieldName);
+    WriteSlab(data, size, file_name, name);
+}
+
 #endif

include/solver.h

@@ -567,62 +567,40 @@ void Solver<howmany, n_str>::postprocess(Reader &reader, const char resultsFileN
         }
     }
 
-    // Write results to results h5 file
-    auto writeData = [&](const char *resultName, const char *resultPrefix, auto *data, hsize_t *dims, int ndims) {
-        if (std::find(reader.resultsToWrite.begin(), reader.resultsToWrite.end(), resultName) != reader.resultsToWrite.end()) {
-            char name[5096];
-            sprintf(name, "%s/load%i/time_step%i/%s", dataset_name, load_idx, time_idx, resultPrefix);
-            reader.WriteData(data, resultsFileName, name, dims, ndims);
-        }
-    };
-
-    auto writeSlab = [&](const char *resultName, const char *resultPrefix, auto *data, int size) {
-        if (std::find(reader.resultsToWrite.begin(), reader.resultsToWrite.end(), resultName) != reader.resultsToWrite.end()) {
-            char name[5096];
-            sprintf(name, "%s/load%i/time_step%i/%s", dataset_name, load_idx, time_idx, resultPrefix);
-            reader.WriteSlab(data, size, resultsFileName, name);
-        }
-    };
-
-    for (int i = 0; i < world_size; ++i) {
-        if (i == world_rank) {
-            if (world_rank == 0) {
-                hsize_t dims[1] = {static_cast<hsize_t>(n_str)};
-                writeData("stress_average", "stress_average", stress_average.data(), dims, 1);
-                writeData("strain_average", "strain_average", strain_average.data(), dims, 1);
-
-                for (int mat_index = 0; mat_index < n_mat; ++mat_index) {
-                    char stress_name[512];
-                    char strain_name[512];
-                    sprintf(stress_name, "phase_stress_average_phase%d", mat_index);
-                    sprintf(strain_name, "phase_strain_average_phase%d", mat_index);
-                    writeData("phase_stress_average", stress_name, phase_stress_average[mat_index].data(), dims, 1);
-                    writeData("phase_strain_average", strain_name, phase_strain_average[mat_index].data(), dims, 1);
-                }
-                dims[0] = iter + 1;
-                writeData("absolute_error", "absolute_error", err_all.data(), dims, 1);
-            }
-            writeSlab("microstructure", "microstructure", ms, 1);
-            writeSlab("displacement_fluctuation", "displacement_fluctuation", v_u, howmany);
-            writeSlab("displacement", "displacement", u_total.data(), howmany);
-            writeSlab("residual", "residual", v_r, howmany);
-            writeSlab("strain", "strain", strain.data(), n_str);
-            writeSlab("stress", "stress", stress.data(), n_str);
-        }
-        MPI_Barrier(MPI_COMM_WORLD);
+    hsize_t dims[1] = {static_cast<hsize_t>(n_str)};
+    reader.writeData("stress_average", resultsFileName, dataset_name, load_idx, time_idx, stress_average.data(), dims, 1);
+    reader.writeData("strain_average", resultsFileName, dataset_name, load_idx, time_idx, strain_average.data(), dims, 1);
+    for (int mat_index = 0; mat_index < n_mat; ++mat_index) {
+        char stress_name[512];
+        char strain_name[512];
+        sprintf(stress_name, "phase_stress_average_phase%d", mat_index);
+        sprintf(strain_name, "phase_strain_average_phase%d", mat_index);
+        reader.writeData(stress_name, resultsFileName, dataset_name, load_idx, time_idx, phase_stress_average[mat_index].data(), dims, 1);
+        reader.writeData(strain_name, resultsFileName, dataset_name, load_idx, time_idx, phase_strain_average[mat_index].data(), dims, 1);
     }
+    dims[0] = iter + 1;
+    reader.writeData("absolute_error", resultsFileName, dataset_name, load_idx, time_idx, err_all.data(), dims, 1);
+
+    vector<int> rank_field(local_n0 * n_y * n_z, world_rank);
+    reader.writeSlab("mpi_rank", resultsFileName, dataset_name, load_idx, time_idx, rank_field.data(), 1);
+    reader.writeSlab("microstructure", resultsFileName, dataset_name, load_idx, time_idx, ms, 1);
+    reader.writeSlab("displacement_fluctuation", resultsFileName, dataset_name, load_idx, time_idx, v_u, howmany);
+    reader.writeSlab("displacement", resultsFileName, dataset_name, load_idx, time_idx, u_total.data(), howmany);
+    reader.writeSlab("residual", resultsFileName, dataset_name, load_idx, time_idx, v_r, howmany);
+    reader.writeSlab("strain", resultsFileName, dataset_name, load_idx, time_idx, strain.data(), n_str);
+    reader.writeSlab("stress", resultsFileName, dataset_name, load_idx, time_idx, stress.data(), n_str);
 
     matmodel->postprocess(*this, reader, resultsFileName, load_idx, time_idx);
 
-    // Compute homogenized tangent
+    // Compute homogenized tangent only if requested
     if (find(reader.resultsToWrite.begin(), reader.resultsToWrite.end(), "homogenized_tangent") != reader.resultsToWrite.end()) {
         homogenized_tangent = get_homogenized_tangent(1e-6);
         hsize_t dims[2]     = {static_cast<hsize_t>(n_str), static_cast<hsize_t>(n_str)};
         if (world_rank == 0) {
             cout << "# Homogenized tangent: " << endl
                  << setprecision(12) << homogenized_tangent << endl
                  << endl;
-            writeData("homogenized_tangent", "homogenized_tangent", homogenized_tangent.data(), dims, 2);
+            reader.writeData("homogenized_tangent", resultsFileName, dataset_name, load_idx, time_idx, homogenized_tangent.data(), dims, 2);
         }
     }
 }

src/reader.cpp

@@ -393,14 +393,15 @@ void Reader ::ReadMS(int hm)
                                      static_cast<size_t>(dims[2]));
 
     if (is_zyx) {
-        /* tmp =  [z][y][x] , we need ms = [x][y][z] */
-        for (size_t z = 0; z < dims[2]; ++z)
-            for (size_t y = 0; y < dims[1]; ++y)
-                for (size_t x = 0; x < static_cast<size_t>(local_n0); ++x) {
-                    size_t idx_tmp = (z * dims[1] + y) * local_n0 + x; // z-major
-                    size_t idx_ms  = (x * dims[1] + y) * dims[2] + z;  // x-major
-                    ms[idx_ms]     = tmp[idx_tmp];
-                }
+        const Eigen::Index Nx = static_cast<Eigen::Index>(local_n0);
+        const Eigen::Index Ny = static_cast<Eigen::Index>(dims[1]);
+        const Eigen::Index Nz = static_cast<Eigen::Index>(dims[2]);
+
+        Eigen::TensorMap<Eigen::Tensor<const unsigned short, 3, Eigen::RowMajor>>
+            input_tensor(tmp, Nz, Ny, Nx); // [Z][Y][X] in file
+        Eigen::TensorMap<Eigen::Tensor<unsigned short, 3, Eigen::RowMajor>>
+            output_tensor(ms, Nx, Ny, Nz); // [X][Y][Z] in memory
+        output_tensor = input_tensor.shuffle(Eigen::array<Eigen::Index, 3>{2, 1, 0});
         FANS_free(tmp);
     } else {
         /* XYZ case: the slab is already in correct order */