Skip to content

Commit

Permalink
fix computation of dt when we have multiple levels (erf-model#1806)
Browse files Browse the repository at this point in the history
  • Loading branch information
@asalmgren
asalmgren authored Sep 13, 2024
1 parent b8b4551 commit 7fd7810
Show file tree
Hide file tree
Showing 3 changed files with 48 additions and 38 deletions.
4 changes: 2 additions & 2 deletions Source/ERF.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -835,8 +835,8 @@ private:
static amrex::Real change_max;

// Fixed dt for level 0 timesteps (only used if positive)
static amrex::Real fixed_dt;
static amrex::Real fixed_fast_dt;
amrex::Vector<amrex::Real> fixed_dt;
amrex::Vector<amrex::Real> fixed_fast_dt;
static int fixed_mri_dt_ratio;

// how often each level regrids the higher levels of refinement
Expand Down
58 changes: 34 additions & 24 deletions Source/ERF.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,8 +30,6 @@ SolverChoice ERF::solverChoice;

// Time step control
Real ERF::cfl = 0.8;
Real ERF::fixed_dt = -1.0;
Real ERF::fixed_fast_dt = -1.0;
Real ERF::init_shrink = 1.0;
Real ERF::change_max = 1.1;
int ERF::fixed_mri_dt_ratio = 0;
Expand Down Expand Up @@ -1371,8 +1369,18 @@ ERF::ReadParameters ()
pp.query("init_shrink", init_shrink);
pp.query("change_max", change_max);

pp.query("fixed_dt", fixed_dt);
pp.query("fixed_fast_dt", fixed_fast_dt);
fixed_dt.resize(max_level+1,-1.);
fixed_fast_dt.resize(max_level+1,-1.);

pp.query("fixed_dt", fixed_dt[0]);
pp.query("fixed_fast_dt", fixed_fast_dt[0]);

for (int lev = 1; lev <= max_level; lev++)
{
fixed_dt[lev] = fixed_dt[lev-1] / static_cast<Real>(MaxRefRatio(lev-1));
fixed_fast_dt[lev] = fixed_fast_dt[lev-1] / static_cast<Real>(MaxRefRatio(lev-1));
}

pp.query("fixed_mri_dt_ratio", fixed_mri_dt_ratio);

// How to initialize
Expand Down Expand Up @@ -1514,7 +1522,7 @@ ERF::ReadParameters ()
void
ERF::ParameterSanityChecks ()
{
AMREX_ALWAYS_ASSERT(cfl > 0. || fixed_dt > 0.);
AMREX_ALWAYS_ASSERT(cfl > 0. || fixed_dt[0] > 0.);

// We don't allow use_real_bcs to be true if init_type is not either real or metgrid
AMREX_ALWAYS_ASSERT(!use_real_bcs || ((init_type == "real") || (init_type == "metgrid")) );
Expand Down Expand Up @@ -1561,31 +1569,33 @@ ERF::ParameterSanityChecks ()
Abort("If you specify fixed_mri_dt_ratio, it must be even");
}

// We ignore fixed_fast_dt if not substepping
if (solverChoice.no_substepping && fixed_fast_dt > 0.) {
fixed_fast_dt = -1.0;
Warning("fixed_fast_dt will be ignored since we are not substepping");


// If both fixed_dt and fast_dt are specified, their ratio must be an even integer
} else if (fixed_dt > 0. && fixed_fast_dt > 0. && fixed_mri_dt_ratio <= 0)
for (int lev = 0; lev <= max_level; lev++)
{
Real eps = 1.e-12;
int ratio = static_cast<int>( ( (1.0+eps) * fixed_dt ) / fixed_fast_dt );
if (fixed_dt / fixed_fast_dt != ratio)
// We ignore fixed_fast_dt if not substepping
if (solverChoice.no_substepping) {
fixed_fast_dt[lev] = -1.0;
}

// If both fixed_dt and fast_dt are specified, their ratio must be an even integer
if (fixed_dt[lev] > 0. && fixed_fast_dt[lev] > 0. && fixed_mri_dt_ratio <= 0)
{
Abort("Ratio of fixed_dt to fixed_fast_dt must be an even integer");
Real eps = 1.e-12;
int ratio = static_cast<int>( ( (1.0+eps) * fixed_dt[lev] ) / fixed_fast_dt[lev] );
if (fixed_dt[lev] / fixed_fast_dt[lev] != ratio)
{
Abort("Ratio of fixed_dt to fixed_fast_dt must be an even integer");
}
}
}

// If all three are specified, they must be consistent
if (fixed_dt > 0. && fixed_fast_dt > 0. && fixed_mri_dt_ratio > 0)
{
if (fixed_dt / fixed_fast_dt != fixed_mri_dt_ratio)
// If all three are specified, they must be consistent
if (fixed_dt[lev] > 0. && fixed_fast_dt[lev] > 0. && fixed_mri_dt_ratio > 0)
{
Abort("Dt is over-specfied");
if (fixed_dt[lev] / fixed_fast_dt[lev] != fixed_mri_dt_ratio)
{
Abort("Dt is over-specfied");
}
}
}
} // lev

