diff --git a/CHANGELOG.md b/CHANGELOG.md index c37f63c6d2d3..69a5b4fc816a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -3,11 +3,11 @@ ## Current develop ### Added (new features/APIs/variables/...) -- [[PR 1171]](https://github.com/parthenon-hpc-lab/parthenon/pull/1171) Add PARTHENON_USE_SYSTEM_PACKAGES build option - +- [[PR 1161]](https://github.com/parthenon-hpc-lab/parthenon/pull/1161) Make flux field Metadata accessible, add Metadata::CellMemAligned flag, small perfomance upgrades ### Changed (changing behavior/API/variables/...) -- [[PR1172]](https://github.com/parthenon-hpc-lab/parthenon/pull/1172) Make parthenon manager robust against external MPI init and finalize calls +- [[PR 1171]](https://github.com/parthenon-hpc-lab/parthenon/pull/1171) Add PARTHENON_USE_SYSTEM_PACKAGES build option +- [[PR 1172]](https://github.com/parthenon-hpc-lab/parthenon/pull/1172) Make parthenon manager robust against external MPI init and finalize calls ### Fixed (not changing behavior/API/variables/...) diff --git a/doc/sphinx/src/interface/metadata.rst b/doc/sphinx/src/interface/metadata.rst index 4f08b7ede310..0ccfe73fe4a2 100644 --- a/doc/sphinx/src/interface/metadata.rst +++ b/doc/sphinx/src/interface/metadata.rst @@ -158,9 +158,22 @@ classes may be allocated. The behaviours are the following: sense (e.g. if the ``WithFluxes`` variable has ``Metadata::Cell`` set the new variable will have ``Metadata::Face``) will be created in the package with the name ``bnd_flux::name_of_original_variable`` and - ``Metadata::Flux`` and ``Metadata::OneCopy``. When creating packs that - include fluxes, the new flux field will be included in the flux portion - of the pack if the parent field is in the pack. + ``Metadata::Flux`` and ``Metadata::OneCopy``. Additionally, the flags + ``Metadata::Sparse``, ``Metadata::Vector``, and ``Metadata::Tensor`` + will propagate to the flux ``Metadata`` if they are set in the base field + ``Metadata``. By default, a flux field for a cell-centered field is + built with ``Metadata::CellMemAligned`` flag set for backwards + compatability. A shared pointer to the ``Metadata`` object for the flux + field can be accessed from the base ``Metadata`` with the method + ``Metadata::GetSPtrFluxMetadata()``. This can be used to set flags other + than the defaults or set custom prolongation/restriction operations for + the fluxes. Note that calling `Metadata::RegisterRefinementOps<...>()` + on the base field propagates the registered refinement operations through + to the flux `Metadata` for backward compatibility. If separate operations + are desired for the fluxes, the ordering of calls to `RegisterRefinementOps` + on the base field and the flux field matters. When creating packs that + include fluxes, the new flux field will be included in the flux portion of + the pack if the parent field is in the pack. - If ``Metadata::Flux`` is set, this field is exchanged on shared elements across fine-coarse boundaries when the flux correction tasks are called. diff --git a/src/bvals/comms/boundary_communication.cpp b/src/bvals/comms/boundary_communication.cpp index 1f238b1fdbf7..78121cd3fac9 100644 --- a/src/bvals/comms/boundary_communication.cpp +++ b/src/bvals/comms/boundary_communication.cpp @@ -213,7 +213,7 @@ TaskStatus ReceiveBoundBufs(std::shared_ptr> &md) { [&all_received](auto pbuf) { all_received = pbuf->TryReceive() && all_received; }); int ibound = 0; - if (Globals::sparse_config.enabled) { + if (Globals::sparse_config.enabled && all_received) { ForEachBoundary( md, [&](auto pmb, sp_mbd_t rc, nb_t &nb, const sp_cv_t v) { const std::size_t ibuf = cache.idx_vec[ibound]; diff --git a/src/interface/meshblock_data.hpp b/src/interface/meshblock_data.hpp index c7133317bd76..9fab1f354674 100644 --- a/src/interface/meshblock_data.hpp +++ b/src/interface/meshblock_data.hpp @@ -77,6 +77,9 @@ class MeshBlockData { void SetAllowedDt(const Real dt) const { GetBlockPointer()->SetAllowedDt(dt); } Mesh *GetMeshPointer() const { return GetBlockPointer()->pmy_mesh; } + // This mirrors a MeshBlockData routine + int NumBlocks() const { return 1; } + template IndexRange GetBoundsI(Ts &&...args) const { return GetBlockPointer()->cellbounds.GetBoundsI(std::forward(args)...); diff --git a/src/interface/metadata.cpp b/src/interface/metadata.cpp index fce88334c652..36f37adf8b1c 100644 --- a/src/interface/metadata.cpp +++ b/src/interface/metadata.cpp @@ -101,10 +101,13 @@ MetadataFlag Metadata::GetUserFlag(const std::string &flagname) { } namespace parthenon { -Metadata::Metadata(const std::vector &bits, const std::vector &shape, +Metadata::Metadata(const std::vector &bits, + const std::vector &flux_bits, + const std::vector &shape, const std::vector &component_labels, const std::string &associated, - const refinement::RefinementFunctions_t ref_funcs_) + const refinement::RefinementFunctions_t ref_funcs_, + const refinement::RefinementFunctions_t flux_ref_funcs_) : shape_(shape), component_labels_(component_labels), associated_(associated) { // set flags for (const auto f : bits) { @@ -164,6 +167,39 @@ Metadata::Metadata(const std::vector &bits, const std::vector deallocation_threshold_ = 0.0; default_value_ = 0.0; } + + // Now create the flux metadata if required + if (IsSet(WithFluxes)) { + std::set flux_flags; + for (const auto f : flux_bits) + flux_flags.insert(f); + + // Set some standard defaults for the flux metadata if no + // flags were provided + if (flux_flags.size() == 0) { + flux_flags.insert(OneCopy); + if (IsSet(Fine)) flux_flags.insert(Fine); + if (IsSet(Cell)) flux_flags.insert(CellMemAligned); + if (IsSet(Sparse)) flux_flags.insert(Sparse); + } + + // These flags are automatically propagated for fluxes + flux_flags.insert(Flux); + if (IsSet(Cell)) { + flux_flags.insert(Face); + } else if (IsSet(Face)) { + flux_flags.insert(Edge); + } else if (IsSet(Edge)) { + flux_flags.insert(Node); + } + + if (IsSet(Tensor)) flux_flags.insert(Tensor); + if (IsSet(Vector)) flux_flags.insert(Vector); + + std::vector flux_flags_vec(flux_flags.begin(), flux_flags.end()); + flux_metadata = std::make_shared(flux_flags_vec, shape, component_labels, + std::string(), flux_ref_funcs_); + } } std::ostream &operator<<(std::ostream &os, const parthenon::Metadata &m) { @@ -252,6 +288,14 @@ bool Metadata::IsValid(bool throw_on_fail) const { } } + if (IsSet(FillGhost) && IsSet(CellMemAligned) && (!IsSet(Cell))) { + valid = false; + if (throw_on_fail) { + PARTHENON_THROW( + "Cannot communicate ghosts of non-cell fields that have cell aligned memory."); + } + } + // Prolongation/restriction if (HasRefinementOps()) { if (refinement_funcs_.label().size() == 0) { @@ -306,16 +350,12 @@ Metadata::GetArrayDims(std::weak_ptr wpmb, bool coarse) const { arrDims[i + 3] = 1; if (IsSet(Cell)) { arrDims[MAX_VARIABLE_DIMENSION - 1] = 1; // Only one cell center per cell - } else if (IsSet(Face) && IsSet(Flux)) { - // 3 directions but keep the same ijk shape as cell var for performance - arrDims[MAX_VARIABLE_DIMENSION - 1] = 3; } else if (IsSet(Face) || IsSet(Edge)) { arrDims[MAX_VARIABLE_DIMENSION - 1] = 3; // Three faces and edges per cell - arrDims[0]++; - if (arrDims[1] > 1) arrDims[1]++; - if (arrDims[2] > 1) arrDims[2]++; } else if (IsSet(Node)) { arrDims[MAX_VARIABLE_DIMENSION - 1] = 1; // Only one lower left node per cell + } + if (!IsSet(CellMemAligned) && !IsSet(Cell)) { arrDims[0]++; if (arrDims[1] > 1) arrDims[1]++; if (arrDims[2] > 1) arrDims[2]++; diff --git a/src/interface/metadata.hpp b/src/interface/metadata.hpp index 5770323c21ba..467038040ed7 100644 --- a/src/interface/metadata.hpp +++ b/src/interface/metadata.hpp @@ -120,6 +120,9 @@ PARTHENON_INTERNAL_FOR_FLAG(Fine) \ /** this variable is the flux for another variable **/ \ PARTHENON_INTERNAL_FOR_FLAG(Flux) \ + /** Align memory of fields to cell centered memory \ + (Field will be missing one layer of ghosts if it is not cell centered) **/ \ + PARTHENON_INTERNAL_FOR_FLAG(CellMemAligned) \ /************************************************/ \ /** Vars specifying coordinates for visualization purposes **/ \ /** You can specify a single 3D var **/ \ @@ -325,28 +328,48 @@ class Metadata { // 4 constructors, this is the general constructor called by all other constructors, so // we do some sanity checks here Metadata( - const std::vector &bits, const std::vector &shape = {}, + const std::vector &bits, const std::vector &flux_bits, + const std::vector &shape = {}, const std::vector &component_labels = {}, const std::string &associated = "", const refinement::RefinementFunctions_t ref_funcs_ = + refinement::RefinementFunctions_t::RegisterOps< + refinement_ops::ProlongateSharedMinMod, refinement_ops::RestrictAverage>(), + const refinement::RefinementFunctions_t flux_ref_funcs_ = refinement::RefinementFunctions_t::RegisterOps< refinement_ops::ProlongateSharedMinMod, refinement_ops::RestrictAverage>()); - // 1 constructor + Metadata( + const std::vector &bits, const std::vector &shape = {}, + const std::vector &component_labels = {}, + const std::string &associated = "", + const refinement::RefinementFunctions_t ref_funcs_ = + refinement::RefinementFunctions_t::RegisterOps< + refinement_ops::ProlongateSharedMinMod, refinement_ops::RestrictAverage>()) + : Metadata(bits, {}, shape, component_labels, associated, ref_funcs_, ref_funcs_) {} + Metadata(const std::vector &bits, const std::vector &shape, const std::string &associated) : Metadata(bits, shape, {}, associated) {} - // 2 constructors Metadata(const std::vector &bits, const std::vector component_labels, const std::string &associated = "") : Metadata(bits, {1}, component_labels, associated) {} - // 1 constructor Metadata(const std::vector &bits, const std::string &associated) : Metadata(bits, {1}, {}, associated) {} + std::shared_ptr GetSPtrFluxMetadata() { + PARTHENON_REQUIRE(IsSet(WithFluxes), + "Asking for flux metadata from metadata that doesn't have it."); + return flux_metadata; + } + + private: + std::shared_ptr flux_metadata; + + public: // Static routines static MetadataFlag AddUserFlag(const std::string &name); static bool FlagNameExists(const std::string &flagname); @@ -536,6 +559,12 @@ class Metadata { refinement_funcs_ = refinement::RefinementFunctions_t::RegisterOps(); + // Propagate refinement operations to flux metadata for backward compatibility + if (IsSet(WithFluxes)) { + flux_metadata->refinement_funcs_ = + refinement::RefinementFunctions_t::RegisterOps(); + } } // Operators @@ -688,7 +717,6 @@ Set_t GetByFlag(const Metadata::FlagCollection &flags, NameMap_t &nameMap, return out; } } // namespace MetadataUtils - } // namespace parthenon #endif // INTERFACE_METADATA_HPP_ diff --git a/src/interface/sparse_pack_base.cpp b/src/interface/sparse_pack_base.cpp index d5a66db01498..d4ead84113c8 100644 --- a/src/interface/sparse_pack_base.cpp +++ b/src/interface/sparse_pack_base.cpp @@ -67,19 +67,16 @@ SparsePackBase::GetAllocStatus(T *pmd, const PackDescriptor &desc, const std::vector &include_block) { using mbd_t = MeshBlockData; - int nvar = desc.nvar_groups; - - std::vector astat; + const int nvar = desc.nvar_groups; + const int nblock = pmd->NumBlocks(); + std::vector astat(nblock * desc.nvar_tot); + int idx = 0; ForEachBlock(pmd, include_block, [&](int b, mbd_t *pmbd) { const auto &uid_map = pmbd->GetUidMap(); for (int i = 0; i < nvar; ++i) { for (const auto &[var_name, uid] : desc.var_groups[i]) { - if (uid_map.count(uid) > 0) { - const auto pv = uid_map.at(uid); - astat.push_back(pv->GetAllocationStatus()); - } else { - astat.push_back(-1); - } + astat[idx++] = + uid_map.count(uid) > 0 ? (uid_map.at(uid))->GetAllocationStatus() : -1; } } }); diff --git a/src/interface/sparse_pack_base.hpp b/src/interface/sparse_pack_base.hpp index 53fc4e37d0b1..4b01547d7b04 100644 --- a/src/interface/sparse_pack_base.hpp +++ b/src/interface/sparse_pack_base.hpp @@ -138,7 +138,7 @@ struct PackDescriptor { // default constructor needed for certain use cases PackDescriptor() : nvar_groups(0), var_group_names({}), var_groups({}), with_fluxes(false), - coarse(false), flat(false), identifier("") {} + coarse(false), flat(false), identifier(""), nvar_tot(0) {} template PackDescriptor(StateDescriptor *psd, const std::vector &var_groups_in, @@ -147,7 +147,7 @@ struct PackDescriptor { var_groups(BuildUids(var_groups_in.size(), psd, selector)), with_fluxes(options.count(PDOpt::WithFluxes)), coarse(options.count(PDOpt::Coarse)), flat(options.count(PDOpt::Flatten)), - identifier(GetIdentifier()) { + identifier(GetIdentifier()), nvar_tot(GetNVarsTotal(var_groups)) { PARTHENON_REQUIRE(!(with_fluxes && coarse), "Probably shouldn't be making a coarse pack with fine fluxes."); } @@ -159,8 +159,18 @@ struct PackDescriptor { const bool coarse; const bool flat; const std::string identifier; + const std::size_t nvar_tot; private: + static int GetNVarsTotal(const std::vector &var_groups) { + int nvar_tot = 0; + for (const auto &group : var_groups) { + for (const auto &[a, b] : group) { + nvar_tot++; + } + } + return nvar_tot; + } std::string GetIdentifier() { std::string ident(""); for (const auto &vgroup : var_groups) { diff --git a/src/interface/sparse_pool.cpp b/src/interface/sparse_pool.cpp index 8a6e0ef213ca..7890cff964de 100644 --- a/src/interface/sparse_pool.cpp +++ b/src/interface/sparse_pool.cpp @@ -11,6 +11,8 @@ // the public, perform publicly and display publicly, and to permit others to do so. //======================================================================================== +#include + #include "interface/sparse_pool.hpp" #include "interface/metadata.hpp" @@ -48,34 +50,31 @@ SparsePool::SparsePool(const std::string &base_name, const Metadata &metadata, } } -const Metadata &SparsePool::AddImpl(int sparse_id, const std::vector &shape, - const MetadataFlag *vector_tensor, - const std::vector &component_labels) { - PARTHENON_REQUIRE_THROWS(sparse_id != InvalidSparseID, - "Tried to add InvalidSparseID to sparse pool " + base_name_); - +std::shared_ptr +MakeSparseVarMetadataImpl(Metadata *in, const std::vector &shape, + const MetadataFlag *vector_tensor, + const std::vector &component_labels) { // copy shared metadata - Metadata this_metadata( - shared_metadata_.Flags(), shape.size() > 0 ? shape : shared_metadata_.Shape(), - component_labels.size() > 0 ? component_labels - : shared_metadata_.getComponentLabels(), - shared_metadata_.getAssociated(), shared_metadata_.GetRefinementFunctions()); + auto this_metadata = std::make_shared( + in->Flags(), shape.size() > 0 ? shape : in->Shape(), + component_labels.size() > 0 ? component_labels : in->getComponentLabels(), + in->getAssociated(), in->GetRefinementFunctions()); - this_metadata.SetSparseThresholds(shared_metadata_.GetAllocationThreshold(), - shared_metadata_.GetDeallocationThreshold(), - shared_metadata_.GetDefaultValue()); + this_metadata->SetSparseThresholds(in->GetAllocationThreshold(), + in->GetDeallocationThreshold(), + in->GetDefaultValue()); // if vector_tensor is set, apply it if (vector_tensor != nullptr) { if (*vector_tensor == Metadata::Vector) { - this_metadata.Unset(Metadata::Tensor); - this_metadata.Set(Metadata::Vector); + this_metadata->Unset(Metadata::Tensor); + this_metadata->Set(Metadata::Vector); } else if (*vector_tensor == Metadata::Tensor) { - this_metadata.Unset(Metadata::Vector); - this_metadata.Set(Metadata::Tensor); + this_metadata->Unset(Metadata::Vector); + this_metadata->Set(Metadata::Tensor); } else if (*vector_tensor == Metadata::None) { - this_metadata.Unset(Metadata::Vector); - this_metadata.Unset(Metadata::Tensor); + this_metadata->Unset(Metadata::Vector); + this_metadata->Unset(Metadata::Tensor); } else { PARTHENON_THROW("Expected MetadataFlag Vector, Tensor, or None, but got " + vector_tensor->Name()); @@ -83,9 +82,26 @@ const Metadata &SparsePool::AddImpl(int sparse_id, const std::vector &shape } // just in case - this_metadata.IsValid(true); + this_metadata->IsValid(true); + + return this_metadata; +} + +const Metadata &SparsePool::AddImpl(int sparse_id, const std::vector &shape, + const MetadataFlag *vector_tensor, + const std::vector &component_labels) { + PARTHENON_REQUIRE_THROWS(sparse_id != InvalidSparseID, + "Tried to add InvalidSparseID to sparse pool " + base_name_); + + auto this_metadata = MakeSparseVarMetadataImpl(&shared_metadata_, shape, vector_tensor, + component_labels); + if (this_metadata->IsSet(Metadata::WithFluxes)) { + this_metadata->GetSPtrFluxMetadata() = + MakeSparseVarMetadataImpl(shared_metadata_.GetSPtrFluxMetadata().get(), shape, + vector_tensor, component_labels); + } - const auto ins = pool_.insert({sparse_id, this_metadata}); + const auto ins = pool_.insert({sparse_id, *this_metadata}); PARTHENON_REQUIRE_THROWS(ins.second, "Tried to add sparse ID " + std::to_string(sparse_id) + " to sparse pool '" + base_name_ + diff --git a/src/interface/state_descriptor.cpp b/src/interface/state_descriptor.cpp index af73259ff6aa..b468b0201091 100644 --- a/src/interface/state_descriptor.cpp +++ b/src/interface/state_descriptor.cpp @@ -274,27 +274,9 @@ bool StateDescriptor::AddFieldImpl(const VarID &vid, const Metadata &m_in, return false; // this field has already been added } else { if (m.IsSet(Metadata::WithFluxes) && m.GetFluxName() == "") { - std::vector mFlags = {Metadata::OneCopy, Metadata::Flux}; - if (m.IsSet(Metadata::Sparse)) mFlags.push_back(Metadata::Sparse); - if (m.IsSet(Metadata::Fine)) mFlags.push_back(Metadata::Fine); - if (m.IsSet(Metadata::Cell)) - mFlags.push_back(Metadata::Face); - else if (m.IsSet(Metadata::Face)) - mFlags.push_back(Metadata::Edge); - else if (m.IsSet(Metadata::Edge)) - mFlags.push_back(Metadata::Node); - - Metadata mf; - if (m.GetRefinementFunctions().label().size() > 0) { - // Propagate custom refinement ops to flux field - mf = Metadata(mFlags, m.Shape(), std::vector(), std::string(), - m.GetRefinementFunctions()); - } else { - mf = Metadata(mFlags, m.Shape()); - } auto fId = VarID{internal_fluxname + internal_varname_seperator + vid.base_name, vid.sparse_id}; - AddFieldImpl(fId, mf, control_vid); + AddFieldImpl(fId, *(m.GetSPtrFluxMetadata()), control_vid); m.SetFluxName(fId.label()); } metadataMap_.insert({vid, m});