diff --git a/CHANGELOG.md b/CHANGELOG.md index 3abf18c52538..678ec7782046 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 1103]](https://github.com/parthenon-hpc-lab/parthenon/pull/1103) Add sparsity to vector wave equation test - [[PR 1185]](https://github.com/parthenon-hpc-lab/parthenon/pull/1185) Bugfix to particle defragmentation - [[PR 1184]](https://github.com/parthenon-hpc-lab/parthenon/pull/1184) Fix swarm block neighbor indexing in 1D, 2D diff --git a/benchmarks/burgers/burgers_driver.cpp b/benchmarks/burgers/burgers_driver.cpp index cddb255a636b..abb2a64a03d1 100644 --- a/benchmarks/burgers/burgers_driver.cpp +++ b/benchmarks/burgers/burgers_driver.cpp @@ -120,27 +120,14 @@ 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()); - } - } - - 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); - - if (stage == integrator->nstages) { - // Update refinement if (pmesh->adaptive) { - auto tag_refine = tl.AddTask( - set_bc, parthenon::Refinement::Tag>, sc1.get()); + auto tag_refine = + tl.AddTask(set_bc, parthenon::Refinement::Tag>, mc1.get()); } } } 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 diff --git a/example/fine_advection/advection_package.cpp b/example/fine_advection/advection_package.cpp index 428d91b9245f..a134d3202ef1 100644 --- a/example/fine_advection/advection_package.cpp +++ b/example/fine_advection/advection_package.cpp @@ -118,50 +118,57 @@ std::shared_ptr Initialize(ParameterInput *pin) { Metadata({Metadata::Cell, Metadata::Derived, Metadata::OneCopy})); } - pkg->CheckRefinementBlock = CheckRefinement; + pkg->CheckRefinementMesh = CheckRefinementMesh; pkg->EstimateTimestepMesh = EstimateTimestep; pkg->FillDerivedMesh = FillDerived; return pkg; } -AmrTag CheckRefinement(MeshBlockData *rc) { +void CheckRefinementMesh(MeshData *md, parthenon::ParArray1D &amr_tags) { std::shared_ptr pkg = - rc->GetMeshPointer()->packages.Get("advection_package"); + md->GetMeshPointer()->packages.Get("advection_package"); auto do_regular_advection = pkg->Param("do_regular_advection"); if (do_regular_advection) { // 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); + static auto desc = parthenon::MakePackDescriptor(md); + auto pack = desc.GetPack(md); - typename Kokkos::MinMax::value_type minmax; - parthenon::par_reduce( - parthenon::loop_pattern_mdrange_tag, PARTHENON_AUTO_LABEL, DevExecSpace(), 0, - pack.GetNBlocks() - 1, // Runs from [0, 0] since pack built from MeshBlockData - pack.GetLowerBoundHost(0), pack.GetUpperBoundHost(0), 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, - typename Kokkos::MinMax::value_type &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(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; + auto ib = md->GetBoundsI(IndexDomain::entire); + auto jb = md->GetBoundsJ(IndexDomain::entire); + auto kb = md->GetBoundsK(IndexDomain::entire); + 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, + KOKKOS_LAMBDA(parthenon::team_mbr_t team_member, const int b, const int n, + const int k) { + 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, + [&](const int j, const int i, + typename Kokkos::MinMax::value_type &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(minmax)); + + 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; + tags_access(b).update(flag); + }); + amr_tags.ContributeScatter(scatter_tags); } - return AmrTag::same; } Real EstimateTimestep(MeshData *md) { diff --git a/example/fine_advection/advection_package.hpp b/example/fine_advection/advection_package.hpp index 12b6b9b65aad..cd1242b01928 100644 --- a/example/fine_advection/advection_package.hpp +++ b/example/fine_advection/advection_package.hpp @@ -49,6 +49,7 @@ VARIABLE(advection, divD); std::shared_ptr Initialize(ParameterInput *pin); AmrTag CheckRefinement(MeshBlockData *rc); +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 41b2237d96ff..a0625b81e652 100644 --- a/src/amr_criteria/amr_criteria.cpp +++ b/src/amr_criteria/amr_criteria.cpp @@ -78,31 +78,44 @@ 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; +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 = md->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]; } - 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); + const int idx = comp4 + n4 * (comp5 + n5 * comp6); + Refinement::FirstDerivative(bnds, md, field, idx, amr_tags, refine_criteria, + derefine_criteria, max_level); } -AmrTag AMRSecondDerivative::operator()(const MeshBlockData *rc) const { - if (!rc->HasVariable(field) || !rc->IsAllocated(field)) { - return AmrTag::same; +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 = md->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]; } - 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); + const int idx = comp4 + n4 * (comp5 + n5 * comp6); + Refinement::SecondDerivative(bnds, md, field, idx, amr_tags, refine_criteria, + derefine_criteria, max_level); } } // namespace parthenon diff --git a/src/amr_criteria/amr_criteria.hpp b/src/amr_criteria/amr_criteria.hpp index e31bd8b381d8..e69aebe84733 100644 --- a/src/amr_criteria/amr_criteria.hpp +++ b/src/amr_criteria/amr_criteria.hpp @@ -18,6 +18,7 @@ #include "defs.hpp" #include "mesh/domain.hpp" +#include "mesh/mesh.hpp" namespace parthenon { @@ -35,26 +36,25 @@ 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; int max_level; 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; }; } // namespace parthenon diff --git a/src/amr_criteria/refinement_package.cpp b/src/amr_criteria/refinement_package.cpp index 81342b70695a..f4523358a692 100644 --- a/src/amr_criteria/refinement_package.cpp +++ b/src/amr_criteria/refinement_package.cpp @@ -20,13 +20,16 @@ #include #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" #include "parameter_input.hpp" +#include "utils/instrument.hpp" namespace parthenon { namespace Refinement { @@ -48,7 +51,25 @@ std::shared_ptr Initialize(ParameterInput *pin) { return ref; } -AmrTag CheckAllRefinement(MeshBlockData *rc) { +ParArray1D CheckAllRefinement(MeshData *md) { + const int nblocks = md->NumBlocks(); + Mesh *pm = md->GetMeshPointer(); + auto amr_tags = pm->GetAmrTags(); + Kokkos::deep_copy(amr_tags.KokkosView(), AmrTag::derefine); + + for (auto &pkg : pm->packages.AllPackages()) { + auto &desc = pkg.second; + desc->CheckRefinement(md, amr_tags); + + for (auto &amr : desc->amr_criteria) { + (*amr)(md, amr_tags); + } + } + + return amr_tags; +} + +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 +83,9 @@ AmrTag CheckAllRefinement(MeshBlockData *rc) { // neighboring blocks. Similarly for "do nothing" PARTHENON_INSTRUMENT MeshBlock *pmb = rc->GetBlockPointer(); - // delta_level holds the max over all criteria. default to derefining. - AmrTag delta_level = AmrTag::derefine; + // 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; delta_level = std::max(delta_level, desc->CheckRefinement(rc)); @@ -71,88 +93,121 @@ AmrTag CheckAllRefinement(MeshBlockData *rc) { // since 1 is the max, we can return without having to look at anything else return AmrTag::refine; } - // 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; } -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; +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(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; + // 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, + 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 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; + // 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); } -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 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(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; + // 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, + 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 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; + // 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); } -void SetRefinement_(MeshBlockData *rc) { +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 <> @@ -163,10 +218,13 @@ TaskStatus Tag(MeshBlockData *rc) { } template <> -TaskStatus Tag(MeshData *rc) { +TaskStatus Tag(MeshData *md) { PARTHENON_INSTRUMENT - for (int i = 0; i < rc->NumBlocks(); i++) { - SetRefinement_(rc->GetBlockData(i).get()); + ParArray1D amr_tags = CheckAllRefinement(md); + auto amr_tags_h = amr_tags.GetHostMirrorAndCopy(); + + 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 40e9cd39fb14..82439f947df1 100644 --- a/src/amr_criteria/refinement_package.hpp +++ b/src/amr_criteria/refinement_package.hpp @@ -37,13 +37,18 @@ std::shared_ptr Initialize(ParameterInput *pin); template TaskStatus Tag(T *rc); -AmrTag CheckAllRefinement(MeshBlockData *rc); - -AmrTag FirstDerivative(const AMRBounds &bnds, const ParArray3D &q, - const Real refine_criteria, const Real derefine_criteria); - -AmrTag SecondDerivative(const AMRBounds &bnds, const ParArray3D &q, - const Real refine_criteria, const Real derefine_criteria); +AmrTag CheckAllRefinement(MeshBlockData *rc, const AmrTag &level); +ParArray1D CheckAllRefinement(MeshData *md); + +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_); + +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_); } // namespace Refinement diff --git a/src/interface/state_descriptor.hpp b/src/interface/state_descriptor.hpp index 8f8898170f76..a516f6be8821 100644 --- a/src/interface/state_descriptor.hpp +++ b/src/interface/state_descriptor.hpp @@ -340,6 +340,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); } @@ -389,6 +393,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; diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp index 86d727bae46b..b4d48ae5bccd 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" @@ -814,8 +815,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 @@ -905,6 +907,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 ed26be1c6a1b..4a3bba6da38d 100644 --- a/src/mesh/mesh.hpp +++ b/src/mesh/mesh.hpp @@ -109,6 +109,8 @@ class Mesh { } const IndexShape GetLeafBlockCellBounds(CellLevel level = CellLevel::same) const; + ParArray1D &GetAmrTags(); + const forest::Forest &Forest() const { return forest; } // data @@ -281,6 +283,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; 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)); 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() diff --git a/src/parthenon_array_generic.hpp b/src/parthenon_array_generic.hpp index ac38b6cb5e5f..0840b46f8d11 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();