if (solverChoice.coupling_type == CouplingType::TwoWay && cf_width > 0) {
Abort("For two-way coupling you must set cf_width = 0");
Expand Down
24 changes: 12 additions & 12 deletions Source/TimeIntegration/ERF_ComputeTimestep.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,7 @@ ERF::ComputeDt (int step)
dt_0 = amrex::min(dt_0, n_factor*dt_tmp[lev]);
if (step == 0){
dt_0 *= init_shrink;
if (verbose) {
if (verbose && init_shrink != 1.0) {
Print() << "Timestep 0: shrink initial dt at level " << lev << " by " << init_shrink << std::endl;
}
}
Expand Down Expand Up @@ -148,7 +148,7 @@ ERF::estTimeStep (int level, long& dt_fast_ratio) const
estdt_lowM = cfl / estdt_lowM_inv;

if (verbose) {
if (fixed_dt <= 0.0) {
if (fixed_dt[level] <= 0.0) {
Print() << "Using cfl = " << cfl << std::endl;
Print() << "Compressible dt at level " << level << ": " << estdt_comp << std::endl;
if (estdt_lowM_inv > 0.0_rt) {
Expand All @@ -158,26 +158,26 @@ ERF::estTimeStep (int level, long& dt_fast_ratio) const
}
}

if (fixed_dt > 0.0) {
if (fixed_dt[level] > 0.0) {
Print() << "Based on cfl of 1.0 " << std::endl;
Print() << "Compressible dt at level " << level << " would be: " << estdt_comp/cfl << std::endl;
if (estdt_lowM_inv > 0.0_rt) {
Print() << "Incompressible dt at level " << level << " would be: " << estdt_lowM/cfl << std::endl;
} else {
Print() << "Incompressible dt at level " << level << " would be undefined " << std::endl;
}
Print() << "Fixed dt at level " << level << " is: " << fixed_dt << std::endl;
if (fixed_fast_dt > 0.0) {
Print() << "Fixed fast dt at level " << level << " is: " << fixed_fast_dt << std::endl;
Print() << "Fixed dt at level " << level << " is: " << fixed_dt[level] << std::endl;
if (fixed_fast_dt[level] > 0.0) {
Print() << "Fixed fast dt at level " << level << " is: " << fixed_fast_dt[level] << std::endl;
}
}
}

if (!l_no_substepping) {
if (fixed_dt > 0. && fixed_fast_dt > 0.) {
dt_fast_ratio = static_cast<long>( fixed_dt / fixed_fast_dt );
} else if (fixed_dt > 0.) {
dt_fast_ratio = static_cast<long>( std::ceil((fixed_dt/estdt_comp)) );
if (fixed_dt[level] > 0. && fixed_fast_dt[level] > 0.) {
dt_fast_ratio = static_cast<long>( fixed_dt[level] / fixed_fast_dt[level] );
} else if (fixed_dt[level] > 0.) {
dt_fast_ratio = static_cast<long>( std::ceil((fixed_dt[level]/estdt_comp)) );
} else {
dt_fast_ratio = (estdt_lowM_inv > 0.0) ? static_cast<long>( std::ceil((estdt_lowM/estdt_comp)) ) : 1;
}
Expand All @@ -198,8 +198,8 @@ ERF::estTimeStep (int level, long& dt_fast_ratio) const
}
} // if substepping

if (fixed_dt > 0.0) {
return fixed_dt;
if (fixed_dt[level] > 0.0) {
return fixed_dt[level];
} else {
// Incompressible (substepping is not allowed)
if (l_incompressible) {
Expand Down

0 comments on commit 7fd7810

Please sign in to comment.