diff --git a/integration/BackwardEuler/be_type.H b/integration/BackwardEuler/be_type.H index 04d7d8f19..f7b5bbc67 100644 --- a/integration/BackwardEuler/be_type.H +++ b/integration/BackwardEuler/be_type.H @@ -51,6 +51,8 @@ struct be_t { ArrayUtil::MathArray2D jac; short jacobian_type; + bool clip_species; + }; #endif diff --git a/integration/RKC/rkc_type.H b/integration/RKC/rkc_type.H index f4a20a0cd..48cc3ec79 100644 --- a/integration/RKC/rkc_type.H +++ b/integration/RKC/rkc_type.H @@ -98,6 +98,9 @@ struct rkc_t { // not used here, but needed for compatibility with other integrators short jacobian_type; + // do we clip the species before evaluating the rates + bool clip_species; + }; #ifdef SDC diff --git a/integration/VODE/vode_type.H b/integration/VODE/vode_type.H index a1df1f9bf..505119b1c 100644 --- a/integration/VODE/vode_type.H +++ b/integration/VODE/vode_type.H @@ -121,6 +121,9 @@ struct dvode_t // jacobian_type = the type of Jacobian to use (1 = analytic, 2 = numerical) short jacobian_type; + // do we clip the species before evaluating the rates + bool clip_species; + // EL = Real array of integration coefficients. See DVSET amrex::Array1D el; diff --git a/integration/_parameters b/integration/_parameters index 56a4a5522..a3a9d149d 100644 --- a/integration/_parameters +++ b/integration/_parameters @@ -73,6 +73,9 @@ use_burn_retry bool 0 # a retry? retry_swap_jacobian bool 1 +# when we retry, do we clip the species (or if we were, stop clipping) +retry_swap_clipping bool 0 + # Tolerances for the solver (relative and absolute), for the # species and energy equations. If set to < 0, then the same # value as the first attempt is used. diff --git a/integration/integrator_setup_sdc.H b/integration/integrator_setup_sdc.H index 2eac92bef..16f8ced70 100644 --- a/integration/integrator_setup_sdc.H +++ b/integration/integrator_setup_sdc.H @@ -50,6 +50,13 @@ IntegratorT integrator_setup (BurnT& state, amrex::Real dt, bool is_retry) int_state.jacobian_type = integrator_rp::jacobian; } + // set whether we clip species + if (is_retry && integrator_rp::retry_swap_clipping) { + int_state.clip_species = (integrator_rp::do_species_clip) ? false : true; + } else { + int_state.clip_species = integrator_rp::do_species_clip; + } + // Fill in the initial integration state. burn_to_int(state, int_state); diff --git a/integration/integrator_type_sdc.H b/integration/integrator_type_sdc.H index d0c699b50..6f6bbbd3b 100644 --- a/integration/integrator_type_sdc.H +++ b/integration/integrator_type_sdc.H @@ -18,7 +18,7 @@ void clean_state(const amrex::Real time, BurnT& state, T& int_state) { // Ensure that mass fractions always stay positive. - if (integrator_rp::do_species_clip) { + if (int_state.clip_species) { for (int n = 1; n <= NumSpec; ++n) { // we use 1-based indexing, so we need to offset SFS int_state.y(SFS+n) = amrex::Clamp(int_state.y(SFS+n), diff --git a/integration/integrator_type_strang.H b/integration/integrator_type_strang.H index 7eb3c6d60..7e71f7279 100644 --- a/integration/integrator_type_strang.H +++ b/integration/integrator_type_strang.H @@ -53,7 +53,7 @@ void clean_state (const amrex::Real time, BurnT& state, I& int_state) // Ensure that mass fractions always stay positive and less than or // equal to 1. - if (integrator_rp::do_species_clip) { + if (int_state.clip_species) { for (int n = 1; n <= NumSpec; ++n) { int_state.y(n) = amrex::Clamp(int_state.y(n), integrator_rp::SMALL_X_SAFE, 1.0_rt); }