Boundary bugfix(es) and exposing FluxDiv_ interface

example/advection/parthenon_app_inputs.cpp

@@ -47,9 +47,9 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) {
   const auto &profile = pkg->Param<std::string>("profile");
 
   auto cellbounds = pmb->cellbounds;
-  IndexRange ib = cellbounds.GetBoundsI(IndexDomain::entire);
-  IndexRange jb = cellbounds.GetBoundsJ(IndexDomain::entire);
-  IndexRange kb = cellbounds.GetBoundsK(IndexDomain::entire);
+  IndexRange ib = cellbounds.GetBoundsI(IndexDomain::interior);
+  IndexRange jb = cellbounds.GetBoundsJ(IndexDomain::interior);
+  IndexRange kb = cellbounds.GetBoundsK(IndexDomain::interior);
 
   auto coords = pmb->coords;
   auto q = data->PackVariables(std::vector<MetadataFlag>{Metadata::Independent});

src/bvals/boundary_conditions.cpp

@@ -135,7 +135,7 @@ void OutflowOuterX2(std::shared_ptr<MeshBlockData<Real>> &rc, bool coarse) {
 void OutflowInnerX3(std::shared_ptr<MeshBlockData<Real>> &rc, bool coarse) {
   std::shared_ptr<MeshBlock> pmb = rc->GetBlockPointer();
   auto bounds = coarse ? pmb->c_cellbounds : pmb->cellbounds;
-  int ref = bounds.GetBoundsJ(IndexDomain::interior).s;
+  int ref = bounds.GetBoundsK(IndexDomain::interior).s;
   auto q = rc->PackVariables(std::vector<MetadataFlag>{Metadata::FillGhost}, coarse);
   auto nb = IndexRange{0, q.GetDim(4) - 1};
   pmb->par_for_bndry(
@@ -148,7 +148,7 @@ void OutflowInnerX3(std::shared_ptr<MeshBlockData<Real>> &rc, bool coarse) {
 void OutflowOuterX3(std::shared_ptr<MeshBlockData<Real>> &rc, bool coarse) {
   std::shared_ptr<MeshBlock> pmb = rc->GetBlockPointer();
   auto bounds = coarse ? pmb->c_cellbounds : pmb->cellbounds;
-  int ref = bounds.GetBoundsJ(IndexDomain::interior).e;
+  int ref = bounds.GetBoundsK(IndexDomain::interior).e;
   auto q = rc->PackVariables(std::vector<MetadataFlag>{Metadata::FillGhost}, coarse);
   auto nb = IndexRange{0, q.GetDim(4) - 1};
   pmb->par_for_bndry(

src/driver/driver.cpp

@@ -156,7 +156,9 @@ void EvolutionDriver::InitializeBlockTimeSteps() {
 void EvolutionDriver::SetGlobalTimeStep() {
   // don't allow dt to grow by more than 2x
   // consider making this configurable in the input
-  tm.dt *= 2.0;
+  if (tm.dt < 0.1 * std::numeric_limits<Real>::max()) {
+    tm.dt *= 2.0;
+  }
   Real big = std::numeric_limits<Real>::max();
   for (auto const &pmb : pmesh->block_list) {
     tm.dt = std::min(tm.dt, pmb->NewDt());

src/interface/update.cpp

@@ -29,22 +29,6 @@ namespace parthenon {
 
 namespace Update {
 
-KOKKOS_FORCEINLINE_FUNCTION
-Real FluxDiv_(const int l, const int k, const int j, const int i, const int ndim,
-              const Coordinates_t &coords, const VariableFluxPack<Real> &v) {
-  Real du = (coords.Area(X1DIR, k, j, i + 1) * v.flux(X1DIR, l, k, j, i + 1) -
-             coords.Area(X1DIR, k, j, i) * v.flux(X1DIR, l, k, j, i));
-  if (ndim >= 2) {
-    du += (coords.Area(X2DIR, k, j + 1, i) * v.flux(X2DIR, l, k, j + 1, i) -
-           coords.Area(X2DIR, k, j, i) * v.flux(X2DIR, l, k, j, i));
-  }
-  if (ndim == 3) {
-    du += (coords.Area(X3DIR, k + 1, j, i) * v.flux(X3DIR, l, k + 1, j, i) -
-           coords.Area(X3DIR, k, j, i) * v.flux(X3DIR, l, k, j, i));
-  }
-  return -du / coords.Volume(k, j, i);
-}
-
 template <>
 TaskStatus FluxDivergence(MeshBlockData<Real> *in, MeshBlockData<Real> *dudt_cont) {
   std::shared_ptr<MeshBlock> pmb = in->GetBlockPointer();
@@ -62,7 +46,7 @@ TaskStatus FluxDivergence(MeshBlockData<Real> *in, MeshBlockData<Real> *dudt_con
   pmb->par_for(
       "FluxDivergenceBlock", 0, vin.GetDim(4) - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
       KOKKOS_LAMBDA(const int l, const int k, const int j, const int i) {
-        dudt(l, k, j, i) = FluxDiv_(l, k, j, i, ndim, coords, vin);
+        dudt(l, k, j, i) = FluxDivHelper(l, k, j, i, ndim, coords, vin);
       });
 
   return TaskStatus::complete;
@@ -86,7 +70,7 @@ TaskStatus FluxDivergence(MeshData<Real> *in_obj, MeshData<Real> *dudt_obj) {
       KOKKOS_LAMBDA(const int m, const int l, const int k, const int j, const int i) {
         const auto &coords = vin.coords(m);
         const auto &v = vin(m);
-        dudt(m, l, k, j, i) = FluxDiv_(l, k, j, i, ndim, coords, v);
+        dudt(m, l, k, j, i) = FluxDivHelper(l, k, j, i, ndim, coords, v);
       });
   return TaskStatus::complete;
 }
@@ -111,7 +95,7 @@ TaskStatus UpdateWithFluxDivergence(MeshBlockData<Real> *u0_data,
       "UpdateWithFluxDivergenceBlock", 0, u0.GetDim(4) - 1, kb.s, kb.e, jb.s, jb.e, ib.s,
       ib.e, KOKKOS_LAMBDA(const int l, const int k, const int j, const int i) {
         u0(l, k, j, i) = gam0 * u0(l, k, j, i) + gam1 * u1(l, k, j, i) +
-                         beta_dt * FluxDiv_(l, k, j, i, ndim, coords, u0);
+                         beta_dt * FluxDivHelper(l, k, j, i, ndim, coords, u0);
       });
 
   return TaskStatus::complete;
@@ -138,7 +122,7 @@ TaskStatus UpdateWithFluxDivergence(MeshData<Real> *u0_data, MeshData<Real> *u1_
         const auto &coords = u0_pack.coords(m);
         const auto &u0 = u0_pack(m);
         u0_pack(m, l, k, j, i) = gam0 * u0(l, k, j, i) + gam1 * u1_pack(m, l, k, j, i) +
-                                 beta_dt * FluxDiv_(l, k, j, i, ndim, coords, u0);
+                                 beta_dt * FluxDivHelper(l, k, j, i, ndim, coords, u0);
       });
   return TaskStatus::complete;
 }

src/interface/update.hpp

@@ -34,6 +34,23 @@ namespace parthenon {
 
 namespace Update {
 
+// Calculate the flux divergence for a specific component l of a variable v
+KOKKOS_FORCEINLINE_FUNCTION
+Real FluxDivHelper(const int l, const int k, const int j, const int i, const int ndim,
+                   const Coordinates_t &coords, const VariableFluxPack<Real> &v) {
+  Real du = (coords.Area(X1DIR, k, j, i + 1) * v.flux(X1DIR, l, k, j, i + 1) -
+             coords.Area(X1DIR, k, j, i) * v.flux(X1DIR, l, k, j, i));
+  if (ndim >= 2) {
+    du += (coords.Area(X2DIR, k, j + 1, i) * v.flux(X2DIR, l, k, j + 1, i) -
+           coords.Area(X2DIR, k, j, i) * v.flux(X2DIR, l, k, j, i));
+  }
+  if (ndim == 3) {
+    du += (coords.Area(X3DIR, k + 1, j, i) * v.flux(X3DIR, l, k + 1, j, i) -
+           coords.Area(X3DIR, k, j, i) * v.flux(X3DIR, l, k, j, i));
+  }
+  return -du / coords.Volume(k, j, i);
+}
+
 template <typename T>
 TaskStatus FluxDivergence(T *in, T *dudt_obj);
 

src/mesh/domain.hpp

@@ -177,19 +177,14 @@ class IndexShape {
     }
   }
 
-  // note ks and ke are different than the others for boundaries.
-  // This is because extending i and j over extra ghost zones
-  // is sufficient to fill whole ghost layer
   KOKKOS_INLINE_FUNCTION int ks(const IndexDomain &domain) const noexcept {
     switch (domain) {
-    case IndexDomain::entire:
-      return 0;
-    case IndexDomain::inner_x3:
-      return 0;
+    case IndexDomain::interior:
+      return x_[2].s;
     case IndexDomain::outer_x3:
       return entire_ncells_[2] == 1 ? 0 : x_[2].e + 1;
-    default: // interior+
-      return x_[2].s;
+    default:
+      return 0;
     }
   }
 
@@ -217,14 +212,12 @@ class IndexShape {
 
   KOKKOS_INLINE_FUNCTION int ke(const IndexDomain &domain) const noexcept {
     switch (domain) {
-    case IndexDomain::entire:
-      return entire_ncells_[2] - 1;
+    case IndexDomain::interior:
+      return x_[2].e;
     case IndexDomain::inner_x3:
       return x_[2].s == 0 ? 0 : x_[2].s - 1;
-    case IndexDomain::outer_x3:
+    default:
       return entire_ncells_[2] - 1;
-    default: // interior+
-      return x_[2].e;
     }
   }
 

tst/unit/test_unit_domain.cpp

@@ -199,33 +199,33 @@ TEST_CASE("Checking IndexShape indices", "[IndexShape]") {
     REQUIRE(shape.ie(inner_x1) == 0);
     REQUIRE(shape.js(inner_x1) == 0);
     REQUIRE(shape.je(inner_x1) == 2);
-    REQUIRE(shape.ks(inner_x1) == 1);
-    REQUIRE(shape.ke(inner_x1) == 4);
-    REQUIRE(shape.GetTotal(inner_x1) == 12);
+    REQUIRE(shape.ks(inner_x1) == 0);
+    REQUIRE(shape.ke(inner_x1) == 5);
+    REQUIRE(shape.GetTotal(inner_x1) == 18);
 
     REQUIRE(shape.is(outer_x1) == 7);
     REQUIRE(shape.ie(outer_x1) == 7);
     REQUIRE(shape.js(outer_x1) == 0);
     REQUIRE(shape.je(outer_x1) == 2);
-    REQUIRE(shape.ks(outer_x1) == 1);
-    REQUIRE(shape.ke(outer_x1) == 4);
-    REQUIRE(shape.GetTotal(outer_x1) == 12);
+    REQUIRE(shape.ks(outer_x1) == 0);
+    REQUIRE(shape.ke(outer_x1) == 5);
+    REQUIRE(shape.GetTotal(outer_x1) == 18);
 
     REQUIRE(shape.is(inner_x2) == 0);
     REQUIRE(shape.ie(inner_x2) == 7);
     REQUIRE(shape.js(inner_x2) == 0);
     REQUIRE(shape.je(inner_x2) == 0);
-    REQUIRE(shape.ks(inner_x2) == 1);
-    REQUIRE(shape.ke(inner_x2) == 4);
-    REQUIRE(shape.GetTotal(inner_x2) == 32);
+    REQUIRE(shape.ks(inner_x2) == 0);
+    REQUIRE(shape.ke(inner_x2) == 5);
+    REQUIRE(shape.GetTotal(inner_x2) == 48);
 
     REQUIRE(shape.is(outer_x2) == 0);
     REQUIRE(shape.ie(outer_x2) == 7);
     REQUIRE(shape.js(outer_x2) == 2);
     REQUIRE(shape.je(outer_x2) == 2);
-    REQUIRE(shape.ks(outer_x2) == 1);
-    REQUIRE(shape.ke(outer_x2) == 4);
-    REQUIRE(shape.GetTotal(outer_x2) == 32);
+    REQUIRE(shape.ks(outer_x2) == 0);
+    REQUIRE(shape.ke(outer_x2) == 5);
+    REQUIRE(shape.GetTotal(outer_x2) == 48);
 
     REQUIRE(shape.is(inner_x3) == 0);
     REQUIRE(shape.ie(inner_x3) == 7);
@@ -398,19 +398,19 @@ TEST_CASE("Checking IndexShape cell counts", "[IndexShape]") {
 
     REQUIRE(shape.ncellsi(inner_x1) == 1);
     REQUIRE(shape.ncellsj(inner_x1) == 3);
-    REQUIRE(shape.ncellsk(inner_x1) == 4);
+    REQUIRE(shape.ncellsk(inner_x1) == 6);
 
     REQUIRE(shape.ncellsi(outer_x1) == 1);
     REQUIRE(shape.ncellsj(outer_x1) == 3);
-    REQUIRE(shape.ncellsk(outer_x1) == 4);
+    REQUIRE(shape.ncellsk(outer_x1) == 6);
 
     REQUIRE(shape.ncellsi(inner_x2) == 8);
     REQUIRE(shape.ncellsj(inner_x2) == 1);
-    REQUIRE(shape.ncellsk(inner_x2) == 4);
+    REQUIRE(shape.ncellsk(inner_x2) == 6);
 
     REQUIRE(shape.ncellsi(outer_x2) == 8);
     REQUIRE(shape.ncellsj(outer_x2) == 1);
-    REQUIRE(shape.ncellsk(outer_x2) == 4);
+    REQUIRE(shape.ncellsk(outer_x2) == 6);
 
     REQUIRE(shape.ncellsi(inner_x3) == 8);
     REQUIRE(shape.ncellsj(inner_x3) == 3);