From 68141738fcd6c8c154be717c18bdd8d6c7647a46 Mon Sep 17 00:00:00 2001 From: adam reyes Date: Thu, 26 Sep 2024 16:19:29 +0200 Subject: [PATCH 01/49] meshdata version of refinement tagging --- src/amr_criteria/amr_criteria.cpp | 23 ++++ src/amr_criteria/amr_criteria.hpp | 8 ++ src/amr_criteria/refinement_package.cpp | 167 +++++++++++++++++++++++- src/amr_criteria/refinement_package.hpp | 15 ++- src/interface/state_descriptor.hpp | 6 + 5 files changed, 213 insertions(+), 6 deletions(-) diff --git a/src/amr_criteria/amr_criteria.cpp b/src/amr_criteria/amr_criteria.cpp index 41b2237d96ff..e215f3f7d4c6 100644 --- a/src/amr_criteria/amr_criteria.cpp +++ b/src/amr_criteria/amr_criteria.cpp @@ -95,6 +95,17 @@ AmrTag AMRFirstDerivative::operator()(const MeshBlockData *rc) const { return Refinement::FirstDerivative(bnds, q, refine_criteria, derefine_criteria); } +void AMRFirstDerivative::operator()(MeshData *mc, + const std::vector &fields, + ParArray1D &delta_levels) const { + auto ib = mc->GetBoundsI(IndexDomain::interior); + auto jb = mc->GetBoundsJ(IndexDomain::interior); + auto kb = mc->GetBoundsK(IndexDomain::interior); + auto bnds = AMRBounds(ib, jb, kb); + Refinement::FirstDerivative(bnds, mc, fields, delta_levels, refine_criteria, + derefine_criteria); +} + AmrTag AMRSecondDerivative::operator()(const MeshBlockData *rc) const { if (!rc->HasVariable(field) || !rc->IsAllocated(field)) { return AmrTag::same; @@ -102,7 +113,19 @@ AmrTag AMRSecondDerivative::operator()(const MeshBlockData *rc) const { auto bnds = GetBounds(rc); auto q = Kokkos::subview(rc->Get(field).data, comp6, comp5, comp4, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()); + printf("first deriv woo\n"); return Refinement::SecondDerivative(bnds, q, refine_criteria, derefine_criteria); } +void AMRSecondDerivative::operator()(MeshData *mc, + const std::vector &fields, + ParArray1D &delta_levels) const { + auto ib = mc->GetBoundsI(IndexDomain::interior); + auto jb = mc->GetBoundsJ(IndexDomain::interior); + auto kb = mc->GetBoundsK(IndexDomain::interior); + auto bnds = AMRBounds(ib, jb, kb); + Refinement::SecondDerivative(bnds, mc, fields, delta_levels, refine_criteria, + derefine_criteria); +} + } // namespace parthenon diff --git a/src/amr_criteria/amr_criteria.hpp b/src/amr_criteria/amr_criteria.hpp index e31bd8b381d8..73509a1d1b0f 100644 --- a/src/amr_criteria/amr_criteria.hpp +++ b/src/amr_criteria/amr_criteria.hpp @@ -15,9 +15,11 @@ #include #include +#include #include "defs.hpp" #include "mesh/domain.hpp" +#include "mesh/mesh.hpp" namespace parthenon { @@ -36,6 +38,8 @@ struct AMRCriteria { AMRCriteria(ParameterInput *pin, std::string &block_name); virtual ~AMRCriteria() {} virtual AmrTag operator()(const MeshBlockData *rc) const = 0; + virtual void operator()(MeshData *mc, const std::vector &fields, + ParArray1D &delta_level) const = 0; std::string field; Real refine_criteria, derefine_criteria; int max_level; @@ -49,12 +53,16 @@ struct AMRFirstDerivative : public AMRCriteria { AMRFirstDerivative(ParameterInput *pin, std::string &block_name) : AMRCriteria(pin, block_name) {} AmrTag operator()(const MeshBlockData *rc) const override; + void operator()(MeshData *mc, const std::vector &fields, + ParArray1D &delta_level) const override; }; struct AMRSecondDerivative : public AMRCriteria { AMRSecondDerivative(ParameterInput *pin, std::string &block_name) : AMRCriteria(pin, block_name) {} AmrTag operator()(const MeshBlockData *rc) const override; + void operator()(MeshData *mc, const std::vector &fields, + ParArray1D &delta_level) const override; }; } // namespace parthenon diff --git a/src/amr_criteria/refinement_package.cpp b/src/amr_criteria/refinement_package.cpp index 81342b70695a..d9df19bd5b6d 100644 --- a/src/amr_criteria/refinement_package.cpp +++ b/src/amr_criteria/refinement_package.cpp @@ -18,11 +18,15 @@ #include #include #include +#include +#include "Kokkos_Macros.hpp" #include "amr_criteria/amr_criteria.hpp" +#include "interface/make_pack_descriptor.hpp" #include "interface/mesh_data.hpp" #include "interface/meshblock_data.hpp" #include "interface/state_descriptor.hpp" +#include "kokkos_abstraction.hpp" #include "mesh/mesh.hpp" #include "mesh/mesh_refinement.hpp" #include "mesh/meshblock.hpp" @@ -33,6 +37,9 @@ namespace Refinement { std::shared_ptr Initialize(ParameterInput *pin) { auto ref = std::make_shared("Refinement"); + bool check_refine_mesh = + pin->GetOrAddBoolean("parthenon/mesh", "CheckRefineMesh", false); + ref->AddParam("check_refine_mesh", check_refine_mesh); int numcrit = 0; while (true) { @@ -45,10 +52,43 @@ std::shared_ptr Initialize(ParameterInput *pin) { ref->amr_criteria.push_back(AMRCriteria::MakeAMRCriteria(method, pin, block_name)); numcrit++; } + std::vector fields; + for (auto &amr : ref->amr_criteria) { + fields.push_back(amr->field); + } + ref->AddParam("refine_fields", fields); return ref; } -AmrTag CheckAllRefinement(MeshBlockData *rc) { +ParArray1D CheckAllRefinement(MeshData *mc) { + printf("allrefine!!\n"); + const int nblocks = mc->NumBlocks(); + // maybe not great to allocate this all the time + auto delta_levels = ParArray1D(Kokkos::View( + Kokkos::view_alloc(Kokkos::WithoutInitializing, "delta_levels"), nblocks)); + Kokkos::deep_copy(delta_levels.KokkosView(), AmrTag::derefine); + + Mesh *pm = mc->GetMeshPointer(); + static const bool check_refine_mesh = + pm->packages.Get("Refinement")->Param("check_refine_mesh"); + static const std::vector refine_fields = + pm->packages.Get("Refinement")->Param>("refine_fields"); + + for (auto &pkg : pm->packages.AllPackages()) { + auto &desc = pkg.second; + desc->CheckRefinement(mc, delta_levels); + + if (check_refine_mesh) { + for (auto &amr : desc->amr_criteria) { + (*amr)(mc, refine_fields, delta_levels); + } + } + } + + return delta_levels; +} + +AmrTag CheckAllRefinement(MeshBlockData *rc, const AmrTag &level) { // Check all refinement criteria and return the maximum recommended change in // refinement level: // delta_level = -1 => recommend derefinement @@ -62,8 +102,10 @@ AmrTag CheckAllRefinement(MeshBlockData *rc) { // neighboring blocks. Similarly for "do nothing" PARTHENON_INSTRUMENT MeshBlock *pmb = rc->GetBlockPointer(); + static const bool check_refine_mesh = + pmb->packages.Get("Refinement")->Param("check_refine_mesh"); // delta_level holds the max over all criteria. default to derefining. - AmrTag delta_level = AmrTag::derefine; + AmrTag delta_level = level; for (auto &pkg : pmb->packages.AllPackages()) { auto &desc = pkg.second; delta_level = std::max(delta_level, desc->CheckRefinement(rc)); @@ -71,6 +113,7 @@ AmrTag CheckAllRefinement(MeshBlockData *rc) { // since 1 is the max, we can return without having to look at anything else return AmrTag::refine; } + if (check_refine_mesh) continue; // call parthenon criteria that were registered for (auto &amr : desc->amr_criteria) { // get the recommended change in refinement level from this criteria @@ -150,9 +193,120 @@ AmrTag SecondDerivative(const AMRBounds &bnds, const ParArray3D &q, return AmrTag::same; } -void SetRefinement_(MeshBlockData *rc) { +void FirstDerivative(const AMRBounds &bnds, MeshData *mc, + const std::vector &fields, + ParArray1D &delta_levels_, const Real refine_criteria_, + const Real derefine_criteria_) { + std::vector> var_regex; + for (std::string var : fields) { + var_regex.push_back(std::pair(var, false)); + } + + static auto desc = + MakePackDescriptor(mc->GetMeshPointer()->resolved_packages.get(), var_regex); + auto pack = desc.GetPack(mc); + const int ndim = mc->GetMeshPointer()->ndim; + const int nvars = pack.GetMaxNumberOfVars(); + + const Real refine_criteria = refine_criteria_; + const Real derefine_criteria = derefine_criteria_; + auto delta_levels = delta_levels_; + par_for_outer( + PARTHENON_AUTO_LABEL, 0, 0, 0, pack.GetNBlocks() - 1, + KOKKOS_LAMBDA(team_mbr_t team_member, const int b) { + if (delta_levels(b) == AmrTag::refine) return; + Real maxd = 0.; + for (int var = 0; var < nvars; var++) { + par_reduce_inner( + inner_loop_pattern_ttr_tag, team_member, bnds.ks, bnds.ke, bnds.js, bnds.je, + bnds.is, bnds.ie, + [&](const int k, const int j, const int i, Real &maxder) { + Real scale = std::abs(pack(b, var, k, j, i)); + Real d = + 0.5 * + std::abs((pack(b, var, k, j, i + 1) - pack(b, var, k, j, i - 1))) / + (scale + TINY_NUMBER); + maxd = (d > maxd ? d : maxd); + if (ndim > 1) { + d = 0.5 * + std::abs((pack(b, var, k, j + 1, i) - pack(b, var, k, j - 1, i))) / + (scale + TINY_NUMBER); + maxd = (d > maxd ? d : maxd); + } + if (ndim > 2) { + d = 0.5 * + std::abs((pack(b, var, k + 1, j, i) - pack(b, var, k - 1, j, i))) / + (scale + TINY_NUMBER); + maxd = (d > maxd ? d : maxd); + } + }, + Kokkos::Max(maxd)); + + AmrTag flag = AmrTag::same; + if (maxd > refine_criteria) flag = AmrTag::refine; + if (maxd < derefine_criteria) flag = AmrTag::derefine; + delta_levels(b) = std::max(delta_levels(b), flag); + } + }); +} + +void SecondDerivative(const AMRBounds &bnds, MeshData *mc, + const std::vector &fields, + ParArray1D &delta_levels_, const Real refine_criteria_, + const Real derefine_criteria_) { + std::vector> var_regex; + for (std::string var : fields) { + var_regex.push_back(std::pair(var, false)); + } + + static auto desc = + MakePackDescriptor(mc->GetMeshPointer()->resolved_packages.get(), var_regex); + auto pack = desc.GetPack(mc); + const int ndim = mc->GetMeshPointer()->ndim; + const int nvars = pack.GetMaxNumberOfVars(); + + const Real refine_criteria = refine_criteria_; + const Real derefine_criteria = derefine_criteria_; + auto delta_levels = delta_levels_; + par_for_outer( + PARTHENON_AUTO_LABEL, 0, 0, 0, pack.GetNBlocks() - 1, + KOKKOS_LAMBDA(team_mbr_t team_member, const int b) { + if (delta_levels(b) == AmrTag::refine) return; + Real maxd = 0.; + for (int var = 0; var < nvars; var++) { + par_reduce_inner( + inner_loop_pattern_ttr_tag, team_member, bnds.ks, bnds.ke, bnds.js, bnds.je, + bnds.is, bnds.ie, + [&](const int k, const int j, const int i, Real &maxder) { + Real aqt = std::abs(pack(b, var, k, j, i)) + TINY_NUMBER; + Real qavg = 0.5 * (pack(b, var, k, j, i + 1) + pack(b, var, k, j, i - 1)); + Real d = std::abs(qavg - pack(b, var, k, j, i)) / (std::abs(qavg) + aqt); + maxd = (d > maxd ? d : maxd); + if (ndim > 1) { + qavg = 0.5 * (pack(b, var, k, j + 1, i) + pack(b, var, k, j - 1, i)); + d = std::abs(qavg - pack(b, var, k, j, i)) / (std::abs(qavg) + aqt); + maxd = (d > maxd ? d : maxd); + } + if (ndim > 2) { + qavg = 0.5 * (pack(b, var, k + 1, j, i) + pack(b, var, k - 1, j, i)); + d = std::abs(qavg - pack(b, var, k, j, i)) / (std::abs(qavg) + aqt); + maxd = (d > maxd ? d : maxd); + } + }, + Kokkos::Max(maxd)); + + AmrTag flag = AmrTag::same; + if (maxd > refine_criteria) flag = AmrTag::refine; + if (maxd < derefine_criteria) flag = AmrTag::derefine; + delta_levels(b) = std::max(delta_levels(b), flag); + } + }); +} + +void SetRefinement_(MeshBlockData *rc, + const AmrTag &delta_level = AmrTag::derefine) { auto pmb = rc->GetBlockPointer(); - pmb->pmr->SetRefinement(CheckAllRefinement(rc)); + pmb->pmr->SetRefinement(CheckAllRefinement(rc, delta_level)); } template <> @@ -165,8 +319,11 @@ TaskStatus Tag(MeshBlockData *rc) { template <> TaskStatus Tag(MeshData *rc) { PARTHENON_INSTRUMENT + ParArray1D delta_levels = CheckAllRefinement(rc); + auto delta_levels_h = delta_levels.GetHostMirrorAndCopy(); + for (int i = 0; i < rc->NumBlocks(); i++) { - SetRefinement_(rc->GetBlockData(i).get()); + SetRefinement_(rc->GetBlockData(i).get(), delta_levels_h(i)); } return TaskStatus::complete; } diff --git a/src/amr_criteria/refinement_package.hpp b/src/amr_criteria/refinement_package.hpp index 40e9cd39fb14..aaaba6c6b787 100644 --- a/src/amr_criteria/refinement_package.hpp +++ b/src/amr_criteria/refinement_package.hpp @@ -16,6 +16,7 @@ #include #include +#include #include "defs.hpp" #include "parthenon_arrays.hpp" @@ -37,14 +38,26 @@ std::shared_ptr Initialize(ParameterInput *pin); template TaskStatus Tag(T *rc); -AmrTag CheckAllRefinement(MeshBlockData *rc); +AmrTag CheckAllRefinement(MeshBlockData *rc, + const AmrTag &level = AmrTag::derefine); +void CheckAllRefinement(MeshData *mc, ParArray1D &delta_level); AmrTag FirstDerivative(const AMRBounds &bnds, const ParArray3D &q, const Real refine_criteria, const Real derefine_criteria); +void FirstDerivative(const AMRBounds &bnds, MeshData *mc, + const std::vector &fields, + ParArray1D &delta_levels_, const Real refine_criteria_, + const Real derefine_criteria_); + AmrTag SecondDerivative(const AMRBounds &bnds, const ParArray3D &q, const Real refine_criteria, const Real derefine_criteria); +void SecondDerivative(const AMRBounds &bnds, MeshData *mc, + const std::vector &fields, + ParArray1D &delta_levels_, const Real refine_criteria_, + const Real derefine_criteria_); + } // namespace Refinement } // namespace parthenon diff --git a/src/interface/state_descriptor.hpp b/src/interface/state_descriptor.hpp index 106c3ea6e525..c3a034bb377a 100644 --- a/src/interface/state_descriptor.hpp +++ b/src/interface/state_descriptor.hpp @@ -326,6 +326,10 @@ class StateDescriptor { return AmrTag::derefine; } + void CheckRefinement(MeshData *mc, ParArray1D &delta_level) const { + if (CheckRefinementMesh != nullptr) CheckRefinementMesh(mc, delta_level); + } + void InitNewlyAllocatedVars(MeshData *rc) const { if (InitNewlyAllocatedVarsMesh != nullptr) return InitNewlyAllocatedVarsMesh(rc); } @@ -375,6 +379,8 @@ class StateDescriptor { std::function *rc)> EstimateTimestepMesh = nullptr; std::function *rc)> CheckRefinementBlock = nullptr; + std::function *rc, ParArray1D &delta_level)> + CheckRefinementMesh = nullptr; std::function *rc)> InitNewlyAllocatedVarsMesh = nullptr; std::function *rc)> InitNewlyAllocatedVarsBlock = nullptr; From 6df75e7f8263776bb589c67ee597f234e71dfca8 Mon Sep 17 00:00:00 2001 From: adam reyes Date: Fri, 27 Sep 2024 00:29:09 +0200 Subject: [PATCH 02/49] got initialization and vector indices working --- src/amr_criteria/amr_criteria.cpp | 11 +- src/amr_criteria/amr_criteria.hpp | 9 +- src/amr_criteria/refinement_package.cpp | 151 ++++++++++-------------- src/amr_criteria/refinement_package.hpp | 16 ++- src/mesh/mesh.cpp | 6 +- 5 files changed, 83 insertions(+), 110 deletions(-) diff --git a/src/amr_criteria/amr_criteria.cpp b/src/amr_criteria/amr_criteria.cpp index e215f3f7d4c6..a7597d7ffa07 100644 --- a/src/amr_criteria/amr_criteria.cpp +++ b/src/amr_criteria/amr_criteria.cpp @@ -96,14 +96,13 @@ AmrTag AMRFirstDerivative::operator()(const MeshBlockData *rc) const { } void AMRFirstDerivative::operator()(MeshData *mc, - const std::vector &fields, ParArray1D &delta_levels) const { auto ib = mc->GetBoundsI(IndexDomain::interior); auto jb = mc->GetBoundsJ(IndexDomain::interior); auto kb = mc->GetBoundsK(IndexDomain::interior); auto bnds = AMRBounds(ib, jb, kb); - Refinement::FirstDerivative(bnds, mc, fields, delta_levels, refine_criteria, - derefine_criteria); + Refinement::FirstDerivative(bnds, mc, field, {comp6, comp5, comp4}, delta_levels, + refine_criteria, derefine_criteria); } AmrTag AMRSecondDerivative::operator()(const MeshBlockData *rc) const { @@ -113,19 +112,17 @@ AmrTag AMRSecondDerivative::operator()(const MeshBlockData *rc) const { auto bnds = GetBounds(rc); auto q = Kokkos::subview(rc->Get(field).data, comp6, comp5, comp4, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()); - printf("first deriv woo\n"); return Refinement::SecondDerivative(bnds, q, refine_criteria, derefine_criteria); } void AMRSecondDerivative::operator()(MeshData *mc, - const std::vector &fields, ParArray1D &delta_levels) const { auto ib = mc->GetBoundsI(IndexDomain::interior); auto jb = mc->GetBoundsJ(IndexDomain::interior); auto kb = mc->GetBoundsK(IndexDomain::interior); auto bnds = AMRBounds(ib, jb, kb); - Refinement::SecondDerivative(bnds, mc, fields, delta_levels, refine_criteria, - derefine_criteria); + Refinement::SecondDerivative(bnds, mc, field, {comp6, comp5, comp4}, delta_levels, + refine_criteria, derefine_criteria); } } // namespace parthenon diff --git a/src/amr_criteria/amr_criteria.hpp b/src/amr_criteria/amr_criteria.hpp index 73509a1d1b0f..411646e9ac9e 100644 --- a/src/amr_criteria/amr_criteria.hpp +++ b/src/amr_criteria/amr_criteria.hpp @@ -38,8 +38,7 @@ struct AMRCriteria { AMRCriteria(ParameterInput *pin, std::string &block_name); virtual ~AMRCriteria() {} virtual AmrTag operator()(const MeshBlockData *rc) const = 0; - virtual void operator()(MeshData *mc, const std::vector &fields, - ParArray1D &delta_level) const = 0; + virtual void operator()(MeshData *mc, ParArray1D &delta_level) const = 0; std::string field; Real refine_criteria, derefine_criteria; int max_level; @@ -53,16 +52,14 @@ struct AMRFirstDerivative : public AMRCriteria { AMRFirstDerivative(ParameterInput *pin, std::string &block_name) : AMRCriteria(pin, block_name) {} AmrTag operator()(const MeshBlockData *rc) const override; - void operator()(MeshData *mc, const std::vector &fields, - ParArray1D &delta_level) const override; + void operator()(MeshData *mc, ParArray1D &delta_level) const override; }; struct AMRSecondDerivative : public AMRCriteria { AMRSecondDerivative(ParameterInput *pin, std::string &block_name) : AMRCriteria(pin, block_name) {} AmrTag operator()(const MeshBlockData *rc) const override; - void operator()(MeshData *mc, const std::vector &fields, - ParArray1D &delta_level) const override; + void operator()(MeshData *mc, ParArray1D &delta_level) const override; }; } // namespace parthenon diff --git a/src/amr_criteria/refinement_package.cpp b/src/amr_criteria/refinement_package.cpp index d9df19bd5b6d..dc4efd854f63 100644 --- a/src/amr_criteria/refinement_package.cpp +++ b/src/amr_criteria/refinement_package.cpp @@ -52,16 +52,10 @@ std::shared_ptr Initialize(ParameterInput *pin) { ref->amr_criteria.push_back(AMRCriteria::MakeAMRCriteria(method, pin, block_name)); numcrit++; } - std::vector fields; - for (auto &amr : ref->amr_criteria) { - fields.push_back(amr->field); - } - ref->AddParam("refine_fields", fields); return ref; } ParArray1D CheckAllRefinement(MeshData *mc) { - printf("allrefine!!\n"); const int nblocks = mc->NumBlocks(); // maybe not great to allocate this all the time auto delta_levels = ParArray1D(Kokkos::View( @@ -71,8 +65,6 @@ ParArray1D CheckAllRefinement(MeshData *mc) { Mesh *pm = mc->GetMeshPointer(); static const bool check_refine_mesh = pm->packages.Get("Refinement")->Param("check_refine_mesh"); - static const std::vector refine_fields = - pm->packages.Get("Refinement")->Param>("refine_fields"); for (auto &pkg : pm->packages.AllPackages()) { auto &desc = pkg.second; @@ -80,7 +72,7 @@ ParArray1D CheckAllRefinement(MeshData *mc) { if (check_refine_mesh) { for (auto &amr : desc->amr_criteria) { - (*amr)(mc, refine_fields, delta_levels); + (*amr)(mc, delta_levels); } } } @@ -193,17 +185,11 @@ AmrTag SecondDerivative(const AMRBounds &bnds, const ParArray3D &q, return AmrTag::same; } -void FirstDerivative(const AMRBounds &bnds, MeshData *mc, - const std::vector &fields, - ParArray1D &delta_levels_, const Real refine_criteria_, - const Real derefine_criteria_) { - std::vector> var_regex; - for (std::string var : fields) { - var_regex.push_back(std::pair(var, false)); - } - - static auto desc = - MakePackDescriptor(mc->GetMeshPointer()->resolved_packages.get(), var_regex); +void FirstDerivative(const AMRBounds &bnds, MeshData *mc, const std::string &field, + Kokkos::Array index, ParArray1D &delta_levels_, + const Real refine_criteria_, const Real derefine_criteria_) { + const auto desc = + MakePackDescriptor(mc->GetMeshPointer()->resolved_packages.get(), {field}); auto pack = desc.GetPack(mc); const int ndim = mc->GetMeshPointer()->ndim; const int nvars = pack.GetMaxNumberOfVars(); @@ -211,56 +197,49 @@ void FirstDerivative(const AMRBounds &bnds, MeshData *mc, const Real refine_criteria = refine_criteria_; const Real derefine_criteria = derefine_criteria_; auto delta_levels = delta_levels_; + // TODO(?): I have no idea how to get the index from the tensor indices... + const int var = index[2]; par_for_outer( PARTHENON_AUTO_LABEL, 0, 0, 0, pack.GetNBlocks() - 1, KOKKOS_LAMBDA(team_mbr_t team_member, const int b) { if (delta_levels(b) == AmrTag::refine) return; Real maxd = 0.; - for (int var = 0; var < nvars; var++) { - par_reduce_inner( - inner_loop_pattern_ttr_tag, team_member, bnds.ks, bnds.ke, bnds.js, bnds.je, - bnds.is, bnds.ie, - [&](const int k, const int j, const int i, Real &maxder) { - Real scale = std::abs(pack(b, var, k, j, i)); - Real d = - 0.5 * - std::abs((pack(b, var, k, j, i + 1) - pack(b, var, k, j, i - 1))) / + par_reduce_inner( + inner_loop_pattern_ttr_tag, team_member, bnds.ks, bnds.ke, bnds.js, bnds.je, + bnds.is, bnds.ie, + [&](const int k, const int j, const int i, Real &maxder) { + Real scale = std::abs(pack(b, var, k, j, i)); + Real d = 0.5 * + std::abs((pack(b, var, k, j, i + 1) - pack(b, var, k, j, i - 1))) / + (scale + TINY_NUMBER); + maxder = (d > maxder ? d : maxder); + if (ndim > 1) { + d = 0.5 * + std::abs((pack(b, var, k, j + 1, i) - pack(b, var, k, j - 1, i))) / + (scale + TINY_NUMBER); + maxder = (d > maxder ? d : maxder); + } + if (ndim > 2) { + d = 0.5 * + std::abs((pack(b, var, k + 1, j, i) - pack(b, var, k - 1, j, i))) / (scale + TINY_NUMBER); - maxd = (d > maxd ? d : maxd); - if (ndim > 1) { - d = 0.5 * - std::abs((pack(b, var, k, j + 1, i) - pack(b, var, k, j - 1, i))) / - (scale + TINY_NUMBER); - maxd = (d > maxd ? d : maxd); - } - if (ndim > 2) { - d = 0.5 * - std::abs((pack(b, var, k + 1, j, i) - pack(b, var, k - 1, j, i))) / - (scale + TINY_NUMBER); - maxd = (d > maxd ? d : maxd); - } - }, - Kokkos::Max(maxd)); + maxder = (d > maxder ? d : maxder); + } + }, + Kokkos::Max(maxd)); - AmrTag flag = AmrTag::same; - if (maxd > refine_criteria) flag = AmrTag::refine; - if (maxd < derefine_criteria) flag = AmrTag::derefine; - delta_levels(b) = std::max(delta_levels(b), flag); - } + AmrTag flag = AmrTag::same; + if (maxd > refine_criteria) flag = AmrTag::refine; + if (maxd < derefine_criteria) flag = AmrTag::derefine; + delta_levels(b) = std::max(delta_levels(b), flag); }); } -void SecondDerivative(const AMRBounds &bnds, MeshData *mc, - const std::vector &fields, - ParArray1D &delta_levels_, const Real refine_criteria_, - const Real derefine_criteria_) { - std::vector> var_regex; - for (std::string var : fields) { - var_regex.push_back(std::pair(var, false)); - } - - static auto desc = - MakePackDescriptor(mc->GetMeshPointer()->resolved_packages.get(), var_regex); +void SecondDerivative(const AMRBounds &bnds, MeshData *mc, const std::string &field, + Kokkos::Array index, ParArray1D &delta_levels_, + const Real refine_criteria_, const Real derefine_criteria_) { + const auto desc = + MakePackDescriptor(mc->GetMeshPointer()->resolved_packages.get(), {field}); auto pack = desc.GetPack(mc); const int ndim = mc->GetMeshPointer()->ndim; const int nvars = pack.GetMaxNumberOfVars(); @@ -268,38 +247,38 @@ void SecondDerivative(const AMRBounds &bnds, MeshData *mc, const Real refine_criteria = refine_criteria_; const Real derefine_criteria = derefine_criteria_; auto delta_levels = delta_levels_; + // TODO(?): I have no idea how to get the index from the tensor indices... + const int var = index[2]; par_for_outer( PARTHENON_AUTO_LABEL, 0, 0, 0, pack.GetNBlocks() - 1, KOKKOS_LAMBDA(team_mbr_t team_member, const int b) { if (delta_levels(b) == AmrTag::refine) return; Real maxd = 0.; - for (int var = 0; var < nvars; var++) { - par_reduce_inner( - inner_loop_pattern_ttr_tag, team_member, bnds.ks, bnds.ke, bnds.js, bnds.je, - bnds.is, bnds.ie, - [&](const int k, const int j, const int i, Real &maxder) { - Real aqt = std::abs(pack(b, var, k, j, i)) + TINY_NUMBER; - Real qavg = 0.5 * (pack(b, var, k, j, i + 1) + pack(b, var, k, j, i - 1)); - Real d = std::abs(qavg - pack(b, var, k, j, i)) / (std::abs(qavg) + aqt); - maxd = (d > maxd ? d : maxd); - if (ndim > 1) { - qavg = 0.5 * (pack(b, var, k, j + 1, i) + pack(b, var, k, j - 1, i)); - d = std::abs(qavg - pack(b, var, k, j, i)) / (std::abs(qavg) + aqt); - maxd = (d > maxd ? d : maxd); - } - if (ndim > 2) { - qavg = 0.5 * (pack(b, var, k + 1, j, i) + pack(b, var, k - 1, j, i)); - d = std::abs(qavg - pack(b, var, k, j, i)) / (std::abs(qavg) + aqt); - maxd = (d > maxd ? d : maxd); - } - }, - Kokkos::Max(maxd)); + par_reduce_inner( + inner_loop_pattern_ttr_tag, team_member, bnds.ks, bnds.ke, bnds.js, bnds.je, + bnds.is, bnds.ie, + [&](const int k, const int j, const int i, Real &maxder) { + Real aqt = std::abs(pack(b, var, k, j, i)) + TINY_NUMBER; + Real qavg = 0.5 * (pack(b, var, k, j, i + 1) + pack(b, var, k, j, i - 1)); + Real d = std::abs(qavg - pack(b, var, k, j, i)) / (std::abs(qavg) + aqt); + maxder = (d > maxder ? d : maxder); + if (ndim > 1) { + qavg = 0.5 * (pack(b, var, k, j + 1, i) + pack(b, var, k, j - 1, i)); + d = std::abs(qavg - pack(b, var, k, j, i)) / (std::abs(qavg) + aqt); + maxder = (d > maxder ? d : maxder); + } + if (ndim > 2) { + qavg = 0.5 * (pack(b, var, k + 1, j, i) + pack(b, var, k - 1, j, i)); + d = std::abs(qavg - pack(b, var, k, j, i)) / (std::abs(qavg) + aqt); + maxder = (d > maxder ? d : maxder); + } + }, + Kokkos::Max(maxd)); - AmrTag flag = AmrTag::same; - if (maxd > refine_criteria) flag = AmrTag::refine; - if (maxd < derefine_criteria) flag = AmrTag::derefine; - delta_levels(b) = std::max(delta_levels(b), flag); - } + AmrTag flag = AmrTag::same; + if (maxd > refine_criteria) flag = AmrTag::refine; + if (maxd < derefine_criteria) flag = AmrTag::derefine; + delta_levels(b) = std::max(delta_levels(b), flag); }); } diff --git a/src/amr_criteria/refinement_package.hpp b/src/amr_criteria/refinement_package.hpp index aaaba6c6b787..e0c1b2329283 100644 --- a/src/amr_criteria/refinement_package.hpp +++ b/src/amr_criteria/refinement_package.hpp @@ -40,23 +40,21 @@ TaskStatus Tag(T *rc); AmrTag CheckAllRefinement(MeshBlockData *rc, const AmrTag &level = AmrTag::derefine); -void CheckAllRefinement(MeshData *mc, ParArray1D &delta_level); +ParArray1D CheckAllRefinement(MeshData *mc); AmrTag FirstDerivative(const AMRBounds &bnds, const ParArray3D &q, const Real refine_criteria, const Real derefine_criteria); -void FirstDerivative(const AMRBounds &bnds, MeshData *mc, - const std::vector &fields, - ParArray1D &delta_levels_, const Real refine_criteria_, - const Real derefine_criteria_); +void FirstDerivative(const AMRBounds &bnds, MeshData *mc, const std::string &field, + Kokkos::Array index, ParArray1D &delta_levels_, + const Real refine_criteria_, const Real derefine_criteria_); AmrTag SecondDerivative(const AMRBounds &bnds, const ParArray3D &q, const Real refine_criteria, const Real derefine_criteria); -void SecondDerivative(const AMRBounds &bnds, MeshData *mc, - const std::vector &fields, - ParArray1D &delta_levels_, const Real refine_criteria_, - const Real derefine_criteria_); +void SecondDerivative(const AMRBounds &bnds, MeshData *mc, const std::string &field, + Kokkos::Array index, ParArray1D &delta_levels_, + const Real refine_criteria_, const Real derefine_criteria_); } // namespace Refinement diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp index d513ec6e0b52..ea838580c2b1 100644 --- a/src/mesh/mesh.cpp +++ b/src/mesh/mesh.cpp @@ -39,6 +39,7 @@ #include "bvals/comms/bvals_in_one.hpp" #include "parthenon_mpi.hpp" +#include "amr_criteria/refinement_package.hpp" #include "application_input.hpp" #include "bvals/boundary_conditions.hpp" #include "bvals/bvals.hpp" @@ -820,8 +821,9 @@ void Mesh::Initialize(bool init_problem, ParameterInput *pin, ApplicationInput * FillDerived(); if (init_problem && adaptive) { - for (int i = 0; i < nmb; ++i) { - block_list[i]->pmr->CheckRefinementCondition(); + for (auto &partition : GetDefaultBlockPartitions(GridIdentifier::leaf())) { + auto &md = mesh_data.Add("base", partition); + Refinement::Tag(md.get()); } init_done = false; // caching nbtotal the private variable my be updated in the following function From 2fdc2ed93101d49b6df574fcf7f3074171bb4892 Mon Sep 17 00:00:00 2001 From: adam reyes Date: Fri, 27 Sep 2024 13:53:47 +0200 Subject: [PATCH 03/49] fix tensor indices --- src/amr_criteria/amr_criteria.cpp | 28 +++++++++++++++++++++---- src/amr_criteria/refinement_package.cpp | 10 ++++----- src/amr_criteria/refinement_package.hpp | 4 ++-- 3 files changed, 30 insertions(+), 12 deletions(-) diff --git a/src/amr_criteria/amr_criteria.cpp b/src/amr_criteria/amr_criteria.cpp index a7597d7ffa07..d867ff39604e 100644 --- a/src/amr_criteria/amr_criteria.cpp +++ b/src/amr_criteria/amr_criteria.cpp @@ -101,8 +101,18 @@ void AMRFirstDerivative::operator()(MeshData *mc, auto jb = mc->GetBoundsJ(IndexDomain::interior); auto kb = mc->GetBoundsK(IndexDomain::interior); auto bnds = AMRBounds(ib, jb, kb); - Refinement::FirstDerivative(bnds, mc, field, {comp6, comp5, comp4}, delta_levels, - refine_criteria, derefine_criteria); + auto dims = mc->GetMeshPointer()->resolved_packages->FieldMetadata(field).Shape(); + int n5(0), n4(0); + if (dims.size() > 2) { + n5 = dims[1]; + n4 = dims[2]; + } else if (dims.size() > 1) { + n5 = dims[0]; + n4 = dims[1]; + } + const int idx = comp4 + n4 * (comp5 + n5 * comp6); + Refinement::FirstDerivative(bnds, mc, field, idx, delta_levels, refine_criteria, + derefine_criteria); } AmrTag AMRSecondDerivative::operator()(const MeshBlockData *rc) const { @@ -121,8 +131,18 @@ void AMRSecondDerivative::operator()(MeshData *mc, auto jb = mc->GetBoundsJ(IndexDomain::interior); auto kb = mc->GetBoundsK(IndexDomain::interior); auto bnds = AMRBounds(ib, jb, kb); - Refinement::SecondDerivative(bnds, mc, field, {comp6, comp5, comp4}, delta_levels, - refine_criteria, derefine_criteria); + auto dims = mc->GetMeshPointer()->resolved_packages->FieldMetadata(field).Shape(); + int n5(0), n4(0); + if (dims.size() > 2) { + n5 = dims[1]; + n4 = dims[2]; + } else if (dims.size() > 1) { + n5 = dims[0]; + n4 = dims[1]; + } + const int idx = comp4 + n4 * (comp5 + n5 * comp6); + Refinement::SecondDerivative(bnds, mc, field, idx, delta_levels, refine_criteria, + derefine_criteria); } } // namespace parthenon diff --git a/src/amr_criteria/refinement_package.cpp b/src/amr_criteria/refinement_package.cpp index dc4efd854f63..7f57f857a3dd 100644 --- a/src/amr_criteria/refinement_package.cpp +++ b/src/amr_criteria/refinement_package.cpp @@ -186,7 +186,7 @@ AmrTag SecondDerivative(const AMRBounds &bnds, const ParArray3D &q, } void FirstDerivative(const AMRBounds &bnds, MeshData *mc, const std::string &field, - Kokkos::Array index, ParArray1D &delta_levels_, + const int &idx, ParArray1D &delta_levels_, const Real refine_criteria_, const Real derefine_criteria_) { const auto desc = MakePackDescriptor(mc->GetMeshPointer()->resolved_packages.get(), {field}); @@ -197,8 +197,7 @@ void FirstDerivative(const AMRBounds &bnds, MeshData *mc, const std::strin const Real refine_criteria = refine_criteria_; const Real derefine_criteria = derefine_criteria_; auto delta_levels = delta_levels_; - // TODO(?): I have no idea how to get the index from the tensor indices... - const int var = index[2]; + const int var = idx; par_for_outer( PARTHENON_AUTO_LABEL, 0, 0, 0, pack.GetNBlocks() - 1, KOKKOS_LAMBDA(team_mbr_t team_member, const int b) { @@ -236,7 +235,7 @@ void FirstDerivative(const AMRBounds &bnds, MeshData *mc, const std::strin } void SecondDerivative(const AMRBounds &bnds, MeshData *mc, const std::string &field, - Kokkos::Array index, ParArray1D &delta_levels_, + const int &idx, ParArray1D &delta_levels_, const Real refine_criteria_, const Real derefine_criteria_) { const auto desc = MakePackDescriptor(mc->GetMeshPointer()->resolved_packages.get(), {field}); @@ -247,8 +246,7 @@ void SecondDerivative(const AMRBounds &bnds, MeshData *mc, const std::stri const Real refine_criteria = refine_criteria_; const Real derefine_criteria = derefine_criteria_; auto delta_levels = delta_levels_; - // TODO(?): I have no idea how to get the index from the tensor indices... - const int var = index[2]; + const int var = idx; par_for_outer( PARTHENON_AUTO_LABEL, 0, 0, 0, pack.GetNBlocks() - 1, KOKKOS_LAMBDA(team_mbr_t team_member, const int b) { diff --git a/src/amr_criteria/refinement_package.hpp b/src/amr_criteria/refinement_package.hpp index e0c1b2329283..781157758d80 100644 --- a/src/amr_criteria/refinement_package.hpp +++ b/src/amr_criteria/refinement_package.hpp @@ -46,14 +46,14 @@ AmrTag FirstDerivative(const AMRBounds &bnds, const ParArray3D &q, const Real refine_criteria, const Real derefine_criteria); void FirstDerivative(const AMRBounds &bnds, MeshData *mc, const std::string &field, - Kokkos::Array index, ParArray1D &delta_levels_, + const int &idx, ParArray1D &delta_levels_, const Real refine_criteria_, const Real derefine_criteria_); AmrTag SecondDerivative(const AMRBounds &bnds, const ParArray3D &q, const Real refine_criteria, const Real derefine_criteria); void SecondDerivative(const AMRBounds &bnds, MeshData *mc, const std::string &field, - Kokkos::Array index, ParArray1D &delta_levels_, + const int &idx, ParArray1D &delta_levels_, const Real refine_criteria_, const Real derefine_criteria_); } // namespace Refinement From 5cb25b44222c72456666e09d67e48874e3e9cdaa Mon Sep 17 00:00:00 2001 From: adam reyes Date: Fri, 27 Sep 2024 23:45:29 +0200 Subject: [PATCH 04/49] added CheckRefinementMesh to fine-advection example --- example/fine_advection/advection_package.cpp | 56 +++++++++++++++++++- example/fine_advection/advection_package.hpp | 1 + 2 files changed, 56 insertions(+), 1 deletion(-) diff --git a/example/fine_advection/advection_package.cpp b/example/fine_advection/advection_package.cpp index 8ca6c6458b20..3a50fbd8f1dd 100644 --- a/example/fine_advection/advection_package.cpp +++ b/example/fine_advection/advection_package.cpp @@ -94,12 +94,66 @@ std::shared_ptr Initialize(ParameterInput *pin) { pkg->AddField( Metadata({Metadata::Cell, Metadata::Derived, Metadata::OneCopy})); - pkg->CheckRefinementBlock = CheckRefinement; + bool check_refine_mesh = + pin->GetOrAddBoolean("parthenon/mesh", "CheckRefineMesh", false); + if (check_refine_mesh) { + pkg->CheckRefinementMesh = CheckRefinementMesh; + } else { + pkg->CheckRefinementBlock = CheckRefinement; + } pkg->EstimateTimestepMesh = EstimateTimestep; pkg->FillDerivedMesh = FillDerived; return pkg; } +void CheckRefinementMesh(MeshData *md, parthenon::ParArray1D &data_levels) { + // refine on advected, for example. could also be a derived quantity + static auto desc = parthenon::MakePackDescriptor(md); + auto pack = desc.GetPack(md); + + auto pkg = md->GetMeshPointer()->packages.Get("advection_package"); + const auto &refine_tol = pkg->Param("refine_tol"); + const auto &derefine_tol = pkg->Param("derefine_tol"); + + auto ib = md->GetBoundsI(IndexDomain::entire); + auto jb = md->GetBoundsJ(IndexDomain::entire); + auto kb = md->GetBoundsK(IndexDomain::entire); + parthenon::par_for_outer( + PARTHENON_AUTO_LABEL, 0, 0, 0, pack.GetNBlocks() - 1, + KOKKOS_LAMBDA(parthenon::team_mbr_t team_member, const int b) { + if (data_levels(b) == AmrTag::refine) return; + using MinMax_t = typename Kokkos::MinMax::value_type; + MinMax_t minmax, temp_minmax; + minmax.max_val = 0.; + minmax.min_val = 0.; + for (int n = pack.GetLowerBound(b); n <= pack.GetUpperBound(b); n++) { + parthenon::par_reduce_inner( + parthenon::inner_loop_pattern_ttr_tag, team_member, kb.s, kb.e, jb.s, jb.e, + ib.s, ib.e, + [&](const int k, const int j, const int i, MinMax_t &lminmax) { + lminmax.min_val = + (pack(b, n, k, j, i) < lminmax.min_val ? pack(b, n, k, j, i) + : lminmax.min_val); + lminmax.max_val = + (pack(b, n, k, j, i) > lminmax.max_val ? pack(b, n, k, j, i) + : lminmax.max_val); + }, + Kokkos::MinMax(temp_minmax)); + minmax.min_val = + temp_minmax.min_val < minmax.min_val ? temp_minmax.min_val : minmax.min_val; + minmax.max_val = + temp_minmax.max_val > minmax.max_val ? temp_minmax.max_val : minmax.max_val; + } + if (minmax.max_val > refine_tol && minmax.min_val < derefine_tol) { + data_levels(b) = AmrTag::refine; + } else if (minmax.max_val < derefine_tol) { + data_levels(b) = AmrTag::derefine; + } else { + data_levels(b) = AmrTag::same; + } + }); +} + AmrTag CheckRefinement(MeshBlockData *rc) { // refine on advected, for example. could also be a derived quantity static auto desc = parthenon::MakePackDescriptor(rc); diff --git a/example/fine_advection/advection_package.hpp b/example/fine_advection/advection_package.hpp index 87e0c0b2edab..42113d2fc090 100644 --- a/example/fine_advection/advection_package.hpp +++ b/example/fine_advection/advection_package.hpp @@ -48,6 +48,7 @@ VARIABLE(advection, divD); std::shared_ptr Initialize(ParameterInput *pin); AmrTag CheckRefinement(MeshBlockData *rc); +void CheckRefinementMesh(MeshData *md, parthenon::ParArray1D &data_levels); Real EstimateTimestep(MeshData *md); TaskStatus FillDerived(MeshData *md); From 3ad74e17baa0e1e6834aef781416cec7454a9adc Mon Sep 17 00:00:00 2001 From: adam reyes Date: Sat, 28 Sep 2024 14:40:45 +0200 Subject: [PATCH 05/49] cleaning up var names --- example/fine_advection/advection_package.cpp | 17 ++++++++--------- example/fine_advection/advection_package.hpp | 2 +- 2 files changed, 9 insertions(+), 10 deletions(-) diff --git a/example/fine_advection/advection_package.cpp b/example/fine_advection/advection_package.cpp index 3a50fbd8f1dd..f790dbedbc58 100644 --- a/example/fine_advection/advection_package.cpp +++ b/example/fine_advection/advection_package.cpp @@ -106,7 +106,8 @@ std::shared_ptr Initialize(ParameterInput *pin) { return pkg; } -void CheckRefinementMesh(MeshData *md, parthenon::ParArray1D &data_levels) { +void CheckRefinementMesh(MeshData *md, + parthenon::ParArray1D &delta_levels) { // refine on advected, for example. could also be a derived quantity static auto desc = parthenon::MakePackDescriptor(md); auto pack = desc.GetPack(md); @@ -121,7 +122,7 @@ void CheckRefinementMesh(MeshData *md, parthenon::ParArray1D &data parthenon::par_for_outer( PARTHENON_AUTO_LABEL, 0, 0, 0, pack.GetNBlocks() - 1, KOKKOS_LAMBDA(parthenon::team_mbr_t team_member, const int b) { - if (data_levels(b) == AmrTag::refine) return; + if (delta_levels(b) == AmrTag::refine) return; using MinMax_t = typename Kokkos::MinMax::value_type; MinMax_t minmax, temp_minmax; minmax.max_val = 0.; @@ -144,13 +145,11 @@ void CheckRefinementMesh(MeshData *md, parthenon::ParArray1D &data minmax.max_val = temp_minmax.max_val > minmax.max_val ? temp_minmax.max_val : minmax.max_val; } - if (minmax.max_val > refine_tol && minmax.min_val < derefine_tol) { - data_levels(b) = AmrTag::refine; - } else if (minmax.max_val < derefine_tol) { - data_levels(b) = AmrTag::derefine; - } else { - data_levels(b) = AmrTag::same; - } + auto flag = AmrTag::same; + if (minmax.max_val > refine_tol && minmax.min_val < derefine_tol) + flag = AmrTag::refine; + if (minmax.max_val < derefine_tol) flag = AmrTag::derefine; + delta_levels(b) = std::max(delta_levels(b), flag); }); } diff --git a/example/fine_advection/advection_package.hpp b/example/fine_advection/advection_package.hpp index 42113d2fc090..d444cf29382a 100644 --- a/example/fine_advection/advection_package.hpp +++ b/example/fine_advection/advection_package.hpp @@ -48,7 +48,7 @@ VARIABLE(advection, divD); std::shared_ptr Initialize(ParameterInput *pin); AmrTag CheckRefinement(MeshBlockData *rc); -void CheckRefinementMesh(MeshData *md, parthenon::ParArray1D &data_levels); +void CheckRefinementMesh(MeshData *md, parthenon::ParArray1D &delta_levels); Real EstimateTimestep(MeshData *md); TaskStatus FillDerived(MeshData *md); From 16a1c5895e237089254fe04ee65510535e701870 Mon Sep 17 00:00:00 2001 From: adam reyes Date: Sat, 28 Sep 2024 14:55:08 +0200 Subject: [PATCH 06/49] cleanup --- src/amr_criteria/amr_criteria.hpp | 1 - src/amr_criteria/refinement_package.cpp | 7 ++++--- src/amr_criteria/refinement_package.hpp | 1 - 3 files changed, 4 insertions(+), 5 deletions(-) diff --git a/src/amr_criteria/amr_criteria.hpp b/src/amr_criteria/amr_criteria.hpp index 411646e9ac9e..a57ebb1ebb9d 100644 --- a/src/amr_criteria/amr_criteria.hpp +++ b/src/amr_criteria/amr_criteria.hpp @@ -15,7 +15,6 @@ #include #include -#include #include "defs.hpp" #include "mesh/domain.hpp" diff --git a/src/amr_criteria/refinement_package.cpp b/src/amr_criteria/refinement_package.cpp index 7f57f857a3dd..eddd8ff4e2d4 100644 --- a/src/amr_criteria/refinement_package.cpp +++ b/src/amr_criteria/refinement_package.cpp @@ -18,9 +18,9 @@ #include #include #include -#include -#include "Kokkos_Macros.hpp" +#include + #include "amr_criteria/amr_criteria.hpp" #include "interface/make_pack_descriptor.hpp" #include "interface/mesh_data.hpp" @@ -96,7 +96,8 @@ AmrTag CheckAllRefinement(MeshBlockData *rc, const AmrTag &level) { MeshBlock *pmb = rc->GetBlockPointer(); static const bool check_refine_mesh = pmb->packages.Get("Refinement")->Param("check_refine_mesh"); - // delta_level holds the max over all criteria. default to derefining. + // delta_level holds the max over all criteria. default to derefining, or level from + // MeshData check. AmrTag delta_level = level; for (auto &pkg : pmb->packages.AllPackages()) { auto &desc = pkg.second; diff --git a/src/amr_criteria/refinement_package.hpp b/src/amr_criteria/refinement_package.hpp index 781157758d80..e8f100d01686 100644 --- a/src/amr_criteria/refinement_package.hpp +++ b/src/amr_criteria/refinement_package.hpp @@ -16,7 +16,6 @@ #include #include -#include #include "defs.hpp" #include "parthenon_arrays.hpp" From 6e1b9056651fd5015fbec2cac0a6def128c28752 Mon Sep 17 00:00:00 2001 From: adam reyes Date: Sat, 28 Sep 2024 23:24:23 +0200 Subject: [PATCH 07/49] adding scatter view utilities --- src/parthenon_array_generic.hpp | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/src/parthenon_array_generic.hpp b/src/parthenon_array_generic.hpp index d527707f9070..5328ae7b73f6 100644 --- a/src/parthenon_array_generic.hpp +++ b/src/parthenon_array_generic.hpp @@ -21,6 +21,8 @@ #include #include +#include + #include "utils/concepts_lite.hpp" namespace parthenon { @@ -214,6 +216,24 @@ class ParArrayGeneric : public State { KOKKOS_INLINE_FUNCTION auto size() const { return data_.size(); } + // utilities for scatter views + template + auto ToScatterView() { + using view_type = std::remove_cv_t>; + using data_type = typename view_type::data_type; + using exec_space = typename view_type::execution_space; + using layout = typename view_type::array_layout; + return Kokkos::Experimental::ScatterView(data_); + } + + template + void ContributeScatter(ScatterView_t scatter) { + static_assert( + is_specialization_of::value, + "Need to provide a Kokkos::Experimental::ScatterView"); + Kokkos::Experimental::contribute(data_, scatter); + } + // a function to get the total size of the array KOKKOS_INLINE_FUNCTION int GetSize() const { return data_.size(); From 2b6545c65ec524418a5bf67ae347beabe2b2b3bf Mon Sep 17 00:00:00 2001 From: adam reyes Date: Sat, 28 Sep 2024 23:24:34 +0200 Subject: [PATCH 08/49] scatterview version of refinement --- example/fine_advection/advection_package.cpp | 44 +++---- src/amr_criteria/refinement_package.cpp | 121 +++++++++---------- 2 files changed, 68 insertions(+), 97 deletions(-) diff --git a/example/fine_advection/advection_package.cpp b/example/fine_advection/advection_package.cpp index f790dbedbc58..ffae5d11938b 100644 --- a/example/fine_advection/advection_package.cpp +++ b/example/fine_advection/advection_package.cpp @@ -17,8 +17,10 @@ #include #include #include +#include #include +#include #include #include #include @@ -28,6 +30,7 @@ #include "kokkos_abstraction.hpp" #include "reconstruct/dc_inline.hpp" #include "utils/error_checking.hpp" +#include "utils/instrument.hpp" using namespace parthenon::package::prelude; @@ -119,38 +122,19 @@ void CheckRefinementMesh(MeshData *md, auto ib = md->GetBoundsI(IndexDomain::entire); auto jb = md->GetBoundsJ(IndexDomain::entire); auto kb = md->GetBoundsK(IndexDomain::entire); - parthenon::par_for_outer( - PARTHENON_AUTO_LABEL, 0, 0, 0, pack.GetNBlocks() - 1, - KOKKOS_LAMBDA(parthenon::team_mbr_t team_member, const int b) { - if (delta_levels(b) == AmrTag::refine) return; - using MinMax_t = typename Kokkos::MinMax::value_type; - MinMax_t minmax, temp_minmax; - minmax.max_val = 0.; - minmax.min_val = 0.; - for (int n = pack.GetLowerBound(b); n <= pack.GetUpperBound(b); n++) { - parthenon::par_reduce_inner( - parthenon::inner_loop_pattern_ttr_tag, team_member, kb.s, kb.e, jb.s, jb.e, - ib.s, ib.e, - [&](const int k, const int j, const int i, MinMax_t &lminmax) { - lminmax.min_val = - (pack(b, n, k, j, i) < lminmax.min_val ? pack(b, n, k, j, i) - : lminmax.min_val); - lminmax.max_val = - (pack(b, n, k, j, i) > lminmax.max_val ? pack(b, n, k, j, i) - : lminmax.max_val); - }, - Kokkos::MinMax(temp_minmax)); - minmax.min_val = - temp_minmax.min_val < minmax.min_val ? temp_minmax.min_val : minmax.min_val; - minmax.max_val = - temp_minmax.max_val > minmax.max_val ? temp_minmax.max_val : minmax.max_val; - } + auto scatter_levels = delta_levels.ToScatterView(); + parthenon::par_for( + PARTHENON_AUTO_LABEL, 0, pack.GetNBlocks() - 1, 0, pack.GetMaxNumberOfVars() - 1, + kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int b, const int n, const int k, const int j, const int i) { + auto levels_access = scatter_levels.access(); auto flag = AmrTag::same; - if (minmax.max_val > refine_tol && minmax.min_val < derefine_tol) - flag = AmrTag::refine; - if (minmax.max_val < derefine_tol) flag = AmrTag::derefine; - delta_levels(b) = std::max(delta_levels(b), flag); + const auto val = pack(b, n, k, j, i); + if (val > refine_tol) flag = AmrTag::refine; + if (val < derefine_tol) flag = AmrTag::derefine; + levels_access(b).update(flag); }); + delta_levels.ContributeScatter(scatter_levels); } AmrTag CheckRefinement(MeshBlockData *rc) { diff --git a/src/amr_criteria/refinement_package.cpp b/src/amr_criteria/refinement_package.cpp index eddd8ff4e2d4..20160ba5c2e4 100644 --- a/src/amr_criteria/refinement_package.cpp +++ b/src/amr_criteria/refinement_package.cpp @@ -21,6 +21,7 @@ #include +#include "Kokkos_ScatterView.hpp" #include "amr_criteria/amr_criteria.hpp" #include "interface/make_pack_descriptor.hpp" #include "interface/mesh_data.hpp" @@ -31,6 +32,7 @@ #include "mesh/mesh_refinement.hpp" #include "mesh/meshblock.hpp" #include "parameter_input.hpp" +#include "utils/instrument.hpp" namespace parthenon { namespace Refinement { @@ -187,7 +189,7 @@ AmrTag SecondDerivative(const AMRBounds &bnds, const ParArray3D &q, } void FirstDerivative(const AMRBounds &bnds, MeshData *mc, const std::string &field, - const int &idx, ParArray1D &delta_levels_, + const int &idx, ParArray1D &delta_levels, const Real refine_criteria_, const Real derefine_criteria_) { const auto desc = MakePackDescriptor(mc->GetMeshPointer()->resolved_packages.get(), {field}); @@ -197,46 +199,37 @@ void FirstDerivative(const AMRBounds &bnds, MeshData *mc, const std::strin const Real refine_criteria = refine_criteria_; const Real derefine_criteria = derefine_criteria_; - auto delta_levels = delta_levels_; const int var = idx; - par_for_outer( - PARTHENON_AUTO_LABEL, 0, 0, 0, pack.GetNBlocks() - 1, - KOKKOS_LAMBDA(team_mbr_t team_member, const int b) { - if (delta_levels(b) == AmrTag::refine) return; - Real maxd = 0.; - par_reduce_inner( - inner_loop_pattern_ttr_tag, team_member, bnds.ks, bnds.ke, bnds.js, bnds.je, - bnds.is, bnds.ie, - [&](const int k, const int j, const int i, Real &maxder) { - Real scale = std::abs(pack(b, var, k, j, i)); - Real d = 0.5 * - std::abs((pack(b, var, k, j, i + 1) - pack(b, var, k, j, i - 1))) / - (scale + TINY_NUMBER); - maxder = (d > maxder ? d : maxder); - if (ndim > 1) { - d = 0.5 * - std::abs((pack(b, var, k, j + 1, i) - pack(b, var, k, j - 1, i))) / - (scale + TINY_NUMBER); - maxder = (d > maxder ? d : maxder); - } - if (ndim > 2) { - d = 0.5 * - std::abs((pack(b, var, k + 1, j, i) - pack(b, var, k - 1, j, i))) / - (scale + TINY_NUMBER); - maxder = (d > maxder ? d : maxder); - } - }, - Kokkos::Max(maxd)); - - AmrTag flag = AmrTag::same; - if (maxd > refine_criteria) flag = AmrTag::refine; - if (maxd < derefine_criteria) flag = AmrTag::derefine; - delta_levels(b) = std::max(delta_levels(b), flag); + auto scatter_levels = delta_levels.ToScatterView(); + par_for( + PARTHENON_AUTO_LABEL, 0, pack.GetNBlocks() - 1, bnds.ks, bnds.ke, bnds.js, bnds.je, + bnds.is, bnds.ie, + KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) { + Real scale = std::abs(pack(b, var, k, j, i)); + Real d = 0.5 * std::abs((pack(b, var, k, j, i + 1) - pack(b, var, k, j, i - 1))) / + (scale + TINY_NUMBER); + Real maxder = d; + if (ndim > 1) { + d = 0.5 * std::abs((pack(b, var, k, j + 1, i) - pack(b, var, k, j - 1, i))) / + (scale + TINY_NUMBER); + maxder = (d > maxder ? d : maxder); + } + if (ndim > 2) { + d = 0.5 * std::abs((pack(b, var, k + 1, j, i) - pack(b, var, k - 1, j, i))) / + (scale + TINY_NUMBER); + maxder = (d > maxder ? d : maxder); + } + auto levels_access = scatter_levels.access(); + auto flag = AmrTag::same; + if (maxder > refine_criteria) flag = AmrTag::refine; + if (maxder < derefine_criteria) flag = AmrTag::derefine; + levels_access(b).update(flag); }); + delta_levels.ContributeScatter(scatter_levels); } void SecondDerivative(const AMRBounds &bnds, MeshData *mc, const std::string &field, - const int &idx, ParArray1D &delta_levels_, + const int &idx, ParArray1D &delta_levels, const Real refine_criteria_, const Real derefine_criteria_) { const auto desc = MakePackDescriptor(mc->GetMeshPointer()->resolved_packages.get(), {field}); @@ -246,39 +239,33 @@ void SecondDerivative(const AMRBounds &bnds, MeshData *mc, const std::stri const Real refine_criteria = refine_criteria_; const Real derefine_criteria = derefine_criteria_; - auto delta_levels = delta_levels_; const int var = idx; - par_for_outer( - PARTHENON_AUTO_LABEL, 0, 0, 0, pack.GetNBlocks() - 1, - KOKKOS_LAMBDA(team_mbr_t team_member, const int b) { - if (delta_levels(b) == AmrTag::refine) return; - Real maxd = 0.; - par_reduce_inner( - inner_loop_pattern_ttr_tag, team_member, bnds.ks, bnds.ke, bnds.js, bnds.je, - bnds.is, bnds.ie, - [&](const int k, const int j, const int i, Real &maxder) { - Real aqt = std::abs(pack(b, var, k, j, i)) + TINY_NUMBER; - Real qavg = 0.5 * (pack(b, var, k, j, i + 1) + pack(b, var, k, j, i - 1)); - Real d = std::abs(qavg - pack(b, var, k, j, i)) / (std::abs(qavg) + aqt); - maxder = (d > maxder ? d : maxder); - if (ndim > 1) { - qavg = 0.5 * (pack(b, var, k, j + 1, i) + pack(b, var, k, j - 1, i)); - d = std::abs(qavg - pack(b, var, k, j, i)) / (std::abs(qavg) + aqt); - maxder = (d > maxder ? d : maxder); - } - if (ndim > 2) { - qavg = 0.5 * (pack(b, var, k + 1, j, i) + pack(b, var, k - 1, j, i)); - d = std::abs(qavg - pack(b, var, k, j, i)) / (std::abs(qavg) + aqt); - maxder = (d > maxder ? d : maxder); - } - }, - Kokkos::Max(maxd)); - - AmrTag flag = AmrTag::same; - if (maxd > refine_criteria) flag = AmrTag::refine; - if (maxd < derefine_criteria) flag = AmrTag::derefine; - delta_levels(b) = std::max(delta_levels(b), flag); + auto scatter_levels = delta_levels.ToScatterView(); + par_for( + PARTHENON_AUTO_LABEL, 0, pack.GetNBlocks() - 1, bnds.ks, bnds.ke, bnds.js, bnds.je, + bnds.is, bnds.ie, + KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) { + Real aqt = std::abs(pack(b, var, k, j, i)) + TINY_NUMBER; + Real qavg = 0.5 * (pack(b, var, k, j, i + 1) + pack(b, var, k, j, i - 1)); + Real d = std::abs(qavg - pack(b, var, k, j, i)) / (std::abs(qavg) + aqt); + Real maxder = d; + if (ndim > 1) { + qavg = 0.5 * (pack(b, var, k, j + 1, i) + pack(b, var, k, j - 1, i)); + d = std::abs(qavg - pack(b, var, k, j, i)) / (std::abs(qavg) + aqt); + maxder = (d > maxder ? d : maxder); + } + if (ndim > 2) { + qavg = 0.5 * (pack(b, var, k + 1, j, i) + pack(b, var, k - 1, j, i)); + d = std::abs(qavg - pack(b, var, k, j, i)) / (std::abs(qavg) + aqt); + maxder = (d > maxder ? d : maxder); + } + auto levels_access = scatter_levels.access(); + auto flag = AmrTag::same; + if (maxder > refine_criteria) flag = AmrTag::refine; + if (maxder < derefine_criteria) flag = AmrTag::derefine; + levels_access(b).update(flag); }); + delta_levels.ContributeScatter(scatter_levels); } void SetRefinement_(MeshBlockData *rc, From ccf72b0472c8e8fac5dfa6bb7f62e49a0eed61ff Mon Sep 17 00:00:00 2001 From: adam reyes Date: Sun, 29 Sep 2024 13:00:44 +0200 Subject: [PATCH 09/49] burgers-benchmark uses Tag --- benchmarks/burgers/burgers_driver.cpp | 10 ++-------- 1 file changed, 2 insertions(+), 8 deletions(-) diff --git a/benchmarks/burgers/burgers_driver.cpp b/benchmarks/burgers/burgers_driver.cpp index cddb255a636b..210f0276c99f 100644 --- a/benchmarks/burgers/burgers_driver.cpp +++ b/benchmarks/burgers/burgers_driver.cpp @@ -123,6 +123,8 @@ TaskCollection BurgersDriver::MakeTaskCollection(BlockList_t &blocks, const int // estimate next time step if (stage == integrator->nstages) { auto new_dt = tl.AddTask(update, EstimateTimestep>, mc1.get()); + auto tag_refine = + tl.AddTask(update, parthenon::Refinement::Tag>, mc1.get()); } } @@ -135,14 +137,6 @@ TaskCollection BurgersDriver::MakeTaskCollection(BlockList_t &blocks, const int // set physical boundaries auto set_bc = tl.AddTask(none, parthenon::ApplyBoundaryConditions, sc1); - - if (stage == integrator->nstages) { - // Update refinement - if (pmesh->adaptive) { - auto tag_refine = tl.AddTask( - set_bc, parthenon::Refinement::Tag>, sc1.get()); - } - } } return tc; } From d7d536fcb2dbef96ca20365aec1154e661254f99 Mon Sep 17 00:00:00 2001 From: adam reyes Date: Sun, 29 Sep 2024 14:29:47 +0200 Subject: [PATCH 10/49] add hierarchial par --- example/fine_advection/advection_package.cpp | 27 ++++-- src/amr_criteria/refinement_package.cpp | 95 +++++++++++--------- 2 files changed, 75 insertions(+), 47 deletions(-) diff --git a/example/fine_advection/advection_package.cpp b/example/fine_advection/advection_package.cpp index ffae5d11938b..e2bcecf84f72 100644 --- a/example/fine_advection/advection_package.cpp +++ b/example/fine_advection/advection_package.cpp @@ -123,15 +123,28 @@ void CheckRefinementMesh(MeshData *md, auto jb = md->GetBoundsJ(IndexDomain::entire); auto kb = md->GetBoundsK(IndexDomain::entire); auto scatter_levels = delta_levels.ToScatterView(); - parthenon::par_for( - PARTHENON_AUTO_LABEL, 0, pack.GetNBlocks() - 1, 0, pack.GetMaxNumberOfVars() - 1, - kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, - KOKKOS_LAMBDA(const int b, const int n, const int k, const int j, const int i) { + parthenon::par_for_outer( + PARTHENON_AUTO_LABEL, 0, 0, 0, pack.GetNBlocks() - 1, 0, + pack.GetMaxNumberOfVars() - 1, kb.s, kb.e, + KOKKOS_LAMBDA(parthenon::team_mbr_t team_member, const int b, const int n, + const int k) { + typename Kokkos::MinMax::value_type minmax; + par_reduce_inner( + parthenon::inner_loop_pattern_ttr_tag, team_member, jb.s, jb.e, ib.s, ib.e, + [&](const int j, const int i, + typename Kokkos::MinMax::value_type &lminmax) { + lminmax.min_val = (pack(n, k, j, i) < lminmax.min_val ? pack(n, k, j, i) + : lminmax.min_val); + lminmax.max_val = (pack(n, k, j, i) > lminmax.max_val ? pack(n, k, j, i) + : lminmax.max_val); + }, + Kokkos::MinMax(minmax)); + auto levels_access = scatter_levels.access(); auto flag = AmrTag::same; - const auto val = pack(b, n, k, j, i); - if (val > refine_tol) flag = AmrTag::refine; - if (val < derefine_tol) flag = AmrTag::derefine; + if (minmax.max_val > refine_tol && minmax.min_val < derefine_tol) + flag = AmrTag::refine; + if (minmax.max_val < derefine_tol) flag = AmrTag::derefine; levels_access(b).update(flag); }); delta_levels.ContributeScatter(scatter_levels); diff --git a/src/amr_criteria/refinement_package.cpp b/src/amr_criteria/refinement_package.cpp index 20160ba5c2e4..4648a94280cd 100644 --- a/src/amr_criteria/refinement_package.cpp +++ b/src/amr_criteria/refinement_package.cpp @@ -201,28 +201,37 @@ void FirstDerivative(const AMRBounds &bnds, MeshData *mc, const std::strin const Real derefine_criteria = derefine_criteria_; const int var = idx; auto scatter_levels = delta_levels.ToScatterView(); - par_for( - PARTHENON_AUTO_LABEL, 0, pack.GetNBlocks() - 1, bnds.ks, bnds.ke, bnds.js, bnds.je, - bnds.is, bnds.ie, - KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) { - Real scale = std::abs(pack(b, var, k, j, i)); - Real d = 0.5 * std::abs((pack(b, var, k, j, i + 1) - pack(b, var, k, j, i - 1))) / - (scale + TINY_NUMBER); - Real maxder = d; - if (ndim > 1) { - d = 0.5 * std::abs((pack(b, var, k, j + 1, i) - pack(b, var, k, j - 1, i))) / - (scale + TINY_NUMBER); - maxder = (d > maxder ? d : maxder); - } - if (ndim > 2) { - d = 0.5 * std::abs((pack(b, var, k + 1, j, i) - pack(b, var, k - 1, j, i))) / - (scale + TINY_NUMBER); - maxder = (d > maxder ? d : maxder); - } + par_for_outer( + PARTHENON_AUTO_LABEL, 0, 0, 0, pack.GetNBlocks() - 1, bnds.ks, bnds.ke, bnds.js, + bnds.je, + KOKKOS_LAMBDA(team_mbr_t team_member, const int b, const int k, const int j) { + Real maxd = 0.; + par_reduce_inner( + inner_loop_pattern_ttr_tag, team_member, bnds.is, bnds.ie, + [&](const int i, Real &maxder) { + Real scale = std::abs(pack(b, var, k, j, i)); + Real d = 0.5 * + std::abs((pack(b, var, k, j, i + 1) - pack(b, var, k, j, i - 1))) / + (scale + TINY_NUMBER); + maxder = (d > maxder ? d : maxder); + if (ndim > 1) { + d = 0.5 * + std::abs((pack(b, var, k, j + 1, i) - pack(b, var, k, j - 1, i))) / + (scale + TINY_NUMBER); + maxder = (d > maxder ? d : maxder); + } + if (ndim > 2) { + d = 0.5 * + std::abs((pack(b, var, k + 1, j, i) - pack(b, var, k - 1, j, i))) / + (scale + TINY_NUMBER); + maxder = (d > maxder ? d : maxder); + } + }, + Kokkos::Max(maxd)); auto levels_access = scatter_levels.access(); auto flag = AmrTag::same; - if (maxder > refine_criteria) flag = AmrTag::refine; - if (maxder < derefine_criteria) flag = AmrTag::derefine; + if (maxd > refine_criteria) flag = AmrTag::refine; + if (maxd < derefine_criteria) flag = AmrTag::derefine; levels_access(b).update(flag); }); delta_levels.ContributeScatter(scatter_levels); @@ -241,28 +250,34 @@ void SecondDerivative(const AMRBounds &bnds, MeshData *mc, const std::stri const Real derefine_criteria = derefine_criteria_; const int var = idx; auto scatter_levels = delta_levels.ToScatterView(); - par_for( - PARTHENON_AUTO_LABEL, 0, pack.GetNBlocks() - 1, bnds.ks, bnds.ke, bnds.js, bnds.je, - bnds.is, bnds.ie, - KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) { - Real aqt = std::abs(pack(b, var, k, j, i)) + TINY_NUMBER; - Real qavg = 0.5 * (pack(b, var, k, j, i + 1) + pack(b, var, k, j, i - 1)); - Real d = std::abs(qavg - pack(b, var, k, j, i)) / (std::abs(qavg) + aqt); - Real maxder = d; - if (ndim > 1) { - qavg = 0.5 * (pack(b, var, k, j + 1, i) + pack(b, var, k, j - 1, i)); - d = std::abs(qavg - pack(b, var, k, j, i)) / (std::abs(qavg) + aqt); - maxder = (d > maxder ? d : maxder); - } - if (ndim > 2) { - qavg = 0.5 * (pack(b, var, k + 1, j, i) + pack(b, var, k - 1, j, i)); - d = std::abs(qavg - pack(b, var, k, j, i)) / (std::abs(qavg) + aqt); - maxder = (d > maxder ? d : maxder); - } + par_for_outer( + PARTHENON_AUTO_LABEL, 0, 0, 0, pack.GetNBlocks() - 1, bnds.ks, bnds.ke, bnds.js, + bnds.je, + KOKKOS_LAMBDA(team_mbr_t team_member, const int b, const int k, const int j) { + Real maxd = 0.; + par_reduce_inner( + inner_loop_pattern_ttr_tag, team_member, bnds.is, bnds.ie, + [&](const int i, Real &maxder) { + Real aqt = std::abs(pack(b, var, k, j, i)) + TINY_NUMBER; + Real qavg = 0.5 * (pack(b, var, k, j, i + 1) + pack(b, var, k, j, i - 1)); + Real d = std::abs(qavg - pack(b, var, k, j, i)) / (std::abs(qavg) + aqt); + maxder = (d > maxder ? d : maxder); + if (ndim > 1) { + qavg = 0.5 * (pack(b, var, k, j + 1, i) + pack(b, var, k, j - 1, i)); + d = std::abs(qavg - pack(b, var, k, j, i)) / (std::abs(qavg) + aqt); + maxder = (d > maxder ? d : maxder); + } + if (ndim > 2) { + qavg = 0.5 * (pack(b, var, k + 1, j, i) + pack(b, var, k - 1, j, i)); + d = std::abs(qavg - pack(b, var, k, j, i)) / (std::abs(qavg) + aqt); + maxder = (d > maxder ? d : maxder); + } + }, + Kokkos::Max(maxd)); auto levels_access = scatter_levels.access(); auto flag = AmrTag::same; - if (maxder > refine_criteria) flag = AmrTag::refine; - if (maxder < derefine_criteria) flag = AmrTag::derefine; + if (maxd > refine_criteria) flag = AmrTag::refine; + if (maxd < derefine_criteria) flag = AmrTag::derefine; levels_access(b).update(flag); }); delta_levels.ContributeScatter(scatter_levels); From fa37fb92de7ab29c654926e02cce5066ec497975 Mon Sep 17 00:00:00 2001 From: adam reyes Date: Sun, 29 Sep 2024 16:16:12 +0200 Subject: [PATCH 11/49] cleanup includes --- example/fine_advection/advection_package.cpp | 1 - src/amr_criteria/refinement_package.cpp | 3 --- 2 files changed, 4 deletions(-) diff --git a/example/fine_advection/advection_package.cpp b/example/fine_advection/advection_package.cpp index e2bcecf84f72..a2b30757002f 100644 --- a/example/fine_advection/advection_package.cpp +++ b/example/fine_advection/advection_package.cpp @@ -17,7 +17,6 @@ #include #include #include -#include #include #include diff --git a/src/amr_criteria/refinement_package.cpp b/src/amr_criteria/refinement_package.cpp index 4648a94280cd..c8a12e2ed657 100644 --- a/src/amr_criteria/refinement_package.cpp +++ b/src/amr_criteria/refinement_package.cpp @@ -19,9 +19,6 @@ #include #include -#include - -#include "Kokkos_ScatterView.hpp" #include "amr_criteria/amr_criteria.hpp" #include "interface/make_pack_descriptor.hpp" #include "interface/mesh_data.hpp" From 401f4123951e8d9e8bee360def659cc76e2a3e26 Mon Sep 17 00:00:00 2001 From: adam reyes Date: Sun, 29 Sep 2024 16:19:21 +0200 Subject: [PATCH 12/49] missed one --- example/fine_advection/advection_package.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/example/fine_advection/advection_package.cpp b/example/fine_advection/advection_package.cpp index a2b30757002f..94a40f6cb130 100644 --- a/example/fine_advection/advection_package.cpp +++ b/example/fine_advection/advection_package.cpp @@ -19,7 +19,6 @@ #include #include -#include #include #include #include @@ -29,7 +28,6 @@ #include "kokkos_abstraction.hpp" #include "reconstruct/dc_inline.hpp" #include "utils/error_checking.hpp" -#include "utils/instrument.hpp" using namespace parthenon::package::prelude; From 638e0773ef4a8385824fa2acd7bfced425dd0e51 Mon Sep 17 00:00:00 2001 From: Adam C Reyes Date: Sun, 29 Sep 2024 18:05:45 +0200 Subject: [PATCH 13/49] Update CHANGELOG.md --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index a7373b19036b..efce2aa0cbec 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -3,6 +3,7 @@ ## Current develop ### Added (new features/APIs/variables/...) +- [[PR 1182]](https://github.com/parthenon-hpc-lab/parthenon/pull/1182) Add a MeshData variant for refinement tagging - [[PR 1179]](https://github.com/parthenon-hpc-lab/parthenon/pull/1179) Make a global variable for whether simulation is a restart - [[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 From b6e0e6ba7cb72cc68ad9bdfcc0bdb5eb5165cefd Mon Sep 17 00:00:00 2001 From: adam reyes Date: Sun, 29 Sep 2024 20:37:26 +0200 Subject: [PATCH 14/49] remove default level tag --- src/amr_criteria/refinement_package.cpp | 2 +- src/amr_criteria/refinement_package.hpp | 3 +-- src/mesh/mesh_refinement.cpp | 11 ----------- 3 files changed, 2 insertions(+), 14 deletions(-) diff --git a/src/amr_criteria/refinement_package.cpp b/src/amr_criteria/refinement_package.cpp index c8a12e2ed657..78d15c689e02 100644 --- a/src/amr_criteria/refinement_package.cpp +++ b/src/amr_criteria/refinement_package.cpp @@ -296,7 +296,7 @@ TaskStatus Tag(MeshBlockData *rc) { template <> TaskStatus Tag(MeshData *rc) { PARTHENON_INSTRUMENT - ParArray1D delta_levels = CheckAllRefinement(rc); + ParArray1D delta_levels = CheckAllRefinement(rc, AmrTag::derefine); auto delta_levels_h = delta_levels.GetHostMirrorAndCopy(); for (int i = 0; i < rc->NumBlocks(); i++) { diff --git a/src/amr_criteria/refinement_package.hpp b/src/amr_criteria/refinement_package.hpp index e8f100d01686..444e1c5e3e32 100644 --- a/src/amr_criteria/refinement_package.hpp +++ b/src/amr_criteria/refinement_package.hpp @@ -37,8 +37,7 @@ std::shared_ptr Initialize(ParameterInput *pin); template TaskStatus Tag(T *rc); -AmrTag CheckAllRefinement(MeshBlockData *rc, - const AmrTag &level = AmrTag::derefine); +AmrTag CheckAllRefinement(MeshBlockData *rc, const AmrTag &level); ParArray1D CheckAllRefinement(MeshData *mc); AmrTag FirstDerivative(const AMRBounds &bnds, const ParArray3D &q, diff --git a/src/mesh/mesh_refinement.cpp b/src/mesh/mesh_refinement.cpp index 75c9964e6a6b..23cd952794d7 100644 --- a/src/mesh/mesh_refinement.cpp +++ b/src/mesh/mesh_refinement.cpp @@ -68,17 +68,6 @@ MeshRefinement::MeshRefinement(std::weak_ptr pmb, ParameterInput *pin } } -//---------------------------------------------------------------------------------------- -//! \fn void MeshRefinement::CheckRefinementCondition() -// \brief Check refinement criteria - -void MeshRefinement::CheckRefinementCondition() { - std::shared_ptr pmb = GetBlockPointer(); - auto &rc = pmb->meshblock_data.Get(); - AmrTag ret = Refinement::CheckAllRefinement(rc.get()); - SetRefinement(ret); -} - void MeshRefinement::SetRefinement(AmrTag flag) { std::shared_ptr pmb = GetBlockPointer(); int aret = std::max(-1, static_cast(flag)); From 2b3f32a8751e91ac514187050bab73ad35a52f67 Mon Sep 17 00:00:00 2001 From: adam reyes Date: Sun, 29 Sep 2024 20:38:50 +0200 Subject: [PATCH 15/49] default CheckRefineMesh to true --- example/fine_advection/advection_package.cpp | 2 +- src/amr_criteria/refinement_package.cpp | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/example/fine_advection/advection_package.cpp b/example/fine_advection/advection_package.cpp index 94a40f6cb130..84d86dec8f67 100644 --- a/example/fine_advection/advection_package.cpp +++ b/example/fine_advection/advection_package.cpp @@ -95,7 +95,7 @@ std::shared_ptr Initialize(ParameterInput *pin) { Metadata({Metadata::Cell, Metadata::Derived, Metadata::OneCopy})); bool check_refine_mesh = - pin->GetOrAddBoolean("parthenon/mesh", "CheckRefineMesh", false); + pin->GetOrAddBoolean("parthenon/mesh", "CheckRefineMesh", true); if (check_refine_mesh) { pkg->CheckRefinementMesh = CheckRefinementMesh; } else { diff --git a/src/amr_criteria/refinement_package.cpp b/src/amr_criteria/refinement_package.cpp index 78d15c689e02..5d947015f646 100644 --- a/src/amr_criteria/refinement_package.cpp +++ b/src/amr_criteria/refinement_package.cpp @@ -37,7 +37,7 @@ namespace Refinement { std::shared_ptr Initialize(ParameterInput *pin) { auto ref = std::make_shared("Refinement"); bool check_refine_mesh = - pin->GetOrAddBoolean("parthenon/mesh", "CheckRefineMesh", false); + pin->GetOrAddBoolean("parthenon/mesh", "CheckRefineMesh", true); ref->AddParam("check_refine_mesh", check_refine_mesh); int numcrit = 0; @@ -296,7 +296,7 @@ TaskStatus Tag(MeshBlockData *rc) { template <> TaskStatus Tag(MeshData *rc) { PARTHENON_INSTRUMENT - ParArray1D delta_levels = CheckAllRefinement(rc, AmrTag::derefine); + ParArray1D delta_levels = CheckAllRefinement(rc); auto delta_levels_h = delta_levels.GetHostMirrorAndCopy(); for (int i = 0; i < rc->NumBlocks(); i++) { From 622882d9d9a24debcf348f9683a93e742ff1b662 Mon Sep 17 00:00:00 2001 From: adam reyes Date: Sun, 29 Sep 2024 21:27:40 +0200 Subject: [PATCH 16/49] respect amr_criteria max_level --- src/amr_criteria/amr_criteria.cpp | 4 ++-- src/amr_criteria/refinement_package.cpp | 14 ++++++++++---- src/amr_criteria/refinement_package.hpp | 6 ++++-- 3 files changed, 16 insertions(+), 8 deletions(-) diff --git a/src/amr_criteria/amr_criteria.cpp b/src/amr_criteria/amr_criteria.cpp index d867ff39604e..aecb766ec5d7 100644 --- a/src/amr_criteria/amr_criteria.cpp +++ b/src/amr_criteria/amr_criteria.cpp @@ -112,7 +112,7 @@ void AMRFirstDerivative::operator()(MeshData *mc, } const int idx = comp4 + n4 * (comp5 + n5 * comp6); Refinement::FirstDerivative(bnds, mc, field, idx, delta_levels, refine_criteria, - derefine_criteria); + derefine_criteria, max_level); } AmrTag AMRSecondDerivative::operator()(const MeshBlockData *rc) const { @@ -142,7 +142,7 @@ void AMRSecondDerivative::operator()(MeshData *mc, } const int idx = comp4 + n4 * (comp5 + n5 * comp6); Refinement::SecondDerivative(bnds, mc, field, idx, delta_levels, refine_criteria, - derefine_criteria); + derefine_criteria, max_level); } } // namespace parthenon diff --git a/src/amr_criteria/refinement_package.cpp b/src/amr_criteria/refinement_package.cpp index 5d947015f646..14f839256d6a 100644 --- a/src/amr_criteria/refinement_package.cpp +++ b/src/amr_criteria/refinement_package.cpp @@ -187,7 +187,8 @@ AmrTag SecondDerivative(const AMRBounds &bnds, const ParArray3D &q, void FirstDerivative(const AMRBounds &bnds, MeshData *mc, const std::string &field, const int &idx, ParArray1D &delta_levels, - const Real refine_criteria_, const Real derefine_criteria_) { + const Real refine_criteria_, const Real derefine_criteria_, + const int max_level_) { const auto desc = MakePackDescriptor(mc->GetMeshPointer()->resolved_packages.get(), {field}); auto pack = desc.GetPack(mc); @@ -196,6 +197,7 @@ void FirstDerivative(const AMRBounds &bnds, MeshData *mc, const std::strin const Real refine_criteria = refine_criteria_; const Real derefine_criteria = derefine_criteria_; + const int max_level = max_level_; const int var = idx; auto scatter_levels = delta_levels.ToScatterView(); par_for_outer( @@ -227,7 +229,8 @@ void FirstDerivative(const AMRBounds &bnds, MeshData *mc, const std::strin Kokkos::Max(maxd)); auto levels_access = scatter_levels.access(); auto flag = AmrTag::same; - if (maxd > refine_criteria) flag = AmrTag::refine; + if (maxd > refine_criteria && pack.GetLevel(b, 0, 0, 0) < max_level) + flag = AmrTag::refine; if (maxd < derefine_criteria) flag = AmrTag::derefine; levels_access(b).update(flag); }); @@ -236,7 +239,8 @@ void FirstDerivative(const AMRBounds &bnds, MeshData *mc, const std::strin void SecondDerivative(const AMRBounds &bnds, MeshData *mc, const std::string &field, const int &idx, ParArray1D &delta_levels, - const Real refine_criteria_, const Real derefine_criteria_) { + const Real refine_criteria_, const Real derefine_criteria_, + const int max_level_) { const auto desc = MakePackDescriptor(mc->GetMeshPointer()->resolved_packages.get(), {field}); auto pack = desc.GetPack(mc); @@ -245,6 +249,7 @@ void SecondDerivative(const AMRBounds &bnds, MeshData *mc, const std::stri const Real refine_criteria = refine_criteria_; const Real derefine_criteria = derefine_criteria_; + const int max_level = max_level_; const int var = idx; auto scatter_levels = delta_levels.ToScatterView(); par_for_outer( @@ -273,7 +278,8 @@ void SecondDerivative(const AMRBounds &bnds, MeshData *mc, const std::stri Kokkos::Max(maxd)); auto levels_access = scatter_levels.access(); auto flag = AmrTag::same; - if (maxd > refine_criteria) flag = AmrTag::refine; + if (maxd > refine_criteria && pack.GetLevel(b, 0, 0, 0) < max_level) + flag = AmrTag::refine; if (maxd < derefine_criteria) flag = AmrTag::derefine; levels_access(b).update(flag); }); diff --git a/src/amr_criteria/refinement_package.hpp b/src/amr_criteria/refinement_package.hpp index 444e1c5e3e32..3b1f92d9a96c 100644 --- a/src/amr_criteria/refinement_package.hpp +++ b/src/amr_criteria/refinement_package.hpp @@ -45,14 +45,16 @@ AmrTag FirstDerivative(const AMRBounds &bnds, const ParArray3D &q, void FirstDerivative(const AMRBounds &bnds, MeshData *mc, const std::string &field, const int &idx, ParArray1D &delta_levels_, - const Real refine_criteria_, const Real derefine_criteria_); + const Real refine_criteria_, const Real derefine_criteria_, + const int max_level_); AmrTag SecondDerivative(const AMRBounds &bnds, const ParArray3D &q, const Real refine_criteria, const Real derefine_criteria); void SecondDerivative(const AMRBounds &bnds, MeshData *mc, const std::string &field, const int &idx, ParArray1D &delta_levels_, - const Real refine_criteria_, const Real derefine_criteria_); + const Real refine_criteria_, const Real derefine_criteria_, + const int max_level_); } // namespace Refinement From 5b7863b3e5a921481b03f947e2de98f0c1acf48b Mon Sep 17 00:00:00 2001 From: adam reyes Date: Tue, 1 Oct 2024 01:27:04 +0200 Subject: [PATCH 17/49] renaming delta_levels->amr_tags, mc->md --- example/fine_advection/advection_package.cpp | 11 ++-- example/fine_advection/advection_package.hpp | 2 +- src/amr_criteria/amr_criteria.cpp | 28 ++++----- src/amr_criteria/amr_criteria.hpp | 6 +- src/amr_criteria/refinement_package.cpp | 64 ++++++++++---------- src/amr_criteria/refinement_package.hpp | 10 +-- 6 files changed, 60 insertions(+), 61 deletions(-) diff --git a/example/fine_advection/advection_package.cpp b/example/fine_advection/advection_package.cpp index 84d86dec8f67..51c86a82bcab 100644 --- a/example/fine_advection/advection_package.cpp +++ b/example/fine_advection/advection_package.cpp @@ -106,8 +106,7 @@ std::shared_ptr Initialize(ParameterInput *pin) { return pkg; } -void CheckRefinementMesh(MeshData *md, - parthenon::ParArray1D &delta_levels) { +void CheckRefinementMesh(MeshData *md, parthenon::ParArray1D &amr_tags) { // refine on advected, for example. could also be a derived quantity static auto desc = parthenon::MakePackDescriptor(md); auto pack = desc.GetPack(md); @@ -119,7 +118,7 @@ void CheckRefinementMesh(MeshData *md, auto ib = md->GetBoundsI(IndexDomain::entire); auto jb = md->GetBoundsJ(IndexDomain::entire); auto kb = md->GetBoundsK(IndexDomain::entire); - auto scatter_levels = delta_levels.ToScatterView(); + auto scatter_tags = amr_tags.ToScatterView(); parthenon::par_for_outer( PARTHENON_AUTO_LABEL, 0, 0, 0, pack.GetNBlocks() - 1, 0, pack.GetMaxNumberOfVars() - 1, kb.s, kb.e, @@ -137,14 +136,14 @@ void CheckRefinementMesh(MeshData *md, }, Kokkos::MinMax(minmax)); - auto levels_access = scatter_levels.access(); + auto tags_access = scatter_tags.access(); auto flag = AmrTag::same; if (minmax.max_val > refine_tol && minmax.min_val < derefine_tol) flag = AmrTag::refine; if (minmax.max_val < derefine_tol) flag = AmrTag::derefine; - levels_access(b).update(flag); + tags_access(b).update(flag); }); - delta_levels.ContributeScatter(scatter_levels); + amr_tags.ContributeScatter(scatter_tags); } AmrTag CheckRefinement(MeshBlockData *rc) { diff --git a/example/fine_advection/advection_package.hpp b/example/fine_advection/advection_package.hpp index d444cf29382a..6344d8320119 100644 --- a/example/fine_advection/advection_package.hpp +++ b/example/fine_advection/advection_package.hpp @@ -48,7 +48,7 @@ VARIABLE(advection, divD); std::shared_ptr Initialize(ParameterInput *pin); AmrTag CheckRefinement(MeshBlockData *rc); -void CheckRefinementMesh(MeshData *md, parthenon::ParArray1D &delta_levels); +void CheckRefinementMesh(MeshData *md, parthenon::ParArray1D &amr_tags); Real EstimateTimestep(MeshData *md); TaskStatus FillDerived(MeshData *md); diff --git a/src/amr_criteria/amr_criteria.cpp b/src/amr_criteria/amr_criteria.cpp index aecb766ec5d7..76f00d7a8905 100644 --- a/src/amr_criteria/amr_criteria.cpp +++ b/src/amr_criteria/amr_criteria.cpp @@ -95,13 +95,13 @@ AmrTag AMRFirstDerivative::operator()(const MeshBlockData *rc) const { return Refinement::FirstDerivative(bnds, q, refine_criteria, derefine_criteria); } -void AMRFirstDerivative::operator()(MeshData *mc, - ParArray1D &delta_levels) const { - auto ib = mc->GetBoundsI(IndexDomain::interior); - auto jb = mc->GetBoundsJ(IndexDomain::interior); - auto kb = mc->GetBoundsK(IndexDomain::interior); +void AMRFirstDerivative::operator()(MeshData *md, + ParArray1D &amr_tags) const { + auto ib = md->GetBoundsI(IndexDomain::interior); + auto jb = md->GetBoundsJ(IndexDomain::interior); + auto kb = md->GetBoundsK(IndexDomain::interior); auto bnds = AMRBounds(ib, jb, kb); - auto dims = mc->GetMeshPointer()->resolved_packages->FieldMetadata(field).Shape(); + auto dims = md->GetMeshPointer()->resolved_packages->FieldMetadata(field).Shape(); int n5(0), n4(0); if (dims.size() > 2) { n5 = dims[1]; @@ -111,7 +111,7 @@ void AMRFirstDerivative::operator()(MeshData *mc, n4 = dims[1]; } const int idx = comp4 + n4 * (comp5 + n5 * comp6); - Refinement::FirstDerivative(bnds, mc, field, idx, delta_levels, refine_criteria, + Refinement::FirstDerivative(bnds, md, field, idx, amr_tags, refine_criteria, derefine_criteria, max_level); } @@ -125,13 +125,13 @@ AmrTag AMRSecondDerivative::operator()(const MeshBlockData *rc) const { return Refinement::SecondDerivative(bnds, q, refine_criteria, derefine_criteria); } -void AMRSecondDerivative::operator()(MeshData *mc, - ParArray1D &delta_levels) const { - auto ib = mc->GetBoundsI(IndexDomain::interior); - auto jb = mc->GetBoundsJ(IndexDomain::interior); - auto kb = mc->GetBoundsK(IndexDomain::interior); +void AMRSecondDerivative::operator()(MeshData *md, + ParArray1D &amr_tags) const { + auto ib = md->GetBoundsI(IndexDomain::interior); + auto jb = md->GetBoundsJ(IndexDomain::interior); + auto kb = md->GetBoundsK(IndexDomain::interior); auto bnds = AMRBounds(ib, jb, kb); - auto dims = mc->GetMeshPointer()->resolved_packages->FieldMetadata(field).Shape(); + auto dims = md->GetMeshPointer()->resolved_packages->FieldMetadata(field).Shape(); int n5(0), n4(0); if (dims.size() > 2) { n5 = dims[1]; @@ -141,7 +141,7 @@ void AMRSecondDerivative::operator()(MeshData *mc, n4 = dims[1]; } const int idx = comp4 + n4 * (comp5 + n5 * comp6); - Refinement::SecondDerivative(bnds, mc, field, idx, delta_levels, refine_criteria, + Refinement::SecondDerivative(bnds, md, field, idx, amr_tags, refine_criteria, derefine_criteria, max_level); } diff --git a/src/amr_criteria/amr_criteria.hpp b/src/amr_criteria/amr_criteria.hpp index a57ebb1ebb9d..380639aa3e5d 100644 --- a/src/amr_criteria/amr_criteria.hpp +++ b/src/amr_criteria/amr_criteria.hpp @@ -37,7 +37,7 @@ struct AMRCriteria { AMRCriteria(ParameterInput *pin, std::string &block_name); virtual ~AMRCriteria() {} virtual AmrTag operator()(const MeshBlockData *rc) const = 0; - virtual void operator()(MeshData *mc, ParArray1D &delta_level) const = 0; + virtual void operator()(MeshData *md, ParArray1D &delta_level) const = 0; std::string field; Real refine_criteria, derefine_criteria; int max_level; @@ -51,14 +51,14 @@ struct AMRFirstDerivative : public AMRCriteria { AMRFirstDerivative(ParameterInput *pin, std::string &block_name) : AMRCriteria(pin, block_name) {} AmrTag operator()(const MeshBlockData *rc) const override; - void operator()(MeshData *mc, ParArray1D &delta_level) const override; + void operator()(MeshData *md, ParArray1D &delta_level) const override; }; struct AMRSecondDerivative : public AMRCriteria { AMRSecondDerivative(ParameterInput *pin, std::string &block_name) : AMRCriteria(pin, block_name) {} AmrTag operator()(const MeshBlockData *rc) const override; - void operator()(MeshData *mc, ParArray1D &delta_level) const override; + void operator()(MeshData *md, ParArray1D &delta_level) const override; }; } // namespace parthenon diff --git a/src/amr_criteria/refinement_package.cpp b/src/amr_criteria/refinement_package.cpp index 14f839256d6a..335d05125295 100644 --- a/src/amr_criteria/refinement_package.cpp +++ b/src/amr_criteria/refinement_package.cpp @@ -54,29 +54,29 @@ std::shared_ptr Initialize(ParameterInput *pin) { return ref; } -ParArray1D CheckAllRefinement(MeshData *mc) { - const int nblocks = mc->NumBlocks(); +ParArray1D CheckAllRefinement(MeshData *md) { + const int nblocks = md->NumBlocks(); // maybe not great to allocate this all the time - auto delta_levels = ParArray1D(Kokkos::View( - Kokkos::view_alloc(Kokkos::WithoutInitializing, "delta_levels"), nblocks)); - Kokkos::deep_copy(delta_levels.KokkosView(), AmrTag::derefine); + auto amr_tags = ParArray1D(Kokkos::View( + Kokkos::view_alloc(Kokkos::WithoutInitializing, "amr_tags"), nblocks)); + Kokkos::deep_copy(amr_tags.KokkosView(), AmrTag::derefine); - Mesh *pm = mc->GetMeshPointer(); + Mesh *pm = md->GetMeshPointer(); static const bool check_refine_mesh = pm->packages.Get("Refinement")->Param("check_refine_mesh"); for (auto &pkg : pm->packages.AllPackages()) { auto &desc = pkg.second; - desc->CheckRefinement(mc, delta_levels); + desc->CheckRefinement(md, amr_tags); if (check_refine_mesh) { for (auto &amr : desc->amr_criteria) { - (*amr)(mc, delta_levels); + (*amr)(md, amr_tags); } } } - return delta_levels; + return amr_tags; } AmrTag CheckAllRefinement(MeshBlockData *rc, const AmrTag &level) { @@ -185,21 +185,21 @@ AmrTag SecondDerivative(const AMRBounds &bnds, const ParArray3D &q, return AmrTag::same; } -void FirstDerivative(const AMRBounds &bnds, MeshData *mc, const std::string &field, - const int &idx, ParArray1D &delta_levels, +void FirstDerivative(const AMRBounds &bnds, MeshData *md, const std::string &field, + const int &idx, ParArray1D &amr_tags, const Real refine_criteria_, const Real derefine_criteria_, const int max_level_) { const auto desc = - MakePackDescriptor(mc->GetMeshPointer()->resolved_packages.get(), {field}); - auto pack = desc.GetPack(mc); - const int ndim = mc->GetMeshPointer()->ndim; + MakePackDescriptor(md->GetMeshPointer()->resolved_packages.get(), {field}); + auto pack = desc.GetPack(md); + const int ndim = md->GetMeshPointer()->ndim; const int nvars = pack.GetMaxNumberOfVars(); const Real refine_criteria = refine_criteria_; const Real derefine_criteria = derefine_criteria_; const int max_level = max_level_; const int var = idx; - auto scatter_levels = delta_levels.ToScatterView(); + auto scatter_tags = amr_tags.ToScatterView(); par_for_outer( PARTHENON_AUTO_LABEL, 0, 0, 0, pack.GetNBlocks() - 1, bnds.ks, bnds.ke, bnds.js, bnds.je, @@ -227,31 +227,31 @@ void FirstDerivative(const AMRBounds &bnds, MeshData *mc, const std::strin } }, Kokkos::Max(maxd)); - auto levels_access = scatter_levels.access(); + auto tags_access = scatter_tags.access(); auto flag = AmrTag::same; if (maxd > refine_criteria && pack.GetLevel(b, 0, 0, 0) < max_level) flag = AmrTag::refine; if (maxd < derefine_criteria) flag = AmrTag::derefine; - levels_access(b).update(flag); + tags_access(b).update(flag); }); - delta_levels.ContributeScatter(scatter_levels); + amr_tags.ContributeScatter(scatter_tags); } -void SecondDerivative(const AMRBounds &bnds, MeshData *mc, const std::string &field, - const int &idx, ParArray1D &delta_levels, +void SecondDerivative(const AMRBounds &bnds, MeshData *md, const std::string &field, + const int &idx, ParArray1D &amr_tags, const Real refine_criteria_, const Real derefine_criteria_, const int max_level_) { const auto desc = - MakePackDescriptor(mc->GetMeshPointer()->resolved_packages.get(), {field}); - auto pack = desc.GetPack(mc); - const int ndim = mc->GetMeshPointer()->ndim; + MakePackDescriptor(md->GetMeshPointer()->resolved_packages.get(), {field}); + auto pack = desc.GetPack(md); + const int ndim = md->GetMeshPointer()->ndim; const int nvars = pack.GetMaxNumberOfVars(); const Real refine_criteria = refine_criteria_; const Real derefine_criteria = derefine_criteria_; const int max_level = max_level_; const int var = idx; - auto scatter_levels = delta_levels.ToScatterView(); + auto scatter_tags = amr_tags.ToScatterView(); par_for_outer( PARTHENON_AUTO_LABEL, 0, 0, 0, pack.GetNBlocks() - 1, bnds.ks, bnds.ke, bnds.js, bnds.je, @@ -276,14 +276,14 @@ void SecondDerivative(const AMRBounds &bnds, MeshData *mc, const std::stri } }, Kokkos::Max(maxd)); - auto levels_access = scatter_levels.access(); + auto tags_access = scatter_tags.access(); auto flag = AmrTag::same; if (maxd > refine_criteria && pack.GetLevel(b, 0, 0, 0) < max_level) flag = AmrTag::refine; if (maxd < derefine_criteria) flag = AmrTag::derefine; - levels_access(b).update(flag); + tags_access(b).update(flag); }); - delta_levels.ContributeScatter(scatter_levels); + amr_tags.ContributeScatter(scatter_tags); } void SetRefinement_(MeshBlockData *rc, @@ -300,13 +300,13 @@ TaskStatus Tag(MeshBlockData *rc) { } template <> -TaskStatus Tag(MeshData *rc) { +TaskStatus Tag(MeshData *md) { PARTHENON_INSTRUMENT - ParArray1D delta_levels = CheckAllRefinement(rc); - auto delta_levels_h = delta_levels.GetHostMirrorAndCopy(); + ParArray1D amr_tags = CheckAllRefinement(md); + auto amr_tags_h = amr_tags.GetHostMirrorAndCopy(); - for (int i = 0; i < rc->NumBlocks(); i++) { - SetRefinement_(rc->GetBlockData(i).get(), delta_levels_h(i)); + for (int i = 0; i < md->NumBlocks(); i++) { + SetRefinement_(md->GetBlockData(i).get(), amr_tags_h(i)); } return TaskStatus::complete; } diff --git a/src/amr_criteria/refinement_package.hpp b/src/amr_criteria/refinement_package.hpp index 3b1f92d9a96c..cc71566fd92e 100644 --- a/src/amr_criteria/refinement_package.hpp +++ b/src/amr_criteria/refinement_package.hpp @@ -38,21 +38,21 @@ template TaskStatus Tag(T *rc); AmrTag CheckAllRefinement(MeshBlockData *rc, const AmrTag &level); -ParArray1D CheckAllRefinement(MeshData *mc); +ParArray1D CheckAllRefinement(MeshData *md); AmrTag FirstDerivative(const AMRBounds &bnds, const ParArray3D &q, const Real refine_criteria, const Real derefine_criteria); -void FirstDerivative(const AMRBounds &bnds, MeshData *mc, const std::string &field, - const int &idx, ParArray1D &delta_levels_, +void FirstDerivative(const AMRBounds &bnds, MeshData *md, const std::string &field, + const int &idx, ParArray1D &amr_tags, const Real refine_criteria_, const Real derefine_criteria_, const int max_level_); AmrTag SecondDerivative(const AMRBounds &bnds, const ParArray3D &q, const Real refine_criteria, const Real derefine_criteria); -void SecondDerivative(const AMRBounds &bnds, MeshData *mc, const std::string &field, - const int &idx, ParArray1D &delta_levels_, +void SecondDerivative(const AMRBounds &bnds, MeshData *md, const std::string &field, + const int &idx, ParArray1D &amr_tags, const Real refine_criteria_, const Real derefine_criteria_, const int max_level_); From 06e0aded5e9bcde3a08bc833a6bc6a883fe2674c Mon Sep 17 00:00:00 2001 From: adam reyes Date: Tue, 1 Oct 2024 01:33:12 +0200 Subject: [PATCH 18/49] fix refinement/bc order --- benchmarks/burgers/burgers_driver.cpp | 19 ++++++------------- 1 file changed, 6 insertions(+), 13 deletions(-) diff --git a/benchmarks/burgers/burgers_driver.cpp b/benchmarks/burgers/burgers_driver.cpp index 210f0276c99f..abb2a64a03d1 100644 --- a/benchmarks/burgers/burgers_driver.cpp +++ b/benchmarks/burgers/burgers_driver.cpp @@ -120,24 +120,17 @@ TaskCollection BurgersDriver::MakeTaskCollection(BlockList_t &blocks, const int auto fill_deriv = tl.AddTask(update, FillDerived>, mc1.get()); + auto set_bc = tl.AddTask(update, parthenon::ApplyBoundaryConditionsMD, mc1); + // estimate next time step if (stage == integrator->nstages) { auto new_dt = tl.AddTask(update, EstimateTimestep>, mc1.get()); - auto tag_refine = - tl.AddTask(update, parthenon::Refinement::Tag>, mc1.get()); + if (pmesh->adaptive) { + auto tag_refine = + tl.AddTask(set_bc, parthenon::Refinement::Tag>, mc1.get()); + } } } - - TaskRegion &async_region2 = tc.AddRegion(blocks.size()); - assert(blocks.size() == async_region2.size()); - for (int i = 0; i < blocks.size(); i++) { - auto &pmb = blocks[i]; - auto &tl = async_region2[i]; - auto &sc1 = pmb->meshblock_data.Get(stage_name[stage]); - - // set physical boundaries - auto set_bc = tl.AddTask(none, parthenon::ApplyBoundaryConditions, sc1); - } return tc; } From 3cf589da5425ed3a517a114490b792852a5f9254 Mon Sep 17 00:00:00 2001 From: adam reyes Date: Tue, 1 Oct 2024 01:40:43 +0200 Subject: [PATCH 19/49] docs for CheckRefinementMesh --- doc/sphinx/src/interface/state.rst | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/doc/sphinx/src/interface/state.rst b/doc/sphinx/src/interface/state.rst index b87832fe78a3..38ff841f2d68 100644 --- a/doc/sphinx/src/interface/state.rst +++ b/doc/sphinx/src/interface/state.rst @@ -104,6 +104,10 @@ several useful features and functions. ``std::function`` member ``CheckRefinementBlock`` if set (defaults to ``nullptr`` and therefore a no-op) that allows an application to define an application-specific refinement/de-refinement tagging function. +- ``void CheckRefinement(MeshData* md)`` delegates to the + ``std::function`` member ``CheckRefinementMesh`` if set (defaults to + ``nullptr`` and therefore a no-op) that allows an application to define + an application-specific refinement/de-refinement tagging function. - ``void PreStepDiagnostics(SimTime const &simtime, MeshData *rc)`` deletgates to the ``std::function`` member ``PreStepDiagnosticsMesh`` if set (defaults to ``nullptr`` an therefore a no-op) to print diagnostics From 7ea604cdc600b57c7a5ad6053c2b4b4681056deb Mon Sep 17 00:00:00 2001 From: adam reyes Date: Tue, 1 Oct 2024 10:41:43 +0200 Subject: [PATCH 20/49] move amr_tags array to mesh --- src/amr_criteria/refinement_package.cpp | 6 ++---- src/mesh/mesh.cpp | 12 ++++++++++++ src/mesh/mesh.hpp | 6 ++++++ 3 files changed, 20 insertions(+), 4 deletions(-) diff --git a/src/amr_criteria/refinement_package.cpp b/src/amr_criteria/refinement_package.cpp index 335d05125295..cd1b1502e90e 100644 --- a/src/amr_criteria/refinement_package.cpp +++ b/src/amr_criteria/refinement_package.cpp @@ -56,12 +56,10 @@ std::shared_ptr Initialize(ParameterInput *pin) { ParArray1D CheckAllRefinement(MeshData *md) { const int nblocks = md->NumBlocks(); - // maybe not great to allocate this all the time - auto amr_tags = ParArray1D(Kokkos::View( - Kokkos::view_alloc(Kokkos::WithoutInitializing, "amr_tags"), nblocks)); + Mesh *pm = md->GetMeshPointer(); + auto amr_tags = pm->GetAmrTags(); Kokkos::deep_copy(amr_tags.KokkosView(), AmrTag::derefine); - Mesh *pm = md->GetMeshPointer(); static const bool check_refine_mesh = pm->packages.Get("Refinement")->Param("check_refine_mesh"); diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp index ea838580c2b1..c4e3096f897b 100644 --- a/src/mesh/mesh.cpp +++ b/src/mesh/mesh.cpp @@ -913,6 +913,18 @@ const IndexShape Mesh::GetLeafBlockCellBounds(CellLevel level) const { } } +ParArray1D &Mesh::GetAmrTags() { + const int nblocks = GetNumMeshBlocksThisRank(); + if (!amr_tags.KokkosView().is_allocated()) { + amr_tags.KokkosView() = Kokkos::View( + Kokkos::view_alloc(Kokkos::WithoutInitializing, "amr_tags"), nblocks); + } + if (amr_tags.KokkosView().size() != nblocks) { + Kokkos::realloc(amr_tags.KokkosView(), nblocks); + } + return amr_tags; +} + // 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 684a897aad56..7a33a4ad43d3 100644 --- a/src/mesh/mesh.hpp +++ b/src/mesh/mesh.hpp @@ -34,6 +34,8 @@ #include #include +#include "Kokkos_Core.hpp" +#include "basic_types.hpp" #include "bvals/boundary_conditions.hpp" #include "bvals/comms/tag_map.hpp" #include "config.hpp" @@ -109,6 +111,8 @@ class Mesh { } const IndexShape GetLeafBlockCellBounds(CellLevel level = CellLevel::same) const; + ParArray1D &GetAmrTags(); + const forest::Forest &Forest() const { return forest; } // data @@ -276,6 +280,8 @@ class Mesh { std::vector bnref, bnderef; std::vector brdisp, bddisp; // the last 4x should be std::size_t, but are limited to int by MPI + // Refinement tags used by MeshData checks + ParArray1D amr_tags; std::vector loclist; From dd4a652dd8eb6eb93b9b8c61a37881c762fb2960 Mon Sep 17 00:00:00 2001 From: adam reyes Date: Wed, 2 Oct 2024 20:51:49 +0200 Subject: [PATCH 21/49] adding comments for ScatterMax view --- src/amr_criteria/refinement_package.cpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/amr_criteria/refinement_package.cpp b/src/amr_criteria/refinement_package.cpp index cd1b1502e90e..27d4dad213bb 100644 --- a/src/amr_criteria/refinement_package.cpp +++ b/src/amr_criteria/refinement_package.cpp @@ -197,6 +197,7 @@ void FirstDerivative(const AMRBounds &bnds, MeshData *md, const std::strin const Real derefine_criteria = derefine_criteria_; const int max_level = max_level_; const int var = idx; + // get a scatterview for the tags that will use Kokkos::Max as the reduction operation auto scatter_tags = amr_tags.ToScatterView(); par_for_outer( PARTHENON_AUTO_LABEL, 0, 0, 0, pack.GetNBlocks() - 1, bnds.ks, bnds.ke, bnds.js, @@ -230,6 +231,8 @@ void FirstDerivative(const AMRBounds &bnds, MeshData *md, const std::strin if (maxd > refine_criteria && pack.GetLevel(b, 0, 0, 0) < max_level) flag = AmrTag::refine; if (maxd < derefine_criteria) flag = AmrTag::derefine; + // ScatterMax view will use an atomic_max to prevent race condition across k,j + // indices tags_access(b).update(flag); }); amr_tags.ContributeScatter(scatter_tags); @@ -249,6 +252,7 @@ void SecondDerivative(const AMRBounds &bnds, MeshData *md, const std::stri const Real derefine_criteria = derefine_criteria_; const int max_level = max_level_; const int var = idx; + // get a scatterview for the tags that will use Kokkos::Max as the reduction operation auto scatter_tags = amr_tags.ToScatterView(); par_for_outer( PARTHENON_AUTO_LABEL, 0, 0, 0, pack.GetNBlocks() - 1, bnds.ks, bnds.ke, bnds.js, @@ -279,6 +283,8 @@ void SecondDerivative(const AMRBounds &bnds, MeshData *md, const std::stri if (maxd > refine_criteria && pack.GetLevel(b, 0, 0, 0) < max_level) flag = AmrTag::refine; if (maxd < derefine_criteria) flag = AmrTag::derefine; + // ScatterMax view will use an atomic_max to prevent race condition across k,j + // indices tags_access(b).update(flag); }); amr_tags.ContributeScatter(scatter_tags); From fe95b47146d73e10347b233692d626bc8c122943 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Thu, 12 Sep 2024 19:15:06 -0600 Subject: [PATCH 22/49] it compiles at least --- .../boundary_exchange_driver.cpp | 16 +-- example/particles/main.cpp | 14 +-- src/CMakeLists.txt | 3 + src/application_input.cpp | 97 +++++++++++++++++++ src/application_input.hpp | 39 ++++++-- src/bvals/boundary_conditions.cpp | 4 +- src/bvals/boundary_conditions.hpp | 6 +- src/bvals/boundary_flag.cpp | 48 +-------- src/bvals/bvals.hpp | 1 - src/defs.hpp | 6 +- src/mesh/forest/forest.hpp | 13 +-- src/mesh/forest/tree.cpp | 73 +++----------- src/mesh/forest/tree.hpp | 12 ++- src/mesh/mesh.cpp | 50 +++++++--- src/mesh/mesh.hpp | 13 +++ src/outputs/parthenon_hdf5.cpp | 8 +- src/outputs/parthenon_hdf5.hpp | 4 + src/outputs/parthenon_hdf5_attributes.cpp | 11 +++ src/outputs/restart.hpp | 1 - src/outputs/restart_hdf5.cpp | 2 - tst/unit/test_swarm.cpp | 2 +- 21 files changed, 256 insertions(+), 167 deletions(-) create mode 100644 src/application_input.cpp diff --git a/example/boundary_exchange/boundary_exchange_driver.cpp b/example/boundary_exchange/boundary_exchange_driver.cpp index bdb43d7e66ca..640705ca77ef 100644 --- a/example/boundary_exchange/boundary_exchange_driver.cpp +++ b/example/boundary_exchange/boundary_exchange_driver.cpp @@ -76,23 +76,23 @@ int main(int argc, char *argv[]) { ar3_t{1.0, 1.0, 1.0}); // forest_def.AddFace(0, {n[0], n[1], n[3], n[2]}, ar3_t{0.0, 0.0, 0.0}, // ar3_t{1.0, 1.0, 1.0}); - forest_def.AddBC(edge_t({n[0], n[1]}), parthenon::BoundaryFlag::outflow); - forest_def.AddBC(edge_t({n[0], n[3]}), parthenon::BoundaryFlag::outflow); + forest_def.AddBC(edge_t({n[0], n[1]})); + forest_def.AddBC(edge_t({n[0], n[3]})); forest_def.AddFace(1, {n[1], n[4], n[2], n[5]}, ar3_t{2.0, 0.0, 0.0}, ar3_t{3.0, 1.0, 1.0}); - forest_def.AddBC(edge_t({n[1], n[4]}), parthenon::BoundaryFlag::outflow); - forest_def.AddBC(edge_t({n[4], n[5]}), parthenon::BoundaryFlag::outflow); + forest_def.AddBC(edge_t({n[1], n[4]})); + forest_def.AddBC(edge_t({n[4], n[5]})); forest_def.AddFace(3, {n[3], n[2], n[6], n[7]}, ar3_t{0.0, 2.0, 0.0}, ar3_t{1.0, 3.0, 1.0}); - forest_def.AddBC(edge_t({n[6], n[7]}), parthenon::BoundaryFlag::outflow); - forest_def.AddBC(edge_t({n[3], n[6]}), parthenon::BoundaryFlag::outflow); + forest_def.AddBC(edge_t({n[6], n[7]})); + forest_def.AddBC(edge_t({n[3], n[6]})); forest_def.AddFace(4, {n[2], n[5], n[7], n[8]}, ar3_t{2.0, 2.0, 0.0}, ar3_t{3.0, 3.0, 1.0}); - forest_def.AddBC(edge_t({n[5], n[8]}), parthenon::BoundaryFlag::outflow); - forest_def.AddBC(edge_t({n[7], n[8]}), parthenon::BoundaryFlag::outflow); + forest_def.AddBC(edge_t({n[5], n[8]})); + forest_def.AddBC(edge_t({n[7], n[8]})); forest_def.AddInitialRefinement(parthenon::LogicalLocation(0, 1, 0, 0, 0)); pman.ParthenonInitPackagesAndMesh(forest_def); diff --git a/example/particles/main.cpp b/example/particles/main.cpp index 45b84fcfe823..e733033608c3 100644 --- a/example/particles/main.cpp +++ b/example/particles/main.cpp @@ -49,12 +49,14 @@ int main(int argc, char *argv[]) { // Redefine parthenon defaults pman.app_input->ProcessPackages = particles_example::ProcessPackages; - pman.app_input->boundary_conditions[BoundaryFace::inner_x1] = - parthenon::BoundaryFunction::OutflowInnerX1; - pman.app_input->boundary_conditions[BoundaryFace::outer_x1] = - parthenon::BoundaryFunction::OutflowOuterX1; - pman.app_input->swarm_boundary_conditions[BoundaryFace::inner_x1] = SwarmUserInnerX1; - pman.app_input->swarm_boundary_conditions[BoundaryFace::outer_x1] = SwarmUserOuterX1; + pman.app_input->RegisterBoundaryCondition(BoundaryFace::inner_x1, + parthenon::BoundaryFunction::OutflowInnerX1); + pman.app_input->RegisterBoundaryCondition(BoundaryFace::outer_x1, + parthenon::BoundaryFunction::OutflowOuterX1); + pman.app_input->RegisterSwarmBoundaryCondition(BoundaryFace::inner_x1, + SwarmUserInnerX1); + pman.app_input->RegisterSwarmBoundaryCondition(BoundaryFace::outer_x1, + SwarmUserOuterX1); // Note that this example does not use a ProblemGenerator pman.ParthenonInitPackagesAndMesh(); diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 133279563126..37883ddcb7af 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -95,6 +95,9 @@ set(COORDINATE_TYPE UniformCartesian) # TODO: Make this an option when more are configure_file(config.hpp.in generated/config.hpp @ONLY) add_library(parthenon + application_input.cpp + application_input.hpp + bvals/comms/bvals_in_one.hpp bvals/comms/bvals_utils.hpp bvals/comms/build_boundary_buffers.cpp diff --git a/src/application_input.cpp b/src/application_input.cpp new file mode 100644 index 000000000000..87721ef36072 --- /dev/null +++ b/src/application_input.cpp @@ -0,0 +1,97 @@ +//======================================================================================== +// (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 +// for the U.S. Department of Energy/National Nuclear Security Administration. All rights +// in the program are reserved by Triad National Security, LLC, and the U.S. Department +// of Energy/National Nuclear Security Administration. The Government is granted for +// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide +// license in this material to reproduce, prepare derivative works, distribute copies to +// the public, perform publicly and display publicly, and to permit others to do so. +//======================================================================================== + +#include + +#include "application_input.hpp" +#include "basic_types.hpp" +#include "bvals/boundary_conditions.hpp" +#include "utils/error_checking.hpp" + +namespace parthenon { +ApplicationInput::ApplicationInput() { + using namespace BoundaryFunction; + const std::string OUTFLOW = "outflow"; + const std::string PERIODIC = "periodic"; + + // Periodic handled separately for mesh + RegisterBoundaryCondition(BoundaryFace::inner_x1, OUTFLOW, &OutflowInnerX1); + RegisterBoundaryCondition(BoundaryFace::outer_x1, OUTFLOW, &OutflowOuterX1); + RegisterBoundaryCondition(BoundaryFace::inner_x2, OUTFLOW, &OutflowInnerX2); + RegisterBoundaryCondition(BoundaryFace::outer_x2, OUTFLOW, &OutflowOuterX2); + RegisterBoundaryCondition(BoundaryFace::inner_x3, OUTFLOW, &OutflowInnerX3); + RegisterBoundaryCondition(BoundaryFace::outer_x3, OUTFLOW, &OutflowOuterX3); + + // Periodic is explicit function for swarms + RegisterSwarmBoundaryCondition(BoundaryFace::inner_x1, OUTFLOW, &SwarmOutflowInnerX1); + RegisterSwarmBoundaryCondition(BoundaryFace::outer_x1, OUTFLOW, &SwarmOutflowOuterX1); + RegisterSwarmBoundaryCondition(BoundaryFace::inner_x2, OUTFLOW, &SwarmOutflowInnerX2); + RegisterSwarmBoundaryCondition(BoundaryFace::outer_x2, OUTFLOW, &SwarmOutflowOuterX2); + RegisterSwarmBoundaryCondition(BoundaryFace::inner_x3, OUTFLOW, &SwarmOutflowInnerX3); + RegisterSwarmBoundaryCondition(BoundaryFace::outer_x3, OUTFLOW, &SwarmOutflowOuterX3); + RegisterSwarmBoundaryCondition(BoundaryFace::inner_x1, PERIODIC, &SwarmPeriodicInnerX1); + RegisterSwarmBoundaryCondition(BoundaryFace::outer_x1, PERIODIC, &SwarmPeriodicOuterX1); + RegisterSwarmBoundaryCondition(BoundaryFace::inner_x2, PERIODIC, &SwarmPeriodicInnerX2); + RegisterSwarmBoundaryCondition(BoundaryFace::outer_x2, PERIODIC, &SwarmPeriodicOuterX2); + RegisterSwarmBoundaryCondition(BoundaryFace::inner_x3, PERIODIC, &SwarmPeriodicInnerX3); + RegisterSwarmBoundaryCondition(BoundaryFace::outer_x3, PERIODIC, &SwarmPeriodicOuterX3); +} + +void ApplicationInput::RegisterBoundaryCondition(BoundaryFace face, + const std::string &name, + BValFunc condition) { + if (boundary_conditions_[face].count(name) > 0) { + PARTHENON_THROW("Boundary condition " + name + " at face " + std::to_string(face) + + "already registered."); + } + boundary_conditions_[face][name] = condition; +} +void ApplicationInput::RegisterSwarmBoundaryCondition(BoundaryFace face, + const std::string &name, + SBValFunc condition) { + if (swarm_boundary_conditions_[face].count(name) > 0) { + PARTHENON_THROW("Swarm boundary condition " + name + " at face " + + std::to_string(face) + "already registered."); + } + swarm_boundary_conditions_[face][name] = condition; +} + +BValFunc ApplicationInput::GetBoundaryCondition(BoundaryFace face, + const std::string &name) const { + if (boundary_conditions_[face].count(name) == 0) { + std::stringstream msg; + msg << "Boundary condition " << name << " at face " << face << "not registered!\n" + << "Available conditions for this face are:\n"; + for (const auto &[name, func] : boundary_conditions_[face]) { + msg << name << "\n"; + } + PARTHENON_THROW(msg); + } + return boundary_conditions_[face].at(name); +} +SBValFunc ApplicationInput::GetSwarmBoundaryCondition(BoundaryFace face, + const std::string &name) const { + if (swarm_boundary_conditions_[face].count(name) == 0) { + std::stringstream msg; + msg << "Swarm boundary condition " << name << " at face " << face + << "not registered!\n" + << "Available conditions for this face are:\n"; + for (const auto &[name, func] : swarm_boundary_conditions_[face]) { + msg << name << "\n"; + } + PARTHENON_THROW(msg); + } + return swarm_boundary_conditions_[face].at(name); +} + +} // namespace parthenon diff --git a/src/application_input.hpp b/src/application_input.hpp index 3a8c9bbbee06..7abd6990e5be 100644 --- a/src/application_input.hpp +++ b/src/application_input.hpp @@ -19,17 +19,25 @@ #include #include +#include "basic_types.hpp" #include "bvals/boundary_conditions.hpp" #include "defs.hpp" -#include "interface/state_descriptor.hpp" -#include "outputs/output_parameters.hpp" -#include "parameter_input.hpp" #include "parthenon_arrays.hpp" namespace parthenon { +class Mesh; +class MeshBlock; +class MeshBlockApplicationData; +class OutputParameters; +class Packages_t; +class ParameterInput; +template +class MeshData; -struct ApplicationInput { +class ApplicationInput { public: + ApplicationInput(); + // ParthenonManager functions std::function &)> ProcessPackages = nullptr; @@ -57,8 +65,6 @@ struct ApplicationInput { std::function UserWorkAfterLoop = nullptr; std::function UserWorkBeforeLoop = nullptr; - BValFunc boundary_conditions[BOUNDARY_NFACES] = {nullptr}; - SBValFunc swarm_boundary_conditions[BOUNDARY_NFACES] = {nullptr}; // MeshBlock functions std::function(MeshBlock *, ParameterInput *)> @@ -68,6 +74,27 @@ struct ApplicationInput { std::function PostInitialization = nullptr; std::function MeshBlockUserWorkBeforeOutput = nullptr; + + // Physical boundaries + void RegisterBoundaryCondition(BoundaryFace face, const std::string &name, + BValFunc condition); + void RegisterBoundaryCondition(BoundaryFace face, BValFunc condition) { + RegisterBoundaryCondition(face, "user", condition); + } + + void RegisterSwarmBoundaryCondition(BoundaryFace face, const std::string &name, + SBValFunc condition); + void RegisterSwarmBoundaryCondition(BoundaryFace face, SBValFunc condition) { + RegisterSwarmBoundaryCondition(face, "user", condition); + } + + // Getters + BValFunc GetBoundaryCondition(BoundaryFace face, const std::string &name) const; + SBValFunc GetSwarmBoundaryCondition(BoundaryFace face, const std::string &name) const; + + private: + Dictionary boundary_conditions_[BOUNDARY_NFACES]; + Dictionary swarm_boundary_conditions_[BOUNDARY_NFACES]; }; } // namespace parthenon diff --git a/src/bvals/boundary_conditions.cpp b/src/bvals/boundary_conditions.cpp index 8f6777f5e705..648f04ac7478 100644 --- a/src/bvals/boundary_conditions.cpp +++ b/src/bvals/boundary_conditions.cpp @@ -207,7 +207,7 @@ bool DoPhysicalBoundary_(const BoundaryFlag flag, const BoundaryFace face, return false; } // ndim always at least 1 - return true; // reflect, outflow, user, dims correct + return true; // user, dims correct } bool DoPhysicalSwarmBoundary_(const BoundaryFlag flag, const BoundaryFace face, @@ -215,7 +215,7 @@ bool DoPhysicalSwarmBoundary_(const BoundaryFlag flag, const BoundaryFace face, if (flag == BoundaryFlag::undef) return false; if (flag == BoundaryFlag::block) return false; - return true; // outflow, periodic, user, dims (particles always 3D) correct + return true; // periodic, user, dims (particles always 3D) correct } } // namespace boundary_cond_impl diff --git a/src/bvals/boundary_conditions.hpp b/src/bvals/boundary_conditions.hpp index 020c20bcb720..b14ad1955e7f 100644 --- a/src/bvals/boundary_conditions.hpp +++ b/src/bvals/boundary_conditions.hpp @@ -14,11 +14,13 @@ #ifndef BVALS_BOUNDARY_CONDITIONS_HPP_ #define BVALS_BOUNDARY_CONDITIONS_HPP_ +#include #include #include -#include +#include #include "basic_types.hpp" +#include "defs.hpp" namespace parthenon { @@ -32,6 +34,8 @@ class Swarm; // Physical boundary conditions using BValFunc = std::function> &, bool)>; using SBValFunc = std::function &)>; +using BValFuncArray_t = std::array, BOUNDARY_NFACES>; +using SBValFuncArray_t = std::array, BOUNDARY_NFACES>; TaskStatus ApplyBoundaryConditionsOnCoarseOrFine(std::shared_ptr> &rc, bool coarse); diff --git a/src/bvals/boundary_flag.cpp b/src/bvals/boundary_flag.cpp index 649ba7bd7c5b..37dfa229939a 100644 --- a/src/bvals/boundary_flag.cpp +++ b/src/bvals/boundary_flag.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. 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 @@ -35,54 +35,14 @@ namespace parthenon { // condition. Typically called in Mesh() ctor and in pgen/*.cpp files. BoundaryFlag GetBoundaryFlag(const std::string &input_string) { - if (input_string == "reflecting") { - return BoundaryFlag::reflect; - } else if (input_string == "outflow") { - return BoundaryFlag::outflow; - } else if (input_string == "periodic") { + if (input_string == "periodic") { return BoundaryFlag::periodic; } else if (input_string == "none") { return BoundaryFlag::undef; } else if (input_string == "block") { return BoundaryFlag::block; - } else if (input_string == "user") { - return BoundaryFlag::user; } else { - std::stringstream msg; - msg << "### FATAL ERROR in GetBoundaryFlag" << std::endl - << "Input string=" << input_string << "\n" - << "is an invalid boundary type" << std::endl; - PARTHENON_FAIL(msg); - } -} - -//---------------------------------------------------------------------------------------- -//! \fn GetBoundaryString(BoundaryFlag input_flag) -// \brief Parses enumerated type BoundaryFlag internal integer representation to return -// string describing the boundary condition. Typicall used to format descriptive errors -// or diagnostics. Inverse of GetBoundaryFlag(). - -std::string GetBoundaryString(BoundaryFlag input_flag) { - switch (input_flag) { - case BoundaryFlag::block: // -1 - return "block"; - case BoundaryFlag::undef: // 0 - return "none"; - case BoundaryFlag::reflect: - return "reflecting"; - case BoundaryFlag::outflow: - return "outflow"; - case BoundaryFlag::periodic: - return "periodic"; - case BoundaryFlag::user: - return "user"; - default: - std::stringstream msg; - msg << "### FATAL ERROR in GetBoundaryString" << std::endl - << "Input enum class BoundaryFlag=" << static_cast(input_flag) << "\n" - << "is an invalid boundary type" << std::endl; - PARTHENON_FAIL(msg); - break; + return BoundaryFlag::user; } } @@ -96,7 +56,7 @@ std::string GetBoundaryString(BoundaryFlag input_flag) { void CheckBoundaryFlag(BoundaryFlag block_flag, CoordinateDirection dir) { std::stringstream msg; msg << "### FATAL ERROR in CheckBoundaryFlag" << std::endl - << "Attempting to set invalid MeshBlock boundary= " << GetBoundaryString(block_flag) + << "Attempting to set invalid MeshBlock boundary" << "\nin x" << dir << " direction" << std::endl; switch (dir) { case CoordinateDirection::X1DIR: diff --git a/src/bvals/bvals.hpp b/src/bvals/bvals.hpp index 0c93ae858f5b..bc6ed0ebed56 100644 --- a/src/bvals/bvals.hpp +++ b/src/bvals/bvals.hpp @@ -50,7 +50,6 @@ struct RegionSize; // free functions to return boundary flag given input string, and vice versa BoundaryFlag GetBoundaryFlag(const std::string &input_string); -std::string GetBoundaryString(BoundaryFlag input_flag); // + confirming that the MeshBlock's boundaries are all valid selections void CheckBoundaryFlag(BoundaryFlag block_flag, CoordinateDirection dir); diff --git a/src/defs.hpp b/src/defs.hpp index d058b14abfc6..cd656d05e037 100644 --- a/src/defs.hpp +++ b/src/defs.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 @@ -19,6 +19,7 @@ //! \file defs.hpp // \brief contains Athena++ general purpose types, structures, enums, etc. +#include #include #include #include @@ -121,7 +122,7 @@ struct RegionSize { // tasks.hpp, ??? // identifiers for boundary conditions -enum class BoundaryFlag { block = -1, undef, reflect, outflow, periodic, user }; +enum class BoundaryFlag { block = -1, undef, periodic, user }; // identifiers for all 6 faces of a MeshBlock constexpr int BOUNDARY_NFACES = 6; @@ -134,6 +135,7 @@ enum BoundaryFace { inner_x3 = 4, outer_x3 = 5 }; +using BValNames_t = std::array; inline BoundaryFace GetInnerBoundaryFace(CoordinateDirection dir) { if (dir == X1DIR) { diff --git a/src/mesh/forest/forest.hpp b/src/mesh/forest/forest.hpp index cf5f2780cbc1..db9de2e0025b 100644 --- a/src/mesh/forest/forest.hpp +++ b/src/mesh/forest/forest.hpp @@ -62,7 +62,8 @@ class ForestDefinition { face_sizes.emplace_back(xmin, xmax, ar3_t{1.0, 1.0, 1.0}, ai3_t{1, 1, 1}); } - void AddBC(Edge edge, BoundaryFlag bf, std::optional periodic_connection = {}) { + void AddBC(Edge edge, BoundaryFlag bf = BoundaryFlag::user, + std::optional periodic_connection = {}) { if (bf == BoundaryFlag::periodic) PARTHENON_REQUIRE(periodic_connection, "Must specify another edge for periodic boundary conditions."); @@ -134,12 +135,12 @@ class Forest { return trees.at(loc.tree())->GetBlockBCs(loc); } - void EnrollBndryFncts( - ApplicationInput *app_in, - std::array, BOUNDARY_NFACES> UserBoundaryFunctions_in, - std::array, BOUNDARY_NFACES> UserSwarmBoundaryFunctions_in) { + void EnrollBndryFncts(ApplicationInput *app_in, const BValNames_t &names, + const BValNames_t &swarm_names, + const BValFuncArray_t &UserBoundaryFunctions_in, + const SBValFuncArray_t &UserSwarmBoundaryFunctions_in) { for (auto &[id, ptree] : trees) - ptree->EnrollBndryFncts(app_in, UserBoundaryFunctions_in, + ptree->EnrollBndryFncts(app_in, names, swarm_names, UserBoundaryFunctions_in, UserSwarmBoundaryFunctions_in); } diff --git a/src/mesh/forest/tree.cpp b/src/mesh/forest/tree.cpp index fe1803882401..d1d019ee5bd4 100644 --- a/src/mesh/forest/tree.cpp +++ b/src/mesh/forest/tree.cpp @@ -24,12 +24,14 @@ #include #include +#include "application_input.hpp" #include "basic_types.hpp" #include "defs.hpp" #include "mesh/forest/forest_topology.hpp" #include "mesh/forest/logical_coordinate_transformation.hpp" #include "mesh/forest/logical_location.hpp" #include "mesh/forest/tree.hpp" +#include "parameter_input.hpp" #include "utils/bit_hacks.hpp" #include "utils/indexer.hpp" @@ -383,70 +385,21 @@ std::int64_t Tree::GetOldGid(const LogicalLocation &loc) const { return -1; } -void Tree::EnrollBndryFncts( - ApplicationInput *app_in, - std::array, BOUNDARY_NFACES> UserBoundaryFunctions_in, - std::array, BOUNDARY_NFACES> UserSwarmBoundaryFunctions_in) { +void Tree::EnrollBndryFncts(ApplicationInput *app_in, const BValNames_t &names, + const BValNames_t &swarm_names, + const BValFuncArray_t &UserBoundaryFunctions_in, + const SBValFuncArray_t &UserSwarmBoundaryFunctions_in) { UserBoundaryFunctions = UserBoundaryFunctions_in; UserSwarmBoundaryFunctions = UserSwarmBoundaryFunctions_in; - static const BValFunc outflow[6] = { - BoundaryFunction::OutflowInnerX1, BoundaryFunction::OutflowOuterX1, - BoundaryFunction::OutflowInnerX2, BoundaryFunction::OutflowOuterX2, - BoundaryFunction::OutflowInnerX3, BoundaryFunction::OutflowOuterX3}; - static const BValFunc reflect[6] = { - BoundaryFunction::ReflectInnerX1, BoundaryFunction::ReflectOuterX1, - BoundaryFunction::ReflectInnerX2, BoundaryFunction::ReflectOuterX2, - BoundaryFunction::ReflectInnerX3, BoundaryFunction::ReflectOuterX3}; - static const SBValFunc soutflow[6] = { - BoundaryFunction::SwarmOutflowInnerX1, BoundaryFunction::SwarmOutflowOuterX1, - BoundaryFunction::SwarmOutflowInnerX2, BoundaryFunction::SwarmOutflowOuterX2, - BoundaryFunction::SwarmOutflowInnerX3, BoundaryFunction::SwarmOutflowOuterX3}; - static const SBValFunc speriodic[6] = { - BoundaryFunction::SwarmPeriodicInnerX1, BoundaryFunction::SwarmPeriodicOuterX1, - BoundaryFunction::SwarmPeriodicInnerX2, BoundaryFunction::SwarmPeriodicOuterX2, - BoundaryFunction::SwarmPeriodicInnerX3, BoundaryFunction::SwarmPeriodicOuterX3}; for (int f = 0; f < BOUNDARY_NFACES; f++) { - switch (boundary_conditions[f]) { - case BoundaryFlag::reflect: - MeshBndryFnctn[f] = reflect[f]; - break; - case BoundaryFlag::outflow: - MeshBndryFnctn[f] = outflow[f]; - SwarmBndryFnctn[f] = soutflow[f]; - break; - case BoundaryFlag::user: - if (app_in->boundary_conditions[f] != nullptr) { - MeshBndryFnctn[f] = app_in->boundary_conditions[f]; - } else { - std::stringstream msg; - msg << "A user boundary condition for face " << f - << " was requested. but no condition was enrolled." << std::endl; - PARTHENON_THROW(msg); - } - break; - default: // periodic/block BCs handled elsewhere. - break; - } - - switch (boundary_conditions[f]) { - case BoundaryFlag::outflow: - SwarmBndryFnctn[f] = soutflow[f]; - break; - case BoundaryFlag::periodic: - SwarmBndryFnctn[f] = speriodic[f]; - break; - case BoundaryFlag::reflect: - // Default "reflect" boundaries not available for swarms; catch later on if swarms - // are present - break; - case BoundaryFlag::user: - if (app_in->swarm_boundary_conditions[f] != nullptr) { - SwarmBndryFnctn[f] = app_in->swarm_boundary_conditions[f]; - } - break; - default: // Default BCs handled elsewhere - break; + auto flag = boundary_conditions[f]; + auto face = static_cast(f); + if (flag == BoundaryFlag::user) { + MeshBndryFnctn[f] = app_in->GetBoundaryCondition(face, names[f]); + SwarmBndryFnctn[f] = app_in->GetSwarmBoundaryCondition(face, swarm_names[f]); + } else if (flag == BoundaryFlag::periodic) { + SwarmBndryFnctn[f] = app_in->GetSwarmBoundaryCondition(face, swarm_names[f]); } } } diff --git a/src/mesh/forest/tree.hpp b/src/mesh/forest/tree.hpp index 7b53073ac92b..f423b0180f78 100644 --- a/src/mesh/forest/tree.hpp +++ b/src/mesh/forest/tree.hpp @@ -23,7 +23,6 @@ #include #include -#include "application_input.hpp" #include "basic_types.hpp" #include "bvals/boundary_conditions.hpp" #include "defs.hpp" @@ -35,6 +34,9 @@ #include "utils/indexer.hpp" namespace parthenon { +class ApplicationInput; +class ParameterInput; + namespace forest { // We don't explicitly allow for periodic boundaries, since we can encode periodicity @@ -117,10 +119,10 @@ class Tree : public std::enable_shared_from_this { std::vector> forest_nodes; // Boundary Functions - void EnrollBndryFncts( - ApplicationInput *app_in, - std::array, BOUNDARY_NFACES> UserBoundaryFunctions_in, - std::array, BOUNDARY_NFACES> UserSwarmBoundaryFunctions_in); + void EnrollBndryFncts(ApplicationInput *app_in, const BValNames_t &names, + const BValNames_t &swarm_names, + const BValFuncArray_t &UserBoundaryFunctions_in, + const SBValFuncArray_t &UserSwarmBoundaryFunctions_in); BValFunc MeshBndryFnctn[BOUNDARY_NFACES]; SBValFunc SwarmBndryFnctn[BOUNDARY_NFACES]; std::array, BOUNDARY_NFACES> UserBoundaryFunctions; diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp index c4e3096f897b..3068c2c2c4ed 100644 --- a/src/mesh/mesh.cpp +++ b/src/mesh/mesh.cpp @@ -152,13 +152,8 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, pin->GetInteger("parthenon/mesh", "nx3")}, {false, pin->GetInteger("parthenon/mesh", "nx2") == 1, pin->GetInteger("parthenon/mesh", "nx3") == 1}); - mesh_bcs = { - GetBoundaryFlag(pin->GetOrAddString("parthenon/mesh", "ix1_bc", "reflecting")), - GetBoundaryFlag(pin->GetOrAddString("parthenon/mesh", "ox1_bc", "reflecting")), - GetBoundaryFlag(pin->GetOrAddString("parthenon/mesh", "ix2_bc", "reflecting")), - GetBoundaryFlag(pin->GetOrAddString("parthenon/mesh", "ox2_bc", "reflecting")), - GetBoundaryFlag(pin->GetOrAddString("parthenon/mesh", "ix3_bc", "reflecting")), - GetBoundaryFlag(pin->GetOrAddString("parthenon/mesh", "ox3_bc", "reflecting"))}; + SetBCNames_(pin); + mesh_bcs = GetBCsFromNames_(mesh_bc_names); ndim = (mesh_size.nx(X3DIR) > 1) ? 3 : ((mesh_size.nx(X2DIR) > 1) ? 2 : 1); for (auto &[dir, label] : std::vector>{ @@ -176,7 +171,8 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, // Load balancing flag and parameters forest = forest::Forest::HyperRectangular(mesh_size, base_block_size, mesh_bcs); root_level = forest.root_level; - forest.EnrollBndryFncts(app_in, resolved_packages->UserBoundaryFunctions, + forest.EnrollBndryFncts(app_in, mesh_bc_names, mesh_swarm_bc_names, + resolved_packages->UserBoundaryFunctions, resolved_packages->UserSwarmBoundaryFunctions); // SMR / AMR: @@ -196,13 +192,9 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, : Mesh(pin, app_in, packages, base_constructor_selector_t()) { mesh_size = RegionSize({0, 0, 0}, {1, 1, 0}, {1, 1, 1}, {1, 1, 1}, {false, false, true}); - mesh_bcs = { - GetBoundaryFlag(pin->GetOrAddString("parthenon/mesh", "ix1_bc", "reflecting")), - GetBoundaryFlag(pin->GetOrAddString("parthenon/mesh", "ox1_bc", "reflecting")), - GetBoundaryFlag(pin->GetOrAddString("parthenon/mesh", "ix2_bc", "reflecting")), - GetBoundaryFlag(pin->GetOrAddString("parthenon/mesh", "ox2_bc", "reflecting")), - GetBoundaryFlag(pin->GetOrAddString("parthenon/mesh", "ix3_bc", "reflecting")), - GetBoundaryFlag(pin->GetOrAddString("parthenon/mesh", "ox3_bc", "reflecting"))}; + + SetBCNames_(pin); + mesh_bcs = GetBCsFromNames_(mesh_bc_names); for (auto &[dir, label] : std::vector>{ {X1DIR, "nx1"}, {X2DIR, "nx2"}, {X3DIR, "nx3"}}) { base_block_size.xrat(dir) = mesh_size.xrat(dir); @@ -220,7 +212,8 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, // Load balancing flag and parameters forest = forest::Forest::Make2D(forest_def); root_level = forest.root_level; - forest.EnrollBndryFncts(app_in, resolved_packages->UserBoundaryFunctions, + forest.EnrollBndryFncts(app_in, mesh_bc_names, mesh_swarm_bc_names, + resolved_packages->UserBoundaryFunctions, resolved_packages->UserSwarmBoundaryFunctions); BuildBlockList(pin, app_in, packages, 0); } @@ -1196,4 +1189,29 @@ Mesh::GetLevelsAndLogicalLocationsFlat() const noexcept { return std::make_pair(levels, logicalLocations); } +void Mesh::SetBCNames_(ParameterInput *pin) { + mesh_bc_names = {pin->GetOrAddString("parthenon/mesh", "ix1_bc", "outflow"), + pin->GetOrAddString("parthenon/mesh", "ox1_bc", "outflow"), + pin->GetOrAddString("parthenon/mesh", "ix2_bc", "outflow"), + pin->GetOrAddString("parthenon/mesh", "ox2_bc", "outflow"), + pin->GetOrAddString("parthenon/mesh", "ix3_bc", "outflow"), + pin->GetOrAddString("parthenon/mesh", "ox3_bc", "outflow")}; + mesh_swarm_bc_names = { + pin->GetOrAddString("parthenon/mesh", "ix1_bc", mesh_bc_names[0]), + pin->GetOrAddString("parthenon/mesh", "ox1_bc", mesh_bc_names[1]), + pin->GetOrAddString("parthenon/mesh", "ix2_bc", mesh_bc_names[2]), + pin->GetOrAddString("parthenon/mesh", "ox2_bc", mesh_bc_names[3]), + pin->GetOrAddString("parthenon/mesh", "ix3_bc", mesh_bc_names[4]), + pin->GetOrAddString("parthenon/mesh", "ox3_bc", mesh_bc_names[5])}; +} + +std::array +Mesh::GetBCsFromNames_(const std::array &names) const { + std::array out; + for (int f = 0; f < BOUNDARY_NFACES; ++f) { + out[f] = GetBoundaryFlag(names[f]); + } + return out; +} + } // namespace parthenon diff --git a/src/mesh/mesh.hpp b/src/mesh/mesh.hpp index 7a33a4ad43d3..072db069d2f7 100644 --- a/src/mesh/mesh.hpp +++ b/src/mesh/mesh.hpp @@ -120,6 +120,15 @@ class Mesh { bool is_restart; RegionSize mesh_size; RegionSize base_block_size; + + BValNames_t mesh_bc_names; + BValNames_t mesh_swarm_bc_names; + + // these are flags not boundary functions + // JMM: A consequence of having only one boundary flag array but + // multiple boundary function arrays is that swarms *must* be + // periodic if the mesh is periodic but otherwise mesh and swarm + // boundaries are decoupled. std::array mesh_bcs; int ndim; // number of dimensions const bool adaptive, multilevel, multigrid; @@ -303,6 +312,10 @@ class Mesh { std::unordered_map mpi_comm_map_; #endif + void SetBCNames_(ParameterInput *pin); + std::array + GetBCsFromNames_(const BValNames_t &names) const; + // functions void CheckMeshValidity() const; void BuildBlockList(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, diff --git a/src/outputs/parthenon_hdf5.cpp b/src/outputs/parthenon_hdf5.cpp index 7d65d366cf77..cd5009d490d9 100644 --- a/src/outputs/parthenon_hdf5.cpp +++ b/src/outputs/parthenon_hdf5.cpp @@ -190,12 +190,8 @@ void PHDF5Output::WriteOutputFileImpl(Mesh *pm, ParameterInput *pin, SimTime *tm info_group); // Boundary conditions - std::vector boundary_condition_str(BOUNDARY_NFACES); - for (size_t i = 0; i < boundary_condition_str.size(); i++) { - boundary_condition_str[i] = GetBoundaryString(pm->mesh_bcs[i]); - } - - HDF5WriteAttribute("BoundaryConditions", boundary_condition_str, info_group); + HDF5WriteAttribute("BoundaryConditions", pm->mesh_bc_names, info_group); + HDF5WriteAttribute("SwarmBoundaryConditions", pm->mesh_swarm_bc_names, info_group); Kokkos::Profiling::popRegion(); // write Info } // Info section diff --git a/src/outputs/parthenon_hdf5.hpp b/src/outputs/parthenon_hdf5.hpp index 72deaab90681..e99c58b14815 100644 --- a/src/outputs/parthenon_hdf5.hpp +++ b/src/outputs/parthenon_hdf5.hpp @@ -18,6 +18,7 @@ #define OUTPUTS_PARTHENON_HDF5_HPP_ #include "config.hpp" +#include "defs.hpp" #include "kokkos_abstraction.hpp" #include "parthenon_arrays.hpp" @@ -122,6 +123,9 @@ void HDF5WriteAttribute(const std::string &name, size_t num_values, const T *dat // In CPP file void HDF5WriteAttribute(const std::string &name, const std::string &value, hid_t location); +void HDF5WriteAttribute(const std::string &name, + const std::array &values, + hid_t location); template void HDF5WriteAttribute(const std::string &name, const std::vector &values, diff --git a/src/outputs/parthenon_hdf5_attributes.cpp b/src/outputs/parthenon_hdf5_attributes.cpp index 45cece79881c..ef1d3c74c462 100644 --- a/src/outputs/parthenon_hdf5_attributes.cpp +++ b/src/outputs/parthenon_hdf5_attributes.cpp @@ -33,6 +33,7 @@ #include #include +#include "defs.hpp" #include "kokkos_abstraction.hpp" #include "utils/concepts_lite.hpp" #include "utils/error_checking.hpp" @@ -94,6 +95,16 @@ void HDF5WriteAttribute(const std::string &name, const std::vector HDF5WriteAttribute(name, char_ptrs, location); } +void HDF5WriteAttribute(const std::string &name, + const std::array &values, + hid_t location) { + std::vector char_ptrs(values.size()); + for (size_t i = 0; i < values.size(); ++i) { + char_ptrs[i] = values[i].c_str(); + } + HDF5WriteAttribute(name, char_ptrs, location); +} + template <> std::vector HDF5ReadAttributeVec(hid_t location, const std::string &name) { // get strings as char pointers, HDF5 will allocate the memory and we need to free it diff --git a/src/outputs/restart.hpp b/src/outputs/restart.hpp index 22b059fc3c5e..138949dd1667 100644 --- a/src/outputs/restart.hpp +++ b/src/outputs/restart.hpp @@ -86,7 +86,6 @@ class RestartReader { struct MeshInfo { int nbnew, nbdel, nbtotal, root_level, includes_ghost, n_ghost; - std::vector bound_cond; std::vector block_size; std::vector grid_dim; std::vector lx123; diff --git a/src/outputs/restart_hdf5.cpp b/src/outputs/restart_hdf5.cpp index c7584f954b8d..cababddbfafe 100644 --- a/src/outputs/restart_hdf5.cpp +++ b/src/outputs/restart_hdf5.cpp @@ -126,8 +126,6 @@ RestartReaderHDF5::MeshInfo RestartReaderHDF5::GetMeshInfo() const { mesh_info.nbtotal = GetAttr("Info", "NumMeshBlocks"); mesh_info.root_level = GetAttr("Info", "RootLevel"); - mesh_info.bound_cond = GetAttrVec("Info", "BoundaryConditions"); - mesh_info.block_size = GetAttrVec("Info", "MeshBlockSize"); mesh_info.includes_ghost = GetAttr("Info", "IncludesGhost"); mesh_info.n_ghost = GetAttr("Info", "NGhost"); diff --git a/tst/unit/test_swarm.cpp b/tst/unit/test_swarm.cpp index 22f6b1f4d735..ada71b112537 100644 --- a/tst/unit/test_swarm.cpp +++ b/tst/unit/test_swarm.cpp @@ -94,7 +94,7 @@ TEST_CASE("Swarm memory management", "[Swarm]") { mesh->mesh_bcs[0] = BoundaryFlag::user; meshblock->boundary_flag[0] = BoundaryFlag::user; for (int i = 1; i < 6; i++) { - mesh->mesh_bcs[i] = BoundaryFlag::outflow; + mesh->mesh_bcs[i] = BoundaryFlag::user; meshblock->boundary_flag[i] = BoundaryFlag::user; } meshblock->pmy_mesh = mesh.get(); From 554a1d42d90c54ad35f38357df841c16d55e3e74 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Thu, 12 Sep 2024 19:29:42 -0600 Subject: [PATCH 23/49] add easy machinery to register reflecting BCs --- doc/sphinx/src/boundary_conditions.rst | 31 +++++++++++++++++--------- src/application_input.cpp | 10 +++++++++ src/application_input.hpp | 1 + 3 files changed, 32 insertions(+), 10 deletions(-) diff --git a/doc/sphinx/src/boundary_conditions.rst b/doc/sphinx/src/boundary_conditions.rst index cb0e08532492..2f23bfe1f744 100644 --- a/doc/sphinx/src/boundary_conditions.rst +++ b/doc/sphinx/src/boundary_conditions.rst @@ -10,7 +10,6 @@ Natively, Parthenon supports three kinds of boundary conditions: - ``periodic`` - ``outflow`` -- ``reflecting`` which are all imposed on variables with the ``Metadata::FillGhost`` metadata flag. To set the boundaries in each direction, set the @@ -22,8 +21,8 @@ metadata flag. To set the boundaries in each direction, set the ix1_bc = outflow ox1_bc = outflow - ix2_bc = reflecting - ox2_bc = reflecting + ix2_bc = outflow + ox2_bc = outflow ix3_bc = periodic ox3_bc = periodic @@ -40,7 +39,9 @@ for your ``parthenon_manager``. e.g., .. code:: c++ - pman.app_input->boundary_conditions[parthenon::BoundaryFace::inner_x1] = MyBoundaryInnerX1; + pman.app_input->RegisterBoundaryCondition( + parthenon::BoundaryFace::inner_x1, + "my_bc_name", MyBoundaryInnerX1); where ``BoundaryFace`` is an enum defined in ``defs.hpp`` as @@ -58,13 +59,13 @@ where ``BoundaryFace`` is an enum defined in ``defs.hpp`` as outer_x3 = 5 }; -You can then set this boundary condition via the ``user`` flag in the -input file: +You can then set this boundary condition by using the name you +registered in the input file: :: - ix1_bc = user + ix1_bc = my_bc_name Boundary conditions so defined should look roughly like @@ -100,9 +101,19 @@ Other than these requirements, the ``Boundary`` object can do whatever you like. Reference implementations of the standard boundary conditions are available `here `__. +.. note:: -Per package user-defined boundary conditions --------------------------------------------- + A per-variable reflecting boundary condition is available, but you + must register it manually. To do so, simply call + ``app_in->RegisterDefaultReflectingBoundaryConditions()`` and it + will be available as a mesh boundary with the name ``reflecting``. + The reason manual registration is required is to support custom + reflecting boundary conditions int he case where a single variable + is used as the state vector. + + +Per package user-defined boundary conditions. +--------------------------------- In addition to user defined *global* boundary conditions, Parthenon also supports registration of boundary conditions at the *per package* level. These per package @@ -133,4 +144,4 @@ for a more complete example): pkg->UserBoundaryFunctions[BF::inner_x1].push_back(GetMyBC()); pkg->UserBoundaryFunctions[BF::inner_x2].push_back(GetMyBC()); ... - } \ No newline at end of file + } diff --git a/src/application_input.cpp b/src/application_input.cpp index 87721ef36072..42f7eb116272 100644 --- a/src/application_input.cpp +++ b/src/application_input.cpp @@ -47,6 +47,16 @@ ApplicationInput::ApplicationInput() { RegisterSwarmBoundaryCondition(BoundaryFace::outer_x3, PERIODIC, &SwarmPeriodicOuterX3); } +void ApplicationInput::RegisterDefaultReflectingBoundaryConditions() { + using namespace BoundaryFunction; + const std::string REFLECTING = "reflecting"; + RegisterBoundaryCondition(BoundaryFace::inner_x1, REFLECTING, &ReflectInnerX1); + RegisterBoundaryCondition(BoundaryFace::outer_x1, REFLECTING, &ReflectOuterX1); + RegisterBoundaryCondition(BoundaryFace::inner_x2, REFLECTING, &ReflectInnerX2); + RegisterBoundaryCondition(BoundaryFace::outer_x2, REFLECTING, &ReflectOuterX2); + RegisterBoundaryCondition(BoundaryFace::inner_x3, REFLECTING, &ReflectInnerX3); + RegisterBoundaryCondition(BoundaryFace::outer_x3, REFLECTING, &ReflectOuterX3); +} void ApplicationInput::RegisterBoundaryCondition(BoundaryFace face, const std::string &name, BValFunc condition) { diff --git a/src/application_input.hpp b/src/application_input.hpp index 7abd6990e5be..9c91585bfb3e 100644 --- a/src/application_input.hpp +++ b/src/application_input.hpp @@ -81,6 +81,7 @@ class ApplicationInput { void RegisterBoundaryCondition(BoundaryFace face, BValFunc condition) { RegisterBoundaryCondition(face, "user", condition); } + void RegisterDefaultReflectingBoundaryConditions(); void RegisterSwarmBoundaryCondition(BoundaryFace face, const std::string &name, SBValFunc condition); From 781d032f7eb4dc9309b96373d57b1c1b18411fc5 Mon Sep 17 00:00:00 2001 From: adam reyes Date: Fri, 8 Nov 2024 08:59:14 -0800 Subject: [PATCH 24/49] changelog --- CHANGELOG.md | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 8476cc25b96c..4be253c63f25 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,6 +13,7 @@ ### Changed (changing behavior/API/variables/...) - [[PR 1186]](https://github.com/parthenon-hpc-lab/parthenon/pull/1186) Bump Kokkos submodule to 4.4.1 +- [[PR1177]](https://github.com/parthenon-hpc-lab/parthenon/pull/1177) Make mesh-level boundary conditions usable without the "user" flag - [[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 @@ -27,7 +28,7 @@ ### Incompatibilities (i.e. breaking changes) - +- [[PR1177]](https://github.com/parthenon-hpc-lab/parthenon/pull/1177) Make mesh-level boundary conditions usable without the "user" flag ## Release 24.08 Date: 2024-08-30 From c1beb9ed60aef68b79accd269416a4c630a4ef1a Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Fri, 13 Sep 2024 08:33:34 -0600 Subject: [PATCH 25/49] swarm bcs differetn from mesh bcs --- src/mesh/mesh.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp index 3068c2c2c4ed..9784298c7efd 100644 --- a/src/mesh/mesh.cpp +++ b/src/mesh/mesh.cpp @@ -1197,12 +1197,12 @@ void Mesh::SetBCNames_(ParameterInput *pin) { pin->GetOrAddString("parthenon/mesh", "ix3_bc", "outflow"), pin->GetOrAddString("parthenon/mesh", "ox3_bc", "outflow")}; mesh_swarm_bc_names = { - pin->GetOrAddString("parthenon/mesh", "ix1_bc", mesh_bc_names[0]), - pin->GetOrAddString("parthenon/mesh", "ox1_bc", mesh_bc_names[1]), - pin->GetOrAddString("parthenon/mesh", "ix2_bc", mesh_bc_names[2]), - pin->GetOrAddString("parthenon/mesh", "ox2_bc", mesh_bc_names[3]), - pin->GetOrAddString("parthenon/mesh", "ix3_bc", mesh_bc_names[4]), - pin->GetOrAddString("parthenon/mesh", "ox3_bc", mesh_bc_names[5])}; + pin->GetOrAddString("parthenon/mesh", "ix1_swarm_bc", mesh_bc_names[0]), + pin->GetOrAddString("parthenon/mesh", "ox1_swarm_bc", mesh_bc_names[1]), + pin->GetOrAddString("parthenon/mesh", "ix2_swarm_bc", mesh_bc_names[2]), + pin->GetOrAddString("parthenon/mesh", "ox2_swarm_bc", mesh_bc_names[3]), + pin->GetOrAddString("parthenon/mesh", "ix3_swarm_bc", mesh_bc_names[4]), + pin->GetOrAddString("parthenon/mesh", "ox3_swarm_bc", mesh_bc_names[5])}; } std::array From 9e8b17fefbb0acec91759ddd48025f016e8a2e62 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Fri, 13 Sep 2024 08:33:58 -0600 Subject: [PATCH 26/49] use new input block --- src/mesh/mesh.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp index 9784298c7efd..7c0355933f43 100644 --- a/src/mesh/mesh.cpp +++ b/src/mesh/mesh.cpp @@ -1197,12 +1197,12 @@ void Mesh::SetBCNames_(ParameterInput *pin) { pin->GetOrAddString("parthenon/mesh", "ix3_bc", "outflow"), pin->GetOrAddString("parthenon/mesh", "ox3_bc", "outflow")}; mesh_swarm_bc_names = { - pin->GetOrAddString("parthenon/mesh", "ix1_swarm_bc", mesh_bc_names[0]), - pin->GetOrAddString("parthenon/mesh", "ox1_swarm_bc", mesh_bc_names[1]), - pin->GetOrAddString("parthenon/mesh", "ix2_swarm_bc", mesh_bc_names[2]), - pin->GetOrAddString("parthenon/mesh", "ox2_swarm_bc", mesh_bc_names[3]), - pin->GetOrAddString("parthenon/mesh", "ix3_swarm_bc", mesh_bc_names[4]), - pin->GetOrAddString("parthenon/mesh", "ox3_swarm_bc", mesh_bc_names[5])}; + pin->GetOrAddString("parthenon/swarm", "ix1_bc", mesh_bc_names[0]), + pin->GetOrAddString("parthenon/swarm", "ox1_bc", mesh_bc_names[1]), + pin->GetOrAddString("parthenon/swarm", "ix2_bc", mesh_bc_names[2]), + pin->GetOrAddString("parthenon/swarm", "ox2_bc", mesh_bc_names[3]), + pin->GetOrAddString("parthenon/swarm", "ix3_bc", mesh_bc_names[4]), + pin->GetOrAddString("parthenon/swarm", "ox3_bc", mesh_bc_names[5])}; } std::array From bcbc1c7ce3423c7e00726bdb292c340ce82983a7 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Fri, 27 Sep 2024 17:02:24 -0600 Subject: [PATCH 27/49] typo --- doc/sphinx/src/boundary_conditions.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/sphinx/src/boundary_conditions.rst b/doc/sphinx/src/boundary_conditions.rst index 2f23bfe1f744..9537f0c65728 100644 --- a/doc/sphinx/src/boundary_conditions.rst +++ b/doc/sphinx/src/boundary_conditions.rst @@ -1,4 +1,4 @@ -.. _sphinx-doc: +.. _boundary-conds: Boundary Conditions =================== @@ -108,7 +108,7 @@ are available `here RegisterDefaultReflectingBoundaryConditions()`` and it will be available as a mesh boundary with the name ``reflecting``. The reason manual registration is required is to support custom - reflecting boundary conditions int he case where a single variable + reflecting boundary conditions in the case where a single variable is used as the state vector. From b2a5e52eeb9a98d86b90530bdb57873546649f37 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Fri, 27 Sep 2024 17:18:37 -0600 Subject: [PATCH 28/49] add error checking for swarm/mesh BC consistency --- src/mesh/mesh.cpp | 10 ++++++++++ src/mesh/mesh.hpp | 4 ---- 2 files changed, 10 insertions(+), 4 deletions(-) diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp index 7c0355933f43..98123483e91e 100644 --- a/src/mesh/mesh.cpp +++ b/src/mesh/mesh.cpp @@ -1203,6 +1203,16 @@ void Mesh::SetBCNames_(ParameterInput *pin) { pin->GetOrAddString("parthenon/swarm", "ox2_bc", mesh_bc_names[3]), pin->GetOrAddString("parthenon/swarm", "ix3_bc", mesh_bc_names[4]), pin->GetOrAddString("parthenon/swarm", "ox3_bc", mesh_bc_names[5])}; + // JMM: A consequence of having only one boundary flag array but + // multiple boundary function arrays is that swarms *must* be + // periodic if the mesh is periodic but otherwise mesh and swarm + // boundaries are decoupled. + for (int i = 0; i < BOUNDARY_NFACES; ++i) { + if (mesh_bc_names[i] == "periodic") { + PARTHENON_REQUIRE(mesh_swarm_bc_names == "periodic", + "If the mesh is periodic, swarms must be also."); + } + } } std::array diff --git a/src/mesh/mesh.hpp b/src/mesh/mesh.hpp index 072db069d2f7..7e500eb3cf59 100644 --- a/src/mesh/mesh.hpp +++ b/src/mesh/mesh.hpp @@ -125,10 +125,6 @@ class Mesh { BValNames_t mesh_swarm_bc_names; // these are flags not boundary functions - // JMM: A consequence of having only one boundary flag array but - // multiple boundary function arrays is that swarms *must* be - // periodic if the mesh is periodic but otherwise mesh and swarm - // boundaries are decoupled. std::array mesh_bcs; int ndim; // number of dimensions const bool adaptive, multilevel, multigrid; From 46517ab169bc1f1642db9817074e04e0fce2c773 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Fri, 27 Sep 2024 17:19:11 -0600 Subject: [PATCH 29/49] typo --- src/mesh/mesh.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp index 98123483e91e..8f5a1d3ed908 100644 --- a/src/mesh/mesh.cpp +++ b/src/mesh/mesh.cpp @@ -1209,7 +1209,7 @@ void Mesh::SetBCNames_(ParameterInput *pin) { // boundaries are decoupled. for (int i = 0; i < BOUNDARY_NFACES; ++i) { if (mesh_bc_names[i] == "periodic") { - PARTHENON_REQUIRE(mesh_swarm_bc_names == "periodic", + PARTHENON_REQUIRE(mesh_swarm_bc_names[i] == "periodic", "If the mesh is periodic, swarms must be also."); } } From f9313d334aae363f41a7b25fd2e244e708e831cb Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Fri, 27 Sep 2024 18:01:56 -0600 Subject: [PATCH 30/49] phdf diff --- .../packages/parthenon_tools/parthenon_tools/phdf_diff.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/scripts/python/packages/parthenon_tools/parthenon_tools/phdf_diff.py b/scripts/python/packages/parthenon_tools/parthenon_tools/phdf_diff.py index b7cb8b1d74c4..dd1ff9051d81 100644 --- a/scripts/python/packages/parthenon_tools/parthenon_tools/phdf_diff.py +++ b/scripts/python/packages/parthenon_tools/parthenon_tools/phdf_diff.py @@ -211,7 +211,8 @@ def compare_metadata(f0, f1, quiet=False, one=False, check_input=False, tol=1.0e if one: return ret_code - # Compare the names of attributes in /Info, except "Time" + # Compare the names of attributes in /Info, except those we know + # may vary safely f0_Info = { key: value for key, value in f0.Info.items() @@ -219,6 +220,8 @@ def compare_metadata(f0, f1, quiet=False, one=False, check_input=False, tol=1.0e and key != "BlocksPerPE" and key != "WallTime" and key != "OutputFormatVersion" + and key != "BoundaryConditions" + and key != "SwarmBoundaryConditions" } f1_Info = { key: value @@ -227,6 +230,8 @@ def compare_metadata(f0, f1, quiet=False, one=False, check_input=False, tol=1.0e and key != "BlocksPerPE" and key != "WallTime" and key != "OutputFormatVersion" + and key != "BoundaryConditions" + and key != "SwarmBoundaryConditions" } if sorted(f0_Info.keys()) != sorted(f1_Info.keys()): print("Names of attributes in '/Info' of differ") From 20abf5a8f677a0e1a92593d9a6306c2fcaa9e895 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Sat, 28 Sep 2024 10:55:57 -0600 Subject: [PATCH 31/49] Register reflecting BCs for advection examples --- example/advection/main.cpp | 1 + example/sparse_advection/main.cpp | 1 + 2 files changed, 2 insertions(+) diff --git a/example/advection/main.cpp b/example/advection/main.cpp index 146093bc22eb..43764c35f774 100644 --- a/example/advection/main.cpp +++ b/example/advection/main.cpp @@ -25,6 +25,7 @@ int main(int argc, char *argv[]) { pman.app_input->ProblemGenerator = advection_example::ProblemGenerator; pman.app_input->UserWorkAfterLoop = advection_example::UserWorkAfterLoop; pman.app_input->UserMeshWorkBeforeOutput = advection_example::UserMeshWorkBeforeOutput; + mpan.app_input->RegisterDefaultReflectingBoundaryConditions(); // call ParthenonInit to initialize MPI and Kokkos, parse the input deck, and set up auto manager_status = pman.ParthenonInitEnv(argc, argv); diff --git a/example/sparse_advection/main.cpp b/example/sparse_advection/main.cpp index 585d1a35ab1e..4e19a6ac8034 100644 --- a/example/sparse_advection/main.cpp +++ b/example/sparse_advection/main.cpp @@ -27,6 +27,7 @@ int main(int argc, char *argv[]) { pman.app_input->ProblemGenerator = sparse_advection_example::ProblemGenerator; pman.app_input->PostStepDiagnosticsInLoop = sparse_advection_example::PostStepDiagnosticsInLoop; + mpan.app_input->RegisterDefaultReflectingBoundaryConditions(); // call ParthenonInit to initialize MPI and Kokkos, parse the input deck, and set up auto manager_status = pman.ParthenonInitEnv(argc, argv); From 1408bf08ff32039d84e78f83f0672d62b9fcb096 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Sat, 28 Sep 2024 13:03:18 -0600 Subject: [PATCH 32/49] typo --- example/advection/main.cpp | 2 +- example/sparse_advection/main.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/example/advection/main.cpp b/example/advection/main.cpp index 43764c35f774..5ed943b07474 100644 --- a/example/advection/main.cpp +++ b/example/advection/main.cpp @@ -25,7 +25,7 @@ int main(int argc, char *argv[]) { pman.app_input->ProblemGenerator = advection_example::ProblemGenerator; pman.app_input->UserWorkAfterLoop = advection_example::UserWorkAfterLoop; pman.app_input->UserMeshWorkBeforeOutput = advection_example::UserMeshWorkBeforeOutput; - mpan.app_input->RegisterDefaultReflectingBoundaryConditions(); + pman.app_input->RegisterDefaultReflectingBoundaryConditions(); // call ParthenonInit to initialize MPI and Kokkos, parse the input deck, and set up auto manager_status = pman.ParthenonInitEnv(argc, argv); diff --git a/example/sparse_advection/main.cpp b/example/sparse_advection/main.cpp index 4e19a6ac8034..ce6b197e335a 100644 --- a/example/sparse_advection/main.cpp +++ b/example/sparse_advection/main.cpp @@ -27,7 +27,7 @@ int main(int argc, char *argv[]) { pman.app_input->ProblemGenerator = sparse_advection_example::ProblemGenerator; pman.app_input->PostStepDiagnosticsInLoop = sparse_advection_example::PostStepDiagnosticsInLoop; - mpan.app_input->RegisterDefaultReflectingBoundaryConditions(); + pman.app_input->RegisterDefaultReflectingBoundaryConditions(); // call ParthenonInit to initialize MPI and Kokkos, parse the input deck, and set up auto manager_status = pman.ParthenonInitEnv(argc, argv); From 5b91a7ad648fa9867cb3e05001db6cd38f33acde Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Sat, 28 Sep 2024 15:03:01 -0600 Subject: [PATCH 33/49] silly backwards compatibility thing to make it so you don't have to specify swarm BCs if you're not using swarms --- src/mesh/mesh.cpp | 17 +++++++++++------ 1 file changed, 11 insertions(+), 6 deletions(-) diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp index 8f5a1d3ed908..27a0c778cb33 100644 --- a/src/mesh/mesh.cpp +++ b/src/mesh/mesh.cpp @@ -1196,13 +1196,18 @@ void Mesh::SetBCNames_(ParameterInput *pin) { pin->GetOrAddString("parthenon/mesh", "ox2_bc", "outflow"), pin->GetOrAddString("parthenon/mesh", "ix3_bc", "outflow"), pin->GetOrAddString("parthenon/mesh", "ox3_bc", "outflow")}; + // JMM: This is needed because not all BCs are necessarily + // implemented for swarms + auto maybe = [](const std::string &s) { + return ((s == "outflow") || (s == "periodic")) ? s : "outflow"; + }; mesh_swarm_bc_names = { - pin->GetOrAddString("parthenon/swarm", "ix1_bc", mesh_bc_names[0]), - pin->GetOrAddString("parthenon/swarm", "ox1_bc", mesh_bc_names[1]), - pin->GetOrAddString("parthenon/swarm", "ix2_bc", mesh_bc_names[2]), - pin->GetOrAddString("parthenon/swarm", "ox2_bc", mesh_bc_names[3]), - pin->GetOrAddString("parthenon/swarm", "ix3_bc", mesh_bc_names[4]), - pin->GetOrAddString("parthenon/swarm", "ox3_bc", mesh_bc_names[5])}; + pin->GetOrAddString("parthenon/swarm", "ix1_bc", maybe(mesh_bc_names[0])), + pin->GetOrAddString("parthenon/swarm", "ox1_bc", maybe(mesh_bc_names[1])), + pin->GetOrAddString("parthenon/swarm", "ix2_bc", maybe(mesh_bc_names[2])), + pin->GetOrAddString("parthenon/swarm", "ox2_bc", maybe(mesh_bc_names[3])), + pin->GetOrAddString("parthenon/swarm", "ix3_bc", maybe(mesh_bc_names[4])), + pin->GetOrAddString("parthenon/swarm", "ox3_bc", maybe(mesh_bc_names[5]))}; // JMM: A consequence of having only one boundary flag array but // multiple boundary function arrays is that swarms *must* be // periodic if the mesh is periodic but otherwise mesh and swarm From 382a6a7d977ba786bc0295041a2c1404e851a92f Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Wed, 9 Oct 2024 20:54:11 -0600 Subject: [PATCH 34/49] working --- src/interface/data_collection.hpp | 2 +- src/interface/mesh_data.hpp | 8 +++++ src/interface/meshblock_data.hpp | 49 ++++++++++++++++++++++++++++-- src/interface/state_descriptor.cpp | 1 + src/interface/state_descriptor.hpp | 8 +++++ src/mesh/mesh-amr_loadbalance.cpp | 5 +-- src/mesh/mesh.cpp | 11 ++++--- src/mesh/mesh.hpp | 3 +- 8 files changed, 76 insertions(+), 11 deletions(-) diff --git a/src/interface/data_collection.hpp b/src/interface/data_collection.hpp index 043bf878f8b8..485b642498f9 100644 --- a/src/interface/data_collection.hpp +++ b/src/interface/data_collection.hpp @@ -58,7 +58,7 @@ class DataCollection { auto key = GetKey(name, src); auto it = containers_.find(key); if (it != containers_.end()) { - if (fields.size() && !(it->second)->Contains(fields)) { + if (fields.size() && !(it->second)->CreatedFrom(fields)) { PARTHENON_THROW(key + " already exists in collection but fields do not match."); } return it->second; diff --git a/src/interface/mesh_data.hpp b/src/interface/mesh_data.hpp index 621ffce9763f..14ae3959c32a 100644 --- a/src/interface/mesh_data.hpp +++ b/src/interface/mesh_data.hpp @@ -476,6 +476,14 @@ class MeshData { [this, vars](const auto &b) { return b->ContainsExactly(vars); }); } + // Checks that the same set of variables was requested to create this container + // (which may be different than the set of variables in the container because of fluxes) + template + bool CreatedFrom(const Vars_t &vars) const noexcept { + return std::all_of(block_data_.begin(), block_data_.end(), + [this, vars](const auto &b) { return b->CreatedFrom(vars); }); + } + std::shared_ptr GetSwarmData(int n) { PARTHENON_REQUIRE(n >= 0 && n < block_data_.size(), "MeshData::GetSwarmData requires n within [0, block_data_.size()]"); diff --git a/src/interface/meshblock_data.hpp b/src/interface/meshblock_data.hpp index 1ead0e5f957d..2e95800c2e79 100644 --- a/src/interface/meshblock_data.hpp +++ b/src/interface/meshblock_data.hpp @@ -141,6 +141,17 @@ class MeshBlockData { resolved_packages = resolved_packages_in; is_shallow_ = shallow_copy; + // Store the list of variables used to create this container + // so we can compare to it when searching the cache + varUidIn_.clear(); + if constexpr (std::is_same_v) { + for (const auto &var : vars) + varUidIn_.insert(Variable::GetUniqueID(var)); + } else { + for (const auto &var : vars) + varUidIn_.insert(var); + } + // clear all variables, maps, and pack caches varVector_.clear(); varMap_.clear(); @@ -185,9 +196,29 @@ class MeshBlockData { if (!found) add_var(src->GetVarPtr(flx_name)); } } + } else if constexpr (std::is_same_v && + std::is_same_v) { + for (const auto &v : vars) { + const auto &vid = resolved_packages->GetFieldVarID(v); + const auto &md = resolved_packages->GetFieldMetadata(v); + AddField(vid.base_name, md, vid.sparse_id); + // Add the associated flux as well if not explicitly + // asked for + if (md.IsSet(Metadata::WithFluxes)) { + auto flx_name = md.GetFluxName(); + bool found = false; + for (const auto &v2 : vars) + if (v2 == flx_name) found = true; + if (!found) { + const auto &vid = resolved_packages->GetFieldVarID(flx_name); + const auto &md = resolved_packages->GetFieldMetadata(flx_name); + AddField(vid.base_name, md, vid.sparse_id); + } + } + } } else { - PARTHENON_FAIL( - "Variable subset selection not yet implemented for MeshBlock input."); + PARTHENON_FAIL("Variable subset selection not yet implemented for MeshBlock " + "input with unique ids."); } } @@ -525,6 +556,18 @@ class MeshBlockData { return Contains(vars) && (vars.size() == varVector_.size()); } + bool CreatedFrom(const std::vector &vars) { + return (vars.size() == varUidIn_.size()) && + std::all_of(vars.begin(), vars.end(), + [this](const auto &v) { return this->varUidIn_.count(v); }); + } + bool CreatedFrom(const std::vector &vars) { + return (vars.size() == varUidIn_.size()) && + std::all_of(vars.begin(), vars.end(), [this](const auto &v) { + return this->varUidIn_.count(Variable::GetUniqueID(v)); + }); + } + void SetAllVariablesToInitialized() { std::for_each(varVector_.begin(), varVector_.end(), [](auto &sp_var) { sp_var->data.initialized = true; }); @@ -561,6 +604,8 @@ class MeshBlockData { VariableVector varVector_; ///< the saved variable array std::map>> varUidMap_; + std::set varUidIn_; // Uid list from which this MeshBlockData was created, + // empty implies all variables were included MapToVars varMap_; MetadataFlagToVariableMap flagsToVars_; diff --git a/src/interface/state_descriptor.cpp b/src/interface/state_descriptor.cpp index 9116298e8bd5..d970887e50ae 100644 --- a/src/interface/state_descriptor.cpp +++ b/src/interface/state_descriptor.cpp @@ -316,6 +316,7 @@ bool StateDescriptor::AddFieldImpl_(const VarID &vid, const Metadata &m_in, AddFieldImpl_(fId, *(m.GetSPtrFluxMetadata()), control_vid); m.SetFluxName(fId.label()); } + labelToVidMap_.insert({vid.label(), vid}); metadataMap_.insert({vid, m}); refinementFuncMaps_.Register(m, vid.label()); allocControllerReverseMap_.insert({vid, control_vid}); diff --git a/src/interface/state_descriptor.hpp b/src/interface/state_descriptor.hpp index c3a034bb377a..ff143605b2a1 100644 --- a/src/interface/state_descriptor.hpp +++ b/src/interface/state_descriptor.hpp @@ -203,6 +203,13 @@ class StateDescriptor { // retrieve all swarm names std::vector Swarms() noexcept; + const auto &GetFieldVarID(const std::string &label) const { + return labelToVidMap_.at(label); + } + const auto &GetFieldMetadata(const std::string &label) const { + return metadataMap_.at(labelToVidMap_.at(label)); + } + const auto &GetFieldMetadata(const VarID &id) const { return metadataMap_.at(id); } const auto &AllFields() const noexcept { return metadataMap_; } const auto &AllSparsePools() const noexcept { return sparsePoolMap_; } const auto &AllSwarms() const noexcept { return swarmMetadataMap_; } @@ -403,6 +410,7 @@ class StateDescriptor { const std::string label_; // for each variable label (full label for sparse variables) hold metadata + std::unordered_map labelToVidMap_; std::unordered_map metadataMap_; std::unordered_map allocControllerReverseMap_; std::unordered_map> allocControllerMap_; diff --git a/src/mesh/mesh-amr_loadbalance.cpp b/src/mesh/mesh-amr_loadbalance.cpp index d54167026066..00d768720988 100644 --- a/src/mesh/mesh-amr_loadbalance.cpp +++ b/src/mesh/mesh-amr_loadbalance.cpp @@ -979,8 +979,9 @@ void Mesh::RedistributeAndRefineMeshBlocks(ParameterInput *pin, ApplicationInput auto &md_noncc = mesh_data.AddShallow(noncc, md, noncc_names); } - CommunicateBoundaries(noncc); // Called to make sure shared values are correct, - // ghosts of non-cell centered vars may get some junk + CommunicateBoundaries( + noncc, noncc_names); // Called to make sure shared values are correct, + // ghosts of non-cell centered vars may get some junk // Now there is the correct data for prolongating on un-shared topological elements // on the new fine blocks if (nprolong > 0) { diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp index 27a0c778cb33..b4d48ae5bccd 100644 --- a/src/mesh/mesh.cpp +++ b/src/mesh/mesh.cpp @@ -638,7 +638,8 @@ void Mesh::BuildTagMapAndBoundaryBuffers() { } } -void Mesh::CommunicateBoundaries(std::string md_name) { +void Mesh::CommunicateBoundaries(std::string md_name, + const std::vector &fields) { const int num_partitions = DefaultNumPartitions(); const int nmb = GetNumMeshBlocksThisRank(Globals::my_rank); constexpr std::int64_t max_it = 1e10; @@ -650,7 +651,7 @@ void Mesh::CommunicateBoundaries(std::string md_name) { do { all_sent = true; for (int i = 0; i < partitions.size(); ++i) { - auto &md = mesh_data.Add(md_name, partitions[i]); + auto &md = mesh_data.Add(md_name, partitions[i], fields); if (!sent[i]) { if (SendBoundaryBuffers(md) != TaskStatus::complete) { all_sent = false; @@ -674,7 +675,7 @@ void Mesh::CommunicateBoundaries(std::string md_name) { do { all_received = true; for (int i = 0; i < partitions.size(); ++i) { - auto &md = mesh_data.Add(md_name, partitions[i]); + auto &md = mesh_data.Add(md_name, partitions[i], fields); if (!received[i]) { if (ReceiveBoundaryBuffers(md) != TaskStatus::complete) { all_received = false; @@ -690,14 +691,14 @@ void Mesh::CommunicateBoundaries(std::string md_name) { "Too many iterations waiting to receive boundary communication buffers."); for (auto &partition : partitions) { - auto &md = mesh_data.Add(md_name, partition); + auto &md = mesh_data.Add(md_name, partition, fields); // unpack FillGhost variables SetBoundaries(md); } // Now do prolongation, compute primitives, apply BCs for (auto &partition : partitions) { - auto &md = mesh_data.Add(md_name, partition); + auto &md = mesh_data.Add(md_name, partition, fields); if (multilevel) { ApplyBoundaryConditionsOnCoarseOrFineMD(md, true); ProlongateBoundaries(md); diff --git a/src/mesh/mesh.hpp b/src/mesh/mesh.hpp index 7e500eb3cf59..86b46f461e52 100644 --- a/src/mesh/mesh.hpp +++ b/src/mesh/mesh.hpp @@ -344,7 +344,8 @@ class Mesh { void SetupMPIComms(); void BuildTagMapAndBoundaryBuffers(); - void CommunicateBoundaries(std::string md_name = "base"); + void CommunicateBoundaries(std::string md_name = "base", + const std::vector &fields = {}); void PreCommFillDerived(); void FillDerived(); From 4ee02c5fec1ae07162a5d0a11017fe03f3c406cc Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Wed, 9 Oct 2024 21:05:48 -0600 Subject: [PATCH 35/49] changelog --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 4be253c63f25..1a56e5c221cc 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,6 +12,7 @@ - [[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/...) +- [[PR 1187]](https://github.com/parthenon-hpc-lab/parthenon/pull/1187) Make DataCollection::Add safer and generalize MeshBlockData::Initialize - [[PR 1186]](https://github.com/parthenon-hpc-lab/parthenon/pull/1186) Bump Kokkos submodule to 4.4.1 - [[PR1177]](https://github.com/parthenon-hpc-lab/parthenon/pull/1177) Make mesh-level boundary conditions usable without the "user" flag - [[PR 1171]](https://github.com/parthenon-hpc-lab/parthenon/pull/1171) Add PARTHENON_USE_SYSTEM_PACKAGES build option From 862c2a98f0f4e112e5c9ef422cd3127a9a93ed75 Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Thu, 10 Oct 2024 09:26:18 -0600 Subject: [PATCH 36/49] Make everything work --- src/interface/meshblock_data.hpp | 15 +++++---------- src/interface/state_descriptor.hpp | 6 ++++++ 2 files changed, 11 insertions(+), 10 deletions(-) diff --git a/src/interface/meshblock_data.hpp b/src/interface/meshblock_data.hpp index 2e95800c2e79..a6b681baf142 100644 --- a/src/interface/meshblock_data.hpp +++ b/src/interface/meshblock_data.hpp @@ -196,8 +196,7 @@ class MeshBlockData { if (!found) add_var(src->GetVarPtr(flx_name)); } } - } else if constexpr (std::is_same_v && - std::is_same_v) { + } else if constexpr (std::is_same_v) { for (const auto &v : vars) { const auto &vid = resolved_packages->GetFieldVarID(v); const auto &md = resolved_packages->GetFieldMetadata(v); @@ -205,20 +204,16 @@ class MeshBlockData { // Add the associated flux as well if not explicitly // asked for if (md.IsSet(Metadata::WithFluxes)) { - auto flx_name = md.GetFluxName(); + auto flx_vid = resolved_packages->GetFieldVarID(md.GetFluxName()); bool found = false; for (const auto &v2 : vars) - if (v2 == flx_name) found = true; + if (resolved_packages->GetFieldVarID(v2) == flx_vid) found = true; if (!found) { - const auto &vid = resolved_packages->GetFieldVarID(flx_name); - const auto &md = resolved_packages->GetFieldMetadata(flx_name); - AddField(vid.base_name, md, vid.sparse_id); + const auto &flx_md = resolved_packages->GetFieldMetadata(flx_vid); + AddField(flx_vid.base_name, flx_md, flx_vid.sparse_id); } } } - } else { - PARTHENON_FAIL("Variable subset selection not yet implemented for MeshBlock " - "input with unique ids."); } } diff --git a/src/interface/state_descriptor.hpp b/src/interface/state_descriptor.hpp index ff143605b2a1..d1af2b85f2dd 100644 --- a/src/interface/state_descriptor.hpp +++ b/src/interface/state_descriptor.hpp @@ -203,6 +203,12 @@ class StateDescriptor { // retrieve all swarm names std::vector Swarms() noexcept; + const auto GetFieldVarID(const VarID &id) const { + PARTHENON_REQUIRE_THROWS(metadataMap_.count(id), + "Asking for a variable that is not in this StateDescriptor."); + return id; + } + const auto &GetFieldVarID(const std::string &label) const { return labelToVidMap_.at(label); } From bd50d2447de4324f2d4ee5289d832040dd75a4f5 Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Thu, 10 Oct 2024 09:28:45 -0600 Subject: [PATCH 37/49] format --- src/interface/state_descriptor.hpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/interface/state_descriptor.hpp b/src/interface/state_descriptor.hpp index d1af2b85f2dd..a516f6be8821 100644 --- a/src/interface/state_descriptor.hpp +++ b/src/interface/state_descriptor.hpp @@ -204,8 +204,9 @@ class StateDescriptor { std::vector Swarms() noexcept; const auto GetFieldVarID(const VarID &id) const { - PARTHENON_REQUIRE_THROWS(metadataMap_.count(id), - "Asking for a variable that is not in this StateDescriptor."); + PARTHENON_REQUIRE_THROWS( + metadataMap_.count(id), + "Asking for a variable that is not in this StateDescriptor."); return id; } From 09464c25efbc32dea26b4a677aaeb30fc6a7ee52 Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Thu, 10 Oct 2024 10:47:18 -0600 Subject: [PATCH 38/49] maybe fix doc issue? --- .github/workflows/docs.yml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml index 3ada901a78ac..4b95511d82f5 100644 --- a/.github/workflows/docs.yml +++ b/.github/workflows/docs.yml @@ -23,9 +23,9 @@ jobs: run: export DEBIAN_FRONTEND=noninteractive - name: install dependencies run: | - pip install sphinx - pip install sphinx-rtd-theme - pip install sphinx-multiversion + pip install --break-system-packages sphinx + pip install --break-system-packages sphinx-rtd-theme + pip install --break-system-packages sphinx-multiversion - name: build docs run: | echo "Repo = ${GITHUB_REPOSITORY}" From df06a508cdbdded0f007fc06a60e55340fa7675e Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Wed, 16 Oct 2024 01:30:55 +0200 Subject: [PATCH 39/49] Address CUDA MPI/ICP issue with Kokkos <=4.4.1 (#1189) --- CHANGELOG.md | 2 ++ CMakeLists.txt | 7 +++++++ 2 files changed, 9 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 1a56e5c221cc..0682db869c5a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -15,10 +15,12 @@ - [[PR 1187]](https://github.com/parthenon-hpc-lab/parthenon/pull/1187) Make DataCollection::Add safer and generalize MeshBlockData::Initialize - [[PR 1186]](https://github.com/parthenon-hpc-lab/parthenon/pull/1186) Bump Kokkos submodule to 4.4.1 - [[PR1177]](https://github.com/parthenon-hpc-lab/parthenon/pull/1177) Make mesh-level boundary conditions usable without the "user" flag +- [[Issue 1165]](https://github.com/parthenon-hpc-lab/parthenon/issues/1165) Bump Kokkos submodule to 4.4.1 - [[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/...) +- [[PR 1189]](https://github.com/parthenon-hpc-lab/parthenon/pull/1189) Address CUDA MPI/ICP issue with Kokkos <=4.4.1 - [[PR 1178]](https://github.com/parthenon-hpc-lab/parthenon/pull/1178) Fix issue with mesh pointer when using relative residual tolerance in BiCGSTAB solver. - [[PR 1173]](https://github.com/parthenon-hpc-lab/parthenon/pull/1173) Make debugging easier by making parthenon throw an error if ParameterInput is different on multiple MPI ranks. diff --git a/CMakeLists.txt b/CMakeLists.txt index 24a7a3b2ebfb..61fd4fd3cfff 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -291,6 +291,13 @@ if (Kokkos_ENABLE_CUDA) if(CHECK_REGISTRY_PRESSURE) add_compile_options(-Xptxas=-v) endif() + + # Async malloc are enabled by default until Kokkos <= 4.4.1 but + # cause issues with IPCs and MPI, + # see https://github.com/parthenon-hpc-lab/athenapk/pull/114#issuecomment-2358981857 + # Given that it's also not expected to be a significant performance gain + # it's hard disabled now. + set(Kokkos_ENABLE_IMPL_CUDA_MALLOC_ASYNC OFF CACHE BOOL "Disable Cuda async malloc (to address MPI/IPC issues)") endif() # Note that these options may not play nice with nvcc wrapper if (TEST_INTEL_OPTIMIZATION) From 511c77cc07077b1df7c09011fc00da62bc783fb7 Mon Sep 17 00:00:00 2001 From: Ben Ryan Date: Thu, 31 Oct 2024 10:54:09 -0600 Subject: [PATCH 40/49] Jonah's fix for this CI issue --- .github/workflows/docs.yml | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml index 4b95511d82f5..9e1cbcfad76f 100644 --- a/.github/workflows/docs.yml +++ b/.github/workflows/docs.yml @@ -12,9 +12,12 @@ on: jobs: docs: name: build and deploy docs - runs-on: ubuntu-latest + # not using latest due to issues with pip user packages, see + # https://github.com/actions/runner-images/issues/10781 and + # https://github.com/actions/runner-images/issues/10636 + runs-on: ubuntu-22.04 - steps: + steps: - name: Checkout code uses: actions/checkout@v2 with: @@ -23,9 +26,9 @@ jobs: run: export DEBIAN_FRONTEND=noninteractive - name: install dependencies run: | - pip install --break-system-packages sphinx - pip install --break-system-packages sphinx-rtd-theme - pip install --break-system-packages sphinx-multiversion + pip install sphinx + pip install sphinx-rtd-theme + pip install sphinx-multiversion - name: build docs run: | echo "Repo = ${GITHUB_REPOSITORY}" From 328dc33ef67d369860f375d60a8f239604ffbdd5 Mon Sep 17 00:00:00 2001 From: Ben Ryan Date: Thu, 31 Oct 2024 10:55:43 -0600 Subject: [PATCH 41/49] CHANGELOG --- CHANGELOG.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 0682db869c5a..5466bbb4f44e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,6 +12,8 @@ - [[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/...) +- [[PR1203]](https://github.com/parthenon-hpc-lab/parthenon/pull/1203) Pin Ubuntu CI image +- [[PR1177]](https://github.com/parthenon-hpc-lab/parthenon/pull/1177) Make mesh-level boundary conditions usable without the "user" flag - [[PR 1187]](https://github.com/parthenon-hpc-lab/parthenon/pull/1187) Make DataCollection::Add safer and generalize MeshBlockData::Initialize - [[PR 1186]](https://github.com/parthenon-hpc-lab/parthenon/pull/1186) Bump Kokkos submodule to 4.4.1 - [[PR1177]](https://github.com/parthenon-hpc-lab/parthenon/pull/1177) Make mesh-level boundary conditions usable without the "user" flag From 5e657c0baa7b6a8ed49c0d6b0edba256124df200 Mon Sep 17 00:00:00 2001 From: Adam Dempsey Date: Thu, 31 Oct 2024 07:38:34 -0600 Subject: [PATCH 42/49] Remove else from if constexpr when there are returns --- src/utils/hash.hpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/utils/hash.hpp b/src/utils/hash.hpp index 77642a64dba3..4f1ff6557c74 100644 --- a/src/utils/hash.hpp +++ b/src/utils/hash.hpp @@ -31,9 +31,8 @@ std::size_t hash_combine(std::size_t lhs, const T &v, Rest &&...rest) { lhs ^= rhs + 0x9e3779b9 + (lhs << 6) + (lhs >> 2); if constexpr (sizeof...(Rest) > 0) { return hash_combine(lhs, std::forward(rest)...); - } else { - return lhs; } + return lhs; } template ::value - 1> From dd928abef0f456ca39a11d066905407ea451b200 Mon Sep 17 00:00:00 2001 From: Alex Long Date: Thu, 31 Oct 2024 22:35:00 -0600 Subject: [PATCH 43/49] Consolidate buffer packing functions with less atomics (#1199) * Consolidate buffer packing functions with less atomics + Fold in buffer sizing to the load buffer function + Use a sort and a discontinuity check instead of an atomic_fetch_add inside of LoadBuffer_ to get particle index into buffer + Reduce transcendental functions call in particle sourcing in particle example * Address PR comments + Call the member variable in SwarmKey the sort_key + Remove CountParticlesInBuffer function + Add buffer_start and buffer_sorted as swarm member variables * Update example/particles/particles.cpp Co-authored-by: Ben Ryan * Update src/interface/swarm_comms.cpp Co-authored-by: Ben Ryan --------- Co-authored-by: Ben Ryan Co-authored-by: Ben Ryan --- example/particles/particles.cpp | 30 ++-- src/interface/swarm.cpp | 23 +-- src/interface/swarm.hpp | 5 +- src/interface/swarm_comms.cpp | 190 ++++++++++++------------- src/interface/swarm_device_context.hpp | 7 +- 5 files changed, 131 insertions(+), 124 deletions(-) diff --git a/example/particles/particles.cpp b/example/particles/particles.cpp index a01a121f75fa..0f250e9ab02f 100644 --- a/example/particles/particles.cpp +++ b/example/particles/particles.cpp @@ -288,17 +288,18 @@ TaskStatus CreateSomeParticles(MeshBlock *pmb, const double t0) { } while (r > 0.5); // Randomly sample direction perpendicular to origin - Real theta = acos(2. * rng_gen.drand() - 1.); - Real phi = 2. * M_PI * rng_gen.drand(); - v(0, n) = sin(theta) * cos(phi); - v(1, n) = sin(theta) * sin(phi); - v(2, n) = cos(theta); + const Real mu = 2.0 * rng_gen.drand() - 1.0; + const Real phi = 2. * M_PI * rng_gen.drand(); + const Real stheta = std::sqrt(1.0 - mu * mu); + v(0, n) = stheta * cos(phi); + v(1, n) = stheta * sin(phi); + v(2, n) = mu; // Project v onto plane normal to sphere Real vdN = v(0, n) * x(n) + v(1, n) * y(n) + v(2, n) * z(n); - Real NdN = r * r; - v(0, n) = v(0, n) - vdN / NdN * x(n); - v(1, n) = v(1, n) - vdN / NdN * y(n); - v(2, n) = v(2, n) - vdN / NdN * z(n); + Real inverse_NdN = 1. / (r * r); + v(0, n) = v(0, n) - vdN * inverse_NdN * x(n); + v(1, n) = v(1, n) - vdN * inverse_NdN * y(n); + v(2, n) = v(2, n) - vdN * inverse_NdN * z(n); // Normalize Real v_tmp = sqrt(v(0, n) * v(0, n) + v(1, n) * v(1, n) + v(2, n) * v(2, n)); @@ -335,11 +336,12 @@ TaskStatus CreateSomeParticles(MeshBlock *pmb, const double t0) { z(n) = minx_k + nx_k * dx_k * rng_gen.drand(); // Randomly sample direction on the unit sphere, fixing speed - Real theta = acos(2. * rng_gen.drand() - 1.); - Real phi = 2. * M_PI * rng_gen.drand(); - v(0, n) = vel * sin(theta) * cos(phi); - v(1, n) = vel * sin(theta) * sin(phi); - v(2, n) = vel * cos(theta); + const Real mu = 2.0 * rng_gen.drand() - 1.0; + const Real phi = 2. * M_PI * rng_gen.drand(); + const Real stheta = std::sqrt(1.0 - mu * mu); + v(0, n) = vel * stheta * cos(phi); + v(1, n) = vel * stheta * sin(phi); + v(2, n) = vel * mu; // Create particles at the beginning of the timestep t(n) = t0; diff --git a/src/interface/swarm.cpp b/src/interface/swarm.cpp index 294338dee0d2..39f1e8a28123 100644 --- a/src/interface/swarm.cpp +++ b/src/interface/swarm.cpp @@ -33,6 +33,7 @@ SwarmDeviceContext Swarm::GetDeviceContext() const { context.block_index_ = block_index_; context.neighbor_indices_ = neighbor_indices_; context.cell_sorted_ = cell_sorted_; + context.buffer_sorted_ = buffer_sorted_; context.cell_sorted_begin_ = cell_sorted_begin_; context.cell_sorted_number_ = cell_sorted_number_; @@ -73,9 +74,10 @@ Swarm::Swarm(const std::string &label, const Metadata &metadata, const int nmax_ new_indices_("new_indices_", nmax_pool_), scratch_a_("scratch_a_", nmax_pool_), scratch_b_("scratch_b_", nmax_pool_), num_particles_to_send_("num_particles_to_send_", NMAX_NEIGHBORS), - buffer_counters_("buffer_counters_", NMAX_NEIGHBORS), + buffer_start_("buffer_start_", NMAX_NEIGHBORS), neighbor_received_particles_("neighbor_received_particles_", NMAX_NEIGHBORS), - cell_sorted_("cell_sorted_", nmax_pool_), mpiStatus(true) { + cell_sorted_("cell_sorted_", nmax_pool_), + buffer_sorted_("buffer_sorted_", nmax_pool_), mpiStatus(true) { PARTHENON_REQUIRE_THROWS(typeid(Coordinates_t) == typeid(UniformCartesian), "SwarmDeviceContext only supports a uniform Cartesian mesh!"); @@ -209,6 +211,9 @@ void Swarm::SetPoolMax(const std::int64_t nmax_pool) { Kokkos::resize(cell_sorted_, nmax_pool); pmb->LogMemUsage(n_new * sizeof(SwarmKey)); + Kokkos::resize(buffer_sorted_, nmax_pool); + pmb->LogMemUsage(n_new * sizeof(SwarmKey)); + block_index_.Resize(nmax_pool); pmb->LogMemUsage(n_new * sizeof(int)); @@ -490,35 +495,35 @@ void Swarm::SortParticlesByCell() { break; } - if (cell_sorted(start_index).cell_idx_1d_ == cell_idx_1d) { + if (cell_sorted(start_index).sort_idx_ == cell_idx_1d) { if (start_index == 0) { break; - } else if (cell_sorted(start_index - 1).cell_idx_1d_ != cell_idx_1d) { + } else if (cell_sorted(start_index - 1).sort_idx_ != cell_idx_1d) { break; } else { start_index--; continue; } } - if (cell_sorted(start_index).cell_idx_1d_ >= cell_idx_1d) { + if (cell_sorted(start_index).sort_idx_ >= cell_idx_1d) { start_index--; if (start_index < 0) { start_index = -1; break; } - if (cell_sorted(start_index).cell_idx_1d_ < cell_idx_1d) { + if (cell_sorted(start_index).sort_idx_ < cell_idx_1d) { start_index = -1; break; } continue; } - if (cell_sorted(start_index).cell_idx_1d_ < cell_idx_1d) { + if (cell_sorted(start_index).sort_idx_ < cell_idx_1d) { start_index++; if (start_index > max_active_index) { start_index = -1; break; } - if (cell_sorted(start_index).cell_idx_1d_ > cell_idx_1d) { + if (cell_sorted(start_index).sort_idx_ > cell_idx_1d) { start_index = -1; break; } @@ -532,7 +537,7 @@ void Swarm::SortParticlesByCell() { int number = 0; int current_index = start_index; while (current_index <= max_active_index && - cell_sorted(current_index).cell_idx_1d_ == cell_idx_1d) { + cell_sorted(current_index).sort_idx_ == cell_idx_1d) { current_index++; number++; cell_sorted_number(k, j, i) = number; diff --git a/src/interface/swarm.hpp b/src/interface/swarm.hpp index 2058dc58f894..49cd906c2a15 100644 --- a/src/interface/swarm.hpp +++ b/src/interface/swarm.hpp @@ -289,7 +289,7 @@ class Swarm { constexpr static int unset_index_ = -1; ParArray1D num_particles_to_send_; - ParArray1D buffer_counters_; + ParArray1D buffer_start_; ParArray1D neighbor_received_particles_; int total_received_particles_; @@ -298,6 +298,9 @@ class Swarm { ParArray1D cell_sorted_; // 1D per-cell sorted array of key-value swarm memory indices + ParArray1D + buffer_sorted_; // 1D per-buffer sorted array of key-value swarm memory indices + ParArrayND cell_sorted_begin_; // Per-cell array of starting indices in cell_sorted_ diff --git a/src/interface/swarm_comms.cpp b/src/interface/swarm_comms.cpp index 3ee4e2af7512..7b236b6e4e1f 100644 --- a/src/interface/swarm_comms.cpp +++ b/src/interface/swarm_comms.cpp @@ -156,72 +156,14 @@ void Swarm::SetupPersistentMPI() { } } -void Swarm::CountParticlesToSend_() { - auto mask_h = Kokkos::create_mirror_view_and_copy(HostMemSpace(), mask_); - auto swarm_d = GetDeviceContext(); - auto pmb = GetBlockPointer(); - const int nbmax = vbswarm->bd_var_.nbmax; - - // Fence to make sure particles aren't currently being transported locally - // TODO(BRR) do this operation on device. - pmb->exec_space.fence(); - const int particle_size = GetParticleDataSize(); - vbswarm->particle_size = particle_size; - - // TODO(BRR) This kernel launch should be folded into the subsequent logic once we - // convert that to kernel-based reductions - auto &x = Get(swarm_position::x::name()).Get(); - auto &y = Get(swarm_position::y::name()).Get(); - auto &z = Get(swarm_position::z::name()).Get(); - const int max_active_index = GetMaxActiveIndex(); - pmb->par_for( - PARTHENON_AUTO_LABEL, 0, max_active_index, KOKKOS_LAMBDA(const int n) { - if (swarm_d.IsActive(n)) { - bool on_current_mesh_block = true; - swarm_d.GetNeighborBlockIndex(n, x(n), y(n), z(n), on_current_mesh_block); - } - }); - - // Facilitate lambda captures - auto &block_index = block_index_; - auto &num_particles_to_send = num_particles_to_send_; - - // Zero out number of particles to send before accumulating - pmb->par_for( - PARTHENON_AUTO_LABEL, 0, NMAX_NEIGHBORS - 1, - KOKKOS_LAMBDA(const int n) { num_particles_to_send[n] = 0; }); - - parthenon::par_for( - PARTHENON_AUTO_LABEL, 0, max_active_index, KOKKOS_LAMBDA(const int n) { - if (swarm_d.IsActive(n)) { - bool on_current_mesh_block = true; - swarm_d.GetNeighborBlockIndex(n, x(n), y(n), z(n), on_current_mesh_block); - - if (block_index(n) >= 0) { - Kokkos::atomic_add(&num_particles_to_send(block_index(n)), 1); - } - } - }); - - auto num_particles_to_send_h = num_particles_to_send_.GetHostMirrorAndCopy(); - - // Resize send buffers if too small - for (int n = 0; n < pmb->neighbors.size(); n++) { - const int bufid = pmb->neighbors[n].bufid; - auto sendbuf = vbswarm->bd_var_.send[bufid]; - if (sendbuf.extent(0) < num_particles_to_send_h(n) * particle_size) { - sendbuf = BufArray1D("Buffer", num_particles_to_send_h(n) * particle_size); - vbswarm->bd_var_.send[bufid] = sendbuf; - } - vbswarm->send_size[bufid] = num_particles_to_send_h(n) * particle_size; - } -} - void Swarm::LoadBuffers_() { auto swarm_d = GetDeviceContext(); auto pmb = GetBlockPointer(); const int particle_size = GetParticleDataSize(); + vbswarm->particle_size = particle_size; const int nneighbor = pmb->neighbors.size(); + // Fence to make sure particles aren't currently being transported locally + pmb->exec_space.fence(); auto &int_vector = std::get()>(vectors_); auto &real_vector = std::get()>(vectors_); @@ -236,42 +178,98 @@ void Swarm::LoadBuffers_() { auto &y = Get(swarm_position::y::name()).Get(); auto &z = Get(swarm_position::z::name()).Get(); - // Zero buffer index counters - auto &buffer_counters = buffer_counters_; - pmb->par_for( - PARTHENON_AUTO_LABEL, 0, NMAX_NEIGHBORS - 1, - KOKKOS_LAMBDA(const int n) { buffer_counters[n] = 0; }); + if (max_active_index_ >= 0) { + auto &buffer_sorted = buffer_sorted_; + auto &buffer_start = buffer_start_; - auto &bdvar = vbswarm->bd_var_; - auto neighbor_buffer_index = neighbor_buffer_index_; - // Loop over active particles and use atomic operations to find indices into buffers if - // this particle will be sent. - pmb->par_for( - PARTHENON_AUTO_LABEL, 0, max_active_index_, KOKKOS_LAMBDA(const int n) { - if (swarm_d.IsActive(n)) { - bool on_current_mesh_block = true; - const int m = - swarm_d.GetNeighborBlockIndex(n, x(n), y(n), z(n), on_current_mesh_block); - const int bufid = neighbor_buffer_index(m); - - if (m >= 0) { - const int bid = Kokkos::atomic_fetch_add(&buffer_counters(m), 1); - int buffer_index = bid * particle_size; - swarm_d.MarkParticleForRemoval(n); - for (int i = 0; i < realPackDim; i++) { - bdvar.send[bufid](buffer_index) = vreal(i, n); - buffer_index++; - } - for (int i = 0; i < intPackDim; i++) { - bdvar.send[bufid](buffer_index) = static_cast(vint(i, n)); - buffer_index++; + pmb->par_for( + PARTHENON_AUTO_LABEL, 0, max_active_index_, KOKKOS_LAMBDA(const int n) { + if (swarm_d.IsActive(n)) { + bool on_current_mesh_block = true; + const int m = + swarm_d.GetNeighborBlockIndex(n, x(n), y(n), z(n), on_current_mesh_block); + buffer_sorted(n) = SwarmKey(m, n); + } else { + buffer_sorted(n) = SwarmKey(this_block_, n); + } + }); + + // sort by buffer index + sort(buffer_sorted, SwarmKeyComparator(), 0, max_active_index_); + + // use discontinuity check to determine start of a buffer in buffer_sorted array + auto &num_particles_to_send = num_particles_to_send_; + auto max_active_index = max_active_index_; + + // Zero out number of particles to send before accumulating + pmb->par_for( + PARTHENON_AUTO_LABEL, 0, NMAX_NEIGHBORS - 1, KOKKOS_LAMBDA(const int n) { + num_particles_to_send[n] = 0; + buffer_start[n] = 0; + }); + + pmb->par_for( + PARTHENON_AUTO_LABEL, 0, max_active_index_, KOKKOS_LAMBDA(const int n) { + auto m = buffer_sorted(n).sort_idx_; + // start checks (used for index of particle in buffer) + if (m >= 0 && n == 0) { + buffer_start(m) = 0; + } else if (m >= 0 && m != buffer_sorted(n - 1).sort_idx_) { + buffer_start(m) = n; + } + // end checks (used to to size particle buffers) + if (m >= 0 && n == max_active_index) { + num_particles_to_send(m) = n + 1; + } else if (m >= 0 && m != buffer_sorted(n + 1).sort_idx_) { + num_particles_to_send(m) = n + 1; + } + }); + + // copy values back to host for buffer sizing + auto num_particles_to_send_h = num_particles_to_send_.GetHostMirrorAndCopy(); + auto buffer_start_h = buffer_start.GetHostMirrorAndCopy(); + + // Resize send buffers if too small + for (int n = 0; n < pmb->neighbors.size(); n++) { + num_particles_to_send_h(n) -= buffer_start_h(n); + const int bufid = pmb->neighbors[n].bufid; + auto sendbuf = vbswarm->bd_var_.send[bufid]; + if (sendbuf.extent(0) < num_particles_to_send_h(n) * particle_size) { + sendbuf = BufArray1D("Buffer", num_particles_to_send_h(n) * particle_size); + vbswarm->bd_var_.send[bufid] = sendbuf; + } + vbswarm->send_size[bufid] = num_particles_to_send_h(n) * particle_size; + } + + auto &bdvar = vbswarm->bd_var_; + auto neighbor_buffer_index = neighbor_buffer_index_; + // Loop over active particles buffer_sorted, use index n and buffer_start to + /// get index in buffer m, pack that particle in buffer + pmb->par_for( + PARTHENON_AUTO_LABEL, 0, max_active_index_, KOKKOS_LAMBDA(const int n) { + auto p_index = buffer_sorted(n).swarm_idx_; + if (swarm_d.IsActive(p_index)) { + const int m = buffer_sorted(n).sort_idx_; + const int bufid = neighbor_buffer_index(m); + if (m >= 0) { + const int bid = n - buffer_start[m]; + int buffer_index = bid * particle_size; + swarm_d.MarkParticleForRemoval(p_index); + for (int i = 0; i < realPackDim; i++) { + bdvar.send[bufid](buffer_index) = vreal(i, p_index); + buffer_index++; + } + for (int i = 0; i < intPackDim; i++) { + bdvar.send[bufid](buffer_index) = static_cast(vint(i, p_index)); + buffer_index++; + } } } - } - }); + }); - // Remove particles that were loaded to send to another block from this block - RemoveMarkedParticles(); + // Remove particles that were loaded to send to another block from this block + RemoveMarkedParticles(); + } } void Swarm::Send(BoundaryCommSubset phase) { @@ -279,10 +277,8 @@ void Swarm::Send(BoundaryCommSubset phase) { const int nneighbor = pmb->neighbors.size(); auto swarm_d = GetDeviceContext(); - // Query particles for those to be sent - CountParticlesToSend_(); - - // Prepare buffers for send operations + // Potentially resize buffer, get consistent index from particle array, get ready to + // send LoadBuffers_(); // Send buffer data diff --git a/src/interface/swarm_device_context.hpp b/src/interface/swarm_device_context.hpp index 936d2d56ad35..b16daf0d0cdf 100644 --- a/src/interface/swarm_device_context.hpp +++ b/src/interface/swarm_device_context.hpp @@ -25,16 +25,16 @@ struct SwarmKey { SwarmKey() {} KOKKOS_INLINE_FUNCTION SwarmKey(const int cell_idx_1d, const int swarm_idx_1d) - : cell_idx_1d_(cell_idx_1d), swarm_idx_(swarm_idx_1d) {} + : sort_idx_(cell_idx_1d), swarm_idx_(swarm_idx_1d) {} - int cell_idx_1d_; + int sort_idx_; int swarm_idx_; }; struct SwarmKeyComparator { KOKKOS_INLINE_FUNCTION bool operator()(const SwarmKey &s1, const SwarmKey &s2) { - return s1.cell_idx_1d_ < s2.cell_idx_1d_; + return s1.sort_idx_ < s2.sort_idx_; } }; @@ -139,6 +139,7 @@ class SwarmDeviceContext { ParArrayND block_index_; ParArrayND neighbor_indices_; // 4x4x4 array of possible block AMR regions ParArray1D cell_sorted_; + ParArray1D buffer_sorted_; ParArrayND cell_sorted_begin_; ParArrayND cell_sorted_number_; int ndim_; From 0cbbdacac5c234e57d52636db64e3180e79edee7 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Mon, 4 Nov 2024 18:55:01 +0100 Subject: [PATCH 44/49] [Trivial] Fix type used for array init (#1170) * Fix type used for array init * init array with constexpr expression * CC --------- Co-authored-by: Jonah Miller --- CHANGELOG.md | 1 + src/outputs/output_utils.hpp | 2 +- src/outputs/restart_hdf5.cpp | 4 ++-- 3 files changed, 4 insertions(+), 3 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 5466bbb4f44e..57095203cc3a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -22,6 +22,7 @@ - [[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/...) +- [[PR 1170]](https://github.com/parthenon-hpc-lab/parthenon/pull/1170) Fixed incorrect initialization of array by a const not constexpr - [[PR 1189]](https://github.com/parthenon-hpc-lab/parthenon/pull/1189) Address CUDA MPI/ICP issue with Kokkos <=4.4.1 - [[PR 1178]](https://github.com/parthenon-hpc-lab/parthenon/pull/1178) Fix issue with mesh pointer when using relative residual tolerance in BiCGSTAB solver. - [[PR 1173]](https://github.com/parthenon-hpc-lab/parthenon/pull/1173) Make debugging easier by making parthenon throw an error if ParameterInput is different on multiple MPI ranks. diff --git a/src/outputs/output_utils.hpp b/src/outputs/output_utils.hpp index 5f95fbbaedc2..5c094ef688d7 100644 --- a/src/outputs/output_utils.hpp +++ b/src/outputs/output_utils.hpp @@ -50,7 +50,7 @@ namespace OutputUtils { // Helper struct containing some information about a variable struct VarInfo { public: - static constexpr int VNDIM = MAX_VARIABLE_DIMENSION; + static constexpr const int VNDIM = MAX_VARIABLE_DIMENSION; std::string label; int num_components; int tensor_rank; // 0- to 3-D for cell-centered variables, 0- to 6-D for arbitrary shape diff --git a/src/outputs/restart_hdf5.cpp b/src/outputs/restart_hdf5.cpp index cababddbfafe..8078fbd5fa9a 100644 --- a/src/outputs/restart_hdf5.cpp +++ b/src/outputs/restart_hdf5.cpp @@ -3,7 +3,7 @@ // Copyright(C) 2020-2024 The Parthenon collaboration // Licensed under the 3-clause BSD License, see LICENSE file for details //======================================================================================== -// (C) (or copyright) 2020-2021. 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 @@ -213,7 +213,7 @@ void RestartReaderHDF5::ReadBlocks(const std::string &name, IndexRange range, #else // HDF5 enabled auto hdl = OpenDataset(name); - const int VNDIM = info.VNDIM; + constexpr int VNDIM = OutputUtils::VarInfo::VNDIM; /** Select hyperslab in dataset **/ int total_dim = 0; From 01517f7f2a0be3857a3eed7461a35c21749e64b9 Mon Sep 17 00:00:00 2001 From: Ben Ryan Date: Thu, 7 Nov 2024 20:43:27 -0700 Subject: [PATCH 45/49] Leapfrog fix (#1206) * Missing send size init * cleanup, CHANGELOG * verbose CI * further CI debugging * This should be working... * This should be fixed... but I get a segfault on GPU * Is it my AMD GPU thats wrong? * Missing a return statement * retest * Oops missing statement * Revert test * revert workflow --- CHANGELOG.md | 1 + src/interface/swarm_comms.cpp | 5 +++++ src/utils/sort.hpp | 17 +++++++++++++++-- 3 files changed, 21 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 57095203cc3a..43671a3b06b9 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,6 +12,7 @@ - [[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/...) +- [[PR 1206]](https://github.com/parthenon-hpc-lab/parthenon/pull/1206) Leapfrog fix - [[PR1203]](https://github.com/parthenon-hpc-lab/parthenon/pull/1203) Pin Ubuntu CI image - [[PR1177]](https://github.com/parthenon-hpc-lab/parthenon/pull/1177) Make mesh-level boundary conditions usable without the "user" flag - [[PR 1187]](https://github.com/parthenon-hpc-lab/parthenon/pull/1187) Make DataCollection::Add safer and generalize MeshBlockData::Initialize diff --git a/src/interface/swarm_comms.cpp b/src/interface/swarm_comms.cpp index 7b236b6e4e1f..054ded34fabb 100644 --- a/src/interface/swarm_comms.cpp +++ b/src/interface/swarm_comms.cpp @@ -269,6 +269,11 @@ void Swarm::LoadBuffers_() { // Remove particles that were loaded to send to another block from this block RemoveMarkedParticles(); + } else { + for (int n = 0; n < pmb->neighbors.size(); n++) { + const int bufid = pmb->neighbors[n].bufid; + vbswarm->send_size[bufid] = 0; + } } } diff --git a/src/utils/sort.hpp b/src/utils/sort.hpp index 97e9c77a88e4..a4ab8139dab3 100644 --- a/src/utils/sort.hpp +++ b/src/utils/sort.hpp @@ -61,7 +61,7 @@ void sort(ParArray1D data, KeyComparator comparator, size_t min_idx, size_t max_idx) { PARTHENON_DEBUG_REQUIRE(min_idx < data.extent(0), "Invalid minimum sort index!"); PARTHENON_DEBUG_REQUIRE(max_idx < data.extent(0), "Invalid maximum sort index!"); -#ifdef KOKKOS_ENABLE_CUDA +#if defined(KOKKOS_ENABLE_CUDA) #ifdef __clang__ PARTHENON_FAIL("sort is using thrust and there exists an incompatibility with clang, " "see https://github.com/lanl/parthenon/issues/647 for more details. We " @@ -74,6 +74,13 @@ void sort(ParArray1D data, KeyComparator comparator, size_t min_idx, thrust::device_ptr last_d = thrust::device_pointer_cast(data.data()) + max_idx + 1; thrust::sort(first_d, last_d, comparator); #endif +#elif defined(KOKKOS_ENABLE_HIP) + auto data_h = Kokkos::create_mirror_view_and_copy(HostMemSpace(), data); + std::sort(data_h.data() + min_idx, data_h.data() + max_idx + 1, comparator); + Kokkos::deep_copy(data, data_h); + // TODO(BRR) With Kokkos 4.4, switch to Kokkos::sort + // auto sub_data = Kokkos::subview(data, std::make_pair(min_idx, max_idx + 1)); + // Kokkos::sort(sub_data, comparator); #else if (std::is_same::value) { std::sort(data.data() + min_idx, data.data() + max_idx + 1, comparator); @@ -89,7 +96,7 @@ template void sort(ParArray1D data, size_t min_idx, size_t max_idx) { PARTHENON_DEBUG_REQUIRE(min_idx < data.extent(0), "Invalid minimum sort index!"); PARTHENON_DEBUG_REQUIRE(max_idx < data.extent(0), "Invalid maximum sort index!"); -#ifdef KOKKOS_ENABLE_CUDA +#if defined(KOKKOS_ENABLE_CUDA) || defined(KOKKOS_ENABLE_HIP) #ifdef __clang__ PARTHENON_FAIL("sort is using thrust and there exists an incompatibility with clang, " "see https://github.com/lanl/parthenon/issues/647 for more details. We " @@ -102,6 +109,12 @@ void sort(ParArray1D data, size_t min_idx, size_t max_idx) { thrust::device_ptr last_d = thrust::device_pointer_cast(data.data()) + max_idx + 1; thrust::sort(first_d, last_d); #endif + auto data_h = Kokkos::create_mirror_view_and_copy(HostMemSpace(), data); + std::sort(data_h.data() + min_idx, data_h.data() + max_idx + 1); + Kokkos::deep_copy(data, data_h); + // TODO(BRR) With Kokkos 4.4, switch to Kokkos::sort + // auto sub_data = Kokkos::subview(data, std::make_pair(min_idx, max_idx + 1)); + // Kokkos::sort(sub_data); #else if (std::is_same::value) { std::sort(data.data() + min_idx, data.data() + max_idx + 1); From 3ddbb9f84a54bc65b7e8ea2b248bd00a7fc9554e Mon Sep 17 00:00:00 2001 From: adam reyes Date: Fri, 8 Nov 2024 15:38:52 -0800 Subject: [PATCH 46/49] removing meshblock amr_criteria --- example/fine_advection/advection_package.cpp | 51 +++----------------- src/amr_criteria/amr_criteria.cpp | 27 ----------- src/amr_criteria/amr_criteria.hpp | 4 -- src/amr_criteria/refinement_package.cpp | 30 +----------- src/amr_criteria/refinement_package.hpp | 6 --- src/mesh/mesh_refinement.hpp | 1 - 6 files changed, 9 insertions(+), 110 deletions(-) diff --git a/example/fine_advection/advection_package.cpp b/example/fine_advection/advection_package.cpp index 51c86a82bcab..6762f2c7e19e 100644 --- a/example/fine_advection/advection_package.cpp +++ b/example/fine_advection/advection_package.cpp @@ -94,13 +94,7 @@ std::shared_ptr Initialize(ParameterInput *pin) { pkg->AddField( Metadata({Metadata::Cell, Metadata::Derived, Metadata::OneCopy})); - bool check_refine_mesh = - pin->GetOrAddBoolean("parthenon/mesh", "CheckRefineMesh", true); - if (check_refine_mesh) { - pkg->CheckRefinementMesh = CheckRefinementMesh; - } else { - pkg->CheckRefinementBlock = CheckRefinement; - } + pkg->CheckRefinementMesh = CheckRefinementMesh; pkg->EstimateTimestepMesh = EstimateTimestep; pkg->FillDerivedMesh = FillDerived; return pkg; @@ -129,10 +123,12 @@ void CheckRefinementMesh(MeshData *md, parthenon::ParArray1D &amr_ parthenon::inner_loop_pattern_ttr_tag, team_member, jb.s, jb.e, ib.s, ib.e, [&](const int j, const int i, typename Kokkos::MinMax::value_type &lminmax) { - lminmax.min_val = (pack(n, k, j, i) < lminmax.min_val ? pack(n, k, j, i) - : lminmax.min_val); - lminmax.max_val = (pack(n, k, j, i) > lminmax.max_val ? pack(n, k, j, i) - : lminmax.max_val); + lminmax.min_val = + (pack(b, n, k, j, i) < lminmax.min_val ? pack(b, n, k, j, i) + : lminmax.min_val); + lminmax.max_val = + (pack(b, n, k, j, i) > lminmax.max_val ? pack(b, n, k, j, i) + : lminmax.max_val); }, Kokkos::MinMax(minmax)); @@ -146,39 +142,6 @@ void CheckRefinementMesh(MeshData *md, parthenon::ParArray1D &amr_ amr_tags.ContributeScatter(scatter_tags); } -AmrTag CheckRefinement(MeshBlockData *rc) { - // refine on advected, for example. could also be a derived quantity - static auto desc = parthenon::MakePackDescriptor(rc); - auto pack = desc.GetPack(rc); - - auto pmb = rc->GetBlockPointer(); - IndexRange ib = pmb->cellbounds.GetBoundsI(IndexDomain::entire); - IndexRange jb = pmb->cellbounds.GetBoundsJ(IndexDomain::entire); - IndexRange kb = pmb->cellbounds.GetBoundsK(IndexDomain::entire); - - typename Kokkos::MinMax::value_type minmax; - parthenon::par_reduce( - parthenon::loop_pattern_mdrange_tag, PARTHENON_AUTO_LABEL, DevExecSpace(), - pack.GetLowerBoundHost(0), pack.GetUpperBoundHost(0), kb.s, kb.e, jb.s, jb.e, ib.s, - ib.e, - KOKKOS_LAMBDA(const int n, const int k, const int j, const int i, - typename Kokkos::MinMax::value_type &lminmax) { - lminmax.min_val = - (pack(n, k, j, i) < lminmax.min_val ? pack(n, k, j, i) : lminmax.min_val); - lminmax.max_val = - (pack(n, k, j, i) > lminmax.max_val ? pack(n, k, j, i) : lminmax.max_val); - }, - Kokkos::MinMax(minmax)); - - auto pkg = pmb->packages.Get("advection_package"); - const auto &refine_tol = pkg->Param("refine_tol"); - const auto &derefine_tol = pkg->Param("derefine_tol"); - - if (minmax.max_val > refine_tol && minmax.min_val < derefine_tol) return AmrTag::refine; - if (minmax.max_val < derefine_tol) return AmrTag::derefine; - return AmrTag::same; -} - Real EstimateTimestep(MeshData *md) { std::shared_ptr pkg = md->GetMeshPointer()->packages.Get("advection_package"); diff --git a/src/amr_criteria/amr_criteria.cpp b/src/amr_criteria/amr_criteria.cpp index 76f00d7a8905..a0625b81e652 100644 --- a/src/amr_criteria/amr_criteria.cpp +++ b/src/amr_criteria/amr_criteria.cpp @@ -78,23 +78,6 @@ std::shared_ptr AMRCriteria::MakeAMRCriteria(std::string &criteria, block_name + ": " + criteria); } -AMRBounds AMRCriteria::GetBounds(const MeshBlockData *rc) const { - auto ib = rc->GetBoundsI(IndexDomain::interior); - auto jb = rc->GetBoundsJ(IndexDomain::interior); - auto kb = rc->GetBoundsK(IndexDomain::interior); - return AMRBounds(ib, jb, kb); -} - -AmrTag AMRFirstDerivative::operator()(const MeshBlockData *rc) const { - if (!rc->HasVariable(field) || !rc->IsAllocated(field)) { - return AmrTag::same; - } - auto bnds = GetBounds(rc); - auto q = Kokkos::subview(rc->Get(field).data, comp6, comp5, comp4, Kokkos::ALL(), - Kokkos::ALL(), Kokkos::ALL()); - return Refinement::FirstDerivative(bnds, q, refine_criteria, derefine_criteria); -} - void AMRFirstDerivative::operator()(MeshData *md, ParArray1D &amr_tags) const { auto ib = md->GetBoundsI(IndexDomain::interior); @@ -115,16 +98,6 @@ void AMRFirstDerivative::operator()(MeshData *md, derefine_criteria, max_level); } -AmrTag AMRSecondDerivative::operator()(const MeshBlockData *rc) const { - if (!rc->HasVariable(field) || !rc->IsAllocated(field)) { - return AmrTag::same; - } - auto bnds = GetBounds(rc); - auto q = Kokkos::subview(rc->Get(field).data, comp6, comp5, comp4, Kokkos::ALL(), - Kokkos::ALL(), Kokkos::ALL()); - return Refinement::SecondDerivative(bnds, q, refine_criteria, derefine_criteria); -} - void AMRSecondDerivative::operator()(MeshData *md, ParArray1D &amr_tags) const { auto ib = md->GetBoundsI(IndexDomain::interior); diff --git a/src/amr_criteria/amr_criteria.hpp b/src/amr_criteria/amr_criteria.hpp index 380639aa3e5d..e69aebe84733 100644 --- a/src/amr_criteria/amr_criteria.hpp +++ b/src/amr_criteria/amr_criteria.hpp @@ -36,7 +36,6 @@ struct AMRBounds { struct AMRCriteria { AMRCriteria(ParameterInput *pin, std::string &block_name); virtual ~AMRCriteria() {} - virtual AmrTag operator()(const MeshBlockData *rc) const = 0; virtual void operator()(MeshData *md, ParArray1D &delta_level) const = 0; std::string field; Real refine_criteria, derefine_criteria; @@ -44,20 +43,17 @@ struct AMRCriteria { int comp6, comp5, comp4; static std::shared_ptr MakeAMRCriteria(std::string &criteria, ParameterInput *pin, std::string &block_name); - AMRBounds GetBounds(const MeshBlockData *rc) const; }; struct AMRFirstDerivative : public AMRCriteria { AMRFirstDerivative(ParameterInput *pin, std::string &block_name) : AMRCriteria(pin, block_name) {} - AmrTag operator()(const MeshBlockData *rc) const override; void operator()(MeshData *md, ParArray1D &delta_level) const override; }; struct AMRSecondDerivative : public AMRCriteria { AMRSecondDerivative(ParameterInput *pin, std::string &block_name) : AMRCriteria(pin, block_name) {} - AmrTag operator()(const MeshBlockData *rc) const override; void operator()(MeshData *md, ParArray1D &delta_level) const override; }; diff --git a/src/amr_criteria/refinement_package.cpp b/src/amr_criteria/refinement_package.cpp index 27d4dad213bb..76d963b49be7 100644 --- a/src/amr_criteria/refinement_package.cpp +++ b/src/amr_criteria/refinement_package.cpp @@ -36,9 +36,6 @@ namespace Refinement { std::shared_ptr Initialize(ParameterInput *pin) { auto ref = std::make_shared("Refinement"); - bool check_refine_mesh = - pin->GetOrAddBoolean("parthenon/mesh", "CheckRefineMesh", true); - ref->AddParam("check_refine_mesh", check_refine_mesh); int numcrit = 0; while (true) { @@ -60,17 +57,12 @@ ParArray1D CheckAllRefinement(MeshData *md) { auto amr_tags = pm->GetAmrTags(); Kokkos::deep_copy(amr_tags.KokkosView(), AmrTag::derefine); - static const bool check_refine_mesh = - pm->packages.Get("Refinement")->Param("check_refine_mesh"); - for (auto &pkg : pm->packages.AllPackages()) { auto &desc = pkg.second; desc->CheckRefinement(md, amr_tags); - if (check_refine_mesh) { - for (auto &amr : desc->amr_criteria) { - (*amr)(md, amr_tags); - } + for (auto &amr : desc->amr_criteria) { + (*amr)(md, amr_tags); } } @@ -91,8 +83,6 @@ AmrTag CheckAllRefinement(MeshBlockData *rc, const AmrTag &level) { // neighboring blocks. Similarly for "do nothing" PARTHENON_INSTRUMENT MeshBlock *pmb = rc->GetBlockPointer(); - static const bool check_refine_mesh = - pmb->packages.Get("Refinement")->Param("check_refine_mesh"); // delta_level holds the max over all criteria. default to derefining, or level from // MeshData check. AmrTag delta_level = level; @@ -103,22 +93,6 @@ AmrTag CheckAllRefinement(MeshBlockData *rc, const AmrTag &level) { // since 1 is the max, we can return without having to look at anything else return AmrTag::refine; } - if (check_refine_mesh) continue; - // call parthenon criteria that were registered - for (auto &amr : desc->amr_criteria) { - // get the recommended change in refinement level from this criteria - AmrTag temp_delta = (*amr)(rc); - if ((temp_delta == AmrTag::refine) && pmb->loc.level() >= amr->max_level) { - // don't refine if we're at the max level - temp_delta = AmrTag::same; - } - // maintain the max across all criteria - delta_level = std::max(delta_level, temp_delta); - if (delta_level == AmrTag::refine) { - // 1 is the max, so just return - return AmrTag::refine; - } - } } return delta_level; } diff --git a/src/amr_criteria/refinement_package.hpp b/src/amr_criteria/refinement_package.hpp index cc71566fd92e..82439f947df1 100644 --- a/src/amr_criteria/refinement_package.hpp +++ b/src/amr_criteria/refinement_package.hpp @@ -40,17 +40,11 @@ TaskStatus Tag(T *rc); AmrTag CheckAllRefinement(MeshBlockData *rc, const AmrTag &level); ParArray1D CheckAllRefinement(MeshData *md); -AmrTag FirstDerivative(const AMRBounds &bnds, const ParArray3D &q, - const Real refine_criteria, const Real derefine_criteria); - void FirstDerivative(const AMRBounds &bnds, MeshData *md, const std::string &field, const int &idx, ParArray1D &amr_tags, const Real refine_criteria_, const Real derefine_criteria_, const int max_level_); -AmrTag SecondDerivative(const AMRBounds &bnds, const ParArray3D &q, - const Real refine_criteria, const Real derefine_criteria); - void SecondDerivative(const AMRBounds &bnds, MeshData *md, const std::string &field, const int &idx, ParArray1D &amr_tags, const Real refine_criteria_, const Real derefine_criteria_, diff --git a/src/mesh/mesh_refinement.hpp b/src/mesh/mesh_refinement.hpp index 800beb7d2a35..2f41effe41ec 100644 --- a/src/mesh/mesh_refinement.hpp +++ b/src/mesh/mesh_refinement.hpp @@ -53,7 +53,6 @@ class MeshRefinement { public: MeshRefinement(std::weak_ptr pmb, ParameterInput *pin); - void CheckRefinementCondition(); void SetRefinement(AmrTag flag); // setter functions for "enrolling" variable arrays in refinement via Mesh::AMR() From ed3e16f0790c358fa596c8d2ffde6d77e43f839c Mon Sep 17 00:00:00 2001 From: adam reyes Date: Fri, 15 Nov 2024 13:45:01 -0800 Subject: [PATCH 47/49] use pack.UpperBound(b) to check allocation --- example/fine_advection/advection_package.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/example/fine_advection/advection_package.cpp b/example/fine_advection/advection_package.cpp index f1d64a64f732..a134d3202ef1 100644 --- a/example/fine_advection/advection_package.cpp +++ b/example/fine_advection/advection_package.cpp @@ -145,7 +145,7 @@ void CheckRefinementMesh(MeshData *md, parthenon::ParArray1D &amr_ pack.GetMaxNumberOfVars() - 1, kb.s, kb.e, KOKKOS_LAMBDA(parthenon::team_mbr_t team_member, const int b, const int n, const int k) { - if (pack.GetIndex(b, Conserved::phi()) < 0) return; + if (n > pack.GetUpperBound(b)) return; typename Kokkos::MinMax::value_type minmax; par_reduce_inner( parthenon::inner_loop_pattern_ttr_tag, team_member, jb.s, jb.e, ib.s, ib.e, From cc97120166a2cccb94f0b0af096698ac8533f888 Mon Sep 17 00:00:00 2001 From: adam reyes Date: Mon, 25 Nov 2024 09:40:19 -0800 Subject: [PATCH 48/49] remove meshblock first/second derivative from cpp --- src/amr_criteria/refinement_package.cpp | 60 ------------------------- 1 file changed, 60 deletions(-) diff --git a/src/amr_criteria/refinement_package.cpp b/src/amr_criteria/refinement_package.cpp index 76d963b49be7..f4523358a692 100644 --- a/src/amr_criteria/refinement_package.cpp +++ b/src/amr_criteria/refinement_package.cpp @@ -97,66 +97,6 @@ AmrTag CheckAllRefinement(MeshBlockData *rc, const AmrTag &level) { return delta_level; } -AmrTag FirstDerivative(const AMRBounds &bnds, const ParArray3D &q, - const Real refine_criteria, const Real derefine_criteria) { - PARTHENON_INSTRUMENT - const int ndim = 1 + (bnds.je > bnds.js) + (bnds.ke > bnds.ks); - Real maxd = 0.0; - par_reduce( - loop_pattern_mdrange_tag, PARTHENON_AUTO_LABEL, DevExecSpace(), bnds.ks, bnds.ke, - bnds.js, bnds.je, bnds.is, bnds.ie, - KOKKOS_LAMBDA(int k, int j, int i, Real &maxd) { - Real scale = std::abs(q(k, j, i)); - Real d = - 0.5 * std::abs((q(k, j, i + 1) - q(k, j, i - 1))) / (scale + TINY_NUMBER); - maxd = (d > maxd ? d : maxd); - if (ndim > 1) { - d = 0.5 * std::abs((q(k, j + 1, i) - q(k, j - 1, i))) / (scale + TINY_NUMBER); - maxd = (d > maxd ? d : maxd); - } - if (ndim > 2) { - d = 0.5 * std::abs((q(k + 1, j, i) - q(k - 1, j, i))) / (scale + TINY_NUMBER); - maxd = (d > maxd ? d : maxd); - } - }, - Kokkos::Max(maxd)); - - if (maxd > refine_criteria) return AmrTag::refine; - if (maxd < derefine_criteria) return AmrTag::derefine; - return AmrTag::same; -} - -AmrTag SecondDerivative(const AMRBounds &bnds, const ParArray3D &q, - const Real refine_criteria, const Real derefine_criteria) { - PARTHENON_INSTRUMENT - const int ndim = 1 + (bnds.je > bnds.js) + (bnds.ke > bnds.ks); - Real maxd = 0.0; - par_reduce( - loop_pattern_mdrange_tag, PARTHENON_AUTO_LABEL, DevExecSpace(), bnds.ks, bnds.ke, - bnds.js, bnds.je, bnds.is, bnds.ie, - KOKKOS_LAMBDA(int k, int j, int i, Real &maxd) { - Real aqt = std::abs(q(k, j, i)) + TINY_NUMBER; - Real qavg = 0.5 * (q(k, j, i + 1) + q(k, j, i - 1)); - Real d = std::abs(qavg - q(k, j, i)) / (std::abs(qavg) + aqt); - maxd = (d > maxd ? d : maxd); - if (ndim > 1) { - qavg = 0.5 * (q(k, j + 1, i) + q(k, j - 1, i)); - d = std::abs(qavg - q(k, j, i)) / (std::abs(qavg) + aqt); - maxd = (d > maxd ? d : maxd); - } - if (ndim > 2) { - qavg = 0.5 * (q(k + 1, j, i) + q(k - 1, j, i)); - d = std::abs(qavg - q(k, j, i)) / (std::abs(qavg) + aqt); - maxd = (d > maxd ? d : maxd); - } - }, - Kokkos::Max(maxd)); - - if (maxd > refine_criteria) return AmrTag::refine; - if (maxd < derefine_criteria) return AmrTag::derefine; - return AmrTag::same; -} - void FirstDerivative(const AMRBounds &bnds, MeshData *md, const std::string &field, const int &idx, ParArray1D &amr_tags, const Real refine_criteria_, const Real derefine_criteria_, From 9338c54fe8cf19b0f04365e47ab77574ba1c8728 Mon Sep 17 00:00:00 2001 From: adam reyes Date: Mon, 25 Nov 2024 09:46:45 -0800 Subject: [PATCH 49/49] linting --- src/mesh/mesh.hpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/mesh/mesh.hpp b/src/mesh/mesh.hpp index 86b46f461e52..4a3bba6da38d 100644 --- a/src/mesh/mesh.hpp +++ b/src/mesh/mesh.hpp @@ -34,8 +34,6 @@ #include #include -#include "Kokkos_Core.hpp" -#include "basic_types.hpp" #include "bvals/boundary_conditions.hpp" #include "bvals/comms/tag_map.hpp" #include "config.hpp"