add NSE bailout to RKC (#1544)

also clean up some VODE stuff

integration/RKC/rkc.H

@@ -15,6 +15,14 @@
 #include <circle_theorem.H>
 #include <integrator_data.H>
 
+#ifdef NSE_TABLE
+#include <nse_table_check.H>
+#endif
+#ifdef NSE_NET
+#include <nse_check.H>
+#endif
+
+
 template <typename BurnT, typename RkcT>
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
 void step (BurnT& state, RkcT& rstate, const amrex::Real h, const int m)
@@ -225,6 +233,35 @@ int rkclow (BurnT& state, RkcT& rstate)
 
         err = std::sqrt(err / int_neqs);
 
+        // before we accept or reject the step, let's check if we've entered
+        // NSE
+
+#ifdef NSE
+       // check if, during the course of integration, we hit NSE, and
+       // if so, bail out we rely on the state being consistent after
+       // the call to dvstep, even if the step failed.
+
+       // we only do this after MIN_NSE_BAILOUT_STEPS to prevent us
+       // from hitting this right at the start when VODE 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.nsteps > MIN_NSE_BAILOUT_STEPS && rstate.t <= rstate.tout) {
+           // first we need to make the burn_t in sync
+
+#ifdef STRANG
+           update_thermodynamics(state, rstate);
+#endif
+#ifdef SDC
+           int_to_burn(rstate.t, rstate, state);
+#endif
+
+           if (in_nse(state)) {
+               return IERR_ENTERED_NSE;
+           }
+       }
+#endif
+
         if (err > 1.0_rt) {
             // Step is rejected.
             rstate.nrejct++;

integration/VODE/vode_type.H

@@ -9,30 +9,26 @@
 
 #include <integrator_data.H>
 
-const amrex::Real UROUND = std::numeric_limits<amrex::Real>::epsilon();
+constexpr amrex::Real UROUND = std::numeric_limits<amrex::Real>::epsilon();
 
 // CCMXJ  = Threshold on DRC for updating the Jacobian
-const amrex::Real CCMXJ = 0.2e0_rt;
+constexpr amrex::Real CCMXJ = 0.2e0_rt;
 
-const amrex::Real HMIN = 0.0_rt;
+constexpr amrex::Real HMIN = 0.0_rt;
 
 // For each species whose abundance is above a certain threshold, we do
 // not allow its mass fraction to change by more than a certain amount
 // in any integration step.
-const amrex::Real vode_increase_change_factor = 4.0_rt;
-const amrex::Real vode_decrease_change_factor = 0.25_rt;
+constexpr amrex::Real vode_increase_change_factor = 4.0_rt;
+constexpr amrex::Real vode_decrease_change_factor = 0.25_rt;
 
 // For the backward differentiation formula (BDF) integration
 // the maximum order should be no greater than 5.
-const int VODE_MAXORD = 5;
-const int VODE_LMAX = VODE_MAXORD + 1;
+constexpr int VODE_MAXORD = 5;
+constexpr int VODE_LMAX = VODE_MAXORD + 1;
 
 // How many timesteps should pass before refreshing the Jacobian
-const int max_steps_between_jacobian_evals = 50;
-
-#ifdef NSE
-const int MIN_NSE_BAILOUT_STEPS = 10;
-#endif
+constexpr int max_steps_between_jacobian_evals = 50;
 
 // Type dvode_t contains the integration solution and control variables
 template<int int_neqs>

integration/integrator_data.H

@@ -5,12 +5,17 @@
 
 // Define the size of the ODE system that VODE will integrate
 
-const int INT_NEQS = NumSpec + 1;
+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).
-const amrex::Real species_failure_tolerance = 1.e-2_rt;
+constexpr amrex::Real species_failure_tolerance = 1.e-2_rt;
+
+#ifdef NSE
+constexpr int MIN_NSE_BAILOUT_STEPS = 10;
+#endif
+
 
 enum integrator_errors : std::int8_t {
     IERR_SUCCESS = 1,