Merge branch 'mfournier01/cubicAMR' into mfournier01/sloshing

Merge cubicAMR into sloshing setup

Author: fournier2

src/hydro/hydro.cpp

@@ -816,6 +816,13 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
     pkg->AddParam<Real>("refinement/maxdensity_refine_above", refine_above);
   } else if (refine_str == "user") {
     pkg->CheckRefinementBlock = Hydro::ProblemCheckRefinementBlock;
+  } else if (refine_str == "cubic") {
+    pkg->CheckRefinementBlock = refinement::other::Cubic;
+    const auto active = pin->GetOrAddBoolean("refinement", "active", false);
+    const auto refinement_width =
+        pin->GetOrAddReal("refinement", "refinement_width", 0.0);
+    pkg->AddParam<>("refinement/active", active);
+    pkg->AddParam<>("refinement/refinement_width", refinement_width);
   }
 
   if (ProblemInitPackageData != nullptr) {

src/refinement/other.cpp

@@ -40,5 +40,60 @@ parthenon::AmrTag MaxDensity(MeshBlockData<Real> *rc) {
   return parthenon::AmrTag::same;
 }
 
+// refinement condition: cubic refinement with activation check
+parthenon::AmrTag Cubic(MeshBlockData<Real> *rc) {
+
+  auto pmb = rc->GetBlockPointer();
+  auto &coords = pmb->coords;
+
+  // Check if refinement is active
+  const bool active = pmb->packages.Get("Hydro")->Param<bool>("refinement/active");
+
+  if (!active) {
+    return parthenon::AmrTag::same; // Skip refinement if inactive
+  }
+
+  const Real refinement_width =
+      pmb->packages.Get("Hydro")->Param<Real>("refinement/refinement_width");
+  const Real half_width = refinement_width / 2.0;
+
+  // Retrieve bounds
+  IndexRange ib = pmb->cellbounds.GetBoundsI(IndexDomain::interior);
+  IndexRange jb = pmb->cellbounds.GetBoundsJ(IndexDomain::interior);
+  IndexRange kb = pmb->cellbounds.GetBoundsK(IndexDomain::interior);
+
+  // Fast bounding box check (block-level)
+  const Real x_min = coords.Xc<1>(ib.s);
+  const Real x_max = coords.Xc<1>(ib.e);
+  const Real y_min = coords.Xc<2>(jb.s);
+  const Real y_max = coords.Xc<2>(jb.e);
+  const Real z_min = coords.Xc<3>(kb.s);
+  const Real z_max = coords.Xc<3>(kb.e);
+
+  if (x_min > half_width || x_max < -half_width || y_min > half_width ||
+      y_max < -half_width || z_min > half_width || z_max < -half_width) {
+    return parthenon::AmrTag::same; // Fully outside, no refinement needed
+  }
+
+  // If the block intersects the cubic region, perform detailed check
+  bool inside_cubic_region = false;
+
+  pmb->par_reduce(
+      "cubic check refinement", kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
+      KOKKOS_LAMBDA(const int k, const int j, const int i, bool &inside) {
+        const Real x = coords.Xc<1>(i);
+        const Real y = coords.Xc<2>(j);
+        const Real z = coords.Xc<3>(k);
+
+        if (std::abs(x) <= half_width && std::abs(y) <= half_width &&
+            std::abs(z) <= half_width) {
+          inside = true;
+        }
+      },
+      Kokkos::LOr<bool>(inside_cubic_region));
+
+  return inside_cubic_region ? parthenon::AmrTag::refine : parthenon::AmrTag::same;
+}
+
 } // namespace other
 } // namespace refinement

src/refinement/refinement.hpp

@@ -20,7 +20,8 @@ AmrTag VelocityGradient(MeshBlockData<Real> *rc);
 } // namespace gradient
 namespace other {
 parthenon::AmrTag MaxDensity(MeshBlockData<Real> *rc);
-}
+parthenon::AmrTag Cubic(MeshBlockData<Real> *rc);
+} // namespace other
 
 } // namespace refinement