From bcaa2766a78f2e027aac11b12842b5ed62bcf294 Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Fri, 15 Dec 2023 12:03:27 -0700 Subject: [PATCH 01/25] allow for disabling of XDMF --- src/outputs/outputs.cpp | 1 + src/outputs/outputs.hpp | 3 ++- src/outputs/parthenon_hdf5.cpp | 10 ++++++---- 3 files changed, 9 insertions(+), 5 deletions(-) diff --git a/src/outputs/outputs.cpp b/src/outputs/outputs.cpp index 33b35a24f284..125b17511fe7 100644 --- a/src/outputs/outputs.cpp +++ b/src/outputs/outputs.cpp @@ -263,6 +263,7 @@ Outputs::Outputs(Mesh *pm, ParameterInput *pin, SimTime *tm) { num_rst_outputs++; } #ifdef ENABLE_HDF5 + op.write_xdmf = pin->GetOrAddBoolean(op.block_name, "write_xdmf", true); pnew_type = new PHDF5Output(op, restart); #else msg << "### FATAL ERROR in Outputs constructor" << std::endl diff --git a/src/outputs/outputs.hpp b/src/outputs/outputs.hpp index 4fd64236a19c..dfa2cfe63f75 100644 --- a/src/outputs/outputs.hpp +++ b/src/outputs/outputs.hpp @@ -64,12 +64,13 @@ struct OutputParameters { bool single_precision_output; bool sparse_seed_nans; int hdf5_compression_level; + bool write_xdmf; // TODO(felker): some of the parameters in this class are not initialized in constructor OutputParameters() : block_number(0), next_time(0.0), dt(-1.0), file_number(0), include_ghost_zones(false), cartesian_vector(false), single_precision_output(false), sparse_seed_nans(false), - hdf5_compression_level(5) {} + hdf5_compression_level(5), write_xdmf(true) {} }; //---------------------------------------------------------------------------------------- diff --git a/src/outputs/parthenon_hdf5.cpp b/src/outputs/parthenon_hdf5.cpp index 004054d8421a..0253aabd2cb4 100644 --- a/src/outputs/parthenon_hdf5.cpp +++ b/src/outputs/parthenon_hdf5.cpp @@ -676,10 +676,12 @@ void PHDF5Output::WriteOutputFileImpl(Mesh *pm, ParameterInput *pin, SimTime *tm } Kokkos::Profiling::popRegion(); // write particle data - Kokkos::Profiling::pushRegion("genXDMF"); - // generate XDMF companion file - XDMF::genXDMF(filename, pm, tm, nx1, nx2, nx3, all_vars_info, swarm_info); - Kokkos::Profiling::popRegion(); // genXDMF + if (output_params.write_xdmf) { + Kokkos::Profiling::pushRegion("genXDMF"); + // generate XDMF companion file + XDMF::genXDMF(filename, pm, tm, nx1, nx2, nx3, all_vars_info, swarm_info); + Kokkos::Profiling::popRegion(); // genXDMF + } Kokkos::Profiling::popRegion(); // WriteOutputFile???Prec } From 56574c80fe71e9af77ec21506c9240fa6d2addcc Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Fri, 15 Dec 2023 18:06:45 -0700 Subject: [PATCH 02/25] Move output loops into PackOrUnpackVar --- src/mesh/domain.hpp | 1 + src/outputs/output_utils.cpp | 35 ++++++++++++++++ src/outputs/output_utils.hpp | 52 ++++++++++++++++++++++++ src/outputs/parthenon_hdf5.cpp | 73 +++++----------------------------- src/parthenon_manager.cpp | 20 ++++------ 5 files changed, 104 insertions(+), 77 deletions(-) diff --git a/src/mesh/domain.hpp b/src/mesh/domain.hpp index 39fe63aea704..e394837be38c 100644 --- a/src/mesh/domain.hpp +++ b/src/mesh/domain.hpp @@ -23,6 +23,7 @@ #include #include +#include "basic_types.hpp" #include "defs.hpp" namespace parthenon { diff --git a/src/outputs/output_utils.cpp b/src/outputs/output_utils.cpp index c6d05451ea29..00aea74153b7 100644 --- a/src/outputs/output_utils.cpp +++ b/src/outputs/output_utils.cpp @@ -100,6 +100,41 @@ AllSwarmInfo::AllSwarmInfo(BlockList_t &block_list, } } +// Tools that can be shared accross Output types +std::vector ComputeXminBlocks(Mesh *pm) { + return FlattenBlockInfo(pm, 3, + [=](MeshBlock *pmb, std::vector &data, int &i) { + auto xmin = pmb->coords.GetXmin(); + data[i++] = xmin[0]; + if (pm->ndim > 1) { + data[i++] = xmin[1]; + } + if (pm->ndim > 2) { + data[i++] = xmin[2]; + } + }); +} + +std::vector ComputeLocs(Mesh *pm) { + return FlattenBlockInfo( + pm, 3, [=](MeshBlock *pmb, std::vector &locs, int &i) { + locs[i++] = pmb->loc.lx1(); + locs[i++] = pmb->loc.lx2(); + locs[i++] = pmb->loc.lx3(); + }); +} + +std::vector ComputeIDsAndFlags(Mesh *pm) { + return FlattenBlockInfo(pm, 5, + [=](MeshBlock *pmb, std::vector &data, int &i) { + data[i++] = pmb->loc.level(); + data[i++] = pmb->gid; + data[i++] = pmb->lid; + data[i++] = pmb->cnghost; + data[i++] = pmb->gflag; + }); +} + // TODO(JMM): may need to generalize this std::size_t MPIPrefixSum(std::size_t local, std::size_t &tot_count) { std::size_t out = 0; diff --git a/src/outputs/output_utils.hpp b/src/outputs/output_utils.hpp index c42a235d6a85..a3b8df2051e1 100644 --- a/src/outputs/output_utils.hpp +++ b/src/outputs/output_utils.hpp @@ -29,8 +29,10 @@ #include // Parthenon +#include "basic_types.hpp" #include "interface/metadata.hpp" #include "interface/variable.hpp" +#include "mesh/domain.hpp" #include "mesh/mesh.hpp" #include "mesh/meshblock.hpp" #include "utils/error_checking.hpp" @@ -200,6 +202,56 @@ struct AllSwarmInfo { bool is_restart); }; +template +std::vector FlattenBlockInfo(Mesh *pm, int shape, Function_t f) { + const int num_blocks_local = static_cast(pm->block_list.size()); + std::vector data(shape * num_blocks_local); + int i = 0; + for (auto &pmb : pm->block_list) { + f(pmb.get(), data, i); + } + return data; +} + +template +void PackOrUnpackVar(MeshBlock *pmb, Variable *pvar, bool do_ghosts, idx_t &idx, + std::vector &data, Function_t f) { + auto v_h = pvar->data.GetHostMirrorAndCopy(); + const auto &Nt = pvar->GetDim(6); + const auto &Nu = pvar->GetDim(5); + const auto &Nv = pvar->GetDim(4); + const IndexDomain theDomain = (do_ghosts ? IndexDomain::entire : IndexDomain::interior); + IndexRange kb, jb, ib; + if (pvar->metadata().Where() == MetadataFlag(Metadata::Cell)) { + kb = pmb->cellbounds.GetBoundsK(theDomain); + jb = pmb->cellbounds.GetBoundsJ(theDomain); + ib = pmb->cellbounds.GetBoundsI(theDomain); + // TODO(JMM): Add topological elements here + } else { // metadata none + kb = {0, pvar->GetDim(3) - 1}; + jb = {0, pvar->GetDim(2) - 1}; + ib = {0, pvar->GetDim(1) - 1}; + } + for (int t = 0; t < Nt; ++t) { + for (int u = 0; u < Nu; ++u) { + for (int v = 0; v < Nv; ++v) { + for (int k = kb.s; k <= kb.e; ++k) { + for (int j = jb.s; j <= jb.e; ++j) { + for (int i = ib.s; i <= ib.e; ++i) { + f(v_h, idx, t, u, v, k, j, i); + idx++; + } + } + } + } + } + } +} + +std::vector ComputeXminBlocks(Mesh *pm); +std::vector ComputeLocs(Mesh *pm); +std::vector ComputeIDsAndFlags(Mesh *pm); + // TODO(JMM): Potentially unsafe if MPI_UNSIGNED_LONG_LONG isn't a size_t // however I think it's probably safe to assume we'll be on systems // where this is the case? diff --git a/src/outputs/parthenon_hdf5.cpp b/src/outputs/parthenon_hdf5.cpp index 0253aabd2cb4..0e733248f08a 100644 --- a/src/outputs/parthenon_hdf5.cpp +++ b/src/outputs/parthenon_hdf5.cpp @@ -247,8 +247,7 @@ void PHDF5Output::WriteOutputFileImpl(Mesh *pm, ParameterInput *pin, SimTime *tm // write Xmin[ndim] for blocks { - std::vector tmpData(num_blocks_local * 3); - ComputeXminBlocks_(pm, tmpData); + std::vector tmpData = OutputUtils::ComputeXminBlocks(pm); local_count[1] = global_count[1] = pm->ndim; HDF5Write2D(gBlocks, "xmin", tmpData.data(), p_loc_offset, p_loc_cnt, p_glob_cnt, pl_xfer); @@ -258,17 +257,15 @@ void PHDF5Output::WriteOutputFileImpl(Mesh *pm, ParameterInput *pin, SimTime *tm { // LOC.lx1,2,3 hsize_t n = 3; - std::vector tmpLoc(num_blocks_local * n); + std::vector tmpLoc = OutputUtils::ComputeLocs(pm); local_count[1] = global_count[1] = n; - ComputeLocs_(pm, tmpLoc); HDF5Write2D(gBlocks, "loc.lx123", tmpLoc.data(), p_loc_offset, p_loc_cnt, p_glob_cnt, pl_xfer); // (LOC.)level, GID, LID, cnghost, gflag n = 5; // this is NOT H5_NDIM - std::vector tmpID(num_blocks_local * n); + std::vector tmpID = OutputUtils::ComputeIDsAndFlags(pm); local_count[1] = global_count[1] = n; - ComputeIDsAndFlags_(pm, tmpID); HDF5Write2D(gBlocks, "loc.level-gid-lid-cnghost-gflag", tmpID.data(), p_loc_offset, p_loc_cnt, p_glob_cnt, pl_xfer); } @@ -495,30 +492,11 @@ void PHDF5Output::WriteOutputFileImpl(Mesh *pm, ParameterInput *pin, SimTime *tm // For reference, if we update the logic here, there's also // a similar block in parthenon_manager.cpp if (v->IsAllocated() && (var_name == v->label())) { - auto v_h = v->data.GetHostMirrorAndCopy(); - for (int t = 0; t < nx6; ++t) { - for (int u = 0; u < nx5; ++u) { - for (int v = 0; v < nx4; ++v) { - if (vinfo.where == MetadataFlag(Metadata::Cell)) { - for (int k = out_kb.s; k <= out_kb.e; ++k) { - for (int j = out_jb.s; j <= out_jb.e; ++j) { - for (int i = out_ib.s; i <= out_ib.e; ++i) { - tmpData[index++] = static_cast(v_h(t, u, v, k, j, i)); - } - } - } - } else { - for (int k = 0; k < vinfo.nx3; ++k) { - for (int j = 0; j < vinfo.nx2; ++j) { - for (int i = 0; i < vinfo.nx1; ++i) { - tmpData[index++] = static_cast(v_h(t, u, v, k, j, i)); - } - } - } - } - } - } - } + OutputUtils::PackOrUnpackVar( + pmb.get(), v.get(), output_params.include_ghost_zones, index, tmpData, + [&](auto v_h, auto index, int t, int u, int v, int k, int j, int i) { + tmpData[index] = static_cast(v_h(t, u, v, k, j, i)); + }); is_allocated = true; break; @@ -720,40 +698,7 @@ std::string PHDF5Output::GenerateFilename_(ParameterInput *pin, SimTime *tm, } return filename; } -// TODO(JMM): Should this live in the base class or output_utils? -void PHDF5Output::ComputeXminBlocks_(Mesh *pm, std::vector &data) { - int i = 0; - for (auto &pmb : pm->block_list) { - auto xmin = pmb->coords.GetXmin(); - data[i++] = xmin[0]; - if (pm->ndim > 1) { - data[i++] = xmin[1]; - } - if (pm->ndim > 2) { - data[i++] = xmin[2]; - } - } -} -// TODO(JMM): Should this live in the base class or output_utils? -void PHDF5Output::ComputeLocs_(Mesh *pm, std::vector &locs) { - int i = 0; - for (auto &pmb : pm->block_list) { - locs[i++] = pmb->loc.lx1(); - locs[i++] = pmb->loc.lx2(); - locs[i++] = pmb->loc.lx3(); - } -} -// TODO(JMM): Should this live in the base class or output_utils? -void PHDF5Output::ComputeIDsAndFlags_(Mesh *pm, std::vector &data) { - int i = 0; - for (auto &pmb : pm->block_list) { - data[i++] = pmb->loc.level(); - data[i++] = pmb->gid; - data[i++] = pmb->lid; - data[i++] = pmb->cnghost; - data[i++] = pmb->gflag; - } -} + // TODO(JMM): Should this live in the base class or output_utils? void PHDF5Output::ComputeCoords_(Mesh *pm, bool face, const IndexRange &ib, const IndexRange &jb, const IndexRange &kb, diff --git a/src/parthenon_manager.cpp b/src/parthenon_manager.cpp index 2329553bfe76..821158ea8f1c 100644 --- a/src/parthenon_manager.cpp +++ b/src/parthenon_manager.cpp @@ -336,6 +336,9 @@ void ParthenonManager::RestartPackages(Mesh &rm, RestartReader &resfile) { // Double note that this also needs to be update in case // we update the HDF5 infrastructure! + // JMM: Version -1 is the old version. We are now on version 3. + // TODO(JMM): Do we need to maintain backwards compatibility + // here? It would be ncie to remove this. if (file_output_format_ver == -1) { for (int k = out_kb.s; k <= out_kb.e; ++k) { for (int j = out_jb.s; j <= out_jb.e; ++j) { @@ -348,19 +351,10 @@ void ParthenonManager::RestartPackages(Mesh &rm, RestartReader &resfile) { } } else if (file_output_format_ver == 2 || file_output_format_ver == HDF5::OUTPUT_VERSION_FORMAT) { - for (int t = 0; t < Nt; ++t) { - for (int u = 0; u < Nu; ++u) { - for (int v = 0; v < Nv; ++v) { - for (int k = out_kb.s; k <= out_kb.e; ++k) { - for (int j = out_jb.s; j <= out_jb.e; ++j) { - for (int i = out_ib.s; i <= out_ib.e; ++i) { - v_h(t, u, v, k, j, i) = tmp[index++]; - } - } - } - } - } - } + OutputUtils::PackOrUnpackVar(pmb.get(), v.get(), resfile.hasGhost, index, tmp, + [&](auto v_h, auto index, int t, int u, int v, int k, + int j, + int i) { v_h(t, u, v, k, j, i) = tmp[index]; }); } else { PARTHENON_THROW("Unknown output format version in restart file.") } From 52da2adcb68c4d3cf25419c838015f26b0c4888e Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Tue, 19 Dec 2023 16:29:50 -0700 Subject: [PATCH 03/25] tiny cleanup with size method in VarInfo struct --- src/outputs/output_utils.hpp | 1 + src/outputs/parthenon_hdf5.cpp | 3 +-- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/outputs/output_utils.hpp b/src/outputs/output_utils.hpp index a3b8df2051e1..d8614f574f83 100644 --- a/src/outputs/output_utils.hpp +++ b/src/outputs/output_utils.hpp @@ -55,6 +55,7 @@ struct VarInfo { bool is_sparse; bool is_vector; std::vector component_labels; + int Size() const { return nx6 * nx5 * nx4 * nx3 * nx2 * nx1; } VarInfo() = delete; diff --git a/src/outputs/parthenon_hdf5.cpp b/src/outputs/parthenon_hdf5.cpp index 0e733248f08a..dcbfeadb6963 100644 --- a/src/outputs/parthenon_hdf5.cpp +++ b/src/outputs/parthenon_hdf5.cpp @@ -394,8 +394,7 @@ void PHDF5Output::WriteOutputFileImpl(Mesh *pm, ParameterInput *pin, SimTime *tm // allocate space for largest size variable int varSize_max = 0; for (auto &vinfo : all_vars_info) { - const int varSize = - vinfo.nx6 * vinfo.nx5 * vinfo.nx4 * vinfo.nx3 * vinfo.nx2 * vinfo.nx1; + const int varSize = vinfo.Size(); varSize_max = std::max(varSize_max, varSize); } From ba333c5f6c037ba1c324f2dd718f6bf668f68906 Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Tue, 19 Dec 2023 16:43:35 -0700 Subject: [PATCH 04/25] move some logic into the mesh --- src/mesh/mesh.hpp | 14 ++++++++++++++ src/outputs/parthenon_hdf5.cpp | 15 +-------------- 2 files changed, 15 insertions(+), 14 deletions(-) diff --git a/src/mesh/mesh.hpp b/src/mesh/mesh.hpp index bf21853d3a82..e8f29f6a35e2 100644 --- a/src/mesh/mesh.hpp +++ b/src/mesh/mesh.hpp @@ -194,6 +194,20 @@ class Mesh { std::vector GetNbList() const noexcept { return nblist; } std::vector GetLocList() const noexcept { return loclist; } + // TODO(JMM): Put in implementation file? + auto GetLevelsAndLogicalLocationsFlat() const noexcept { + std::vector levels, logicalLocations; + levels.reserve(nbtotal); + logicalLocations.reserve(nbtotal * 3); + for (const auto &loc : loclist) { + levels.push_back(loc.level() - GetRootLevel()); + logicalLocations.push_back(loc.lx1()); + logicalLocations.push_back(loc.lx2()); + logicalLocations.push_back(loc.lx3()); + } + return std::make_pair(levels, logicalLocations); + } + void OutputMeshStructure(const int dim, const bool dump_mesh_structure = true); // Ordering here is important to prevent deallocation of pools before boundary diff --git a/src/outputs/parthenon_hdf5.cpp b/src/outputs/parthenon_hdf5.cpp index dcbfeadb6963..e60c68847f05 100644 --- a/src/outputs/parthenon_hdf5.cpp +++ b/src/outputs/parthenon_hdf5.cpp @@ -302,20 +302,7 @@ void PHDF5Output::WriteOutputFileImpl(Mesh *pm, ParameterInput *pin, SimTime *tm // levels and logical locations for all meshblocks on all ranks { Kokkos::Profiling::pushRegion("write levels and locations"); - const auto &loclist = pm->GetLocList(); - - std::vector levels; - levels.reserve(pm->nbtotal); - - std::vector logicalLocations; - logicalLocations.reserve(pm->nbtotal * 3); - - for (const auto &loc : loclist) { - levels.push_back(loc.level() - pm->GetRootLevel()); - logicalLocations.push_back(loc.lx1()); - logicalLocations.push_back(loc.lx2()); - logicalLocations.push_back(loc.lx3()); - } + auto [levels, logicalLocations] = pm->GetLevelsAndLogicalLocationsFlat(); // Only write levels on rank 0 since it has data for all ranks local_count[0] = (Globals::my_rank == 0) ? pm->nbtotal : 0; From 38affd4ebe9b885ded5beb29ca3b0d61d10198b0 Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Wed, 20 Dec 2023 12:09:00 -0700 Subject: [PATCH 05/25] move compute coords into I/O utils --- src/basic_types.hpp | 2 +- src/mesh/domain.hpp | 11 +++++++++++ src/mesh/mesh.cpp | 16 ++++++++++++++++ src/mesh/mesh.hpp | 1 + src/outputs/output_utils.cpp | 30 ++++++++++++++++++++++++++++++ src/outputs/output_utils.hpp | 3 +++ src/outputs/outputs.hpp | 3 --- src/outputs/parthenon_hdf5.cpp | 28 ++-------------------------- 8 files changed, 64 insertions(+), 30 deletions(-) diff --git a/src/basic_types.hpp b/src/basic_types.hpp index 9898140f82b4..7e8c0eabea62 100644 --- a/src/basic_types.hpp +++ b/src/basic_types.hpp @@ -49,7 +49,7 @@ enum class TaskStatus { fail, complete, incomplete, iterate, skip, waiting }; enum class AmrTag : int { derefine = -1, same = 0, refine = 1 }; enum class RefinementOp_t { Prolongation, Restriction, None }; - +enum class CellLevel : int { coarse = -1, same = 0, fine = 1 }; // JMM: Not clear this is the best place for this but it minimizes // circular dependency nonsense. constexpr int NUM_BNDRY_TYPES = 10; diff --git a/src/mesh/domain.hpp b/src/mesh/domain.hpp index e394837be38c..5cf9a79bfbbb 100644 --- a/src/mesh/domain.hpp +++ b/src/mesh/domain.hpp @@ -183,6 +183,17 @@ class IndexShape { : IndexRange{ks(domain, el), ke(domain, el)}; } + template + KOKKOS_INLINE_FUNCTION const IndexRange GetBounds(const IndexDomain &domain, + TE el = TE::CC) const noexcept { + if constexpr (dim == 0) + return GetBoundsI(domain, el); + else if constexpr (dim == 1) + return GetBoundsJ(domain, el); + else + return GetBoundsK(domain, el); + } + KOKKOS_INLINE_FUNCTION int is(const IndexDomain &domain, TE el = TE::CC) const noexcept { switch (domain) { diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp index 33d38bcdeee8..0f46524e82e9 100644 --- a/src/mesh/mesh.cpp +++ b/src/mesh/mesh.cpp @@ -1204,6 +1204,22 @@ int Mesh::GetNumberOfMeshBlockCells() const { } const RegionSize &Mesh::GetBlockSize() const { return base_block_size; } +const IndexShape &Mesh::GetLeafBlockCellBounds(CellLevel level) const { + // TODO(JMM): Luke this is for your Metadata::fine stuff. + PARTHENON_DEBUG_REQUIRE(level != CellLevel::fine, + "Currently no access to finer cellbounds"); + MeshBlock *pmb = block_list[0].get(); + if (level == CellLevel::same) { + return pmb->cellbounds; + // TODO(JMM): + // } else if (level == CellLevel::fine) { + // return pmb->fine_cellbounds; + // } + } else { // if (level == CellLevel::coarse) { + return pmb->c_cellbounds; + } +} + // Functionality re-used in mesh constructor void Mesh::RegisterLoadBalancing_(ParameterInput *pin) { #ifdef MPI_PARALLEL // JMM: Not sure this ifdef is needed diff --git a/src/mesh/mesh.hpp b/src/mesh/mesh.hpp index e8f29f6a35e2..c7d5d9e55e9e 100644 --- a/src/mesh/mesh.hpp +++ b/src/mesh/mesh.hpp @@ -98,6 +98,7 @@ class Mesh { int GetNumberOfMeshBlockCells() const; const RegionSize &GetBlockSize() const; RegionSize GetBlockSize(const LogicalLocation &loc) const; + const IndexShape &GetLeafBlockCellBounds(CellLevel level = CellLevel::same) const; // data bool modified; diff --git a/src/outputs/output_utils.cpp b/src/outputs/output_utils.cpp index 00aea74153b7..42e2aaf7f88e 100644 --- a/src/outputs/output_utils.cpp +++ b/src/outputs/output_utils.cpp @@ -101,6 +101,7 @@ AllSwarmInfo::AllSwarmInfo(BlockList_t &block_list, } // Tools that can be shared accross Output types + std::vector ComputeXminBlocks(Mesh *pm) { return FlattenBlockInfo(pm, 3, [=](MeshBlock *pmb, std::vector &data, int &i) { @@ -135,6 +136,35 @@ std::vector ComputeIDsAndFlags(Mesh *pm) { }); } +// TODO(JMM): I could make this use the other loop +// functionality/high-order functions. but it was more code than this +// for, I think, little benefit. +void ComputeCoords(Mesh *pm, bool face, const IndexRange &ib, const IndexRange &jb, + const IndexRange &kb, std::vector &x, std::vector &y, + std::vector &z) { + const int nx1 = ib.e - ib.s + 1; + const int nx2 = jb.e - jb.s + 1; + const int nx3 = kb.e - kb.s + 1; + const int num_blocks = pm->block_list.size(); + x.resize((nx1 + face) * num_blocks); + y.resize((nx2 + face) * num_blocks); + z.resize((nx3 + face) * num_blocks); + std::size_t idx_x = 0, idx_y = 0, idx_z = 0; + + // note relies on casting of bool to int + for (auto &pmb : pm->block_list) { + for (int i = ib.s; i <= ib.e + face; ++i) { + x[idx_x++] = face ? pmb->coords.Xf<1>(i) : pmb->coords.Xc<1>(i); + } + for (int j = jb.s; j <= jb.e + face; ++j) { + y[idx_y++] = face ? pmb->coords.Xf<2>(j) : pmb->coords.Xc<2>(j); + } + for (int k = kb.s; k <= kb.e + face; ++k) { + z[idx_z++] = face ? pmb->coords.Xf<3>(k) : pmb->coords.Xc<3>(k); + } + } +} + // TODO(JMM): may need to generalize this std::size_t MPIPrefixSum(std::size_t local, std::size_t &tot_count) { std::size_t out = 0; diff --git a/src/outputs/output_utils.hpp b/src/outputs/output_utils.hpp index d8614f574f83..749d6b333a1f 100644 --- a/src/outputs/output_utils.hpp +++ b/src/outputs/output_utils.hpp @@ -249,6 +249,9 @@ void PackOrUnpackVar(MeshBlock *pmb, Variable *pvar, bool do_ghosts, idx_t } } +void ComputeCoords(Mesh *pm, bool face, const IndexRange &ib, const IndexRange &jb, + const IndexRange &kb, std::vector &x, std::vector &y, + std::vector &z); std::vector ComputeXminBlocks(Mesh *pm); std::vector ComputeLocs(Mesh *pm); std::vector ComputeIDsAndFlags(Mesh *pm); diff --git a/src/outputs/outputs.hpp b/src/outputs/outputs.hpp index dfa2cfe63f75..f3c3d2f3436d 100644 --- a/src/outputs/outputs.hpp +++ b/src/outputs/outputs.hpp @@ -214,9 +214,6 @@ class PHDF5Output : public OutputType { void ComputeXminBlocks_(Mesh *pm, std::vector &data); void ComputeLocs_(Mesh *pm, std::vector &locs); void ComputeIDsAndFlags_(Mesh *pm, std::vector &data); - void ComputeCoords_(Mesh *pm, bool face, const IndexRange &ib, const IndexRange &jb, - const IndexRange &kb, std::vector &x, std::vector &y, - std::vector &z); }; //---------------------------------------------------------------------------------------- diff --git a/src/outputs/parthenon_hdf5.cpp b/src/outputs/parthenon_hdf5.cpp index e60c68847f05..ecd8041a42b9 100644 --- a/src/outputs/parthenon_hdf5.cpp +++ b/src/outputs/parthenon_hdf5.cpp @@ -278,11 +278,8 @@ void PHDF5Output::WriteOutputFileImpl(Mesh *pm, ParameterInput *pin, SimTime *tm const H5G gLocations = MakeGroup(file, face ? "/Locations" : "/VolumeLocations"); // write X coordinates - std::vector loc_x((nx1 + face) * num_blocks_local); - std::vector loc_y((nx2 + face) * num_blocks_local); - std::vector loc_z((nx3 + face) * num_blocks_local); - - ComputeCoords_(pm, face, out_ib, out_jb, out_kb, loc_x, loc_y, loc_z); + std::vector loc_x, loc_y, loc_z; + OutputUtils::ComputeCoords(pm, face, out_ib, out_jb, out_kb, loc_x, loc_y, loc_z); local_count[1] = global_count[1] = nx1 + face; HDF5Write2D(gLocations, "x", loc_x.data(), p_loc_offset, p_loc_cnt, p_glob_cnt, @@ -685,27 +682,6 @@ std::string PHDF5Output::GenerateFilename_(ParameterInput *pin, SimTime *tm, return filename; } -// TODO(JMM): Should this live in the base class or output_utils? -void PHDF5Output::ComputeCoords_(Mesh *pm, bool face, const IndexRange &ib, - const IndexRange &jb, const IndexRange &kb, - std::vector &x, std::vector &y, - std::vector &z) { - std::size_t idx_x = 0, idx_y = 0, idx_z = 0; - - // note relies on casting of bool to int - for (auto &pmb : pm->block_list) { - for (int i = ib.s; i <= ib.e + face; ++i) { - x[idx_x++] = face ? pmb->coords.Xf<1>(i) : pmb->coords.Xc<1>(i); - } - for (int j = jb.s; j <= jb.e + face; ++j) { - y[idx_y++] = face ? pmb->coords.Xf<2>(j) : pmb->coords.Xc<2>(j); - } - for (int k = kb.s; k <= kb.e + face; ++k) { - z[idx_z++] = face ? pmb->coords.Xf<3>(k) : pmb->coords.Xc<3>(k); - } - } -} - // explicit template instantiation template void PHDF5Output::WriteOutputFileImpl(Mesh *, ParameterInput *, SimTime *, SignalHandler::OutputSignal); From 90be429ec51a19ded924e59abaa0a099ac3ac5fe Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Wed, 20 Dec 2023 12:15:50 -0700 Subject: [PATCH 06/25] Remove unneeded prototypes from PHDF5Output class --- src/outputs/outputs.hpp | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/outputs/outputs.hpp b/src/outputs/outputs.hpp index f3c3d2f3436d..d76fb03de595 100644 --- a/src/outputs/outputs.hpp +++ b/src/outputs/outputs.hpp @@ -211,9 +211,6 @@ class PHDF5Output : public OutputType { const SignalHandler::OutputSignal signal); const bool restart_; // true if we write a restart file, false for regular output files // TODO(JMM): these methods might want to live in the base class or in output_utils.hpp - void ComputeXminBlocks_(Mesh *pm, std::vector &data); - void ComputeLocs_(Mesh *pm, std::vector &locs); - void ComputeIDsAndFlags_(Mesh *pm, std::vector &data); }; //---------------------------------------------------------------------------------------- From c519b0f80e3f1486a6285b7bcb010a3764ab0537 Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Wed, 20 Dec 2023 16:33:51 -0700 Subject: [PATCH 07/25] Mark TopologicalOffset constexpr so that we can simplify the coords logic --- src/basic_types.hpp | 7 ++++--- src/coordinates/uniform_cartesian.hpp | 13 +++---------- 2 files changed, 7 insertions(+), 13 deletions(-) diff --git a/src/basic_types.hpp b/src/basic_types.hpp index 7e8c0eabea62..db409d6a51f9 100644 --- a/src/basic_types.hpp +++ b/src/basic_types.hpp @@ -152,16 +152,17 @@ inline std::vector GetTopologicalElements(TopologicalType tt if (tt == TT::Face) return {TE::F1, TE::F2, TE::F3}; return {TE::CC}; } + using TE = TopologicalElement; // Returns one if the I coordinate of el is offset from the zone center coordinates, // and zero otherwise -KOKKOS_INLINE_FUNCTION int TopologicalOffsetI(TE el) noexcept { +constexpr int TopologicalOffsetI(TE el) noexcept { return (el == TE::F1 || el == TE::E2 || el == TE::E3 || el == TE::NN); } -KOKKOS_INLINE_FUNCTION int TopologicalOffsetJ(TE el) noexcept { +constexpr int TopologicalOffsetJ(TE el) noexcept { return (el == TE::F2 || el == TE::E3 || el == TE::E1 || el == TE::NN); } -KOKKOS_INLINE_FUNCTION int TopologicalOffsetK(TE el) noexcept { +constexpr int TopologicalOffsetK(TE el) noexcept { return (el == TE::F3 || el == TE::E2 || el == TE::E1 || el == TE::NN); } diff --git a/src/coordinates/uniform_cartesian.hpp b/src/coordinates/uniform_cartesian.hpp index b10cbe0ad06c..af2bb088432e 100644 --- a/src/coordinates/uniform_cartesian.hpp +++ b/src/coordinates/uniform_cartesian.hpp @@ -164,18 +164,11 @@ class UniformCartesian { template KOKKOS_FORCEINLINE_FUNCTION Real X(const int idx) const { - using TE = TopologicalElement; - [[maybe_unused]] bool constexpr X1EDGE = - el == TE::F1 || el == TE::E2 || el == TE::E3 || el == TE::NN; - [[maybe_unused]] bool constexpr X2EDGE = - el == TE::F2 || el == TE::E3 || el == TE::E1 || el == TE::NN; - [[maybe_unused]] bool constexpr X3EDGE = - el == TE::F3 || el == TE::E1 || el == TE::E2 || el == TE::NN; - if constexpr (dir == X1DIR && X1EDGE) { + if constexpr (dir == X1DIR && TopologicalOIffsetI(el)) { return xmin_[dir - 1] + idx * dx_[dir - 1]; // idx - 1/2 - } else if constexpr (dir == X2DIR && X2EDGE) { + } else if constexpr (dir == X2DIR && TopologicalOffsetJ(el)) { return xmin_[dir - 1] + idx * dx_[dir - 1]; // idx - 1/2 - } else if constexpr (dir == X3DIR && X3EDGE) { + } else if constexpr (dir == X3DIR && TopologicalOffsetK(el)) { return xmin_[dir - 1] + idx * dx_[dir - 1]; // idx - 1/2 } else { return xmin_[dir - 1] + (idx + 0.5) * dx_[dir - 1]; // idx From 34a4aaf89174ebf4bf32500ff112649cb51d6b1b Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Wed, 20 Dec 2023 16:37:18 -0700 Subject: [PATCH 08/25] move WriteBlocksMetadata into a helper function --- src/outputs/outputs.hpp | 3 +- src/outputs/parthenon_hdf5.cpp | 81 +++++++++++++++++++--------------- 2 files changed, 48 insertions(+), 36 deletions(-) diff --git a/src/outputs/outputs.hpp b/src/outputs/outputs.hpp index d76fb03de595..cccde00a27b3 100644 --- a/src/outputs/outputs.hpp +++ b/src/outputs/outputs.hpp @@ -209,8 +209,9 @@ class PHDF5Output : public OutputType { private: std::string GenerateFilename_(ParameterInput *pin, SimTime *tm, const SignalHandler::OutputSignal signal); + void WriteBlocksMetadata_(Mesh *pm, hid_t file, const HDF5::H5P &pl, hsize_t offset, + hsize_t max_blocks_global) const; const bool restart_; // true if we write a restart file, false for regular output files - // TODO(JMM): these methods might want to live in the base class or in output_utils.hpp }; //---------------------------------------------------------------------------------------- diff --git a/src/outputs/parthenon_hdf5.cpp b/src/outputs/parthenon_hdf5.cpp index ecd8041a42b9..a43a784fe261 100644 --- a/src/outputs/parthenon_hdf5.cpp +++ b/src/outputs/parthenon_hdf5.cpp @@ -241,36 +241,7 @@ void PHDF5Output::WriteOutputFileImpl(Mesh *pm, ParameterInput *pin, SimTime *tm #endif // write Blocks metadata - { - Kokkos::Profiling::pushRegion("write block metadata"); - const H5G gBlocks = MakeGroup(file, "/Blocks"); - - // write Xmin[ndim] for blocks - { - std::vector tmpData = OutputUtils::ComputeXminBlocks(pm); - local_count[1] = global_count[1] = pm->ndim; - HDF5Write2D(gBlocks, "xmin", tmpData.data(), p_loc_offset, p_loc_cnt, p_glob_cnt, - pl_xfer); - } - - // write Block ID - { - // LOC.lx1,2,3 - hsize_t n = 3; - std::vector tmpLoc = OutputUtils::ComputeLocs(pm); - local_count[1] = global_count[1] = n; - HDF5Write2D(gBlocks, "loc.lx123", tmpLoc.data(), p_loc_offset, p_loc_cnt, - p_glob_cnt, pl_xfer); - - // (LOC.)level, GID, LID, cnghost, gflag - n = 5; // this is NOT H5_NDIM - std::vector tmpID = OutputUtils::ComputeIDsAndFlags(pm); - local_count[1] = global_count[1] = n; - HDF5Write2D(gBlocks, "loc.level-gid-lid-cnghost-gflag", tmpID.data(), p_loc_offset, - p_loc_cnt, p_glob_cnt, pl_xfer); - } - Kokkos::Profiling::popRegion(); // write block metadata - } // Block section + WriteBlocksMetadata_(pm, file, pl_xfer, my_offset, max_blocks_global); // Write mesh coordinates to file Kokkos::Profiling::pushRegion("write mesh coords"); @@ -646,6 +617,11 @@ void PHDF5Output::WriteOutputFileImpl(Mesh *pm, ParameterInput *pin, SimTime *tm Kokkos::Profiling::popRegion(); // WriteOutputFile???Prec } +// explicit template instantiation +template void PHDF5Output::WriteOutputFileImpl(Mesh *, ParameterInput *, SimTime *, + SignalHandler::OutputSignal); +template void PHDF5Output::WriteOutputFileImpl(Mesh *, ParameterInput *, SimTime *, + SignalHandler::OutputSignal); std::string PHDF5Output::GenerateFilename_(ParameterInput *pin, SimTime *tm, const SignalHandler::OutputSignal signal) { @@ -682,11 +658,46 @@ std::string PHDF5Output::GenerateFilename_(ParameterInput *pin, SimTime *tm, return filename; } -// explicit template instantiation -template void PHDF5Output::WriteOutputFileImpl(Mesh *, ParameterInput *, SimTime *, - SignalHandler::OutputSignal); -template void PHDF5Output::WriteOutputFileImpl(Mesh *, ParameterInput *, SimTime *, - SignalHandler::OutputSignal); +void PHDF5Output::WriteBlocksMetadata_(Mesh *pm, hid_t file, const HDF5::H5P &pl, + hsize_t offset, hsize_t max_blocks_global) const { + using namespace HDF5; + Kokkos::Profiling::pushRegion("I/O HDF5: write block metadata"); + const H5G gBlocks = MakeGroup(file, "/Blocks"); + const hsize_t num_blocks_local = pm->block_list.size(); + const hsize_t ndim = pm->ndim; + const hsize_t loc_offset[2] = {offset, 0}; + + // write Xmin[ndim] for blocks + { + // JMM: These arrays chould be shared, but I think this is clearer + // as to what's going on. + hsize_t loc_cnt[2] = {num_blocks_local, ndim}; + hsize_t glob_cnt[2] = {max_blocks_global, ndim}; + + std::vector tmpData = OutputUtils::ComputeXminBlocks(pm); + HDF5Write2D(gBlocks, "xmin", tmpData.data(), &loc_offset[0], &loc_cnt[0], + &glob_cnt[0], pl); + } + + { + // LOC.lx1,2,3 + hsize_t loc_cnt[2] = {num_blocks_local, 3}; + hsize_t glob_cnt[2] = {max_blocks_global, 3}; + std::vector tmpLoc = OutputUtils::ComputeLocs(pm); + HDF5Write2D(gBlocks, "loc.lx123", tmpLoc.data(), &loc_offset[0], &loc_cnt[0], + &glob_cnt[0], pl); + } + + { + // (LOC.)level, GID, LID, cnghost, gflag + hsize_t loc_cnt[2] = {num_blocks_local, 5}; + hsize_t glob_cnt[2] = {max_blocks_global, 5}; + std::vector tmpID = OutputUtils::ComputeIDsAndFlags(pm); + HDF5Write2D(gBlocks, "loc.level-gid-lid-cnghost-gflag", tmpID.data(), &loc_offset[0], + &loc_cnt[0], &glob_cnt[0], pl); + } + Kokkos::Profiling::popRegion(); // write block metadata +} // Utility functions implemented namespace HDF5 { From 93d8e2ee7823083ddd087ccee4a66017162805c3 Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Wed, 20 Dec 2023 16:43:10 -0700 Subject: [PATCH 09/25] oops fix typo --- src/basic_types.hpp | 6 +++--- src/coordinates/uniform_cartesian.hpp | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/basic_types.hpp b/src/basic_types.hpp index db409d6a51f9..919781c6ac17 100644 --- a/src/basic_types.hpp +++ b/src/basic_types.hpp @@ -156,13 +156,13 @@ inline std::vector GetTopologicalElements(TopologicalType tt using TE = TopologicalElement; // Returns one if the I coordinate of el is offset from the zone center coordinates, // and zero otherwise -constexpr int TopologicalOffsetI(TE el) noexcept { +inline constexpr int TopologicalOffsetI(TE el) { return (el == TE::F1 || el == TE::E2 || el == TE::E3 || el == TE::NN); } -constexpr int TopologicalOffsetJ(TE el) noexcept { +inline constexpr int TopologicalOffsetJ(TE el) { return (el == TE::F2 || el == TE::E3 || el == TE::E1 || el == TE::NN); } -constexpr int TopologicalOffsetK(TE el) noexcept { +inline constexpr int TopologicalOffsetK(TE el) { return (el == TE::F3 || el == TE::E2 || el == TE::E1 || el == TE::NN); } diff --git a/src/coordinates/uniform_cartesian.hpp b/src/coordinates/uniform_cartesian.hpp index af2bb088432e..e52e71fa614d 100644 --- a/src/coordinates/uniform_cartesian.hpp +++ b/src/coordinates/uniform_cartesian.hpp @@ -164,7 +164,7 @@ class UniformCartesian { template KOKKOS_FORCEINLINE_FUNCTION Real X(const int idx) const { - if constexpr (dir == X1DIR && TopologicalOIffsetI(el)) { + if constexpr (dir == X1DIR && TopologicalOffsetI(el)) { return xmin_[dir - 1] + idx * dx_[dir - 1]; // idx - 1/2 } else if constexpr (dir == X2DIR && TopologicalOffsetJ(el)) { return xmin_[dir - 1] + idx * dx_[dir - 1]; // idx - 1/2 From 870050893d97a72e18ff67409329af4aae8089a7 Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Wed, 20 Dec 2023 16:57:00 -0700 Subject: [PATCH 10/25] locations into subfunction --- src/outputs/outputs.hpp | 3 ++ src/outputs/parthenon_hdf5.cpp | 58 ++++++++++++++++++++++------------ 2 files changed, 40 insertions(+), 21 deletions(-) diff --git a/src/outputs/outputs.hpp b/src/outputs/outputs.hpp index cccde00a27b3..d87cd658af00 100644 --- a/src/outputs/outputs.hpp +++ b/src/outputs/outputs.hpp @@ -211,6 +211,9 @@ class PHDF5Output : public OutputType { const SignalHandler::OutputSignal signal); void WriteBlocksMetadata_(Mesh *pm, hid_t file, const HDF5::H5P &pl, hsize_t offset, hsize_t max_blocks_global) const; + void WriteCoordinates_(Mesh *pm, const IndexDomain &domain, hid_t file, + const HDF5::H5P &pl, hsize_t offset, + hsize_t max_blocks_global) const; const bool restart_; // true if we write a restart file, false for regular output files }; diff --git a/src/outputs/parthenon_hdf5.cpp b/src/outputs/parthenon_hdf5.cpp index a43a784fe261..4ed199a7f406 100644 --- a/src/outputs/parthenon_hdf5.cpp +++ b/src/outputs/parthenon_hdf5.cpp @@ -244,27 +244,7 @@ void PHDF5Output::WriteOutputFileImpl(Mesh *pm, ParameterInput *pin, SimTime *tm WriteBlocksMetadata_(pm, file, pl_xfer, my_offset, max_blocks_global); // Write mesh coordinates to file - Kokkos::Profiling::pushRegion("write mesh coords"); - for (const bool face : {true, false}) { - const H5G gLocations = MakeGroup(file, face ? "/Locations" : "/VolumeLocations"); - - // write X coordinates - std::vector loc_x, loc_y, loc_z; - OutputUtils::ComputeCoords(pm, face, out_ib, out_jb, out_kb, loc_x, loc_y, loc_z); - - local_count[1] = global_count[1] = nx1 + face; - HDF5Write2D(gLocations, "x", loc_x.data(), p_loc_offset, p_loc_cnt, p_glob_cnt, - pl_xfer); - - local_count[1] = global_count[1] = nx2 + face; - HDF5Write2D(gLocations, "y", loc_y.data(), p_loc_offset, p_loc_cnt, p_glob_cnt, - pl_xfer); - - local_count[1] = global_count[1] = nx3 + face; - HDF5Write2D(gLocations, "z", loc_z.data(), p_loc_offset, p_loc_cnt, p_glob_cnt, - pl_xfer); - } - Kokkos::Profiling::popRegion(); // write mesh coords + WriteCoordinates_(pm, theDomain, file, pl_xfer, my_offset, max_blocks_global); // Write Levels and Logical Locations with the level for each Meshblock loclist contains // levels and logical locations for all meshblocks on all ranks @@ -699,6 +679,42 @@ void PHDF5Output::WriteBlocksMetadata_(Mesh *pm, hid_t file, const HDF5::H5P &pl Kokkos::Profiling::popRegion(); // write block metadata } +void PHDF5Output::WriteCoordinates_(Mesh *pm, const IndexDomain &domain, hid_t file, + const HDF5::H5P &pl, hsize_t offset, + hsize_t max_blocks_global) const { + using namespace HDF5; + Kokkos::Profiling::pushRegion("write mesh coords"); + const IndexShape &shape = pm->GetLeafBlockCellBounds(); + const IndexRange ib = shape.GetBoundsI(domain); + const IndexRange jb = shape.GetBoundsJ(domain); + const IndexRange kb = shape.GetBoundsK(domain); + + const hsize_t num_blocks_local = pm->block_list.size(); + const hsize_t loc_offset[2] = {offset, 0}; + hsize_t loc_cnt[2] = {num_blocks_local, 1}; + hsize_t glob_cnt[2] = {max_blocks_global, 1}; + + for (const bool face : {true, false}) { + const H5G gLocations = MakeGroup(file, face ? "/Locations" : "/VolumeLocations"); + + std::vector loc_x, loc_y, loc_z; + OutputUtils::ComputeCoords(pm, face, ib, jb, kb, loc_x, loc_y, loc_z); + + loc_cnt[1] = glob_cnt[1] = (ib.e - ib.s + 1) + face; + HDF5Write2D(gLocations, "x", loc_x.data(), &loc_offset[0], &loc_cnt[0], &glob_cnt[0], + pl); + + loc_cnt[1] = glob_cnt[1] = (jb.e - jb.s + 1) + face; + HDF5Write2D(gLocations, "y", loc_y.data(), &loc_offset[0], &loc_cnt[0], &glob_cnt[0], + pl); + + loc_cnt[1] = glob_cnt[1] = (kb.e - kb.s + 1) + face; + HDF5Write2D(gLocations, "z", loc_z.data(), &loc_offset[0], &loc_cnt[0], &glob_cnt[0], + pl); + } + Kokkos::Profiling::popRegion(); // write mesh coords +} + // Utility functions implemented namespace HDF5 { std::tuple, std::size_t> From a370038f719c9bb77fd414a640e3ff3cbe489793 Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Wed, 20 Dec 2023 17:21:33 -0700 Subject: [PATCH 11/25] Move WriteLevelsAndLocs to subroutine --- src/mesh/domain.hpp | 11 ---------- src/outputs/outputs.hpp | 2 ++ src/outputs/parthenon_hdf5.cpp | 38 +++++++++++++++++++--------------- 3 files changed, 23 insertions(+), 28 deletions(-) diff --git a/src/mesh/domain.hpp b/src/mesh/domain.hpp index 5cf9a79bfbbb..e394837be38c 100644 --- a/src/mesh/domain.hpp +++ b/src/mesh/domain.hpp @@ -183,17 +183,6 @@ class IndexShape { : IndexRange{ks(domain, el), ke(domain, el)}; } - template - KOKKOS_INLINE_FUNCTION const IndexRange GetBounds(const IndexDomain &domain, - TE el = TE::CC) const noexcept { - if constexpr (dim == 0) - return GetBoundsI(domain, el); - else if constexpr (dim == 1) - return GetBoundsJ(domain, el); - else - return GetBoundsK(domain, el); - } - KOKKOS_INLINE_FUNCTION int is(const IndexDomain &domain, TE el = TE::CC) const noexcept { switch (domain) { diff --git a/src/outputs/outputs.hpp b/src/outputs/outputs.hpp index d87cd658af00..6153d6a0ec96 100644 --- a/src/outputs/outputs.hpp +++ b/src/outputs/outputs.hpp @@ -214,6 +214,8 @@ class PHDF5Output : public OutputType { void WriteCoordinates_(Mesh *pm, const IndexDomain &domain, hid_t file, const HDF5::H5P &pl, hsize_t offset, hsize_t max_blocks_global) const; + void WriteLevelsAndLocs_(Mesh *pm, hid_t file, const HDF5::H5P &pl, hsize_t offset, + hsize_t max_blocks_global) const; const bool restart_; // true if we write a restart file, false for regular output files }; diff --git a/src/outputs/parthenon_hdf5.cpp b/src/outputs/parthenon_hdf5.cpp index 4ed199a7f406..0162024ea66b 100644 --- a/src/outputs/parthenon_hdf5.cpp +++ b/src/outputs/parthenon_hdf5.cpp @@ -248,23 +248,7 @@ void PHDF5Output::WriteOutputFileImpl(Mesh *pm, ParameterInput *pin, SimTime *tm // Write Levels and Logical Locations with the level for each Meshblock loclist contains // levels and logical locations for all meshblocks on all ranks - { - Kokkos::Profiling::pushRegion("write levels and locations"); - auto [levels, logicalLocations] = pm->GetLevelsAndLogicalLocationsFlat(); - - // Only write levels on rank 0 since it has data for all ranks - local_count[0] = (Globals::my_rank == 0) ? pm->nbtotal : 0; - HDF5WriteND(file, "Levels", levels.data(), 1, local_offset.data(), local_count.data(), - global_count.data(), pl_xfer, H5P_DEFAULT); - - local_count[1] = global_count[1] = 3; - HDF5Write2D(file, "LogicalLocations", logicalLocations.data(), local_offset.data(), - local_count.data(), global_count.data(), pl_xfer); - - // reset for collective output - local_count[0] = num_blocks_local; - Kokkos::Profiling::popRegion(); // write levels and locations - } + WriteLevelsAndLocs_(pm, file, pl_xfer, my_offset, max_blocks_global); // -------------------------------------------------------------------------------- // // WRITING VARIABLES DATA // @@ -715,6 +699,26 @@ void PHDF5Output::WriteCoordinates_(Mesh *pm, const IndexDomain &domain, hid_t f Kokkos::Profiling::popRegion(); // write mesh coords } +void PHDF5Output::WriteLevelsAndLocs_(Mesh *pm, hid_t file, const HDF5::H5P &pl, + hsize_t offset, hsize_t max_blocks_global) const { + using namespace HDF5; + Kokkos::Profiling::pushRegion("write levels and locations"); + auto [levels, logicalLocations] = pm->GetLevelsAndLogicalLocationsFlat(); + + // Only write levels on rank 0 since it has data for all ranks + const hsize_t num_blocks_local = pm->block_list.size(); + const hsize_t loc_offset[2] = {offset, 0}; + const hsize_t loc_cnt[2] = {(Globals::my_rank == 0) ? max_blocks_global : 0, 3}; + const hsize_t glob_cnt[2] = {max_blocks_global, 3}; + + HDF5Write1D(file, "Levels", levels.data(), &loc_offset[0], &loc_cnt[0], &glob_cnt[0], + pl); + HDF5Write2D(file, "LogicalLocations", logicalLocations.data(), &loc_offset[0], + &loc_cnt[0], &glob_cnt[0], pl); + + Kokkos::Profiling::popRegion(); // write levels and locations +} + // Utility functions implemented namespace HDF5 { std::tuple, std::size_t> From ad1c5ebb2b04d1f8b6bc6324209e9057e1bf2f06 Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Wed, 20 Dec 2023 17:43:10 -0700 Subject: [PATCH 12/25] Fully disentangle the HDF5 shapes in the file --- src/outputs/outputs.hpp | 4 ++ src/outputs/parthenon_hdf5.cpp | 91 +++++++++++++++++----------------- 2 files changed, 49 insertions(+), 46 deletions(-) diff --git a/src/outputs/outputs.hpp b/src/outputs/outputs.hpp index 6153d6a0ec96..afc5583c003c 100644 --- a/src/outputs/outputs.hpp +++ b/src/outputs/outputs.hpp @@ -216,6 +216,10 @@ class PHDF5Output : public OutputType { hsize_t max_blocks_global) const; void WriteLevelsAndLocs_(Mesh *pm, hid_t file, const HDF5::H5P &pl, hsize_t offset, hsize_t max_blocks_global) const; + void WriteSparseInfo_(Mesh *pm, hbool_t *sparse_allocated, + const std::vector &sparse_names, hsize_t num_sparse, + hid_t file, const HDF5::H5P &pl, size_t offset, + hsize_t max_blocks_global) const; const bool restart_; // true if we write a restart file, false for regular output files }; diff --git a/src/outputs/parthenon_hdf5.cpp b/src/outputs/parthenon_hdf5.cpp index 0162024ea66b..74434fae4b57 100644 --- a/src/outputs/parthenon_hdf5.cpp +++ b/src/outputs/parthenon_hdf5.cpp @@ -217,19 +217,6 @@ void PHDF5Output::WriteOutputFileImpl(Mesh *pm, ParameterInput *pin, SimTime *tm my_offset += nblist[i]; } - const std::array local_offset({my_offset, 0, 0, 0, 0, 0, 0}); - - // these can vary by data set, except index 0 is always the same - std::array local_count( - {static_cast(num_blocks_local), 1, 1, 1, 1, 1, 1}); - std::array global_count( - {static_cast(max_blocks_global), 1, 1, 1, 1, 1, 1}); - - // for convenience - const hsize_t *const p_loc_offset = local_offset.data(); - const hsize_t *const p_loc_cnt = local_count.data(); - const hsize_t *const p_glob_cnt = global_count.data(); - H5P const pl_xfer = H5P::FromHIDCheck(H5Pcreate(H5P_DATASET_XFER)); H5P const pl_dcreate = H5P::FromHIDCheck(H5Pcreate(H5P_DATASET_CREATE)); @@ -240,14 +227,8 @@ void PHDF5Output::WriteOutputFileImpl(Mesh *pm, ParameterInput *pin, SimTime *tm PARTHENON_HDF5_CHECK(H5Pset_dxpl_mpio(pl_xfer, H5FD_MPIO_COLLECTIVE)); #endif - // write Blocks metadata WriteBlocksMetadata_(pm, file, pl_xfer, my_offset, max_blocks_global); - - // Write mesh coordinates to file WriteCoordinates_(pm, theDomain, file, pl_xfer, my_offset, max_blocks_global); - - // Write Levels and Logical Locations with the level for each Meshblock loclist contains - // levels and logical locations for all meshblocks on all ranks WriteLevelsAndLocs_(pm, file, pl_xfer, my_offset, max_blocks_global); // -------------------------------------------------------------------------------- // @@ -320,13 +301,6 @@ void PHDF5Output::WriteOutputFileImpl(Mesh *pm, ParameterInput *pin, SimTime *tm using OutT = typename std::conditional::type; std::vector tmpData(varSize_max * num_blocks_local); - // create persistent spaces - local_count[0] = num_blocks_local; - global_count[0] = max_blocks_global; - local_count[4] = global_count[4] = nx3; - local_count[5] = global_count[5] = nx2; - local_count[6] = global_count[6] = nx1; - // for each variable we write for (auto &vinfo : all_vars_info) { Kokkos::Profiling::pushRegion("write variable loop"); @@ -338,9 +312,21 @@ void PHDF5Output::WriteOutputFileImpl(Mesh *pm, ParameterInput *pin, SimTime *tm const hsize_t nx5 = vinfo.nx5; const hsize_t nx4 = vinfo.nx4; - local_count[1] = global_count[1] = nx6; - local_count[2] = global_count[2] = nx5; - local_count[3] = global_count[3] = nx4; + hsize_t local_offset[H5_NDIM] = {my_offset, 0, 0, 0, 0, 0, 0}; + hsize_t local_count[H5_NDIM] = {static_cast(num_blocks_local), + static_cast(nx6), + static_cast(nx5), + static_cast(nx4), + static_cast(nx3), + static_cast(nx2), + static_cast(nx1)}; + hsize_t global_count[H5_NDIM] = {static_cast(max_blocks_global), + static_cast(nx6), + static_cast(nx5), + static_cast(nx4), + static_cast(nx3), + static_cast(nx2), + static_cast(nx1)}; std::vector alldims({nx6, nx5, nx4, static_cast(vinfo.nx3), static_cast(vinfo.nx2), @@ -451,8 +437,8 @@ void PHDF5Output::WriteOutputFileImpl(Mesh *pm, ParameterInput *pin, SimTime *tm Kokkos::Profiling::pushRegion("write variable data"); // write data to file - HDF5WriteND(file, var_name, tmpData.data(), ndim, p_loc_offset, p_loc_cnt, p_glob_cnt, - pl_xfer, pl_dcreate); + HDF5WriteND(file, var_name, tmpData.data(), ndim, &local_offset[0], &local_count[0], + &global_count[0], pl_xfer, pl_dcreate); Kokkos::Profiling::popRegion(); // write variable data Kokkos::Profiling::popRegion(); // write variable loop } @@ -489,21 +475,9 @@ void PHDF5Output::WriteOutputFileImpl(Mesh *pm, ParameterInput *pin, SimTime *tm // write SparseInfo and SparseFields (we can't write a zero-size dataset, so only write // this if we have sparse fields) if (num_sparse > 0) { - Kokkos::Profiling::pushRegion("write sparse info"); - local_count[1] = global_count[1] = num_sparse; - - HDF5Write2D(file, "SparseInfo", sparse_allocated.get(), p_loc_offset, p_loc_cnt, - p_glob_cnt, pl_xfer); - - // write names of sparse fields as attribute, first convert to vector of const char* - std::vector names(num_sparse); - for (size_t i = 0; i < num_sparse; ++i) - names[i] = sparse_names[i].c_str(); - - const H5D dset = H5D::FromHIDCheck(H5Dopen2(file, "SparseInfo", H5P_DEFAULT)); - HDF5WriteAttribute("SparseFields", names, dset); - Kokkos::Profiling::popRegion(); // write sparse info - } // SparseInfo and SparseFields sections + WriteSparseInfo_(pm, sparse_allocated.get(), sparse_names, num_sparse, file, pl_xfer, + my_offset, max_blocks_global); + } // SparseInfo and SparseFields sections // -------------------------------------------------------------------------------- // // WRITING PARTICLE DATA // @@ -719,6 +693,31 @@ void PHDF5Output::WriteLevelsAndLocs_(Mesh *pm, hid_t file, const HDF5::H5P &pl, Kokkos::Profiling::popRegion(); // write levels and locations } +void PHDF5Output::WriteSparseInfo_(Mesh *pm, hbool_t *sparse_allocated, + const std::vector &sparse_names, + hsize_t num_sparse, hid_t file, const HDF5::H5P &pl, + size_t offset, hsize_t max_blocks_global) const { + using namespace HDF5; + Kokkos::Profiling::pushRegion("write sparse info"); + + const hsize_t num_blocks_local = pm->block_list.size(); + const hsize_t loc_offset[2] = {offset, num_sparse}; + const hsize_t loc_cnt[2] = {num_blocks_local, num_sparse}; + const hsize_t glob_cnt[2] = {max_blocks_global, num_sparse}; + + HDF5Write2D(file, "SparseInfo", sparse_allocated, &loc_offset[0], &loc_cnt[0], + &glob_cnt[0], pl); + + // write names of sparse fields as attribute, first convert to vector of const char* + std::vector names(num_sparse); + for (size_t i = 0; i < num_sparse; ++i) + names[i] = sparse_names[i].c_str(); + + const H5D dset = H5D::FromHIDCheck(H5Dopen2(file, "SparseInfo", H5P_DEFAULT)); + HDF5WriteAttribute("SparseFields", names, dset); + Kokkos::Profiling::popRegion(); // write sparse info +} + // Utility functions implemented namespace HDF5 { std::tuple, std::size_t> From 3b3f8aecfa44310d3652dc015bef5de74edabf02 Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Wed, 20 Dec 2023 17:48:52 -0700 Subject: [PATCH 13/25] offset in 1 direction basically always 0 --- src/outputs/parthenon_hdf5.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/outputs/parthenon_hdf5.cpp b/src/outputs/parthenon_hdf5.cpp index 74434fae4b57..7d6cf7bfee25 100644 --- a/src/outputs/parthenon_hdf5.cpp +++ b/src/outputs/parthenon_hdf5.cpp @@ -701,7 +701,7 @@ void PHDF5Output::WriteSparseInfo_(Mesh *pm, hbool_t *sparse_allocated, Kokkos::Profiling::pushRegion("write sparse info"); const hsize_t num_blocks_local = pm->block_list.size(); - const hsize_t loc_offset[2] = {offset, num_sparse}; + const hsize_t loc_offset[2] = {offset, 0}; const hsize_t loc_cnt[2] = {num_blocks_local, num_sparse}; const hsize_t glob_cnt[2] = {max_blocks_global, num_sparse}; From 872eb4bcc871ffc2f18f7b79c8ec5145beae616b Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Tue, 9 Jan 2024 17:12:18 -0700 Subject: [PATCH 14/25] pgrete comments 1 --- doc/sphinx/src/outputs.rst | 1 + src/outputs/output_utils.hpp | 8 ++++---- src/outputs/outputs.hpp | 2 +- src/parthenon_manager.cpp | 2 ++ 4 files changed, 8 insertions(+), 5 deletions(-) diff --git a/doc/sphinx/src/outputs.rst b/doc/sphinx/src/outputs.rst index c01c957fdcf5..170c83ebac7e 100644 --- a/doc/sphinx/src/outputs.rst +++ b/doc/sphinx/src/outputs.rst @@ -46,6 +46,7 @@ look like file_type = hdf5 + write_xdmf = true # Determines whether xdmf annotations are output # nonexistent variables/swarms are ignored variables = density, velocity, & # comments are still ok energy # notice the & continuation character diff --git a/src/outputs/output_utils.hpp b/src/outputs/output_utils.hpp index 749d6b333a1f..22723e96b5ac 100644 --- a/src/outputs/output_utils.hpp +++ b/src/outputs/output_utils.hpp @@ -221,12 +221,12 @@ void PackOrUnpackVar(MeshBlock *pmb, Variable *pvar, bool do_ghosts, idx_t const auto &Nt = pvar->GetDim(6); const auto &Nu = pvar->GetDim(5); const auto &Nv = pvar->GetDim(4); - const IndexDomain theDomain = (do_ghosts ? IndexDomain::entire : IndexDomain::interior); + const IndexDomain domain = (do_ghosts ? IndexDomain::entire : IndexDomain::interior); IndexRange kb, jb, ib; if (pvar->metadata().Where() == MetadataFlag(Metadata::Cell)) { - kb = pmb->cellbounds.GetBoundsK(theDomain); - jb = pmb->cellbounds.GetBoundsJ(theDomain); - ib = pmb->cellbounds.GetBoundsI(theDomain); + kb = pmb->cellbounds.GetBoundsK(domain); + jb = pmb->cellbounds.GetBoundsJ(domain); + ib = pmb->cellbounds.GetBoundsI(domain); // TODO(JMM): Add topological elements here } else { // metadata none kb = {0, pvar->GetDim(3) - 1}; diff --git a/src/outputs/outputs.hpp b/src/outputs/outputs.hpp index afc5583c003c..0bc85ee909d0 100644 --- a/src/outputs/outputs.hpp +++ b/src/outputs/outputs.hpp @@ -70,7 +70,7 @@ struct OutputParameters { : block_number(0), next_time(0.0), dt(-1.0), file_number(0), include_ghost_zones(false), cartesian_vector(false), single_precision_output(false), sparse_seed_nans(false), - hdf5_compression_level(5), write_xdmf(true) {} + hdf5_compression_level(5), write_xdmf(false) {} }; //---------------------------------------------------------------------------------------- diff --git a/src/parthenon_manager.cpp b/src/parthenon_manager.cpp index 821158ea8f1c..53a4887faec1 100644 --- a/src/parthenon_manager.cpp +++ b/src/parthenon_manager.cpp @@ -340,6 +340,8 @@ void ParthenonManager::RestartPackages(Mesh &rm, RestartReader &resfile) { // TODO(JMM): Do we need to maintain backwards compatibility // here? It would be ncie to remove this. if (file_output_format_ver == -1) { + PARTHENON_WARN("This file output format version is deprecrated and will be " + "removed in a future release."); for (int k = out_kb.s; k <= out_kb.e; ++k) { for (int j = out_jb.s; j <= out_jb.e; ++j) { for (int i = out_ib.s; i <= out_ib.e; ++i) { From 553f19eb01010e050c6fef19951c8ae9c0ba87f3 Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Tue, 9 Jan 2024 17:16:59 -0700 Subject: [PATCH 15/25] 3->ndim --- src/outputs/output_utils.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/outputs/output_utils.cpp b/src/outputs/output_utils.cpp index 42e2aaf7f88e..2869a5c42578 100644 --- a/src/outputs/output_utils.cpp +++ b/src/outputs/output_utils.cpp @@ -103,7 +103,7 @@ AllSwarmInfo::AllSwarmInfo(BlockList_t &block_list, // Tools that can be shared accross Output types std::vector ComputeXminBlocks(Mesh *pm) { - return FlattenBlockInfo(pm, 3, + return FlattenBlockInfo(pm, pm->ndim, [=](MeshBlock *pmb, std::vector &data, int &i) { auto xmin = pmb->coords.GetXmin(); data[i++] = xmin[0]; From 802e4d255adae055c135d607dc41c1e0b87215f2 Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Tue, 9 Jan 2024 17:19:07 -0700 Subject: [PATCH 16/25] changelog --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 417305e6bba0..897acba89ea8 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -17,6 +17,7 @@ - [[PR978]](https://github.com/parthenon-hpc-lab/parthenon/pull/978) remove erroneous sparse check ### Infrastructure (changes irrelevant to downstream codes) +- [[PR 990]](https://github.com/parthenon-hpc-lab/parthenon/pull/990) Partial refactor of HDF5 I/O code for readability/extendability - [[PR 982]](https://github.com/parthenon-hpc-lab/parthenon/pull/982) add some gut check testing for parthenon-VIBE ### Removed (removing behavior/API/varaibles/...) From 11d832c85e826577c655ab804f7b73e04e6d8042 Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Tue, 9 Jan 2024 17:21:11 -0700 Subject: [PATCH 17/25] CC --- src/basic_types.hpp | 2 +- src/coordinates/uniform_cartesian.hpp | 2 +- src/mesh/domain.hpp | 2 +- src/mesh/mesh.cpp | 2 +- src/mesh/mesh.hpp | 2 +- src/outputs/output_utils.cpp | 2 +- src/outputs/output_utils.hpp | 2 +- src/outputs/outputs.cpp | 2 +- src/outputs/outputs.hpp | 2 +- src/outputs/parthenon_hdf5.cpp | 2 +- src/parthenon_manager.cpp | 2 +- 11 files changed, 11 insertions(+), 11 deletions(-) diff --git a/src/basic_types.hpp b/src/basic_types.hpp index 919781c6ac17..c2883286f320 100644 --- a/src/basic_types.hpp +++ b/src/basic_types.hpp @@ -1,5 +1,5 @@ //======================================================================================== -// (C) (or copyright) 2021-2023. Triad National Security, LLC. All rights reserved. +// (C) (or copyright) 2021-2024. Triad National Security, LLC. All rights reserved. // // This program was produced under U.S. Government contract 89233218CNA000001 for Los // Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC diff --git a/src/coordinates/uniform_cartesian.hpp b/src/coordinates/uniform_cartesian.hpp index e52e71fa614d..51b06467fae3 100644 --- a/src/coordinates/uniform_cartesian.hpp +++ b/src/coordinates/uniform_cartesian.hpp @@ -1,5 +1,5 @@ //======================================================================================== -// (C) (or copyright) 2023. Triad National Security, LLC. All rights reserved. +// (C) (or copyright) 2023-2024. Triad National Security, LLC. All rights reserved. // // This program was produced under U.S. Government contract 89233218CNA000001 for Los // Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC diff --git a/src/mesh/domain.hpp b/src/mesh/domain.hpp index e394837be38c..49ec1a37af67 100644 --- a/src/mesh/domain.hpp +++ b/src/mesh/domain.hpp @@ -3,7 +3,7 @@ // Copyright(C) 2014 James M. Stone and other code contributors // Licensed under the 3-clause BSD License, see LICENSE file for details //======================================================================================== -// (C) (or copyright) 2020. Triad National Security, LLC. All rights reserved. +// (C) (or copyright) 2020-2024. Triad National Security, LLC. All rights reserved. // // This program was produced under U.S. Government contract 89233218CNA000001 for Los // Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp index 0f46524e82e9..35c56ac7505c 100644 --- a/src/mesh/mesh.cpp +++ b/src/mesh/mesh.cpp @@ -3,7 +3,7 @@ // Copyright(C) 2014 James M. Stone and other code contributors // Licensed under the 3-clause BSD License, see LICENSE file for details //======================================================================================== -// (C) (or copyright) 2020-2023. Triad National Security, LLC. All rights reserved. +// (C) (or copyright) 2020-2024. Triad National Security, LLC. All rights reserved. // // This program was produced under U.S. Government contract 89233218CNA000001 for Los // Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC diff --git a/src/mesh/mesh.hpp b/src/mesh/mesh.hpp index c7d5d9e55e9e..58b1dfc99094 100644 --- a/src/mesh/mesh.hpp +++ b/src/mesh/mesh.hpp @@ -3,7 +3,7 @@ // Copyright(C) 2014 James M. Stone and other code contributors // Licensed under the 3-clause BSD License, see LICENSE file for details //======================================================================================== -// (C) (or copyright) 2020-2023. Triad National Security, LLC. All rights reserved. +// (C) (or copyright) 2020-2024. Triad National Security, LLC. All rights reserved. // // This program was produced under U.S. Government contract 89233218CNA000001 for Los // Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC diff --git a/src/outputs/output_utils.cpp b/src/outputs/output_utils.cpp index 2869a5c42578..64c099d85617 100644 --- a/src/outputs/output_utils.cpp +++ b/src/outputs/output_utils.cpp @@ -3,7 +3,7 @@ // Copyright(C) 2023 The Parthenon collaboration // Licensed under the 3-clause BSD License, see LICENSE file for details //======================================================================================== -// (C) (or copyright) 2020-2023. Triad National Security, LLC. All rights reserved. +// (C) (or copyright) 2020-2024. Triad National Security, LLC. All rights reserved. // // This program was produced under U.S. Government contract 89233218CNA000001 for Los // Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC diff --git a/src/outputs/output_utils.hpp b/src/outputs/output_utils.hpp index 22723e96b5ac..a768ead308e7 100644 --- a/src/outputs/output_utils.hpp +++ b/src/outputs/output_utils.hpp @@ -3,7 +3,7 @@ // Copyright(C) 2023 The Parthenon collaboration // Licensed under the 3-clause BSD License, see LICENSE file for details //======================================================================================== -// (C) (or copyright) 2020-2023. Triad National Security, LLC. All rights reserved. +// (C) (or copyright) 2020-2024. Triad National Security, LLC. All rights reserved. // // This program was produced under U.S. Government contract 89233218CNA000001 for Los // Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC diff --git a/src/outputs/outputs.cpp b/src/outputs/outputs.cpp index 125b17511fe7..2df7a9beb3c9 100644 --- a/src/outputs/outputs.cpp +++ b/src/outputs/outputs.cpp @@ -1,6 +1,6 @@ //======================================================================================== // Parthenon performance portable AMR framework -// Copyright(C) 2020-2023 The Parthenon collaboration +// Copyright(C) 2020-2024 The Parthenon collaboration // Licensed under the 3-clause BSD License, see LICENSE file for details //======================================================================================== // Athena++ astrophysical MHD code diff --git a/src/outputs/outputs.hpp b/src/outputs/outputs.hpp index 0bc85ee909d0..ff268926d7e1 100644 --- a/src/outputs/outputs.hpp +++ b/src/outputs/outputs.hpp @@ -3,7 +3,7 @@ // Copyright(C) 2014 James M. Stone and other code contributors // Licensed under the 3-clause BSD License, see LICENSE file for details //======================================================================================== -// (C) (or copyright) 2020-2023. Triad National Security, LLC. All rights reserved. +// (C) (or copyright) 2020-2024. Triad National Security, LLC. All rights reserved. // // This program was produced under U.S. Government contract 89233218CNA000001 for Los // Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC diff --git a/src/outputs/parthenon_hdf5.cpp b/src/outputs/parthenon_hdf5.cpp index 7d6cf7bfee25..58f80eefdf1b 100644 --- a/src/outputs/parthenon_hdf5.cpp +++ b/src/outputs/parthenon_hdf5.cpp @@ -3,7 +3,7 @@ // Copyright(C) 2020-2023 The Parthenon collaboration // Licensed under the 3-clause BSD License, see LICENSE file for details //======================================================================================== -// (C) (or copyright) 2020-2023. Triad National Security, LLC. All rights reserved. +// (C) (or copyright) 2020-2024. Triad National Security, LLC. All rights reserved. // // This program was produced under U.S. Government contract 89233218CNA000001 for Los // Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC diff --git a/src/parthenon_manager.cpp b/src/parthenon_manager.cpp index 53a4887faec1..bdb7b6a0b56e 100644 --- a/src/parthenon_manager.cpp +++ b/src/parthenon_manager.cpp @@ -3,7 +3,7 @@ // Copyright(C) 2020-2023 The Parthenon collaboration // Licensed under the 3-clause BSD License, see LICENSE file for details //======================================================================================== -// (C) (or copyright) 2020-2023. Triad National Security, LLC. All rights reserved. +// (C) (or copyright) 2020-2024. Triad National Security, LLC. All rights reserved. // // This program was produced under U.S. Government contract 89233218CNA000001 for Los // Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC From 9945b5c0aa3bb3413f1a2fe88a2cd4fdc7d4a6db Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Tue, 9 Jan 2024 17:29:31 -0700 Subject: [PATCH 18/25] deprecation --- src/parthenon_manager.cpp | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/parthenon_manager.cpp b/src/parthenon_manager.cpp index bdb7b6a0b56e..5cb5d3ce29e9 100644 --- a/src/parthenon_manager.cpp +++ b/src/parthenon_manager.cpp @@ -336,9 +336,6 @@ void ParthenonManager::RestartPackages(Mesh &rm, RestartReader &resfile) { // Double note that this also needs to be update in case // we update the HDF5 infrastructure! - // JMM: Version -1 is the old version. We are now on version 3. - // TODO(JMM): Do we need to maintain backwards compatibility - // here? It would be ncie to remove this. if (file_output_format_ver == -1) { PARTHENON_WARN("This file output format version is deprecrated and will be " "removed in a future release."); From abf16afce27661159a66c1025398b2f20d099567 Mon Sep 17 00:00:00 2001 From: Patrick Mullen Date: Wed, 24 Jan 2024 17:55:03 -0700 Subject: [PATCH 19/25] Add a post-initialization hook --- src/application_input.hpp | 4 +++- src/mesh/mesh.cpp | 27 ++++++++++++++++++++++++++- src/mesh/mesh.hpp | 3 ++- src/mesh/meshblock.cpp | 8 +++++++- src/mesh/meshblock.hpp | 4 +++- src/pgen/default_pgen.cpp | 11 +++++++++-- 6 files changed, 50 insertions(+), 7 deletions(-) diff --git a/src/application_input.hpp b/src/application_input.hpp index a9ec96c1551e..a7cbaf4fe9c7 100644 --- a/src/application_input.hpp +++ b/src/application_input.hpp @@ -1,5 +1,5 @@ //======================================================================================== -// (C) (or copyright) 2020-2022. Triad National Security, LLC. All rights reserved. +// (C) (or copyright) 2020-2024. Triad National Security, LLC. All rights reserved. // // This program was produced under U.S. Government contract 89233218CNA000001 for Los // Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC @@ -36,6 +36,7 @@ struct ApplicationInput { std::function InitUserMeshData = nullptr; std::function *)> MeshProblemGenerator = nullptr; + std::function *)> MeshPostInitialization = nullptr; std::function PreStepMeshUserWorkInLoop = nullptr; @@ -57,6 +58,7 @@ struct ApplicationInput { InitApplicationMeshBlockData = nullptr; std::function InitMeshBlockUserData = nullptr; std::function ProblemGenerator = nullptr; + std::function *)> PostInitialization = nullptr; std::function MeshBlockUserWorkBeforeOutput = nullptr; }; diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp index 7115a0cde719..7fbc4b9f8e47 100644 --- a/src/mesh/mesh.cpp +++ b/src/mesh/mesh.cpp @@ -3,7 +3,7 @@ // Copyright(C) 2014 James M. Stone and other code contributors // Licensed under the 3-clause BSD License, see LICENSE file for details //======================================================================================== -// (C) (or copyright) 2020-2023. Triad National Security, LLC. All rights reserved. +// (C) (or copyright) 2020-2024. Triad National Security, LLC. All rights reserved. // // This program was produced under U.S. Government contract 89233218CNA000001 for Los // Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC @@ -156,6 +156,9 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, if (app_in->MeshProblemGenerator != nullptr) { ProblemGenerator = app_in->MeshProblemGenerator; } + if (app_in->MeshPostInitialization != nullptr) { + PostInitialization = app_in->MeshPostInitialization; + } if (app_in->PreStepMeshUserWorkInLoop != nullptr) { PreStepUserWorkInLoop = app_in->PreStepMeshUserWorkInLoop; } @@ -930,6 +933,10 @@ void Mesh::Initialize(bool init_problem, ParameterInput *pin, ApplicationInput * PARTHENON_REQUIRE_THROWS( !(ProblemGenerator != nullptr && block_list[0]->ProblemGenerator != nullptr), "Mesh and MeshBlock ProblemGenerators are defined. Please use only one."); + PARTHENON_REQUIRE_THROWS( + !(PostInitialization != nullptr && + block_list[0]->PostInitialization != nullptr), + "Mesh and MeshBlock PostInitializations are defined. Please use only one."); // Call Mesh ProblemGenerator if (ProblemGenerator != nullptr) { @@ -946,6 +953,24 @@ void Mesh::Initialize(bool init_problem, ParameterInput *pin, ApplicationInput * pmb->ProblemGenerator(pmb.get(), pin); } } + + // Call Mesh PostInitialization + if (PostInitialization != nullptr) { + PARTHENON_REQUIRE(num_partitions == 1, + "Mesh PostInitialization requires parthenon/mesh/pack_size=-1 " + "during first initialization."); + + auto &md = mesh_data.GetOrAdd("base", 0); + PostInitialization(this, md.get()); + // Call individual MeshBlock PostInitialization + } else { + for (int i = 0; i < nmb; ++i) { + auto &pmb = block_list[i]; + auto &mbd = pmb->meshblock_data.Get(); + pmb->PostInitialization(mbd.get()); + } + } + std::for_each(block_list.begin(), block_list.end(), [](auto &sp_block) { sp_block->SetAllVariablesToInitialized(); }); } diff --git a/src/mesh/mesh.hpp b/src/mesh/mesh.hpp index bf21853d3a82..37d555201598 100644 --- a/src/mesh/mesh.hpp +++ b/src/mesh/mesh.hpp @@ -3,7 +3,7 @@ // Copyright(C) 2014 James M. Stone and other code contributors // Licensed under the 3-clause BSD License, see LICENSE file for details //======================================================================================== -// (C) (or copyright) 2020-2023. Triad National Security, LLC. All rights reserved. +// (C) (or copyright) 2020-2024. Triad National Security, LLC. All rights reserved. // // This program was produced under U.S. Government contract 89233218CNA000001 for Los // Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC @@ -160,6 +160,7 @@ class Mesh { // defined in either the prob file or default_pgen.cpp in ../pgen/ std::function *)> ProblemGenerator = nullptr; + std::function *)> PostInitialization = nullptr; static void UserWorkAfterLoopDefault(Mesh *mesh, ParameterInput *pin, SimTime &tm); // called in main loop std::function UserWorkAfterLoop = diff --git a/src/mesh/meshblock.cpp b/src/mesh/meshblock.cpp index e2d7af7c6ac0..409f29b1cd57 100644 --- a/src/mesh/meshblock.cpp +++ b/src/mesh/meshblock.cpp @@ -3,7 +3,7 @@ // Copyright(C) 2014 James M. Stone and other code contributors // Licensed under the 3-clause BSD License, see LICENSE file for details //======================================================================================== -// (C) (or copyright) 2020-2023. Triad National Security, LLC. All rights reserved. +// (C) (or copyright) 2020-2024. Triad National Security, LLC. All rights reserved. // // This program was produced under U.S. Government contract 89233218CNA000001 for Los // Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC @@ -115,6 +115,12 @@ void MeshBlock::Initialize(int igid, int ilid, LogicalLocation iloc, } else if (app_in->MeshProblemGenerator == nullptr) { ProblemGenerator = &ProblemGeneratorDefault; } + if (app_in->PostInitialization != nullptr) { + PostInitialization = app_in->PostInitialization; + // Only set default post-init when no mesh post-init is set + } else if (app_in->MeshPostInitialization == nullptr) { + PostInitialization = &PostInitializationDefault; + } if (app_in->MeshBlockUserWorkBeforeOutput != nullptr) { UserWorkBeforeOutput = app_in->MeshBlockUserWorkBeforeOutput; } diff --git a/src/mesh/meshblock.hpp b/src/mesh/meshblock.hpp index 3df86d7f73e7..9ad9ca05cae4 100644 --- a/src/mesh/meshblock.hpp +++ b/src/mesh/meshblock.hpp @@ -3,7 +3,7 @@ // Copyright(C) 2014 James M. Stone and other code contributors // Licensed under the 3-clause BSD License, see LICENSE file for details //======================================================================================== -// (C) (or copyright) 2020-2023. Triad National Security, LLC. All rights reserved. +// (C) (or copyright) 2020-2024. Triad National Security, LLC. All rights reserved. // // This program was produced under U.S. Government contract 89233218CNA000001 for Los // Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC @@ -440,7 +440,9 @@ class MeshBlock : public std::enable_shared_from_this { // defined in either the prob file or default_pgen.cpp in ../pgen/ static void ProblemGeneratorDefault(MeshBlock *pmb, ParameterInput *pin); + static void PostInitializationDefault(MeshBlockData *mbd); std::function ProblemGenerator = nullptr; + std::function *)> PostInitialization = nullptr; static pMeshBlockApplicationData_t InitApplicationMeshBlockDataDefault(MeshBlock *, ParameterInput *pin); std::function diff --git a/src/pgen/default_pgen.cpp b/src/pgen/default_pgen.cpp index fc745dd008d8..e8038946a543 100644 --- a/src/pgen/default_pgen.cpp +++ b/src/pgen/default_pgen.cpp @@ -3,7 +3,7 @@ // Copyright(C) 2014 James M. Stone and other code contributors // Licensed under the 3-clause BSD License, see LICENSE file for details //======================================================================================== -// (C) (or copyright) 2020-2022. Triad National Security, LLC. All rights reserved. +// (C) (or copyright) 2020-2024. Triad National Security, LLC. All rights reserved. // // This program was produced under U.S. Government contract 89233218CNA000001 for Los // Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC @@ -110,7 +110,7 @@ void MeshBlock::InitMeshBlockUserDataDefault(MeshBlock *pmb, ParameterInput *pin } //======================================================================================== -//! \fn void MeshBlock::ProblemGeneratorDefault(ParameterInput *pin) +//! \fn void MeshBlock::ProblemGeneratorDefault(MeshBlock *pmb, ParameterInput *pin) // \brief Should be used to set initial conditions. //======================================================================================== @@ -120,6 +120,13 @@ void MeshBlock::ProblemGeneratorDefault(MeshBlock *pmb, ParameterInput *pin) { return; } +//======================================================================================== +//! \fn void MeshBlock::PostInitializationDefault(MeshBlockData *mbd) +// \brief Should be used to perform post initialization ops. +//======================================================================================== + +void MeshBlock::PostInitializationDefault(MeshBlockData *mbd) { return; } + //======================================================================================== //! \fn void MeshBlock::UserWorkBeforeOutputDefault(MeshBlock *pmb, ParameterInput *pin) // \brief Function called before generating output files From f02a89443ea0f08d308587439a0707b970cc4666 Mon Sep 17 00:00:00 2001 From: Patrick Mullen Date: Wed, 24 Jan 2024 18:03:23 -0700 Subject: [PATCH 20/25] Update changelog --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 38e742cf236b..04110c635def 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -3,6 +3,7 @@ ## Current develop ### Added (new features/APIs/variables/...) +- [[PR 999]](https://github.com/parthenon-hpc-lab/parthenon/pull/999) Add a post-initialization hook - [[PR 987]](https://github.com/parthenon-hpc-lab/parthenon/pull/987) New tasking infrastructure and capabilities - [[PR 969]](https://github.com/parthenon-hpc-lab/parthenon/pull/969) New macro-based auto-naming of profiling regions and kernels - [[PR 981]](https://github.com/parthenon-hpc-lab/parthenon/pull/981) Add IndexSplit From 0fc37853d511bcd67ebb9230af019e65d8eef129 Mon Sep 17 00:00:00 2001 From: Patrick Mullen Date: Thu, 25 Jan 2024 14:49:50 -0700 Subject: [PATCH 21/25] Add minimal documentation --- doc/sphinx/src/README.rst | 15 +++++++++------ doc/sphinx/src/parthenon_manager.rst | 13 +++++++++---- 2 files changed, 18 insertions(+), 10 deletions(-) diff --git a/doc/sphinx/src/README.rst b/doc/sphinx/src/README.rst index 521cabf3bc03..96f7258bf33a 100644 --- a/doc/sphinx/src/README.rst +++ b/doc/sphinx/src/README.rst @@ -178,6 +178,8 @@ Mesh ^^^^ - ``InitUserMeshData`` +- ``ProblemGenerator`` +- ``PostInitialization`` - ``PreStepUserWorkInLoop`` - ``PostStepUserWorkInLoop`` - ``UserWorkAfterLoop`` @@ -188,6 +190,7 @@ MeshBlock - ``InitApplicationMeshBlockData`` - ``InitMeshBlockUserData`` - ``ProblemGenerator`` +- ``PostInitialization`` - ``UserWorkBeforeOutput`` To redefine these functions, the user sets the respective function @@ -195,12 +198,12 @@ pointers in the ApplicationInput member app_input of the ParthenonManager class prior to calling ``ParthenonInit``. This is demonstrated in the ``main()`` functions in the examples. -Note that the ``ProblemGenerator``\ s of ``Mesh`` and ``MeshBlock`` are -mutually exclusive. Moreover, the ``Mesh`` one requires -``parthenon/mesh/pack_size=-1`` during initialization, i.e., all blocks -on a rank need to be in a single pack. This allows to use MPI reductions -inside the function, for example, to globally normalize quantities. The -``parthenon/mesh/pack_size=-1`` exists only during problem +Note that the ``ProblemGenerator``\ s (and ``PostInitialization``\ s) of +``Mesh`` and ``MeshBlock`` are mutually exclusive. Moreover, the ``Mesh`` +ones requires ``parthenon/mesh/pack_size=-1`` during initialization, i.e., +all blocks on a rank need to be in a single pack. This allows to use MPI +reductions inside the function, for example, to globally normalize quantities. +The ``parthenon/mesh/pack_size=-1`` exists only during problem inititalization, i.e., simulations can be restarted with an arbitrary ``pack_size``. For an example of the ``Mesh`` version, see the `Poisson example `__. diff --git a/doc/sphinx/src/parthenon_manager.rst b/doc/sphinx/src/parthenon_manager.rst index 15211e8bccdb..e0bac60476d3 100644 --- a/doc/sphinx/src/parthenon_manager.rst +++ b/doc/sphinx/src/parthenon_manager.rst @@ -31,12 +31,16 @@ runtimes. The function Calls the ``Initialize(ParameterInput *pin)`` function of all packages to be utilized and creates the grid hierarchy, including the ``Mesh`` -and ``MeshBlock`` objects, and calls the ``ProblemGenerator`` -initialization routines. +and ``MeshBlock`` objects, and calls the ``ProblemGenerator`` (and +``PostInitialization``) routines. The reason these functions are split out is to enable decisions to be made by the application between reading the input deck and setting up -the grid. For example, a common use-case is: +the grid. For example, during problem initialization, ``ProblemGenerator`` +may be used to be the user-facing API to describe initial conditions, +whereas, ``PostInitialization`` could use those user-specified fields +to sync *all* fields prior to entering communication routines. A common +use-case is: .. code:: cpp @@ -53,13 +57,14 @@ the grid. For example, a common use-case is: if (manager_status == ParthenonStatus::error) { pman.ParthenonFinalize(); return 1; - } + } // Redefine parthenon defaults pman.app_input->ProcessPackages = MyProcessPackages; std::string prob = pman.pin->GetString("app", "problem"); if (prob == "problem1") { pman.app_input->ProblemGenerator = Problem1Generator; + pman.app_input->PostInitialization = Problem1PostInitialization; } else { pman.app_input->ProblemGenerator = Problem2Generator; } From a7276ce423e4504e3a923411989060874f4eca09 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Fri, 26 Jan 2024 13:59:19 -0700 Subject: [PATCH 22/25] fix restart bug. --- src/outputs/output_utils.hpp | 7 ++++--- src/outputs/parthenon_hdf5.cpp | 6 +++--- src/parthenon_manager.cpp | 6 ++++-- 3 files changed, 11 insertions(+), 8 deletions(-) diff --git a/src/outputs/output_utils.hpp b/src/outputs/output_utils.hpp index a768ead308e7..f3f85fbcc6ac 100644 --- a/src/outputs/output_utils.hpp +++ b/src/outputs/output_utils.hpp @@ -214,10 +214,11 @@ std::vector FlattenBlockInfo(Mesh *pm, int shape, Function_t f) { return data; } +// mirror must be provided because copying done externally template -void PackOrUnpackVar(MeshBlock *pmb, Variable *pvar, bool do_ghosts, idx_t &idx, +void PackOrUnpackVar(MeshBlock *pmb, Variable *pvar, + bool do_ghosts, idx_t &idx, std::vector &data, Function_t f) { - auto v_h = pvar->data.GetHostMirrorAndCopy(); const auto &Nt = pvar->GetDim(6); const auto &Nu = pvar->GetDim(5); const auto &Nv = pvar->GetDim(4); @@ -239,7 +240,7 @@ void PackOrUnpackVar(MeshBlock *pmb, Variable *pvar, bool do_ghosts, idx_t for (int k = kb.s; k <= kb.e; ++k) { for (int j = jb.s; j <= jb.e; ++j) { for (int i = ib.s; i <= ib.e; ++i) { - f(v_h, idx, t, u, v, k, j, i); + f(idx, t, u, v, k, j, i); idx++; } } diff --git a/src/outputs/parthenon_hdf5.cpp b/src/outputs/parthenon_hdf5.cpp index 58f80eefdf1b..5e6dc1fa8cf0 100644 --- a/src/outputs/parthenon_hdf5.cpp +++ b/src/outputs/parthenon_hdf5.cpp @@ -396,9 +396,9 @@ void PHDF5Output::WriteOutputFileImpl(Mesh *pm, ParameterInput *pin, SimTime *tm // For reference, if we update the logic here, there's also // a similar block in parthenon_manager.cpp if (v->IsAllocated() && (var_name == v->label())) { - OutputUtils::PackOrUnpackVar( - pmb.get(), v.get(), output_params.include_ghost_zones, index, tmpData, - [&](auto v_h, auto index, int t, int u, int v, int k, int j, int i) { + auto v_h = v->data.GetHostMirrorAndCopy(); + OutputUtils::PackOrUnpackVar(pmb.get(), v.get(), output_params.include_ghost_zones, index, tmpData, + [&](auto index, int t, int u, int v, int k, int j, int i) { tmpData[index] = static_cast(v_h(t, u, v, k, j, i)); }); diff --git a/src/parthenon_manager.cpp b/src/parthenon_manager.cpp index 5cb5d3ce29e9..2d4d063a9ad6 100644 --- a/src/parthenon_manager.cpp +++ b/src/parthenon_manager.cpp @@ -351,9 +351,11 @@ void ParthenonManager::RestartPackages(Mesh &rm, RestartReader &resfile) { } else if (file_output_format_ver == 2 || file_output_format_ver == HDF5::OUTPUT_VERSION_FORMAT) { OutputUtils::PackOrUnpackVar(pmb.get(), v.get(), resfile.hasGhost, index, tmp, - [&](auto v_h, auto index, int t, int u, int v, int k, + [&](auto index, int t, int u, int v, int k, int j, - int i) { v_h(t, u, v, k, j, i) = tmp[index]; }); + int i) { + v_h(t, u, v, k, j, i) = tmp[index]; + }); } else { PARTHENON_THROW("Unknown output format version in restart file.") } From 9abc4012ddeefe94acf3b342baa3135b56c29311 Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Fri, 26 Jan 2024 14:03:05 -0700 Subject: [PATCH 23/25] formatting --- src/outputs/output_utils.hpp | 3 +-- src/outputs/parthenon_hdf5.cpp | 5 +++-- src/parthenon_manager.cpp | 7 ++----- 3 files changed, 6 insertions(+), 9 deletions(-) diff --git a/src/outputs/output_utils.hpp b/src/outputs/output_utils.hpp index f3f85fbcc6ac..db6353090cbf 100644 --- a/src/outputs/output_utils.hpp +++ b/src/outputs/output_utils.hpp @@ -216,8 +216,7 @@ std::vector FlattenBlockInfo(Mesh *pm, int shape, Function_t f) { // mirror must be provided because copying done externally template -void PackOrUnpackVar(MeshBlock *pmb, Variable *pvar, - bool do_ghosts, idx_t &idx, +void PackOrUnpackVar(MeshBlock *pmb, Variable *pvar, bool do_ghosts, idx_t &idx, std::vector &data, Function_t f) { const auto &Nt = pvar->GetDim(6); const auto &Nu = pvar->GetDim(5); diff --git a/src/outputs/parthenon_hdf5.cpp b/src/outputs/parthenon_hdf5.cpp index 5e6dc1fa8cf0..03c14cc540b9 100644 --- a/src/outputs/parthenon_hdf5.cpp +++ b/src/outputs/parthenon_hdf5.cpp @@ -396,8 +396,9 @@ void PHDF5Output::WriteOutputFileImpl(Mesh *pm, ParameterInput *pin, SimTime *tm // For reference, if we update the logic here, there's also // a similar block in parthenon_manager.cpp if (v->IsAllocated() && (var_name == v->label())) { - auto v_h = v->data.GetHostMirrorAndCopy(); - OutputUtils::PackOrUnpackVar(pmb.get(), v.get(), output_params.include_ghost_zones, index, tmpData, + auto v_h = v->data.GetHostMirrorAndCopy(); + OutputUtils::PackOrUnpackVar( + pmb.get(), v.get(), output_params.include_ghost_zones, index, tmpData, [&](auto index, int t, int u, int v, int k, int j, int i) { tmpData[index] = static_cast(v_h(t, u, v, k, j, i)); }); diff --git a/src/parthenon_manager.cpp b/src/parthenon_manager.cpp index 2d4d063a9ad6..6815adfbd2b5 100644 --- a/src/parthenon_manager.cpp +++ b/src/parthenon_manager.cpp @@ -351,11 +351,8 @@ void ParthenonManager::RestartPackages(Mesh &rm, RestartReader &resfile) { } else if (file_output_format_ver == 2 || file_output_format_ver == HDF5::OUTPUT_VERSION_FORMAT) { OutputUtils::PackOrUnpackVar(pmb.get(), v.get(), resfile.hasGhost, index, tmp, - [&](auto index, int t, int u, int v, int k, - int j, - int i) { - v_h(t, u, v, k, j, i) = tmp[index]; - }); + [&](auto index, int t, int u, int v, int k, int j, + int i) { v_h(t, u, v, k, j, i) = tmp[index]; }); } else { PARTHENON_THROW("Unknown output format version in restart file.") } From 1ef878e0727e98483b71b29e52cd37809aadad44 Mon Sep 17 00:00:00 2001 From: Patrick Mullen Date: Tue, 30 Jan 2024 10:54:41 -0700 Subject: [PATCH 24/25] Address @pgrete's comments --- doc/sphinx/src/interface/state.rst | 19 +++++++++++++------ src/application_input.hpp | 6 ++++-- src/mesh/mesh.cpp | 4 ++-- src/mesh/mesh.hpp | 3 ++- src/mesh/meshblock.hpp | 5 +++-- src/pgen/default_pgen.cpp | 7 +++++-- 6 files changed, 29 insertions(+), 15 deletions(-) diff --git a/doc/sphinx/src/interface/state.rst b/doc/sphinx/src/interface/state.rst index e0760d4b8ff1..6a4cfaf840e1 100644 --- a/doc/sphinx/src/interface/state.rst +++ b/doc/sphinx/src/interface/state.rst @@ -113,12 +113,19 @@ several useful features and functions. if set (defaults to ``nullptr`` an therefore a no-op) to print diagnostics after the time-integration advance - ``void UserWorkBeforeLoopMesh(Mesh *, ParameterInput *pin, SimTime - &tm)`` performs a per-package, mesh-wide calculation after the mesh - has been generated, and problem generators called, but before any - time evolution. This work is done both on first initialization and - on restart. If you would like to avoid doing the work upon restart, - you can check for the const ``is_restart`` member field of the ``Mesh`` - object. + &tm)`` performs a per-package, mesh-wide calculation after (1) the mesh + has been generated, (2) problem generators are called, and (3) comms + are executed, but before any time evolution. This work is done both on + first initialization and on restart. If you would like to avoid doing the + work upon restart, you can check for the const ``is_restart`` member + field of the ``Mesh`` object. It is worth making a clear distinction + between ``UserWorkBeforeLoopMesh`` and ``ApplicationInput``s + ``PostInitialization``. ``PostInitialization`` is very much so tied to + initialization, and will not be called upon restarts. ``PostInitialization`` + is also carefully positioned after ``ProblemGenerator`` and before + ``PreCommFillDerived`` (and hence communications). In practice, when + additional granularity is required inbetween initialization and communication, + ``PostInitialization`` may be the desired hook. The reasoning for providing ``FillDerived*`` and ``EstimateTimestep*`` function pointers appropriate for usage with both ``MeshData`` and diff --git a/src/application_input.hpp b/src/application_input.hpp index a7cbaf4fe9c7..bffa9fd44715 100644 --- a/src/application_input.hpp +++ b/src/application_input.hpp @@ -36,7 +36,8 @@ struct ApplicationInput { std::function InitUserMeshData = nullptr; std::function *)> MeshProblemGenerator = nullptr; - std::function *)> MeshPostInitialization = nullptr; + std::function *)> MeshPostInitialization = + nullptr; std::function PreStepMeshUserWorkInLoop = nullptr; @@ -58,7 +59,8 @@ struct ApplicationInput { InitApplicationMeshBlockData = nullptr; std::function InitMeshBlockUserData = nullptr; std::function ProblemGenerator = nullptr; - std::function *)> PostInitialization = nullptr; + std::function *, ParameterInput *)> PostInitialization = + nullptr; std::function MeshBlockUserWorkBeforeOutput = nullptr; }; diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp index 7fbc4b9f8e47..b79fb626718b 100644 --- a/src/mesh/mesh.cpp +++ b/src/mesh/mesh.cpp @@ -961,13 +961,13 @@ void Mesh::Initialize(bool init_problem, ParameterInput *pin, ApplicationInput * "during first initialization."); auto &md = mesh_data.GetOrAdd("base", 0); - PostInitialization(this, md.get()); + PostInitialization(this, pin, md.get()); // Call individual MeshBlock PostInitialization } else { for (int i = 0; i < nmb; ++i) { auto &pmb = block_list[i]; auto &mbd = pmb->meshblock_data.Get(); - pmb->PostInitialization(mbd.get()); + pmb->PostInitialization(mbd.get(), pin); } } diff --git a/src/mesh/mesh.hpp b/src/mesh/mesh.hpp index 37d555201598..48bba457d643 100644 --- a/src/mesh/mesh.hpp +++ b/src/mesh/mesh.hpp @@ -160,7 +160,8 @@ class Mesh { // defined in either the prob file or default_pgen.cpp in ../pgen/ std::function *)> ProblemGenerator = nullptr; - std::function *)> PostInitialization = nullptr; + std::function *)> PostInitialization = + nullptr; static void UserWorkAfterLoopDefault(Mesh *mesh, ParameterInput *pin, SimTime &tm); // called in main loop std::function UserWorkAfterLoop = diff --git a/src/mesh/meshblock.hpp b/src/mesh/meshblock.hpp index 9ad9ca05cae4..b5b6f4df7bd6 100644 --- a/src/mesh/meshblock.hpp +++ b/src/mesh/meshblock.hpp @@ -440,9 +440,10 @@ class MeshBlock : public std::enable_shared_from_this { // defined in either the prob file or default_pgen.cpp in ../pgen/ static void ProblemGeneratorDefault(MeshBlock *pmb, ParameterInput *pin); - static void PostInitializationDefault(MeshBlockData *mbd); + static void PostInitializationDefault(MeshBlockData *mbd, ParameterInput *pin); std::function ProblemGenerator = nullptr; - std::function *)> PostInitialization = nullptr; + std::function *, ParameterInput *)> PostInitialization = + nullptr; static pMeshBlockApplicationData_t InitApplicationMeshBlockDataDefault(MeshBlock *, ParameterInput *pin); std::function diff --git a/src/pgen/default_pgen.cpp b/src/pgen/default_pgen.cpp index e8038946a543..76dc2a530c6f 100644 --- a/src/pgen/default_pgen.cpp +++ b/src/pgen/default_pgen.cpp @@ -121,11 +121,14 @@ void MeshBlock::ProblemGeneratorDefault(MeshBlock *pmb, ParameterInput *pin) { } //======================================================================================== -//! \fn void MeshBlock::PostInitializationDefault(MeshBlockData *mbd) +//! \fn void MeshBlock::PostInitializationDefault(MeshBlockData *mbd, +//! ParameterInput *pin) // \brief Should be used to perform post initialization ops. //======================================================================================== -void MeshBlock::PostInitializationDefault(MeshBlockData *mbd) { return; } +void MeshBlock::PostInitializationDefault(MeshBlockData *mbd, ParameterInput *pin) { + return; +} //======================================================================================== //! \fn void MeshBlock::UserWorkBeforeOutputDefault(MeshBlock *pmb, ParameterInput *pin) From 9bf8094bbd40837821f519d160f54546fe2dffbd Mon Sep 17 00:00:00 2001 From: Patrick Mullen Date: Tue, 30 Jan 2024 11:26:22 -0700 Subject: [PATCH 25/25] Alright... let's just be consistent with ProblemGenerator for now --- src/application_input.hpp | 3 +-- src/mesh/mesh.cpp | 3 +-- src/mesh/meshblock.hpp | 5 ++--- src/pgen/default_pgen.cpp | 7 ++----- 4 files changed, 6 insertions(+), 12 deletions(-) diff --git a/src/application_input.hpp b/src/application_input.hpp index bffa9fd44715..57e1c1962752 100644 --- a/src/application_input.hpp +++ b/src/application_input.hpp @@ -59,8 +59,7 @@ struct ApplicationInput { InitApplicationMeshBlockData = nullptr; std::function InitMeshBlockUserData = nullptr; std::function ProblemGenerator = nullptr; - std::function *, ParameterInput *)> PostInitialization = - nullptr; + std::function PostInitialization = nullptr; std::function MeshBlockUserWorkBeforeOutput = nullptr; }; diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp index b79fb626718b..5ee2f92daf80 100644 --- a/src/mesh/mesh.cpp +++ b/src/mesh/mesh.cpp @@ -966,8 +966,7 @@ void Mesh::Initialize(bool init_problem, ParameterInput *pin, ApplicationInput * } else { for (int i = 0; i < nmb; ++i) { auto &pmb = block_list[i]; - auto &mbd = pmb->meshblock_data.Get(); - pmb->PostInitialization(mbd.get(), pin); + pmb->PostInitialization(pmb.get(), pin); } } diff --git a/src/mesh/meshblock.hpp b/src/mesh/meshblock.hpp index b5b6f4df7bd6..a7b16b9955ac 100644 --- a/src/mesh/meshblock.hpp +++ b/src/mesh/meshblock.hpp @@ -440,10 +440,9 @@ class MeshBlock : public std::enable_shared_from_this { // defined in either the prob file or default_pgen.cpp in ../pgen/ static void ProblemGeneratorDefault(MeshBlock *pmb, ParameterInput *pin); - static void PostInitializationDefault(MeshBlockData *mbd, ParameterInput *pin); + static void PostInitializationDefault(MeshBlock *pmb, ParameterInput *pin); std::function ProblemGenerator = nullptr; - std::function *, ParameterInput *)> PostInitialization = - nullptr; + std::function PostInitialization = nullptr; static pMeshBlockApplicationData_t InitApplicationMeshBlockDataDefault(MeshBlock *, ParameterInput *pin); std::function diff --git a/src/pgen/default_pgen.cpp b/src/pgen/default_pgen.cpp index 76dc2a530c6f..d3c0ed50d164 100644 --- a/src/pgen/default_pgen.cpp +++ b/src/pgen/default_pgen.cpp @@ -121,14 +121,11 @@ void MeshBlock::ProblemGeneratorDefault(MeshBlock *pmb, ParameterInput *pin) { } //======================================================================================== -//! \fn void MeshBlock::PostInitializationDefault(MeshBlockData *mbd, -//! ParameterInput *pin) +//! \fn void MeshBlock::PostInitializationDefault(MeshBlock *pmb, ParameterInput *pin) // \brief Should be used to perform post initialization ops. //======================================================================================== -void MeshBlock::PostInitializationDefault(MeshBlockData *mbd, ParameterInput *pin) { - return; -} +void MeshBlock::PostInitializationDefault(MeshBlock *pmb, ParameterInput *pin) { return; } //======================================================================================== //! \fn void MeshBlock::UserWorkBeforeOutputDefault(MeshBlock *pmb, ParameterInput *pin)