add check interval for checking nse inside integrator

Docs/source/nse_net.rst

@@ -354,3 +354,7 @@ to the self-consistent nse check:
   the subsequent NSE checks. This is mainly to avoid unnecessary computations
   of computing the NSE mass fractions when the current temperature is too low.
   This is set to 4.0e9 by default.
+
+* ``nse.nse_check_interval`` is the interval at which we perform the
+  nse check in the integrator. This is used to minimize the
+  expensive check overhead. This is set to 5 by default.

integration/RKC/rkc.H

@@ -246,15 +246,17 @@ int rkclow (BurnT& state, RkcT& rstate)
        // wild exploration.  Also ensure we are not working > tmax,
        // so we don't need to worry about extrapolating back in time.
 
-       if (rstate.nsteps > MIN_NSE_BAILOUT_STEPS && rstate.t <= rstate.tout) {
-           // first we need to make the burn_t in sync
+       // For NSE_NET, we only check for some interval to prevent
+       // expensive check overhead
 
-#ifdef STRANG
-           update_thermodynamics(state, rstate);
+       if (rstate.nsteps > MIN_NSE_BAILOUT_STEPS && rstate.t <= rstate.tout
+#ifdef NSE_NET
+           && (vstate.n_step % nse_check_interval == 0 || state.T > 6.e9_rt)
 #endif
-#ifdef SDC
+           ) {
+           // first we need to make the burn_t in sync
+
            int_to_burn(rstate.t, rstate, state);
-#endif
 
            if (in_nse(state)) {
                return IERR_ENTERED_NSE;

integration/VODE/vode_dvode.H

@@ -233,7 +233,14 @@ int dvode (BurnT& state, DvodeT& vstate)
        // wild exploration.  Also ensure we are not working > tmax,
        // so we don't need to worry about extrapolating back in time.
 
-       if (vstate.n_step > MIN_NSE_BAILOUT_STEPS && vstate.tn <= vstate.tout) {
+       // For NSE_NET, we only check for some interval to prevent
+       // expensive check overhead
+
+       if (vstate.n_step > MIN_NSE_BAILOUT_STEPS && vstate.tn <= vstate.tout
+#ifdef NSE_NET
+           && (vstate.n_step % nse_check_interval == 0 || state.T > 6.e9_rt)
+#endif
+           ) {
            // first we need to make the burn_t in sync
 
            int_to_burn(vstate.tn, vstate, state);

interfaces/burner.H

@@ -50,23 +50,27 @@ void burner (BurnT& state, amrex::Real dt)
         // left after the failure
         amrex::Real dt_remaining = amrex::max(dt - state.time, 0.0_rt);
 
-        // we use a relaxed NSE criteria now to catch states that are
-        // right on the edge of being in NSE
+        if (!state.success && dt_remaining > 0.0_rt) {
+
+            // Proceed with NSE burn if NSE is detected within integrator
+            // If we allow a relaxed NSE criteria to catch states that are
+            // right on the edge of being in NSE, then evaluate that here.
 #ifdef NSE_TABLE
-        if (in_nse(state, true) && !state.success && dt_remaining > 0.0) {
+            if (state.error_code == IERR_ENTERED_NSE || in_nse(state, true)) {
 #else
-        if (in_nse(state, nse_skip_molar) && !state.success && dt_remaining > 0.0) {
+            if (state.error_code == IERR_ENTERED_NSE || (nse_skip_molar && in_nse(state, nse_skip_molar))) {
 #endif
 
 #ifndef AMREX_USE_GPU
-            std::cout << "recovering burn failure in NSE, zone = (" << state.i << ", " << state.j << ", " << state.k << ")" << std::endl;
+                std::cout << "recovering burn failure in NSE, zone = (" << state.i << ", " << state.j << ", " << state.k << ")" << std::endl;
 #endif
 
-            // This will append to state.e the amount additional
-            // energy released from adjusting to the new NSE state
+                // This will append to state.e the amount additional
+                // energy released from adjusting to the new NSE state
 #ifdef SDC
-            sdc_nse_burn(state, dt_remaining);
+                sdc_nse_burn(state, dt_remaining);
 #endif
+            }
         }
     }
 #else

nse_solver/_parameters

@@ -43,3 +43,6 @@ T_nse_net               real       -1.0
 
 # Minimum Temperature required for NSE
 T_min_nse               real       4.0e9
+
+# NSE check interval in integrator
+nse_check_interval      int        5

nse_solver/nse_check.H

@@ -555,15 +555,6 @@ bool in_nse(burn_t& current_state, bool skip_molar_check=false) {
         return current_state.nse;
     }
 
-    // Check if any of the input rhoX is negative
-    // This indicates a burn failure, so we should return false
-
-    for (int n = 0; n < NumSpec; ++n) {
-        if (current_state.y[SFS+n] < 0.0_rt) {
-            return current_state.nse;
-        }
-    }
-
     // Get the NSE state using (rho, T, Ye) as input
     // Then compare input molar fractions to the NSE molar fractions.