add retry framework to burn

integration/BackwardEuler/actual_integrator.H

@@ -6,22 +6,34 @@
 
 template <typename BurnT>
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-void actual_integrator (BurnT& state, const Real dt)
+void actual_integrator (BurnT& state, const Real dt, bool is_retry=false)
 {
     constexpr int int_neqs = integrator_neqs<BurnT>();
 
     be_t<int_neqs> be;
 
     // Set the tolerances.
 
-    be.atol_spec = atol_spec;  // mass fractions
-    be.atol_enuc = atol_enuc;  // energy generated
+    if (!is_retry) {
+        be.atol_spec = atol_spec;  // mass fractions
+        be.atol_enuc = atol_enuc;  // energy generated
 
-    be.rtol_spec = rtol_spec;  // mass fractions
-    be.rtol_enuc = rtol_enuc;  // energy generated
+        be.rtol_spec = rtol_spec;  // mass fractions
+        be.rtol_enuc = rtol_enuc;  // energy generated
+    } else {
+        be.atol_spec = retry_atol_spec; // mass fractions
+        be.atol_enuc = retry_atol_enuc; // energy generated
+
+        be.rtol_spec = retry_rtol_spec; // mass fractions
+        be.rtol_enuc = retry_rtol_enuc; // energy generated
+    }
 
     // set the Jacobian type
-    be.jacobian_type = jacobian;
+    if (is_retry && retry_swap_jacobian) {
+        be.jacobian_type = (jacobian == 1) ? 2 : 1;
+    } else {
+        be.jacobian_type = jacobian;
+    }
 
     // Start off by assuming a successful burn.
 

integration/BackwardEuler/actual_integrator_simplified_sdc.H

@@ -10,15 +10,12 @@
 
 template <typename BurnT>
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-void actual_integrator (BurnT& state, const Real dt)
+void actual_integrator (BurnT& state, const Real dt, bool is_retry=false)
 {
     constexpr int int_neqs = integrator_neqs<BurnT>();
 
     be_t<int_neqs> be;
 
-    // set the Jacobian type
-    be.jacobian_type = jacobian;
-
     // Start off by assuming a successful burn.
 
     state.success = true;
@@ -28,6 +25,13 @@ void actual_integrator (BurnT& state, const Real dt)
     be.t = 0.0;
     be.tout = dt;
 
+    // set the Jacobian type
+    if (is_retry && retry_swap_jacobian) {
+        be.jacobian_type = (jacobian == 1) ? 2 : 1;
+    } else {
+        be.jacobian_type = jacobian;
+    }
+
     // Fill in the initial integration state.
 
     burn_to_int(state, be);
@@ -60,22 +64,38 @@ void actual_integrator (BurnT& state, const Real dt)
 
     Real sdc_min_density = amrex::min(state.rho, state.rho_orig + state.ydot_a[SRHO] * dt);
 
-    be.atol_enuc = sdc_min_density * atol_enuc * sdc_tol_fac;
-    be.rtol_enuc = rtol_enuc * sdc_tol_fac;
+    if (!is_retry) {
+
+        be.atol_enuc = sdc_min_density * atol_enuc * sdc_tol_fac;
+        be.rtol_enuc = rtol_enuc * sdc_tol_fac;
+
+        // Note: we define the input atol for species to refer only to the
+        // mass fraction part, and we multiply by a representative density
+        // so that atol becomes an absolutely tolerance on (rho X)
+
+        be.atol_spec = sdc_min_density * atol_spec * sdc_tol_fac;
+        be.rtol_spec = rtol_spec * sdc_tol_fac;
+
+    } else {
+
+        be.atol_enuc = sdc_min_density * retry_atol_enuc * sdc_tol_fac;
+        be.rtol_enuc = retry_rtol_enuc * sdc_tol_fac;
+
+        // Note: we define the input atol for species to refer only to the
+        // mass fraction part, and we multiply by a representative density
+        // so that atol becomes an absolutely tolerance on (rho X)
+
+        be.atol_spec = sdc_min_density * retry_atol_spec * sdc_tol_fac;
+        be.rtol_spec = retry_rtol_spec * sdc_tol_fac;
+
+    }
 
     if (scale_system) {
         // the absolute tol for energy needs to reflect the scaled
         // energy the integrator sees
         be.atol_enuc /= state.e_scale;
     }
 
-    // Note: we define the input atol for species to refer only to the
-    // mass fraction part, and we multiply by a representative density
-    // so that atol becomes an absolutely tolerance on (rho X)
-
-    be.atol_spec = sdc_min_density * atol_spec * sdc_tol_fac;
-    be.rtol_spec = rtol_spec * sdc_tol_fac;
-
     // Call the integration routine.
 
     int istate = be_integrator(state, be);

integration/ForwardEuler/actual_integrator.H

@@ -157,7 +157,7 @@ void initialize_int_state (IntT& int_state)
 
 template <typename BurnT>
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-void actual_integrator (BurnT& state, Real dt)
+void actual_integrator (BurnT& state, Real dt, bool is_retry=false)
 {
     using namespace microphysics::forward_euler;
 

integration/QSS/actual_integrator.H

@@ -197,7 +197,7 @@ bool corrector (const BurnT& state_0, const Array1D<Real, 1, NumSpec>& f_minus_0
 
 template <typename BurnT>
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-void actual_integrator (BurnT& state, Real dt)
+void actual_integrator (BurnT& state, Real dt, const bool is_retry=false)
 {
     initialize_state(state);
 

integration/RKC/actual_integrator.H

@@ -13,14 +13,28 @@
 
 template <typename BurnT>
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-void actual_integrator (BurnT& state, Real dt)
+void actual_integrator (BurnT& state, Real dt, bool is_retry=false)
 {
     constexpr int int_neqs = integrator_neqs<BurnT>();
 
     rkc_t<int_neqs> rkc_state{};
 
     // Set the tolerances.
 
+    if (!is_retry) {
+        rkc_state.atol_spec = atol_spec;  // mass fractions
+        rkc_state.atol_enuc = atol_enuc;  // energy generated
+
+        rkc_state.rtol_spec = rtol_spec;  // mass fractions
+        rkc_state.rtol_enuc = rtol_enuc;  // energy generated
+    } else {
+        rkc_state.atol_spec = retry_atol_spec; // mass fractions
+        rkc_state.atol_enuc = retry_atol_enuc; // energy generated
+
+        rkc_state.rtol_spec = retry_rtol_spec; // mass fractions
+        rkc_state.rtol_enuc = retry_rtol_enuc; // energy generated
+    }
+
     rkc_state.atol_spec = atol_spec; // mass fractions
     rkc_state.atol_enuc = atol_enuc; // energy generated
 

integration/VODE/actual_integrator.H

@@ -18,22 +18,34 @@ using namespace integrator_rp;
 
 template <typename BurnT>
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-void actual_integrator (BurnT& state, Real dt)
+void actual_integrator (BurnT& state, Real dt, bool is_retry=false)
 {
     constexpr int int_neqs = integrator_neqs<BurnT>();
 
     dvode_t<int_neqs> vode_state{};
 
     // Set the tolerances.
 
-    vode_state.atol_spec = atol_spec; // mass fractions
-    vode_state.atol_enuc = atol_enuc; // energy generated
+    if (!is_retry) {
+        vode_state.atol_spec = atol_spec; // mass fractions
+        vode_state.atol_enuc = atol_enuc; // energy generated
 
-    vode_state.rtol_spec = rtol_spec; // mass fractions
-    vode_state.rtol_enuc = rtol_enuc; // energy generated
+        vode_state.rtol_spec = rtol_spec; // mass fractions
+        vode_state.rtol_enuc = rtol_enuc; // energy generated
+    } else {
+        vode_state.atol_spec = retry_atol_spec; // mass fractions
+        vode_state.atol_enuc = retry_atol_enuc; // energy generated
+
+        vode_state.rtol_spec = retry_rtol_spec; // mass fractions
+        vode_state.rtol_enuc = retry_rtol_enuc; // energy generated
+    }
 
     // set the Jacobian type
-    vode_state.jacobian_type = static_cast<short>(jacobian);
+    if (is_retry && retry_swap_jacobian) {
+        vode_state.jacobian_type = (jacobian == 1) ? 2 : 1;
+    } else {
+        vode_state.jacobian_type = static_cast<short>(jacobian);
+    }
 
     // Start off by assuming a successful burn.
 

integration/VODE/actual_integrator_simplified_sdc.H

@@ -18,15 +18,12 @@ using namespace integrator_rp;
 
 template <typename BurnT>
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-void actual_integrator (BurnT& state, Real dt)
+void actual_integrator (BurnT& state, Real dt, bool is_retry=false)
 {
     constexpr int int_neqs = integrator_neqs<BurnT>();
 
     dvode_t<int_neqs> vode_state{};
 
-    // set the Jacobian type
-    vode_state.jacobian_type = jacobian;
-
     // Start off by assuming a successful burn.
 
     state.success = true;
@@ -40,6 +37,13 @@ void actual_integrator (BurnT& state, Real dt)
 
     vode_state.HMXI = 1.0_rt / ode_max_dt;
 
+    // set the Jacobian type
+    if (is_retry && retry_swap_jacobian) {
+        vode_state.jacobian_type = (jacobian == 1) ? 2 : 1;
+    } else {
+        vode_state.jacobian_type = jacobian;
+    }
+
     // Fill in the initial integration state.
 
     burn_to_int(state, vode_state);
@@ -72,22 +76,38 @@ void actual_integrator (BurnT& state, Real dt)
 
     Real sdc_min_density = amrex::min(state.rho, state.rho_orig + state.ydot_a[SRHO] * dt);
 
-    vode_state.atol_enuc = sdc_min_density * atol_enuc * sdc_tol_fac;
-    vode_state.rtol_enuc = rtol_enuc * sdc_tol_fac;
+    if (!is_retry) {
+
+        vode_state.atol_enuc = sdc_min_density * atol_enuc * sdc_tol_fac;
+        vode_state.rtol_enuc = rtol_enuc * sdc_tol_fac;
+
+        // Note: we define the input atol for species to refer only to the
+        // mass fraction part, and we multiply by a representative density
+        // so that atol becomes an absolutely tolerance on (rho X)
+
+        vode_state.atol_spec = sdc_min_density * atol_spec * sdc_tol_fac;
+        vode_state.rtol_spec = rtol_spec * sdc_tol_fac;
+
+    } else {
+
+        vode_state.atol_enuc = sdc_min_density * retry_atol_enuc * sdc_tol_fac;
+        vode_state.rtol_enuc = retry_rtol_enuc * sdc_tol_fac;
+
+        // Note: we define the input atol for species to refer only to the
+        // mass fraction part, and we multiply by a representative density
+        // so that atol becomes an absolutely tolerance on (rho X)
+
+        vode_state.atol_spec = sdc_min_density * retry_atol_spec * sdc_tol_fac;
+        vode_state.rtol_spec = retry_rtol_spec * sdc_tol_fac;
+
+    }
 
     if (scale_system) {
         // the absolute tol for energy needs to reflect the scaled
         // energy the integrator sees
         vode_state.atol_enuc /= state.e_scale;
     }
 
-    // Note: we define the input atol for species to refer only to the
-    // mass fraction part, and we multiply by a representative density
-    // so that atol becomes an absolutely tolerance on (rho X)
-
-    vode_state.atol_spec = sdc_min_density * atol_spec * sdc_tol_fac;
-    vode_state.rtol_spec = rtol_spec * sdc_tol_fac;
-
     // Call the integration routine.
 
     int istate = dvode(state, vode_state);

integration/_parameters

@@ -60,6 +60,21 @@ nonaka_k                integer           0
 nonaka_level            integer           0
 nonaka_file             character         "nonaka_plot.dat"
 
+# do we retry a failed burn with different parameters?
+use_burn_retry            integer   0
+
+# do we swap the Jacobian (from analytic to numerical or vice versa) on
+# a retry?
+retry_swap_jacobian       integer     1
+
+# Tolerances for the solver (relative and absolute), for the
+# species and energy equations.
+retry_rtol_spec                real      1.d-12
+retry_rtol_enuc                real      1.d-6
+
+retry_atol_spec                real      1.d-8
+retry_atol_enuc                real      1.d-6
+
 # in the clean_state process, do we clip the species such that they
 # are in [0, 1]?
 do_species_clip              integer          1

integration/integrator.H

@@ -7,11 +7,41 @@
 #include <actual_integrator.H>
 #endif
 
+template <typename BurnT, bool enable_retry>
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+void integrator_wrapper (BurnT& state, Real dt)
+{
+
+    if constexpr (enable_retry) {
+        burn_t old_state{state};
+
+        actual_integrator(state, dt);
+
+        if (!state.success) {
+            state = old_state;
+            const bool is_retry = true;
+            actual_integrator(state, dt, is_retry);
+        }
+    } else {
+        actual_integrator(state, dt);
+    }
+
+}
+
+
 template <typename BurnT>
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
 void integrator (BurnT& state, Real dt)
 {
-    actual_integrator(state, dt);
+
+    if (integrator_rp::use_burn_retry) {
+        constexpr bool enable_retry{true};
+        integrator_wrapper<BurnT, enable_retry>(state, dt);
+    } else {
+        constexpr bool enable_retry{false};
+        integrator_wrapper<BurnT, enable_retry>(state, dt);
+
+    }
 }
 
 #endif