add a retry_clip_species option

integration/BackwardEuler/be_type.H

@@ -51,6 +51,8 @@ struct be_t {
     ArrayUtil::MathArray2D<jac_t, 1, int_neqs, 1, int_neqs> jac;
 
     short jacobian_type;
+    bool clip_species;
+
 };
 
 #endif

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

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<amrex::Real, 1, VODE_LMAX> el;
 

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.

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);

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),

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);
         }