diff --git a/integration/BackwardEuler/be_integrator.H b/integration/BackwardEuler/be_integrator.H index 8e9ea227e..da0980f57 100644 --- a/integration/BackwardEuler/be_integrator.H +++ b/integration/BackwardEuler/be_integrator.H @@ -28,19 +28,25 @@ /// template 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(); 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 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; @@ -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; @@ -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; @@ -183,6 +190,8 @@ int single_step (BurnT& state, BeT& be, const amrex::Real dt) } + be.t = time_save; + return ierr; } @@ -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 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 @@ -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 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 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 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) { @@ -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 { @@ -317,13 +376,52 @@ 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; } @@ -331,6 +429,15 @@ int be_integrator (BurnT& state, BeT& be) 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; } diff --git a/integration/integrator_rhs_sdc.H b/integration/integrator_rhs_sdc.H index dead4a654..7d5a1895f 100644 --- a/integration/integrator_rhs_sdc.H +++ b/integration/integrator_rhs_sdc.H @@ -23,7 +23,7 @@ template 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 @@ -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 @@ -86,7 +88,8 @@ void rhs(const amrex::Real time, BurnT& state, T& int_state, RArray1D& ydot, template 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 @@ -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 diff --git a/integration/integrator_rhs_strang.H b/integration/integrator_rhs_strang.H index 90b35a167..5a17471bc 100644 --- a/integration/integrator_rhs_strang.H +++ b/integration/integrator_rhs_strang.H @@ -22,7 +22,8 @@ template 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) +void rhs (const amrex::Real time, BurnT& state, T& int_state, RArray1D& ydot, + [[maybe_unused]] const bool in_jacobian=false, const bool do_clean_state=true) { // We are integrating a system of @@ -38,7 +39,9 @@ void rhs (const amrex::Real time, BurnT& state, T& int_state, RArray1D& ydot, [[ // Fix the state as necessary -- this ensures that the mass // fractions that enter are valid (and optionally normalized) - clean_state(time, state, int_state); + if (do_clean_state) { + clean_state(time, state, int_state); + } // Update the thermodynamics as necessary -- this primarily takes // the information in int_state (X and e), copies it to the @@ -112,7 +115,8 @@ void rhs (const amrex::Real time, BurnT& state, T& int_state, RArray1D& ydot, [[ // Analytical Jacobian template AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void jac ([[maybe_unused]] const amrex::Real time, BurnT& state, T& int_state, MatrixType& pd) +void jac ([[maybe_unused]] const amrex::Real time, BurnT& state, T& int_state, MatrixType& pd, + [[maybe_unused]] const bool do_clean_state=true) { // Only do the burn if the incoming temperature is within the temperature diff --git a/integration/utils/initial_timestep.H b/integration/utils/initial_timestep.H index c53fa026f..4fb206ee1 100644 --- a/integration/utils/initial_timestep.H +++ b/integration/utils/initial_timestep.H @@ -62,7 +62,8 @@ amrex::Real initial_react_dt (BurnT& burn_state, IntT& int_state, // Call the RHS, then estimate the finite difference. - rhs(int_state.t, burn_state, int_state, ydot_temp); + rhs(int_state.t, burn_state, int_state, ydot_temp, false, + ! integrator_rp::use_number_densities); for (int ii = 1; ii <= INT_NEQS; ii++) { ddydtt(ii) = (ydot_temp(ii) - ydot(ii)) / h; @@ -89,7 +90,7 @@ amrex::Real initial_react_dt (BurnT& burn_state, IntT& int_state, // Save the final timestep, with a bias factor. amrex::Real dt = h / 2.0_rt; - dt = amrex::Clamp(h, hL, hU); + dt = amrex::Clamp(dt, hL, hU); dt = amrex::min(dt, integrator_rp::ode_max_dt); diff --git a/unit_test/burn_cell/inputs_vode_example b/unit_test/burn_cell/inputs_vode_example index d4c895b08..6b33307ac 100644 --- a/unit_test/burn_cell/inputs_vode_example +++ b/unit_test/burn_cell/inputs_vode_example @@ -15,9 +15,10 @@ unit_test.X2 = 0.0e0 unit_test.X3 = 0.0e0 -integrator.rtol_spec = 1.0e-8 -integrator.rtol_enuc = 1.0e-8 +integrator.rtol_spec = 1.0e-3 +integrator.rtol_enuc = 1.0e-3 integrator.atol_spec = 1.0e-14 integrator.atol_enuc = 1.0e-8 +integrator.use_number_densities = 1