Merge branch 'development' into remove_divide_in_screening

Author: zhichen3

EOS/helmholtz/actual_eos.H

@@ -1315,10 +1315,11 @@ void actual_eos_init ()
     for (int j = 0; j < jmax; ++j) {
         amrex::Real tsav = tlo + j * tstp;
         t[j] = std::pow(10.0e0_rt, tsav);
-        for (int i = 0; i < imax; ++i) {
-            amrex::Real dsav = dlo + i * dstp;
-            d[i] = std::pow(10.0e0_rt, dsav);
-        }
+    }
+
+    for (int i = 0; i < imax; ++i) {
+        amrex::Real dsav = dlo + i * dstp;
+        d[i] = std::pow(10.0e0_rt, dsav);
     }
 
     // it does not work on all machines (for GPUs) broadcast to other

integration/BackwardEuler/be_integrator.H

@@ -57,6 +57,27 @@ int single_step (BurnT& state, BeT& be, const amrex::Real dt)
 
         // work with the current guess
 
+        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;
+#else
+            const amrex::Real thresh = species_failure_tolerance;
+#endif
+            bool fail{};
+            for (int i = 1; i <= NumSpec; ++i) {
+                if (be.y(i) < -thresh) {
+                    fail = true;
+                    break;
+                }
+            }
+
+            if (fail) {
+                ierr = IERR_BAD_STATE_IN_CORRECTOR;
+                break;
+            }
+        }
+
         // get the ydots for our current guess of y
 
         rhs(be.t, state, be, ydot);

integration/integrator_data.H

@@ -32,6 +32,7 @@ enum integrator_errors : std::int8_t {
     IERR_TOO_MUCH_ACCURACY_REQUESTED = -5,
     IERR_CORRECTOR_CONVERGENCE = -6,
     IERR_LU_DECOMPOSITION_ERROR = -7,
+    IERR_BAD_STATE_IN_CORRECTOR = -8,
     IERR_ENTERED_NSE = -100
 };
 

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 * 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 + species_failure_tolerance)) {
             state.success = false;
         }
     }