diff --git a/src/algorithms.jl b/src/algorithms.jl index b44516cd18..68c6b77a49 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -3085,20 +3085,23 @@ for Alg in [ :Rodas5Pe, :Rodas5Pr] @eval begin - struct $Alg{CS, AD, F, P, FDT, ST, CJ, StepLimiter} <: + struct $Alg{CS, AD, F, P, FDT, ST, CJ, StepLimiter, StageLimiter} <: OrdinaryDiffEqRosenbrockAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} linsolve::F precs::P step_limiter!::StepLimiter + stage_limiter!::StageLimiter end function $Alg(; chunk_size = Val{0}(), autodiff = Val{true}(), standardtag = Val{true}(), concrete_jac = nothing, diff_type = Val{:forward}, linsolve = nothing, - precs = DEFAULT_PRECS, step_limiter! = trivial_limiter!) + precs = DEFAULT_PRECS, step_limiter! = trivial_limiter!, + stage_limiter! = trivial_limiter!) $Alg{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), typeof(precs), diff_type, _unwrap_val(standardtag), - _unwrap_val(concrete_jac), typeof(step_limiter!)}(linsolve, - precs, step_limiter!) + _unwrap_val(concrete_jac), typeof(step_limiter!), + typeof(stage_limiter!)}(linsolve, precs, step_limiter!, + stage_limiter!) end end diff --git a/src/caches/rosenbrock_caches.jl b/src/caches/rosenbrock_caches.jl index 9f0d3c8db6..e2eb013809 100644 --- a/src/caches/rosenbrock_caches.jl +++ b/src/caches/rosenbrock_caches.jl @@ -5,7 +5,7 @@ abstract type RosenbrockMutableCache <: OrdinaryDiffEqMutableCache end @cache mutable struct Rosenbrock23Cache{uType, rateType, uNoUnitsType, JType, WType, TabType, TFType, UFType, F, JCType, GCType, - RTolType, A, AV, StepLimiter} <: RosenbrockMutableCache + RTolType, A, AV, StepLimiter, StageLimiter} <: RosenbrockMutableCache u::uType uprev::uType k₁::rateType @@ -33,13 +33,14 @@ abstract type RosenbrockMutableCache <: OrdinaryDiffEqMutableCache end alg::A algebraic_vars::AV step_limiter!::StepLimiter + stage_limiter!::StageLimiter end TruncatedStacktraces.@truncate_stacktrace Rosenbrock23Cache 1 @cache mutable struct Rosenbrock32Cache{uType, rateType, uNoUnitsType, JType, WType, TabType, TFType, UFType, F, JCType, GCType, - RTolType, A, AV, StepLimiter} <: RosenbrockMutableCache + RTolType, A, AV, StepLimiter, StageLimiter} <: RosenbrockMutableCache u::uType uprev::uType k₁::rateType @@ -67,6 +68,7 @@ TruncatedStacktraces.@truncate_stacktrace Rosenbrock23Cache 1 alg::A algebraic_vars::AV step_limiter!::StepLimiter + stage_limiter!::StageLimiter end function alg_cache(alg::Rosenbrock23, u, rate_prototype, ::Type{uEltypeNoUnits}, @@ -110,7 +112,8 @@ function alg_cache(alg::Rosenbrock23, u, rate_prototype, ::Type{uEltypeNoUnits}, Rosenbrock23Cache(u, uprev, k₁, k₂, k₃, du1, du2, f₁, fsalfirst, fsallast, dT, J, W, tmp, atmp, weight, tab, tf, uf, linsolve_tmp, - linsolve, jac_config, grad_config, reltol, alg, algebraic_vars, alg.step_limiter!) + linsolve, jac_config, grad_config, reltol, alg, algebraic_vars, alg.step_limiter!, + alg.stage_limiter!) end function alg_cache(alg::Rosenbrock32, u, rate_prototype, ::Type{uEltypeNoUnits}, @@ -153,7 +156,7 @@ function alg_cache(alg::Rosenbrock32, u, rate_prototype, ::Type{uEltypeNoUnits}, Rosenbrock32Cache(u, uprev, k₁, k₂, k₃, du1, du2, f₁, fsalfirst, fsallast, dT, J, W, tmp, atmp, weight, tab, tf, uf, linsolve_tmp, linsolve, jac_config, - grad_config, reltol, alg, algebraic_vars, alg.step_limiter!) + grad_config, reltol, alg, algebraic_vars, alg.step_limiter!, alg.stage_limiter!) end struct Rosenbrock23ConstantCache{T, TF, UF, JType, WType, F, AD} <: @@ -232,7 +235,7 @@ end @cache mutable struct Rosenbrock33Cache{uType, rateType, uNoUnitsType, JType, WType, TabType, TFType, UFType, F, JCType, GCType, - RTolType, A, StepLimiter} <: RosenbrockMutableCache + RTolType, A, StepLimiter, StageLimiter} <: RosenbrockMutableCache u::uType uprev::uType du::rateType @@ -260,6 +263,7 @@ end reltol::RTolType alg::A step_limiter!::StepLimiter + stage_limiter!::StageLimiter end function alg_cache(alg::ROS3P, u, rate_prototype, ::Type{uEltypeNoUnits}, @@ -298,7 +302,8 @@ function alg_cache(alg::ROS3P, u, rate_prototype, ::Type{uEltypeNoUnits}, Rosenbrock33Cache(u, uprev, du, du1, du2, k1, k2, k3, k4, fsalfirst, fsallast, dT, J, W, tmp, atmp, weight, tab, tf, uf, linsolve_tmp, - linsolve, jac_config, grad_config, reltol, alg, alg.step_limiter!) + linsolve, jac_config, grad_config, reltol, alg, alg.step_limiter!, + alg.stage_limiter!) end function alg_cache(alg::ROS3P, u, rate_prototype, ::Type{uEltypeNoUnits}, @@ -316,7 +321,7 @@ function alg_cache(alg::ROS3P, u, rate_prototype, ::Type{uEltypeNoUnits}, end @cache mutable struct Rosenbrock34Cache{uType, rateType, uNoUnitsType, JType, WType, - TabType, TFType, UFType, F, JCType, GCType, StepLimiter} <: + TabType, TFType, UFType, F, JCType, GCType, StepLimiter, StageLimiter} <: RosenbrockMutableCache u::uType uprev::uType @@ -343,6 +348,7 @@ end jac_config::JCType grad_config::GCType step_limiter!::StepLimiter + stage_limiter!::StageLimiter end function alg_cache(alg::Rodas3, u, rate_prototype, ::Type{uEltypeNoUnits}, @@ -382,7 +388,8 @@ function alg_cache(alg::Rodas3, u, rate_prototype, ::Type{uEltypeNoUnits}, Rosenbrock34Cache(u, uprev, du, du1, du2, k1, k2, k3, k4, fsalfirst, fsallast, dT, J, W, tmp, atmp, weight, tab, tf, uf, linsolve_tmp, - linsolve, jac_config, grad_config, alg.step_limiter!) + linsolve, jac_config, grad_config, alg.step_limiter!, + alg.stage_limiter!) end struct Rosenbrock34ConstantCache{TF, UF, Tab, JType, WType, F} <: @@ -460,7 +467,7 @@ struct Rodas3PConstantCache{TF, UF, Tab, JType, WType, F, AD} <: OrdinaryDiffEqC end @cache mutable struct Rodas23WCache{uType, rateType, uNoUnitsType, JType, WType, TabType, - TFType, UFType, F, JCType, GCType, RTolType, A, StepLimiter} <: + TFType, UFType, F, JCType, GCType, RTolType, A, StepLimiter, StageLimiter} <: RosenbrockMutableCache u::uType uprev::uType @@ -493,10 +500,11 @@ end reltol::RTolType alg::A step_limiter!::StepLimiter + stage_limiter!::StageLimiter end @cache mutable struct Rodas3PCache{uType, rateType, uNoUnitsType, JType, WType, TabType, - TFType, UFType, F, JCType, GCType, RTolType, A, StepLimiter} <: + TFType, UFType, F, JCType, GCType, RTolType, A, StepLimiter, StageLimiter} <: RosenbrockMutableCache u::uType uprev::uType @@ -529,6 +537,7 @@ end reltol::RTolType alg::A step_limiter!::StepLimiter + stage_limiter!::StageLimiter end function alg_cache(alg::Rodas23W, u, rate_prototype, ::Type{uEltypeNoUnits}, @@ -571,7 +580,8 @@ function alg_cache(alg::Rodas23W, u, rate_prototype, ::Type{uEltypeNoUnits}, jac_config = build_jac_config(alg, f, uf, du1, uprev, u, tmp, du2) Rodas23WCache(u, uprev, dense1, dense2, dense3, du, du1, du2, k1, k2, k3, k4, k5, fsalfirst, fsallast, dT, J, W, tmp, atmp, weight, tab, tf, uf, linsolve_tmp, - linsolve, jac_config, grad_config, reltol, alg, alg.step_limiter!) + linsolve, jac_config, grad_config, reltol, alg, alg.step_limiter!, + alg.stage_limiter!) end TruncatedStacktraces.@truncate_stacktrace Rodas23WCache 1 @@ -615,7 +625,8 @@ function alg_cache(alg::Rodas3P, u, rate_prototype, ::Type{uEltypeNoUnits}, jac_config = build_jac_config(alg, f, uf, du1, uprev, u, tmp, du2) Rodas3PCache(u, uprev, dense1, dense2, dense3, du, du1, du2, k1, k2, k3, k4, k5, fsalfirst, fsallast, dT, J, W, tmp, atmp, weight, tab, tf, uf, linsolve_tmp, - linsolve, jac_config, grad_config, reltol, alg, alg.step_limiter!) + linsolve, jac_config, grad_config, reltol, alg, alg.step_limiter!, + alg.stage_limiter!) end TruncatedStacktraces.@truncate_stacktrace Rodas3PCache 1 @@ -663,7 +674,7 @@ struct Rodas4ConstantCache{TF, UF, Tab, JType, WType, F, AD} <: OrdinaryDiffEqCo end @cache mutable struct Rodas4Cache{uType, rateType, uNoUnitsType, JType, WType, TabType, - TFType, UFType, F, JCType, GCType, RTolType, A, StepLimiter} <: + TFType, UFType, F, JCType, GCType, RTolType, A, StepLimiter, StageLimiter} <: RosenbrockMutableCache u::uType uprev::uType @@ -696,6 +707,7 @@ end reltol::RTolType alg::A step_limiter!::StepLimiter + stage_limiter!::StageLimiter end function alg_cache(alg::Rodas4, u, rate_prototype, ::Type{uEltypeNoUnits}, @@ -739,7 +751,8 @@ function alg_cache(alg::Rodas4, u, rate_prototype, ::Type{uEltypeNoUnits}, Rodas4Cache(u, uprev, dense1, dense2, du, du1, du2, k1, k2, k3, k4, k5, k6, fsalfirst, fsallast, dT, J, W, tmp, atmp, weight, tab, tf, uf, linsolve_tmp, - linsolve, jac_config, grad_config, reltol, alg, alg.step_limiter!) + linsolve, jac_config, grad_config, reltol, alg, alg.step_limiter!, + alg.stage_limiter!) end TruncatedStacktraces.@truncate_stacktrace Rodas4Cache 1 @@ -800,7 +813,8 @@ function alg_cache(alg::Rodas42, u, rate_prototype, ::Type{uEltypeNoUnits}, Rodas4Cache(u, uprev, dense1, dense2, du, du1, du2, k1, k2, k3, k4, k5, k6, fsalfirst, fsallast, dT, J, W, tmp, atmp, weight, tab, tf, uf, linsolve_tmp, - linsolve, jac_config, grad_config, reltol, alg, alg.step_limiter!) + linsolve, jac_config, grad_config, reltol, alg, alg.step_limiter!, + alg.stage_limiter!) end function alg_cache(alg::Rodas42, u, rate_prototype, ::Type{uEltypeNoUnits}, @@ -859,7 +873,8 @@ function alg_cache(alg::Rodas4P, u, rate_prototype, ::Type{uEltypeNoUnits}, Rodas4Cache(u, uprev, dense1, dense2, du, du1, du2, k1, k2, k3, k4, k5, k6, fsalfirst, fsallast, dT, J, W, tmp, atmp, weight, tab, tf, uf, linsolve_tmp, - linsolve, jac_config, grad_config, reltol, alg, alg.step_limiter!) + linsolve, jac_config, grad_config, reltol, alg, alg.step_limiter!, + alg.stage_limiter!) end function alg_cache(alg::Rodas4P, u, rate_prototype, ::Type{uEltypeNoUnits}, @@ -918,7 +933,8 @@ function alg_cache(alg::Rodas4P2, u, rate_prototype, ::Type{uEltypeNoUnits}, Rodas4Cache(u, uprev, dense1, dense2, du, du1, du2, k1, k2, k3, k4, k5, k6, fsalfirst, fsallast, dT, J, W, tmp, atmp, weight, tab, tf, uf, linsolve_tmp, - linsolve, jac_config, grad_config, reltol, alg, alg.step_limiter!) + linsolve, jac_config, grad_config, reltol, alg, alg.step_limiter!, + alg.stage_limiter!) end function alg_cache(alg::Rodas4P2, u, rate_prototype, ::Type{uEltypeNoUnits}, @@ -951,7 +967,7 @@ end @cache mutable struct Rosenbrock5Cache{ uType, rateType, uNoUnitsType, JType, WType, TabType, - TFType, UFType, F, JCType, GCType, RTolType, A, StepLimiter} <: + TFType, UFType, F, JCType, GCType, RTolType, A, StepLimiter, StageLimiter} <: RosenbrockMutableCache u::uType uprev::uType @@ -987,6 +1003,7 @@ end reltol::RTolType alg::A step_limiter!::StepLimiter + stage_limiter!::StageLimiter end TruncatedStacktraces.@truncate_stacktrace Rosenbrock5Cache 1 @@ -1036,7 +1053,8 @@ function alg_cache(alg::Rodas5, u, rate_prototype, ::Type{uEltypeNoUnits}, k5, k6, k7, k8, fsalfirst, fsallast, dT, J, W, tmp, atmp, weight, tab, tf, uf, linsolve_tmp, - linsolve, jac_config, grad_config, reltol, alg, alg.step_limiter!) + linsolve, jac_config, grad_config, reltol, alg, alg.step_limiter!, + alg.stage_limiter!) end function alg_cache(alg::Rodas5, u, rate_prototype, ::Type{uEltypeNoUnits}, @@ -1099,7 +1117,8 @@ function alg_cache( k5, k6, k7, k8, fsalfirst, fsallast, dT, J, W, tmp, atmp, weight, tab, tf, uf, linsolve_tmp, - linsolve, jac_config, grad_config, reltol, alg, alg.step_limiter!) + linsolve, jac_config, grad_config, reltol, alg, alg.step_limiter!, + alg.stage_limiter!) end function alg_cache( diff --git a/src/perform_step/rosenbrock_perform_step.jl b/src/perform_step/rosenbrock_perform_step.jl index c3236a5b53..0521dc9779 100644 --- a/src/perform_step/rosenbrock_perform_step.jl +++ b/src/perform_step/rosenbrock_perform_step.jl @@ -27,7 +27,7 @@ end @muladd function perform_step!(integrator, cache::Rosenbrock23Cache, repeat_step = false) @unpack t, dt, uprev, u, f, p, opts = integrator - @unpack k₁, k₂, k₃, du1, du2, f₁, fsalfirst, fsallast, dT, J, W, tmp, uf, tf, linsolve_tmp, jac_config, atmp, weight, step_limiter! = cache + @unpack k₁, k₂, k₃, du1, du2, f₁, fsalfirst, fsallast, dT, J, W, tmp, uf, tf, linsolve_tmp, jac_config, atmp, weight, stage_limiter!, step_limiter! = cache @unpack c₃₂, d = cache.tab # Assignments @@ -68,6 +68,7 @@ end integrator.stats.nsolve += 1 @.. broadcast=false u=uprev + dto2 * k₁ + stage_limiter!(u, integrator, p, t + dto2) f(f₁, u, p, t + dto2) integrator.stats.nf += 1 @@ -88,7 +89,7 @@ end @.. broadcast=false k₂+=k₁ @.. broadcast=false u=uprev + dt * k₂ - + stage_limiter!(u, integrator, p, t + dt) step_limiter!(u, integrator, p, t + dt) if integrator.opts.adaptive @@ -138,7 +139,7 @@ end @muladd function perform_step!(integrator, cache::Rosenbrock23Cache{<:Array}, repeat_step = false) @unpack t, dt, uprev, u, f, p, opts = integrator - @unpack k₁, k₂, k₃, du1, du2, f₁, fsalfirst, fsallast, dT, J, W, tmp, uf, tf, linsolve_tmp, jac_config, atmp, weight, step_limiter! = cache + @unpack k₁, k₂, k₃, du1, du2, f₁, fsalfirst, fsallast, dT, J, W, tmp, uf, tf, linsolve_tmp, jac_config, atmp, weight, stage_limiter!, step_limiter! = cache @unpack c₃₂, d = cache.tab # Assignments @@ -182,6 +183,7 @@ end @inbounds @simd ivdep for i in eachindex(u) u[i] = uprev[i] + dto2 * k₁[i] end + stage_limiter!(u, integrator, p, t + dto2) f(f₁, u, p, t + dto2) integrator.stats.nf += 1 @@ -209,7 +211,7 @@ end @inbounds @simd ivdep for i in eachindex(u) u[i] = uprev[i] + dt * k₂[i] end - + stage_limiter!(u, integrator, p, t + dt) step_limiter!(u, integrator, p, t + dt) if integrator.opts.adaptive @@ -278,7 +280,7 @@ end @muladd function perform_step!(integrator, cache::Rosenbrock32Cache, repeat_step = false) @unpack t, dt, uprev, u, f, p, opts = integrator - @unpack k₁, k₂, k₃, du1, du2, f₁, fsalfirst, fsallast, dT, J, W, tmp, uf, tf, linsolve_tmp, jac_config, atmp, weight, step_limiter! = cache + @unpack k₁, k₂, k₃, du1, du2, f₁, fsalfirst, fsallast, dT, J, W, tmp, uf, tf, linsolve_tmp, jac_config, atmp, weight, stage_limiter!, step_limiter! = cache @unpack c₃₂, d = cache.tab # Assignments @@ -319,6 +321,7 @@ end integrator.stats.nsolve += 1 @.. broadcast=false u=uprev + dto2 * k₁ + stage_limiter!(u, integrator, p, t + dto2) f(f₁, u, p, t + dto2) integrator.stats.nf += 1 @@ -339,6 +342,7 @@ end @.. broadcast=false k₂+=k₁ @.. broadcast=false tmp=uprev + dt * k₂ + stage_limiter!(u, integrator, p, t + dt) f(fsallast, tmp, p, t + dt) integrator.stats.nf += 1 @@ -634,7 +638,7 @@ end @muladd function perform_step!(integrator, cache::Rosenbrock33Cache, repeat_step = false) @unpack t, dt, uprev, u, f, p = integrator - @unpack du, du1, du2, fsalfirst, fsallast, k1, k2, k3, dT, J, W, uf, tf, linsolve_tmp, jac_config, atmp, weight, step_limiter! = cache + @unpack du, du1, du2, fsalfirst, fsallast, k1, k2, k3, dT, J, W, uf, tf, linsolve_tmp, jac_config, atmp, weight, stage_limiter!, step_limiter! = cache @unpack a21, a31, a32, C21, C31, C32, b1, b2, b3, btilde1, btilde2, btilde3, gamma, c2, c3, d1, d2, d3 = cache.tab # Assignments @@ -676,6 +680,7 @@ end integrator.stats.nsolve += 1 @.. broadcast=false u=uprev + a21 * k1 + stage_limiter!(u, integrator, p, t + c2 * dt) f(du, u, p, t + c2 * dt) integrator.stats.nf += 1 @@ -696,6 +701,7 @@ end integrator.stats.nsolve += 1 @.. broadcast=false u=uprev + a31 * k1 + a32 * k2 + stage_limiter!(u, integrator, p, t + c3 * dt) f(du, u, p, t + c3 * dt) integrator.stats.nf += 1 @@ -822,7 +828,7 @@ end @muladd function perform_step!(integrator, cache::Rosenbrock34Cache, repeat_step = false) @unpack t, dt, uprev, u, f, p = integrator - @unpack du, du1, du2, fsalfirst, fsallast, k1, k2, k3, k4, dT, J, W, uf, tf, linsolve_tmp, jac_config, atmp, weight, step_limiter! = cache + @unpack du, du1, du2, fsalfirst, fsallast, k1, k2, k3, k4, dT, J, W, uf, tf, linsolve_tmp, jac_config, atmp, weight, stage_limiter!, step_limiter! = cache @unpack a21, a31, a32, a41, a42, a43, C21, C31, C32, C41, C42, C43, b1, b2, b3, b4, btilde1, btilde2, btilde3, btilde4, gamma, c2, c3, d1, d2, d3, d4 = cache.tab # Assignments @@ -890,6 +896,7 @@ end integrator.stats.nsolve += 1 @.. broadcast=false u=uprev + a31 * k1 + a32 * k2 + stage_limiter!(u, integrator, p, t + c3 * dt) f(du, u, p, t + c3 * dt) integrator.stats.nf += 1 @@ -906,6 +913,7 @@ end @.. broadcast=false veck3=-vecu integrator.stats.nsolve += 1 @.. broadcast=false u=uprev + a41 * k1 + a42 * k2 + a43 * k3 + stage_limiter!(u, integrator, p, t + dt) f(du, u, p, t + dt) #-- c4 = 1 integrator.stats.nf += 1 @@ -1123,7 +1131,7 @@ end @muladd function perform_step!( integrator, cache::Union{Rodas23WCache, Rodas3PCache}, repeat_step = false) @unpack t, dt, uprev, u, f, p = integrator - @unpack du, du1, du2, dT, J, W, uf, tf, k1, k2, k3, k4, k5, linsolve_tmp, jac_config, atmp, weight, step_limiter! = cache + @unpack du, du1, du2, dT, J, W, uf, tf, k1, k2, k3, k4, k5, linsolve_tmp, jac_config, atmp, weight, stage_limiter!, step_limiter! = cache @unpack a21, a41, a42, a43, C21, C31, C32, C41, C42, C43, C51, C52, C53, C54, gamma, c2, c3, d1, d2, d3 = cache.tab # Assignments @@ -1173,6 +1181,7 @@ end integrator.stats.nsolve += 1 @.. broadcast=false u=uprev + a21 * k1 + stage_limiter!(u, integrator, p, t + c2 * dt) f(du, u, p, t + c2 * dt) integrator.stats.nf += 1 @@ -1202,6 +1211,7 @@ end integrator.stats.nsolve += 1 @.. broadcast=false u=uprev + a41 * k1 + a42 * k2 + a43 * k3 + stage_limiter!(u, integrator, p, t + c2 * dt) f(du, u, p, t + dt) integrator.stats.nf += 1 @@ -1679,7 +1689,7 @@ end @muladd function perform_step!(integrator, cache::Rodas4Cache, repeat_step = false) @unpack t, dt, uprev, u, f, p = integrator - @unpack du, du1, du2, dT, J, W, uf, tf, k1, k2, k3, k4, k5, k6, linsolve_tmp, jac_config, atmp, weight, step_limiter! = cache + @unpack du, du1, du2, dT, J, W, uf, tf, k1, k2, k3, k4, k5, k6, linsolve_tmp, jac_config, atmp, weight, stage_limiter!, step_limiter! = cache @unpack a21, a31, a32, a41, a42, a43, a51, a52, a53, a54, C21, C31, C32, C41, C42, C43, C51, C52, C53, C54, C61, C62, C63, C64, C65, gamma, c2, c3, c4, d1, d2, d3, d4 = cache.tab # Assignments @@ -1735,6 +1745,7 @@ end integrator.stats.nsolve += 1 @.. broadcast=false u=uprev + a21 * k1 + stage_limiter!(u, integrator, p, t + c2 * dt) f(du, u, p, t + c2 * dt) integrator.stats.nf += 1 @@ -1751,6 +1762,7 @@ end integrator.stats.nsolve += 1 @.. broadcast=false u=uprev + a31 * k1 + a32 * k2 + stage_limiter!(u, integrator, p, t + c3 * dt) f(du, u, p, t + c3 * dt) integrator.stats.nf += 1 @@ -1767,6 +1779,7 @@ end integrator.stats.nsolve += 1 @.. broadcast=false u=uprev + a41 * k1 + a42 * k2 + a43 * k3 + stage_limiter!(u, integrator, p, t + c4 * dt) f(du, u, p, t + c4 * dt) integrator.stats.nf += 1 @@ -1784,6 +1797,7 @@ end integrator.stats.nsolve += 1 @.. broadcast=false u=uprev + a51 * k1 + a52 * k2 + a53 * k3 + a54 * k4 + stage_limiter!(u, integrator, p, t + dt) f(du, u, p, t + dt) integrator.stats.nf += 1 @@ -1840,7 +1854,7 @@ end @muladd function perform_step!(integrator, cache::Rodas4Cache{<:Array}, repeat_step = false) @unpack t, dt, uprev, u, f, p = integrator - @unpack du, du1, du2, dT, J, W, uf, tf, k1, k2, k3, k4, k5, k6, linsolve_tmp, jac_config, atmp, weight, step_limiter! = cache + @unpack du, du1, du2, dT, J, W, uf, tf, k1, k2, k3, k4, k5, k6, linsolve_tmp, jac_config, atmp, weight, stage_limiter!, step_limiter! = cache @unpack a21, a31, a32, a41, a42, a43, a51, a52, a53, a54, C21, C31, C32, C41, C42, C43, C51, C52, C53, C54, C61, C62, C63, C64, C65, gamma, c2, c3, c4, d1, d2, d3, d4 = cache.tab # Assignments @@ -1899,6 +1913,7 @@ end @inbounds @simd ivdep for i in eachindex(u) u[i] = uprev[i] + a21 * k1[i] end + stage_limiter!(u, integrator, p, t + c2 * dt) f(du, u, p, t + c2 * dt) integrator.stats.nf += 1 @@ -1926,6 +1941,7 @@ end @inbounds @simd ivdep for i in eachindex(u) u[i] = uprev[i] + a31 * k1[i] + a32 * k2[i] end + stage_limiter!(u, integrator, p, t + c3 * dt) f(du, u, p, t + c3 * dt) integrator.stats.nf += 1 @@ -1953,6 +1969,7 @@ end @inbounds @simd ivdep for i in eachindex(u) u[i] = uprev[i] + a41 * k1[i] + a42 * k2[i] + a43 * k3[i] end + stage_limiter!(u, integrator, p, t + c4 * dt) f(du, u, p, t + c4 * dt) integrator.stats.nf += 1 @@ -1981,6 +1998,7 @@ end @inbounds @simd ivdep for i in eachindex(u) u[i] = uprev[i] + a51 * k1[i] + a52 * k2[i] + a53 * k3[i] + a54 * k4[i] end + stage_limiter!(u, integrator, p, t + dt) f(du, u, p, t + dt) integrator.stats.nf += 1 @@ -2009,6 +2027,7 @@ end @inbounds @simd ivdep for i in eachindex(u) u[i] += k5[i] end + stage_limiter!(u, integrator, p, t + dt) f(du, u, p, t + dt) integrator.stats.nf += 1 @@ -2284,7 +2303,7 @@ end @muladd function perform_step!(integrator, cache::Rosenbrock5Cache, repeat_step = false) @unpack t, dt, uprev, u, f, p = integrator - @unpack du, du1, du2, k1, k2, k3, k4, k5, k6, k7, k8, dT, J, W, uf, tf, linsolve_tmp, jac_config, atmp, weight, step_limiter! = cache + @unpack du, du1, du2, k1, k2, k3, k4, k5, k6, k7, k8, dT, J, W, uf, tf, linsolve_tmp, jac_config, atmp, weight, stage_limiter!, step_limiter! = cache @unpack a21, a31, a32, a41, a42, a43, a51, a52, a53, a54, a61, a62, a63, a64, a65, C21, C31, C32, C41, C42, C43, C51, C52, C53, C54, C61, C62, C63, C64, C65, C71, C72, C73, C74, C75, C76, C81, C82, C83, C84, C85, C86, C87, gamma, d1, d2, d3, d4, d5, c2, c3, c4, c5 = cache.tab # Assignments @@ -2356,6 +2375,7 @@ end integrator.stats.nsolve += 1 @.. broadcast=false u=uprev + a21 * k1 + stage_limiter!(u, integrator, p, t + c2 * dt) f(du, u, p, t + c2 * dt) integrator.stats.nf += 1 @@ -2373,6 +2393,7 @@ end integrator.stats.nsolve += 1 @.. broadcast=false u=uprev + a31 * k1 + a32 * k2 + stage_limiter!(u, integrator, p, t + c3 * dt) f(du, u, p, t + c3 * dt) integrator.stats.nf += 1 @@ -2390,6 +2411,7 @@ end integrator.stats.nsolve += 1 @.. broadcast=false u=uprev + a41 * k1 + a42 * k2 + a43 * k3 + stage_limiter!(u, integrator, p, t + c4 * dt) f(du, u, p, t + c4 * dt) integrator.stats.nf += 1 @@ -2408,6 +2430,7 @@ end integrator.stats.nsolve += 1 @.. broadcast=false u=uprev + a51 * k1 + a52 * k2 + a53 * k3 + a54 * k4 + stage_limiter!(u, integrator, p, t + c5 * dt) f(du, u, p, t + c5 * dt) integrator.stats.nf += 1 @@ -2426,6 +2449,7 @@ end integrator.stats.nsolve += 1 @.. broadcast=false u=uprev + a61 * k1 + a62 * k2 + a63 * k3 + a64 * k4 + a65 * k5 + stage_limiter!(u, integrator, p, t + dt) f(du, u, p, t + dt) integrator.stats.nf += 1 @@ -2445,6 +2469,7 @@ end integrator.stats.nsolve += 1 u .+= k6 + stage_limiter!(u, integrator, p, t + dt) f(du, u, p, t + dt) integrator.stats.nf += 1 @@ -2530,7 +2555,7 @@ end @muladd function perform_step!(integrator, cache::Rosenbrock5Cache{<:Array}, repeat_step = false) @unpack t, dt, uprev, u, f, p = integrator - @unpack du, du1, du2, k1, k2, k3, k4, k5, k6, k7, k8, dT, J, W, uf, tf, linsolve_tmp, jac_config, atmp, weight, step_limiter! = cache + @unpack du, du1, du2, k1, k2, k3, k4, k5, k6, k7, k8, dT, J, W, uf, tf, linsolve_tmp, jac_config, atmp, weight, stage_limiter!, step_limiter! = cache @unpack a21, a31, a32, a41, a42, a43, a51, a52, a53, a54, a61, a62, a63, a64, a65, C21, C31, C32, C41, C42, C43, C51, C52, C53, C54, C61, C62, C63, C64, C65, C71, C72, C73, C74, C75, C76, C81, C82, C83, C84, C85, C86, C87, gamma, d1, d2, d3, d4, d5, c2, c3, c4, c5 = cache.tab # Assignments @@ -2607,6 +2632,7 @@ end @inbounds @simd ivdep for i in eachindex(u) u[i] = uprev[i] + a21 * k1[i] end + stage_limiter!(u, integrator, p, t + c2 * dt) f(du, u, p, t + c2 * dt) integrator.stats.nf += 1 @@ -2635,6 +2661,7 @@ end @inbounds @simd ivdep for i in eachindex(u) u[i] = uprev[i] + a31 * k1[i] + a32 * k2[i] end + stage_limiter!(u, integrator, p, t + c3 * dt) f(du, u, p, t + c3 * dt) integrator.stats.nf += 1 @@ -2662,6 +2689,7 @@ end @inbounds @simd ivdep for i in eachindex(u) u[i] = uprev[i] + a41 * k1[i] + a42 * k2[i] + a43 * k3[i] end + stage_limiter!(u, integrator, p, t + c4 * dt) f(du, u, p, t + c4 * dt) integrator.stats.nf += 1 @@ -2690,6 +2718,7 @@ end @inbounds @simd ivdep for i in eachindex(u) u[i] = uprev[i] + a51 * k1[i] + a52 * k2[i] + a53 * k3[i] + a54 * k4[i] end + stage_limiter!(u, integrator, p, t + c5 * dt) f(du, u, p, t + c5 * dt) integrator.stats.nf += 1 @@ -2720,6 +2749,7 @@ end u[i] = uprev[i] + a61 * k1[i] + a62 * k2[i] + a63 * k3[i] + a64 * k4[i] + a65 * k5[i] end + stage_limiter!(u, integrator, p, t + dt) f(du, u, p, t + dt) integrator.stats.nf += 1 @@ -2749,6 +2779,7 @@ end @inbounds @simd ivdep for i in eachindex(u) u[i] += k6[i] end + stage_limiter!(u, integrator, p, t + dt) f(du, u, p, t + dt) integrator.stats.nf += 1 @@ -2778,6 +2809,7 @@ end @inbounds @simd ivdep for i in eachindex(u) u[i] += k7[i] end + stage_limiter!(u, integrator, p, t + dt) f(du, u, p, t + dt) integrator.stats.nf += 1