diff --git a/docs/src/massmatrixdae/Rosenbrock.md b/docs/src/massmatrixdae/Rosenbrock.md index d62013d132..aadbf41b66 100644 --- a/docs/src/massmatrixdae/Rosenbrock.md +++ b/docs/src/massmatrixdae/Rosenbrock.md @@ -11,7 +11,7 @@ the former is more efficient, but the latter is more reliable. For larger systems look at multistep methods. !!! warn - + In order to use OrdinaryDiffEqRosenbrock with DAEs that require a non-trivial consistent initialization, a nonlinear solver is required and thus `using OrdinaryDiffEqNonlinearSolve` is required or you must pass an `initializealg` diff --git a/lib/OrdinaryDiffEqBDF/src/bdf_perform_step.jl b/lib/OrdinaryDiffEqBDF/src/bdf_perform_step.jl index 4f527bac7a..bca93c4d01 100644 --- a/lib/OrdinaryDiffEqBDF/src/bdf_perform_step.jl +++ b/lib/OrdinaryDiffEqBDF/src/bdf_perform_step.jl @@ -1113,7 +1113,8 @@ function perform_step!(integrator, cache::FBDFConstantCache{max_order}, end tmp = -uprev * bdf_coeffs[k, 2] for i in 1:(k - 1) - tmp = @.. tmp - $(_reshape(view(u_corrector, :, i), axes(u))) * bdf_coeffs[k, i + 2] + tmp = @.. tmp - + $(_reshape(view(u_corrector, :, i), axes(u))) * bdf_coeffs[k, i + 2] end end @@ -1170,7 +1171,9 @@ function perform_step!(integrator, cache::FBDFConstantCache{max_order}, terk *= abs(dt^(k)) else for i in 2:(k + 1) - terk = @.. terk + fd_weights[i, k + 1] * $(_reshape(view(u_history, :, i - 1), axes(u))) + terk = @.. terk + + fd_weights[i, k + 1] * + $(_reshape(view(u_history, :, i - 1), axes(u))) end terk *= abs(dt^(k)) end diff --git a/lib/OrdinaryDiffEqBDF/src/controllers.jl b/lib/OrdinaryDiffEqBDF/src/controllers.jl index e992aaeb1d..cd65919803 100644 --- a/lib/OrdinaryDiffEqBDF/src/controllers.jl +++ b/lib/OrdinaryDiffEqBDF/src/controllers.jl @@ -217,7 +217,8 @@ function choose_order!(alg::FBDF, integrator, terk_tmp = similar(u) @.. terk_tmp = fd_weights[k - 2, 1] * _vec(u) for i in 2:(k - 2) - @.. terk_tmp += fd_weights[i, k - 2] * $(_reshape(view(u_history, :, i - 1), axes(u))) + @.. terk_tmp += fd_weights[i, k - 2] * + $(_reshape(view(u_history, :, i - 1), axes(u))) end @.. terk_tmp *= abs(dt^(k - 2)) end diff --git a/lib/OrdinaryDiffEqCore/src/caches/basic_caches.jl b/lib/OrdinaryDiffEqCore/src/caches/basic_caches.jl index bd4edd7d15..f836d61b9b 100644 --- a/lib/OrdinaryDiffEqCore/src/caches/basic_caches.jl +++ b/lib/OrdinaryDiffEqCore/src/caches/basic_caches.jl @@ -51,7 +51,8 @@ function get_fsalfirstlast(cache::DefaultCache, u) (cache.u, cache.u) # will be overwritten by the cache choice end -function ismutablecache(cache::DefaultCache{T1, T2, T3, T4, T5, T6, A, F, uType}) where {T1, T2, T3, T4, T5, T6, A, F, uType} +function ismutablecache(cache::DefaultCache{ + T1, T2, T3, T4, T5, T6, A, F, uType}) where {T1, T2, T3, T4, T5, T6, A, F, uType} T1 <: OrdinaryDiffEqMutableCache end diff --git a/lib/OrdinaryDiffEqCore/src/initialize_dae.jl b/lib/OrdinaryDiffEqCore/src/initialize_dae.jl index 7a8fc14b28..c352cb3ddf 100644 --- a/lib/OrdinaryDiffEqCore/src/initialize_dae.jl +++ b/lib/OrdinaryDiffEqCore/src/initialize_dae.jl @@ -121,20 +121,23 @@ function default_nlsolve( end function default_nlsolve( - ::Nothing, isinplace::Val{false}, u::Nothing, ::NonlinearProblem, autodiff = false) -nothing + ::Nothing, isinplace::Val{false}, u::Nothing, ::NonlinearProblem, autodiff = false) + nothing end function default_nlsolve( - ::Nothing, isinplace::Val{false}, u::Nothing, ::NonlinearLeastSquaresProblem, autodiff = false) -nothing + ::Nothing, isinplace::Val{false}, u::Nothing, + ::NonlinearLeastSquaresProblem, autodiff = false) + nothing end -function OrdinaryDiffEqCore.default_nlsolve(::Nothing, isinplace, u, ::NonlinearProblem, autodiff = false) +function OrdinaryDiffEqCore.default_nlsolve( + ::Nothing, isinplace, u, ::NonlinearProblem, autodiff = false) error("This ODE requires a DAE initialization and thus a nonlinear solve but no nonlinear solve has been loaded. To solve this problem, do `using OrdinaryDiffEqNonlinearSolve` or pass a custom `nlsolve` choice into the `initializealg`.") end -function OrdinaryDiffEqCore.default_nlsolve(::Nothing, isinplace, u, ::NonlinearLeastSquaresProblem, autodiff = false) +function OrdinaryDiffEqCore.default_nlsolve( + ::Nothing, isinplace, u, ::NonlinearLeastSquaresProblem, autodiff = false) error("This ODE requires a DAE initialization and thus a nonlinear solve but no nonlinear solve has been loaded. To solve this problem, do `using OrdinaryDiffEqNonlinearSolve` or pass a custom `nlsolve` choice into the `initializealg`.") end @@ -179,12 +182,13 @@ end ## CheckInit struct CheckInitFailureError <: Exception - normresid - abstol + normresid::Any + abstol::Any end function Base.showerror(io::IO, e::CheckInitFailureError) - print(io, "CheckInit specified but initialization not satisifed. normresid = $(e.normresid) > abstol = $(e.abstol)") + print(io, + "CheckInit specified but initialization not satisifed. normresid = $(e.normresid) > abstol = $(e.abstol)") end function _initialize_dae!(integrator, prob::ODEProblem, alg::CheckInit, @@ -202,7 +206,7 @@ function _initialize_dae!(integrator, prob::ODEProblem, alg::CheckInit, tmp .= ArrayInterface.restructure(tmp, algebraic_eqs .* _vec(tmp)) normresid = integrator.opts.internalnorm(tmp, t) - if normresid > integrator.opts.abstol + if normresid > integrator.opts.abstol throw(CheckInitFailureError(normresid, integrator.opts.abstol)) end end @@ -221,7 +225,7 @@ function _initialize_dae!(integrator, prob::ODEProblem, alg::CheckInit, resid = _vec(du)[algebraic_eqs] normresid = integrator.opts.internalnorm(resid, t) - if normresid > integrator.opts.abstol + if normresid > integrator.opts.abstol throw(CheckInitFailureError(normresid, integrator.opts.abstol)) end end @@ -234,7 +238,7 @@ function _initialize_dae!(integrator, prob::DAEProblem, f(resid, integrator.du, u0, p, t) normresid = integrator.opts.internalnorm(resid, t) - if normresid > integrator.opts.abstol + if normresid > integrator.opts.abstol throw(CheckInitFailureError(normresid, integrator.opts.abstol)) end end @@ -252,7 +256,7 @@ function _initialize_dae!(integrator, prob::DAEProblem, resid = f(integrator.du, u0, p, t) normresid = integrator.opts.internalnorm(resid, t) - if normresid > integrator.opts.abstol + if normresid > integrator.opts.abstol throw(CheckInitFailureError(normresid, integrator.opts.abstol)) end -end \ No newline at end of file +end diff --git a/lib/OrdinaryDiffEqCore/src/integrators/integrator_interface.jl b/lib/OrdinaryDiffEqCore/src/integrators/integrator_interface.jl index e7fa75fbd8..b629a81141 100644 --- a/lib/OrdinaryDiffEqCore/src/integrators/integrator_interface.jl +++ b/lib/OrdinaryDiffEqCore/src/integrators/integrator_interface.jl @@ -38,9 +38,10 @@ end function DiffEqBase.reeval_internals_due_to_modification!( integrator::ODEIntegrator, continuous_modification = true; callback_initializealg = nothing) - if integrator.isdae - DiffEqBase.initialize_dae!(integrator, isnothing(callback_initializealg) ? integrator.initializealg : callback_initializealg) + DiffEqBase.initialize_dae!(integrator, + isnothing(callback_initializealg) ? integrator.initializealg : + callback_initializealg) update_uprev!(integrator) end diff --git a/lib/OrdinaryDiffEqCore/src/interp_func.jl b/lib/OrdinaryDiffEqCore/src/interp_func.jl index 3d737af0b5..ed989a15b9 100644 --- a/lib/OrdinaryDiffEqCore/src/interp_func.jl +++ b/lib/OrdinaryDiffEqCore/src/interp_func.jl @@ -76,7 +76,8 @@ end function strip_cache(cache) if hasfield(typeof(cache), :jac_config) || hasfield(typeof(cache), :grad_config) || - hasfield(typeof(cache), :nlsolver) || hasfield(typeof(cache), :tf) || hasfield(typeof(cache), :uf) + hasfield(typeof(cache), :nlsolver) || hasfield(typeof(cache), :tf) || + hasfield(typeof(cache), :uf) fieldnums = length(fieldnames(typeof(cache))) noth_list = fill(nothing, fieldnums) cache_type_name = Base.typename(typeof(cache)).wrapper diff --git a/lib/OrdinaryDiffEqDefault/test/default_solver_tests.jl b/lib/OrdinaryDiffEqDefault/test/default_solver_tests.jl index 4191dbf4d8..8676b4a57a 100644 --- a/lib/OrdinaryDiffEqDefault/test/default_solver_tests.jl +++ b/lib/OrdinaryDiffEqDefault/test/default_solver_tests.jl @@ -1,4 +1,5 @@ -using OrdinaryDiffEqDefault, OrdinaryDiffEqTsit5, OrdinaryDiffEqVerner, OrdinaryDiffEqRosenbrock, OrdinaryDiffEqBDF +using OrdinaryDiffEqDefault, OrdinaryDiffEqTsit5, OrdinaryDiffEqVerner, + OrdinaryDiffEqRosenbrock, OrdinaryDiffEqBDF using Test, LinearSolve, LinearAlgebra, SparseArrays, StaticArrays f_2dlinear = (du, u, p, t) -> (@. du = p * u) diff --git a/lib/OrdinaryDiffEqDifferentiation/src/alg_utils.jl b/lib/OrdinaryDiffEqDifferentiation/src/alg_utils.jl index cfc14219c5..50e605e409 100644 --- a/lib/OrdinaryDiffEqDifferentiation/src/alg_utils.jl +++ b/lib/OrdinaryDiffEqDifferentiation/src/alg_utils.jl @@ -76,4 +76,4 @@ end @generated function pick_static_chunksize(::Val{chunksize}) where {chunksize} x = ForwardDiff.pickchunksize(chunksize) :(Val{$x}()) -end \ No newline at end of file +end diff --git a/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl b/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl index e351ab54c3..7951f1fc72 100644 --- a/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl +++ b/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl @@ -12,7 +12,7 @@ struct StaticWOperator{isinv, T, F} <: AbstractSciMLOperator{T} # doing to how StaticArrays and StaticArraysCore are split up StaticArrays.LU(LowerTriangular(W), UpperTriangular(W), SVector{n}(1:n)) else - lu(W, check=false) + lu(W, check = false) end # when constructing W for the first time for the type # inv(W) can be singular @@ -938,28 +938,28 @@ function LinearSolve.init_cacheval( end for alg in [LinearSolve.AppleAccelerateLUFactorization, - LinearSolve.BunchKaufmanFactorization, - LinearSolve.CHOLMODFactorization, - LinearSolve.CholeskyFactorization, - LinearSolve.CudaOffloadFactorization, - LinearSolve.DiagonalFactorization, - LinearSolve.FastLUFactorization, - LinearSolve.FastQRFactorization, - LinearSolve.GenericFactorization, - LinearSolve.GenericLUFactorization, - LinearSolve.KLUFactorization, - LinearSolve.LDLtFactorization, - LinearSolve.LUFactorization, - LinearSolve.MKLLUFactorization, - LinearSolve.MetalLUFactorization, - LinearSolve.NormalBunchKaufmanFactorization, - LinearSolve.NormalCholeskyFactorization, - LinearSolve.QRFactorization, - LinearSolve.RFLUFactorization, - LinearSolve.SVDFactorization, - LinearSolve.SimpleLUFactorization, - LinearSolve.SparspakFactorization, - LinearSolve.UMFPACKFactorization] + LinearSolve.BunchKaufmanFactorization, + LinearSolve.CHOLMODFactorization, + LinearSolve.CholeskyFactorization, + LinearSolve.CudaOffloadFactorization, + LinearSolve.DiagonalFactorization, + LinearSolve.FastLUFactorization, + LinearSolve.FastQRFactorization, + LinearSolve.GenericFactorization, + LinearSolve.GenericLUFactorization, + LinearSolve.KLUFactorization, + LinearSolve.LDLtFactorization, + LinearSolve.LUFactorization, + LinearSolve.MKLLUFactorization, + LinearSolve.MetalLUFactorization, + LinearSolve.NormalBunchKaufmanFactorization, + LinearSolve.NormalCholeskyFactorization, + LinearSolve.QRFactorization, + LinearSolve.RFLUFactorization, + LinearSolve.SVDFactorization, + LinearSolve.SimpleLUFactorization, + LinearSolve.SparspakFactorization, + LinearSolve.UMFPACKFactorization] @eval function LinearSolve.init_cacheval(alg::$alg, A::WOperator, b, u, Pl, Pr, maxiters::Int, abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) @@ -1003,4 +1003,4 @@ function resize_J_W!(cache, integrator, i) end nothing -end \ No newline at end of file +end diff --git a/lib/OrdinaryDiffEqExtrapolation/src/extrapolation_perform_step.jl b/lib/OrdinaryDiffEqExtrapolation/src/extrapolation_perform_step.jl index 5c33ce5ed3..2c8be5cd49 100644 --- a/lib/OrdinaryDiffEqExtrapolation/src/extrapolation_perform_step.jl +++ b/lib/OrdinaryDiffEqExtrapolation/src/extrapolation_perform_step.jl @@ -1193,7 +1193,7 @@ function perform_step!(integrator, cache::ImplicitDeuflhardExtrapolationCache, for j in 2:j_int f(k, cache.u_temp1, p, t + (j - 1) * dt_int) OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1) - @.. broadcast=false linsolve_tmps[1]=k - (u_temp1 - u_temp2)/dt_int + @.. broadcast=false linsolve_tmps[1]=k - (u_temp1 - u_temp2) / dt_int linsolve = cache.linsolve[1] @@ -1270,7 +1270,8 @@ function perform_step!(integrator, cache::ImplicitDeuflhardExtrapolationCache, p, t + (j - 1) * dt_int_temp) @.. broadcast=false linsolve_tmps[Threads.threadid()]=k_tmps[Threads.threadid()] - (u_temp3[Threads.threadid()] - - u_temp4[Threads.threadid()])/dt_int_temp + u_temp4[Threads.threadid()]) / + dt_int_temp linsolve = cache.linsolve[Threads.threadid()] @@ -1354,7 +1355,8 @@ function perform_step!(integrator, cache::ImplicitDeuflhardExtrapolationCache, p, t + (j - 1) * dt_int_temp) @.. broadcast=false linsolve_tmps[Threads.threadid()]=k_tmps[Threads.threadid()] - (u_temp3[Threads.threadid()] - - u_temp4[Threads.threadid()])/dt_int_temp + u_temp4[Threads.threadid()]) / + dt_int_temp linsolve = cache.linsolve[Threads.threadid()] @@ -2555,7 +2557,7 @@ function perform_step!(integrator, cache::ImplicitHairerWannerExtrapolationCache for j in 2:(j_int + 1) f(k, cache.u_temp1, p, t + (j - 1) * dt_int) OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1) - @.. broadcast=false linsolve_tmps[1]=k - (u_temp1 - u_temp2)/dt_int + @.. broadcast=false linsolve_tmps[1]=k - (u_temp1 - u_temp2) / dt_int linsolve = cache.linsolve[1] @@ -2634,9 +2636,10 @@ function perform_step!(integrator, cache::ImplicitHairerWannerExtrapolationCache f(k_tmps[Threads.threadid()], cache.u_temp3[Threads.threadid()], p, t + (j - 1) * dt_int_temp) - @.. broadcast=false linsolve_tmps[Threads.threadid()]= k_tmps[Threads.threadid()] - + @.. broadcast=false linsolve_tmps[Threads.threadid()]=k_tmps[Threads.threadid()] - (u_temp3[Threads.threadid()] - - u_temp4[Threads.threadid()]) / dt_int_temp + u_temp4[Threads.threadid()]) / + dt_int_temp linsolve = cache.linsolve[Threads.threadid()] if !repeat_step && j == 1 @@ -2717,7 +2720,8 @@ function perform_step!(integrator, cache::ImplicitHairerWannerExtrapolationCache for j in 2:(j_int_temp + 1) f(ktmp, cache.u_temp3[tid], p, t + (j - 1) * dt_int_temp) @.. broadcast=false linsolvetmp=ktmp - - (u_temp3[tid] - u_temp4[tid])/dt_int_temp + (u_temp3[tid] - u_temp4[tid]) / + dt_int_temp linsolve = cache.linsolve[tid] @@ -2832,7 +2836,7 @@ function perform_step!(integrator, cache::ImplicitHairerWannerExtrapolationCache for j in 2:(j_int + 1) f(k, cache.u_temp1, p, t + (j - 1) * dt_int) OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1) - @.. broadcast=false linsolve_tmps[1]=k - (u_temp1 - u_temp2)/dt_int + @.. broadcast=false linsolve_tmps[1]=k - (u_temp1 - u_temp2) / dt_int linsolve = cache.linsolve[1] diff --git a/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl b/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl index 0432ae1103..6789b307f6 100644 --- a/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl +++ b/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl @@ -902,7 +902,7 @@ end rhs3 = @.. broadcast=false fw3 - β1dt * Mw2-α1dt * Mw3 rhs4 = @.. broadcast=false fw4 - α2dt * Mw4+β2dt * Mw5 rhs5 = @.. broadcast=false fw5 - β2dt * Mw4-α2dt * Mw5 - dw1 = _reshape(LU1 \ _vec(rhs1), axes(u)) + dw1 = _reshape(LU1 \ _vec(rhs1), axes(u)) dw23 = _reshape(LU2 \ _vec(@.. broadcast=false rhs2+rhs3 * im), axes(u)) dw45 = _reshape(LU3 \ _vec(@.. broadcast=false rhs4+rhs5 * im), axes(u)) integrator.stats.nsolve += 3 diff --git a/lib/OrdinaryDiffEqNonlinearSolve/src/initialize_dae.jl b/lib/OrdinaryDiffEqNonlinearSolve/src/initialize_dae.jl index 347621b6de..67c3fb177e 100644 --- a/lib/OrdinaryDiffEqNonlinearSolve/src/initialize_dae.jl +++ b/lib/OrdinaryDiffEqNonlinearSolve/src/initialize_dae.jl @@ -1,21 +1,23 @@ -function default_nlsolve(::Nothing, isinplace::Val{true}, u, ::NonlinearProblem, autodiff = false) +function default_nlsolve( + ::Nothing, isinplace::Val{true}, u, ::NonlinearProblem, autodiff = false) FastShortcutNonlinearPolyalg(; autodiff = autodiff ? AutoForwardDiff() : AutoFiniteDiff()) end function default_nlsolve( - ::Nothing, isinplace::Val{true}, u, ::NonlinearLeastSquaresProblem, autodiff = false) + ::Nothing, isinplace::Val{true}, u, ::NonlinearLeastSquaresProblem, autodiff = false) FastShortcutNLLSPolyalg(; autodiff = autodiff ? AutoForwardDiff() : AutoFiniteDiff()) end -function default_nlsolve(::Nothing, isinplace::Val{false}, u, ::NonlinearProblem, autodiff = false) +function default_nlsolve( + ::Nothing, isinplace::Val{false}, u, ::NonlinearProblem, autodiff = false) FastShortcutNonlinearPolyalg(; autodiff = autodiff ? AutoForwardDiff() : AutoFiniteDiff()) end function default_nlsolve( - ::Nothing, isinplace::Val{false}, u, ::NonlinearLeastSquaresProblem, autodiff = false) + ::Nothing, isinplace::Val{false}, u, ::NonlinearLeastSquaresProblem, autodiff = false) FastShortcutNLLSPolyalg(; autodiff = autodiff ? AutoForwardDiff() : AutoFiniteDiff()) end function default_nlsolve(::Nothing, isinplace::Val{false}, u::StaticArray, - ::NonlinearProblem, autodiff = false) + ::NonlinearProblem, autodiff = false) SimpleTrustRegion(autodiff = autodiff ? AutoForwardDiff() : AutoFiniteDiff()) end function default_nlsolve(::Nothing, isinplace::Val{false}, u::StaticArray, diff --git a/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_caches.jl b/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_caches.jl index 2b44500e74..5d6abd9f42 100644 --- a/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_caches.jl +++ b/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_caches.jl @@ -13,7 +13,8 @@ end # Shampine's Low-order Rosenbrocks mutable struct RosenbrockCache{uType, rateType, uNoUnitsType, JType, WType, TabType, - TFType, UFType, F, JCType, GCType, RTolType, A, StepLimiter, StageLimiter} <: RosenbrockMutableCache + TFType, UFType, F, JCType, GCType, RTolType, A, StepLimiter, StageLimiter} <: + RosenbrockMutableCache u::uType uprev::uType dense::Vector{rateType} @@ -43,7 +44,7 @@ mutable struct RosenbrockCache{uType, rateType, uNoUnitsType, JType, WType, TabT end function full_cache(c::RosenbrockCache) return [c.u, c.uprev, c.dense..., c.du, c.du1, c.du2, - c.ks..., c.fsalfirst, c.fsallast, c.dT, c.tmp, c.atmp, c.weight, c.linsolve_tmp] + c.ks..., c.fsalfirst, c.fsallast, c.dT, c.tmp, c.atmp, c.weight, c.linsolve_tmp] end @cache mutable struct Rosenbrock23Cache{uType, rateType, uNoUnitsType, JType, WType, @@ -716,7 +717,8 @@ tabtype(::Rodas42) = Rodas42Tableau tabtype(::Rodas4P) = Rodas4PTableau tabtype(::Rodas4P2) = Rodas4P2Tableau -function alg_cache(alg::Union{Rodas4, Rodas42, Rodas4P, Rodas4P2}, u, rate_prototype, ::Type{uEltypeNoUnits}, +function alg_cache(alg::Union{Rodas4, Rodas42, Rodas4P, Rodas4P2}, + u, rate_prototype, ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, dt, reltol, p, calck, ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} @@ -731,9 +733,11 @@ function alg_cache(alg::Union{Rodas4, Rodas42, Rodas4P, Rodas4P2}, u, rate_proto alg_autodiff(alg)) end -function alg_cache(alg::Union{Rodas4, Rodas42, Rodas4P, Rodas4P2}, u, rate_prototype, ::Type{uEltypeNoUnits}, - ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, - dt, reltol, p, calck, ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} +function alg_cache(alg::Union{Rodas4, Rodas42, Rodas4P, Rodas4P2}, + u, rate_prototype, ::Type{uEltypeNoUnits}, + ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, + dt, reltol, p, calck, + ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} # Initialize vectors dense = [zero(rate_prototype) for _ in 1:2] @@ -765,11 +769,11 @@ function alg_cache(alg::Union{Rodas4, Rodas42, Rodas4P, Rodas4P2}, u, rate_proto Pl, Pr = wrapprecs( alg.precs(W, nothing, u, p, t, nothing, nothing, nothing, - nothing)..., weight, tmp) + nothing)..., weight, tmp) linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true, - Pl = Pl, Pr = Pr, - assumptions = LinearSolve.OperatorAssumptions(true)) + Pl = Pl, Pr = Pr, + assumptions = LinearSolve.OperatorAssumptions(true)) grad_config = build_grad_config(alg, f, tf, du1, t) jac_config = build_jac_config(alg, f, uf, du1, uprev, u, tmp, du2) diff --git a/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_perform_step.jl b/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_perform_step.jl index 49c6999279..84d93854be 100644 --- a/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_perform_step.jl +++ b/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_perform_step.jl @@ -66,7 +66,7 @@ end @.. veck₁ = vecu * neginvdtγ integrator.stats.nsolve += 1 - @.. u=uprev + dto2 * k₁ + @.. u = uprev + dto2 * k₁ stage_limiter!(u, integrator, p, t + dto2) f(f₁, u, p, t + dto2) OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1) @@ -231,7 +231,7 @@ end if mass_matrix !== I invatol = inv(integrator.opts.abstol) - @.. atmp=ifelse(cache.algebraic_vars, fsallast, false) * invatol + @.. atmp = ifelse(cache.algebraic_vars, fsallast, false) * invatol integrator.EEst += integrator.opts.internalnorm(atmp, t) end end @@ -286,11 +286,11 @@ end if mass_matrix === I linsolve_tmp = @.. (integrator.fsallast - c₃₂ * (k₂ - f₁) - - 2 * (k₁ - integrator.fsalfirst) + dt * dT) + 2 * (k₁ - integrator.fsalfirst) + dt * dT) else linsolve_tmp = mass_matrix * (@.. c₃₂ * k₂ + 2 * k₁) linsolve_tmp = @.. (integrator.fsallast - linsolve_tmp + - c₃₂ * f₁ + 2 * integrator.fsalfirst + dt * dT) + c₃₂ * f₁ + 2 * integrator.fsalfirst + dt * dT) end k₃ = _reshape(W \ _vec(linsolve_tmp), axes(uprev)) * neginvdtγ integrator.stats.nsolve += 1 @@ -306,7 +306,8 @@ end if mass_matrix !== I invatol = inv(integrator.opts.abstol) - atmp = @. ifelse(integrator.differential_vars, false, integrator.fsallast) * invatol + atmp = @. ifelse(integrator.differential_vars, false, integrator.fsallast) * + invatol integrator.EEst += integrator.opts.internalnorm(atmp, t) end end @@ -343,7 +344,7 @@ end return nothing end - k₁ = _reshape(W \ -_vec((integrator.fsalfirst + dtγ * dT)), axes(uprev))/dtγ + k₁ = _reshape(W \ -_vec((integrator.fsalfirst + dtγ * dT)), axes(uprev)) / dtγ integrator.stats.nsolve += 1 tmp = @.. uprev + dto2 * k₁ f₁ = f(tmp, p, t + dto2) @@ -364,11 +365,11 @@ end if mass_matrix === I linsolve_tmp = @.. (integrator.fsallast - c₃₂ * (k₂ - f₁) - - 2(k₁ - integrator.fsalfirst) + dt * dT) + 2(k₁ - integrator.fsalfirst) + dt * dT) else linsolve_tmp = mass_matrix * (@.. c₃₂ * k₂ + 2 * k₁) linsolve_tmp = @.. (integrator.fsallast - linsolve_tmp + - c₃₂ * f₁ + 2 * integrator.fsalfirst + dt * dT) + c₃₂ * f₁ + 2 * integrator.fsalfirst + dt * dT) end k₃ = _reshape(W \ _vec(linsolve_tmp), axes(uprev)) * neginvdtγ integrator.stats.nsolve += 1 @@ -382,7 +383,8 @@ end if mass_matrix !== I invatol = inv(integrator.opts.abstol) - atmp = ifelse(integrator.differential_vars, false, integrator.fsallast) .* invatol + atmp = ifelse(integrator.differential_vars, false, integrator.fsallast) .* + invatol integrator.EEst += integrator.opts.internalnorm(atmp, t) end end @@ -1209,9 +1211,9 @@ function initialize!(integrator, cache::Rodas4ConstantCache) end @muladd function perform_step!(integrator, cache::Rodas4ConstantCache, repeat_step = false) - (;t, dt, uprev, u, f, p) = integrator - (;tf, uf) = cache - (;A, C, gamma, c, d, H) = cache.tab + (; t, dt, uprev, u, f, p) = integrator + (; tf, uf) = cache + (; A, C, gamma, c, d, H) = cache.tab # Precalculations dtC = C ./ dt @@ -1231,7 +1233,7 @@ end end # Initialize ks - num_stages = size(A,1) + num_stages = size(A, 1) du = f(uprev, p, t) linsolve_tmp = @.. du + dtd[1] * dT k1 = _reshape(W \ -_vec(linsolve_tmp), axes(uprev)) @@ -1240,7 +1242,7 @@ end # Loop for stages for stage in 2:num_stages u = uprev - for i in 1:stage-1 + for i in 1:(stage - 1) u = @.. u + A[stage, i] * ks[i] end @@ -1250,11 +1252,11 @@ end # Compute linsolve_tmp for current stage linsolve_tmp = zero(du) if mass_matrix === I - for i in 1:stage-1 + for i in 1:(stage - 1) linsolve_tmp = @.. linsolve_tmp + dtC[stage, i] * ks[i] end else - for i in 1:stage-1 + for i in 1:(stage - 1) linsolve_tmp = @.. linsolve_tmp + dtC[stage, i] * ks[i] end linsolve_tmp = mass_matrix * linsolve_tmp @@ -1297,11 +1299,10 @@ function initialize!(integrator, cache::RosenbrockCache) integrator.k[2] = dense2 end - @muladd function perform_step!(integrator, cache::RosenbrockCache, repeat_step = false) - (;t, dt, uprev, u, f, p) = integrator - (;du, du1, du2, dT, J, W, uf, tf, ks, linsolve_tmp, jac_config, atmp, weight, stage_limiter!, step_limiter!) = cache - (;A, C, gamma, c, d, H) = cache.tab + (; t, dt, uprev, u, f, p) = integrator + (; du, du1, du2, dT, J, W, uf, tf, ks, linsolve_tmp, jac_config, atmp, weight, stage_limiter!, step_limiter!) = cache + (; A, C, gamma, c, d, H) = cache.tab # Assignments sizeu = size(u) @@ -1333,12 +1334,12 @@ end solverdata = (; gamma = dtgamma)) end - @.. $(_vec(ks[1]))=-linres.u + @.. $(_vec(ks[1])) = -linres.u integrator.stats.nsolve += 1 for stage in 2:length(ks) u .= uprev - for i in 1:stage-1 + for i in 1:(stage - 1) @.. u += A[stage, i] * ks[i] end @@ -1348,11 +1349,11 @@ end du1 .= 0 if mass_matrix === I - for i in 1:stage-1 + for i in 1:(stage - 1) @.. du1 += dtC[stage, i] * ks[i] end else - for i in 1:stage-1 + for i in 1:(stage - 1) @.. du1 += dtC[stage, i] * ks[i] end mul!(_vec(du2), mass_matrix, _vec(du1)) @@ -1361,7 +1362,7 @@ end @.. linsolve_tmp = du + dtd[stage] * dT + du1 linres = dolinsolve(integrator, linres.cache; b = _vec(linsolve_tmp)) - @.. $(_vec(ks[stage]))=-linres.u + @.. $(_vec(ks[stage])) = -linres.u integrator.stats.nsolve += 1 end u .+= ks[end] diff --git a/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_tableaus.jl b/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_tableaus.jl index 19e7a989d7..f64d5e1165 100644 --- a/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_tableaus.jl +++ b/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_tableaus.jl @@ -231,28 +231,23 @@ struct RodasTableau{T, T2} H::Matrix{T} end - function Rodas4Tableau(T, T2) gamma = convert(T, 1 // 4) #BET2P=0.0317D0 #BET3P=0.0635D0 #BET4P=0.3438D0 - A = T[ - 0 0 0 0 0 0 - 1.544 0 0 0 0 0 - 0.9466785280815826 0.2557011698983284 0 0 0 0 - 3.314825187068521 2.896124015972201 0.9986419139977817 0 0 0 - 1.221224509226641 6.019134481288629 12.53708332932087 -0.6878860361058950 0 0 - 1.221224509226641 6.019134481288629 12.53708332932087 -0.6878860361058950 1 0 - ] - C = T[ - 0 0 0 0 0 - -5.6688 0 0 0 0 - -2.430093356833875 -0.2063599157091915 0 0 0 - -0.1073529058151375 -9.594562251023355 -20.47028614809616 0 0 - 7.496443313967647 -10.24680431464352 -33.99990352819905 11.70890893206160 0 - 8.083246795921522 -7.981132988064893 -31.52159432874371 16.31930543123136 -6.058818238834054 - ] + A = T[0 0 0 0 0 0 + 1.544 0 0 0 0 0 + 0.9466785280815826 0.2557011698983284 0 0 0 0 + 3.314825187068521 2.896124015972201 0.9986419139977817 0 0 0 + 1.221224509226641 6.019134481288629 12.53708332932087 -0.6878860361058950 0 0 + 1.221224509226641 6.019134481288629 12.53708332932087 -0.6878860361058950 1 0] + C = T[0 0 0 0 0 + -5.6688 0 0 0 0 + -2.430093356833875 -0.2063599157091915 0 0 0 + -0.1073529058151375 -9.594562251023355 -20.47028614809616 0 0 + 7.496443313967647 -10.24680431464352 -33.99990352819905 11.70890893206160 0 + 8.083246795921522 -7.981132988064893 -31.52159432874371 16.31930543123136 -6.058818238834054] c = T2[0, 0.386, 0.21, 0.63, 1, 1] d = T[0.25, -0.1043, 0.1035, -0.0362, 0, 0] H = T[10.12623508344586 -7.487995877610167 -34.80091861555747 -7.992771707568823 1.025137723295662 0 @@ -262,22 +257,18 @@ end function Rodas42Tableau(T, T2) gamma = convert(T, 1 // 4) - A = T[ - 0.0 0 0 0 0 0 - 1.4028884 0 0 0 0 0 - 0.6581212688557198 -1.320936088384301 0 0 0 0 - 7.131197445744498 16.02964143958207 -5.561572550509766 0 0 0 - 22.73885722420363 67.38147284535289 -31.21877493038560 0.7285641833203814 0 0 - 22.73885722420363 67.38147284535289 -31.21877493038560 0.7285641833203814 1 0 - ] - C = T[ - 0 0 0 0 0 - -5.1043536 0 0 0 0 - -2.899967805418783 4.040399359702244 0 0 0 - -32.64449927841361 -99.35311008728094 49.99119122405989 0 0 - -76.46023087151691 -278.5942120829058 153.9294840910643 10.97101866258358 0 - -76.29701586804983 -294.2795630511232 162.0029695867566 23.65166903095270 -7.652977706771382 - ] + A = T[0.0 0 0 0 0 0 + 1.4028884 0 0 0 0 0 + 0.6581212688557198 -1.320936088384301 0 0 0 0 + 7.131197445744498 16.02964143958207 -5.561572550509766 0 0 0 + 22.73885722420363 67.38147284535289 -31.21877493038560 0.7285641833203814 0 0 + 22.73885722420363 67.38147284535289 -31.21877493038560 0.7285641833203814 1 0] + C = T[0 0 0 0 0 + -5.1043536 0 0 0 0 + -2.899967805418783 4.040399359702244 0 0 0 + -32.64449927841361 -99.35311008728094 49.99119122405989 0 0 + -76.46023087151691 -278.5942120829058 153.9294840910643 10.97101866258358 0 + -76.29701586804983 -294.2795630511232 162.0029695867566 23.65166903095270 -7.652977706771382] c = T2[0, 0.3507221, 0.2557041, 0.681779, 1, 1] d = T[0.25, -0.0690221, -0.0009672, -0.087979, 0, 0] H = T[-38.71940424117216 -135.8025833007622 64.51068857505875 -4.192663174613162 -2.531932050335060 0 @@ -290,22 +281,18 @@ function Rodas4PTableau(T, T2) #BET2P=0.D0 #BET3P=c3*c3*(c3/6.d0-GAMMA/2.d0)/(GAMMA*GAMMA) #BET4P=0.3438D0 - A = T[ - 0 0 0 0 0 0 - 3 0 0 0 0 0 - 1.831036793486759 0.4955183967433795 0 0 0 0 - 2.304376582692669 -0.05249275245743001 -1.176798761832782 0 0 0 - -7.170454962423024 -4.741636671481785 -16.31002631330971 -1.062004044111401 0 0 - -7.170454962423024 -4.741636671481785 -16.31002631330971 -1.062004044111401 1 0 - ] - C = T[ - 0 0 0 0 0 - -12 0 0 0 0 - -8.791795173947035 -2.207865586973518 0 0 0 - 10.81793056857153 6.780270611428266 19.53485944642410 0 0 - 34.19095006749676 15.49671153725963 54.74760875964130 14.16005392148534 0 - 34.62605830930532 15.30084976114473 56.99955578662667 18.40807009793095 -5.714285714285717 - ] + A = T[0 0 0 0 0 0 + 3 0 0 0 0 0 + 1.831036793486759 0.4955183967433795 0 0 0 0 + 2.304376582692669 -0.05249275245743001 -1.176798761832782 0 0 0 + -7.170454962423024 -4.741636671481785 -16.31002631330971 -1.062004044111401 0 0 + -7.170454962423024 -4.741636671481785 -16.31002631330971 -1.062004044111401 1 0] + C = T[0 0 0 0 0 + -12 0 0 0 0 + -8.791795173947035 -2.207865586973518 0 0 0 + 10.81793056857153 6.780270611428266 19.53485944642410 0 0 + 34.19095006749676 15.49671153725963 54.74760875964130 14.16005392148534 0 + 34.62605830930532 15.30084976114473 56.99955578662667 18.40807009793095 -5.714285714285717] c = T2[0, 0.75, 0.21, 0.63, 1, 1] d = T[0.25, -0.5, -0.023504, -0.0362, 0, 0] H = T[25.09876703708589 11.62013104361867 28.49148307714626 -5.664021568594133 0 0 @@ -315,26 +302,22 @@ end function Rodas4P2Tableau(T, T2) gamma = convert(T, 1 // 4) - A = T[ - 0 0 0 0 0 0 - 3 0 0 0 0 0 - 0.906377755268814 -0.189707390391685 0 0 0 0 - 3.758617027739064 1.161741776019525 -0.849258085312803 0 0 0 - 7.089566927282776 4.573591406461604 -8.423496976860259 -0.959280113459775 0 0 - 7.089566927282776 4.573591406461604 -8.423496976860259 -0.959280113459775 1 0 - ] - C = T[ - 0 0 0 0 0 - -12 0 0 0 0 - -6.354581592719008 0.338972550544623 0 0 0 - -8.575016317114033 -7.606483992117508 12.224997650124820 0 0 - -5.888975457523102 -8.157396617841821 24.805546872612922 12.790401512796979 0 - -4.408651676063871 -6.692003137674639 24.625568527593117 16.627521966636085 -5.714285714285718 - ] + A = T[0 0 0 0 0 0 + 3 0 0 0 0 0 + 0.906377755268814 -0.189707390391685 0 0 0 0 + 3.758617027739064 1.161741776019525 -0.849258085312803 0 0 0 + 7.089566927282776 4.573591406461604 -8.423496976860259 -0.959280113459775 0 0 + 7.089566927282776 4.573591406461604 -8.423496976860259 -0.959280113459775 1 0] + C = T[0 0 0 0 0 + -12 0 0 0 0 + -6.354581592719008 0.338972550544623 0 0 0 + -8.575016317114033 -7.606483992117508 12.224997650124820 0 0 + -5.888975457523102 -8.157396617841821 24.805546872612922 12.790401512796979 0 + -4.408651676063871 -6.692003137674639 24.625568527593117 16.627521966636085 -5.714285714285718] c = T2[0, 0.75, 0.321448134013046, 0.519745732277726, 1, 1] d = T[0.25, -0.5, -0.189532918363016, 0.085612108792769, 0, 0] - H = [-5.323528268423303 -10.042123754867493 17.175254928256965 -5.079931171878093 -0.016185991706112 0 - 6.984505741529879 6.914061169603662 -0.849178943070653 18.104410789349338 -3.516963011559032 0] + H = [-5.323528268423303 -10.042123754867493 17.175254928256965 -5.079931171878093 -0.016185991706112 0 + 6.984505741529879 6.914061169603662 -0.849178943070653 18.104410789349338 -3.516963011559032 0] RodasTableau(A, C, gamma, c, d, H) end diff --git a/lib/OrdinaryDiffEqRosenbrock/src/stiff_addsteps.jl b/lib/OrdinaryDiffEqRosenbrock/src/stiff_addsteps.jl index f68c812373..35f9359b07 100644 --- a/lib/OrdinaryDiffEqRosenbrock/src/stiff_addsteps.jl +++ b/lib/OrdinaryDiffEqRosenbrock/src/stiff_addsteps.jl @@ -22,9 +22,9 @@ function _ode_addsteps!(k, t, uprev, u, dt, f, p, else J = ForwardDiff.jacobian(uf, uprev) if mass_matrix isa UniformScaling - W = neginvdtγ*mass_matrix + J + W = neginvdtγ * mass_matrix + J else - W = @.. neginvdtγ*mass_matrix .+ J + W = @.. neginvdtγ * mass_matrix .+ J end end f₀ = f(uprev, p, t) @@ -59,7 +59,7 @@ function _ode_addsteps!(k, t, uprev, u, dt, f, p, neginvdtγ = -inv(dtγ) dto2 = dt / 2 - @.. linsolve_tmp=@muladd fsalfirst + dtγ * dT + @.. linsolve_tmp = @muladd fsalfirst + dtγ * dT ### Jacobian does not need to be re-evaluated after an event ### Since it's unchanged @@ -74,7 +74,7 @@ function _ode_addsteps!(k, t, uprev, u, dt, f, p, veck₁ = _vec(k₁) @.. veck₁ = vecu * neginvdtγ - @.. tmp=uprev + dto2 * k₁ + @.. tmp = uprev + dto2 * k₁ f(f₁, tmp, p, t + dto2) if mass_matrix === I @@ -301,8 +301,8 @@ function _ode_addsteps!(k, t, uprev, u, dt, f, p, cache::Rodas4ConstantCache, always_calc_begin = false, allow_calc_end = true, force_calc_end = false) if length(k) < 2 || always_calc_begin - (;tf, uf) = cache - (;A, C, gamma, c, d, H) = cache.tab + (; tf, uf) = cache + (; A, C, gamma, c, d, H) = cache.tab # Precalculations dtC = C ./ dt @@ -327,18 +327,17 @@ function _ode_addsteps!(k, t, uprev, u, dt, f, p, cache::Rodas4ConstantCache, J = ForwardDiff.derivative(uf, uprev) W = 1 / dtgamma - J end - - - num_stages = size(A,1) + + num_stages = size(A, 1) du = f(u, p, t) linsolve_tmp = @.. du + dtd[1] * dT k1 = _reshape(W \ _vec(linsolve_tmp), axes(uprev)) # constant number for type stability make sure this is greater than num_stages ks = ntuple(Returns(k1), 10) # Last stage doesn't affect ks - for stage in 2:num_stages-1 + for stage in 2:(num_stages - 1) u = uprev - for i in 1:stage-1 + for i in 1:(stage - 1) u = @.. u + A[stage, i] * ks[i] end @@ -347,11 +346,11 @@ function _ode_addsteps!(k, t, uprev, u, dt, f, p, cache::Rodas4ConstantCache, # Compute linsolve_tmp for current stage linsolve_tmp = zero(du) if mass_matrix === I - for i in 1:stage-1 + for i in 1:(stage - 1) linsolve_tmp = @.. linsolve_tmp + dtC[stage, i] * ks[i] end else - for i in 1:stage-1 + for i in 1:(stage - 1) linsolve_tmp = @.. linsolve_tmp + dtC[stage, i] * ks[i] end linsolve_tmp = mass_matrix * linsolve_tmp @@ -364,7 +363,7 @@ function _ode_addsteps!(k, t, uprev, u, dt, f, p, cache::Rodas4ConstantCache, k2 = zero(ks[1]) H = cache.tab.H # Last stage doesn't affect ks - for i in 1:num_stages-1 + for i in 1:(num_stages - 1) k1 = @.. k1 + H[1, i] * ks[i] k2 = @.. k2 + H[2, i] * ks[i] end @@ -379,8 +378,8 @@ function _ode_addsteps!(k, t, uprev, u, dt, f, p, cache::RosenbrockCache, always_calc_begin = false, allow_calc_end = true, force_calc_end = false) if length(k) < 2 || always_calc_begin - (;du, du1, du2, tmp, ks, dT, J, W, uf, tf, linsolve_tmp, jac_config, fsalfirst, weight) = cache - (;A, C, gamma, c, d , H) = cache.tab + (; du, du1, du2, tmp, ks, dT, J, W, uf, tf, linsolve_tmp, jac_config, fsalfirst, weight) = cache + (; A, C, gamma, c, d, H) = cache.tab # Assignments sizeu = size(u) @@ -392,45 +391,47 @@ function _ode_addsteps!(k, t, uprev, u, dt, f, p, cache::RosenbrockCache, dtd = dt .* d dtgamma = dt * gamma - @.. linsolve_tmp=@muladd fsalfirst + dtgamma * dT + @.. linsolve_tmp = @muladd fsalfirst + dtgamma * dT # Jacobian does not need to be re-evaluated after an event since it's unchanged jacobian2W!(W, mass_matrix, dtgamma, J, true) linsolve = cache.linsolve - linres = dolinsolve(cache, linsolve; A = W, b = _vec(linsolve_tmp), reltol = cache.reltol) - @.. $(_vec(ks[1]))=-linres.u + linres = dolinsolve( + cache, linsolve; A = W, b = _vec(linsolve_tmp), reltol = cache.reltol) + @.. $(_vec(ks[1])) = -linres.u # Last stage doesn't affect ks - for stage in 2:length(ks)-1 + for stage in 2:(length(ks) - 1) tmp .= uprev - for i in 1:stage-1 + for i in 1:(stage - 1) @.. tmp += A[stage, i] * _vec(ks[i]) end f(du, tmp, p, t + c[stage] * dt) if mass_matrix === I @.. linsolve_tmp = du + dtd[stage] * dT - for i in 1:stage-1 + for i in 1:(stage - 1) @.. linsolve_tmp += dtC[stage, i] * _vec(ks[i]) end else du1 .= du - for i in 1:stage-1 + for i in 1:(stage - 1) @.. du1 += dtC[stage, i] * _vec(ks[i]) end mul!(_vec(du2), mass_matrix, _vec(du1)) @.. linsolve_tmp = du + dtd[stage] * dT + du2 end - linres = dolinsolve(cache, linres.cache; b = _vec(linsolve_tmp), reltol = cache.reltol) - @.. $(_vec(ks[stage]))=-linres.u + linres = dolinsolve( + cache, linres.cache; b = _vec(linsolve_tmp), reltol = cache.reltol) + @.. $(_vec(ks[stage])) = -linres.u end copyat_or_push!(k, 1, zero(du)) copyat_or_push!(k, 2, zero(du)) # Last stage doesn't affect ks - for i in 1:length(ks)-1 + for i in 1:(length(ks) - 1) @.. k[1] += H[1, i] * _vec(ks[i]) @.. k[2] += H[2, i] * _vec(ks[i]) end diff --git a/test/integrators/callback_allocation_tests.jl b/test/integrators/callback_allocation_tests.jl index 6b9bc02e99..a6d3f47b26 100644 --- a/test/integrators/callback_allocation_tests.jl +++ b/test/integrators/callback_allocation_tests.jl @@ -32,8 +32,8 @@ cbs = CallbackSet(ContinuousCallback(cond_1, cb_affect!), ContinuousCallback(cond_9, cb_affect!)) integrator = init( - ODEProblem{true, SciMLBase.FullSpecialize}(f!, [0.8, 1.0], - (0.0, 100.0), [0, 0]), Tsit5(), callback = cbs, + ODEProblem{true, SciMLBase.FullSpecialize}(f!, [0.8, 1.0], + (0.0, 100.0), [0, 0]), Tsit5(), callback = cbs, save_on = false); # Force a callback event to occur so we can call handle_callbacks! directly. # Step to a point where u[1] is still > 0.5, so we can force it below 0.5 and diff --git a/test/interface/checkinit_tests.jl b/test/interface/checkinit_tests.jl index e50a55a9b1..8320449151 100644 --- a/test/interface/checkinit_tests.jl +++ b/test/interface/checkinit_tests.jl @@ -24,8 +24,11 @@ roberf_oop = ODEFunction{false}(rober, mass_matrix = M) prob_mm = ODEProblem(roberf, [1.0, 0.0, 0.2], (0.0, 1e5), (0.04, 3e7, 1e4)) prob_mm_oop = ODEProblem(roberf_oop, [1.0, 0.0, 0.2], (0.0, 1e5), (0.04, 3e7, 1e4)) -@test_throws OrdinaryDiffEqCore.CheckInitFailureError solve(prob_mm, Rodas5P(), reltol = 1e-8, abstol = 1e-8, initializealg = SciMLBase.CheckInit()) -@test_throws OrdinaryDiffEqCore.CheckInitFailureError solve(prob_mm_oop, Rodas5P(), reltol = 1e-8, abstol = 1e-8, initializealg = SciMLBase.CheckInit()) +@test_throws OrdinaryDiffEqCore.CheckInitFailureError solve( + prob_mm, Rodas5P(), reltol = 1e-8, abstol = 1e-8, initializealg = SciMLBase.CheckInit()) +@test_throws OrdinaryDiffEqCore.CheckInitFailureError solve( + prob_mm_oop, Rodas5P(), reltol = 1e-8, abstol = 1e-8, + initializealg = SciMLBase.CheckInit()) f_oop = function (du, u, p, t) out1 = -0.04u[1] + 1e4 * u[2] * u[3] - du[1] @@ -46,5 +49,7 @@ tspan = (0.0, 100000.0) differential_vars = [true, true, false] prob = DAEProblem(f, du₀, u₀, tspan, differential_vars = differential_vars) prob_oop = DAEProblem(f_oop, du₀, u₀, tspan, differential_vars = differential_vars) -@test_throws OrdinaryDiffEqCore.CheckInitFailureError solve(prob, DFBDF(), reltol = 1e-8, abstol = 1e-8, initializealg = SciMLBase.CheckInit()) -@test_throws OrdinaryDiffEqCore.CheckInitFailureError solve(prob_oop, DFBDF(), reltol = 1e-8, abstol = 1e-8, initializealg = SciMLBase.CheckInit()) \ No newline at end of file +@test_throws OrdinaryDiffEqCore.CheckInitFailureError solve( + prob, DFBDF(), reltol = 1e-8, abstol = 1e-8, initializealg = SciMLBase.CheckInit()) +@test_throws OrdinaryDiffEqCore.CheckInitFailureError solve( + prob_oop, DFBDF(), reltol = 1e-8, abstol = 1e-8, initializealg = SciMLBase.CheckInit()) diff --git a/test/interface/composite_algorithm_test.jl b/test/interface/composite_algorithm_test.jl index f3e091e3d4..be3a496fcf 100644 --- a/test/interface/composite_algorithm_test.jl +++ b/test/interface/composite_algorithm_test.jl @@ -95,5 +95,6 @@ M = [1.0 0 0 0 0 0] f = ODEFunction(rober, mass_matrix = M) prob_mm = ODEProblem(f, [1.0, 0.0, 0.0], (0.0, 1e5), (0.04, 3e7, 1e4)) -cb = DiscreteCallback((u,t,integrator)->true, (integrator)->u_modified!(integrator,true)) -sol = solve(prob_mm, DefaultODEAlgorithm(), callback = cb) \ No newline at end of file +cb = DiscreteCallback( + (u, t, integrator) -> true, (integrator) -> u_modified!(integrator, true)) +sol = solve(prob_mm, DefaultODEAlgorithm(), callback = cb)