diff --git a/integration/Rosenbrock/rosenbrock_integrator.H b/integration/Rosenbrock/rosenbrock_integrator.H index c9b5f792e..e7104c2c6 100644 --- a/integration/Rosenbrock/rosenbrock_integrator.H +++ b/integration/Rosenbrock/rosenbrock_integrator.H @@ -194,13 +194,29 @@ amrex::Real yass_change_norm (const rosenbrock_t& rstate) return max_change / integrator_rp::yass_epsilon; } +template +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 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool valid_integrator_state (const amrex::Array1D& y) { for (int n = 1; n <= int_neqs; ++n) { const amrex::Real yn = y(n); - if (yn != yn || std::abs(yn) == std::numeric_limits::infinity()) { + if (std::isnan(yn) || std::isinf(yn)) { return false; } } @@ -236,7 +252,10 @@ void evaluate_jacobian (BurnT& state, rosenbrock_t& rstate, const amre if (rstate.jacobian_type == 1) { jac(time, state, rstate, rstate.jac); rstate.n_jac += 1; - return; + + if (matrix_is_finite(rstate.jac)) { + return; + } } constexpr amrex::Real UROUND = std::numeric_limits::epsilon(); diff --git a/integration/integrator_setup_strang.H b/integration/integrator_setup_strang.H index 03252ab15..60f874ab4 100644 --- a/integration/integrator_setup_strang.H +++ b/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(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.