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 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

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<MeshData<Real>>, mc1.get());
+      auto tag_refine =
+          tl.AddTask(update, parthenon::Refinement::Tag<MeshData<Real>>, 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<MeshBlockData<Real>>, sc1.get());
-      }
-    }
   }
   return tc;
 }

example/fine_advection/advection_package.cpp

@@ -94,12 +94,59 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
   pkg->AddField<Conserved::divD>(
       Metadata({Metadata::Cell, Metadata::Derived, Metadata::OneCopy}));
 
-  pkg->CheckRefinementBlock = CheckRefinement;
+  bool check_refine_mesh =
+      pin->GetOrAddBoolean("parthenon/mesh", "CheckRefineMesh", true);
+  if (check_refine_mesh) {
+    pkg->CheckRefinementMesh = CheckRefinementMesh;
+  } else {
+    pkg->CheckRefinementBlock = CheckRefinement;
+  }
   pkg->EstimateTimestepMesh = EstimateTimestep;
   pkg->FillDerivedMesh = FillDerived;
   return pkg;
 }
 
+void CheckRefinementMesh(MeshData<Real> *md,
+                         parthenon::ParArray1D<AmrTag> &delta_levels) {
+  // refine on advected, for example.  could also be a derived quantity
+  static auto desc = parthenon::MakePackDescriptor<Conserved::phi>(md);
+  auto pack = desc.GetPack(md);
+
+  auto pkg = md->GetMeshPointer()->packages.Get("advection_package");
+  const auto &refine_tol = pkg->Param<Real>("refine_tol");
+  const auto &derefine_tol = pkg->Param<Real>("derefine_tol");
+
+  auto ib = md->GetBoundsI(IndexDomain::entire);
+  auto jb = md->GetBoundsJ(IndexDomain::entire);
+  auto kb = md->GetBoundsK(IndexDomain::entire);
+  auto scatter_levels = delta_levels.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) {
+        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(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<Real>(minmax));
+
+        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;
+        levels_access(b).update(flag);
+      });
+  delta_levels.ContributeScatter(scatter_levels);
+}
+
 AmrTag CheckRefinement(MeshBlockData<Real> *rc) {
   // refine on advected, for example.  could also be a derived quantity
   static auto desc = parthenon::MakePackDescriptor<Conserved::phi>(rc);

example/fine_advection/advection_package.hpp

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

src/amr_criteria/amr_criteria.cpp

@@ -95,6 +95,26 @@ AmrTag AMRFirstDerivative::operator()(const MeshBlockData<Real> *rc) const {
   return Refinement::FirstDerivative(bnds, q, refine_criteria, derefine_criteria);
 }
 
+void AMRFirstDerivative::operator()(MeshData<Real> *mc,
+                                    ParArray1D<AmrTag> &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);
+  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, max_level);
+}
+
 AmrTag AMRSecondDerivative::operator()(const MeshBlockData<Real> *rc) const {
   if (!rc->HasVariable(field) || !rc->IsAllocated(field)) {
     return AmrTag::same;
@@ -105,4 +125,24 @@ AmrTag AMRSecondDerivative::operator()(const MeshBlockData<Real> *rc) const {
   return Refinement::SecondDerivative(bnds, q, refine_criteria, derefine_criteria);
 }
 
+void AMRSecondDerivative::operator()(MeshData<Real> *mc,
+                                     ParArray1D<AmrTag> &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);
+  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, 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 {
 
@@ -36,6 +37,7 @@ 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> *mc, ParArray1D<AmrTag> &delta_level) const = 0;
   std::string field;
   Real refine_criteria, derefine_criteria;
   int max_level;
@@ -49,12 +51,14 @@ 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> *mc, 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> *mc, ParArray1D<AmrTag> &delta_level) const override;
 };
 
 } // namespace parthenon

src/amr_criteria/refinement_package.cpp

@@ -20,19 +20,25 @@
 #include <utility>
 
 #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 {
 
 std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
   auto ref = std::make_shared<StateDescriptor>("Refinement");
+  bool check_refine_mesh =
+      pin->GetOrAddBoolean("parthenon/mesh", "CheckRefineMesh", true);
+  ref->AddParam("check_refine_mesh", check_refine_mesh);
 
   int numcrit = 0;
   while (true) {
@@ -48,7 +54,32 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
   return ref;
 }
 
-AmrTag CheckAllRefinement(MeshBlockData<Real> *rc) {
+ParArray1D<AmrTag> CheckAllRefinement(MeshData<Real> *mc) {
+  const int nblocks = mc->NumBlocks();
+  // maybe not great to allocate this all the time
+  auto delta_levels = ParArray1D<AmrTag>(Kokkos::View<AmrTag *>(
+      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<bool>("check_refine_mesh");
+
+  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, delta_levels);
+      }
+    }
+  }
+
+  return delta_levels;
+}
+
+AmrTag CheckAllRefinement(MeshBlockData<Real> *rc, const AmrTag &level) {
   // Check all refinement criteria and return the maximum recommended change in
   // refinement level:
   //   delta_level = -1 => recommend derefinement
@@ -62,15 +93,19 @@ AmrTag CheckAllRefinement(MeshBlockData<Real> *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;
+  static const bool check_refine_mesh =
+      pmb->packages.Get("Refinement")->Param<bool>("check_refine_mesh");
+  // 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));
     if (delta_level == AmrTag::refine) {
       // 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 +185,111 @@ AmrTag SecondDerivative(const AMRBounds &bnds, const ParArray3D<Real> &q,
   return AmrTag::same;
 }
 
-void SetRefinement_(MeshBlockData<Real> *rc) {
+void FirstDerivative(const AMRBounds &bnds, MeshData<Real> *mc, const std::string &field,
+                     const int &idx, ParArray1D<AmrTag> &delta_levels,
+                     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;
+  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<Kokkos::Experimental::ScatterMax>();
+  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<Real>(maxd));
+        auto levels_access = scatter_levels.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);
+      });
+  delta_levels.ContributeScatter(scatter_levels);
+}
+
+void SecondDerivative(const AMRBounds &bnds, MeshData<Real> *mc, const std::string &field,
+                      const int &idx, ParArray1D<AmrTag> &delta_levels,
+                      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;
+  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<Kokkos::Experimental::ScatterMax>();
+  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<Real>(maxd));
+        auto levels_access = scatter_levels.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);
+      });
+  delta_levels.ContributeScatter(scatter_levels);
+}
+
+void SetRefinement_(MeshBlockData<Real> *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 +302,11 @@ TaskStatus Tag(MeshBlockData<Real> *rc) {
 template <>
 TaskStatus Tag(MeshData<Real> *rc) {
   PARTHENON_INSTRUMENT
+  ParArray1D<AmrTag> 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;
 }