From d149fade99ebd1c340784deb4b2cf2cfbf6600cd Mon Sep 17 00:00:00 2001 From: "Todd A. Oliver" Date: Mon, 29 Jan 2024 18:33:32 -0600 Subject: [PATCH 01/14] Style --- src/cycle_avg_joule_coupling.cpp | 1 - src/tps2Boltzmann.cpp | 5 +++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/cycle_avg_joule_coupling.cpp b/src/cycle_avg_joule_coupling.cpp index f5f5c3a5c..acf78235b 100644 --- a/src/cycle_avg_joule_coupling.cpp +++ b/src/cycle_avg_joule_coupling.cpp @@ -313,7 +313,6 @@ void CycleAvgJouleCoupling::interpElectricFieldFromEMToFlow() { efieldR_->SetFromTrueDofs(interp_vals); efieldR_->HostRead(); - const ParGridFunction *efield_imag_gf = qmsa_solver_->getElectricFieldimag(); interp_em_to_flow_->Interpolate(vxyz, *efield_imag_gf, interp_vals); efieldI_->SetFromTrueDofs(interp_vals); diff --git a/src/tps2Boltzmann.cpp b/src/tps2Boltzmann.cpp index 171313f0d..59f0d18ca 100644 --- a/src/tps2Boltzmann.cpp +++ b/src/tps2Boltzmann.cpp @@ -41,9 +41,10 @@ #include #endif +#include + #include #include -#include namespace TPS { @@ -85,7 +86,7 @@ Tps2Boltzmann::Tps2Boltzmann(Tps *tps) : NIndexes(7), tps_(tps), all_fes_(nullpt assert(basis_type_ == 0 || basis_type_ == 1); tps->getRequiredInput("em/current_frequency", EfieldAngularFreq_); - EfieldAngularFreq_ *= 2.*M_PI; + EfieldAngularFreq_ *= 2. * M_PI; offsets.SetSize(NIndexes + 1); ncomps.SetSize(NIndexes + 1); From 2b25303aae141bea6a2925e93758643877523a1d Mon Sep 17 00:00:00 2001 From: "Todd A. Oliver" Date: Mon, 29 Jan 2024 19:10:06 -0600 Subject: [PATCH 02/14] Eliminate M2ulPhyS::writeHistoryFile Based on kevin's comments, we don't use or maintain this anymore, so better get rid of it. --- src/M2ulPhyS.cpp | 6 ------ src/M2ulPhyS.hpp | 5 ----- src/io.cpp | 25 ------------------------- 3 files changed, 36 deletions(-) diff --git a/src/M2ulPhyS.cpp b/src/M2ulPhyS.cpp index 822a2eed9..4e9b6d687 100644 --- a/src/M2ulPhyS.cpp +++ b/src/M2ulPhyS.cpp @@ -1988,12 +1988,6 @@ void M2ulPhyS::solveStep() { const int vis_steps = config.GetNumItersOutput(); if (iter % vis_steps == 0) { - // dump history - // NOTE(kevin): this routine is currently obsolete. - // It computes `dof`-averaged state and time-derivative, which are useless at this point. - // This will not be supported. - writeHistoryFile(); - #ifdef HAVE_MASA if (config.use_mms_) { if (config.mmsSaveDetails_) { diff --git a/src/M2ulPhyS.hpp b/src/M2ulPhyS.hpp index 934131487..5bfb3f33a 100644 --- a/src/M2ulPhyS.hpp +++ b/src/M2ulPhyS.hpp @@ -343,11 +343,6 @@ class M2ulPhyS : public TPS::Solver { void initGradUp(); void initilizeSpeciesFromLTE(); - // NOTE(kevin): this routine is currently obsolete. - // It computes `dof`-averaged state and time-derivative, which are useless at this point. - // This will not be supported. - void writeHistoryFile(); - // i/o routines void read_partitioned_soln_data(hid_t file, string varName, size_t index, double *data); void read_serialized_soln_data(hid_t file, string varName, int numDof, int varOffset, double *data, IOFamily &fam); diff --git a/src/io.cpp b/src/io.cpp index b2cd3032e..e6f98f017 100644 --- a/src/io.cpp +++ b/src/io.cpp @@ -669,31 +669,6 @@ void M2ulPhyS::write_soln_data(hid_t group, string varName, hid_t dataspace, dou return; } -void M2ulPhyS::writeHistoryFile() { - MPI_Comm TPSCommWorld = this->groupsMPI->getTPSCommWorld(); - double global_dUdt[5]; - MPI_Allreduce(rhsOperator->getLocalTimeDerivatives(), &global_dUdt, 5, MPI_DOUBLE, MPI_SUM, TPSCommWorld); - - if (rank0_) { - histFile << time << "," << iter; - for (int eq = 0; eq < 5; eq++) { - histFile << "," << global_dUdt[eq] / static_cast(nprocs_); - } - } - - if (average->ComputeMean()) { - double global_meanData[5 + 6]; - MPI_Allreduce(average->getLocalSums(), &global_meanData, 5 + 6, MPI_DOUBLE, MPI_SUM, TPSCommWorld); - - if (rank0_) { - histFile << "," << average->GetSamplesMean(); - for (int n = 0; n < 5 + 6; n++) histFile << "," << global_meanData[n] / static_cast(nprocs_); - } - } - - if (rank0_) histFile << endl; -} - void M2ulPhyS::readTable(const std::string &inputPath, TableInput &result) { MPI_Comm TPSCommWorld = this->groupsMPI->getTPSCommWorld(); tpsP->getInput((inputPath + "/x_log").c_str(), result.xLogScale, false); From 4335a39f02d042902bbe133bd5f194e84d101b60 Mon Sep 17 00:00:00 2001 From: "Todd A. Oliver" Date: Mon, 29 Jan 2024 18:34:45 -0600 Subject: [PATCH 03/14] Move some M2ulPhyS IO support into standalone fcns (#244) On the path to making the IO capabilities more reusable, in this commit we extract 5 functions out of M2ulPhyS and into standalone functions. The refactored functions are write_soln_data read_partitioned_soln_data read_serialized_soln_data partitioning_file_hdf5 serialize_soln_for_write --- src/M2ulPhyS.cpp | 6 +- src/M2ulPhyS.hpp | 5 -- src/io.cpp | 132 ++++++++++++++++++-------------------- src/io.hpp | 13 ++++ src/run_configuration.hpp | 2 +- 5 files changed, 81 insertions(+), 77 deletions(-) diff --git a/src/M2ulPhyS.cpp b/src/M2ulPhyS.cpp index 4e9b6d687..de7be70ec 100644 --- a/src/M2ulPhyS.cpp +++ b/src/M2ulPhyS.cpp @@ -329,9 +329,9 @@ void M2ulPhyS::initVariables() { if (config.isRestartSerialized("read")) { assert(serial_mesh->Conforming()); partitioning_ = Array(serial_mesh->GeneratePartitioning(nprocs_, defaultPartMethod), nelemGlobal_); - partitioning_file_hdf5("write"); + partitioning_file_hdf5("write", config, groupsMPI, nelemGlobal_, partitioning_); } else { - partitioning_file_hdf5("read"); + partitioning_file_hdf5("read", config, groupsMPI, nelemGlobal_, partitioning_); } } @@ -359,7 +359,7 @@ void M2ulPhyS::initVariables() { if (nprocs_ > 1) { assert(serial_mesh->Conforming()); partitioning_ = Array(serial_mesh->GeneratePartitioning(nprocs_, defaultPartMethod), nelemGlobal_); - if (rank0_) partitioning_file_hdf5("write"); + if (rank0_) partitioning_file_hdf5("write", config, groupsMPI, nelemGlobal_, partitioning_); MPI_Barrier(groupsMPI->getTPSCommWorld()); } diff --git a/src/M2ulPhyS.hpp b/src/M2ulPhyS.hpp index 5bfb3f33a..05e47db78 100644 --- a/src/M2ulPhyS.hpp +++ b/src/M2ulPhyS.hpp @@ -344,12 +344,7 @@ class M2ulPhyS : public TPS::Solver { void initilizeSpeciesFromLTE(); // i/o routines - void read_partitioned_soln_data(hid_t file, string varName, size_t index, double *data); - void read_serialized_soln_data(hid_t file, string varName, int numDof, int varOffset, double *data, IOFamily &fam); void restart_files_hdf5(string mode, string inputFileName = std::string()); - void partitioning_file_hdf5(string mode); - void serialize_soln_for_write(IOFamily &fam); - void write_soln_data(hid_t group, string varName, hid_t dataspace, double *data); void Check_NAN(); bool Check_JobResubmit(); diff --git a/src/io.cpp b/src/io.cpp index e6f98f017..c97a40f7e 100644 --- a/src/io.cpp +++ b/src/io.cpp @@ -265,7 +265,7 @@ void M2ulPhyS::restart_files_hdf5(string mode, string inputFileName) { double *data = fam.pfunc_->HostReadWrite(); // special case if writing a serial restart if ((config.isRestartSerialized(mode)) && (nprocs_ > 1)) { - serialize_soln_for_write(fam); + serialize_soln_for_write(fam, groupsMPI, mesh->GetNE(), nelemGlobal_, locToGlobElem, partitioning_); if (rank0_) data = fam.serial_sol->HostReadWrite(); } @@ -274,7 +274,7 @@ void M2ulPhyS::restart_files_hdf5(string mode, string inputFileName) { // save raw data if (rank0_ || config.isRestartPartitioned(mode)) { - for (auto var : vars) write_soln_data(group, var.varName_, dataspace, data + var.index_ * dims[0]); + for (auto var : vars) write_soln_data(group, var.varName_, dataspace, data + var.index_ * dims[0], rank0_); } if (group >= 0) H5Gclose(group); @@ -326,10 +326,14 @@ void M2ulPhyS::restart_files_hdf5(string mode, string inputFileName) { if (var.inRestartFile_) { std::string h5Path = fam.group_ + "/" + var.varName_; if (rank0_) grvy_printf(ginfo, "--> Reading h5 path = %s\n", h5Path.c_str()); - if (config.isRestartPartitioned(mode)) + if (config.isRestartPartitioned(mode)) { read_partitioned_soln_data(file, h5Path.c_str(), var.index_ * numInSoln, data); - else - read_serialized_soln_data(file, h5Path.c_str(), dof, var.index_, data, fam); + } else { + assert(partitioning_ != NULL); + assert(config.isRestartSerialized("read")); + read_serialized_soln_data(file, h5Path.c_str(), dof, var.index_, data, fam, groupsMPI, partitioning_, + nelemGlobal_); + } } } } @@ -381,23 +385,27 @@ void M2ulPhyS::restart_files_hdf5(string mode, string inputFileName) { return; } -void M2ulPhyS::partitioning_file_hdf5(std::string mode) { +void partitioning_file_hdf5(std::string mode, const RunConfiguration &config, MPI_Groups *groupsMPI, int nelemGlobal, + Array &partitioning) { MPI_Comm TPSCommWorld = groupsMPI->getTPSCommWorld(); + const bool rank0 = groupsMPI->isWorldRoot(); + const int nprocs = groupsMPI->getTPSWorldSize(); + grvy_timer_begin(__func__); // only rank 0 writes partitioning file - if (!rank0_ && (mode == "write")) return; + if (!rank0 && (mode == "write")) return; // hid_t file, dataspace, data_soln; hid_t file = -1, dataspace; herr_t status; std::string fileName = config.GetPartitionBaseName(); - fileName += "." + std::to_string(nprocs_) + "p.h5"; + fileName += "." + std::to_string(nprocs) + "p.h5"; assert((mode == "read") || (mode == "write")); if (mode == "write") { - assert(partitioning_.Size() > 0); + assert(partitioning.Size() > 0); if (file_exists(fileName)) { grvy_printf(gwarn, "Removing existing partition file: %s\n", fileName.c_str()); @@ -421,26 +429,26 @@ void M2ulPhyS::partitioning_file_hdf5(std::string mode) { if (mode == "write") { // Attributes - h5_save_attribute(file, "numProcs", nprocs_); + h5_save_attribute(file, "numProcs", nprocs); // Raw partition info hsize_t dims[1]; hid_t data; - dims[0] = partitioning_.Size(); + dims[0] = partitioning.Size(); dataspace = H5Screate_simple(1, dims, NULL); assert(dataspace >= 0); data = H5Dcreate2(file, "partitioning", H5T_NATIVE_INT, dataspace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); assert(data >= 0); - status = H5Dwrite(data, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, partitioning_.GetData()); + status = H5Dwrite(data, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, partitioning.HostRead()); assert(status >= 0); H5Dclose(data); H5Fclose(file); } else if (mode == "read") { - partitioning_.SetSize(nelemGlobal_); + partitioning.SetSize(nelemGlobal); - if (rank0_) { + if (rank0) { grvy_printf(ginfo, "Reading original domain decomposition partition file: %s\n", fileName.c_str()); // verify partition info matches current proc count @@ -448,8 +456,8 @@ void M2ulPhyS::partitioning_file_hdf5(std::string mode) { int np = 0; h5_read_attribute(file, "numProcs", np); grvy_printf(ginfo, "--> # partitions defined = %i\n", np); - if (np != nprocs_) { - grvy_printf(gerror, "[ERROR]: Partition info does not match current processor count -> %i\n", nprocs_); + if (np != nprocs) { + grvy_printf(gerror, "[ERROR]: Partition info does not match current processor count -> %i\n", nprocs); exit(ERROR); } } @@ -462,55 +470,50 @@ void M2ulPhyS::partitioning_file_hdf5(std::string mode) { dataspace = H5Dget_space(data); numInFile = H5Sget_simple_extent_npoints(dataspace); - assert((int)numInFile == nelemGlobal_); + assert((int)numInFile == nelemGlobal); - status = H5Dread(data, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, partitioning_.GetData()); + status = H5Dread(data, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, partitioning.HostWrite()); assert(status >= 0); H5Dclose(data); } // <-- end rank0_ #if 0 // distribute partition vector to all procs - MPI_Bcast(partitioning_.GetData(), nelemGlobal_, MPI_INT, 0, TPSCommWorld); + MPI_Bcast(partitioning_.GetData(), nelemGlobal, MPI_INT, 0, TPSCommWorld); #endif // distribute partition vector to all procs (serialzed per process variant) int tag = 21; - if (rank0_) { - for (int rank = 1; rank < nprocs_; rank++) { - MPI_Send(partitioning_.GetData(), nelemGlobal_, MPI_INT, rank, tag, TPSCommWorld); + if (rank0) { + for (int rank = 1; rank < nprocs; rank++) { + MPI_Send(partitioning.HostRead(), nelemGlobal, MPI_INT, rank, tag, TPSCommWorld); grvy_printf(gdebug, "Sent partitioning data to rank %i\n", rank); } } else { - MPI_Recv(partitioning_.GetData(), nelemGlobal_, MPI_INT, 0, tag, TPSCommWorld, MPI_STATUS_IGNORE); + MPI_Recv(partitioning.HostWrite(), nelemGlobal, MPI_INT, 0, tag, TPSCommWorld, MPI_STATUS_IGNORE); } - if (rank0_) grvy_printf(ginfo, "--> partition file read complete\n"); + if (rank0) grvy_printf(ginfo, "--> partition file read complete\n"); } grvy_timer_end(__func__); } -void M2ulPhyS::serialize_soln_for_write(IOFamily &fam) { +void serialize_soln_for_write(IOFamily &fam, MPI_Groups *groupsMPI, int local_ne, int global_ne, + const int *locToGlobElem, const Array &partitioning) { MPI_Comm TPSCommWorld = groupsMPI->getTPSCommWorld(); - // Get total number of elements - const int local_ne = mesh->GetNE(); - int global_ne; - MPI_Reduce(&local_ne, &global_ne, 1, MPI_INT, MPI_SUM, 0, TPSCommWorld); - - // ks note: need to double check this routine when additional solution - // families are added and then remove next assert - // assert(ioFamily == "/solution"); - - // get pargridfunction for this IO family - // int iFamily = ioData.getIOFamilyIndex(ioFamily); - // assert(iFamily >= 0); - ParGridFunction *pfunc = fam.pfunc_; - - if (rank0_) assert(global_ne == nelemGlobal_); - - assert((locToGlobElem != NULL) && (partitioning_ != NULL)); + const bool rank0 = groupsMPI->isWorldRoot(); + + // Ensure consistency + int global_ne_check; + MPI_Reduce(&local_ne, &global_ne_check, 1, MPI_INT, MPI_SUM, 0, TPSCommWorld); + if (rank0) { + assert(global_ne_check == global_ne); + assert(partitioning.Size() == global_ne); + } + assert(locToGlobElem != NULL); - if (rank0_) { + ParGridFunction *pfunc = fam.pfunc_; + if (rank0) { grvy_printf(ginfo, "Generating serialized restart file (group %s...)\n", fam.group_.c_str()); // copy my own data Array lvdofs, gvdofs; @@ -525,7 +528,7 @@ void M2ulPhyS::serialize_soln_for_write(IOFamily &fam) { // have rank 0 receive data from other tasks and copy its own for (int gelem = 0; gelem < global_ne; gelem++) { - int from_rank = partitioning_[gelem]; + int from_rank = partitioning[gelem]; if (from_rank != 0) { fam.serial_fes->GetElementVDofs(gelem, gvdofs); lsoln.SetSize(gvdofs.Size()); @@ -553,9 +556,7 @@ void M2ulPhyS::serialize_soln_for_write(IOFamily &fam) { } // end function: serialize_soln_for_write() // convenience function to read solution data for parallel restarts -void M2ulPhyS::read_partitioned_soln_data(hid_t file, string varName, size_t index, double *data) { - assert(config.isRestartPartitioned("read")); - +void read_partitioned_soln_data(hid_t file, string varName, size_t index, double *data) { hid_t data_soln; herr_t status; @@ -567,10 +568,12 @@ void M2ulPhyS::read_partitioned_soln_data(hid_t file, string varName, size_t ind } // convenience function to read and distribute solution data for serialized restarts -void M2ulPhyS::read_serialized_soln_data(hid_t file, string varName, int numDof, int varOffset, double *data, - IOFamily &fam) { - MPI_Comm TPSCommWorld = this->groupsMPI->getTPSCommWorld(); - assert(config.isRestartSerialized("read")); +void read_serialized_soln_data(hid_t file, string varName, int numDof, int varOffset, double *data, IOFamily &fam, + MPI_Groups *groupsMPI, Array partitioning, int nelemGlobal) { + MPI_Comm TPSCommWorld = groupsMPI->getTPSCommWorld(); + const bool rank0 = groupsMPI->isWorldRoot(); + const int nprocs = groupsMPI->getTPSWorldSize(); + const int rank = groupsMPI->getTPSWorldRank(); hid_t data_soln; herr_t status; @@ -580,13 +583,7 @@ void M2ulPhyS::read_serialized_soln_data(hid_t file, string varName, int numDof, const int tag = 20; - // assert( (dim == 2) || (dim == 3) ); - // if(dim == 2) - // numStateVars = 4; - // else - // numStateVars = 5; - - if (rank0_) { + if (rank0) { grvy_printf(ginfo, "[RestartSerial]: Reading %s for distribution\n", varName.c_str()); Vector data_serial; @@ -598,15 +595,15 @@ void M2ulPhyS::read_serialized_soln_data(hid_t file, string varName, int numDof, H5Dclose(data_soln); // assign solution owned by rank 0 - assert(partitioning_ != NULL); + assert(partitioning != NULL); Array lvdofs, gvdofs; Vector lnodes; int counter = 0; // int ndof_per_elem; - for (int gelem = 0; gelem < nelemGlobal_; gelem++) { - if (partitioning_[gelem] == 0) { + for (int gelem = 0; gelem < nelemGlobal; gelem++) { + if (partitioning[gelem] == 0) { // cull out subset of local vdof vector to use for this solution var fam.pfunc_->ParFESpace()->GetElementVDofs(counter, lvdofs); int numDof_per_this_elem = lvdofs.Size() / numStateVars; @@ -625,10 +622,10 @@ void M2ulPhyS::read_serialized_soln_data(hid_t file, string varName, int numDof, } // pack remaining data and send to other processors - for (int rank = 1; rank < nprocs_; rank++) { + for (int rank = 1; rank < nprocs; rank++) { std::vector packedData; - for (int gelem = 0; gelem < nelemGlobal_; gelem++) - if (partitioning_[gelem] == rank) { + for (int gelem = 0; gelem < nelemGlobal; gelem++) + if (partitioning[gelem] == rank) { fam.serial_fes->GetElementVDofs(gelem, gvdofs); int numDof_per_this_elem = gvdofs.Size() / numStateVars; for (int i = 0; i < numDof_per_this_elem; i++) packedData.push_back(data_serial[gvdofs[i]]); @@ -638,8 +635,7 @@ void M2ulPhyS::read_serialized_soln_data(hid_t file, string varName, int numDof, } } else { // <-- end rank 0 int numlDofs = fam.pfunc_->ParFESpace()->GetNDofs(); - grvy_printf(gdebug, "[%i]: local number of state vars to receive = %i (var=%s)\n", rank_, numlDofs, - varName.c_str()); + grvy_printf(gdebug, "[%i]: local number of state vars to receive = %i (var=%s)\n", rank, numlDofs, varName.c_str()); std::vector packedData(numlDofs); @@ -652,12 +648,12 @@ void M2ulPhyS::read_serialized_soln_data(hid_t file, string varName, int numDof, } // convenience function to write HDF5 data -void M2ulPhyS::write_soln_data(hid_t group, string varName, hid_t dataspace, double *data) { +void write_soln_data(hid_t group, string varName, hid_t dataspace, double *data, bool rank0) { hid_t data_soln; herr_t status; assert(group >= 0); - if (rank0_) grvy_printf(ginfo, " --> Saving (%s)\n", varName.c_str()); + if (rank0) grvy_printf(ginfo, " --> Saving (%s)\n", varName.c_str()); data_soln = H5Dcreate2(group, varName.c_str(), H5T_NATIVE_DOUBLE, dataspace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); assert(data_soln >= 0); diff --git a/src/io.hpp b/src/io.hpp index 2e632dacb..0e7637f12 100644 --- a/src/io.hpp +++ b/src/io.hpp @@ -32,8 +32,13 @@ #ifndef IO_HPP_ #define IO_HPP_ +#include + +#include #include +#include "mpi_groups.hpp" +#include "run_configuration.hpp" #include "tps_mfem_wrap.hpp" class IOFamily { @@ -73,4 +78,12 @@ class IODataOrganizer { void initializeSerial(bool root, bool serial, mfem::Mesh *serial_mesh); }; +void read_partitioned_soln_data(hid_t file, std::string varName, size_t index, double *data); +void read_serialized_soln_data(hid_t file, std::string varName, int numDof, int varOffset, double *data, IOFamily &fam, + MPI_Groups *groupsMPI, mfem::Array partitioning, int nelemGlobal); +void write_soln_data(hid_t group, std::string varName, hid_t dataspace, double *data, bool rank0); +void partitioning_file_hdf5(std::string mode, const RunConfiguration &config, MPI_Groups *groupsMPI, int nelemGlobal, + mfem::Array &partitioning); +void serialize_soln_for_write(IOFamily &fam, MPI_Groups *groupsMPI, int local_ne, int global_ne, + const int *locToGlobElem, const mfem::Array &partitioning); #endif // IO_HPP_ diff --git a/src/run_configuration.hpp b/src/run_configuration.hpp index 6cd5e0e5f..c56359d57 100644 --- a/src/run_configuration.hpp +++ b/src/run_configuration.hpp @@ -298,7 +298,7 @@ class RunConfiguration { string GetMeshFileName() { return meshFile; } string GetOutputName() { return outputFile; } - string GetPartitionBaseName() { return partFile; } + string GetPartitionBaseName() const { return partFile; } int GetUniformRefLevels() { return ref_levels; } int GetTimeIntegratorType() { return timeIntegratorType; } From d6b962cbef27b3e8bb903f7d548e84ebe9ff19b0 Mon Sep 17 00:00:00 2001 From: "Todd A. Oliver" Date: Tue, 30 Jan 2024 10:17:26 -0600 Subject: [PATCH 04/14] Split M2ulPhyS::restart_files_hdf5 into multiple functions (#244) This is an intermediate commit. The idea is to extract the most of the read and write functionality from M2ulPhyS::restart_files_hdf5 into separate functions. In the previous state the read and write code were intertwined, which made it difficult to understand. The eventual goal is to move the parts that are generic into IODataOrganizer (and related classes), in order to minimize the restart read/write code within M2ulPhyS, which should make it easier to reuse these capabilities for other Solver types down the road. --- src/M2ulPhyS.hpp | 2 + src/io.cpp | 438 ++++++++++++++++++++++++----------------------- 2 files changed, 229 insertions(+), 211 deletions(-) diff --git a/src/M2ulPhyS.hpp b/src/M2ulPhyS.hpp index 05e47db78..56ed34e1c 100644 --- a/src/M2ulPhyS.hpp +++ b/src/M2ulPhyS.hpp @@ -345,6 +345,8 @@ class M2ulPhyS : public TPS::Solver { // i/o routines void restart_files_hdf5(string mode, string inputFileName = std::string()); + void write_restart_files_hdf5(hid_t file, bool serialized_write); + void read_restart_files_hdf5(hid_t file, bool serialized_read); void Check_NAN(); bool Check_JobResubmit(); diff --git a/src/io.cpp b/src/io.cpp index c97a40f7e..3901af758 100644 --- a/src/io.cpp +++ b/src/io.cpp @@ -36,255 +36,199 @@ #include "M2ulPhyS.hpp" #include "utils.hpp" -void M2ulPhyS::restart_files_hdf5(string mode, string inputFileName) { +void M2ulPhyS::write_restart_files_hdf5(hid_t file, bool serialized_write) { MPI_Comm TPSCommWorld = this->groupsMPI->getTPSCommWorld(); -#ifdef HAVE_GRVY - grvy_timer_begin(__func__); -#endif - hid_t file = -1; - hid_t data_soln; - herr_t status; - Vector dataSerial; + // ------------------------------------------------------------------- + // Attributes - relevant solution metadata saved as attributes + // ------------------------------------------------------------------- - string serialName; - if (inputFileName.length() > 0) { - if (inputFileName.substr(inputFileName.length() - 3) != ".h5") { - grvy_printf(gerror, "[ERROR]: M2ulPhyS::restart_files_hdf5 - input file name has a wrong format -> %s\n", - inputFileName.c_str()); - grvy_printf(GRVY_INFO, "format: %s\n", (inputFileName.substr(inputFileName.length() - 3)).c_str()); - exit(ERROR); + // note: all tasks save unless we are writing a serial restart file + if (rank0_ || !serialized_write) { + // current iteration count + h5_save_attribute(file, "iteration", iter); + // total time + h5_save_attribute(file, "time", time); + // timestep + h5_save_attribute(file, "dt", dt); + // solution order + h5_save_attribute(file, "order", order); + // spatial dimension + h5_save_attribute(file, "dimension", dim); + + if (average->ComputeMean()) { + // samples meanUp + h5_save_attribute(file, "samplesMean", average->GetSamplesMean()); + h5_save_attribute(file, "samplesInterval", average->GetSamplesInterval()); } - serialName = inputFileName; - } else { - serialName = "restart_"; - serialName.append(config.GetOutputName()); - serialName.append(".sol.h5"); - } - string fileName; + // code revision +#ifdef BUILD_VERSION + { + hid_t ctype = H5Tcopy(H5T_C_S1); + int shaLength = strlen(BUILD_VERSION); + hsize_t dims[1] = {1}; + H5Tset_size(ctype, shaLength); - assert((mode == "read") || (mode == "write")); + hid_t dspace1dim = H5Screate_simple(1, dims, NULL); - // determine restart file name (either serial or partitioned) + hid_t attr = H5Acreate(file, "revision", ctype, dspace1dim, H5P_DEFAULT, H5P_DEFAULT); + assert(attr >= 0); + herr_t status = H5Awrite(attr, ctype, BUILD_VERSION); + assert(status >= 0); + H5Sclose(dspace1dim); + H5Aclose(attr); + } + } +#endif - if (config.isRestartSerialized(mode)) { - fileName = serialName; - } else { - fileName = groupsMPI->getParallelName(serialName); + // included total dofs for partitioned files + if (!serialized_write) { + int ldofs = vfes->GetNDofs(); + int gdofs; + MPI_Allreduce(&ldofs, &gdofs, 1, MPI_INT, MPI_SUM, TPSCommWorld); + h5_save_attribute(file, "dofs_global", gdofs); } - // Variables used if (and only if) restarting from different order - FiniteElementCollection *aux_fec = NULL; - ParFiniteElementSpace *aux_vfes = NULL; - ParGridFunction *aux_U = NULL; - double *aux_U_data = NULL; - int auxOrder = -1; + // ------------------------------------------------------------------- + // Write solution data defined by IO families + // ------------------------------------------------------------------- + hsize_t dims[1]; - if (rank0_) cout << "HDF5 restart files mode: " << mode << endl; + //------------------------------------------------------- + // Loop over defined IO families to save desired output + //------------------------------------------------------- + for (size_t n = 0; n < ioData.families_.size(); n++) { + IOFamily &fam = ioData.families_[n]; - // open restart files - if (mode == "write") { - if (rank0_ || config.isRestartPartitioned(mode)) { - file = H5Fcreate(fileName.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); - assert(file >= 0); - } - } else if (mode == "read") { - if (config.isRestartSerialized(mode)) { - if (rank0_) { - file = H5Fopen(fileName.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT); - assert(file >= 0); - } + if (serialized_write && (nprocs_ > 1) && (rank0_)) { + assert((locToGlobElem != NULL) && (partitioning_ != NULL)); + assert(fam.serial_fes != NULL); + dims[0] = fam.serial_fes->GetNDofs(); } else { - // verify we have all desired files and open on each process - int gstatus; - int status = static_cast(file_exists(fileName)); - MPI_Allreduce(&status, &gstatus, 1, MPI_INT, MPI_MIN, TPSCommWorld); + dims[0] = fam.pfunc_->ParFESpace()->GetNDofs(); + } + // define groups based on defined IO families + if (rank0_) { + grvy_printf(ginfo, "\nCreating HDF5 group for defined IO families\n"); + grvy_printf(ginfo, "--> %s : %s\n", fam.group_.c_str(), fam.description_.c_str()); + } - if (gstatus == 0) { - grvy_printf(gerror, "[ERROR]: Unable to access desired restart file -> %s\n", fileName.c_str()); - exit(ERROR); - } + hid_t group = -1; + hid_t dataspace = -1; - file = H5Fopen(fileName.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT); - assert(file >= 0); + if (rank0_ || !serialized_write) { + dataspace = H5Screate_simple(1, dims, NULL); + assert(dataspace >= 0); + group = H5Gcreate(file, fam.group_.c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + assert(group >= 0); } - } - - // ------------------------------------------------------------------- - // Attributes - relevant solution metadata saved as attributes - // ------------------------------------------------------------------- - // hid_t aid, attr; - hid_t attr; - if (mode == "write") { - // note: all tasks save unless we are writing a serial restart file - if (rank0_ || config.isRestartPartitioned(mode)) { - // current iteration count - h5_save_attribute(file, "iteration", iter); - // total time - h5_save_attribute(file, "time", time); - // timestep - h5_save_attribute(file, "dt", dt); - // solution order - h5_save_attribute(file, "order", order); - // spatial dimension - h5_save_attribute(file, "dimension", dim); - - if (average->ComputeMean()) { - // samples meanUp - h5_save_attribute(file, "samplesMean", average->GetSamplesMean()); - h5_save_attribute(file, "samplesInterval", average->GetSamplesInterval()); - } - // code revision -#ifdef BUILD_VERSION - { - hid_t ctype = H5Tcopy(H5T_C_S1); - int shaLength = strlen(BUILD_VERSION); - hsize_t dims[1] = {1}; - H5Tset_size(ctype, shaLength); - - hid_t dspace1dim = H5Screate_simple(1, dims, NULL); - - attr = H5Acreate(file, "revision", ctype, dspace1dim, H5P_DEFAULT, H5P_DEFAULT); - assert(attr >= 0); - status = H5Awrite(attr, ctype, BUILD_VERSION); - assert(status >= 0); - H5Sclose(dspace1dim); - H5Aclose(attr); - } + // get pointer to raw data + double *data = fam.pfunc_->HostReadWrite(); + // special case if writing a serial restart + if (serialized_write && (nprocs_ > 1)) { + serialize_soln_for_write(fam, groupsMPI, mesh->GetNE(), nelemGlobal_, locToGlobElem, partitioning_); + if (rank0_) data = fam.serial_sol->HostReadWrite(); } -#endif - // included total dofs for partitioned files - if (!config.isRestartSerialized(mode)) { - int ldofs = vfes->GetNDofs(); - int gdofs; - MPI_Allreduce(&ldofs, &gdofs, 1, MPI_INT, MPI_SUM, TPSCommWorld); - h5_save_attribute(file, "dofs_global", gdofs); + // get defined variables for this IO family + vector vars = ioData.vars_[fam.group_]; + + // save raw data + if (rank0_ || !serialized_write) { + for (auto var : vars) write_soln_data(group, var.varName_, dataspace, data + var.index_ * dims[0], rank0_); } - } else { // read - int read_order; - // normal restarts have each process read their own portion of the - // solution; a serial restart only reads on rank 0 and distributes to - // the remaining processes + if (group >= 0) H5Gclose(group); + if (dataspace >= 0) H5Sclose(dataspace); + } +} + +void M2ulPhyS::read_restart_files_hdf5(hid_t file, bool serialized_read) { + MPI_Comm TPSCommWorld = this->groupsMPI->getTPSCommWorld(); - if (rank0_ || config.isRestartPartitioned(mode)) { - h5_read_attribute(file, "iteration", iter); - h5_read_attribute(file, "time", time); - h5_read_attribute(file, "dt", dt); - h5_read_attribute(file, "order", read_order); - if (average->ComputeMean() && config.GetRestartMean()) { - int samplesMean, intervals; - h5_read_attribute(file, "samplesMean", samplesMean); - h5_read_attribute(file, "samplesInterval", intervals); - average->SetSamplesMean(samplesMean); - average->SetSamplesInterval(intervals); - } + int read_order; + + // Variables used if (and only if) restarting from different order + FiniteElementCollection *aux_fec = NULL; + ParFiniteElementSpace *aux_vfes = NULL; + ParGridFunction *aux_U = NULL; + double *aux_U_data = NULL; + int auxOrder = -1; + + // normal restarts have each process read their own portion of the + // solution; a serial restart only reads on rank 0 and distributes to + // the remaining processes + + if (rank0_ || !serialized_read) { + h5_read_attribute(file, "iteration", iter); + h5_read_attribute(file, "time", time); + h5_read_attribute(file, "dt", dt); + h5_read_attribute(file, "order", read_order); + if (average->ComputeMean() && config.GetRestartMean()) { + int samplesMean, intervals; + h5_read_attribute(file, "samplesMean", samplesMean); + h5_read_attribute(file, "samplesInterval", intervals); + average->SetSamplesMean(samplesMean); + average->SetSamplesInterval(intervals); } + } - if (config.isRestartSerialized(mode)) { - MPI_Bcast(&iter, 1, MPI_INT, 0, TPSCommWorld); - MPI_Bcast(&time, 1, MPI_DOUBLE, 0, TPSCommWorld); - MPI_Bcast(&dt, 1, MPI_DOUBLE, 0, TPSCommWorld); - MPI_Bcast(&read_order, 1, MPI_INT, 0, TPSCommWorld); - if (average->ComputeMean() && config.GetRestartMean()) { - int sampMean = average->GetSamplesMean(); - int intervals = average->GetSamplesInterval(); - MPI_Bcast(&sampMean, 1, MPI_INT, 0, TPSCommWorld); - MPI_Bcast(&intervals, 1, MPI_INT, 0, TPSCommWorld); - } + if (serialized_read) { + MPI_Bcast(&iter, 1, MPI_INT, 0, TPSCommWorld); + MPI_Bcast(&time, 1, MPI_DOUBLE, 0, TPSCommWorld); + MPI_Bcast(&dt, 1, MPI_DOUBLE, 0, TPSCommWorld); + MPI_Bcast(&read_order, 1, MPI_INT, 0, TPSCommWorld); + if (average->ComputeMean() && config.GetRestartMean()) { + int sampMean = average->GetSamplesMean(); + int intervals = average->GetSamplesInterval(); + MPI_Bcast(&sampMean, 1, MPI_INT, 0, TPSCommWorld); + MPI_Bcast(&intervals, 1, MPI_INT, 0, TPSCommWorld); } + } - if (loadFromAuxSol) { - assert(config.RestartSerial() == "no"); // koomie, check cis - auxOrder = read_order; + if (loadFromAuxSol) { + assert(!serialized_read); // koomie, check cis + auxOrder = read_order; - if (basisType == 0) - aux_fec = new DG_FECollection(auxOrder, dim, BasisType::GaussLegendre); - else if (basisType == 1) - aux_fec = new DG_FECollection(auxOrder, dim, BasisType::GaussLobatto); + if (basisType == 0) + aux_fec = new DG_FECollection(auxOrder, dim, BasisType::GaussLegendre); + else if (basisType == 1) + aux_fec = new DG_FECollection(auxOrder, dim, BasisType::GaussLobatto); - aux_vfes = new ParFiniteElementSpace(mesh, aux_fec, num_equation, Ordering::byNODES); + aux_vfes = new ParFiniteElementSpace(mesh, aux_fec, num_equation, Ordering::byNODES); - aux_U_data = new double[num_equation * aux_vfes->GetNDofs()]; - aux_U = new ParGridFunction(aux_vfes, aux_U_data); - } else { - assert(read_order == order); - } + aux_U_data = new double[num_equation * aux_vfes->GetNDofs()]; + aux_U = new ParGridFunction(aux_vfes, aux_U_data); + } else { + assert(read_order == order); + } - if (rank0_) { - grvy_printf(ginfo, "Restarting from iteration = %i\n", iter); - grvy_printf(ginfo, "--> time = %e\n", time); - grvy_printf(ginfo, "--> dt = %e\n", dt); - if (average != NULL) { - grvy_printf(ginfo, "Restarting averages with %i\n samples", average->GetSamplesMean()); - } + if (rank0_) { + grvy_printf(ginfo, "Restarting from iteration = %i\n", iter); + grvy_printf(ginfo, "--> time = %e\n", time); + grvy_printf(ginfo, "--> dt = %e\n", dt); + if (average != NULL) { + grvy_printf(ginfo, "Restarting averages with %i\n samples", average->GetSamplesMean()); } } // ------------------------------------------------------------------- // Read/write solution data defined by IO families // ------------------------------------------------------------------- - - hsize_t dims[1]; - // hsize_t maxdims[1]; hsize_t numInSoln = 0; + hid_t data_soln; //------------------------------------------------------- // Loop over defined IO families to save desired output //------------------------------------------------------- for (size_t n = 0; n < ioData.families_.size(); n++) { IOFamily &fam = ioData.families_[n]; - - if (mode == "write") { - if ((config.isRestartSerialized(mode)) && (nprocs_ > 1) && (rank0_)) { - assert((locToGlobElem != NULL) && (partitioning_ != NULL)); - assert(fam.serial_fes != NULL); - dims[0] = fam.serial_fes->GetNDofs(); - } else { - dims[0] = fam.pfunc_->ParFESpace()->GetNDofs(); - } - // define groups based on defined IO families - if (rank0_) { - grvy_printf(ginfo, "\nCreating HDF5 group for defined IO families\n"); - grvy_printf(ginfo, "--> %s : %s\n", fam.group_.c_str(), fam.description_.c_str()); - } - - hid_t group = -1; - hid_t dataspace = -1; - - if (rank0_ || !config.isRestartSerialized(mode)) { - dataspace = H5Screate_simple(1, dims, NULL); - assert(dataspace >= 0); - group = H5Gcreate(file, fam.group_.c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); - assert(group >= 0); - } - - // get pointer to raw data - double *data = fam.pfunc_->HostReadWrite(); - // special case if writing a serial restart - if ((config.isRestartSerialized(mode)) && (nprocs_ > 1)) { - serialize_soln_for_write(fam, groupsMPI, mesh->GetNE(), nelemGlobal_, locToGlobElem, partitioning_); - if (rank0_) data = fam.serial_sol->HostReadWrite(); - } - - // get defined variables for this IO family - vector vars = ioData.vars_[fam.group_]; - - // save raw data - if (rank0_ || config.isRestartPartitioned(mode)) { - for (auto var : vars) write_soln_data(group, var.varName_, dataspace, data + var.index_ * dims[0], rank0_); - } - - if (group >= 0) H5Gclose(group); - if (dataspace >= 0) H5Sclose(dataspace); - - } else if (fam.inRestartFile) { // read mode + if (fam.inRestartFile) { // read mode if (rank0_) cout << "Reading in solutiond data from restart..." << endl; // verify Dofs match expectations with current mesh - if (rank0_ || config.isRestartPartitioned(mode)) { + if (rank0_ || !serialized_read) { hid_t dataspace; vector vars = ioData.vars_[fam.group_]; @@ -300,7 +244,7 @@ void M2ulPhyS::restart_files_hdf5(string mode, string inputFileName) { } int dof = fam.pfunc_->ParFESpace()->GetNDofs(); - if (config.isRestartSerialized(mode)) { + if (serialized_read) { int dof_global; MPI_Reduce(&dof, &dof_global, 1, MPI_INT, MPI_SUM, 0, TPSCommWorld); dof = dof_global; @@ -308,7 +252,7 @@ void M2ulPhyS::restart_files_hdf5(string mode, string inputFileName) { if (loadFromAuxSol && fam.allowsAuxRestart) dof = aux_vfes->GetNDofs(); - if (rank0_ || config.isRestartPartitioned(mode)) assert((int)numInSoln == dof); + if (rank0_ || !serialized_read) assert((int)numInSoln == dof); // get pointer to raw data double *data = fam.pfunc_->HostReadWrite(); @@ -326,11 +270,11 @@ void M2ulPhyS::restart_files_hdf5(string mode, string inputFileName) { if (var.inRestartFile_) { std::string h5Path = fam.group_ + "/" + var.varName_; if (rank0_) grvy_printf(ginfo, "--> Reading h5 path = %s\n", h5Path.c_str()); - if (config.isRestartPartitioned(mode)) { + if (!serialized_read) { read_partitioned_soln_data(file, h5Path.c_str(), var.index_ * numInSoln, data); } else { assert(partitioning_ != NULL); - assert(config.isRestartSerialized("read")); + assert(serialized_read); read_serialized_soln_data(file, h5Path.c_str(), dof, var.index_, data, fam, groupsMPI, partitioning_, nelemGlobal_); } @@ -339,9 +283,7 @@ void M2ulPhyS::restart_files_hdf5(string mode, string inputFileName) { } } - if (file >= 0) H5Fclose(file); - - if (mode == "read" && loadFromAuxSol) { + if (loadFromAuxSol) { if (rank0_) cout << "Interpolating from auxOrder = " << auxOrder << " to order = " << order << endl; // Interpolate from aux to new order @@ -378,6 +320,80 @@ void M2ulPhyS::restart_files_hdf5(string mode, string inputFileName) { delete aux_vfes; delete aux_fec; } +} + +void M2ulPhyS::restart_files_hdf5(string mode, string inputFileName) { + assert((mode == "read") || (mode == "write")); + MPI_Comm TPSCommWorld = this->groupsMPI->getTPSCommWorld(); +#ifdef HAVE_GRVY + grvy_timer_begin(__func__); +#endif + + string serialName; + if (inputFileName.length() > 0) { + if (inputFileName.substr(inputFileName.length() - 3) != ".h5") { + grvy_printf(gerror, "[ERROR]: M2ulPhyS::restart_files_hdf5 - input file name has a wrong format -> %s\n", + inputFileName.c_str()); + grvy_printf(GRVY_INFO, "format: %s\n", (inputFileName.substr(inputFileName.length() - 3)).c_str()); + exit(ERROR); + } + serialName = inputFileName; + } else { + serialName = "restart_"; + serialName.append(config.GetOutputName()); + serialName.append(".sol.h5"); + } + + // determine restart file name (either serial or partitioned) + string fileName; + if (config.isRestartSerialized(mode)) { + fileName = serialName; + } else { + fileName = groupsMPI->getParallelName(serialName); + } + + if (rank0_) cout << "HDF5 restart files mode: " << mode << endl; + + // open restart files + hid_t file = -1; + if (mode == "write") { + if (rank0_ || config.isRestartPartitioned(mode)) { + file = H5Fcreate(fileName.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); + assert(file >= 0); + } + } else if (mode == "read") { + if (config.isRestartSerialized(mode)) { + if (rank0_) { + file = H5Fopen(fileName.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT); + assert(file >= 0); + } + } else { + // verify we have all desired files and open on each process + int gstatus; + int status = static_cast(file_exists(fileName)); + MPI_Allreduce(&status, &gstatus, 1, MPI_INT, MPI_MIN, TPSCommWorld); + + if (gstatus == 0) { + grvy_printf(gerror, "[ERROR]: Unable to access desired restart file -> %s\n", fileName.c_str()); + exit(ERROR); + } + + file = H5Fopen(fileName.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT); + assert(file >= 0); + } + } + + // ------------------------------------------------------------------- + // Attributes - relevant solution metadata saved as attributes + // ------------------------------------------------------------------- + + if (mode == "write") { + write_restart_files_hdf5(file, config.isRestartSerialized(mode)); + } else { // read + read_restart_files_hdf5(file, config.isRestartSerialized(mode)); + } + + if (file >= 0) H5Fclose(file); #ifdef HAVE_GRVY grvy_timer_end(__func__); From 050a14630f7a3295fbe8c02cda3be2e6d4f1da16 Mon Sep 17 00:00:00 2001 From: "Todd A. Oliver" Date: Tue, 30 Jan 2024 11:31:20 -0600 Subject: [PATCH 05/14] Move primary data write support into IODataOrganizer class (#244) This involves 3 functions that now perform the bulk of the work. First, the (previously standalone) function serialize_soln_for_write is moved into the IOFamily class, as serializeForWrite. Second, two new functions (IODataOrganizer::write and IODataOrganizer::writeSerial) are introduced to handle writing field data that has been registered through the IODataOrganizer class. At the moment these functions are separate b/c they require different input data (i.e., writeSerial needs information, such as the partitioning array, that write does not). They could be unified, potentially by extending IODataOrganizer::initializeSerial to obtain this extra information rather than passing it directly to writeSerial. But, this refactor is left for later. --- src/io.cpp | 274 ++++++++++++++++++++++++++++++++--------------------- src/io.hpp | 16 +++- 2 files changed, 176 insertions(+), 114 deletions(-) diff --git a/src/io.cpp b/src/io.cpp index 3901af758..606b16929 100644 --- a/src/io.cpp +++ b/src/io.cpp @@ -89,58 +89,10 @@ void M2ulPhyS::write_restart_files_hdf5(hid_t file, bool serialized_write) { h5_save_attribute(file, "dofs_global", gdofs); } - // ------------------------------------------------------------------- - // Write solution data defined by IO families - // ------------------------------------------------------------------- - hsize_t dims[1]; - - //------------------------------------------------------- - // Loop over defined IO families to save desired output - //------------------------------------------------------- - for (size_t n = 0; n < ioData.families_.size(); n++) { - IOFamily &fam = ioData.families_[n]; - - if (serialized_write && (nprocs_ > 1) && (rank0_)) { - assert((locToGlobElem != NULL) && (partitioning_ != NULL)); - assert(fam.serial_fes != NULL); - dims[0] = fam.serial_fes->GetNDofs(); - } else { - dims[0] = fam.pfunc_->ParFESpace()->GetNDofs(); - } - // define groups based on defined IO families - if (rank0_) { - grvy_printf(ginfo, "\nCreating HDF5 group for defined IO families\n"); - grvy_printf(ginfo, "--> %s : %s\n", fam.group_.c_str(), fam.description_.c_str()); - } - - hid_t group = -1; - hid_t dataspace = -1; - - if (rank0_ || !serialized_write) { - dataspace = H5Screate_simple(1, dims, NULL); - assert(dataspace >= 0); - group = H5Gcreate(file, fam.group_.c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); - assert(group >= 0); - } - - // get pointer to raw data - double *data = fam.pfunc_->HostReadWrite(); - // special case if writing a serial restart - if (serialized_write && (nprocs_ > 1)) { - serialize_soln_for_write(fam, groupsMPI, mesh->GetNE(), nelemGlobal_, locToGlobElem, partitioning_); - if (rank0_) data = fam.serial_sol->HostReadWrite(); - } - - // get defined variables for this IO family - vector vars = ioData.vars_[fam.group_]; - - // save raw data - if (rank0_ || !serialized_write) { - for (auto var : vars) write_soln_data(group, var.varName_, dataspace, data + var.index_ * dims[0], rank0_); - } - - if (group >= 0) H5Gclose(group); - if (dataspace >= 0) H5Sclose(dataspace); + if (serialized_write) { + ioData.writeSerial(file, groupsMPI->getTPSCommWorld(), mesh->GetNE(), nelemGlobal_, locToGlobElem, partitioning_); + } else { + ioData.write(file, rank0_); } } @@ -514,63 +466,6 @@ void partitioning_file_hdf5(std::string mode, const RunConfiguration &config, MP grvy_timer_end(__func__); } -void serialize_soln_for_write(IOFamily &fam, MPI_Groups *groupsMPI, int local_ne, int global_ne, - const int *locToGlobElem, const Array &partitioning) { - MPI_Comm TPSCommWorld = groupsMPI->getTPSCommWorld(); - const bool rank0 = groupsMPI->isWorldRoot(); - - // Ensure consistency - int global_ne_check; - MPI_Reduce(&local_ne, &global_ne_check, 1, MPI_INT, MPI_SUM, 0, TPSCommWorld); - if (rank0) { - assert(global_ne_check == global_ne); - assert(partitioning.Size() == global_ne); - } - assert(locToGlobElem != NULL); - - ParGridFunction *pfunc = fam.pfunc_; - if (rank0) { - grvy_printf(ginfo, "Generating serialized restart file (group %s...)\n", fam.group_.c_str()); - // copy my own data - Array lvdofs, gvdofs; - Vector lsoln; - for (int elem = 0; elem < local_ne; elem++) { - int gelem = locToGlobElem[elem]; - fam.pfunc_->ParFESpace()->GetElementVDofs(elem, lvdofs); - pfunc->GetSubVector(lvdofs, lsoln); - fam.serial_fes->GetElementVDofs(gelem, gvdofs); - fam.serial_sol->SetSubVector(gvdofs, lsoln); - } - - // have rank 0 receive data from other tasks and copy its own - for (int gelem = 0; gelem < global_ne; gelem++) { - int from_rank = partitioning[gelem]; - if (from_rank != 0) { - fam.serial_fes->GetElementVDofs(gelem, gvdofs); - lsoln.SetSize(gvdofs.Size()); - - MPI_Recv(lsoln.HostReadWrite(), gvdofs.Size(), MPI_DOUBLE, from_rank, gelem, TPSCommWorld, MPI_STATUS_IGNORE); - - fam.serial_sol->SetSubVector(gvdofs, lsoln); - } - } - - } else { - // have non-zero ranks send their data to rank 0 - Array lvdofs; - Vector lsoln; - for (int elem = 0; elem < local_ne; elem++) { - int gelem = locToGlobElem[elem]; - // assert(gelem > 0); - fam.pfunc_->ParFESpace()->GetElementVDofs(elem, lvdofs); - pfunc->GetSubVector(lvdofs, lsoln); // work for gpu build? - - // send to task 0 - MPI_Send(lsoln.HostReadWrite(), lsoln.Size(), MPI_DOUBLE, 0, gelem, TPSCommWorld); - } - } -} // end function: serialize_soln_for_write() - // convenience function to read solution data for parallel restarts void read_partitioned_soln_data(hid_t file, string varName, size_t index, double *data) { hid_t data_soln; @@ -728,10 +623,68 @@ void M2ulPhyS::readTable(const std::string &inputPath, TableInput &result) { // Routines for I/O data organizer helper class // --------------------------------------------- +void IOFamily::serializeForWrite(MPI_Comm comm, int local_ne, int global_ne, const int *locToGlobElem, + const Array &partitioning) { + int rank; + MPI_Comm_rank(comm, &rank); + const bool rank0 = (rank == 0); + + // Ensure consistency + int global_ne_check; + MPI_Reduce(&local_ne, &global_ne_check, 1, MPI_INT, MPI_SUM, 0, comm); + if (rank0) { + assert(global_ne_check == global_ne); + assert(partitioning.Size() == global_ne); + } + assert(locToGlobElem != NULL); + + ParGridFunction *pfunc = this->pfunc_; + if (rank0) { + grvy_printf(ginfo, "Generating serialized restart file (group %s...)\n", this->group_.c_str()); + // copy my own data + Array lvdofs, gvdofs; + Vector lsoln; + for (int elem = 0; elem < local_ne; elem++) { + int gelem = locToGlobElem[elem]; + this->pfunc_->ParFESpace()->GetElementVDofs(elem, lvdofs); + pfunc->GetSubVector(lvdofs, lsoln); + this->serial_fes->GetElementVDofs(gelem, gvdofs); + this->serial_sol->SetSubVector(gvdofs, lsoln); + } + + // have rank 0 receive data from other tasks and copy its own + for (int gelem = 0; gelem < global_ne; gelem++) { + int from_rank = partitioning[gelem]; + if (from_rank != 0) { + this->serial_fes->GetElementVDofs(gelem, gvdofs); + lsoln.SetSize(gvdofs.Size()); + + MPI_Recv(lsoln.HostReadWrite(), gvdofs.Size(), MPI_DOUBLE, from_rank, gelem, comm, MPI_STATUS_IGNORE); + + this->serial_sol->SetSubVector(gvdofs, lsoln); + } + } + + } else { + // have non-zero ranks send their data to rank 0 + Array lvdofs; + Vector lsoln; + for (int elem = 0; elem < local_ne; elem++) { + int gelem = locToGlobElem[elem]; + // assert(gelem > 0); + this->pfunc_->ParFESpace()->GetElementVDofs(elem, lvdofs); + pfunc->GetSubVector(lvdofs, lsoln); // work for gpu build? + + // send to task 0 + MPI_Send(lsoln.HostReadWrite(), lsoln.Size(), MPI_DOUBLE, 0, gelem, comm); + } + } +} + // register a new IO family which maps to a ParGridFunction void IODataOrganizer::registerIOFamily(std::string description, std::string group, ParGridFunction *pfunc, bool auxRestart, bool _inRestartFile) { - IOFamily family{description, group, pfunc}; + IOFamily family(description, group, pfunc); std::vector vars; family.allowsAuxRestart = auxRestart; family.inRestartFile = _inRestartFile; @@ -775,6 +728,107 @@ void IODataOrganizer::initializeSerial(bool root, bool serial, Mesh *serial_mesh } } +void IODataOrganizer::write(hid_t file, bool rank0) { + // ------------------------------------------------------------------- + // Write solution data defined by IO families + // ------------------------------------------------------------------- + hsize_t dims[1]; + + //------------------------------------------------------- + // Loop over defined IO families to save desired output + //------------------------------------------------------- + for (size_t n = 0; n < families_.size(); n++) { + IOFamily &fam = families_[n]; + + dims[0] = fam.pfunc_->ParFESpace()->GetNDofs(); + + // define groups based on defined IO families + if (rank0) { + grvy_printf(ginfo, "\nCreating HDF5 group for defined IO families\n"); + grvy_printf(ginfo, "--> %s : %s\n", fam.group_.c_str(), fam.description_.c_str()); + } + + hid_t group = -1; + hid_t dataspace = -1; + + dataspace = H5Screate_simple(1, dims, NULL); + assert(dataspace >= 0); + group = H5Gcreate(file, fam.group_.c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + assert(group >= 0); + + // get pointer to raw data + double *data = fam.pfunc_->HostReadWrite(); + + // get defined variables for this IO family + vector vars = vars_[fam.group_]; + + // save raw data + for (auto var : vars) write_soln_data(group, var.varName_, dataspace, data + var.index_ * dims[0], rank0); + + if (group >= 0) H5Gclose(group); + if (dataspace >= 0) H5Sclose(dataspace); + } +} + +void IODataOrganizer::writeSerial(hid_t file, MPI_Comm comm, int local_ne, int global_ne, const int *locToGlobElem, + const Array &partitioning) { + int rank; + MPI_Comm_rank(comm, &rank); + int nprocs; + MPI_Comm_size(comm, &nprocs); + const bool rank0 = (rank == 0); + + // if running serially, writeSerial and write are logically the + // same, but infrastructure is different. Thus, for serial run, + // just call write() here + if (nprocs == 1) { + this->write(file, rank0); + return; + } + + hsize_t dims[1]; + + // Loop over defined IO families to save desired output + for (size_t n = 0; n < families_.size(); n++) { + IOFamily &fam = families_[n]; + + if (rank0) { + assert((locToGlobElem != NULL) && (partitioning.Size() == global_ne)); + assert(fam.serial_fes != NULL); + dims[0] = fam.serial_fes->GetNDofs(); + grvy_printf(ginfo, "\nCreating HDF5 group for defined IO families\n"); + grvy_printf(ginfo, "--> %s : %s\n", fam.group_.c_str(), fam.description_.c_str()); + } + + hid_t group = -1; + hid_t dataspace = -1; + + if (rank0) { + dataspace = H5Screate_simple(1, dims, NULL); + assert(dataspace >= 0); + group = H5Gcreate(file, fam.group_.c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + assert(group >= 0); + } + + fam.serializeForWrite(comm, local_ne, global_ne, locToGlobElem, partitioning); + + // get pointer to raw data + double *data = nullptr; + if (rank0) data = fam.serial_sol->HostReadWrite(); + + // get defined variables for this IO family + vector vars = vars_[fam.group_]; + + // save raw data + if (rank0) { + for (auto var : vars) write_soln_data(group, var.varName_, dataspace, data + var.index_ * dims[0], rank0); + } + + if (group >= 0) H5Gclose(group); + if (dataspace >= 0) H5Sclose(dataspace); + } +} + IODataOrganizer::~IODataOrganizer() { for (size_t n = 0; n < families_.size(); n++) { IOFamily fam = families_[n]; diff --git a/src/io.hpp b/src/io.hpp index 0e7637f12..faa778147 100644 --- a/src/io.hpp +++ b/src/io.hpp @@ -52,8 +52,14 @@ class IOFamily { // true if the family is to be found in restart files bool inRestartFile; - mfem::FiniteElementSpace *serial_fes; - mfem::GridFunction *serial_sol; + mfem::FiniteElementSpace *serial_fes = nullptr; + mfem::GridFunction *serial_sol = nullptr; + + IOFamily(std::string desc, std::string grp, mfem::ParGridFunction *pf) + : description_(desc), group_(grp), pfunc_(pf) {} + + void serializeForWrite(MPI_Comm comm, int local_ne, int global_ne, const int *locToGlobElem, + const Array &partitioning); }; class IOVar { @@ -76,6 +82,10 @@ class IODataOrganizer { int getIOFamilyIndex(std::string group); void initializeSerial(bool root, bool serial, mfem::Mesh *serial_mesh); + + void write(hid_t file, bool rank0); + void writeSerial(hid_t file, MPI_Comm comm, int local_ne, int global_ne, const int *locToGlobElem, + const Array &partitioning); }; void read_partitioned_soln_data(hid_t file, std::string varName, size_t index, double *data); @@ -84,6 +94,4 @@ void read_serialized_soln_data(hid_t file, std::string varName, int numDof, int void write_soln_data(hid_t group, std::string varName, hid_t dataspace, double *data, bool rank0); void partitioning_file_hdf5(std::string mode, const RunConfiguration &config, MPI_Groups *groupsMPI, int nelemGlobal, mfem::Array &partitioning); -void serialize_soln_for_write(IOFamily &fam, MPI_Groups *groupsMPI, int local_ne, int global_ne, - const int *locToGlobElem, const mfem::Array &partitioning); #endif // IO_HPP_ From bd58df19e1552704ce2da40f331f837b021cf249 Mon Sep 17 00:00:00 2001 From: "Todd A. Oliver" Date: Tue, 30 Jan 2024 13:47:33 -0600 Subject: [PATCH 06/14] Move primary data read support into IODataOrganizer class (#244) This is mostly analogous to the write support from the previous commit, with new read functions being added to IODataOrganizer. But, in addition we have to support the option to read data from a different order solution (to support restarting from lower order fields). This is facilitated by adding "aux" objects (FiniteElementCollection FiniteElementSpace, and ParGridFunction) to the IOFamily class. In the previous implementation, these were instantiated directly within the read function and were specific to the DG discretization used by M2ulPhyS. The new implementation should be generic (although that generality is not yet tested). --- src/M2ulPhyS.cpp | 2 +- src/io.cpp | 316 ++++++++++++++++++++++++++++++----------------- src/io.hpp | 13 +- 3 files changed, 218 insertions(+), 113 deletions(-) diff --git a/src/M2ulPhyS.cpp b/src/M2ulPhyS.cpp index de7be70ec..ac820a276 100644 --- a/src/M2ulPhyS.cpp +++ b/src/M2ulPhyS.cpp @@ -1795,7 +1795,7 @@ void M2ulPhyS::initSolutionAndVisualizationVectors() { } // define solution parameters for i/o - ioData.registerIOFamily("Solution state variables", "/solution", U); + ioData.registerIOFamily("Solution state variables", "/solution", U, true, true, fec); ioData.registerIOVar("/solution", "density", 0); ioData.registerIOVar("/solution", "rho-u", 1); ioData.registerIOVar("/solution", "rho-v", 2); diff --git a/src/io.cpp b/src/io.cpp index 606b16929..b217d18c0 100644 --- a/src/io.cpp +++ b/src/io.cpp @@ -89,6 +89,9 @@ void M2ulPhyS::write_restart_files_hdf5(hid_t file, bool serialized_write) { h5_save_attribute(file, "dofs_global", gdofs); } + // ------------------------------------------------------------------- + // Data - actual data write handled by IODataOrganizer class + // ------------------------------------------------------------------- if (serialized_write) { ioData.writeSerial(file, groupsMPI->getTPSCommWorld(), mesh->GetNE(), nelemGlobal_, locToGlobElem, partitioning_); } else { @@ -98,20 +101,16 @@ void M2ulPhyS::write_restart_files_hdf5(hid_t file, bool serialized_write) { void M2ulPhyS::read_restart_files_hdf5(hid_t file, bool serialized_read) { MPI_Comm TPSCommWorld = this->groupsMPI->getTPSCommWorld(); - int read_order; - // Variables used if (and only if) restarting from different order - FiniteElementCollection *aux_fec = NULL; - ParFiniteElementSpace *aux_vfes = NULL; - ParGridFunction *aux_U = NULL; - double *aux_U_data = NULL; - int auxOrder = -1; - // normal restarts have each process read their own portion of the // solution; a serial restart only reads on rank 0 and distributes to // the remaining processes + // ------------------------------------------------------------------- + // Attributes - read relevant solution metadata for this Solver + // ------------------------------------------------------------------- + if (rank0_ || !serialized_read) { h5_read_attribute(file, "iteration", iter); h5_read_attribute(file, "time", time); @@ -139,23 +138,6 @@ void M2ulPhyS::read_restart_files_hdf5(hid_t file, bool serialized_read) { } } - if (loadFromAuxSol) { - assert(!serialized_read); // koomie, check cis - auxOrder = read_order; - - if (basisType == 0) - aux_fec = new DG_FECollection(auxOrder, dim, BasisType::GaussLegendre); - else if (basisType == 1) - aux_fec = new DG_FECollection(auxOrder, dim, BasisType::GaussLobatto); - - aux_vfes = new ParFiniteElementSpace(mesh, aux_fec, num_equation, Ordering::byNODES); - - aux_U_data = new double[num_equation * aux_vfes->GetNDofs()]; - aux_U = new ParGridFunction(aux_vfes, aux_U_data); - } else { - assert(read_order == order); - } - if (rank0_) { grvy_printf(ginfo, "Restarting from iteration = %i\n", iter); grvy_printf(ginfo, "--> time = %e\n", time); @@ -166,89 +148,22 @@ void M2ulPhyS::read_restart_files_hdf5(hid_t file, bool serialized_read) { } // ------------------------------------------------------------------- - // Read/write solution data defined by IO families + // Data - actual data read handled by IODataOrganizer class // ------------------------------------------------------------------- - hsize_t numInSoln = 0; - hid_t data_soln; - - //------------------------------------------------------- - // Loop over defined IO families to save desired output - //------------------------------------------------------- - for (size_t n = 0; n < ioData.families_.size(); n++) { - IOFamily &fam = ioData.families_[n]; - if (fam.inRestartFile) { // read mode - if (rank0_) cout << "Reading in solutiond data from restart..." << endl; - - // verify Dofs match expectations with current mesh - if (rank0_ || !serialized_read) { - hid_t dataspace; - - vector vars = ioData.vars_[fam.group_]; - string varGroupName = fam.group_; - varGroupName.append("/"); - varGroupName.append(vars[0].varName_); - - data_soln = H5Dopen2(file, varGroupName.c_str(), H5P_DEFAULT); - assert(data_soln >= 0); - dataspace = H5Dget_space(data_soln); - numInSoln = H5Sget_simple_extent_npoints(dataspace); - H5Dclose(data_soln); - } - - int dof = fam.pfunc_->ParFESpace()->GetNDofs(); - if (serialized_read) { - int dof_global; - MPI_Reduce(&dof, &dof_global, 1, MPI_INT, MPI_SUM, 0, TPSCommWorld); - dof = dof_global; - } - - if (loadFromAuxSol && fam.allowsAuxRestart) dof = aux_vfes->GetNDofs(); - - if (rank0_ || !serialized_read) assert((int)numInSoln == dof); - - // get pointer to raw data - double *data = fam.pfunc_->HostReadWrite(); - // special case if starting from aux soln - if (loadFromAuxSol && fam.allowsAuxRestart) { - data = aux_U->HostReadWrite(); - - // ks note: would need to add additional logic to handle read of - // multiple pargridfunctions when changing the solution order - assert(ioData.families_.size() == 1); - } - - vector vars = ioData.vars_[fam.group_]; - for (auto var : vars) { - if (var.inRestartFile_) { - std::string h5Path = fam.group_ + "/" + var.varName_; - if (rank0_) grvy_printf(ginfo, "--> Reading h5 path = %s\n", h5Path.c_str()); - if (!serialized_read) { - read_partitioned_soln_data(file, h5Path.c_str(), var.index_ * numInSoln, data); - } else { - assert(partitioning_ != NULL); - assert(serialized_read); - read_serialized_soln_data(file, h5Path.c_str(), dof, var.index_, data, fam, groupsMPI, partitioning_, - nelemGlobal_); - } - } - } - } + if (!loadFromAuxSol && !serialized_read) { + assert(read_order == order); + ioData.read(file, rank0_); + } else if (!loadFromAuxSol && serialized_read) { + assert(read_order == order); + ioData.readSerial(file, rank0_, groupsMPI, partitioning_, nelemGlobal_); + } else if (loadFromAuxSol && !serialized_read) { + ioData.readChangeOrder(file, rank0_, read_order); + } else { + if (rank0_) grvy_printf(gerror, "[ERROR]: Unsupported read mode requested.\n"); + exit(ERROR); } if (loadFromAuxSol) { - if (rank0_) cout << "Interpolating from auxOrder = " << auxOrder << " to order = " << order << endl; - - // Interpolate from aux to new order - U->ProjectGridFunction(*aux_U); - - // If interpolation was successful, the L2 norm of the - // difference between the auxOrder and order versions should be - // zero (well... to within round-off-induced error) - VectorGridFunctionCoefficient lowOrderCoeff(aux_U); - double err = U->ComputeL2Error(lowOrderCoeff); - - if (rank0_) cout << "|| interpolation error ||_2 = " << err << endl; - // Update primitive variables. This will be done automatically // before taking a time step, but we may write paraview output // prior to that, in which case the primitives are incorrect. @@ -265,12 +180,6 @@ void M2ulPhyS::read_restart_files_hdf5(hid_t file, bool serialized_read) { mixture->GetPrimitivesFromConservatives(conserved, primitive); for (int eq = 0; eq < num_equation; eq++) dataUp[i + eq * vfes->GetNDofs()] = primitive[eq]; } - - // clean up aux data - delete aux_U; - delete aux_U_data; - delete aux_vfes; - delete aux_fec; } } @@ -683,11 +592,15 @@ void IOFamily::serializeForWrite(MPI_Comm comm, int local_ne, int global_ne, con // register a new IO family which maps to a ParGridFunction void IODataOrganizer::registerIOFamily(std::string description, std::string group, ParGridFunction *pfunc, - bool auxRestart, bool _inRestartFile) { + bool auxRestart, bool _inRestartFile, FiniteElementCollection *fec) { IOFamily family(description, group, pfunc); std::vector vars; family.allowsAuxRestart = auxRestart; family.inRestartFile = _inRestartFile; + family.fec_ = fec; + if (family.allowsAuxRestart) { + assert(family.fec_ != NULL); + } families_.push_back(family); vars_[group] = vars; @@ -829,6 +742,187 @@ void IODataOrganizer::writeSerial(hid_t file, MPI_Comm comm, int local_ne, int g } } +void IODataOrganizer::read(hid_t file, bool rank0) { + hsize_t numInSoln = 0; + hid_t data_soln; + + // Loop over defined IO families to save desired output + for (auto fam : families_) { + if (fam.inRestartFile) { // read mode + if (rank0) cout << "Reading in solutiond data from restart..." << endl; + + // verify Dofs match expectations with current mesh + hid_t dataspace; + + vector vars = vars_[fam.group_]; + string varGroupName = fam.group_; + varGroupName.append("/"); + varGroupName.append(vars[0].varName_); + + data_soln = H5Dopen2(file, varGroupName.c_str(), H5P_DEFAULT); + assert(data_soln >= 0); + dataspace = H5Dget_space(data_soln); + numInSoln = H5Sget_simple_extent_npoints(dataspace); + H5Dclose(data_soln); + + int dof = fam.pfunc_->ParFESpace()->GetNDofs(); + assert((int)numInSoln == dof); + + // get pointer to raw data + double *data = fam.pfunc_->HostReadWrite(); + + // read from file into appropriate spot in data + for (auto var : vars) { + if (var.inRestartFile_) { + std::string h5Path = fam.group_ + "/" + var.varName_; + if (rank0) grvy_printf(ginfo, "--> Reading h5 path = %s\n", h5Path.c_str()); + read_partitioned_soln_data(file, h5Path.c_str(), var.index_ * numInSoln, data); + } + } + } + } +} + +void IODataOrganizer::readSerial(hid_t file, bool rank0, MPI_Groups *groupsMPI, Array partitioning, + int global_ne) { + MPI_Comm TPSCommWorld = groupsMPI->getTPSCommWorld(); + // ------------------------------------------------------------------- + // Read/write solution data defined by IO families + // ------------------------------------------------------------------- + hsize_t numInSoln = 0; + hid_t data_soln; + + //------------------------------------------------------- + // Loop over defined IO families to save desired output + //------------------------------------------------------- + for (size_t n = 0; n < families_.size(); n++) { + IOFamily &fam = families_[n]; + if (fam.inRestartFile) { // read mode + if (rank0) cout << "Reading in solutiond data from restart..." << endl; + + // verify Dofs match expectations with current mesh + if (rank0) { + hid_t dataspace; + + vector vars = vars_[fam.group_]; + string varGroupName = fam.group_; + varGroupName.append("/"); + varGroupName.append(vars[0].varName_); + + data_soln = H5Dopen2(file, varGroupName.c_str(), H5P_DEFAULT); + assert(data_soln >= 0); + dataspace = H5Dget_space(data_soln); + numInSoln = H5Sget_simple_extent_npoints(dataspace); + H5Dclose(data_soln); + } + + int dof = fam.pfunc_->ParFESpace()->GetNDofs(); + + int dof_global; + MPI_Reduce(&dof, &dof_global, 1, MPI_INT, MPI_SUM, 0, TPSCommWorld); + dof = dof_global; + + if (rank0) assert((int)numInSoln == dof); + + // get pointer to raw data + double *data = fam.pfunc_->HostReadWrite(); + + vector vars = vars_[fam.group_]; + for (auto var : vars) { + if (var.inRestartFile_) { + std::string h5Path = fam.group_ + "/" + var.varName_; + if (rank0) grvy_printf(ginfo, "--> Reading h5 path = %s\n", h5Path.c_str()); + assert(partitioning.Size() == global_ne); + read_serialized_soln_data(file, h5Path.c_str(), dof, var.index_, data, fam, groupsMPI, partitioning, + global_ne); + } + } + } + } // end loop over families +} + +void IODataOrganizer::readChangeOrder(hid_t file, bool rank0, int read_order) { + hsize_t numInSoln = 0; + hid_t data_soln; + double *aux_U_data = NULL; + + // NOTE(trevilo): this assert is left over from the original + // implementation. I believe the refactored implementation can + // support multiple IOFamily objects and changing order on restart. + // But, since we do not have a use or test case, I am leaving this + // assert for the time being. + assert(families_.size() == 1); + + // Loop over defined IO families to save desired output + for (size_t n = 0; n < families_.size(); n++) { + IOFamily &fam = families_[n]; + if (fam.inRestartFile && fam.allowsAuxRestart) { // read mode + if (rank0) cout << "Reading in solutiond data from restart..." << endl; + + // verify Dofs match expectations with current mesh + hid_t dataspace; + + vector vars = vars_[fam.group_]; + string varGroupName = fam.group_; + varGroupName.append("/"); + varGroupName.append(vars[0].varName_); + + data_soln = H5Dopen2(file, varGroupName.c_str(), H5P_DEFAULT); + assert(data_soln >= 0); + dataspace = H5Dget_space(data_soln); + numInSoln = H5Sget_simple_extent_npoints(dataspace); + H5Dclose(data_soln); + + assert(fam.fec_ != NULL); + fam.aux_fec_ = fam.fec_->Clone(read_order); + + int nvars = vars_[fam.group_].size(); + // fam.aux_pfes_ = new ParFiniteElementSpace(mesh, fam.aux_fec_, nvars, Ordering::byNODES); + fam.aux_pfes_ = + new ParFiniteElementSpace(fam.pfunc_->ParFESpace()->GetParMesh(), fam.aux_fec_, nvars, Ordering::byNODES); + + aux_U_data = new double[nvars * fam.aux_pfes_->GetNDofs()]; + fam.aux_pfunc_ = new ParGridFunction(fam.aux_pfes_, aux_U_data); + + int dof = fam.aux_pfes_->GetNDofs(); + + assert((int)numInSoln == dof); + + // get pointer to raw data + double *data = fam.aux_pfunc_->HostReadWrite(); + + for (auto var : vars) { + if (var.inRestartFile_) { + std::string h5Path = fam.group_ + "/" + var.varName_; + if (rank0) grvy_printf(ginfo, "--> Reading h5 path = %s\n", h5Path.c_str()); + read_partitioned_soln_data(file, h5Path.c_str(), var.index_ * numInSoln, data); + } + } + + if (rank0) cout << "Interpolating from read_order = " << read_order << " to solution space" << endl; + + // Interpolate from aux to new order + fam.pfunc_->ProjectGridFunction(*fam.aux_pfunc_); + + // If interpolation was successful, the L2 norm of the + // difference between the auxOrder and order versions should be + // zero (well... to within round-off-induced error) + VectorGridFunctionCoefficient lowOrderCoeff(fam.aux_pfunc_); + double err = fam.pfunc_->ComputeL2Error(lowOrderCoeff); + + if (rank0) cout << "|| interpolation error ||_2 = " << err << endl; + + // clean up aux data + delete fam.aux_pfunc_; + delete[] aux_U_data; + delete fam.aux_pfes_; + delete fam.aux_fec_; + } else { + if (rank0) cout << "Skipping family " << fam.group_ << " in readChangeOrder" << endl; + } + } // end loop over families +} + IODataOrganizer::~IODataOrganizer() { for (size_t n = 0; n < families_.size(); n++) { IOFamily fam = families_[n]; diff --git a/src/io.hpp b/src/io.hpp index faa778147..38b95ba5a 100644 --- a/src/io.hpp +++ b/src/io.hpp @@ -52,9 +52,16 @@ class IOFamily { // true if the family is to be found in restart files bool inRestartFile; + // space and grid function used for serial read/write mfem::FiniteElementSpace *serial_fes = nullptr; mfem::GridFunction *serial_sol = nullptr; + // Variables used for "auxilliary" read (i.e., restart from different order soln) + mfem::FiniteElementCollection *fec_ = nullptr; // collection used to instantiate pfunc_ + mfem::FiniteElementCollection *aux_fec_ = nullptr; + mfem::ParFiniteElementSpace *aux_pfes_ = nullptr; + mfem::ParGridFunction *aux_pfunc_ = nullptr; + IOFamily(std::string desc, std::string grp, mfem::ParGridFunction *pf) : description_(desc), group_(grp), pfunc_(pf) {} @@ -75,7 +82,7 @@ class IODataOrganizer { std::map> vars_; // solution var info for each IO family void registerIOFamily(std::string description, std::string group, mfem::ParGridFunction *pfunc, - bool auxRestart = true, bool inRestartFile = true); + bool auxRestart = true, bool inRestartFile = true, mfem::FiniteElementCollection *fec = NULL); ~IODataOrganizer(); void registerIOVar(std::string group, std::string varName, int index, bool inRestartFile = true); @@ -86,6 +93,10 @@ class IODataOrganizer { void write(hid_t file, bool rank0); void writeSerial(hid_t file, MPI_Comm comm, int local_ne, int global_ne, const int *locToGlobElem, const Array &partitioning); + + void read(hid_t file, bool rank0); + void readSerial(hid_t file, bool rank0, MPI_Groups *groupsMPI, Array partitioning, int global_ne); + void readChangeOrder(hid_t file, bool rank0, int read_order); }; void read_partitioned_soln_data(hid_t file, std::string varName, size_t index, double *data); From 975a659c4a4f8e8ae44690cabba47d2f51069e2d Mon Sep 17 00:00:00 2001 From: "Todd A. Oliver" Date: Wed, 31 Jan 2024 12:39:07 -0600 Subject: [PATCH 07/14] Fold serial write capability in to IODataOrganizer::write function (#244) The distinctions between write and writeSerial are small, so it is better to handle both cases through a single function, which we do now through IODataOrganizer::write. The IODataOrganizer::writeSerial function is no longer necessary and is therefore eliminated. --- src/M2ulPhyS.cpp | 2 +- src/io.cpp | 166 +++++++++++++++++++++-------------------------- src/io.hpp | 19 ++++-- 3 files changed, 87 insertions(+), 100 deletions(-) diff --git a/src/M2ulPhyS.cpp b/src/M2ulPhyS.cpp index ac820a276..8212d9ebb 100644 --- a/src/M2ulPhyS.cpp +++ b/src/M2ulPhyS.cpp @@ -655,7 +655,7 @@ void M2ulPhyS::initVariables() { ioData.registerIOVar("/rmsData", "vw", 5); } - ioData.initializeSerial(rank0_, (config.RestartSerial() != "no"), serial_mesh); + ioData.initializeSerial(rank0_, (config.RestartSerial() != "no"), serial_mesh, locToGlobElem, &partitioning_); projectInitialSolution(); // Boundary attributes in present partition diff --git a/src/io.cpp b/src/io.cpp index b217d18c0..bd8e56fb9 100644 --- a/src/io.cpp +++ b/src/io.cpp @@ -92,11 +92,7 @@ void M2ulPhyS::write_restart_files_hdf5(hid_t file, bool serialized_write) { // ------------------------------------------------------------------- // Data - actual data write handled by IODataOrganizer class // ------------------------------------------------------------------- - if (serialized_write) { - ioData.writeSerial(file, groupsMPI->getTPSCommWorld(), mesh->GetNE(), nelemGlobal_, locToGlobElem, partitioning_); - } else { - ioData.write(file, rank0_); - } + ioData.write(file, serialized_write); } void M2ulPhyS::read_restart_files_hdf5(hid_t file, bool serialized_read) { @@ -468,7 +464,7 @@ void read_serialized_soln_data(hid_t file, string varName, int numDof, int varOf } // convenience function to write HDF5 data -void write_soln_data(hid_t group, string varName, hid_t dataspace, double *data, bool rank0) { +void write_soln_data(hid_t group, string varName, hid_t dataspace, const double *data, bool rank0) { hid_t data_soln; herr_t status; assert(group >= 0); @@ -532,17 +528,20 @@ void M2ulPhyS::readTable(const std::string &inputPath, TableInput &result) { // Routines for I/O data organizer helper class // --------------------------------------------- -void IOFamily::serializeForWrite(MPI_Comm comm, int local_ne, int global_ne, const int *locToGlobElem, - const Array &partitioning) { - int rank; - MPI_Comm_rank(comm, &rank); +void IOFamily::serializeForWrite() { + MPI_Comm comm = this->pfunc_->ParFESpace()->GetComm(); + const int rank = this->pfunc_->ParFESpace()->GetMyRank(); const bool rank0 = (rank == 0); + const int local_ne = this->pfunc_->ParFESpace()->GetParMesh()->GetNE(); + + const Array &partitioning = *(this->partitioning_); + const int *locToGlobElem = this->local_to_global_elem_; + // Ensure consistency - int global_ne_check; - MPI_Reduce(&local_ne, &global_ne_check, 1, MPI_INT, MPI_SUM, 0, comm); + int global_ne; + MPI_Reduce(&local_ne, &global_ne, 1, MPI_INT, MPI_SUM, 0, comm); if (rank0) { - assert(global_ne_check == global_ne); assert(partitioning.Size() == global_ne); } assert(locToGlobElem != NULL); @@ -624,119 +623,102 @@ int IODataOrganizer::getIOFamilyIndex(std::string group) { return (-1); } -void IODataOrganizer::initializeSerial(bool root, bool serial, Mesh *serial_mesh) { +void IODataOrganizer::initializeSerial(bool root, bool serial, Mesh *serial_mesh, int *locToGlob, Array *part) { + supports_serial_ = serial; + // loop through families for (size_t n = 0; n < families_.size(); n++) { IOFamily &fam = families_[n]; fam.serial_fes = NULL; fam.serial_sol = NULL; - if (root && serial) { - const FiniteElementCollection *fec = fam.pfunc_->ParFESpace()->FEColl(); - int numVars = fam.pfunc_->Size() / fam.pfunc_->ParFESpace()->GetNDofs(); - fam.serial_fes = new FiniteElementSpace(serial_mesh, fec, numVars, Ordering::byNODES); - fam.serial_sol = new GridFunction(fam.serial_fes); - // cout<<"I/O organizer for group "<ParFESpace()->FEColl(); + int numVars = fam.pfunc_->Size() / fam.pfunc_->ParFESpace()->GetNDofs(); + + fam.serial_fes = new FiniteElementSpace(serial_mesh, fec, numVars, Ordering::byNODES); + fam.serial_sol = new GridFunction(fam.serial_fes); + // cout<<"I/O organizer for group "<ParFESpace()->GetNDofs(); +void IODataOrganizer::write(hid_t file, bool serial) { + if (serial) assert(supports_serial_); - // define groups based on defined IO families - if (rank0) { - grvy_printf(ginfo, "\nCreating HDF5 group for defined IO families\n"); - grvy_printf(ginfo, "--> %s : %s\n", fam.group_.c_str(), fam.description_.c_str()); + // NOTE(trevilo): this loop envisions families that have different + // number of mpi ranks, but rest of the code is not general enough + // to support this configuration yet (or at a minimum there are no + // tests of such a set up) + int nprocs_max = 0; + for (auto fam : families_) { + int nprocs_tmp = fam.pfunc_->ParFESpace()->GetNRanks(); + if (nprocs_tmp > nprocs_max) { + nprocs_max = nprocs_tmp; } + } + assert(nprocs_max > 0); - hid_t group = -1; - hid_t dataspace = -1; + // Do we need to serialize the data prior to the write? + const bool require_serialization = (serial && nprocs_max > 1); - dataspace = H5Screate_simple(1, dims, NULL); - assert(dataspace >= 0); - group = H5Gcreate(file, fam.group_.c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); - assert(group >= 0); - - // get pointer to raw data - double *data = fam.pfunc_->HostReadWrite(); + // Loop over defined IO families to save desired output + for (auto fam : families_) { + if (require_serialization) fam.serializeForWrite(); // get defined variables for this IO family vector vars = vars_[fam.group_]; - // save raw data - for (auto var : vars) write_soln_data(group, var.varName_, dataspace, data + var.index_ * dims[0], rank0); + // determine if rank 0 + const int rank = fam.pfunc_->ParFESpace()->GetMyRank(); + const bool rank0 = (rank == 0); - if (group >= 0) H5Gclose(group); - if (dataspace >= 0) H5Sclose(dataspace); - } -} - -void IODataOrganizer::writeSerial(hid_t file, MPI_Comm comm, int local_ne, int global_ne, const int *locToGlobElem, - const Array &partitioning) { - int rank; - MPI_Comm_rank(comm, &rank); - int nprocs; - MPI_Comm_size(comm, &nprocs); - const bool rank0 = (rank == 0); - - // if running serially, writeSerial and write are logically the - // same, but infrastructure is different. Thus, for serial run, - // just call write() here - if (nprocs == 1) { - this->write(file, rank0); - return; - } - - hsize_t dims[1]; - - // Loop over defined IO families to save desired output - for (size_t n = 0; n < families_.size(); n++) { - IOFamily &fam = families_[n]; + hsize_t dims[1]; + hid_t group = -1; + hid_t dataspace = -1; if (rank0) { - assert((locToGlobElem != NULL) && (partitioning.Size() == global_ne)); - assert(fam.serial_fes != NULL); - dims[0] = fam.serial_fes->GetNDofs(); grvy_printf(ginfo, "\nCreating HDF5 group for defined IO families\n"); grvy_printf(ginfo, "--> %s : %s\n", fam.group_.c_str(), fam.description_.c_str()); } - hid_t group = -1; - hid_t dataspace = -1; + if (require_serialization) { + // if require_serialization, then we just populated fam.serial_sol above, and now we use only rank0 to write + if (rank0) { + assert(fam.serial_fes != NULL); + dims[0] = fam.serial_fes->GetNDofs(); - if (rank0) { + dataspace = H5Screate_simple(1, dims, NULL); + assert(dataspace >= 0); + group = H5Gcreate(file, fam.group_.c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + assert(group >= 0); + + // get pointer to raw data + assert(fam.serial_sol != NULL); + const double *data = fam.serial_sol->HostRead(); + + // save raw data + for (auto var : vars) write_soln_data(group, var.varName_, dataspace, data + var.index_ * dims[0], rank0); + } + } else { + // otherwise, use all ranks to write data from fam.pfunc_ + dims[0] = fam.pfunc_->ParFESpace()->GetNDofs(); dataspace = H5Screate_simple(1, dims, NULL); assert(dataspace >= 0); group = H5Gcreate(file, fam.group_.c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); assert(group >= 0); - } - fam.serializeForWrite(comm, local_ne, global_ne, locToGlobElem, partitioning); - - // get pointer to raw data - double *data = nullptr; - if (rank0) data = fam.serial_sol->HostReadWrite(); - - // get defined variables for this IO family - vector vars = vars_[fam.group_]; + // get pointer to raw data + const double *data = fam.pfunc_->HostRead(); - // save raw data - if (rank0) { + // save raw data for (auto var : vars) write_soln_data(group, var.varName_, dataspace, data + var.index_ * dims[0], rank0); } - if (group >= 0) H5Gclose(group); if (dataspace >= 0) H5Sclose(dataspace); } diff --git a/src/io.hpp b/src/io.hpp index 38b95ba5a..048a8218c 100644 --- a/src/io.hpp +++ b/src/io.hpp @@ -56,6 +56,11 @@ class IOFamily { mfem::FiniteElementSpace *serial_fes = nullptr; mfem::GridFunction *serial_sol = nullptr; + // additional data used for serial read/write (needed to properly + // serialize fields in write or properly distribute data in read) + int *local_to_global_elem_ = nullptr; // maps from local elem index to global + mfem::Array *partitioning_ = nullptr; // partitioning[i] mpi rank for element i (in global numbering) + // Variables used for "auxilliary" read (i.e., restart from different order soln) mfem::FiniteElementCollection *fec_ = nullptr; // collection used to instantiate pfunc_ mfem::FiniteElementCollection *aux_fec_ = nullptr; @@ -65,8 +70,7 @@ class IOFamily { IOFamily(std::string desc, std::string grp, mfem::ParGridFunction *pf) : description_(desc), group_(grp), pfunc_(pf) {} - void serializeForWrite(MPI_Comm comm, int local_ne, int global_ne, const int *locToGlobElem, - const Array &partitioning); + void serializeForWrite(); }; class IOVar { @@ -77,6 +81,9 @@ class IOVar { }; class IODataOrganizer { + protected: + bool supports_serial_ = false; + public: std::vector families_; // registered IO families std::map> vars_; // solution var info for each IO family @@ -88,11 +95,9 @@ class IODataOrganizer { void registerIOVar(std::string group, std::string varName, int index, bool inRestartFile = true); int getIOFamilyIndex(std::string group); - void initializeSerial(bool root, bool serial, mfem::Mesh *serial_mesh); + void initializeSerial(bool root, bool serial, mfem::Mesh *serial_mesh, int *locToGlob, mfem::Array *part); - void write(hid_t file, bool rank0); - void writeSerial(hid_t file, MPI_Comm comm, int local_ne, int global_ne, const int *locToGlobElem, - const Array &partitioning); + void write(hid_t file, bool serial); void read(hid_t file, bool rank0); void readSerial(hid_t file, bool rank0, MPI_Groups *groupsMPI, Array partitioning, int global_ne); @@ -102,7 +107,7 @@ class IODataOrganizer { void read_partitioned_soln_data(hid_t file, std::string varName, size_t index, double *data); void read_serialized_soln_data(hid_t file, std::string varName, int numDof, int varOffset, double *data, IOFamily &fam, MPI_Groups *groupsMPI, mfem::Array partitioning, int nelemGlobal); -void write_soln_data(hid_t group, std::string varName, hid_t dataspace, double *data, bool rank0); +void write_soln_data(hid_t group, std::string varName, hid_t dataspace, const double *data, bool rank0); void partitioning_file_hdf5(std::string mode, const RunConfiguration &config, MPI_Groups *groupsMPI, int nelemGlobal, mfem::Array &partitioning); #endif // IO_HPP_ From debec1dc5fb0b7d55cd4404a138118111bd866ce Mon Sep 17 00:00:00 2001 From: "Todd A. Oliver" Date: Wed, 31 Jan 2024 14:30:47 -0600 Subject: [PATCH 08/14] Consolidate IODataOrganizer read capabilities (#244) Refactor functions in IODataOrganizer so that capabilities previously distributed between read, readSerial, and readChangeOrder are all now handled by a single read function. --- src/io.cpp | 331 ++++++++++++++++++++++++++--------------------------- src/io.hpp | 7 +- 2 files changed, 163 insertions(+), 175 deletions(-) diff --git a/src/io.cpp b/src/io.cpp index bd8e56fb9..8c13a3a96 100644 --- a/src/io.cpp +++ b/src/io.cpp @@ -146,18 +146,7 @@ void M2ulPhyS::read_restart_files_hdf5(hid_t file, bool serialized_read) { // ------------------------------------------------------------------- // Data - actual data read handled by IODataOrganizer class // ------------------------------------------------------------------- - if (!loadFromAuxSol && !serialized_read) { - assert(read_order == order); - ioData.read(file, rank0_); - } else if (!loadFromAuxSol && serialized_read) { - assert(read_order == order); - ioData.readSerial(file, rank0_, groupsMPI, partitioning_, nelemGlobal_); - } else if (loadFromAuxSol && !serialized_read) { - ioData.readChangeOrder(file, rank0_, read_order); - } else { - if (rank0_) grvy_printf(gerror, "[ERROR]: Unsupported read mode requested.\n"); - exit(ERROR); - } + ioData.read(file, serialized_read, read_order); if (loadFromAuxSol) { // Update primitive variables. This will be done automatically @@ -385,11 +374,17 @@ void read_partitioned_soln_data(hid_t file, string varName, size_t index, double // convenience function to read and distribute solution data for serialized restarts void read_serialized_soln_data(hid_t file, string varName, int numDof, int varOffset, double *data, IOFamily &fam, - MPI_Groups *groupsMPI, Array partitioning, int nelemGlobal) { - MPI_Comm TPSCommWorld = groupsMPI->getTPSCommWorld(); - const bool rank0 = groupsMPI->isWorldRoot(); - const int nprocs = groupsMPI->getTPSWorldSize(); - const int rank = groupsMPI->getTPSWorldRank(); + Array partitioning) { + MPI_Comm comm = fam.pfunc_->ParFESpace()->GetComm(); + const int rank = fam.pfunc_->ParFESpace()->GetMyRank(); + const bool rank0 = (rank == 0); + const int nprocs = fam.pfunc_->ParFESpace()->GetNRanks(); + + // Ensure number of elements is consistent with partition array + const int local_ne = fam.pfunc_->ParFESpace()->GetParMesh()->GetNE(); + int global_ne; + MPI_Allreduce(&local_ne, &global_ne, 1, MPI_INT, MPI_SUM, comm); + assert(partitioning.Size() == global_ne); hid_t data_soln; herr_t status; @@ -411,14 +406,12 @@ void read_serialized_soln_data(hid_t file, string varName, int numDof, int varOf H5Dclose(data_soln); // assign solution owned by rank 0 - assert(partitioning != NULL); - Array lvdofs, gvdofs; Vector lnodes; int counter = 0; // int ndof_per_elem; - for (int gelem = 0; gelem < nelemGlobal; gelem++) { + for (int gelem = 0; gelem < global_ne; gelem++) { if (partitioning[gelem] == 0) { // cull out subset of local vdof vector to use for this solution var fam.pfunc_->ParFESpace()->GetElementVDofs(counter, lvdofs); @@ -440,14 +433,14 @@ void read_serialized_soln_data(hid_t file, string varName, int numDof, int varOf // pack remaining data and send to other processors for (int rank = 1; rank < nprocs; rank++) { std::vector packedData; - for (int gelem = 0; gelem < nelemGlobal; gelem++) + for (int gelem = 0; gelem < global_ne; gelem++) if (partitioning[gelem] == rank) { fam.serial_fes->GetElementVDofs(gelem, gvdofs); int numDof_per_this_elem = gvdofs.Size() / numStateVars; for (int i = 0; i < numDof_per_this_elem; i++) packedData.push_back(data_serial[gvdofs[i]]); } int tag = 20; - MPI_Send(packedData.data(), packedData.size(), MPI_DOUBLE, rank, tag, TPSCommWorld); + MPI_Send(packedData.data(), packedData.size(), MPI_DOUBLE, rank, tag, comm); } } else { // <-- end rank 0 int numlDofs = fam.pfunc_->ParFESpace()->GetNDofs(); @@ -456,7 +449,7 @@ void read_serialized_soln_data(hid_t file, string varName, int numDof, int varOf std::vector packedData(numlDofs); // receive solution data from rank 0 - MPI_Recv(packedData.data(), numlDofs, MPI_DOUBLE, 0, tag, TPSCommWorld, MPI_STATUS_IGNORE); + MPI_Recv(packedData.data(), numlDofs, MPI_DOUBLE, 0, tag, comm, MPI_STATUS_IGNORE); // update local state vector for (int i = 0; i < numlDofs; i++) data[i + varOffset * numlDofs] = packedData[i]; @@ -724,185 +717,183 @@ void IODataOrganizer::write(hid_t file, bool serial) { } } -void IODataOrganizer::read(hid_t file, bool rank0) { - hsize_t numInSoln = 0; - hid_t data_soln; +void IODataOrganizer::read(hid_t file, bool serial, int read_order) { + if (serial) assert(supports_serial_); - // Loop over defined IO families to save desired output + // NOTE(trevilo): this loop envisions families that have different + // number of mpi ranks, but rest of the code is not general enough + // to support this configuration yet (or at a minimum there are no + // tests of such a set up) + int nprocs_max = 0; for (auto fam : families_) { - if (fam.inRestartFile) { // read mode - if (rank0) cout << "Reading in solutiond data from restart..." << endl; - - // verify Dofs match expectations with current mesh - hid_t dataspace; - - vector vars = vars_[fam.group_]; - string varGroupName = fam.group_; - varGroupName.append("/"); - varGroupName.append(vars[0].varName_); - - data_soln = H5Dopen2(file, varGroupName.c_str(), H5P_DEFAULT); - assert(data_soln >= 0); - dataspace = H5Dget_space(data_soln); - numInSoln = H5Sget_simple_extent_npoints(dataspace); - H5Dclose(data_soln); - - int dof = fam.pfunc_->ParFESpace()->GetNDofs(); - assert((int)numInSoln == dof); - - // get pointer to raw data - double *data = fam.pfunc_->HostReadWrite(); - - // read from file into appropriate spot in data - for (auto var : vars) { - if (var.inRestartFile_) { - std::string h5Path = fam.group_ + "/" + var.varName_; - if (rank0) grvy_printf(ginfo, "--> Reading h5 path = %s\n", h5Path.c_str()); - read_partitioned_soln_data(file, h5Path.c_str(), var.index_ * numInSoln, data); - } - } + int nprocs_tmp = fam.pfunc_->ParFESpace()->GetNRanks(); + if (nprocs_tmp > nprocs_max) { + nprocs_max = nprocs_tmp; } } -} + assert(nprocs_max > 0); -void IODataOrganizer::readSerial(hid_t file, bool rank0, MPI_Groups *groupsMPI, Array partitioning, - int global_ne) { - MPI_Comm TPSCommWorld = groupsMPI->getTPSCommWorld(); - // ------------------------------------------------------------------- - // Read/write solution data defined by IO families - // ------------------------------------------------------------------- hsize_t numInSoln = 0; hid_t data_soln; - //------------------------------------------------------- - // Loop over defined IO families to save desired output - //------------------------------------------------------- - for (size_t n = 0; n < families_.size(); n++) { - IOFamily &fam = families_[n]; + // Loop over defined IO families to load desired input + for (auto fam : families_) { if (fam.inRestartFile) { // read mode + const int rank = fam.pfunc_->ParFESpace()->GetMyRank(); + const bool rank0 = (rank == 0); + if (rank0) cout << "Reading in solutiond data from restart..." << endl; - // verify Dofs match expectations with current mesh - if (rank0) { - hid_t dataspace; - - vector vars = vars_[fam.group_]; - string varGroupName = fam.group_; - varGroupName.append("/"); - varGroupName.append(vars[0].varName_); - - data_soln = H5Dopen2(file, varGroupName.c_str(), H5P_DEFAULT); - assert(data_soln >= 0); - dataspace = H5Dget_space(data_soln); - numInSoln = H5Sget_simple_extent_npoints(dataspace); - H5Dclose(data_soln); + int fam_order; + bool change_order = false; + if (read_order >= 0) { + assert(!fam.pfunc_->ParFESpace()->IsVariableOrder()); + fam_order = fam.pfunc_->ParFESpace()->FEColl()->GetOrder(); + change_order = (fam_order != read_order); } - int dof = fam.pfunc_->ParFESpace()->GetNDofs(); + if (change_order) { + if (serial) { + if (rank0) grvy_printf(gerror, "[ERROR]: Serial read does not support order change.\n"); + exit(ERROR); + } - int dof_global; - MPI_Reduce(&dof, &dof_global, 1, MPI_INT, MPI_SUM, 0, TPSCommWorld); - dof = dof_global; + if (fam.allowsAuxRestart) { // read mode + double *aux_U_data = NULL; - if (rank0) assert((int)numInSoln == dof); + // verify Dofs match expectations with current mesh + hid_t dataspace; - // get pointer to raw data - double *data = fam.pfunc_->HostReadWrite(); - - vector vars = vars_[fam.group_]; - for (auto var : vars) { - if (var.inRestartFile_) { - std::string h5Path = fam.group_ + "/" + var.varName_; - if (rank0) grvy_printf(ginfo, "--> Reading h5 path = %s\n", h5Path.c_str()); - assert(partitioning.Size() == global_ne); - read_serialized_soln_data(file, h5Path.c_str(), dof, var.index_, data, fam, groupsMPI, partitioning, - global_ne); - } - } - } - } // end loop over families -} + vector vars = vars_[fam.group_]; + string varGroupName = fam.group_; + varGroupName.append("/"); + varGroupName.append(vars[0].varName_); -void IODataOrganizer::readChangeOrder(hid_t file, bool rank0, int read_order) { - hsize_t numInSoln = 0; - hid_t data_soln; - double *aux_U_data = NULL; + data_soln = H5Dopen2(file, varGroupName.c_str(), H5P_DEFAULT); + assert(data_soln >= 0); + dataspace = H5Dget_space(data_soln); + numInSoln = H5Sget_simple_extent_npoints(dataspace); + H5Dclose(data_soln); - // NOTE(trevilo): this assert is left over from the original - // implementation. I believe the refactored implementation can - // support multiple IOFamily objects and changing order on restart. - // But, since we do not have a use or test case, I am leaving this - // assert for the time being. - assert(families_.size() == 1); + assert(fam.fec_ != NULL); + fam.aux_fec_ = fam.fec_->Clone(read_order); - // Loop over defined IO families to save desired output - for (size_t n = 0; n < families_.size(); n++) { - IOFamily &fam = families_[n]; - if (fam.inRestartFile && fam.allowsAuxRestart) { // read mode - if (rank0) cout << "Reading in solutiond data from restart..." << endl; + int nvars = vars_[fam.group_].size(); + fam.aux_pfes_ = + new ParFiniteElementSpace(fam.pfunc_->ParFESpace()->GetParMesh(), fam.aux_fec_, nvars, Ordering::byNODES); - // verify Dofs match expectations with current mesh - hid_t dataspace; + aux_U_data = new double[nvars * fam.aux_pfes_->GetNDofs()]; + fam.aux_pfunc_ = new ParGridFunction(fam.aux_pfes_, aux_U_data); - vector vars = vars_[fam.group_]; - string varGroupName = fam.group_; - varGroupName.append("/"); - varGroupName.append(vars[0].varName_); + int dof = fam.aux_pfes_->GetNDofs(); - data_soln = H5Dopen2(file, varGroupName.c_str(), H5P_DEFAULT); - assert(data_soln >= 0); - dataspace = H5Dget_space(data_soln); - numInSoln = H5Sget_simple_extent_npoints(dataspace); - H5Dclose(data_soln); + assert((int)numInSoln == dof); - assert(fam.fec_ != NULL); - fam.aux_fec_ = fam.fec_->Clone(read_order); + // get pointer to raw data + double *data = fam.aux_pfunc_->HostReadWrite(); - int nvars = vars_[fam.group_].size(); - // fam.aux_pfes_ = new ParFiniteElementSpace(mesh, fam.aux_fec_, nvars, Ordering::byNODES); - fam.aux_pfes_ = - new ParFiniteElementSpace(fam.pfunc_->ParFESpace()->GetParMesh(), fam.aux_fec_, nvars, Ordering::byNODES); + for (auto var : vars) { + if (var.inRestartFile_) { + std::string h5Path = fam.group_ + "/" + var.varName_; + if (rank0) grvy_printf(ginfo, "--> Reading h5 path = %s\n", h5Path.c_str()); + read_partitioned_soln_data(file, h5Path.c_str(), var.index_ * numInSoln, data); + } + } - aux_U_data = new double[nvars * fam.aux_pfes_->GetNDofs()]; - fam.aux_pfunc_ = new ParGridFunction(fam.aux_pfes_, aux_U_data); + if (rank0) cout << "Interpolating from read_order = " << read_order << " to solution space" << endl; - int dof = fam.aux_pfes_->GetNDofs(); + // Interpolate from aux to new order + fam.pfunc_->ProjectGridFunction(*fam.aux_pfunc_); - assert((int)numInSoln == dof); + // If interpolation was successful, the L2 norm of the + // difference between the auxOrder and order versions should be + // zero (well... to within round-off-induced error) + VectorGridFunctionCoefficient lowOrderCoeff(fam.aux_pfunc_); + double err = fam.pfunc_->ComputeL2Error(lowOrderCoeff); - // get pointer to raw data - double *data = fam.aux_pfunc_->HostReadWrite(); + if (rank0) cout << "|| interpolation error ||_2 = " << err << endl; - for (auto var : vars) { - if (var.inRestartFile_) { - std::string h5Path = fam.group_ + "/" + var.varName_; - if (rank0) grvy_printf(ginfo, "--> Reading h5 path = %s\n", h5Path.c_str()); - read_partitioned_soln_data(file, h5Path.c_str(), var.index_ * numInSoln, data); + // clean up aux data + delete fam.aux_pfunc_; + delete[] aux_U_data; + delete fam.aux_pfes_; + delete fam.aux_fec_; + } else { + if (rank0) cout << "Skipping family " << fam.group_ << " doesn't allow changing order." << endl; + } + } else { + if (serial && nprocs_max > 1) { + const Array &partitioning = *(fam.partitioning_); + MPI_Comm comm = fam.pfunc_->ParFESpace()->GetComm(); + + // verify Dofs match expectations with current mesh + if (rank0) { + hid_t dataspace; + + vector vars = vars_[fam.group_]; + string varGroupName = fam.group_; + varGroupName.append("/"); + varGroupName.append(vars[0].varName_); + + data_soln = H5Dopen2(file, varGroupName.c_str(), H5P_DEFAULT); + assert(data_soln >= 0); + dataspace = H5Dget_space(data_soln); + numInSoln = H5Sget_simple_extent_npoints(dataspace); + H5Dclose(data_soln); + } + + int dof = fam.pfunc_->ParFESpace()->GetNDofs(); + + int dof_global; + MPI_Reduce(&dof, &dof_global, 1, MPI_INT, MPI_SUM, 0, comm); + dof = dof_global; + + if (rank0) assert((int)numInSoln == dof); + + // get pointer to raw data + double *data = fam.pfunc_->HostReadWrite(); + + vector vars = vars_[fam.group_]; + for (auto var : vars) { + if (var.inRestartFile_) { + std::string h5Path = fam.group_ + "/" + var.varName_; + if (rank0) grvy_printf(ginfo, "--> Reading h5 path = %s\n", h5Path.c_str()); + read_serialized_soln_data(file, h5Path.c_str(), dof, var.index_, data, fam, partitioning); + } + } + } else { + // verify Dofs match expectations with current mesh + hid_t dataspace; + + vector vars = vars_[fam.group_]; + string varGroupName = fam.group_; + varGroupName.append("/"); + varGroupName.append(vars[0].varName_); + + data_soln = H5Dopen2(file, varGroupName.c_str(), H5P_DEFAULT); + assert(data_soln >= 0); + dataspace = H5Dget_space(data_soln); + numInSoln = H5Sget_simple_extent_npoints(dataspace); + H5Dclose(data_soln); + + int dof = fam.pfunc_->ParFESpace()->GetNDofs(); + assert((int)numInSoln == dof); + + // get pointer to raw data + double *data = fam.pfunc_->HostReadWrite(); + + // read from file into appropriate spot in data + for (auto var : vars) { + if (var.inRestartFile_) { + std::string h5Path = fam.group_ + "/" + var.varName_; + if (rank0) grvy_printf(ginfo, "--> Reading h5 path = %s\n", h5Path.c_str()); + read_partitioned_soln_data(file, h5Path.c_str(), var.index_ * numInSoln, data); + } + } } } - - if (rank0) cout << "Interpolating from read_order = " << read_order << " to solution space" << endl; - - // Interpolate from aux to new order - fam.pfunc_->ProjectGridFunction(*fam.aux_pfunc_); - - // If interpolation was successful, the L2 norm of the - // difference between the auxOrder and order versions should be - // zero (well... to within round-off-induced error) - VectorGridFunctionCoefficient lowOrderCoeff(fam.aux_pfunc_); - double err = fam.pfunc_->ComputeL2Error(lowOrderCoeff); - - if (rank0) cout << "|| interpolation error ||_2 = " << err << endl; - - // clean up aux data - delete fam.aux_pfunc_; - delete[] aux_U_data; - delete fam.aux_pfes_; - delete fam.aux_fec_; - } else { - if (rank0) cout << "Skipping family " << fam.group_ << " in readChangeOrder" << endl; } - } // end loop over families + } } IODataOrganizer::~IODataOrganizer() { diff --git a/src/io.hpp b/src/io.hpp index 048a8218c..af8ebce0f 100644 --- a/src/io.hpp +++ b/src/io.hpp @@ -98,15 +98,12 @@ class IODataOrganizer { void initializeSerial(bool root, bool serial, mfem::Mesh *serial_mesh, int *locToGlob, mfem::Array *part); void write(hid_t file, bool serial); - - void read(hid_t file, bool rank0); - void readSerial(hid_t file, bool rank0, MPI_Groups *groupsMPI, Array partitioning, int global_ne); - void readChangeOrder(hid_t file, bool rank0, int read_order); + void read(hid_t file, bool serial, int read_order = -1); }; void read_partitioned_soln_data(hid_t file, std::string varName, size_t index, double *data); void read_serialized_soln_data(hid_t file, std::string varName, int numDof, int varOffset, double *data, IOFamily &fam, - MPI_Groups *groupsMPI, mfem::Array partitioning, int nelemGlobal); + mfem::Array partitioning); void write_soln_data(hid_t group, std::string varName, hid_t dataspace, const double *data, bool rank0); void partitioning_file_hdf5(std::string mode, const RunConfiguration &config, MPI_Groups *groupsMPI, int nelemGlobal, mfem::Array &partitioning); From 6dbe4b3a5693c048bc113ed8567fce69d4a9881f Mon Sep 17 00:00:00 2001 From: "Todd A. Oliver" Date: Wed, 31 Jan 2024 22:01:22 -0600 Subject: [PATCH 09/14] For readability, push some read methods down into IOFamily (#244) The logic in IODataOrganizer::read is hard to follow in big chunks, so move those big chunks into the IOFamily class, which is where they belong. Also, move IOVar vector into IOFamily, which simplifies the code slightly, and just makes more sense than storing the info in IODataOrganizer where we have to map from the family name to the variables. --- src/io.cpp | 306 +++++++++++++++++++++++++++-------------------------- src/io.hpp | 27 +++-- 2 files changed, 172 insertions(+), 161 deletions(-) diff --git a/src/io.cpp b/src/io.cpp index 8c13a3a96..53696f237 100644 --- a/src/io.cpp +++ b/src/io.cpp @@ -521,10 +521,13 @@ void M2ulPhyS::readTable(const std::string &inputPath, TableInput &result) { // Routines for I/O data organizer helper class // --------------------------------------------- +IOFamily::IOFamily(std::string desc, std::string grp, mfem::ParGridFunction *pf) + : description_(desc), group_(grp), pfunc_(pf) { + rank0_ = (pfunc_->ParFESpace()->GetMyRank() == 0); +} + void IOFamily::serializeForWrite() { MPI_Comm comm = this->pfunc_->ParFESpace()->GetComm(); - const int rank = this->pfunc_->ParFESpace()->GetMyRank(); - const bool rank0 = (rank == 0); const int local_ne = this->pfunc_->ParFESpace()->GetParMesh()->GetNE(); @@ -534,13 +537,13 @@ void IOFamily::serializeForWrite() { // Ensure consistency int global_ne; MPI_Reduce(&local_ne, &global_ne, 1, MPI_INT, MPI_SUM, 0, comm); - if (rank0) { + if (rank0_) { assert(partitioning.Size() == global_ne); } assert(locToGlobElem != NULL); ParGridFunction *pfunc = this->pfunc_; - if (rank0) { + if (rank0_) { grvy_printf(ginfo, "Generating serialized restart file (group %s...)\n", this->group_.c_str()); // copy my own data Array lvdofs, gvdofs; @@ -582,11 +585,150 @@ void IOFamily::serializeForWrite() { } } +void IOFamily::readPartitioned(hid_t file) { + hid_t dataspace; + hsize_t numInSoln = 0; + hid_t data_soln; + + string varGroupName = group_; + varGroupName.append("/"); + varGroupName.append(vars_[0].varName_); + + data_soln = H5Dopen2(file, varGroupName.c_str(), H5P_DEFAULT); + assert(data_soln >= 0); + dataspace = H5Dget_space(data_soln); + numInSoln = H5Sget_simple_extent_npoints(dataspace); + H5Dclose(data_soln); + + int dof = pfunc_->ParFESpace()->GetNDofs(); + + // verify Dofs match expectations with current mesh + assert((int)numInSoln == dof); + + // get pointer to raw data + double *data = pfunc_->HostWrite(); + + // read from file into appropriate spot in data + for (auto var : vars_) { + if (var.inRestartFile_) { + std::string h5Path = group_ + "/" + var.varName_; + if (rank0_) grvy_printf(ginfo, "--> Reading h5 path = %s\n", h5Path.c_str()); + read_partitioned_soln_data(file, h5Path.c_str(), var.index_ * numInSoln, data); + } + } +} + +void IOFamily::readSerial(hid_t file) { + const Array &partitioning = *partitioning_; + MPI_Comm comm = pfunc_->ParFESpace()->GetComm(); + + hsize_t numInSoln = 0; + hid_t data_soln; + + // verify Dofs match expectations with current mesh + if (rank0_) { + hid_t dataspace; + + string varGroupName = group_; + varGroupName.append("/"); + varGroupName.append(vars_[0].varName_); + + data_soln = H5Dopen2(file, varGroupName.c_str(), H5P_DEFAULT); + assert(data_soln >= 0); + dataspace = H5Dget_space(data_soln); + numInSoln = H5Sget_simple_extent_npoints(dataspace); + H5Dclose(data_soln); + } + + int dof = pfunc_->ParFESpace()->GetNDofs(); + + int dof_global; + MPI_Reduce(&dof, &dof_global, 1, MPI_INT, MPI_SUM, 0, comm); + dof = dof_global; + + if (rank0_) assert((int)numInSoln == dof); + + // get pointer to raw data + double *data = pfunc_->HostWrite(); + + for (auto var : vars_) { + if (var.inRestartFile_) { + std::string h5Path = group_ + "/" + var.varName_; + if (rank0_) grvy_printf(ginfo, "--> Reading h5 path = %s\n", h5Path.c_str()); + read_serialized_soln_data(file, h5Path.c_str(), dof, var.index_, data, *this, partitioning); + } + } +} + +void IOFamily::readChangeOrder(hid_t file, int read_order) { + if (allowsAuxRestart) { + hid_t dataspace; + hsize_t numInSoln = 0; + hid_t data_soln; + + double *aux_U_data = NULL; + + string varGroupName = group_; + varGroupName.append("/"); + varGroupName.append(vars_[0].varName_); + + data_soln = H5Dopen2(file, varGroupName.c_str(), H5P_DEFAULT); + assert(data_soln >= 0); + dataspace = H5Dget_space(data_soln); + numInSoln = H5Sget_simple_extent_npoints(dataspace); + H5Dclose(data_soln); + + assert(fec_ != NULL); + aux_fec_ = fec_->Clone(read_order); + + int nvars = vars_.size(); + aux_pfes_ = new ParFiniteElementSpace(pfunc_->ParFESpace()->GetParMesh(), aux_fec_, nvars, Ordering::byNODES); + + aux_U_data = new double[nvars * aux_pfes_->GetNDofs()]; + aux_pfunc_ = new ParGridFunction(aux_pfes_, aux_U_data); + + int dof = aux_pfes_->GetNDofs(); + + assert((int)numInSoln == dof); + + // get pointer to raw data + double *data = aux_pfunc_->HostWrite(); + + for (auto var : vars_) { + if (var.inRestartFile_) { + std::string h5Path = group_ + "/" + var.varName_; + if (rank0_) grvy_printf(ginfo, "--> Reading h5 path = %s\n", h5Path.c_str()); + read_partitioned_soln_data(file, h5Path.c_str(), var.index_ * numInSoln, data); + } + } + + if (rank0_) cout << "Interpolating from read_order = " << read_order << " to solution space" << endl; + + // Interpolate from aux to new order + pfunc_->ProjectGridFunction(*aux_pfunc_); + + // If interpolation was successful, the L2 norm of the + // difference between the auxOrder and order versions should be + // zero (well... to within round-off-induced error) + VectorGridFunctionCoefficient lowOrderCoeff(aux_pfunc_); + double err = pfunc_->ComputeL2Error(lowOrderCoeff); + + if (rank0_) cout << "|| interpolation error ||_2 = " << err << endl; + + // clean up aux data + delete aux_pfunc_; + delete[] aux_U_data; + delete aux_pfes_; + delete aux_fec_; + } else { + if (rank0_) cout << "Skipping family " << group_ << " doesn't allow changing order." << endl; + } +} + // register a new IO family which maps to a ParGridFunction void IODataOrganizer::registerIOFamily(std::string description, std::string group, ParGridFunction *pfunc, bool auxRestart, bool _inRestartFile, FiniteElementCollection *fec) { IOFamily family(description, group, pfunc); - std::vector vars; family.allowsAuxRestart = auxRestart; family.inRestartFile = _inRestartFile; family.fec_ = fec; @@ -595,17 +737,14 @@ void IODataOrganizer::registerIOFamily(std::string description, std::string grou } families_.push_back(family); - vars_[group] = vars; - - return; } // register individual variables for IO family void IODataOrganizer::registerIOVar(std::string group, std::string varName, int index, bool inRestartFile) { IOVar newvar{varName, index, inRestartFile}; - vars_[group].push_back(newvar); - - return; + const int f = getIOFamilyIndex(group); + assert(f >= 0); + families_[f].vars_.push_back(newvar); } // return the index to IO family given a group name @@ -664,9 +803,6 @@ void IODataOrganizer::write(hid_t file, bool serial) { for (auto fam : families_) { if (require_serialization) fam.serializeForWrite(); - // get defined variables for this IO family - vector vars = vars_[fam.group_]; - // determine if rank 0 const int rank = fam.pfunc_->ParFESpace()->GetMyRank(); const bool rank0 = (rank == 0); @@ -696,7 +832,7 @@ void IODataOrganizer::write(hid_t file, bool serial) { const double *data = fam.serial_sol->HostRead(); // save raw data - for (auto var : vars) write_soln_data(group, var.varName_, dataspace, data + var.index_ * dims[0], rank0); + for (auto var : fam.vars_) write_soln_data(group, var.varName_, dataspace, data + var.index_ * dims[0], rank0); } } else { // otherwise, use all ranks to write data from fam.pfunc_ @@ -710,7 +846,7 @@ void IODataOrganizer::write(hid_t file, bool serial) { const double *data = fam.pfunc_->HostRead(); // save raw data - for (auto var : vars) write_soln_data(group, var.varName_, dataspace, data + var.index_ * dims[0], rank0); + for (auto var : fam.vars_) write_soln_data(group, var.varName_, dataspace, data + var.index_ * dims[0], rank0); } if (group >= 0) H5Gclose(group); if (dataspace >= 0) H5Sclose(dataspace); @@ -733,16 +869,13 @@ void IODataOrganizer::read(hid_t file, bool serial, int read_order) { } assert(nprocs_max > 0); - hsize_t numInSoln = 0; - hid_t data_soln; - // Loop over defined IO families to load desired input for (auto fam : families_) { if (fam.inRestartFile) { // read mode const int rank = fam.pfunc_->ParFESpace()->GetMyRank(); const bool rank0 = (rank == 0); - if (rank0) cout << "Reading in solutiond data from restart..." << endl; + if (rank0) cout << "Reading in solution data from restart..." << endl; int fam_order; bool change_order = false; @@ -757,139 +890,12 @@ void IODataOrganizer::read(hid_t file, bool serial, int read_order) { if (rank0) grvy_printf(gerror, "[ERROR]: Serial read does not support order change.\n"); exit(ERROR); } - - if (fam.allowsAuxRestart) { // read mode - double *aux_U_data = NULL; - - // verify Dofs match expectations with current mesh - hid_t dataspace; - - vector vars = vars_[fam.group_]; - string varGroupName = fam.group_; - varGroupName.append("/"); - varGroupName.append(vars[0].varName_); - - data_soln = H5Dopen2(file, varGroupName.c_str(), H5P_DEFAULT); - assert(data_soln >= 0); - dataspace = H5Dget_space(data_soln); - numInSoln = H5Sget_simple_extent_npoints(dataspace); - H5Dclose(data_soln); - - assert(fam.fec_ != NULL); - fam.aux_fec_ = fam.fec_->Clone(read_order); - - int nvars = vars_[fam.group_].size(); - fam.aux_pfes_ = - new ParFiniteElementSpace(fam.pfunc_->ParFESpace()->GetParMesh(), fam.aux_fec_, nvars, Ordering::byNODES); - - aux_U_data = new double[nvars * fam.aux_pfes_->GetNDofs()]; - fam.aux_pfunc_ = new ParGridFunction(fam.aux_pfes_, aux_U_data); - - int dof = fam.aux_pfes_->GetNDofs(); - - assert((int)numInSoln == dof); - - // get pointer to raw data - double *data = fam.aux_pfunc_->HostReadWrite(); - - for (auto var : vars) { - if (var.inRestartFile_) { - std::string h5Path = fam.group_ + "/" + var.varName_; - if (rank0) grvy_printf(ginfo, "--> Reading h5 path = %s\n", h5Path.c_str()); - read_partitioned_soln_data(file, h5Path.c_str(), var.index_ * numInSoln, data); - } - } - - if (rank0) cout << "Interpolating from read_order = " << read_order << " to solution space" << endl; - - // Interpolate from aux to new order - fam.pfunc_->ProjectGridFunction(*fam.aux_pfunc_); - - // If interpolation was successful, the L2 norm of the - // difference between the auxOrder and order versions should be - // zero (well... to within round-off-induced error) - VectorGridFunctionCoefficient lowOrderCoeff(fam.aux_pfunc_); - double err = fam.pfunc_->ComputeL2Error(lowOrderCoeff); - - if (rank0) cout << "|| interpolation error ||_2 = " << err << endl; - - // clean up aux data - delete fam.aux_pfunc_; - delete[] aux_U_data; - delete fam.aux_pfes_; - delete fam.aux_fec_; - } else { - if (rank0) cout << "Skipping family " << fam.group_ << " doesn't allow changing order." << endl; - } + fam.readChangeOrder(file, read_order); } else { if (serial && nprocs_max > 1) { - const Array &partitioning = *(fam.partitioning_); - MPI_Comm comm = fam.pfunc_->ParFESpace()->GetComm(); - - // verify Dofs match expectations with current mesh - if (rank0) { - hid_t dataspace; - - vector vars = vars_[fam.group_]; - string varGroupName = fam.group_; - varGroupName.append("/"); - varGroupName.append(vars[0].varName_); - - data_soln = H5Dopen2(file, varGroupName.c_str(), H5P_DEFAULT); - assert(data_soln >= 0); - dataspace = H5Dget_space(data_soln); - numInSoln = H5Sget_simple_extent_npoints(dataspace); - H5Dclose(data_soln); - } - - int dof = fam.pfunc_->ParFESpace()->GetNDofs(); - - int dof_global; - MPI_Reduce(&dof, &dof_global, 1, MPI_INT, MPI_SUM, 0, comm); - dof = dof_global; - - if (rank0) assert((int)numInSoln == dof); - - // get pointer to raw data - double *data = fam.pfunc_->HostReadWrite(); - - vector vars = vars_[fam.group_]; - for (auto var : vars) { - if (var.inRestartFile_) { - std::string h5Path = fam.group_ + "/" + var.varName_; - if (rank0) grvy_printf(ginfo, "--> Reading h5 path = %s\n", h5Path.c_str()); - read_serialized_soln_data(file, h5Path.c_str(), dof, var.index_, data, fam, partitioning); - } - } + fam.readSerial(file); } else { - // verify Dofs match expectations with current mesh - hid_t dataspace; - - vector vars = vars_[fam.group_]; - string varGroupName = fam.group_; - varGroupName.append("/"); - varGroupName.append(vars[0].varName_); - - data_soln = H5Dopen2(file, varGroupName.c_str(), H5P_DEFAULT); - assert(data_soln >= 0); - dataspace = H5Dget_space(data_soln); - numInSoln = H5Sget_simple_extent_npoints(dataspace); - H5Dclose(data_soln); - - int dof = fam.pfunc_->ParFESpace()->GetNDofs(); - assert((int)numInSoln == dof); - - // get pointer to raw data - double *data = fam.pfunc_->HostReadWrite(); - - // read from file into appropriate spot in data - for (auto var : vars) { - if (var.inRestartFile_) { - std::string h5Path = fam.group_ + "/" + var.varName_; - if (rank0) grvy_printf(ginfo, "--> Reading h5 path = %s\n", h5Path.c_str()); - read_partitioned_soln_data(file, h5Path.c_str(), var.index_ * numInSoln, data); - } - } + fam.readPartitioned(file); } } } diff --git a/src/io.hpp b/src/io.hpp index af8ebce0f..fdef4dd8f 100644 --- a/src/io.hpp +++ b/src/io.hpp @@ -41,11 +41,22 @@ #include "run_configuration.hpp" #include "tps_mfem_wrap.hpp" +class IOVar { + public: + std::string varName_; // solution variable + int index_; // variable index in the pargrid function + bool inRestartFile_; // Check if we want to read this variable +}; + class IOFamily { + protected: + bool rank0_; + public: std::string description_; // family description std::string group_; // HDF5 group name mfem::ParGridFunction *pfunc_; // pargrid function owning the data for this IO family + std::vector vars_; // solution var info for this family bool allowsAuxRestart; @@ -67,17 +78,12 @@ class IOFamily { mfem::ParFiniteElementSpace *aux_pfes_ = nullptr; mfem::ParGridFunction *aux_pfunc_ = nullptr; - IOFamily(std::string desc, std::string grp, mfem::ParGridFunction *pf) - : description_(desc), group_(grp), pfunc_(pf) {} + IOFamily(std::string desc, std::string grp, mfem::ParGridFunction *pf); void serializeForWrite(); -}; - -class IOVar { - public: - std::string varName_; // solution variable - int index_; // variable index in the pargrid function - bool inRestartFile_; // Check if we want to read this variable + void readPartitioned(hid_t file); + void readSerial(hid_t file); + void readChangeOrder(hid_t file, int read_order); }; class IODataOrganizer { @@ -85,8 +91,7 @@ class IODataOrganizer { bool supports_serial_ = false; public: - std::vector families_; // registered IO families - std::map> vars_; // solution var info for each IO family + std::vector families_; // registered IO families void registerIOFamily(std::string description, std::string group, mfem::ParGridFunction *pfunc, bool auxRestart = true, bool inRestartFile = true, mfem::FiniteElementCollection *fec = NULL); From d1a5fe5f1a535db5e42798a6f0c14ae0af1f93f2 Mon Sep 17 00:00:00 2001 From: "Todd A. Oliver" Date: Thu, 1 Feb 2024 08:39:55 -0600 Subject: [PATCH 10/14] For readability, move bulk of write code into IOFamily (#244) Now the methods IOFamily::writePartitioned and IOFamily::writeSerial handle the data write, making the logic in IODataOrganizer::write easier to follow. --- src/io.cpp | 94 ++++++++++++++++++++++++++++++++---------------------- src/io.hpp | 5 ++- 2 files changed, 60 insertions(+), 39 deletions(-) diff --git a/src/io.cpp b/src/io.cpp index 53696f237..e7133b4fc 100644 --- a/src/io.cpp +++ b/src/io.cpp @@ -585,6 +585,58 @@ void IOFamily::serializeForWrite() { } } +void IOFamily::writePartitioned(hid_t file) { + hsize_t dims[1]; + hid_t group = -1; + hid_t dataspace = -1; + + // use all ranks to write data from this->.pfunc_ + dims[0] = pfunc_->ParFESpace()->GetNDofs(); + dataspace = H5Screate_simple(1, dims, NULL); + assert(dataspace >= 0); + group = H5Gcreate(file, group_.c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + assert(group >= 0); + + // get pointer to raw data + const double *data = pfunc_->HostRead(); + + // save raw data + for (auto var : vars_) write_soln_data(group, var.varName_, dataspace, data + var.index_ * dims[0], rank0_); + + if (group >= 0) H5Gclose(group); + if (dataspace >= 0) H5Sclose(dataspace); +} + +void IOFamily::writeSerial(hid_t file) { + // Only get here if need to serialize + serializeForWrite(); + + // Only rank 0 writes data (which has just been collected into this->serial_sol) + if (rank0_) { + hsize_t dims[1]; + hid_t group = -1; + hid_t dataspace = -1; + + assert(serial_fes != NULL); + dims[0] = serial_fes->GetNDofs(); + + dataspace = H5Screate_simple(1, dims, NULL); + assert(dataspace >= 0); + group = H5Gcreate(file, group_.c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + assert(group >= 0); + + // get pointer to raw data + assert(serial_sol != NULL); + const double *data = serial_sol->HostRead(); + + // save raw data + for (auto var : vars_) write_soln_data(group, var.varName_, dataspace, data + var.index_ * dims[0], rank0_); + + if (group >= 0) H5Gclose(group); + if (dataspace >= 0) H5Sclose(dataspace); + } +} + void IOFamily::readPartitioned(hid_t file) { hid_t dataspace; hsize_t numInSoln = 0; @@ -801,55 +853,20 @@ void IODataOrganizer::write(hid_t file, bool serial) { // Loop over defined IO families to save desired output for (auto fam : families_) { - if (require_serialization) fam.serializeForWrite(); - - // determine if rank 0 const int rank = fam.pfunc_->ParFESpace()->GetMyRank(); const bool rank0 = (rank == 0); - hsize_t dims[1]; - hid_t group = -1; - hid_t dataspace = -1; - if (rank0) { grvy_printf(ginfo, "\nCreating HDF5 group for defined IO families\n"); grvy_printf(ginfo, "--> %s : %s\n", fam.group_.c_str(), fam.description_.c_str()); } + // Write handled by appropriate method from IOFamily if (require_serialization) { - // if require_serialization, then we just populated fam.serial_sol above, and now we use only rank0 to write - if (rank0) { - assert(fam.serial_fes != NULL); - dims[0] = fam.serial_fes->GetNDofs(); - - dataspace = H5Screate_simple(1, dims, NULL); - assert(dataspace >= 0); - group = H5Gcreate(file, fam.group_.c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); - assert(group >= 0); - - // get pointer to raw data - assert(fam.serial_sol != NULL); - const double *data = fam.serial_sol->HostRead(); - - // save raw data - for (auto var : fam.vars_) write_soln_data(group, var.varName_, dataspace, data + var.index_ * dims[0], rank0); - } + fam.writeSerial(file); } else { - // otherwise, use all ranks to write data from fam.pfunc_ - dims[0] = fam.pfunc_->ParFESpace()->GetNDofs(); - dataspace = H5Screate_simple(1, dims, NULL); - assert(dataspace >= 0); - group = H5Gcreate(file, fam.group_.c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); - assert(group >= 0); - - // get pointer to raw data - const double *data = fam.pfunc_->HostRead(); - - // save raw data - for (auto var : fam.vars_) write_soln_data(group, var.varName_, dataspace, data + var.index_ * dims[0], rank0); + fam.writePartitioned(file); } - if (group >= 0) H5Gclose(group); - if (dataspace >= 0) H5Sclose(dataspace); } } @@ -885,6 +902,7 @@ void IODataOrganizer::read(hid_t file, bool serial, int read_order) { change_order = (fam_order != read_order); } + // Read handled by appropriate method from IOFamily if (change_order) { if (serial) { if (rank0) grvy_printf(gerror, "[ERROR]: Serial read does not support order change.\n"); diff --git a/src/io.hpp b/src/io.hpp index fdef4dd8f..5a8e7a77c 100644 --- a/src/io.hpp +++ b/src/io.hpp @@ -52,6 +52,8 @@ class IOFamily { protected: bool rank0_; + void serializeForWrite(); + public: std::string description_; // family description std::string group_; // HDF5 group name @@ -80,7 +82,8 @@ class IOFamily { IOFamily(std::string desc, std::string grp, mfem::ParGridFunction *pf); - void serializeForWrite(); + void writePartitioned(hid_t file); + void writeSerial(hid_t file); void readPartitioned(hid_t file); void readSerial(hid_t file); void readChangeOrder(hid_t file, int read_order); From caaba320882e11564a9a65daaa25e234f144c7c3 Mon Sep 17 00:00:00 2001 From: "Todd A. Oliver" Date: Thu, 1 Feb 2024 09:22:46 -0600 Subject: [PATCH 11/14] Move read_serialized_soln_data functionality into IOFamily (#244) And rename to readDistributeSerializedVariable. This function needed the IOFamily info, so makes more sense as a class method. --- src/io.cpp | 57 ++++++++++++++++++++++++++---------------------------- src/io.hpp | 5 ++--- 2 files changed, 29 insertions(+), 33 deletions(-) diff --git a/src/io.cpp b/src/io.cpp index e7133b4fc..cff0de1b1 100644 --- a/src/io.cpp +++ b/src/io.cpp @@ -373,37 +373,32 @@ void read_partitioned_soln_data(hid_t file, string varName, size_t index, double } // convenience function to read and distribute solution data for serialized restarts -void read_serialized_soln_data(hid_t file, string varName, int numDof, int varOffset, double *data, IOFamily &fam, - Array partitioning) { - MPI_Comm comm = fam.pfunc_->ParFESpace()->GetComm(); - const int rank = fam.pfunc_->ParFESpace()->GetMyRank(); - const bool rank0 = (rank == 0); - const int nprocs = fam.pfunc_->ParFESpace()->GetNRanks(); +void IOFamily::readDistributeSerializedVariable(hid_t file, string varName, int numDof, int varOffset, double *data) { + MPI_Comm comm = pfunc_->ParFESpace()->GetComm(); + const int myrank = pfunc_->ParFESpace()->GetMyRank(); + const int nprocs = pfunc_->ParFESpace()->GetNRanks(); // Ensure number of elements is consistent with partition array - const int local_ne = fam.pfunc_->ParFESpace()->GetParMesh()->GetNE(); + const int local_ne = pfunc_->ParFESpace()->GetParMesh()->GetNE(); int global_ne; MPI_Allreduce(&local_ne, &global_ne, 1, MPI_INT, MPI_SUM, comm); + + Array &partitioning = *partitioning_; assert(partitioning.Size() == global_ne); - hid_t data_soln; - herr_t status; int numStateVars; - numStateVars = fam.pfunc_->Size() / fam.pfunc_->ParFESpace()->GetNDofs(); + numStateVars = pfunc_->Size() / pfunc_->ParFESpace()->GetNDofs(); const int tag = 20; - if (rank0) { + if (rank0_) { grvy_printf(ginfo, "[RestartSerial]: Reading %s for distribution\n", varName.c_str()); Vector data_serial; data_serial.SetSize(numDof); - data_soln = H5Dopen2(file, varName.c_str(), H5P_DEFAULT); - assert(data_soln >= 0); - status = H5Dread(data_soln, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, data_serial.HostReadWrite()); - assert(status >= 0); - H5Dclose(data_soln); + + read_partitioned_soln_data(file, varName, 0, data_serial.HostWrite()); // assign solution owned by rank 0 Array lvdofs, gvdofs; @@ -414,14 +409,14 @@ void read_serialized_soln_data(hid_t file, string varName, int numDof, int varOf for (int gelem = 0; gelem < global_ne; gelem++) { if (partitioning[gelem] == 0) { // cull out subset of local vdof vector to use for this solution var - fam.pfunc_->ParFESpace()->GetElementVDofs(counter, lvdofs); + pfunc_->ParFESpace()->GetElementVDofs(counter, lvdofs); int numDof_per_this_elem = lvdofs.Size() / numStateVars; int ldof_start_index = varOffset * numDof_per_this_elem; // cull out global vdofs - [subtle note]: we only use the // indexing corresponding to the first state vector reference in // gvdofs() since data_serial() holds the solution for a single state vector - fam.serial_fes->GetElementVDofs(gelem, gvdofs); + serial_fes->GetElementVDofs(gelem, gvdofs); int gdof_start_index = 0; for (int i = 0; i < numDof_per_this_elem; i++) @@ -435,7 +430,7 @@ void read_serialized_soln_data(hid_t file, string varName, int numDof, int varOf std::vector packedData; for (int gelem = 0; gelem < global_ne; gelem++) if (partitioning[gelem] == rank) { - fam.serial_fes->GetElementVDofs(gelem, gvdofs); + serial_fes->GetElementVDofs(gelem, gvdofs); int numDof_per_this_elem = gvdofs.Size() / numStateVars; for (int i = 0; i < numDof_per_this_elem; i++) packedData.push_back(data_serial[gvdofs[i]]); } @@ -443,8 +438,9 @@ void read_serialized_soln_data(hid_t file, string varName, int numDof, int varOf MPI_Send(packedData.data(), packedData.size(), MPI_DOUBLE, rank, tag, comm); } } else { // <-- end rank 0 - int numlDofs = fam.pfunc_->ParFESpace()->GetNDofs(); - grvy_printf(gdebug, "[%i]: local number of state vars to receive = %i (var=%s)\n", rank, numlDofs, varName.c_str()); + int numlDofs = pfunc_->ParFESpace()->GetNDofs(); + grvy_printf(gdebug, "[%i]: local number of state vars to receive = %i (var=%s)\n", myrank, numlDofs, + varName.c_str()); std::vector packedData(numlDofs); @@ -457,21 +453,17 @@ void read_serialized_soln_data(hid_t file, string varName, int numDof, int varOf } // convenience function to write HDF5 data -void write_soln_data(hid_t group, string varName, hid_t dataspace, const double *data, bool rank0) { +void write_soln_data(hid_t group, string varName, hid_t dataspace, const double *data) { hid_t data_soln; herr_t status; assert(group >= 0); - if (rank0) grvy_printf(ginfo, " --> Saving (%s)\n", varName.c_str()); - data_soln = H5Dcreate2(group, varName.c_str(), H5T_NATIVE_DOUBLE, dataspace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); assert(data_soln >= 0); status = H5Dwrite(data_soln, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, data); assert(status >= 0); H5Dclose(data_soln); - - return; } void M2ulPhyS::readTable(const std::string &inputPath, TableInput &result) { @@ -601,7 +593,10 @@ void IOFamily::writePartitioned(hid_t file) { const double *data = pfunc_->HostRead(); // save raw data - for (auto var : vars_) write_soln_data(group, var.varName_, dataspace, data + var.index_ * dims[0], rank0_); + for (auto var : vars_) { + if (rank0_) grvy_printf(ginfo, " --> Saving (%s)\n", var.varName_.c_str()); + write_soln_data(group, var.varName_, dataspace, data + var.index_ * dims[0]); + } if (group >= 0) H5Gclose(group); if (dataspace >= 0) H5Sclose(dataspace); @@ -630,7 +625,10 @@ void IOFamily::writeSerial(hid_t file) { const double *data = serial_sol->HostRead(); // save raw data - for (auto var : vars_) write_soln_data(group, var.varName_, dataspace, data + var.index_ * dims[0], rank0_); + for (auto var : vars_) { + grvy_printf(ginfo, " --> Saving (%s)\n", var.varName_.c_str()); + write_soln_data(group, var.varName_, dataspace, data + var.index_ * dims[0]); + } if (group >= 0) H5Gclose(group); if (dataspace >= 0) H5Sclose(dataspace); @@ -671,7 +669,6 @@ void IOFamily::readPartitioned(hid_t file) { } void IOFamily::readSerial(hid_t file) { - const Array &partitioning = *partitioning_; MPI_Comm comm = pfunc_->ParFESpace()->GetComm(); hsize_t numInSoln = 0; @@ -707,7 +704,7 @@ void IOFamily::readSerial(hid_t file) { if (var.inRestartFile_) { std::string h5Path = group_ + "/" + var.varName_; if (rank0_) grvy_printf(ginfo, "--> Reading h5 path = %s\n", h5Path.c_str()); - read_serialized_soln_data(file, h5Path.c_str(), dof, var.index_, data, *this, partitioning); + readDistributeSerializedVariable(file, h5Path.c_str(), dof, var.index_, data); } } } diff --git a/src/io.hpp b/src/io.hpp index 5a8e7a77c..f5076885f 100644 --- a/src/io.hpp +++ b/src/io.hpp @@ -53,6 +53,7 @@ class IOFamily { bool rank0_; void serializeForWrite(); + void readDistributeSerializedVariable(hid_t file, string varName, int numDof, int varOffset, double *data); public: std::string description_; // family description @@ -110,9 +111,7 @@ class IODataOrganizer { }; void read_partitioned_soln_data(hid_t file, std::string varName, size_t index, double *data); -void read_serialized_soln_data(hid_t file, std::string varName, int numDof, int varOffset, double *data, IOFamily &fam, - mfem::Array partitioning); -void write_soln_data(hid_t group, std::string varName, hid_t dataspace, const double *data, bool rank0); +void write_soln_data(hid_t group, std::string varName, hid_t dataspace, const double *data); void partitioning_file_hdf5(std::string mode, const RunConfiguration &config, MPI_Groups *groupsMPI, int nelemGlobal, mfem::Array &partitioning); #endif // IO_HPP_ From 9ac2fa84b770fd23c34e7f61b98919eb77a850da Mon Sep 17 00:00:00 2001 From: "Todd A. Oliver" Date: Thu, 1 Feb 2024 11:54:06 -0600 Subject: [PATCH 12/14] Add get_variable_size_hdf5 (#244) To reduce duplicated code in size consistency checks. Also, some other minor cleaning up. --- src/io.cpp | 156 +++++++++++++++++++++-------------------------------- src/io.hpp | 19 +++++-- 2 files changed, 75 insertions(+), 100 deletions(-) diff --git a/src/io.cpp b/src/io.cpp index cff0de1b1..8791d89bc 100644 --- a/src/io.cpp +++ b/src/io.cpp @@ -360,6 +360,16 @@ void partitioning_file_hdf5(std::string mode, const RunConfiguration &config, MP grvy_timer_end(__func__); } +// convenience function to check size of variable in hdf5 file +hsize_t get_variable_size_hdf5(hid_t file, std::string name) { + hid_t data_soln = H5Dopen2(file, name.c_str(), H5P_DEFAULT); + assert(data_soln >= 0); + hid_t dataspace = H5Dget_space(data_soln); + hsize_t numInSoln = H5Sget_simple_extent_npoints(dataspace); + H5Dclose(data_soln); + return numInSoln; +} + // convenience function to read solution data for parallel restarts void read_partitioned_soln_data(hid_t file, string varName, size_t index, double *data) { hid_t data_soln; @@ -373,22 +383,21 @@ void read_partitioned_soln_data(hid_t file, string varName, size_t index, double } // convenience function to read and distribute solution data for serialized restarts -void IOFamily::readDistributeSerializedVariable(hid_t file, string varName, int numDof, int varOffset, double *data) { +void IOFamily::readDistributeSerializedVariable(hid_t file, const IOVar &var, int numDof, double *data) { + std::string varName = group_ + "/" + var.varName_; + if (rank0_) grvy_printf(ginfo, "--> Reading h5 path = %s\n", varName.c_str()); + + const int varOffset = var.index_; + MPI_Comm comm = pfunc_->ParFESpace()->GetComm(); const int myrank = pfunc_->ParFESpace()->GetMyRank(); const int nprocs = pfunc_->ParFESpace()->GetNRanks(); - // Ensure number of elements is consistent with partition array - const int local_ne = pfunc_->ParFESpace()->GetParMesh()->GetNE(); - int global_ne; - MPI_Allreduce(&local_ne, &global_ne, 1, MPI_INT, MPI_SUM, comm); - Array &partitioning = *partitioning_; - assert(partitioning.Size() == global_ne); - - int numStateVars; + assert(partitioning.Size() == global_ne_); - numStateVars = pfunc_->Size() / pfunc_->ParFESpace()->GetNDofs(); + const unsigned int numStateVars = pfunc_->Size() / pfunc_->ParFESpace()->GetNDofs(); + assert(numStateVars == vars_.size()); const int tag = 20; @@ -406,7 +415,7 @@ void IOFamily::readDistributeSerializedVariable(hid_t file, string varName, int int counter = 0; // int ndof_per_elem; - for (int gelem = 0; gelem < global_ne; gelem++) { + for (int gelem = 0; gelem < global_ne_; gelem++) { if (partitioning[gelem] == 0) { // cull out subset of local vdof vector to use for this solution var pfunc_->ParFESpace()->GetElementVDofs(counter, lvdofs); @@ -428,7 +437,7 @@ void IOFamily::readDistributeSerializedVariable(hid_t file, string varName, int // pack remaining data and send to other processors for (int rank = 1; rank < nprocs; rank++) { std::vector packedData; - for (int gelem = 0; gelem < global_ne; gelem++) + for (int gelem = 0; gelem < global_ne_; gelem++) if (partitioning[gelem] == rank) { serial_fes->GetElementVDofs(gelem, gvdofs); int numDof_per_this_elem = gvdofs.Size() / numStateVars; @@ -516,22 +525,25 @@ void M2ulPhyS::readTable(const std::string &inputPath, TableInput &result) { IOFamily::IOFamily(std::string desc, std::string grp, mfem::ParGridFunction *pf) : description_(desc), group_(grp), pfunc_(pf) { rank0_ = (pfunc_->ParFESpace()->GetMyRank() == 0); + + MPI_Comm comm = this->pfunc_->ParFESpace()->GetComm(); + + // Set local and global number of elements + local_ne_ = this->pfunc_->ParFESpace()->GetParMesh()->GetNE(); + MPI_Allreduce(&local_ne_, &global_ne_, 1, MPI_INT, MPI_SUM, comm); + + // Set local and global number of dofs (per scalar field) + local_ndofs_ = pfunc_->ParFESpace()->GetNDofs(); + MPI_Allreduce(&local_ndofs_, &global_ndofs_, 1, MPI_INT, MPI_SUM, comm); } void IOFamily::serializeForWrite() { MPI_Comm comm = this->pfunc_->ParFESpace()->GetComm(); - const int local_ne = this->pfunc_->ParFESpace()->GetParMesh()->GetNE(); - const Array &partitioning = *(this->partitioning_); - const int *locToGlobElem = this->local_to_global_elem_; + assert(partitioning.Size() == global_ne_); - // Ensure consistency - int global_ne; - MPI_Reduce(&local_ne, &global_ne, 1, MPI_INT, MPI_SUM, 0, comm); - if (rank0_) { - assert(partitioning.Size() == global_ne); - } + const int *locToGlobElem = this->local_to_global_elem_; assert(locToGlobElem != NULL); ParGridFunction *pfunc = this->pfunc_; @@ -540,7 +552,7 @@ void IOFamily::serializeForWrite() { // copy my own data Array lvdofs, gvdofs; Vector lsoln; - for (int elem = 0; elem < local_ne; elem++) { + for (int elem = 0; elem < local_ne_; elem++) { int gelem = locToGlobElem[elem]; this->pfunc_->ParFESpace()->GetElementVDofs(elem, lvdofs); pfunc->GetSubVector(lvdofs, lsoln); @@ -549,7 +561,7 @@ void IOFamily::serializeForWrite() { } // have rank 0 receive data from other tasks and copy its own - for (int gelem = 0; gelem < global_ne; gelem++) { + for (int gelem = 0; gelem < global_ne_; gelem++) { int from_rank = partitioning[gelem]; if (from_rank != 0) { this->serial_fes->GetElementVDofs(gelem, gvdofs); @@ -565,7 +577,7 @@ void IOFamily::serializeForWrite() { // have non-zero ranks send their data to rank 0 Array lvdofs; Vector lsoln; - for (int elem = 0; elem < local_ne; elem++) { + for (int elem = 0; elem < local_ne_; elem++) { int gelem = locToGlobElem[elem]; // assert(gelem > 0); this->pfunc_->ParFESpace()->GetElementVDofs(elem, lvdofs); @@ -636,24 +648,10 @@ void IOFamily::writeSerial(hid_t file) { } void IOFamily::readPartitioned(hid_t file) { - hid_t dataspace; - hsize_t numInSoln = 0; - hid_t data_soln; - - string varGroupName = group_; - varGroupName.append("/"); - varGroupName.append(vars_[0].varName_); - - data_soln = H5Dopen2(file, varGroupName.c_str(), H5P_DEFAULT); - assert(data_soln >= 0); - dataspace = H5Dget_space(data_soln); - numInSoln = H5Sget_simple_extent_npoints(dataspace); - H5Dclose(data_soln); - - int dof = pfunc_->ParFESpace()->GetNDofs(); - - // verify Dofs match expectations with current mesh - assert((int)numInSoln == dof); + // Ensure that size of read matches expectations + std::string varGroupName = group_ + "/" + vars_[0].varName_; + const hsize_t numInSoln = get_variable_size_hdf5(file, varGroupName); + assert((int)numInSoln == local_ndofs_); // get pointer to raw data double *data = pfunc_->HostWrite(); @@ -669,76 +667,41 @@ void IOFamily::readPartitioned(hid_t file) { } void IOFamily::readSerial(hid_t file) { - MPI_Comm comm = pfunc_->ParFESpace()->GetComm(); - - hsize_t numInSoln = 0; - hid_t data_soln; - - // verify Dofs match expectations with current mesh + // Ensure that size of read matches expectations if (rank0_) { - hid_t dataspace; - - string varGroupName = group_; - varGroupName.append("/"); - varGroupName.append(vars_[0].varName_); - - data_soln = H5Dopen2(file, varGroupName.c_str(), H5P_DEFAULT); - assert(data_soln >= 0); - dataspace = H5Dget_space(data_soln); - numInSoln = H5Sget_simple_extent_npoints(dataspace); - H5Dclose(data_soln); + std::string varGroupName = group_ + "/" + vars_[0].varName_; + const hsize_t numInSoln = get_variable_size_hdf5(file, varGroupName); + assert((int)numInSoln == global_ndofs_); } - int dof = pfunc_->ParFESpace()->GetNDofs(); - - int dof_global; - MPI_Reduce(&dof, &dof_global, 1, MPI_INT, MPI_SUM, 0, comm); - dof = dof_global; - - if (rank0_) assert((int)numInSoln == dof); - // get pointer to raw data double *data = pfunc_->HostWrite(); + // read on rank 0 and distribute into data for (auto var : vars_) { if (var.inRestartFile_) { - std::string h5Path = group_ + "/" + var.varName_; - if (rank0_) grvy_printf(ginfo, "--> Reading h5 path = %s\n", h5Path.c_str()); - readDistributeSerializedVariable(file, h5Path.c_str(), dof, var.index_, data); + readDistributeSerializedVariable(file, var, global_ndofs_, data); } } } void IOFamily::readChangeOrder(hid_t file, int read_order) { if (allowsAuxRestart) { - hid_t dataspace; - hsize_t numInSoln = 0; - hid_t data_soln; - - double *aux_U_data = NULL; - - string varGroupName = group_; - varGroupName.append("/"); - varGroupName.append(vars_[0].varName_); - - data_soln = H5Dopen2(file, varGroupName.c_str(), H5P_DEFAULT); - assert(data_soln >= 0); - dataspace = H5Dget_space(data_soln); - numInSoln = H5Sget_simple_extent_npoints(dataspace); - H5Dclose(data_soln); - + // Set up "auxilliary" FiniteElementCollection, which matches original, but with different order assert(fec_ != NULL); aux_fec_ = fec_->Clone(read_order); - int nvars = vars_.size(); + // Set up "auxilliary" ParGridFunction, to read data into before order change + const int nvars = vars_.size(); aux_pfes_ = new ParFiniteElementSpace(pfunc_->ParFESpace()->GetParMesh(), aux_fec_, nvars, Ordering::byNODES); - - aux_U_data = new double[nvars * aux_pfes_->GetNDofs()]; + double *aux_U_data = new double[nvars * aux_pfes_->GetNDofs()]; aux_pfunc_ = new ParGridFunction(aux_pfes_, aux_U_data); + const int aux_dof = aux_pfes_->GetNDofs(); - int dof = aux_pfes_->GetNDofs(); - - assert((int)numInSoln == dof); + // Ensure size of read matches expectations + std::string varGroupName = group_ + "/" + vars_[0].varName_; + const hsize_t numInSoln = get_variable_size_hdf5(file, varGroupName); + assert((int)numInSoln == aux_dof); // get pointer to raw data double *data = aux_pfunc_->HostWrite(); @@ -747,15 +710,15 @@ void IOFamily::readChangeOrder(hid_t file, int read_order) { if (var.inRestartFile_) { std::string h5Path = group_ + "/" + var.varName_; if (rank0_) grvy_printf(ginfo, "--> Reading h5 path = %s\n", h5Path.c_str()); - read_partitioned_soln_data(file, h5Path.c_str(), var.index_ * numInSoln, data); + read_partitioned_soln_data(file, h5Path.c_str(), var.index_ * aux_dof, data); } } + // Interpolate from "auxilliary" to new order if (rank0_) cout << "Interpolating from read_order = " << read_order << " to solution space" << endl; - - // Interpolate from aux to new order pfunc_->ProjectGridFunction(*aux_pfunc_); +#if 0 // If interpolation was successful, the L2 norm of the // difference between the auxOrder and order versions should be // zero (well... to within round-off-induced error) @@ -763,6 +726,7 @@ void IOFamily::readChangeOrder(hid_t file, int read_order) { double err = pfunc_->ComputeL2Error(lowOrderCoeff); if (rank0_) cout << "|| interpolation error ||_2 = " << err << endl; +#endif // clean up aux data delete aux_pfunc_; @@ -770,7 +734,7 @@ void IOFamily::readChangeOrder(hid_t file, int read_order) { delete aux_pfes_; delete aux_fec_; } else { - if (rank0_) cout << "Skipping family " << group_ << " doesn't allow changing order." << endl; + if (rank0_) cout << "Skipping family " << group_ << " because it doesn't allow changing order." << endl; } } diff --git a/src/io.hpp b/src/io.hpp index f5076885f..fc414f6ae 100644 --- a/src/io.hpp +++ b/src/io.hpp @@ -41,6 +41,8 @@ #include "run_configuration.hpp" #include "tps_mfem_wrap.hpp" +class IODataOrganizer; + class IOVar { public: std::string varName_; // solution variable @@ -49,13 +51,17 @@ class IOVar { }; class IOFamily { + friend class IODataOrganizer; + protected: - bool rank0_; + bool rank0_ = false; - void serializeForWrite(); - void readDistributeSerializedVariable(hid_t file, string varName, int numDof, int varOffset, double *data); + int local_ne_ = -1; + int global_ne_ = -1; + + int local_ndofs_ = -1; + int global_ndofs_ = -1; - public: std::string description_; // family description std::string group_; // HDF5 group name mfem::ParGridFunction *pfunc_; // pargrid function owning the data for this IO family @@ -81,6 +87,10 @@ class IOFamily { mfem::ParFiniteElementSpace *aux_pfes_ = nullptr; mfem::ParGridFunction *aux_pfunc_ = nullptr; + void serializeForWrite(); + void readDistributeSerializedVariable(hid_t file, const IOVar &var, int numDof, double *data); + + public: IOFamily(std::string desc, std::string grp, mfem::ParGridFunction *pf); void writePartitioned(hid_t file); @@ -110,6 +120,7 @@ class IODataOrganizer { void read(hid_t file, bool serial, int read_order = -1); }; +hsize_t get_variable_size_hdf5(hid_t file, std::string name); void read_partitioned_soln_data(hid_t file, std::string varName, size_t index, double *data); void write_soln_data(hid_t group, std::string varName, hid_t dataspace, const double *data); void partitioning_file_hdf5(std::string mode, const RunConfiguration &config, MPI_Groups *groupsMPI, int nelemGlobal, From fbfce11e0b01f407ede76fe5e825bddc53979b67 Mon Sep 17 00:00:00 2001 From: "Todd A. Oliver" Date: Thu, 1 Feb 2024 14:27:36 -0600 Subject: [PATCH 13/14] Add doxygen comments for IO classes (#244) Plus a couple name changes to be a tiny bit more descriptive --- src/io.cpp | 19 ++-- src/io.hpp | 251 ++++++++++++++++++++++++++++++++++++++++++++++++----- 2 files changed, 242 insertions(+), 28 deletions(-) diff --git a/src/io.cpp b/src/io.cpp index 8791d89bc..37b2cc085 100644 --- a/src/io.cpp +++ b/src/io.cpp @@ -29,6 +29,9 @@ // OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. // -----------------------------------------------------------------------------------el- +/** @file + * @copydoc io.hpp + */ #include "io.hpp" #include @@ -371,7 +374,7 @@ hsize_t get_variable_size_hdf5(hid_t file, std::string name) { } // convenience function to read solution data for parallel restarts -void read_partitioned_soln_data(hid_t file, string varName, size_t index, double *data) { +void read_variable_data_hdf5(hid_t file, string varName, size_t index, double *data) { hid_t data_soln; herr_t status; @@ -407,7 +410,7 @@ void IOFamily::readDistributeSerializedVariable(hid_t file, const IOVar &var, in Vector data_serial; data_serial.SetSize(numDof); - read_partitioned_soln_data(file, varName, 0, data_serial.HostWrite()); + read_variable_data_hdf5(file, varName, 0, data_serial.HostWrite()); // assign solution owned by rank 0 Array lvdofs, gvdofs; @@ -462,7 +465,7 @@ void IOFamily::readDistributeSerializedVariable(hid_t file, const IOVar &var, in } // convenience function to write HDF5 data -void write_soln_data(hid_t group, string varName, hid_t dataspace, const double *data) { +void write_variable_data_hdf5(hid_t group, string varName, hid_t dataspace, const double *data) { hid_t data_soln; herr_t status; assert(group >= 0); @@ -607,7 +610,7 @@ void IOFamily::writePartitioned(hid_t file) { // save raw data for (auto var : vars_) { if (rank0_) grvy_printf(ginfo, " --> Saving (%s)\n", var.varName_.c_str()); - write_soln_data(group, var.varName_, dataspace, data + var.index_ * dims[0]); + write_variable_data_hdf5(group, var.varName_, dataspace, data + var.index_ * dims[0]); } if (group >= 0) H5Gclose(group); @@ -639,7 +642,7 @@ void IOFamily::writeSerial(hid_t file) { // save raw data for (auto var : vars_) { grvy_printf(ginfo, " --> Saving (%s)\n", var.varName_.c_str()); - write_soln_data(group, var.varName_, dataspace, data + var.index_ * dims[0]); + write_variable_data_hdf5(group, var.varName_, dataspace, data + var.index_ * dims[0]); } if (group >= 0) H5Gclose(group); @@ -661,7 +664,7 @@ void IOFamily::readPartitioned(hid_t file) { if (var.inRestartFile_) { std::string h5Path = group_ + "/" + var.varName_; if (rank0_) grvy_printf(ginfo, "--> Reading h5 path = %s\n", h5Path.c_str()); - read_partitioned_soln_data(file, h5Path.c_str(), var.index_ * numInSoln, data); + read_variable_data_hdf5(file, h5Path.c_str(), var.index_ * numInSoln, data); } } } @@ -710,7 +713,7 @@ void IOFamily::readChangeOrder(hid_t file, int read_order) { if (var.inRestartFile_) { std::string h5Path = group_ + "/" + var.varName_; if (rank0_) grvy_printf(ginfo, "--> Reading h5 path = %s\n", h5Path.c_str()); - read_partitioned_soln_data(file, h5Path.c_str(), var.index_ * aux_dof, data); + read_variable_data_hdf5(file, h5Path.c_str(), var.index_ * aux_dof, data); } } @@ -761,7 +764,7 @@ void IODataOrganizer::registerIOVar(std::string group, std::string varName, int } // return the index to IO family given a group name -int IODataOrganizer::getIOFamilyIndex(std::string group) { +int IODataOrganizer::getIOFamilyIndex(std::string group) const { for (size_t i = 0; i < families_.size(); i++) if (families_[i].group_ == group) return (i); diff --git a/src/io.hpp b/src/io.hpp index fc414f6ae..556f5cd30 100644 --- a/src/io.hpp +++ b/src/io.hpp @@ -31,6 +31,11 @@ // -----------------------------------------------------------------------------------el- #ifndef IO_HPP_ #define IO_HPP_ +/** @file + * @brief Contains utilities for abstracting field IO tasks + * + * The primary class is IODataOrganizer. For example usage, see M2ulPhyS. + */ #include @@ -43,86 +48,292 @@ class IODataOrganizer; +/** + * @brief Stores info about variables to read/write. + * + * This info is used in IOFamily and IODataOrganizer and is not + * intended to be needed separately. + */ class IOVar { public: - std::string varName_; // solution variable - int index_; // variable index in the pargrid function - bool inRestartFile_; // Check if we want to read this variable + /** Variable name, which combines with family name to set path in hdf5 file */ + std::string varName_; + + /** Variable index within variables stored in corresponding ParGridFunction */ + int index_; // variable index in the pargrid function + + /** Whether variable is to be read + * + * This allows variables in the state to be skipped, e.g., to enable + * restarts from models that do not have all the same state variables. + */ + bool inRestartFile_; }; +/** + * @brief Stores info about a "family" of variables to read/write. + * + * Here a "family" refers to data stored in a single + * mfem::ParGridFunction object. Thus, the family has a single mesh + * and parallel decomposition. For the moment, the underlying mesh + * and decomposition is assumed fixed after initialization. + * + * This class is intended to be used through IODataOrganizer, which + * can track multiple IOFamily objects. + */ class IOFamily { friend class IODataOrganizer; protected: + /** Indicates whether current mpi rank is rank 0*/ bool rank0_ = false; + /** Local number of elements in this mpi rank */ int local_ne_ = -1; + + /** Global number of elements in whole mesh */ int global_ne_ = -1; + /** Local number of degrees of freedom on this mpi rank, as returned by ParGridFunction::GetNDofs */ int local_ndofs_ = -1; + + /** Global number of degrees of freedom */ int global_ndofs_ = -1; - std::string description_; // family description - std::string group_; // HDF5 group name - mfem::ParGridFunction *pfunc_; // pargrid function owning the data for this IO family - std::vector vars_; // solution var info for this family + /** Description string */ + std::string description_; + + /** Group name, used to set corresponding group in hdf5 restart files */ + std::string group_; + + /** mfem::ParGridFunction owning the data to be written or read into */ + mfem::ParGridFunction *pfunc_ = nullptr; + /** Variables that make up the family */ + std::vector vars_; + + /** Whether this family allows initialization (read) from other order elements */ bool allowsAuxRestart; - // true if the family is to be found in restart files + /** true if the family is to be found in restart files */ bool inRestartFile; - // space and grid function used for serial read/write + /** mfem::FiniteElementSpace for serial mesh (used if serial read and/or write requested) */ mfem::FiniteElementSpace *serial_fes = nullptr; + + /** mfem::GridFunction on serial mesh (used if serial read and/or write requested) */ mfem::GridFunction *serial_sol = nullptr; - // additional data used for serial read/write (needed to properly - // serialize fields in write or properly distribute data in read) - int *local_to_global_elem_ = nullptr; // maps from local elem index to global - mfem::Array *partitioning_ = nullptr; // partitioning[i] mpi rank for element i (in global numbering) + /** Map from local to global element numbering (used for serial read/write) */ + int *local_to_global_elem_ = nullptr; + + /** Partition array; partitioning[i] = mpi rank that owns the ith element (global numbering) */ + mfem::Array *partitioning_ = nullptr; - // Variables used for "auxilliary" read (i.e., restart from different order soln) - mfem::FiniteElementCollection *fec_ = nullptr; // collection used to instantiate pfunc_ + /** mfem::FiniteElementCollection used by pfunc_ (used for variable order read) */ + mfem::FiniteElementCollection *fec_ = nullptr; + + /** mfem::FiniteElementCollection with read order (used for variable order read) */ mfem::FiniteElementCollection *aux_fec_ = nullptr; + + /** mfem::FiniteElementCollection with read order elements (used for variable order read) */ mfem::ParFiniteElementSpace *aux_pfes_ = nullptr; + + /** mfem::ParGridFunction with read order elements (used for variable order read) */ mfem::ParGridFunction *aux_pfunc_ = nullptr; + /** + * @brief Prepare for serial write by collecting data onto rank 0 + * + * Each rank sends its data in pfunc_ to rank 0 and stores the + * results in serial_sol + */ void serializeForWrite(); + + /** + * @brief Read and distribute a single variable in a serialized restart file + * + * Rank 0 reads the data and sends to appropriate rank and location in pfunc_ + * @param file Open HDF5 file handle (e.g., as returned by H5Fcreate or H5Fopen) + * @param var IOVar for the variable being read + */ void readDistributeSerializedVariable(hid_t file, const IOVar &var, int numDof, double *data); public: + /** + * @brief Constructor, must provide description, group name, and ParGridFunction + * + * @param desc Description string + * @param grp Group name + * @param pf Pointer to mfem::ParGridFunction that owns the data + */ IOFamily(std::string desc, std::string grp, mfem::ParGridFunction *pf); + /** @brief "Partitioned" write (each mpi rank writes its local part of the ParGridFunction to a separate file) */ void writePartitioned(hid_t file); + + /** @brief "Serial" write (all data communicated to rank 0, which writes everything to a single file) */ void writeSerial(hid_t file); + + /** @brief "Partitioned" read (inverse of writePartitioned) */ void readPartitioned(hid_t file); + + /** @brief "Serial" read (inverse of readPartitioned) */ void readSerial(hid_t file); + + /** + * @brief "Partitioned" read of a different order solution + * + * which is then interpolated onto the desired space, as defined by pfunc_ + * @param file Open HDF5 file handle (e.g., as returned by H5Fcreate or H5Fopen) + * @param read_order Polynomial order of function in HDF5 file + */ void readChangeOrder(hid_t file, int read_order); }; +/** + * @brief Provides IO interface for multiple "families" + * + * This class is used to read & write the solution data contained in + * the mfem::ParGridFunction object in each IOFamily that is + * "registered". Once families and variables are registered, the + * underlying data can be read or written with a single call. Since + * it is expected that the calling class will have some class-specific + * metadata to read or write in addition to the solution data, this + * class does not handle opening or creating the hdf5 file or reading + * or writing solution metadata. These tasks must be handled outside + * the IODataOrganizer class. For example usage, see + * M2ulPhyS::initSolutionAndVisualizationVectors and + * M2ulPhyS::write_restart_files_hdf5. + */ class IODataOrganizer { protected: + /** Whether serial read/write is supported. To enable serial support, must call initializeSerial */ bool supports_serial_ = false; - public: - std::vector families_; // registered IO families + /** Registered IO families */ + std::vector families_; + /** Get the index of the IOFamily with group name group */ + int getIOFamilyIndex(std::string group) const; + + public: + /** + * @brief Register an IOFamily + * + * Constructs the family using the input arguments. + * + * @param description Description string for the family + * @param group Group name for the family + * @param pfunc mfem::ParGridFunction for the family + * @param auxRestart Allow initialization from different order read (optional, defaults to true) + * @param inRestartFile If false, skip family on read (optional, defaults to true) + * @param fec FiniteElementCollection for grid function (optional, only used for different order read) + */ void registerIOFamily(std::string description, std::string group, mfem::ParGridFunction *pfunc, bool auxRestart = true, bool inRestartFile = true, mfem::FiniteElementCollection *fec = NULL); + /** Destructor */ ~IODataOrganizer(); + /** + * @brief Register an IOVar as part of the corresponding family + * + * @param group Group name for the family to add the variable to + * @param varName Name of the variable + * @param index Variable index within the ParGridFunction of the family + * @param inRestartFile If false, skip variable on read (optional, defaults to true) + */ void registerIOVar(std::string group, std::string varName, int index, bool inRestartFile = true); - int getIOFamilyIndex(std::string group); + /** + * @brief Initialize data supporting serial read/write. + * + * @param root True for mpi rank = 0 + * @param serial True if we want seral read and/or write + * @param serial_mesh Pointer to serial mesh object + * @param locToGlob Map from local element index to global element index + * @param part Map from global element index to mpi rank owning that element + */ void initializeSerial(bool root, bool serial, mfem::Mesh *serial_mesh, int *locToGlob, mfem::Array *part); + /** + * @brief Write data from all families + * + * All families written to the same HDF5 file, which must be already + * open. If serial write requested, a data is written to a single + * HDF5 file by rank 0. Otherwise, all ranks write to separate + * files. + * + * @param file HDF5 file handle (e.g., as returned by H5Fcreate or H5Fopen) + * @param serial True for serial write + */ void write(hid_t file, bool serial); + + /** + * @brief Read data for all families + * + * All families are read from the same HDF5 file, which must be + * already open. If serial read is requested, the data is read from + * a single HDF5 file by rank 0. Otherwise, all ranks read from + * separate files. For variable order (i.e., data in file had + * different order than desired), the order of the file being read + * must be supplied. + * + * @param file HDF5 file handle (e.g., as returned by H5Fopen) + * @param serial True for serial read + * @param read_order Polynomial order for data in file (optional, only required for variable order read) + */ void read(hid_t file, bool serial, int read_order = -1); }; +/** + * @brief Get number of points to be read + * + * This convenience function returns the number of points to be read + * for a given name. This is used frequently to check for consistency + * between the size of the read and the size of the container being + * read into. + * + * @param file Open HDF5 file + * @param name Name to query (corresponds to "family group/variable name") + * @returns Size of corresponding data in file + */ hsize_t get_variable_size_hdf5(hid_t file, std::string name); -void read_partitioned_soln_data(hid_t file, std::string varName, size_t index, double *data); -void write_soln_data(hid_t group, std::string varName, hid_t dataspace, const double *data); + +/** + * @brief Read data from hdf5 file + * + * Reads data from the hdf5 file into the data buffer, starting at + * location index (i.e., fill starting at data[index]) + * + * @param file Open HDF5 file + * @param varName Name of object to read (corresponds to "family group/variable name") + * @param index Starting location within data buffer + * @param data Data buffer + */ +void read_variable_data_hdf5(hid_t file, std::string varName, size_t index, double *data); + +/** + * @brief Write data to an hdf5 file + * + * @param group Group within the hdf5 file to write to (e.g., as returned by H5Gcreate) + * @param varName Name of variable to write + * @param dataspace HDF5 dataspace (e.g., as returned by H5Screate_simple) + * @param data Buffer with data to write + */ +void write_variable_data_hdf5(hid_t group, std::string varName, hid_t dataspace, const double *data); + +/** + * @brief Read/write partitioning information + * + * For normal restarts, this is used to write/read the partitioning + * information to ensure that it remains consistent across the + * restart. For serial restarts, it is also used to properly + * distribute the incoming fields. + * + * @todo Refactor this function to make it more generic and fit better + * into the IODataOrganizer paradigm. + */ void partitioning_file_hdf5(std::string mode, const RunConfiguration &config, MPI_Groups *groupsMPI, int nelemGlobal, mfem::Array &partitioning); #endif // IO_HPP_ From 73e97ef8be2fd00147c5d9442ebef820fcc1463c Mon Sep 17 00:00:00 2001 From: "Todd A. Oliver" Date: Thu, 1 Feb 2024 14:34:35 -0600 Subject: [PATCH 14/14] Make underscore style consistent in IOFamily Trailing underscore indicates class member variable. --- src/io.cpp | 46 +++++++++++++++++++++++----------------------- src/io.hpp | 8 ++++---- 2 files changed, 27 insertions(+), 27 deletions(-) diff --git a/src/io.cpp b/src/io.cpp index 37b2cc085..c385047a4 100644 --- a/src/io.cpp +++ b/src/io.cpp @@ -428,7 +428,7 @@ void IOFamily::readDistributeSerializedVariable(hid_t file, const IOVar &var, in // cull out global vdofs - [subtle note]: we only use the // indexing corresponding to the first state vector reference in // gvdofs() since data_serial() holds the solution for a single state vector - serial_fes->GetElementVDofs(gelem, gvdofs); + serial_fes_->GetElementVDofs(gelem, gvdofs); int gdof_start_index = 0; for (int i = 0; i < numDof_per_this_elem; i++) @@ -442,7 +442,7 @@ void IOFamily::readDistributeSerializedVariable(hid_t file, const IOVar &var, in std::vector packedData; for (int gelem = 0; gelem < global_ne_; gelem++) if (partitioning[gelem] == rank) { - serial_fes->GetElementVDofs(gelem, gvdofs); + serial_fes_->GetElementVDofs(gelem, gvdofs); int numDof_per_this_elem = gvdofs.Size() / numStateVars; for (int i = 0; i < numDof_per_this_elem; i++) packedData.push_back(data_serial[gvdofs[i]]); } @@ -559,20 +559,20 @@ void IOFamily::serializeForWrite() { int gelem = locToGlobElem[elem]; this->pfunc_->ParFESpace()->GetElementVDofs(elem, lvdofs); pfunc->GetSubVector(lvdofs, lsoln); - this->serial_fes->GetElementVDofs(gelem, gvdofs); - this->serial_sol->SetSubVector(gvdofs, lsoln); + this->serial_fes_->GetElementVDofs(gelem, gvdofs); + this->serial_sol_->SetSubVector(gvdofs, lsoln); } // have rank 0 receive data from other tasks and copy its own for (int gelem = 0; gelem < global_ne_; gelem++) { int from_rank = partitioning[gelem]; if (from_rank != 0) { - this->serial_fes->GetElementVDofs(gelem, gvdofs); + this->serial_fes_->GetElementVDofs(gelem, gvdofs); lsoln.SetSize(gvdofs.Size()); MPI_Recv(lsoln.HostReadWrite(), gvdofs.Size(), MPI_DOUBLE, from_rank, gelem, comm, MPI_STATUS_IGNORE); - this->serial_sol->SetSubVector(gvdofs, lsoln); + this->serial_sol_->SetSubVector(gvdofs, lsoln); } } @@ -621,14 +621,14 @@ void IOFamily::writeSerial(hid_t file) { // Only get here if need to serialize serializeForWrite(); - // Only rank 0 writes data (which has just been collected into this->serial_sol) + // Only rank 0 writes data (which has just been collected into this->serial_sol_) if (rank0_) { hsize_t dims[1]; hid_t group = -1; hid_t dataspace = -1; - assert(serial_fes != NULL); - dims[0] = serial_fes->GetNDofs(); + assert(serial_fes_ != NULL); + dims[0] = serial_fes_->GetNDofs(); dataspace = H5Screate_simple(1, dims, NULL); assert(dataspace >= 0); @@ -636,8 +636,8 @@ void IOFamily::writeSerial(hid_t file) { assert(group >= 0); // get pointer to raw data - assert(serial_sol != NULL); - const double *data = serial_sol->HostRead(); + assert(serial_sol_ != NULL); + const double *data = serial_sol_->HostRead(); // save raw data for (auto var : vars_) { @@ -689,7 +689,7 @@ void IOFamily::readSerial(hid_t file) { } void IOFamily::readChangeOrder(hid_t file, int read_order) { - if (allowsAuxRestart) { + if (allowsAuxRestart_) { // Set up "auxilliary" FiniteElementCollection, which matches original, but with different order assert(fec_ != NULL); aux_fec_ = fec_->Clone(read_order); @@ -743,12 +743,12 @@ void IOFamily::readChangeOrder(hid_t file, int read_order) { // register a new IO family which maps to a ParGridFunction void IODataOrganizer::registerIOFamily(std::string description, std::string group, ParGridFunction *pfunc, - bool auxRestart, bool _inRestartFile, FiniteElementCollection *fec) { + bool auxRestart, bool inRestartFile, FiniteElementCollection *fec) { IOFamily family(description, group, pfunc); - family.allowsAuxRestart = auxRestart; - family.inRestartFile = _inRestartFile; + family.allowsAuxRestart_ = auxRestart; + family.inRestartFile_ = inRestartFile; family.fec_ = fec; - if (family.allowsAuxRestart) { + if (family.allowsAuxRestart_) { assert(family.fec_ != NULL); } @@ -777,8 +777,8 @@ void IODataOrganizer::initializeSerial(bool root, bool serial, Mesh *serial_mesh // loop through families for (size_t n = 0; n < families_.size(); n++) { IOFamily &fam = families_[n]; - fam.serial_fes = NULL; - fam.serial_sol = NULL; + fam.serial_fes_ = NULL; + fam.serial_sol_ = NULL; // If serial support is requested, need to initialize required data if (supports_serial_) { @@ -788,8 +788,8 @@ void IODataOrganizer::initializeSerial(bool root, bool serial, Mesh *serial_mesh const FiniteElementCollection *fec = fam.pfunc_->ParFESpace()->FEColl(); int numVars = fam.pfunc_->Size() / fam.pfunc_->ParFESpace()->GetNDofs(); - fam.serial_fes = new FiniteElementSpace(serial_mesh, fec, numVars, Ordering::byNODES); - fam.serial_sol = new GridFunction(fam.serial_fes); + fam.serial_fes_ = new FiniteElementSpace(serial_mesh, fec, numVars, Ordering::byNODES); + fam.serial_sol_ = new GridFunction(fam.serial_fes_); // cout<<"I/O organizer for group "<ParFESpace()->GetMyRank(); const bool rank0 = (rank == 0); @@ -887,7 +887,7 @@ void IODataOrganizer::read(hid_t file, bool serial, int read_order) { IODataOrganizer::~IODataOrganizer() { for (size_t n = 0; n < families_.size(); n++) { IOFamily fam = families_[n]; - if (fam.serial_fes != NULL) delete fam.serial_fes; - if (fam.serial_sol != NULL) delete fam.serial_sol; + if (fam.serial_fes_ != NULL) delete fam.serial_fes_; + if (fam.serial_sol_ != NULL) delete fam.serial_sol_; } } diff --git a/src/io.hpp b/src/io.hpp index 556f5cd30..18cb9b4c2 100644 --- a/src/io.hpp +++ b/src/io.hpp @@ -113,16 +113,16 @@ class IOFamily { std::vector vars_; /** Whether this family allows initialization (read) from other order elements */ - bool allowsAuxRestart; + bool allowsAuxRestart_; /** true if the family is to be found in restart files */ - bool inRestartFile; + bool inRestartFile_; /** mfem::FiniteElementSpace for serial mesh (used if serial read and/or write requested) */ - mfem::FiniteElementSpace *serial_fes = nullptr; + mfem::FiniteElementSpace *serial_fes_ = nullptr; /** mfem::GridFunction on serial mesh (used if serial read and/or write requested) */ - mfem::GridFunction *serial_sol = nullptr; + mfem::GridFunction *serial_sol_ = nullptr; /** Map from local to global element numbering (used for serial read/write) */ int *local_to_global_elem_ = nullptr;