From 6e46fc7a6a2a16bdf75ed4228e68fd27c3ff7a21 Mon Sep 17 00:00:00 2001 From: Ben Wibking Date: Mon, 8 Jun 2026 09:32:17 -0400 Subject: [PATCH 1/4] remove Rosenbrock-Euler tableau --- integration/Rosenbrock/_parameters | 3 +- .../Rosenbrock/rosenbrock_integrator.H | 3 -- integration/Rosenbrock/rosenbrock_tableau.H | 37 ------------------- 3 files changed, 1 insertion(+), 42 deletions(-) diff --git a/integration/Rosenbrock/_parameters b/integration/Rosenbrock/_parameters index 680d7962f..361690286 100644 --- a/integration/Rosenbrock/_parameters +++ b/integration/Rosenbrock/_parameters @@ -8,9 +8,8 @@ X_reject_buffer real 1.0 # 0 = Rodas5P 8-stage method from DifferentialEquations.jl # 1 = Rodas4P 6-stage method from DifferentialEquations.jl # 2 = Rodas3P 5-stage method from DifferentialEquations.jl -# 3 = existing 3-stage ROS2S tableau +# 3 = Hamkar et al. (2012) 3-stage ROS2S tableau # 4 = Verwer et al. (1999) 2-stage ROS2 tableau -# 5 = Rosenbrock-Euler 1-stage method rosenbrock_tableau int 0 # H211b timestep controller parameters diff --git a/integration/Rosenbrock/rosenbrock_integrator.H b/integration/Rosenbrock/rosenbrock_integrator.H index 6c8311808..fc324a9b3 100644 --- a/integration/Rosenbrock/rosenbrock_integrator.H +++ b/integration/Rosenbrock/rosenbrock_integrator.H @@ -600,9 +600,6 @@ int rosenbrock_integrator (BurnT& state, rosenbrock_t()>& if (integrator_rp::rosenbrock_tableau == 4) { return rosenbrock_integrator(state, rstate); } - if (integrator_rp::rosenbrock_tableau == 5) { - return rosenbrock_integrator(state, rstate); - } return rosenbrock_integrator(state, rstate); } diff --git a/integration/Rosenbrock/rosenbrock_tableau.H b/integration/Rosenbrock/rosenbrock_tableau.H index c4a2fd839..3eb47b620 100644 --- a/integration/Rosenbrock/rosenbrock_tableau.H +++ b/integration/Rosenbrock/rosenbrock_tableau.H @@ -530,43 +530,6 @@ struct rodas4p_tableau { } }; -struct rosenbrock_euler_tableau { - static constexpr int stages = 1; - static constexpr amrex::Real gamma = 1.0_rt; - - AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE - static constexpr amrex::Real ct (const int i) { - (void) i; - return 0.0_rt; - } - - AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE - static constexpr amrex::Real a (const int i, const int j) { - (void) i; - (void) j; - return 0.0_rt; - } - - AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE - static constexpr amrex::Real c (const int i, const int j) { - (void) i; - (void) j; - return 0.0_rt; - } - - AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE - static constexpr amrex::Real b (const int i) { - (void) i; - return 1.0_rt; - } - - AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE - static constexpr amrex::Real e (const int i) { - (void) i; - return 1.0_rt; - } -}; - using coefficients = rodas5p_tableau; } From 7afaee90f89718a64c740a0fc7a44bc34d69744a Mon Sep 17 00:00:00 2001 From: Ben Wibking Date: Mon, 8 Jun 2026 09:35:12 -0400 Subject: [PATCH 2/4] trigger CI --- integration/Rosenbrock/_parameters | 1 + 1 file changed, 1 insertion(+) diff --git a/integration/Rosenbrock/_parameters b/integration/Rosenbrock/_parameters index 361690286..9d772d860 100644 --- a/integration/Rosenbrock/_parameters +++ b/integration/Rosenbrock/_parameters @@ -18,3 +18,4 @@ h211b_k real 2.5 h211b_fac_min real 0.2 h211b_fac_max real 6.0 h211b_reduction_fac real 0.5 + From 9d654f3e01ae8e8c264095a4ea17ecb8685e02a2 Mon Sep 17 00:00:00 2001 From: Ben Wibking Date: Mon, 8 Jun 2026 12:10:32 -0400 Subject: [PATCH 3/4] reduce memory usage --- integration/Rosenbrock/actual_integrator.H | 10 + .../Rosenbrock/actual_integrator_sdc.H | 10 + .../Rosenbrock/rosenbrock_integrator.H | 850 +++++++++++++++++- integration/Rosenbrock/rosenbrock_type.H | 31 + 4 files changed, 897 insertions(+), 4 deletions(-) diff --git a/integration/Rosenbrock/actual_integrator.H b/integration/Rosenbrock/actual_integrator.H index 07252876d..da5224c9f 100644 --- a/integration/Rosenbrock/actual_integrator.H +++ b/integration/Rosenbrock/actual_integrator.H @@ -17,6 +17,16 @@ void actual_integrator (BurnT& state, amrex::Real dt, bool is_retry=false) { constexpr int int_neqs = integrator_neqs(); + if (integrator_rp::rosenbrock_tableau == 0) { + auto rodas5p_state = integrator_setup>(state, dt, is_retry); + auto state_save = integrator_backup(state); + + int istate = rodas5p_integrator(state, rodas5p_state); + + integrator_cleanup(rodas5p_state, state, istate, state_save, dt); + return; + } + auto rosenbrock_state = integrator_setup>(state, dt, is_retry); auto state_save = integrator_backup(state); diff --git a/integration/Rosenbrock/actual_integrator_sdc.H b/integration/Rosenbrock/actual_integrator_sdc.H index 32bd901c8..0ec9fcc02 100644 --- a/integration/Rosenbrock/actual_integrator_sdc.H +++ b/integration/Rosenbrock/actual_integrator_sdc.H @@ -17,6 +17,16 @@ void actual_integrator (BurnT& state, amrex::Real dt, bool is_retry=false) { constexpr int int_neqs = integrator_neqs(); + if (integrator_rp::rosenbrock_tableau == 0) { + auto rodas5p_state = integrator_setup>(state, dt, is_retry); + auto state_save = integrator_backup(state); + + int istate = rodas5p_integrator(state, rodas5p_state); + + integrator_cleanup(rodas5p_state, state, istate, state_save, dt); + return; + } + auto rosenbrock_state = integrator_setup>(state, dt, is_retry); auto state_save = integrator_backup(state); diff --git a/integration/Rosenbrock/rosenbrock_integrator.H b/integration/Rosenbrock/rosenbrock_integrator.H index fc324a9b3..b48c7f172 100644 --- a/integration/Rosenbrock/rosenbrock_integrator.H +++ b/integration/Rosenbrock/rosenbrock_integrator.H @@ -217,6 +217,374 @@ bool valid_integrator_state (const amrex::Array1D& y) return true; } +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +bool valid_integrator_component (const int n, const amrex::Real yn) +{ + if (std::isnan(yn) || std::isinf(yn)) { + return false; + } + +#ifdef STRANG + if (n <= NumSpec) { + if (yn < -integrator_rp::atol_spec) { + return false; + } + + if (! integrator_rp::use_number_densities && yn > 1.0_rt) { + return false; + } + } + + if (n == net_ienuc && yn <= 0.0_rt) { + return false; + } +#endif + + return true; +} + +template +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +bool valid_rodas5p_stage5 (const rodas5p_t& rstate) +{ + using C = rodas5p_tableau; + + for (int n = 1; n <= int_neqs; ++n) { + const amrex::Real yn = rstate.y(n) + + C::a(5, 1) * rstate.ak1(n) + + C::a(5, 2) * rstate.ak2(n) + + C::a(5, 3) * rstate.ak3(n) + + C::a(5, 4) * rstate.ak4(n); + if (! valid_integrator_component(n, yn)) { + return false; + } + } + + return true; +} + +template +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +bool valid_rodas5p_update_stage (const rodas5p_t& rstate, + const amrex::Array1D& update) +{ + for (int n = 1; n <= int_neqs; ++n) { + if (! valid_integrator_component(n, rstate.y(n) + update(n))) { + return false; + } + } + + return true; +} + +#ifdef STRANG +template +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +void finish_rodas5p_stage_state (const amrex::Real time, BurnT& state) +{ + if (integrator_rp::do_species_clip) { + for (int n = 1; n <= NumSpec; ++n) { + state.xn[n-1] = amrex::Clamp(state.xn[n-1], integrator_rp::SMALL_X_SAFE, 1.0_rt); + } + } + + if (integrator_rp::use_number_densities) { + for (int n = 1; n <= NumSpec; ++n) { + state.xn[n-1] = amrex::max(state.xn[n-1], integrator_rp::SMALL_X_SAFE); + } + } + + if (integrator_rp::renormalize_abundances) { + amrex::Real sum = 0.0_rt; + + for (int n = 1; n <= NumSpec; ++n) { + sum += state.xn[n-1]; + } + + for (int n = 1; n <= NumSpec; ++n) { + state.xn[n-1] /= sum; + } + } + +#ifdef AUX_THERMO + state.abar = 0.0_rt; + state.y_e = 0.0_rt; + for (int n = 1; n <= NumSpec; ++n) { + state.abar += state.xn[n-1] * aion_inv[n-1]; + state.y_e += state.xn[n-1] * zion[n-1] * aion_inv[n-1]; + } + state.abar = 1.0_rt / state.abar; + state.zbar = state.abar * state.y_e; + + state.aux[AuxZero::iabar] = state.abar; + state.aux[AuxZero::iye] = state.y_e; + state.aux[AuxZero::ibea] = 0.0_rt; +#endif + + if (integrator_rp::scale_system) { + state.e *= state.e_scale; + } + + if (integrator_rp::call_eos_in_rhs) { + eos(eos_input_re, state); + } + + if (state.T_fixed > 0.0_rt) { + state.T = state.T_fixed; + } + + state.time = time; +} +#endif + +#ifdef SDC +template +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +void finish_rodas5p_stage_state (const amrex::Real time, BurnT& state) +{ + update_density_in_time(time, state); + + if (integrator_rp::do_species_clip) { + for (int n = 1; n <= NumSpec; ++n) { + state.y[SFS+n-1] = amrex::Clamp(state.y[SFS+n-1], + state.rho * integrator_rp::SMALL_X_SAFE, + state.rho); + } + } + + if (integrator_rp::renormalize_abundances) { + amrex::Real nspec_sum = 0.0_rt; + for (int n = 1; n <= NumSpec; n++) { + nspec_sum += state.y[SFS+n-1]; + } + nspec_sum /= state.y[SRHO]; + + for (int n = 1; n <= NumSpec; n++) { + state.y[SFS+n-1] /= nspec_sum; + } + } + + if (integrator_rp::scale_system) { + state.y[SEINT] *= state.e_scale; + } + + update_density_in_time(time, state); + + amrex::Real rhoInv = 1.0_rt / state.rho; + + for (int n = 0; n < NumSpec; n++) { + state.xn[n] = state.y[SFS+n] * rhoInv; + } + +#ifdef AUX_THERMO + set_aux_comp_from_X(state); + + for (int n = 0; n < NumAux; n++) { + state.y[SFX+n] = state.rho * state.aux[n]; + } +#endif + + state.e = state.y[SEINT] * rhoInv; + + if (integrator_rp::call_eos_in_rhs) { + eos(eos_input_re, state); + } + + if (state.T_fixed > 0.0_rt) { + state.T = state.T_fixed; + } + + state.time = time; +} +#endif + +template +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +void prepare_rodas5p_stage5_state (const amrex::Real time, BurnT& state, + const rodas5p_t& rstate) +{ + using C = rodas5p_tableau; + +#ifdef STRANG + for (int n = 1; n <= NumSpec; ++n) { + const amrex::Real yn = rstate.y(n) + + C::a(5, 1) * rstate.ak1(n) + + C::a(5, 2) * rstate.ak2(n) + + C::a(5, 3) * rstate.ak3(n) + + C::a(5, 4) * rstate.ak4(n); + state.y[n-1] = yn; + state.xn[n-1] = yn; + } + + state.y[SEINT] = rstate.y(net_ienuc) + + C::a(5, 1) * rstate.ak1(net_ienuc) + + C::a(5, 2) * rstate.ak2(net_ienuc) + + C::a(5, 3) * rstate.ak3(net_ienuc) + + C::a(5, 4) * rstate.ak4(net_ienuc); + state.e = state.y[SEINT]; +#endif + +#ifdef SDC + for (int n = 0; n < SVAR_EVOLVE; n++) { + const int i = n + 1; + state.y[n] = rstate.y(i) + + C::a(5, 1) * rstate.ak1(i) + + C::a(5, 2) * rstate.ak2(i) + + C::a(5, 3) * rstate.ak3(i) + + C::a(5, 4) * rstate.ak4(i); + } +#endif + + finish_rodas5p_stage_state(time, state); +} + +template +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +void prepare_rodas5p_update_stage_state (const amrex::Real time, BurnT& state, + const rodas5p_t& rstate, + const amrex::Array1D& update) +{ +#ifdef STRANG + for (int n = 1; n <= NumSpec; ++n) { + const amrex::Real yn = rstate.y(n) + update(n); + state.y[n-1] = yn; + state.xn[n-1] = yn; + } + + state.y[SEINT] = rstate.y(net_ienuc) + update(net_ienuc); + state.e = state.y[SEINT]; +#endif + +#ifdef SDC + for (int n = 0; n < SVAR_EVOLVE; n++) { + const int i = n + 1; + state.y[n] = rstate.y(i) + update(i); + } +#endif + + finish_rodas5p_stage_state(time, state); +} + +template +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +amrex::Real rodas5p_stage_component (const BurnT& state, const int n) +{ +#ifdef STRANG + if (n <= NumSpec) { + return state.y[n-1]; + } + + return state.y[SEINT]; +#endif + +#ifdef SDC + const int i = n - 1; + if (i == SEINT && integrator_rp::scale_system) { + return state.y[i] / state.e_scale; + } + + return state.y[i]; +#endif +} + +template +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +void rodas5p_rhs_from_prepared_state (const amrex::Real time, BurnT& state, + RArray1D& ydot, + [[maybe_unused]] const bool in_jacobian=false) +{ + amrex::ignore_unused(time); + + if (state.T <= EOSData::mintemp || state.T >= integrator_rp::MAX_TEMP) { + for (int n = 1; n <= INT_NEQS; ++n) { + ydot(n) = 0.0_rt; + } + + return; + } + +#ifdef STRANG +#ifdef NEW_NETWORK_IMPLEMENTATION + RHS::rhs(state, ydot); +#else + actual_rhs(state, ydot); +#endif + +#ifdef NONAKA_PLOT + if (! in_jacobian) { + nonaka_rhs(time, state, ydot); + } +#endif + + if (! integrator_rp::use_number_densities) { + for (int n = 1; n <= NumSpec; ++n) { + ydot(n) *= aion[n-1]; + } + } + + if (integrator_rp::scale_system) { + ydot(net_ienuc) /= state.e_scale; + } + + if (! integrator_rp::integrate_energy) { + ydot(net_ienuc) = 0.0_rt; + } + + if (integrator_rp::react_boost > 0.0_rt) { + for (int n = 1; n <= INT_NEQS; ++n) { + ydot(n) *= integrator_rp::react_boost; + } + } +#endif + +#ifdef SDC + actual_rhs(state, ydot); + +#ifdef NONAKA_PLOT + if (! in_jacobian) { + nonaka_rhs(time, state, ydot); + } +#endif + + if (integrator_rp::react_boost > 0.0_rt) { + for (int n = 1; n <= neqs; ++n) { + ydot(n) *= integrator_rp::react_boost; + } + } + + rhs_to_int(time, state, ydot); +#endif +} + +template +struct rhs_array_ref { + amrex::Array1D* data; + + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + amrex::Real& operator() (const int n) const { + return (*data)(n); + } +}; + +template +struct rhs_state_view { + amrex::Real t; + rhs_array_ref y; +}; + +template +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +void rhs_at_array (const amrex::Real time, BurnT& state, RosenbrockT& rstate, + amrex::Array1D& y, + RArray1D& ydot) +{ + amrex::ignore_unused(rstate); + constexpr int int_neqs = RosenbrockT::neqs; + rhs_state_view rhs_state{time, rhs_array_ref{&y}}; + + rhs(time, state, rhs_state, ydot); +} + template AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void evaluate_jacobian (BurnT& state, rosenbrock_t& rstate, const amrex::Real time) @@ -281,6 +649,72 @@ void evaluate_jacobian (BurnT& state, rosenbrock_t& rstate, const amre rstate.n_jac += 1; } +template +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +void evaluate_jacobian_rodas5p (BurnT& state, rodas5p_t& rstate, + const amrex::Real time) +{ + constexpr bool in_jacobian = true; + + rhs(time, state, rstate, rstate.ak1, in_jacobian); + rstate.n_rhs += 1; + + if (rstate.jacobian_type == 1) { + jac(time, state, rstate, rstate.jac); + rstate.n_jac += 1; + + if (matrix_is_finite(rstate.jac)) { + return; + } + } + + constexpr amrex::Real UROUND = std::numeric_limits::epsilon(); + + // rhs() may clean or renormalize all components of rstate.y. + for (int i = 1; i <= int_neqs; ++i) { + rstate.ak2(i) = rstate.y(i); + } + + amrex::Real fac = 0.0_rt; + for (int i = 1; i <= int_neqs; ++i) { + const amrex::Real ewt = + 1.0_rt / (rtol_for(rstate, i) * std::abs(rstate.y(i)) + atol_for(rstate, i)); + fac += (rstate.ak1(i) * ewt) * (rstate.ak1(i) * ewt); + } + fac = std::sqrt(fac / static_cast(int_neqs)); + + amrex::Real R0 = 1000.0_rt * std::abs(rstate.dt) * + UROUND * static_cast(int_neqs) * fac; + if (R0 == 0.0_rt) { + R0 = 1.0_rt; + } + + for (int j = 1; j <= int_neqs; ++j) { + const amrex::Real yj = rstate.ak2(j); + const amrex::Real ewt_j = + 1.0_rt / (rtol_for(rstate, j) * std::abs(yj) + atol_for(rstate, j)); + const amrex::Real R = amrex::max(std::sqrt(UROUND) * std::abs(yj), R0 / ewt_j); + for (int n = 1; n <= int_neqs; ++n) { + rstate.y(n) = rstate.ak2(n); + } + rstate.y(j) += R; + + rhs(time, state, rstate, rstate.ak3, in_jacobian); + + fac = 1.0_rt / R; + for (int i = 1; i <= int_neqs; ++i) { + rstate.jac.set(i, j, (rstate.ak3(i) - rstate.ak1(i)) * fac); + } + } + + for (int n = 1; n <= int_neqs; ++n) { + rstate.y(n) = rstate.ak2(n); + } + + rstate.n_rhs += int_neqs; + rstate.n_jac += 1; +} + template AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void rhs_at (const amrex::Real time, BurnT& state, rosenbrock_t& rstate, @@ -301,10 +735,12 @@ void rhs_at (const amrex::Real time, BurnT& state, rosenbrock_t& rstat } } -template +template AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -int decompose (rosenbrock_t& rstate, const amrex::Real fac) +int decompose (RosenbrockT& rstate, const amrex::Real fac) { + constexpr int int_neqs = RosenbrockT::neqs; + for (int j = 1; j <= int_neqs; ++j) { for (int i = 1; i <= int_neqs; ++i) { rstate.jac(i, j) = -rstate.jac(i, j); @@ -324,10 +760,12 @@ int decompose (rosenbrock_t& rstate, const amrex::Real fac) return ierr_linpack; } -template +template AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void solve (rosenbrock_t& rstate, RArray1D& b) +void solve (RosenbrockT& rstate, RArray1D& b) { + constexpr int int_neqs = RosenbrockT::neqs; + if (integrator_rp::linalg_do_pivoting == 1) { constexpr bool allow_pivot{true}; dgesl(rstate.jac, rstate.pivot, b); @@ -584,6 +1022,410 @@ int rosenbrock_integrator (BurnT& state, rosenbrock_t()>& } } +namespace rosenbrock { + +template +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +amrex::Real error_norm (const rodas5p_t& rstate, + const amrex::Array1D& ynew, + const amrex::Array1D& error) +{ + amrex::Real err = 0.0_rt; + for (int n = 1; n <= int_neqs; ++n) { + const amrex::Real sk = atol_for(rstate, n) + + rtol_for(rstate, n) * amrex::max(std::abs(rstate.y(n)), std::abs(ynew(n))); + const amrex::Real term = error(n) / sk; + err += term * term; + } + return std::sqrt(err / static_cast(int_neqs)); +} + +template +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +bool species_change_valid (const rodas5p_t& rstate, + const amrex::Array1D& ynew) +{ +#ifdef STRANG + constexpr amrex::Real increase_change_factor = 4.0_rt; + constexpr amrex::Real decrease_change_factor = 0.25_rt; + + for (int i = 1; i <= NumSpec; ++i) { + if (std::abs(rstate.y(i)) > integrator_rp::X_reject_buffer * rstate.atol_spec && + std::abs(ynew(i)) > integrator_rp::X_reject_buffer * rstate.atol_spec && + (std::abs(ynew(i)) > increase_change_factor * std::abs(rstate.y(i)) || + std::abs(ynew(i)) < decrease_change_factor * std::abs(rstate.y(i)))) { + return false; + } + } +#endif + + return true; +} + +} // namespace rosenbrock + +template +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +int rodas5p_integrator (BurnT& state, rodas5p_t()>& rstate) +{ + using C = rosenbrock::rodas5p_tableau; + constexpr int int_neqs = integrator_neqs(); + + if (rstate.tout == rstate.t) { + return IERR_SUCCESS; + } + + for (int n = 1; n <= int_neqs; ++n) { + if (rosenbrock::atol_for(rstate, n) <= 0.0_rt || + rosenbrock::rtol_for(rstate, n) <= 10.0_rt * std::numeric_limits::epsilon()) { + return IERR_TOO_MUCH_ACCURACY_REQUESTED; + } + } + + rstate.n_rhs = 0; + rstate.n_jac = 0; + rstate.n_step = 0; + + rstate.dt = rstate.tout - rstate.t; + + const amrex::Real posneg = rstate.tout >= rstate.t ? 1.0_rt : -1.0_rt; + const amrex::Real hmax = amrex::min(integrator_rp::ode_max_dt, + std::abs(rstate.tout - rstate.t)); + amrex::Real h = amrex::min(std::abs(rstate.dt), hmax) * posneg; + if (std::abs(h) <= 10.0_rt * std::numeric_limits::epsilon()) { + h = 1.e-6_rt * posneg; + } + + bool reject = false; + bool last = false; + int nsing = 0; + int n_reject = 0; + amrex::Real facold = 1.0_rt; + amrex::Real errold = 1.0_rt; + amrex::Real hopt = h; + amrex::Real x = rstate.t; + + while (true) { + + if (rstate.n_step > integrator_rp::ode_max_steps) { + rstate.t = x; + rstate.dt = h; + return IERR_TOO_MANY_STEPS; + } + + if (0.1_rt * std::abs(h) <= std::abs(x) * std::numeric_limits::epsilon()) { + rstate.t = x; + rstate.dt = h; + return IERR_DT_UNDERFLOW; + } + + if (last) { + rstate.t = x; + rstate.dt = hopt; + return IERR_SUCCESS; + } + + hopt = h; + if ((x + h * (1.0_rt + timestep_safety_factor) - rstate.tout) * posneg >= 0.0_rt) { + h = rstate.tout - x; + last = true; + } + +#ifdef NSE + if (rstate.n_step > MIN_NSE_BAILOUT_STEPS && x <= rstate.tout) { + int_to_burn(x, rstate, state); + + if (in_nse(state)) { + rstate.t = x; + return IERR_ENTERED_NSE; + } + } +#endif + + while (true) { + rstate.dt = h; + + if (0.1_rt * std::abs(h) <= std::abs(x) * std::numeric_limits::epsilon()) { + rstate.t = x; + rstate.dt = h; + return IERR_DT_UNDERFLOW; + } + + const amrex::Real fac = 1.0_rt / (h * C::gamma); + rosenbrock::evaluate_jacobian_rodas5p(state, rstate, x); + int ierr_linpack = rosenbrock::decompose(rstate, fac); + + if (ierr_linpack != 0) { + nsing += 1; + if (nsing >= 5) { + rstate.t = x; + rstate.dt = h; + return IERR_LU_DECOMPOSITION_ERROR; + } + h *= 0.5_rt; + reject = true; + last = false; + continue; + } + + rosenbrock::solve(rstate, rstate.ak1); + + bool valid_stage_state = true; + + for (int n = 1; n <= int_neqs; ++n) { + rstate.ak2(n) = rstate.y(n) + C::a(2, 1) * rstate.ak1(n); + } + valid_stage_state = rosenbrock::valid_integrator_state(rstate.ak2); + if (valid_stage_state) { + rosenbrock::rhs_at_array(x + C::ct(2) * h, state, rstate, rstate.ak2, rstate.ak2); + rstate.n_rhs += 1; + for (int n = 1; n <= int_neqs; ++n) { + rstate.ak2(n) += (C::c(2, 1) / h) * rstate.ak1(n); + } + rosenbrock::solve(rstate, rstate.ak2); + } + + if (valid_stage_state) { + for (int n = 1; n <= int_neqs; ++n) { + rstate.ak3(n) = rstate.y(n) + + C::a(3, 1) * rstate.ak1(n) + + C::a(3, 2) * rstate.ak2(n); + } + valid_stage_state = rosenbrock::valid_integrator_state(rstate.ak3); + } + if (valid_stage_state) { + rosenbrock::rhs_at_array(x + C::ct(3) * h, state, rstate, rstate.ak3, rstate.ak3); + rstate.n_rhs += 1; + for (int n = 1; n <= int_neqs; ++n) { + rstate.ak3(n) += + (C::c(3, 1) * rstate.ak1(n) + + C::c(3, 2) * rstate.ak2(n)) / h; + } + rosenbrock::solve(rstate, rstate.ak3); + } + + if (valid_stage_state) { + for (int n = 1; n <= int_neqs; ++n) { + rstate.ak4(n) = rstate.y(n) + + C::a(4, 1) * rstate.ak1(n) + + C::a(4, 2) * rstate.ak2(n) + + C::a(4, 3) * rstate.ak3(n); + } + valid_stage_state = rosenbrock::valid_integrator_state(rstate.ak4); + } + if (valid_stage_state) { + rosenbrock::rhs_at_array(x + C::ct(4) * h, state, rstate, rstate.ak4, rstate.ak4); + rstate.n_rhs += 1; + for (int n = 1; n <= int_neqs; ++n) { + rstate.ak4(n) += + (C::c(4, 1) * rstate.ak1(n) + + C::c(4, 2) * rstate.ak2(n) + + C::c(4, 3) * rstate.ak3(n)) / h; + } + rosenbrock::solve(rstate, rstate.ak4); + } + + if (valid_stage_state) { + valid_stage_state = rosenbrock::valid_rodas5p_stage5(rstate); + } + if (valid_stage_state) { + rosenbrock::prepare_rodas5p_stage5_state(x + C::ct(5) * h, state, rstate); + rosenbrock::rodas5p_rhs_from_prepared_state(x + C::ct(5) * h, state, rstate.ak4); + rstate.n_rhs += 1; + for (int n = 1; n <= int_neqs; ++n) { + const amrex::Real k4 = + (rosenbrock::rodas5p_stage_component(state, n) - + rstate.y(n) - + C::a(5, 1) * rstate.ak1(n) - + C::a(5, 2) * rstate.ak2(n) - + C::a(5, 3) * rstate.ak3(n)) / C::a(5, 4); + + rstate.ak4(n) += + (C::c(5, 1) * rstate.ak1(n) + + C::c(5, 2) * rstate.ak2(n) + + C::c(5, 3) * rstate.ak3(n) + + C::c(5, 4) * k4) / h; + } + rosenbrock::solve(rstate, rstate.ak4); + } + + if (! valid_stage_state) { + h *= 0.25_rt; + reject = true; + n_reject += 1; + last = false; + continue; + } + + for (int n = 1; n <= int_neqs; ++n) { + const amrex::Real k1 = rstate.ak1(n); + const amrex::Real k2 = rstate.ak2(n); + const amrex::Real k3 = rstate.ak3(n); + const amrex::Real k5 = rstate.ak4(n); + const amrex::Real k4 = + (rosenbrock::rodas5p_stage_component(state, n) - + rstate.y(n) - + C::a(5, 1) * k1 - + C::a(5, 2) * k2 - + C::a(5, 3) * k3) / C::a(5, 4); + + rstate.ak1(n) = + C::b(1) * k1 + C::b(2) * k2 + C::b(3) * k3 + + C::b(4) * k4 + C::b(5) * k5; + rstate.ak2(n) = + C::c(6, 1) * k1 + C::c(6, 2) * k2 + C::c(6, 3) * k3 + + C::c(6, 4) * k4 + C::c(6, 5) * k5; + rstate.ak3(n) = + C::c(7, 1) * k1 + C::c(7, 2) * k2 + C::c(7, 3) * k3 + + C::c(7, 4) * k4 + C::c(7, 5) * k5; + rstate.ak4(n) = + C::c(8, 1) * k1 + C::c(8, 2) * k2 + C::c(8, 3) * k3 + + C::c(8, 4) * k4 + C::c(8, 5) * k5; + } + + valid_stage_state = rosenbrock::valid_rodas5p_update_stage(rstate, rstate.ak1); + if (valid_stage_state) { + rosenbrock::prepare_rodas5p_update_stage_state(x + C::ct(6) * h, + state, rstate, rstate.ak1); + rosenbrock::rodas5p_rhs_from_prepared_state(x + C::ct(6) * h, state, rstate.ak1); + rstate.n_rhs += 1; + for (int n = 1; n <= int_neqs; ++n) { + rstate.ak1(n) += rstate.ak2(n) / h; + } + rosenbrock::solve(rstate, rstate.ak1); + } + + if (valid_stage_state) { + for (int n = 1; n <= int_neqs; ++n) { + const amrex::Real update = + rosenbrock::rodas5p_stage_component(state, n) - rstate.y(n); + const amrex::Real k6 = rstate.ak1(n); + rstate.ak1(n) = update + k6; + rstate.ak3(n) += C::c(7, 6) * k6; + rstate.ak4(n) += C::c(8, 6) * k6; + } + } else { + h *= 0.25_rt; + reject = true; + n_reject += 1; + last = false; + continue; + } + + valid_stage_state = rosenbrock::valid_rodas5p_update_stage(rstate, rstate.ak1); + if (valid_stage_state) { + rosenbrock::prepare_rodas5p_update_stage_state(x + C::ct(7) * h, + state, rstate, rstate.ak1); + rosenbrock::rodas5p_rhs_from_prepared_state(x + C::ct(7) * h, state, rstate.ak1); + rstate.n_rhs += 1; + for (int n = 1; n <= int_neqs; ++n) { + rstate.ak1(n) += rstate.ak3(n) / h; + } + rosenbrock::solve(rstate, rstate.ak1); + } + + if (valid_stage_state) { + for (int n = 1; n <= int_neqs; ++n) { + const amrex::Real update = + rosenbrock::rodas5p_stage_component(state, n) - rstate.y(n); + const amrex::Real k7 = rstate.ak1(n); + rstate.ak1(n) = update + k7; + rstate.ak4(n) += C::c(8, 7) * k7; + } + } else { + h *= 0.25_rt; + reject = true; + n_reject += 1; + last = false; + continue; + } + + valid_stage_state = rosenbrock::valid_rodas5p_update_stage(rstate, rstate.ak1); + if (valid_stage_state) { + rosenbrock::prepare_rodas5p_update_stage_state(x + C::ct(8) * h, + state, rstate, rstate.ak1); + rosenbrock::rodas5p_rhs_from_prepared_state(x + C::ct(8) * h, state, rstate.ak1); + rstate.n_rhs += 1; + for (int n = 1; n <= int_neqs; ++n) { + rstate.ak1(n) += rstate.ak4(n) / h; + } + rosenbrock::solve(rstate, rstate.ak1); + } + + if (! valid_stage_state) { + h *= 0.25_rt; + reject = true; + n_reject += 1; + last = false; + continue; + } + + for (int n = 1; n <= int_neqs; ++n) { + const amrex::Real update = + rosenbrock::rodas5p_stage_component(state, n) - rstate.y(n); + rstate.ak2(n) = rstate.y(n) + update + rstate.ak1(n); + } + + rstate.n_step += 1; + + const bool valid_state = rosenbrock::valid_integrator_state(rstate.ak2); + if (! valid_state) { + reject = true; + n_reject += 1; + last = false; + h *= 0.25_rt; + continue; + } + + const bool valid_update = rosenbrock::species_change_valid(rstate, rstate.ak2); + constexpr amrex::Real err_min = 1.e-10_rt; + amrex::Real err = rosenbrock::error_norm(rstate, rstate.ak2, rstate.ak1); + amrex::Real controller_fac = amrex::Clamp( + std::pow(1.0_rt / amrex::max(err, err_min), + 1.0_rt / (integrator_rp::h211b_b * integrator_rp::h211b_k)) * + std::pow(1.0_rt / amrex::max(errold, err_min), + 1.0_rt / (integrator_rp::h211b_b * integrator_rp::h211b_k)) * + std::pow(facold, -1.0_rt / integrator_rp::h211b_b), + integrator_rp::h211b_fac_min, integrator_rp::h211b_fac_max); + facold = controller_fac; + errold = err; + + amrex::Real hnew = h * controller_fac; + + if (err <= 1.0_rt && valid_update) { + for (int n = 1; n <= int_neqs; ++n) { + rstate.y(n) = rstate.ak2(n); + } + x += h; + + if (std::abs(hnew) > hmax) { + hnew = posneg * hmax; + } + if (reject) { + hnew = posneg * amrex::min(std::abs(hnew), std::abs(h)); + } + reject = false; + n_reject = 0; + h = hnew; + break; + } + + reject = true; + n_reject += 1; + last = false; + if (n_reject >= 2) { + hnew *= integrator_rp::h211b_reduction_fac; + } + if (valid_update) { + h = posneg * amrex::min(std::abs(hnew), + integrator_rp::h211b_reduction_fac * std::abs(h)); + } else { + h *= 0.25_rt; + } + } + } +} + template AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE int rosenbrock_integrator (BurnT& state, rosenbrock_t()>& rstate) diff --git a/integration/Rosenbrock/rosenbrock_type.H b/integration/Rosenbrock/rosenbrock_type.H index a7751d9ff..f4dc7729d 100644 --- a/integration/Rosenbrock/rosenbrock_type.H +++ b/integration/Rosenbrock/rosenbrock_type.H @@ -22,6 +22,7 @@ const amrex::Real timestep_safety_factor = 1.0e-4_rt; template struct rosenbrock_t { + static constexpr int neqs = int_neqs; amrex::Real t; // the starting time amrex::Real tout; // the stopping time @@ -55,4 +56,34 @@ struct rosenbrock_t { short jacobian_type; }; +template +struct rodas5p_t { + static constexpr int neqs = int_neqs; + + amrex::Real t; // the starting time + amrex::Real tout; // the stopping time + amrex::Real dt; // next internal timestep + + int n_step; + int n_rhs; + int n_jac; + + amrex::Real atol_spec; + amrex::Real rtol_spec; + + amrex::Real atol_enuc; + amrex::Real rtol_enuc; + + amrex::Array1D y; + amrex::Array1D ak1; + amrex::Array1D ak2; + amrex::Array1D ak3; + amrex::Array1D ak4; + + ArrayUtil::MathArray2D jac; + IArray1D pivot; + + short jacobian_type; +}; + #endif From 1bf7a0e0d66cd131b56e78b8b83f4d1d239647f7 Mon Sep 17 00:00:00 2001 From: Ben Wibking Date: Mon, 8 Jun 2026 12:22:11 -0400 Subject: [PATCH 4/4] trigger CI --- integration/Rosenbrock/_parameters | 1 - 1 file changed, 1 deletion(-) diff --git a/integration/Rosenbrock/_parameters b/integration/Rosenbrock/_parameters index 9d772d860..361690286 100644 --- a/integration/Rosenbrock/_parameters +++ b/integration/Rosenbrock/_parameters @@ -18,4 +18,3 @@ h211b_k real 2.5 h211b_fac_min real 0.2 h211b_fac_max real 6.0 h211b_reduction_fac real 0.5 -