Add a MeshData variant for refinement tagging

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/files) Bugfix to particle defragmentation
 - [[PR 1184]](https://github.com/parthenon-hpc-lab/parthenon/pull/1184) Fix swarm block neighbor indexing in 1D, 2D

benchmarks/burgers/burgers_driver.cpp

@@ -120,27 +120,14 @@ TaskCollection BurgersDriver::MakeTaskCollection(BlockList_t &blocks, const int
 
     auto fill_deriv = tl.AddTask(update, FillDerived<MeshData<Real>>, 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<MeshData<Real>>, 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<MeshBlockData<Real>>, sc1.get());
+        auto tag_refine =
+            tl.AddTask(set_bc, parthenon::Refinement::Tag<MeshData<Real>>, mc1.get());
       }
     }
   }

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<Real>* 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<Real> *rc)``
   deletgates to the ``std::function`` member ``PreStepDiagnosticsMesh`` if
   set (defaults to ``nullptr`` an therefore a no-op) to print diagnostics

example/fine_advection/advection_package.cpp

@@ -118,50 +118,57 @@ std::shared_ptr<StateDescriptor> 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<Real> *rc) {
+void CheckRefinementMesh(MeshData<Real> *md, parthenon::ParArray1D<AmrTag> &amr_tags) {
   std::shared_ptr<StateDescriptor> pkg =
-      rc->GetMeshPointer()->packages.Get("advection_package");
+      md->GetMeshPointer()->packages.Get("advection_package");
   auto do_regular_advection = pkg->Param<bool>("do_regular_advection");
   if (do_regular_advection) {
     // refine on advected, for example.  could also be a derived quantity
-    static auto desc = parthenon::MakePackDescriptor<Conserved::phi>(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<Conserved::phi>(md);
+    auto pack = desc.GetPack(md);
 
-    typename Kokkos::MinMax<Real>::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<Real>::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<Real>(minmax));
-
-    auto pkg = pmb->packages.Get("advection_package");
     const auto &refine_tol = pkg->Param<Real>("refine_tol");
     const auto &derefine_tol = pkg->Param<Real>("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<Kokkos::Experimental::ScatterMax>();
+    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 (pack.GetIndex(b, Conserved::phi()) < 0) return;
+          typename Kokkos::MinMax<Real>::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<Real>::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<Real>(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<Real> *md) {

example/fine_advection/advection_package.hpp

@@ -49,6 +49,7 @@ VARIABLE(advection, divD);
 
 std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin);
 AmrTag CheckRefinement(MeshBlockData<Real> *rc);
+void CheckRefinementMesh(MeshData<Real> *md, parthenon::ParArray1D<AmrTag> &amr_tags);
 Real EstimateTimestep(MeshData<Real> *md);
 TaskStatus FillDerived(MeshData<Real> *md);
 

src/amr_criteria/amr_criteria.cpp

@@ -78,31 +78,44 @@ std::shared_ptr<AMRCriteria> AMRCriteria::MakeAMRCriteria(std::string &criteria,
                               block_name + ": " + criteria);
 }
 
-AMRBounds AMRCriteria::GetBounds(const MeshBlockData<Real> *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<Real> *rc) const {
-  if (!rc->HasVariable(field) || !rc->IsAllocated(field)) {
-    return AmrTag::same;
+void AMRFirstDerivative::operator()(MeshData<Real> *md,
+                                    ParArray1D<AmrTag> &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<Real> *rc) const {
-  if (!rc->HasVariable(field) || !rc->IsAllocated(field)) {
-    return AmrTag::same;
+void AMRSecondDerivative::operator()(MeshData<Real> *md,
+                                     ParArray1D<AmrTag> &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

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<Real> *rc) const = 0;
+  virtual void operator()(MeshData<Real> *md, ParArray1D<AmrTag> &delta_level) const = 0;
   std::string field;
   Real refine_criteria, derefine_criteria;
   int max_level;
   int comp6, comp5, comp4;
   static std::shared_ptr<AMRCriteria>
   MakeAMRCriteria(std::string &criteria, ParameterInput *pin, std::string &block_name);
-  AMRBounds GetBounds(const MeshBlockData<Real> *rc) const;
 };
 
 struct AMRFirstDerivative : public AMRCriteria {
   AMRFirstDerivative(ParameterInput *pin, std::string &block_name)
       : AMRCriteria(pin, block_name) {}
-  AmrTag operator()(const MeshBlockData<Real> *rc) const override;
+  void operator()(MeshData<Real> *md, ParArray1D<AmrTag> &delta_level) const override;
 };
 
 struct AMRSecondDerivative : public AMRCriteria {
   AMRSecondDerivative(ParameterInput *pin, std::string &block_name)
       : AMRCriteria(pin, block_name) {}
-  AmrTag operator()(const MeshBlockData<Real> *rc) const override;
+  void operator()(MeshData<Real> *md, ParArray1D<AmrTag> &delta_level) const override;
 };
 
 } // namespace parthenon