Skip to content

Commit

Permalink
Store and restore dt_X and centralize dt_diff calc
Browse files Browse the repository at this point in the history
  • Loading branch information
@pgrete
pgrete committed Nov 9, 2024
1 parent 580e573 commit 9a7ede7
Show file tree
Hide file tree
Showing 3 changed files with 52 additions and 58 deletions.
3 changes: 0 additions & 3 deletions docs/input.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down
91 changes: 40 additions & 51 deletions src/hydro/hydro.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@
// Licensed under the BSD 3-Clause License (the "LICENSE").
//========================================================================================

#include <algorithm>
#include <limits>
#include <memory>
#include <string>
Expand All @@ -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"
Expand Down Expand Up @@ -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) {
Expand Down Expand Up @@ -252,17 +217,20 @@ 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);
pkg->AddParam<Real>(
"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<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
Expand Down Expand Up @@ -633,12 +601,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) {
Expand Down Expand Up @@ -855,9 +824,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");
Expand All @@ -869,17 +841,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 (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<Real>("dt_diff");
hydro_pkg->UpdateParam("dt_diff", std::min(dt_diff, dt_diff_param));
}

if (ProblemEstimateTimestep != nullptr) {
Expand Down
16 changes: 12 additions & 4 deletions src/hydro/hydro_driver.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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];
Expand All @@ -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());
Expand Down Expand Up @@ -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());
Expand Down

0 comments on commit 9a7ede7

Please sign in to comment.