From 7d5a6634909980575de0e166a40857dee3caeff5 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Wed, 8 May 2024 08:38:28 -0400 Subject: [PATCH 1/2] add NSE bailout to RKC also clean up some VODE stuff --- integration/RKC/rkc.H | 37 +++++++++++++++++++++++++++++++++++ integration/VODE/vode_type.H | 20 ++++++++----------- integration/integrator_data.H | 9 +++++++-- 3 files changed, 52 insertions(+), 14 deletions(-) diff --git a/integration/RKC/rkc.H b/integration/RKC/rkc.H index d81c794384..2ffd0fce4f 100644 --- a/integration/RKC/rkc.H +++ b/integration/RKC/rkc.H @@ -15,6 +15,14 @@ #include #include +#ifdef NSE_TABLE +#include +#endif +#ifdef NSE_NET +#include +#endif + + template 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.NST > MIN_NSE_BAILOUT_STEPS && rstate.t <= vstate.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++; diff --git a/integration/VODE/vode_type.H b/integration/VODE/vode_type.H index cab1171cfc..9878e8d68c 100644 --- a/integration/VODE/vode_type.H +++ b/integration/VODE/vode_type.H @@ -9,30 +9,26 @@ #include -const amrex::Real UROUND = std::numeric_limits::epsilon(); +constexpr amrex::Real UROUND = std::numeric_limits::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 diff --git a/integration/integrator_data.H b/integration/integrator_data.H index 974dbf482d..1988d39a80 100644 --- a/integration/integrator_data.H +++ b/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 { IERR_SUCCESS = 1, From 36da3f63a441f9ca2ef34950243d5ed6d238134d Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Wed, 8 May 2024 09:23:57 -0400 Subject: [PATCH 2/2] fix NSE --- integration/RKC/rkc.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/integration/RKC/rkc.H b/integration/RKC/rkc.H index 2ffd0fce4f..1482a752d6 100644 --- a/integration/RKC/rkc.H +++ b/integration/RKC/rkc.H @@ -246,7 +246,7 @@ 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.NST > MIN_NSE_BAILOUT_STEPS && rstate.t <= vstate.tout) { + if (rstate.nsteps > MIN_NSE_BAILOUT_STEPS && rstate.t <= rstate.tout) { // first we need to make the burn_t in sync #ifdef STRANG