Merge branch 'main' into BenWibking/restart-reproducer

CHANGELOG.md

@@ -7,13 +7,15 @@
 - [[PR 1]](https://github.com/parthenon-hpc-lab/athenapk/pull/1) Add isotropic thermal conduction and RKL2 supertimestepping
 
 ### Changed (changing behavior/API/variables/...)
+- [[PR 122]](https://github.com/parthenon-hpc-lab/athenapk/pull/122) Fixed sqrt(4pi) factor in CGS Gauss unit and add unit doc
 - [[PR 119]](https://github.com/parthenon-hpc-lab/athenapk/pull/119) Fixed Athena++ paper test case for KHI pgen. Added turbulence pgen doc.
 - [[PR 97]](https://github.com/parthenon-hpc-lab/athenapk/pull/97) Fixed Schure cooling curve. Removed SD one. Added description of cooling function conventions.
 - [[PR 84]](https://github.com/parthenon-hpc-lab/athenapk/pull/84) Bump Parthenon to latest develop (2024-02-15)
 
 ### 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
@@ -23,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)

README.md

@@ -131,21 +131,8 @@ the `file_type = hdf5` format, see
 [VisIt](https://wci.llnl.gov/simulation/computer-codes/visit/).
 In ParaView, select the "XDMF Reader" when prompted.
 
-2. With [yt](https://yt-project.org/) -- though currently through a custom frontend
-that is not yet part of the main yt branch and, thus, has to be installed manually, e.g.,
-as follows:
-```bash
-cd ~/src # or any other folder of choice
-git clone https://github.com/forrestglines/yt.git
-cd yt
-git checkout parthenon-frontend
-
-# If you're using conda or virtualenv
-pip install -e .
-# OR alternatively, if you using the plain Python environment
-pip install --user -e .
-```
-Afterwards, `*.phdf` files can be read as usual with `yt.load()`.
+2. With [yt](https://yt-project.org/)
+As of versions >=4.4 `*.phdf` files can be read as usual with `yt.load()`.
 
 3. Using [Ascent](https://github.com/Alpine-DAV/ascent) (for in situ visualization and analysis).
 This requires Ascent to be installed/available at compile time of AthenaPK.

docs/README.md

@@ -13,6 +13,7 @@ The documentation currently includes
   - [Notebooks to calculate cooling tables from literature](cooling)
 - [Brief notes on developing code for AthenaPK](development.md)
 - [How to add a custom/user problem generator](pgen.md)
+- [Units](units.md)
 - Detailed descriptions of more complex problem generators
   - [Galaxy Cluster and Cluster-like Problem Setup](cluster.md)
   - [Driven turbulence](turbulence.md)

docs/input.md

@@ -67,6 +67,10 @@ To control the floors, following parameters can be set in the `<hydro>` block:
 *Note* the pressure floor will take precedence over the temperature floor in the
 conserved to primitive conversion if both are defined.
 
+#### Units
+
+See(here)[units.md].
+
 #### Diffusive processes
 
 Diffusive processes in AthenaPK can be configured in the `<diffusion>` block of the input file.
@@ -114,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.

docs/units.md

@@ -0,0 +1,58 @@
+# AthenaPK units
+
+## General unit system
+
+Internally, all calculations are done in "code units" and there are no conversions between
+code and physical units during runtime (with excetion like the temperature when cooling is
+being used).
+Therefore, in general no units need to be prescribed to run a simulation.
+
+If units are required (e.g., if cooling is used and, thus, a conversion between internal energy
+in code units and physical temperature is required) they are configured in the input block
+as follows:
+
+```
+<units>
+code_length_cgs = 3.085677580962325e+24 # 1 Mpc in cm
+code_mass_cgs = 1.98841586e+47          # 1e14 Msun in g
+code_time_cgs = 3.15576e+16             # 1 Gyr in s
+```
+
+This information will also be used by postprocessing tools (like yt) to convert between
+code units and a physical unit system (like cgs).
+
+Moreover, internally a set of factors from code to cgs units are available to process conversions
+if required (e.g., from the input file).
+
+For example, for an input parameter (in the input file) like
+
+```
+<problem/cloud>
+r0_cgs = 3.085677580962325e+20          # 100 pc
+```
+
+the conversion should happen in the problem generator lik
+
+```c++
+  r_cloud = pin->GetReal("problem/cloud", "r0_cgs") / units.code_length_cgs();
+```
+
+so that the resulting quantity is internally in code units (here code length).
+
+It is highly recommended to be *very* explicit/specific about units everywhere (as it is
+a common source of confusion) like adding the `_cgs` suffix to the parameter in the
+input file above.
+
+## Magnetic units
+
+Internally, AthenaPK (and almost all MHD codes) use
+[Heaviside-Lorentz units](https://en.wikipedia.org/wiki/Heaviside%E2%80%93Lorentz_units),
+where the magnetic field is transformed from $B \rightarrow B / \sqrt{4 \pi}$.
+(See also the note in the
+[Castro documentation](https://amrex-astro.github.io/Castro/docs/mhd.html) about this.)
+
+So when converting from CGS-Gaussian units to code units, it is necessary to divide
+by $\sqrt{4 \pi}$ (in addition to the base dimensional factors).
+This is automatically handled by the `units.code_magnetic_cgs()` factors.
+
+

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;

src/units.hpp

@@ -90,7 +90,8 @@ class Units {
     return code_energy_cgs() / kev_cgs * code_length_cgs() * code_length_cgs();
   }
   parthenon::Real code_magnetic_cgs() const {
-    return sqrt(code_mass_cgs()) / sqrt(code_length_cgs()) / code_time_cgs();
+    return std::sqrt(4.0 * M_PI) * sqrt(code_mass_cgs()) / sqrt(code_length_cgs()) /
+           code_time_cgs();
   }
 
   // Physical Constants in code units