add NSE bailout to Rosenbrock integrator

integration/Rosenbrock/rosenbrock_integrator.H

@@ -16,6 +16,12 @@
 #ifdef SDC
 #include <integrator_rhs_sdc.H>
 #endif
+#ifdef NSE_TABLE
+#include <nse_table_check.H>
+#endif
+#ifdef NSE_NET
+#include <nse_check.H>
+#endif
 
 #include <integrator_data.H>
 #include <rosenbrock_type.H>
@@ -389,7 +395,14 @@ int rosenbrock_integrator (BurnT& state, rosenbrock_t<integrator_neqs<BurnT>()>&
     amrex::Real x = rstate.t;
     amrex::Array1D<amrex::Real, 1, int_neqs> rhs_tmp;
 
+    // main evolution loop
+
+    // within this loop, x will be the current time and h will be the
+    // step size we are attempting.  The solution at the current time
+    // is rstate.y().
+
     while (true) {
+
         if (rstate.n_step > integrator_rp::ode_max_steps) {
             rstate.t = x;
             rstate.dt = h;
@@ -414,6 +427,32 @@ int rosenbrock_integrator (BurnT& state, rosenbrock_t<integrator_neqs<BurnT>()>&
             last = true;
         }
 
+#ifdef NSE
+        // check if, during the course of integration, we hit NSE, and
+        // if so, bail out.
+
+        // we only do this after MIN_NSE_BAILOUT_STEPS to prevent us
+        // from hitting this right at the start when the integrator
+        // might do so wild exploration.  Also ensure we are not
+        // working > tmax, so we don't need to worry about
+        // extrapolating back in time.
+
+        if (rstate.n_step > MIN_NSE_BAILOUT_STEPS && x <= rstate.tout) {
+            // first we need to make the burn_t in sync
+
+            int_to_burn(x, rstate, state);
+
+            if (in_nse(state)) {
+                rstate.t = x;
+                return IERR_ENTERED_NSE;
+            }
+        }
+#endif
+
+
+        // take a step -- this loop will keep trying until a step is
+        // successful or a catastrophic error is encountered.
+
         while (true) {
             rstate.dt = h;
 
@@ -527,6 +566,9 @@ int rosenbrock_integrator (BurnT& state, rosenbrock_t<integrator_neqs<BurnT>()>&
 
             amrex::Real hnew = h * controller_fac;
 
+            // update the starting time and state for the next step
+            // and predict the next timestep
+
             if (err <= 1.0_rt && valid_update) {
                 for (int n = 1; n <= int_neqs; ++n) {
                     rstate.y(n) = rstate.ynew(n);
@@ -545,6 +587,9 @@ int rosenbrock_integrator (BurnT& state, rosenbrock_t<integrator_neqs<BurnT>()>&
                 break;
             }
 
+            // if we made it here, then we've failed
+            // cut the timestep in prep for trying again
+
             reject = true;
             n_reject += 1;
             last = false;