[VODE] Allow small margin for negativity caused by interpolation (#2017)

Author: BenWibking

integration/VODE/actual_integrator.H

@@ -22,7 +22,8 @@ void actual_integrator (BurnT& state, amrex::Real dt, bool is_retry=false)
 
     auto istate = dvode(state, vode_state);
 
-    integrator_cleanup(vode_state, state, istate, state_save, dt);
+    integrator_cleanup(vode_state, state, istate, state_save, dt,
+                       vode_final_state_species_failure_tolerance_factor);
 
 }
 

integration/VODE/actual_integrator_sdc.H

@@ -27,7 +27,8 @@ void actual_integrator (BurnT& state, amrex::Real dt, bool is_retry=false)
     auto istate = dvode(state, vode_state);
     state.error_code = istate;
 
-    integrator_cleanup(vode_state, state, istate, state_save, dt);
+    integrator_cleanup(vode_state, state, istate, state_save, dt,
+                       vode_final_state_species_failure_tolerance_factor);
 
 
 }

integration/VODE/vode_type.H

@@ -22,6 +22,10 @@ constexpr amrex::Real HMIN = 0.0_rt;
 constexpr amrex::Real vode_increase_change_factor = 4.0_rt;
 constexpr amrex::Real vode_decrease_change_factor = 0.25_rt;
 
+// The interpolation back to tout is not monotonic, so the final state
+// allows slightly more negativity than internal integration nodes.
+constexpr amrex::Real vode_final_state_species_failure_tolerance_factor = 1.5_rt;
+
 // For the backward differentiation formula (BDF) integration
 // the maximum order should be no greater than 5.
 constexpr int VODE_MAXORD = 5;

integration/integrator_setup_sdc.H

@@ -135,7 +135,8 @@ state_backup_t integrator_backup (const BurnT& state) {
 template <typename BurnT, typename IntegratorT>
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
 void integrator_cleanup (IntegratorT& int_state, BurnT& state,
-                         int istate, const state_backup_t& state_save, amrex::Real dt)
+                         int istate, const state_backup_t& state_save, amrex::Real dt,
+                         const amrex::Real final_state_species_failure_tolerance_factor = 1.0_rt)
 {
 
     // Copy the integration data back to the burn state.
@@ -185,8 +186,11 @@ void integrator_cleanup (IntegratorT& int_state, BurnT& state,
         state.success = false;
     }
 
+    const amrex::Real final_state_species_failure_tolerance =
+        final_state_species_failure_tolerance_factor * state.rho * integrator_rp::species_failure_tolerance;
+
     for (int n = 0; n < NumSpec; ++n) {
-        if (state.y[SFS+n] < -state.rho * integrator_rp::species_failure_tolerance) {
+        if (state.y[SFS+n] < -final_state_species_failure_tolerance) {
             state.success = false;
         }
 

integration/integrator_setup_strang.H

@@ -119,7 +119,8 @@ state_backup_t integrator_backup (const BurnT& state) {
 template <typename BurnT, typename IntegratorT>
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
 void integrator_cleanup (IntegratorT& int_state, BurnT& state,
-                         int istate, const state_backup_t& state_save, amrex::Real dt)
+                         int istate, const state_backup_t& state_save, amrex::Real dt,
+                         const amrex::Real final_state_species_failure_tolerance_factor = 1.0_rt)
 {
 
     // Copy the integration data back to the burn state.
@@ -156,8 +157,11 @@ void integrator_cleanup (IntegratorT& int_state, BurnT& state,
         state.success = false;
     }
 
+    const amrex::Real final_state_species_failure_tolerance =
+        final_state_species_failure_tolerance_factor * integrator_rp::species_failure_tolerance;
+
     for (int n = 1; n <= NumSpec; ++n) {
-        if (int_state.y(n) < -integrator_rp::species_failure_tolerance) {
+        if (int_state.y(n) < -final_state_species_failure_tolerance) {
             state.success = false;
         }