Merge branch 'development' into fix_atol_sdc

Author: zingale

integration/BackwardEuler/be_integrator.H

@@ -60,9 +60,9 @@ int single_step (BurnT& state, BeT& be, const amrex::Real dt)
         if (integrator_rp::do_corrector_validation) {
 #ifdef SDC
             const amrex::Real rho_current = state.rho_orig + be.t * state.ydot_a[SRHO];
-            const amrex::Real thresh = species_failure_tolerance * rho_current;
+            const amrex::Real thresh = integrator_rp::species_failure_tolerance * rho_current;
 #else
-            const amrex::Real thresh = species_failure_tolerance;
+            const amrex::Real thresh = integrator_rp::species_failure_tolerance;
 #endif
             bool fail{};
             for (int i = 1; i <= NumSpec; ++i) {

integration/VODE/vode_dvnlsd.H

@@ -151,9 +151,9 @@ amrex::Real dvnlsd (int& NFLAG, BurnT& state, DvodeT& vstate)
             if (integrator_rp::do_corrector_validation && !integrator_rp::use_number_densities) {
 #ifdef SDC
                 const amrex::Real rho_current = state.rho_orig + vstate.tn * state.ydot_a[SRHO];
-                const amrex::Real thresh = species_failure_tolerance * rho_current;
+                const amrex::Real thresh = integrator_rp::species_failure_tolerance * rho_current;
 #else
-                const amrex::Real thresh = species_failure_tolerance;
+                const amrex::Real thresh = integrator_rp::species_failure_tolerance;
 #endif
                 bool fail{};
                 for (int i = 1; i <= NumSpec; ++i) {

integration/VODE/vode_dvstep.H

@@ -244,15 +244,15 @@ int dvstep (BurnT& state, DvodeT& vstate)
 
             // Constrain abundances such that they are not negative (within a tolerance)
             // or greater than one (within a tolerance).
-            if (vstate.y(i) < -species_failure_tolerance) {
+            if (vstate.y(i) < -integrator_rp::species_failure_tolerance) {
                 valid_update = false;
                 break;
             }
 
             // Don't enforce the condition below if
             // vstate.y contains number densities
             if (! integrator_rp::use_number_densities) {
-                if (vstate.y(i) > 1.0_rt + species_failure_tolerance) {
+                if (vstate.y(i) > 1.0_rt + integrator_rp::species_failure_tolerance) {
                     valid_update = false;
                     break;
                 }
@@ -286,12 +286,12 @@ int dvstep (BurnT& state, DvodeT& vstate)
                 break;
             }
 
-            if (vstate.y(SFS+i) < -species_failure_tolerance * rho_current) {
+            if (vstate.y(SFS+i) < -integrator_rp::species_failure_tolerance * rho_current) {
                 valid_update = false;
                 break;
             }
 
-            if (vstate.y(SFS+i) > (1.0_rt + species_failure_tolerance) * rho_current) {
+            if (vstate.y(SFS+i) > (1.0_rt + integrator_rp::species_failure_tolerance) * rho_current) {
                 valid_update = false;
                 break;
             }

integration/_parameters

@@ -117,3 +117,9 @@ nse_include_enu_weak       bool        1
 
 # for the linear algebra, do we allow pivoting?
 linalg_do_pivoting         bool        1
+
+# We will use this parameter to determine if a given species
+# abundance is unreasonably small or large (each X must satisfy
+# -failure_tolerance <= X <= 1.0 + failure_tolerance).
+species_failure_tolerance  real       1.e-2_rt
+

integration/integrator_data.H

@@ -9,10 +9,6 @@
 
 constexpr int INT_NEQS = NumSpec + 1;
 
-// We will use this parameter to determine if a given species
-// abundance is unreasonably small or large (each X must satisfy
-// -failure_tolerance <= X <= 1.0 + failure_tolerance).
-constexpr amrex::Real species_failure_tolerance = 1.e-2_rt;
 
 // for outputting the burn state failure, use enough precision
 // to be reproducible (not guaranteed for std::fixed)

integration/integrator_setup_sdc.H

@@ -186,11 +186,11 @@ void integrator_cleanup (IntegratorT& int_state, BurnT& state,
     }
 
     for (int n = 0; n < NumSpec; ++n) {
-        if (state.y[SFS+n] < -state.rho * species_failure_tolerance) {
+        if (state.y[SFS+n] < -state.rho * integrator_rp::species_failure_tolerance) {
             state.success = false;
         }
 
-        if (state.y[SFS+n] > state.rho * (1.0_rt + species_failure_tolerance)) {
+        if (state.y[SFS+n] > state.rho * (1.0_rt + integrator_rp::species_failure_tolerance)) {
             state.success = false;
         }
     }

integration/integrator_setup_strang.H

@@ -155,14 +155,14 @@ void integrator_cleanup (IntegratorT& int_state, BurnT& state,
     }
 
     for (int n = 1; n <= NumSpec; ++n) {
-        if (int_state.y(n) < -species_failure_tolerance) {
+        if (int_state.y(n) < -integrator_rp::species_failure_tolerance) {
             state.success = false;
         }
 
         // Don't enforce a max if we are evolving number densities
 
         if (! integrator_rp::use_number_densities) {
-            if (int_state.y(n) > 1.0_rt + species_failure_tolerance) {
+            if (int_state.y(n) > 1.0_rt + integrator_rp::species_failure_tolerance) {
                 state.success = false;
             }
         }