diff --git a/Docs/source/ode_integrators.rst b/Docs/source/ode_integrators.rst index 4644a347a..26cb3a366 100644 --- a/Docs/source/ode_integrators.rst +++ b/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 diff --git a/integration/VODE/vode_dvnlsd.H b/integration/VODE/vode_dvnlsd.H index 46ede5d86..0af50f9ec 100644 --- a/integration/VODE/vode_dvnlsd.H +++ b/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. diff --git a/integration/_parameters b/integration/_parameters index a11b5f3ea..f6fb6a41f 100644 --- a/integration/_parameters +++ b/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