diff --git a/CHANGELOG.md b/CHANGELOG.md index 78fa0bfd3c01..c97cf4c57d08 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -26,6 +26,7 @@ - [[PR 1004]](https://github.com/parthenon-hpc-lab/parthenon/pull/1004) Allow parameter modification from an input file for restarts ### Fixed (not changing behavior/API/variables/...) +- [[PR 1092]](https://github.com/parthenon-hpc-lab/parthenon/pull/1092) Updates to DataCollection and MeshData to remove requirement of predefining MeshBlockData - [[PR 1113]](https://github.com/parthenon-hpc-lab/parthenon/pull/1113) Prevent division by zero - [[PR 1112]](https://github.com/parthenon-hpc-lab/parthenon/pull/1112) Remove shared_ptr cycle in forest::Tree - [[PR 1104]](https://github.com/parthenon-hpc-lab/parthenon/pull/1104) Fix reading restarts due to hidden ghost var diff --git a/example/fine_advection/advection_driver.cpp b/example/fine_advection/advection_driver.cpp index f755eac8ab61..d7094663a508 100644 --- a/example/fine_advection/advection_driver.cpp +++ b/example/fine_advection/advection_driver.cpp @@ -62,17 +62,6 @@ TaskCollection AdvectionDriver::MakeTaskCollection(BlockList_t &blocks, const in // Build MeshBlockData containers that will be included in MeshData containers. It is // gross that this has to be done by hand. const auto &stage_name = integrator->stage_name; - if (stage == 1) { - for (int i = 0; i < blocks.size(); i++) { - auto &pmb = blocks[i]; - // first make other useful containers - auto &base = pmb->meshblock_data.Get(); - pmb->meshblock_data.Add("dUdt", base); - for (int s = 1; s < integrator->nstages; s++) - pmb->meshblock_data.Add(stage_name[s], base); - } - } - const Real beta = integrator->beta[stage - 1]; const Real dt = integrator->dt; diff --git a/src/interface/data_collection.cpp b/src/interface/data_collection.cpp index 81e6886e29e4..20ad1b9a2836 100644 --- a/src/interface/data_collection.cpp +++ b/src/interface/data_collection.cpp @@ -32,30 +32,26 @@ std::shared_ptr &DataCollection::Add(const std::string &label) { return containers_[label]; } +template <> std::shared_ptr> & -GetOrAdd_impl(Mesh *pmy_mesh_, - std::map>> &containers_, - BlockList_t &block_list, const std::string &mbd_label, - const int &partition_id, const std::optional gmg_level) { - std::string label = mbd_label + "_part-" + std::to_string(partition_id); - if (gmg_level) label = label + "_gmg-" + std::to_string(*gmg_level); +DataCollection>::GetOrAdd_impl(const std::string &mbd_label, + const int &partition_id, + const std::optional gmg_level) { + std::string label = GetKey(mbd_label, partition_id, gmg_level); auto it = containers_.find(label); if (it == containers_.end()) { // TODO(someone) add caching of partitions to Mesh at some point const int pack_size = pmy_mesh_->DefaultPackSize(); + auto &block_list = + gmg_level ? pmy_mesh_->gmg_block_lists[*gmg_level] : pmy_mesh_->block_list; auto partitions = partition::ToSizeN(block_list, pack_size); // Account for possibly empty block_list if (partitions.size() == 0) partitions = std::vector(1); for (auto i = 0; i < partitions.size(); i++) { - std::string md_label = mbd_label + "_part-" + std::to_string(i); - if (gmg_level) md_label = md_label + "_gmg-" + std::to_string(*gmg_level); + std::string md_label = GetKey(mbd_label, partition_id, gmg_level); containers_[md_label] = std::make_shared>(mbd_label); - containers_[md_label]->Set(partitions[i], pmy_mesh_); - if (gmg_level) { - containers_[md_label]->grid = GridIdentifier::two_level_composite(*gmg_level); - } else { - containers_[md_label]->grid = GridIdentifier::leaf(); - } + containers_[md_label]->Initialize(partitions[i], pmy_mesh_, gmg_level); + containers_[md_label]->partition = i; } } return containers_[label]; @@ -65,16 +61,14 @@ template <> std::shared_ptr> & DataCollection>::GetOrAdd(const std::string &mbd_label, const int &partition_id) { - return GetOrAdd_impl(pmy_mesh_, containers_, pmy_mesh_->block_list, mbd_label, - partition_id, {}); + return GetOrAdd_impl(mbd_label, partition_id, {}); } template <> std::shared_ptr> & DataCollection>::GetOrAdd(int gmg_level, const std::string &mbd_label, const int &partition_id) { - return GetOrAdd_impl(pmy_mesh_, containers_, pmy_mesh_->gmg_block_lists[gmg_level], - mbd_label, partition_id, gmg_level); + return GetOrAdd_impl(mbd_label, partition_id, gmg_level); } template class DataCollection>; diff --git a/src/interface/data_collection.hpp b/src/interface/data_collection.hpp index 913fc8439015..29551d35290c 100644 --- a/src/interface/data_collection.hpp +++ b/src/interface/data_collection.hpp @@ -19,10 +19,17 @@ #include #include +#include "basic_types.hpp" +#include "utils/concepts_lite.hpp" #include "utils/error_checking.hpp" namespace parthenon { class Mesh; +class MeshBlock; +template +class MeshData; +template +class MeshBlockData; /// The DataCollection class is an abstract container that contains at least a /// "base" container of some type (e.g., of MeshData or MeshBlockData) plus /// additional containers identified by string labels. @@ -44,31 +51,39 @@ class DataCollection { void SetMeshPointer(Mesh *pmesh) { pmy_mesh_ = pmesh; } - template - std::shared_ptr &Add(const std::string &name, const std::shared_ptr &src, + template + std::shared_ptr &Add(const std::string &name, const std::shared_ptr &src, const std::vector &fields, const bool shallow) { - auto it = containers_.find(name); + if constexpr (!((std::is_same_v && + std::is_same_v>) || + std::is_same_v)) { + // SRC_t and T are incompatible + static_assert(always_false, "Incompatible source and container types."); + } + + auto key = GetKey(name, src); + auto it = containers_.find(key); if (it != containers_.end()) { if (fields.size() && !(it->second)->Contains(fields)) { - PARTHENON_THROW(name + " already exists in collection but fields do not match."); + PARTHENON_THROW(key + " already exists in collection but fields do not match."); } return it->second; } auto c = std::make_shared(name); - c->Initialize(src.get(), fields, shallow); + c->Initialize(src, fields, shallow); - Set(name, c); - - return containers_[name]; + containers_[key] = c; + return containers_[key]; } - template - std::shared_ptr &Add(const std::string &label, const std::shared_ptr &src, + template + std::shared_ptr &Add(const std::string &label, const std::shared_ptr &src, const std::vector &fields = {}) { return Add(label, src, fields, false); } - template - std::shared_ptr &AddShallow(const std::string &label, const std::shared_ptr &src, + template + std::shared_ptr &AddShallow(const std::string &label, + const std::shared_ptr &src, const std::vector &fields = {}) { return Add(label, src, fields, true); } @@ -105,6 +120,28 @@ class DataCollection { } private: + template + std::string GetKey(const std::string &stage_label, const std::shared_ptr &in) { + if constexpr (std::is_same_v>) { + std::string key = stage_label + "_part-" + std::to_string(in->partition); + if (in->grid.type == GridType::two_level_composite) + key = key + "_gmg-" + std::to_string(in->grid.logical_level); + return key; + } else { + return stage_label; + } + } + std::string GetKey(const std::string &stage_label, std::optional partition_id, + std::optional gmg_level) { + std::string key = stage_label; + if (partition_id) key = key + "_part-" + std::to_string(*partition_id); + if (gmg_level) key = key + "_gmg-" + std::to_string(*gmg_level); + return key; + } + + std::shared_ptr &GetOrAdd_impl(const std::string &mbd_label, const int &partition_id, + const std::optional gmg_level); + Mesh *pmy_mesh_; std::map> containers_; }; diff --git a/src/interface/mesh_data.cpp b/src/interface/mesh_data.cpp index e463c58219a0..bc1a2e91556a 100644 --- a/src/interface/mesh_data.cpp +++ b/src/interface/mesh_data.cpp @@ -17,23 +17,30 @@ namespace parthenon { template -void MeshData::Set(BlockList_t blocks, Mesh *pmesh, int ndim) { +void MeshData::Initialize(BlockList_t blocks, Mesh *pmesh, int ndim, + std::optional gmg_level) { const int nblocks = blocks.size(); ndim_ = ndim; block_data_.resize(nblocks); SetMeshPointer(pmesh); for (int i = 0; i < nblocks; i++) { - block_data_[i] = blocks[i]->meshblock_data.Get(stage_name_); + block_data_[i] = blocks[i]->meshblock_data.Add(stage_name_, blocks[i]); + } + if (gmg_level) { + grid = GridIdentifier::two_level_composite(*gmg_level); + } else { + grid = GridIdentifier::leaf(); } } template -void MeshData::Set(BlockList_t blocks, Mesh *pmesh) { +void MeshData::Initialize(BlockList_t blocks, Mesh *pmesh, + std::optional gmg_level) { int ndim; if (pmesh != nullptr) { ndim = pmesh->ndim; } - Set(blocks, pmesh, ndim); + Initialize(blocks, pmesh, ndim, gmg_level); } template diff --git a/src/interface/mesh_data.hpp b/src/interface/mesh_data.hpp index cb1a24c187d5..4bde8c311017 100644 --- a/src/interface/mesh_data.hpp +++ b/src/interface/mesh_data.hpp @@ -190,10 +190,12 @@ const MeshBlockPack

&PackOnMesh(M &map, BlockDataList_t &block_data_, template class MeshData { public: + using parent_t = Mesh; MeshData() = default; explicit MeshData(const std::string &name) : stage_name_(name) {} GridIdentifier grid; + int partition; const auto &StageName() const { return stage_name_; } @@ -241,11 +243,12 @@ class MeshData { } } - void Set(BlockList_t blocks, Mesh *pmesh, int ndim); - void Set(BlockList_t blocks, Mesh *pmesh); + void Initialize(BlockList_t blocks, Mesh *pmesh, int ndim, + std::optional gmg_level = {}); + void Initialize(BlockList_t blocks, Mesh *pmesh, std::optional gmg_level = {}); template - void Initialize(const MeshData *src, const std::vector &vars, + void Initialize(std::shared_ptr> src, const std::vector &vars, const bool shallow) { if (src == nullptr) { PARTHENON_THROW("src points at null"); diff --git a/src/interface/meshblock_data.cpp b/src/interface/meshblock_data.cpp index 3017d5c419d4..9a927062f8aa 100644 --- a/src/interface/meshblock_data.cpp +++ b/src/interface/meshblock_data.cpp @@ -33,37 +33,6 @@ #include "utils/utils.hpp" namespace parthenon { - -template -void MeshBlockData::Initialize( - const std::shared_ptr resolved_packages, - const std::shared_ptr pmb) { - SetBlockPointer(pmb); - resolved_packages_ = resolved_packages; - - // clear all variables, maps, and pack caches - varVector_.clear(); - varMap_.clear(); - varUidMap_.clear(); - flagsToVars_.clear(); - varPackMap_.clear(); - coarseVarPackMap_.clear(); - varFluxPackMap_.clear(); - - for (auto const &q : resolved_packages->AllFields()) { - AddField(q.first.base_name, q.second, q.first.sparse_id); - } - - const auto &swarm_container = GetSwarmData(); - swarm_container->Initialize(resolved_packages, pmb); - - Metadata::FlagCollection flags({Metadata::Sparse, Metadata::ForceAllocOnNewBlocks}); - auto vars = GetVariablesByFlag(flags); - for (auto &v : vars.vars()) { - AllocateSparse(v->label()); - } -} - /// /// The internal routine for adding a new field. This subroutine /// is topology aware and will allocate accordingly. @@ -270,9 +239,9 @@ MeshBlockData::GetVariablesByName(const std::vector &names, } else { var_list.Add(v, sparse_ids_set); } - } else if ((resolved_packages_ != nullptr) && - (resolved_packages_->SparseBaseNamePresent(name))) { - const auto &sparse_pool = resolved_packages_->GetSparsePool(name); + } else if ((resolved_packages != nullptr) && + (resolved_packages->SparseBaseNamePresent(name))) { + const auto &sparse_pool = resolved_packages->GetSparsePool(name); // add all sparse ids of the pool for (const auto iter : sparse_pool.pool()) { diff --git a/src/interface/meshblock_data.hpp b/src/interface/meshblock_data.hpp index 206f9beaf75a..a993bbb07208 100644 --- a/src/interface/meshblock_data.hpp +++ b/src/interface/meshblock_data.hpp @@ -31,11 +31,11 @@ #include "interface/variable.hpp" #include "interface/variable_pack.hpp" #include "mesh/domain.hpp" +#include "utils/concepts_lite.hpp" #include "utils/error_checking.hpp" #include "utils/unique_id.hpp" namespace parthenon { - /// Interface to underlying infrastructure for data declaration and access. /// /// The MeshBlockData class is a container for the variables that make up @@ -66,13 +66,6 @@ class MeshBlockData { MeshBlockData() = default; explicit MeshBlockData(const std::string &name) : stage_name_(name) {} - // Constructors for getting sub-containers - // the variables returned are all shallow copies of the src container. - MeshBlockData(const MeshBlockData &src, const std::vector &names, - const std::vector &sparse_ids = {}); - MeshBlockData(const MeshBlockData &src, const std::vector &flags, - const std::vector &sparse_ids = {}); - std::shared_ptr GetBlockSharedPointer() const { if (pmy_block.expired()) { PARTHENON_THROW("Invalid pointer to MeshBlock!"); @@ -122,18 +115,27 @@ class MeshBlockData { pmy_block = other->GetBlockSharedPointer(); } - void Initialize(const std::shared_ptr resolved_packages, - const std::shared_ptr pmb); - /// Create copy of MeshBlockData, possibly with a subset of named fields, /// and possibly shallow. Note when shallow=false, new storage is allocated /// for non-OneCopy vars, but the data from src is not actually deep copied - template - void Initialize(const MeshBlockData *src, const std::vector &vars, - const bool shallow_copy) { + template + void Initialize(const std::shared_ptr src, const std::vector &vars = {}, + const bool shallow_copy = false) { + Initialize(src->resolved_packages, src, vars, shallow_copy); + } + + template + void Initialize(const std::shared_ptr resolved_packages_in, + const std::shared_ptr src, const std::vector &vars = {}, + const bool shallow_copy = false) { + if constexpr (!(std::is_same_v> || + std::is_same_v)) { + // We don't allow other types + static_assert(always_false, "Bad source type for initialization."); + } PARTHENON_DEBUG_REQUIRE(src != nullptr, "Source data must be non-null."); SetBlockPointer(src); - resolved_packages_ = src->resolved_packages_; + resolved_packages = resolved_packages_in; is_shallow_ = shallow_copy; // clear all variables, maps, and pack caches @@ -155,23 +157,51 @@ class MeshBlockData { // special case when the list of vars is empty, copy everything if (vars.empty()) { - for (auto v : src->GetVariableVector()) { - add_var(v); + if constexpr (std::is_same_v>) { + for (auto v : src->GetVariableVector()) { + add_var(v); + } + } else if constexpr (std::is_same_v) { + for (auto const &q : resolved_packages->AllFields()) { + AddField(q.first.base_name, q.second, q.first.sparse_id); + } } } else { - for (const auto &v : vars) { - auto var = src->GetVarPtr(v); - add_var(var); - // Add the associated flux as well if not explicitly - // asked for - if (var->IsSet(Metadata::WithFluxes)) { - auto flx_name = var->metadata().GetFluxName(); - bool found = false; - for (const auto &v2 : vars) { - if (src->GetVarPtr(v2)->label() == flx_name) found = true; + if constexpr (std::is_same_v>) { + for (const auto &v : vars) { + auto var = src->GetVarPtr(v); + add_var(var); + // Add the associated flux as well if not explicitly + // asked for + if (var->IsSet(Metadata::WithFluxes)) { + auto flx_name = var->metadata().GetFluxName(); + bool found = false; + for (const auto &v2 : vars) { + if (src->GetVarPtr(v2)->label() == flx_name) found = true; + } + if (!found) add_var(src->GetVarPtr(flx_name)); } - if (!found) add_var(src->GetVarPtr(flx_name)); } + } else { + PARTHENON_FAIL( + "Variable subset selection not yet implemented for MeshBlock input."); + } + } + + // TODO(LFR): Not sure why we only do this in the MeshBlock case, but this carries + // over from the previous iteration. + if constexpr (std::is_same_v) { + if (stage_name_ == "base") { + const auto &swarm_container = GetSwarmData(); + swarm_container->Initialize(resolved_packages, GetBlockSharedPointer()); + } + + // This seems to work fine outside the constexpr if, but having it inside is + // consistent with the old code. + Metadata::FlagCollection flags({Metadata::Sparse, Metadata::ForceAllocOnNewBlocks}); + auto alloc_vars = GetVariablesByFlag(flags); + for (auto &v : alloc_vars.vars()) { + AllocateSparse(v->label()); } } } @@ -573,7 +603,7 @@ class MeshBlockData { } std::weak_ptr pmy_block; - std::shared_ptr resolved_packages_; + std::shared_ptr resolved_packages; bool is_shallow_ = false; const std::string stage_name_; diff --git a/src/mesh/mesh-amr_loadbalance.cpp b/src/mesh/mesh-amr_loadbalance.cpp index be5e7ec08d0e..7c7b418e4891 100644 --- a/src/mesh/mesh-amr_loadbalance.cpp +++ b/src/mesh/mesh-amr_loadbalance.cpp @@ -993,7 +993,7 @@ void Mesh::RedistributeAndRefineMeshBlocks(ParameterInput *pin, ApplicationInput FillDerived(); // Initialize the "base" MeshData object - mesh_data.Get()->Set(block_list, this); + mesh_data.Get()->Initialize(block_list, this); } // AMR Recv and unpack data ResetLoadBalanceVariables(); diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp index 12447244d932..450585e63506 100644 --- a/src/mesh/mesh.cpp +++ b/src/mesh/mesh.cpp @@ -855,7 +855,7 @@ void Mesh::Initialize(bool init_problem, ParameterInput *pin, ApplicationInput * } while (!init_done); // Initialize the "base" MeshData object - mesh_data.Get()->Set(block_list, this); + mesh_data.Get()->Initialize(block_list, this); } /// Finds location of a block with ID `tgid`. diff --git a/src/mesh/meshblock.cpp b/src/mesh/meshblock.cpp index dff5a86af08c..b36fb6489e91 100644 --- a/src/mesh/meshblock.cpp +++ b/src/mesh/meshblock.cpp @@ -142,7 +142,7 @@ void MeshBlock::Initialize(int igid, int ilid, LogicalLocation iloc, // Resolve issues. auto &real_container = meshblock_data.Get(); - real_container->Initialize(resolved_packages, shared_from_this()); + real_container->Initialize(shared_from_this()); // Initialize swarm boundary condition flags real_container->GetSwarmData()->InitializeBoundaries(shared_from_this()); diff --git a/src/outputs/history.cpp b/src/outputs/history.cpp index 2bd1f1b4fc3f..90188278d7fb 100644 --- a/src/outputs/history.cpp +++ b/src/outputs/history.cpp @@ -83,13 +83,13 @@ void HistoryOutput::WriteOutputFile(Mesh *pm, ParameterInput *pin, SimTime *tm, auto &md_base = pm->mesh_data.Get(); // Populated with all blocks if (md_base->NumBlocks() == 0) { - md_base->Set(pm->block_list, pm); + md_base->Initialize(pm->block_list, pm); } else if (md_base->NumBlocks() != pm->block_list.size()) { PARTHENON_WARN( "Resetting \"base\" MeshData to contain all blocks. This indicates that " "the \"base\" MeshData container has been modified elsewhere. Double check " "that the modification was intentional and is compatible with this reset.") - md_base->Set(pm->block_list, pm); + md_base->Initialize(pm->block_list, pm); } // Loop over all packages of the application diff --git a/src/utils/concepts_lite.hpp b/src/utils/concepts_lite.hpp index 43d2344cd855..b3b49d756198 100644 --- a/src/utils/concepts_lite.hpp +++ b/src/utils/concepts_lite.hpp @@ -16,6 +16,13 @@ #include #include +// This is a class template that is required for doing something like static_assert(false) +// in constexpr if blocks. Actually writing static_assert(false) will always cause a +// compilation error, even if it is an unchosen constexpr if block. This is fixed in C++23 +// I think. +template +constexpr std::false_type always_false{}; + // These macros are just to make code more readable and self-explanatory, // generally it is best to write template<..., REQUIRES(... && ...)> in the code // but there are some instance where this causes issues. Switching to the construct diff --git a/tst/unit/test_index_split.cpp b/tst/unit/test_index_split.cpp index a327dc8adbb9..c060f19df325 100644 --- a/tst/unit/test_index_split.cpp +++ b/tst/unit/test_index_split.cpp @@ -99,7 +99,7 @@ TEST_CASE("IndexSplit", "[IndexSplit]") { BlockList_t block_list = MakeBlockList(pkg, NBLOCKS, N, NDIM); MeshData mesh_data("base"); - mesh_data.Set(block_list, nullptr, NDIM); + mesh_data.Initialize(block_list, nullptr, NDIM); WHEN("We initialize an IndexSplit with all outer k and no outer j") { IndexSplit sp(&mesh_data, IndexDomain::interior, IndexSplit::all_outer, diff --git a/tst/unit/test_mesh_data.cpp b/tst/unit/test_mesh_data.cpp index a4cba8f30f8a..f9fa78103ff0 100644 --- a/tst/unit/test_mesh_data.cpp +++ b/tst/unit/test_mesh_data.cpp @@ -78,7 +78,7 @@ TEST_CASE("MeshData works as expected for simple packs", "[MeshData]") { BlockList_t block_list = MakeBlockList(pkg, NBLOCKS, N, NDIM); MeshData mesh_data("base"); - mesh_data.Set(block_list, nullptr); + mesh_data.Initialize(block_list, nullptr); THEN("The number of blocks is correct") { REQUIRE(mesh_data.NumBlocks() == NBLOCKS); } diff --git a/tst/unit/test_sparse_pack.cpp b/tst/unit/test_sparse_pack.cpp index 106dc97befc2..33ffd70e9bb5 100644 --- a/tst/unit/test_sparse_pack.cpp +++ b/tst/unit/test_sparse_pack.cpp @@ -100,7 +100,7 @@ TEST_CASE("Test behavior of sparse packs", "[SparsePack]") { BlockList_t block_list = MakeBlockList(pkg, NBLOCKS, N, NDIM); MeshData mesh_data("base"); - mesh_data.Set(block_list, nullptr); + mesh_data.Initialize(block_list, nullptr); WHEN("We initialize the independent variables by hand and deallocate one") { auto ib = block_list[0]->cellbounds.GetBoundsI(IndexDomain::entire); @@ -168,7 +168,7 @@ TEST_CASE("Test behavior of sparse packs", "[SparsePack]") { BlockList_t block_list = MakeBlockList(pkg, NBLOCKS, N, NDIM); MeshData mesh_data("base"); - mesh_data.Set(block_list, nullptr); + mesh_data.Initialize(block_list, nullptr); WHEN("We initialize the independent variables by hand and deallocate one") { auto ib = block_list[0]->cellbounds.GetBoundsI(IndexDomain::entire);