From 580e5739d1119c409ec6ed22878a52bbab1a4bcc Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Sat, 9 Nov 2024 12:15:49 +0100 Subject: [PATCH 1/5] Bump Kokkos 4.4.1, Parth vov-fix --- external/Kokkos | 2 +- external/parthenon | 2 +- inputs/cloud.in | 2 +- src/main.cpp | 4 ++-- 4 files changed, 5 insertions(+), 5 deletions(-) diff --git a/external/Kokkos b/external/Kokkos index 08ceff92..15dc143e 160000 --- a/external/Kokkos +++ b/external/Kokkos @@ -1 +1 @@ -Subproject commit 08ceff92bcf3a828844480bc1e6137eb74028517 +Subproject commit 15dc143e5f39949eece972a798e175c4b463d4b8 diff --git a/external/parthenon b/external/parthenon index ec61c9cb..e00017c6 160000 --- a/external/parthenon +++ b/external/parthenon @@ -1 +1 @@ -Subproject commit ec61c9cb102d6a67a4dd03e7d7fa4e10c82c1032 +Subproject commit e00017c6ece547f3adf9a16b94ddf0bdf8bb4773 diff --git a/inputs/cloud.in b/inputs/cloud.in index 002febd0..e7e05f69 100644 --- a/inputs/cloud.in +++ b/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 diff --git a/src/main.cpp b/src/main.cpp index f239b197..23cdb825 100644 --- a/src/main.cpp +++ b/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; From 9a7ede7423123d4b0d62f44777dbcb2cbfe34c4b Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Sat, 9 Nov 2024 19:09:54 +0100 Subject: [PATCH 2/5] Store and restore dt_X and centralize dt_diff calc --- docs/input.md | 3 -- src/hydro/hydro.cpp | 91 +++++++++++++++++--------------------- src/hydro/hydro_driver.cpp | 16 +++++-- 3 files changed, 52 insertions(+), 58 deletions(-) diff --git a/docs/input.md b/docs/input.md index 73f2b23b..0395a448 100644 --- a/docs/input.md +++ b/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. diff --git a/src/hydro/hydro.cpp b/src/hydro/hydro.cpp index 7bc95be7..3861e53a 100644 --- a/src/hydro/hydro.cpp +++ b/src/hydro/hydro.cpp @@ -4,6 +4,7 @@ // Licensed under the BSD 3-Clause License (the "LICENSE"). //======================================================================================== +#include #include #include #include @@ -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 &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::rkl2) { - auto dt_diff = std::numeric_limits::max(); - if (hydro_pkg->Param("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::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::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("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 Real HydroHst(MeshData *md) { @@ -252,17 +217,20 @@ std::shared_ptr Initialize(ParameterInput *pin) { auto glmmhd_alpha = pin->GetOrAddReal("hydro", "glmmhd_alpha", 0.1); pkg->AddParam("glmmhd_alpha", glmmhd_alpha); calc_c_h = true; - pkg->AddParam("c_h", 0.0, true); // hyperbolic divergence cleaning speed - // global minimum dx (used to calc c_h) - pkg->AddParam("mindx", std::numeric_limits::max(), true); - // hyperbolic timestep constraint - pkg->AddParam("dt_hyp", std::numeric_limits::max(), true); + pkg->AddParam( + "c_h", 0.0, Params::Mutability::Restart); // hyperbolic divergence cleaning speed } else { PARTHENON_FAIL("AthenaPK hydro: Unknown fluid method."); } pkg->AddParam<>("fluid", fluid); pkg->AddParam<>("nhydro", nhydro); pkg->AddParam<>("calc_c_h", calc_c_h); + // global minimum dx (used to calc c_h) + pkg->AddParam("mindx", std::numeric_limits::max(), + Params::Mutability::Restart); + // hyperbolic timestep constraint + pkg->AddParam("dt_hyp", std::numeric_limits::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 +601,13 @@ std::shared_ptr Initialize(ParameterInput *pin) { "Options are: none, unsplit, rkl2"); } if (diffint != DiffInt::none) { - pkg->AddParam("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("cfl")); pkg->AddParam<>("cfl_diff", cfl_diff); } + pkg->AddParam("dt_diff", std::numeric_limits::max(), + Params::Mutability::Restart); // diffusive timestep constraint pkg->AddParam<>("diffint", diffint); if (fluid == Fluid::euler) { @@ -855,9 +824,12 @@ Real EstimateTimestep(MeshData *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::max(); + auto dt_hyp = std::numeric_limits::max(); - if (hydro_pkg->Param("calc_dt_hyp")) { - min_dt = std::min(min_dt, EstimateHyperbolicTimestep(md)); + const auto calc_dt_hyp = hydro_pkg->Param("calc_dt_hyp"); + if (calc_dt_hyp) { + dt_hyp = EstimateHyperbolicTimestep(md); + min_dt = std::min(min_dt, dt_hyp); } const auto &enable_cooling = hydro_pkg->Param("enable_cooling"); @@ -869,17 +841,34 @@ Real EstimateTimestep(MeshData *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::unsplit) { + auto dt_diff = std::numeric_limits::max(); + if (hydro_pkg->Param("diffint") != DiffInt::none) { if (hydro_pkg->Param("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::none) { - min_dt = std::min(min_dt, EstimateViscosityTimestep(md)); + dt_diff = std::min(dt_diff, EstimateViscosityTimestep(md)); } if (hydro_pkg->Param("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::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::rkl2) { + const auto max_dt_ratio = hydro_pkg->Param("rkl2_max_dt_ratio"); + if (calc_dt_hyp && 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("dt_diff"); + hydro_pkg->UpdateParam("dt_diff", std::min(dt_diff, dt_diff_param)); } if (ProblemEstimateTimestep != nullptr) { diff --git a/src/hydro/hydro_driver.cpp b/src/hydro/hydro_driver.cpp index 36bdce91..eeb29d3d 100644 --- a/src/hydro/hydro_driver.cpp +++ b/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("calc_c_h") && (stage == 1)) { + // TODO(PG) move this to PreStepMeshUserWorkInLoop + if ((hydro_pkg->Param("calc_c_h") || + hydro_pkg->Param("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("mindx"); mins[1] = hydro_pkg->Param("dt_hyp"); - PARTHENON_MPI_CHECK(MPI_Allreduce(MPI_IN_PLACE, mins, 2, MPI_PARTHENON_REAL, + mins[2] = hydro_pkg->Param("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,7 +662,9 @@ 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("calc_c_h")) { + if (stage == integrator->nstages && + (hydro_pkg->Param("calc_c_h") || + hydro_pkg->Param("diffint") != DiffInt::none)) { TaskRegion &reset_reduction_vars_region = tc.AddRegion(1); auto &tl = reset_reduction_vars_region[0]; tl.AddTask( @@ -665,6 +672,7 @@ TaskCollection HydroDriver::MakeTaskCollection(BlockList_t &blocks, int stage) { [](StateDescriptor *hydro_pkg) { hydro_pkg->UpdateParam("mindx", std::numeric_limits::max()); hydro_pkg->UpdateParam("dt_hyp", std::numeric_limits::max()); + hydro_pkg->UpdateParam("dt_diff", std::numeric_limits::max()); return TaskStatus::complete; }, hydro_pkg.get()); From 867fb3a2da10e9881be788f37fc177084b211ebe Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Sun, 10 Nov 2024 11:45:30 +0100 Subject: [PATCH 3/5] Enforce max ratio (i.e., limit number of stages) even without dt_hyp --- src/hydro/hydro.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/hydro/hydro.cpp b/src/hydro/hydro.cpp index 3861e53a..4e2e5eec 100644 --- a/src/hydro/hydro.cpp +++ b/src/hydro/hydro.cpp @@ -860,7 +860,7 @@ Real EstimateTimestep(MeshData *md) { // or not constrain limit further (which is why RKL2 is there in first place) } else if (hydro_pkg->Param("diffint") == DiffInt::rkl2) { const auto max_dt_ratio = hydro_pkg->Param("rkl2_max_dt_ratio"); - if (calc_dt_hyp && max_dt_ratio > 0.0 && dt_hyp / dt_diff > 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 { From fcb5f397aca9032410050da03d45985148cd91a3 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Sun, 10 Nov 2024 12:47:59 +0100 Subject: [PATCH 4/5] Ensure c_h is present --- src/hydro/hydro.cpp | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/hydro/hydro.cpp b/src/hydro/hydro.cpp index 4e2e5eec..5720346c 100644 --- a/src/hydro/hydro.cpp +++ b/src/hydro/hydro.cpp @@ -217,14 +217,17 @@ std::shared_ptr Initialize(ParameterInput *pin) { auto glmmhd_alpha = pin->GetOrAddReal("hydro", "glmmhd_alpha", 0.1); pkg->AddParam("glmmhd_alpha", glmmhd_alpha); calc_c_h = true; - pkg->AddParam( - "c_h", 0.0, Params::Mutability::Restart); // hyperbolic divergence cleaning speed } 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( + "c_h", 0.0, Params::Mutability::Restart); // hyperbolic divergence cleaning speed // global minimum dx (used to calc c_h) pkg->AddParam("mindx", std::numeric_limits::max(), Params::Mutability::Restart); From f4d3b824b0a55fd749a013cc84fa125ff6b0eaa3 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Sun, 10 Nov 2024 17:43:58 +0100 Subject: [PATCH 5/5] Add Changelog --- CHANGELOG.md | 3 +++ 1 file changed, 3 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index c95a86ae..b1d49a38 100644 --- a/CHANGELOG.md +++ b/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)