diff --git a/.gitignore b/.gitignore index 250e33146..59044a53e 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,4 @@ + # latex intermediates *.aux *.log diff --git a/Docs/source/ode_integrators.rst b/Docs/source/ode_integrators.rst index 0399416df..1aa3a0765 100644 --- a/Docs/source/ode_integrators.rst +++ b/Docs/source/ode_integrators.rst @@ -72,13 +72,13 @@ Presently, allowed integrators are: * ``4`` : ROS2 method, a 2nd order, L-stable method :cite:`ros2`. - * ``5`` : a first-order method based on the YASS method (described in :cite:`YASS`). + * ``5`` : Rosenbrock-Euler, a first-order method. Here the "P" suffix refers to methods developed to satisfy the stiff accuracy conditions of :cite:`Prothero1974` (ROS2S also satisfies these). - By default, the H211b error-history timestep controller from :cite:`h211b` (see Eq. 31) is used, but the YASS heuristic method can be used instead by setting ``integrator.rosenbrock_timestep_controller=1``. + The H211b error-history timestep controller from :cite:`h211b` (see Eq. 31) is used. * ``VODE``: the VODE :cite:`vode` integration package. We ported this integrator to C++ and removed the non-stiff integration code paths. diff --git a/integration/Rosenbrock/_parameters b/integration/Rosenbrock/_parameters index 4bd9289ef..680d7962f 100644 --- a/integration/Rosenbrock/_parameters +++ b/integration/Rosenbrock/_parameters @@ -13,24 +13,9 @@ X_reject_buffer real 1.0 # 5 = Rosenbrock-Euler 1-stage method rosenbrock_tableau int 0 -# Timestep controller selector: -# 0 = H211b error-history controller -# 1 = Khokhlov/YASS concentration-change controller -rosenbrock_timestep_controller int 0 - # H211b timestep controller parameters h211b_b real 4.0 h211b_k real 2.5 h211b_fac_min real 0.2 h211b_fac_max real 6.0 h211b_reduction_fac real 0.5 - -# Khokhlov/YASS timestep controller parameters. This controller accepts a -# step when each active species changes by no more than yass_epsilon -# relative to its starting value. Species with concentration <= yass_floor -# times the most abundant species do not constrain the step. -yass_epsilon real 0.1 -yass_floor real 1.e-3 -yass_fac_min real 0.2 -yass_fac_max real 6.0 -yass_reduction_fac real 0.5 diff --git a/integration/Rosenbrock/rosenbrock_integrator.H b/integration/Rosenbrock/rosenbrock_integrator.H index e7104c2c6..6c8311808 100644 --- a/integration/Rosenbrock/rosenbrock_integrator.H +++ b/integration/Rosenbrock/rosenbrock_integrator.H @@ -171,29 +171,6 @@ bool species_change_valid (const rosenbrock_t& rstate) return true; } -template -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -amrex::Real yass_change_norm (const rosenbrock_t& rstate) -{ - amrex::Real max_change = 0.0_rt; - amrex::Real max_species = 0.0_rt; - - for (int i = 1; i <= NumSpec && i <= int_neqs; ++i) { - max_species = amrex::max(max_species, std::abs(rstate.y(i))); - } - - const amrex::Real active_floor = integrator_rp::yass_floor * max_species; - - for (int i = 1; i <= NumSpec && i <= int_neqs; ++i) { - const amrex::Real yi = std::abs(rstate.y(i)); - if (yi > active_floor) { - max_change = amrex::max(max_change, std::abs(rstate.ynew(i) - rstate.y(i)) / yi); - } - } - - return max_change / integrator_rp::yass_epsilon; -} - template AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool matrix_is_finite (const MatrixType& matrix) @@ -381,15 +358,6 @@ int rosenbrock_integrator (BurnT& state, rosenbrock_t()>& } } - if (integrator_rp::rosenbrock_timestep_controller == 1 && - (integrator_rp::yass_epsilon <= 0.0_rt || - integrator_rp::yass_floor < 0.0_rt || - integrator_rp::yass_fac_min <= 0.0_rt || - integrator_rp::yass_fac_max < integrator_rp::yass_fac_min || - integrator_rp::yass_reduction_fac <= 0.0_rt)) { - return IERR_BAD_INPUTS; - } - rstate.n_rhs = 0; rstate.n_jac = 0; rstate.n_step = 0; @@ -564,22 +532,13 @@ int rosenbrock_integrator (BurnT& state, rosenbrock_t()>& const bool valid_update = rosenbrock::species_change_valid(rstate); constexpr amrex::Real err_min = 1.e-10_rt; amrex::Real err = rosenbrock::error_norm(rstate); - amrex::Real controller_fac = 1.0_rt; - - if (integrator_rp::rosenbrock_timestep_controller == 1) { - err = rosenbrock::yass_change_norm(rstate); - controller_fac = amrex::Clamp(1.0_rt / amrex::max(err, err_min), - integrator_rp::yass_fac_min, - integrator_rp::yass_fac_max); - } else { - controller_fac = amrex::Clamp( - std::pow(1.0_rt / amrex::max(err, err_min), - 1.0_rt / (integrator_rp::h211b_b * integrator_rp::h211b_k)) * - std::pow(1.0_rt / amrex::max(errold, err_min), - 1.0_rt / (integrator_rp::h211b_b * integrator_rp::h211b_k)) * - std::pow(facold, -1.0_rt / integrator_rp::h211b_b), - integrator_rp::h211b_fac_min, integrator_rp::h211b_fac_max); - } + amrex::Real controller_fac = amrex::Clamp( + std::pow(1.0_rt / amrex::max(err, err_min), + 1.0_rt / (integrator_rp::h211b_b * integrator_rp::h211b_k)) * + std::pow(1.0_rt / amrex::max(errold, err_min), + 1.0_rt / (integrator_rp::h211b_b * integrator_rp::h211b_k)) * + std::pow(facold, -1.0_rt / integrator_rp::h211b_b), + integrator_rp::h211b_fac_min, integrator_rp::h211b_fac_max); facold = controller_fac; errold = err; @@ -613,14 +572,11 @@ int rosenbrock_integrator (BurnT& state, rosenbrock_t()>& n_reject += 1; last = false; if (n_reject >= 2) { - hnew *= integrator_rp::rosenbrock_timestep_controller == 1 ? - integrator_rp::yass_reduction_fac : integrator_rp::h211b_reduction_fac; + hnew *= integrator_rp::h211b_reduction_fac; } if (valid_update) { - const amrex::Real reject_fac = - integrator_rp::rosenbrock_timestep_controller == 1 ? - integrator_rp::yass_reduction_fac : integrator_rp::h211b_reduction_fac; - h = posneg * amrex::min(std::abs(hnew), reject_fac * std::abs(h)); + h = posneg * amrex::min(std::abs(hnew), + integrator_rp::h211b_reduction_fac * std::abs(h)); } else { h *= 0.25_rt; }