Bump Kokkos 4.4.1, Parth vov-fix

CHANGELOG.md

@@ -15,6 +15,7 @@
 ### Fixed (not changing behavior/API/variables/...)
 
 ### Infrastructure
+- [[PR 124]](https://github.com/parthenon-hpc-lab/athenapk/pull/124) Bump Kokkos 4.4.1 (and Parthenon to include view-of-view fix)
 - [[PR 117]](https://github.com/parthenon-hpc-lab/athenapk/pull/117) Update devcontainer.json to latest CI container
 - [[PR 114]](https://github.com/parthenon-hpc-lab/athenapk/pull/114) Bump Parthenon 24.08 and Kokkos to 4.4.00
 - [[PR 112]](https://github.com/parthenon-hpc-lab/athenapk/pull/112) Add dev container configuration
@@ -24,6 +25,8 @@
 ### Removed (removing behavior/API/varaibles/...)
 
 ### Incompatibilities (i.e. breaking changes)
+- [[PR 124]](https://github.com/parthenon-hpc-lab/athenapk/pull/124) Enrolling custom boundary conditions changed
+  - Boundary conditions can now be enrolled using a string that can be subsequently be used in the input file (see, e.g., cloud problem generator)
 - [[PR 114]](https://github.com/parthenon-hpc-lab/athenapk/pull/114) Bump Parthenon 24.08 and Kokkos to 4.4.00
   - Changed signature of `UserWorkBeforeOutput` to include `SimTime` as last paramter
   - Fixes bitwise idential restarts for AMR simulations (the derefinement counter is now included)

docs/input.md

@@ -118,9 +118,6 @@ However, as reported by [^V+17] a large number of stages
 (in combination with anisotropic, limited diffusion) may lead to a loss in accuracy, which
 is why the difference between hyperbolic and parabolic timesteps can be limited by
 `diffusion/rkl2_max_dt_ratio=...` and a warning is shown if the ratio is above 400.
-Note that if this limit is enforced the `dt=` shown on the terminal might not be the actual
-`dt` taken in the code as the limit is currently enforced only after the output
-has been printed.
 
 [^M+14]:
     C. D. Meyer, D. S. Balsara, and T. D. Aslam, “A stabilized Runge–Kutta–Legendre method for explicit super-time-stepping of parabolic and mixed equations,” Journal of Computational Physics, vol. 257, pp. 594–626, 2014, doi: https://doi.org/10.1016/j.jcp.2013.08.021.

inputs/cloud.in

@@ -22,7 +22,7 @@ ox1_bc = outflow
 nx2 = 320
 x2min = -400
 x2max = 2100
-ix2_bc = user
+ix2_bc = cloud_inflow_x2
 ox2_bc = outflow
 
 nx3 = 192

src/hydro/hydro.cpp

@@ -4,6 +4,7 @@
 // Licensed under the BSD 3-Clause License (the "LICENSE").
 //========================================================================================
 
+#include <algorithm>
 #include <limits>
 #include <memory>
 #include <string>
@@ -29,6 +30,7 @@
 #include "diffusion/diffusion.hpp"
 #include "glmmhd/glmmhd.hpp"
 #include "hydro.hpp"
+#include "interface/params.hpp"
 #include "outputs/outputs.hpp"
 #include "prolongation/custom_ops.hpp"
 #include "rsolvers/rsolvers.hpp"
@@ -60,44 +62,7 @@ parthenon::Packages_t ProcessPackages(std::unique_ptr<ParameterInput> &pin) {
 // the task list is constructed (versus when the task list is being executed).
 // TODO(next person touching this function): If more/separate feature are required
 // please separate concerns.
-void PreStepMeshUserWorkInLoop(Mesh *pmesh, ParameterInput *pin, SimTime &tm) {
-  auto hydro_pkg = pmesh->block_list[0]->packages.Get("Hydro");
-  const auto num_partitions = pmesh->DefaultNumPartitions();
-
-  if (hydro_pkg->Param<DiffInt>("diffint") == DiffInt::rkl2) {
-    auto dt_diff = std::numeric_limits<Real>::max();
-    if (hydro_pkg->Param<Conduction>("conduction") != Conduction::none) {
-      for (auto i = 0; i < num_partitions; i++) {
-        auto &md = pmesh->mesh_data.GetOrAdd("base", i);
-
-        dt_diff = std::min(dt_diff, EstimateConductionTimestep(md.get()));
-      }
-    }
-    if (hydro_pkg->Param<Viscosity>("viscosity") != Viscosity::none) {
-      for (auto i = 0; i < num_partitions; i++) {
-        auto &md = pmesh->mesh_data.GetOrAdd("base", i);
-
-        dt_diff = std::min(dt_diff, EstimateViscosityTimestep(md.get()));
-      }
-    }
-    if (hydro_pkg->Param<Resistivity>("resistivity") != Resistivity::none) {
-      for (auto i = 0; i < num_partitions; i++) {
-        auto &md = pmesh->mesh_data.GetOrAdd("base", i);
-
-        dt_diff = std::min(dt_diff, EstimateResistivityTimestep(md.get()));
-      }
-    }
-#ifdef MPI_PARALLEL
-    PARTHENON_MPI_CHECK(MPI_Allreduce(MPI_IN_PLACE, &dt_diff, 1, MPI_PARTHENON_REAL,
-                                      MPI_MIN, MPI_COMM_WORLD));
-#endif
-    hydro_pkg->UpdateParam("dt_diff", dt_diff);
-    const auto max_dt_ratio = hydro_pkg->Param<Real>("rkl2_max_dt_ratio");
-    if (max_dt_ratio > 0.0 && tm.dt / dt_diff > max_dt_ratio) {
-      tm.dt = max_dt_ratio * dt_diff;
-    }
-  }
-}
+void PreStepMeshUserWorkInLoop(Mesh *pmesh, ParameterInput *pin, SimTime &tm) {}
 
 template <Hst hst, int idx = -1>
 Real HydroHst(MeshData<Real> *md) {
@@ -252,17 +217,23 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
     auto glmmhd_alpha = pin->GetOrAddReal("hydro", "glmmhd_alpha", 0.1);
     pkg->AddParam<Real>("glmmhd_alpha", glmmhd_alpha);
     calc_c_h = true;
-    pkg->AddParam<Real>("c_h", 0.0, true); // hyperbolic divergence cleaning speed
-    // global minimum dx (used to calc c_h)
-    pkg->AddParam<Real>("mindx", std::numeric_limits<Real>::max(), true);
-    // hyperbolic timestep constraint
-    pkg->AddParam<Real>("dt_hyp", std::numeric_limits<Real>::max(), true);
   } else {
     PARTHENON_FAIL("AthenaPK hydro: Unknown fluid method.");
   }
   pkg->AddParam<>("fluid", fluid);
   pkg->AddParam<>("nhydro", nhydro);
   pkg->AddParam<>("calc_c_h", calc_c_h);
+  // Following params should (currently) be present independent of solver because
+  // they're all used in the main loop.
+  // TODO(pgrete) think about which approach (selective versus always is preferable)
+  pkg->AddParam<Real>(
+      "c_h", 0.0, Params::Mutability::Restart); // hyperbolic divergence cleaning speed
+  // global minimum dx (used to calc c_h)
+  pkg->AddParam<Real>("mindx", std::numeric_limits<Real>::max(),
+                      Params::Mutability::Restart);
+  // hyperbolic timestep constraint
+  pkg->AddParam<Real>("dt_hyp", std::numeric_limits<Real>::max(),
+                      Params::Mutability::Restart);
 
   const auto recon_str = pin->GetString("hydro", "reconstruction");
   int recon_need_nghost = 3; // largest number for the choices below
@@ -633,12 +604,13 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
                      "Options are: none, unsplit, rkl2");
     }
     if (diffint != DiffInt::none) {
-      pkg->AddParam<Real>("dt_diff", 0.0, true); // diffusive timestep constraint
       // As in Athena++ a cfl safety factor is also applied to the theoretical limit.
       // By default it is equal to the hyperbolic cfl.
       auto cfl_diff = pin->GetOrAddReal("diffusion", "cfl", pkg->Param<Real>("cfl"));
       pkg->AddParam<>("cfl_diff", cfl_diff);
     }
+    pkg->AddParam<Real>("dt_diff", std::numeric_limits<Real>::max(),
+                        Params::Mutability::Restart); // diffusive timestep constraint
     pkg->AddParam<>("diffint", diffint);
 
     if (fluid == Fluid::euler) {
@@ -855,9 +827,12 @@ Real EstimateTimestep(MeshData<Real> *md) {
   // get to package via first block in Meshdata (which exists by construction)
   auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro");
   auto min_dt = std::numeric_limits<Real>::max();
+  auto dt_hyp = std::numeric_limits<Real>::max();
 
-  if (hydro_pkg->Param<bool>("calc_dt_hyp")) {
-    min_dt = std::min(min_dt, EstimateHyperbolicTimestep<fluid>(md));
+  const auto calc_dt_hyp = hydro_pkg->Param<bool>("calc_dt_hyp");
+  if (calc_dt_hyp) {
+    dt_hyp = EstimateHyperbolicTimestep<fluid>(md);
+    min_dt = std::min(min_dt, dt_hyp);
   }
 
   const auto &enable_cooling = hydro_pkg->Param<Cooling>("enable_cooling");
@@ -869,17 +844,34 @@ Real EstimateTimestep(MeshData<Real> *md) {
     min_dt = std::min(min_dt, tabular_cooling.EstimateTimeStep(md));
   }
 
-  // For RKL2 STS, the diffusive timestep is calculated separately in the driver
-  if (hydro_pkg->Param<DiffInt>("diffint") == DiffInt::unsplit) {
+  auto dt_diff = std::numeric_limits<Real>::max();
+  if (hydro_pkg->Param<DiffInt>("diffint") != DiffInt::none) {
     if (hydro_pkg->Param<Conduction>("conduction") != Conduction::none) {
-      min_dt = std::min(min_dt, EstimateConductionTimestep(md));
+      dt_diff = std::min(dt_diff, EstimateConductionTimestep(md));
     }
     if (hydro_pkg->Param<Viscosity>("viscosity") != Viscosity::none) {
-      min_dt = std::min(min_dt, EstimateViscosityTimestep(md));
+      dt_diff = std::min(dt_diff, EstimateViscosityTimestep(md));
     }
     if (hydro_pkg->Param<Resistivity>("resistivity") != Resistivity::none) {
-      min_dt = std::min(min_dt, EstimateResistivityTimestep(md));
+      dt_diff = std::min(dt_diff, EstimateResistivityTimestep(md));
+    }
+
+    // For unsplit ingegration use strict limit
+    if (hydro_pkg->Param<DiffInt>("diffint") == DiffInt::unsplit) {
+      min_dt = std::min(min_dt, dt_diff);
+      // and for RKL2 integration use limit taking into account the maxium ratio
+      // or not constrain limit further (which is why RKL2 is there in first place)
+    } else if (hydro_pkg->Param<DiffInt>("diffint") == DiffInt::rkl2) {
+      const auto max_dt_ratio = hydro_pkg->Param<Real>("rkl2_max_dt_ratio");
+      if (max_dt_ratio > 0.0 && dt_hyp / dt_diff > max_dt_ratio) {
+        min_dt = std::min(min_dt, max_dt_ratio * dt_diff);
+      }
+    } else {
+      PARTHENON_THROW("Looks like a a new diffusion integrator was implemented without "
+                      "taking into accout timestep contstraints yet.");
     }
+    auto dt_diff_param = hydro_pkg->Param<Real>("dt_diff");
+    hydro_pkg->UpdateParam("dt_diff", std::min(dt_diff, dt_diff_param));
   }
 
   if (ProblemEstimateTimestep != nullptr) {

src/hydro/hydro_driver.cpp

@@ -447,7 +447,10 @@ TaskCollection HydroDriver::MakeTaskCollection(BlockList_t &blocks, int stage) {
   // Calculate hyperbolic divergence cleaning speed
   // TODO(pgrete) Calculating mindx is only required after remeshing. Need to find a clean
   // solution for this one-off global reduction.
-  if (hydro_pkg->Param<bool>("calc_c_h") && (stage == 1)) {
+  // TODO(PG) move this to PreStepMeshUserWorkInLoop
+  if ((hydro_pkg->Param<bool>("calc_c_h") ||
+       hydro_pkg->Param<DiffInt>("diffint") != DiffInt::none) &&
+      (stage == 1)) {
     // need to make sure that there's only one region in order to MPI_reduce to work
     TaskRegion &single_task_region = tc.AddRegion(1);
     auto &tl = single_task_region[0];
@@ -468,14 +471,16 @@ TaskCollection HydroDriver::MakeTaskCollection(BlockList_t &blocks, int stage) {
     reduce_c_h = tl.AddTask(
         prev_task,
         [](StateDescriptor *hydro_pkg) {
-          Real mins[2];
+          Real mins[3];
           mins[0] = hydro_pkg->Param<Real>("mindx");
           mins[1] = hydro_pkg->Param<Real>("dt_hyp");
-          PARTHENON_MPI_CHECK(MPI_Allreduce(MPI_IN_PLACE, mins, 2, MPI_PARTHENON_REAL,
+          mins[2] = hydro_pkg->Param<Real>("dt_diff");
+          PARTHENON_MPI_CHECK(MPI_Allreduce(MPI_IN_PLACE, mins, 3, MPI_PARTHENON_REAL,
                                             MPI_MIN, MPI_COMM_WORLD));
 
           hydro_pkg->UpdateParam("mindx", mins[0]);
           hydro_pkg->UpdateParam("dt_hyp", mins[1]);
+          hydro_pkg->UpdateParam("dt_diff", mins[2]);
           return TaskStatus::complete;
         },
         hydro_pkg.get());
@@ -657,14 +662,17 @@ TaskCollection HydroDriver::MakeTaskCollection(BlockList_t &blocks, int stage) {
 
   // Single task in single (serial) region to reset global vars used in reductions in the
   // first stage.
-  if (stage == integrator->nstages && hydro_pkg->Param<bool>("calc_c_h")) {
+  if (stage == integrator->nstages &&
+      (hydro_pkg->Param<bool>("calc_c_h") ||
+       hydro_pkg->Param<DiffInt>("diffint") != DiffInt::none)) {
     TaskRegion &reset_reduction_vars_region = tc.AddRegion(1);
     auto &tl = reset_reduction_vars_region[0];
     tl.AddTask(
         none,
         [](StateDescriptor *hydro_pkg) {
           hydro_pkg->UpdateParam("mindx", std::numeric_limits<Real>::max());
           hydro_pkg->UpdateParam("dt_hyp", std::numeric_limits<Real>::max());
+          hydro_pkg->UpdateParam("dt_diff", std::numeric_limits<Real>::max());
           return TaskStatus::complete;
         },
         hydro_pkg.get());

src/main.cpp

@@ -63,8 +63,8 @@ int main(int argc, char *argv[]) {
   } else if (problem == "cloud") {
     pman.app_input->InitUserMeshData = cloud::InitUserMeshData;
     pman.app_input->ProblemGenerator = cloud::ProblemGenerator;
-    pman.app_input->boundary_conditions[parthenon::BoundaryFace::inner_x2] =
-        cloud::InflowWindX2;
+    pman.app_input->RegisterBoundaryCondition(parthenon::BoundaryFace::inner_x2,
+                                              "cloud_inflow_x2", cloud::InflowWindX2);
     Hydro::ProblemCheckRefinementBlock = cloud::ProblemCheckRefinementBlock;
   } else if (problem == "blast") {
     pman.app_input->InitUserMeshData = blast::InitUserMeshData;