VODE: bail out of the nonlinear solve if we predict negative mass fractions

Docs/source/ode_integrators.rst

@@ -248,7 +248,7 @@ fractions.  Looser than this can produce large errors.
 Controlling Species $\sum_k X_k = 1$
 ====================================
 
-.. index:: integrator.renormalize_abundances, integrator.SMALL_X_SAFE, integrator.do_species_clip
+.. index:: integrator.renormalize_abundances, integrator.SMALL_X_SAFE, integrator.do_species_clip, integrator.do_corrector_validation
 
 The ODE integrators don't know about the constraint that
 
@@ -275,10 +275,18 @@ constraint on the intermediate states during the integration.
   The default is ``1.e-30``.
 
 * ``integrator.do_species_clip`` : this enforces that the mass fractions
-  all in $[\mathtt{SMALL\_X\_SAFE}, 1.0]$.
-
-  This is enabled by default.
-
+  all in $[\mathtt{SMALL\_X\_SAFE}, 1.0]$ before calling the network righthand
+  side function.
+
+  This is off by default.  Turning this on can sometimes make the integrator
+  work a lot harder.
+
+* ``integrator.do_corrector_validation`` : in the nonlinear solve
+  corrector loop, when we get a corrected integration state, do we
+  check to make sure the mass fractions are valid before calling the
+  righthand function?  This is needed in some cases if
+  ``integrator.do_species_clip`` is disabled.  Note: this is not
+  implemented for every integrator.
 
 
 Retry Mechanism

integration/VODE/vode_dvnlsd.H

@@ -145,6 +145,32 @@ amrex::Real dvnlsd (int& NFLAG, BurnT& state, DvodeT& vstate)
                 vstate.y(i) = vstate.yh(i,1) + vstate.acor(i);
             }
 
+            // sometime VODE goes way off tangent.  If these mass fractions
+            // are really bad, then let's just bail now
+
+            if (integrator_rp::do_corrector_validation) {
+#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;
+#else
+                const amrex::Real thresh = species_failure_tolerance;
+#endif
+                bool fail{};
+                for (int i = 1; i <= NumSpec; ++i) {
+                    if (vstate.y(i) < -thresh) {
+                        fail = true;
+                        break;
+                    }
+                }
+
+                if (fail) {
+                    NFLAG = -1;
+                    vstate.ICF = 2;
+                    vstate.IPUP = 1;
+                    return ACNRM;
+                }
+            }
+
             // Test for convergence. If M > 0, an estimate of the convergence
             // rate constant is stored in CRATE, and this is used in the test.
 

integration/_parameters

@@ -91,7 +91,11 @@ retry_atol_enuc                real      -1
 
 # in the clean_state process, do we clip the species such that they
 # are in [0, 1]?
-do_species_clip              bool         1
+do_species_clip              bool         0
+
+# in the corrector loop, do we check if the predicted state is
+# valid (X > 0) before calling the RHS?
+do_corrector_validation      bool         1
 
 # flag for turning on the use of number densities for all species
 use_number_densities     bool     0