add fallback to finite diff Jacobian

integration/Rosenbrock/rosenbrock_integrator.H

@@ -194,13 +194,29 @@ amrex::Real yass_change_norm (const rosenbrock_t<int_neqs>& rstate)
     return max_change / integrator_rp::yass_epsilon;
 }
 
+template <typename MatrixType, int int_neqs>
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+bool matrix_is_finite (const MatrixType& matrix)
+{
+    for (int j = 1; j <= int_neqs; ++j) {
+        for (int i = 1; i <= int_neqs; ++i) {
+            const amrex::Real value = matrix(i, j);
+            if (std::isnan(value) || std::isinf(value)) {
+                return false;
+            }
+        }
+    }
+
+    return true;
+}
+
 template <int int_neqs>
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
 bool valid_integrator_state (const amrex::Array1D<amrex::Real, 1, int_neqs>& y)
 {
     for (int n = 1; n <= int_neqs; ++n) {
         const amrex::Real yn = y(n);
-        if (yn != yn || std::abs(yn) == std::numeric_limits<amrex::Real>::infinity()) {
+        if (std::isnan(yn) || std::isinf(yn)) {
             return false;
         }
     }
@@ -236,7 +252,10 @@ void evaluate_jacobian (BurnT& state, rosenbrock_t<int_neqs>& rstate, const amre
     if (rstate.jacobian_type == 1) {
         jac(time, state, rstate, rstate.jac);
         rstate.n_jac += 1;
-        return;
+
+        if (matrix_is_finite<decltype(rstate.jac), int_neqs>(rstate.jac)) {
+            return;
+        }
     }
 
     constexpr amrex::Real UROUND = std::numeric_limits<amrex::Real>::epsilon();

integration/integrator_setup_strang.H

@@ -146,6 +146,8 @@ void integrator_cleanup (IntegratorT& int_state, BurnT& state,
     state.n_jac = int_state.n_jac;
     state.n_step = int_state.n_step;
 
+    state.error_code = static_cast<short>(istate);
+
     // The integrator may not always fail even though it can lead to
     // unphysical states.  Add some checks that indicate a burn fail
     // even if the integrator thinks the integration was successful.