Skip to content

Improve BackwardEuler #1987

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
BenWibking wants to merge 1 commit into from
Closed

Improve BackwardEuler #1987

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
213 changes: 160 additions & 53 deletions integration/BackwardEuler/be_integrator.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,19 +28,25 @@
///
template <typename BurnT, typename BeT>
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
int single_step (BurnT& state, BeT& be, const amrex::Real dt)
int single_step (BurnT& state, BeT& be, const amrex::Real time_start,
const amrex::Real dt)
{
constexpr int int_neqs = integrator_neqs<BurnT>();

int ierr = IERR_SUCCESS;
bool converged = false;
const amrex::Real time_end = time_start + dt;
const amrex::Real time_save = be.t;

// create our current guess for the solution -- just as a first
// order explicit prediction

amrex::Array1D<amrex::Real, 1, int_neqs> ydot;

rhs(be.t, state, be, ydot);
const bool do_clean_state = ! integrator_rp::use_number_densities;

be.t = time_start;
rhs(time_start, state, be, ydot, false, do_clean_state);

be.n_rhs += 1;

Expand All @@ -57,9 +63,9 @@ int single_step (BurnT& state, BeT& be, const amrex::Real dt)

// work with the current guess

if (integrator_rp::do_corrector_validation) {
if (integrator_rp::do_corrector_validation && !integrator_rp::use_number_densities) {
#ifdef SDC
const amrex::Real rho_current = state.rho_orig + be.t * state.ydot_a[SRHO];
const amrex::Real rho_current = state.rho_orig + time_end * state.ydot_a[SRHO];
const amrex::Real thresh = integrator_rp::species_failure_tolerance * rho_current;
#else
const amrex::Real thresh = integrator_rp::species_failure_tolerance;
Expand All @@ -80,13 +86,14 @@ int single_step (BurnT& state, BeT& be, const amrex::Real dt)

// get the ydots for our current guess of y

rhs(be.t, state, be, ydot);
be.t = time_end;
rhs(time_end, state, be, ydot, false, do_clean_state);
be.n_rhs += 1;

// construct the Jacobian

if (be.jacobian_type == 1) {
jac(be.t, state, be, be.jac);
jac(time_end, state, be, be.jac, do_clean_state);
} else {
jac_info_t jac_info;
jac_info.h = dt;
Expand Down Expand Up @@ -183,6 +190,8 @@ int single_step (BurnT& state, BeT& be, const amrex::Real dt)

}

be.t = time_save;

return ierr;

}
Expand All @@ -202,22 +211,27 @@ int be_integrator (BurnT& state, BeT& be)
if (integrator_rp::do_single_step == 1) {
// Do a single step to the final time
const amrex::Real dt_single_step = be.tout - be.t;
ierr = single_step(state, be, dt_single_step);
ierr = single_step(state, be, be.t, dt_single_step);
if (ierr == IERR_SUCCESS) {
be.t = be.tout;
++be.n_step;
}
++be.n_step;
return ierr;
}

// estimate the timestep

amrex::Array1D<amrex::Real, 1, int_neqs> ydot;
rhs(be.t, state, be, ydot);
rhs(be.t, state, be, ydot, false,
! integrator_rp::use_number_densities);

be.n_rhs += 1;

amrex::Real dt_sub = initial_react_dt(state, be, ydot);
int n_attempt = 0;
int n_reject = 0;
int n_consecutive_reject = 0;
bool previous_step_rejected = false;

// main timestepping loop

Expand All @@ -236,58 +250,70 @@ int be_integrator (BurnT& state, BeT& be)
dt_sub = be.tout - be.t;
}

// our strategy is to take 2 steps at dt/2 and one at dt and
// to compute the error from those

// first do 2 (fine) dt/2 steps
const amrex::Real t_old = be.t;

amrex::Array1D<amrex::Real, 1, int_neqs> y_fine;
// Estimate the local truncation error by comparing the backward
// Euler candidate against a forward Euler candidate over the same dt.

// keep track of whether we have a valid fine solution
// this means no errors from either half step
bool have_fine_solution = false;

ierr = single_step(state, be, dt_sub/2);
if (ierr == IERR_SUCCESS) {
ierr = single_step(state, be, dt_sub/2);

if (ierr == IERR_SUCCESS) {
// store the fine dt solution
for (int n = 1; n <= int_neqs; ++n) {
y_fine(n) = be.y(n);
}
have_fine_solution = true;
amrex::Array1D<amrex::Real, 1, int_neqs> y_forward;
rhs(t_old, state, be, ydot, false,
! integrator_rp::use_number_densities);
be.n_rhs += 1;

// now do a single (coarse) dt step
// first reset the solution
for (int n = 1; n <= int_neqs; ++n) {
be.y(n) = y_old(n);
}
ierr = single_step(state, be, dt_sub);
}
for (int n = 1; n <= int_neqs; ++n) {
y_forward(n) = y_old(n) + dt_sub * ydot(n);
}

ierr = single_step(state, be, t_old, dt_sub);

amrex::Real rel_error = 0.0_rt;
amrex::Real max_error = 0.0_rt;
bool step_success = false;
if (ierr == IERR_SUCCESS && have_fine_solution) {
int err_idx = -1;
const amrex::Real dt_attempt = dt_sub;
const bool have_error_estimate = ierr == IERR_SUCCESS;
if (have_error_estimate) {
// define a weight for each variable to use in checking the error

amrex::Array1D<amrex::Real, 1, int_neqs> w;
for (int n = 1; n <= NumSpec; n++) {
w(n) = 1.0_rt / (be.rtol_spec * std::abs(y_fine(n)) + be.atol_spec);
w(n) = 1.0_rt / (be.rtol_spec * std::abs(be.y(n)) + be.atol_spec);
}
w(net_ienuc) = 1.0_rt / (be.rtol_enuc * std::abs(y_fine(net_ienuc)) + be.atol_enuc);
w(net_ienuc) = 1.0_rt / (be.rtol_enuc * std::abs(be.y(net_ienuc)) + be.atol_enuc);

// now look for w |y_fine - y_coarse| < 1
// Now look for weighted |y_backward - y_forward| < 1. Use the
// same RMS norm style as VODE for the controller, while keeping
// the maximum component around for diagnostics.

for (int n = 1; n <= NumSpec; n++) {
rel_error = amrex::max(rel_error, w(n) * std::abs(y_fine(n) - be.y(n)));
const amrex::Real component_error = w(n) * std::abs(be.y(n) - y_forward(n));
rel_error += component_error * component_error;
if (component_error > max_error) {
max_error = component_error;
err_idx = n;
}
}
rel_error = amrex::max(rel_error, w(net_ienuc) * std::abs(y_fine(net_ienuc) - be.y(net_ienuc)));
const amrex::Real enuc_error = w(net_ienuc) * std::abs(be.y(net_ienuc) - y_forward(net_ienuc));
rel_error += enuc_error * enuc_error;
if (enuc_error > max_error) {
max_error = enuc_error;
err_idx = net_ienuc;
}
constexpr amrex::Real TQ_ORDER1 = 2.0_rt;
rel_error = std::sqrt(rel_error / int_neqs) / TQ_ORDER1;

if (rel_error < 1.0_rt) {
step_success = true;
}

if (step_success && integrator_rp::use_number_densities) {
for (int n = 1; n <= NumSpec; ++n) {
if (be.y(n) < -integrator_rp::species_failure_tolerance) {
step_success = false;
break;
}
}
}
}

if (ierr == IERR_SUCCESS && step_success) {
Expand All @@ -296,19 +322,52 @@ int be_integrator (BurnT& state, BeT& be)
// (linear algebra, etc.) and we met our error
// goals.

// y_fine has the current best solution

be.t += dt_sub;

for (int n = 1; n <= int_neqs; ++n) {
be.y(n) = y_fine(n);
// Can we potentially increase the timestep? Backward-Euler has
// a local truncation error of O(dt**2). Use a VODE-like biased
// growth factor and threshold, but continue to estimate the error
// from the forward Euler / backward Euler comparison.

constexpr amrex::Real ADDON = 1.0e-6_rt;
constexpr amrex::Real BIAS2 = 2.0_rt;
constexpr amrex::Real ETAMX = 10.0_rt;
constexpr amrex::Real THRESH = 1.5_rt;
constexpr amrex::Real tiny_error = 1.0e-30_rt;

amrex::Real dt_factor = 1.0_rt;
const amrex::Real safe_error = amrex::max(rel_error, tiny_error);
const amrex::Real eta = 1.0_rt / (std::sqrt(BIAS2 * safe_error) + ADDON);

if (! previous_step_rejected && eta >= THRESH) {
dt_factor = amrex::min(eta, ETAMX);
}

// can we potentially increase the timestep?
// backward-Euler has a local truncation error of dt**2
dt_sub *= dt_factor;
previous_step_rejected = false;
n_consecutive_reject = 0;

#ifndef AMREX_USE_GPU
if (integrator_rp::burner_verbose &&
(be.n_step < 40 || be.n_step % 1000 == 0)) {
std::cout << "[BE_STEP] accept n=" << be.n_step
<< " t_new=" << be.t
<< " dt=" << dt_attempt
<< " rel_error=" << rel_error
<< " max_error=" << max_error
<< " err_idx=" << err_idx
<< " eta=" << eta
<< " dt_factor=" << dt_factor
<< " dt_next=" << dt_sub
<< " y=";
for (int n = 1; n <= NumSpec; ++n) {
std::cout << be.y(n) << " ";
}
std::cout << std::endl;
}
#endif

amrex::Real dt_new = dt_sub * std::sqrt(1.0_rt / rel_error);
dt_sub = amrex::Clamp(dt_new, dt_sub / 2.0, 2.0 * dt_sub);
++be.n_step;

} else {

Expand All @@ -317,20 +376,68 @@ int be_integrator (BurnT& state, BeT& be)
be.y(n) = y_old(n);
}

// adjust the timestep and try again
dt_sub /= 2;
// Adjust the timestep and try again. This mirrors VODE's
// same-order error-failure step-size update for order 1, using
// the forward Euler / backward Euler error estimate.

constexpr amrex::Real ADDON = 1.0e-6_rt;
constexpr amrex::Real BIAS2 = 2.0_rt;
constexpr amrex::Real ETAMIN = 0.1_rt;
constexpr amrex::Real ETAMXF = 0.2_rt;
constexpr amrex::Real tiny_error = 1.0e-30_rt;

amrex::Real dt_factor = ETAMIN;
if (have_error_estimate) {
const amrex::Real safe_error = amrex::max(rel_error, tiny_error);
dt_factor = 1.0_rt / (std::sqrt(BIAS2 * safe_error) + ADDON);
dt_factor = amrex::max(dt_factor, ETAMIN);
if (n_consecutive_reject >= 1) {
dt_factor = amrex::min(dt_factor, ETAMXF);
}
}

dt_sub *= dt_factor;
previous_step_rejected = true;
++n_consecutive_reject;
++n_reject;

#ifndef AMREX_USE_GPU
if (integrator_rp::burner_verbose &&
(be.n_step < 40 || be.n_step % 1000 == 0)) {
std::cout << "[BE_STEP] reject n=" << be.n_step
<< " t=" << be.t
<< " dt=" << dt_attempt
<< " dt_next=" << dt_sub
<< " ierr=" << ierr
<< " have_error_estimate=" << have_error_estimate
<< " rel_error=" << rel_error
<< " max_error=" << max_error
<< " err_idx=" << err_idx
<< " dt_factor=" << dt_factor
<< " consecutive_rejects=" << n_consecutive_reject
<< std::endl;
}
#endif

}

++be.n_step;
++be.n_step;
++n_attempt;

}

if (be.n_step >= integrator_rp::ode_max_steps) {
ierr = IERR_TOO_MANY_STEPS;
}

#ifndef AMREX_USE_GPU
if (integrator_rp::burner_verbose) {
std::cout << "[BE_SUMMARY] attempts=" << n_attempt
<< " accepts=" << be.n_step
<< " rejects=" << n_reject
<< std::endl;
}
#endif

return ierr;

}
Expand Down
13 changes: 9 additions & 4 deletions integration/integrator_rhs_sdc.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@
template <typename BurnT, typename T>
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
void rhs(const amrex::Real time, BurnT& state, T& int_state, RArray1D& ydot,
[[maybe_unused]] const bool in_jacobian=false)
[[maybe_unused]] const bool in_jacobian=false, const bool do_clean_state=true)
{

// update rho
Expand All @@ -33,7 +33,9 @@ void rhs(const amrex::Real time, BurnT& state, T& int_state, RArray1D& ydot,
// ensure that the mass fractions are valid -- only int_state is
// updated here

clean_state(time, state, int_state);
if (do_clean_state) {
clean_state(time, state, int_state);
}

// convert to the burn_t -- this does an EOS call to get T
// and populates the (burn_t) state
Expand Down Expand Up @@ -86,7 +88,8 @@ void rhs(const amrex::Real time, BurnT& state, T& int_state, RArray1D& ydot,

template<typename BurnT, typename T, class MatrixType>
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
void jac (const amrex::Real time, BurnT& state, T& int_state, MatrixType& pd)
void jac (const amrex::Real time, BurnT& state, T& int_state, MatrixType& pd,
const bool do_clean_state=true)
{

// update rho, rho*u, ... in the burn_t state
Expand All @@ -96,7 +99,9 @@ void jac (const amrex::Real time, BurnT& state, T& int_state, MatrixType& pd)
// ensure that the mass fractions are valid -- only int_state is
// updated here

clean_state(time, state, int_state);
if (do_clean_state) {
clean_state(time, state, int_state);
}

// convert to the burn_t -- this does an EOS call to get T
// and populates the (burn_t) state
Expand Down
Loading
Loading