diff --git a/REQUIRE b/REQUIRE index e8b79cb768..7d5609dbf5 100644 --- a/REQUIRE +++ b/REQUIRE @@ -7,7 +7,6 @@ GenericSVD 0.0.2 NLsolve 0.14.1 RecursiveArrayTools 0.13.0 DiffEqDiffTools 0.4.0 -DataStructures 0.4.6 Reexport MuladdMacro 0.2.1 StaticArrays diff --git a/src/OrdinaryDiffEq.jl b/src/OrdinaryDiffEq.jl index 0872b1ce26..1c81c1db75 100644 --- a/src/OrdinaryDiffEq.jl +++ b/src/OrdinaryDiffEq.jl @@ -157,7 +157,7 @@ module OrdinaryDiffEq export FunctionMap, Euler, Heun, Ralston, Midpoint, SSPRK22, SSPRK33, SSPRK53, SSPRK53_2N1,SSPRK53_2N2, SSPRK63, SSPRK73, SSPRK83, SSPRK432, SSPRK932, SSPRK54, SSPRK104, RK4, ExplicitRK, OwrenZen3, OwrenZen4, OwrenZen5, - LDDRK64, LDDRK46, BS3, BS5, CarpenterKennedy2N54, + LDDRK64, CFRLDDRK64, NDBLSRK144, BS3, BS5, CarpenterKennedy2N54, ORK256, RK46NL, DP5, DP5Threaded, Tsit5, DP8, Vern6, Vern7, Vern8, TanYam7, TsitPap8, Vern9,Feagin10, Feagin12, Feagin14, CompositeAlgorithm, Anas5 diff --git a/src/alg_utils.jl b/src/alg_utils.jl index 6848b6e3dd..9a1f46fca9 100644 --- a/src/alg_utils.jl +++ b/src/alg_utils.jl @@ -186,7 +186,8 @@ alg_order(alg::OwrenZen3) = 3 alg_order(alg::OwrenZen4) = 4 alg_order(alg::OwrenZen5) = 5 alg_order(alg::LDDRK64) = 4 -alg_order(alg::LDDRK46) = 4 +alg_order(alg::CFRLDDRK64) = 4 +alg_order(alg::NDBLSRK144) = 4 alg_order(alg::DP5) = 5 alg_order(alg::DP5Threaded) = 5 diff --git a/src/algorithms.jl b/src/algorithms.jl index 94d37bc799..ba7b0c2336 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -46,7 +46,8 @@ struct OwrenZen3 <: OrdinaryDiffEqAdaptiveAlgorithm end struct OwrenZen4 <: OrdinaryDiffEqAdaptiveAlgorithm end struct OwrenZen5 <: OrdinaryDiffEqAdaptiveAlgorithm end struct LDDRK64 <: OrdinaryDiffEqAlgorithm end -struct LDDRK46 <: OrdinaryDiffEqAlgorithm end +struct CFRLDDRK64 <: OrdinaryDiffEqAlgorithm end +struct NDBLSRK144 <: OrdinaryDiffEqAlgorithm end struct CarpenterKennedy2N54 <: OrdinaryDiffEqAlgorithm end struct ORK256 <: OrdinaryDiffEqAlgorithm end struct SSPRK22{StageLimiter,StepLimiter} <: OrdinaryDiffEqAlgorithm diff --git a/src/caches/low_order_rk_caches.jl b/src/caches/low_order_rk_caches.jl index aae5c93c59..c7d7f60be6 100644 --- a/src/caches/low_order_rk_caches.jl +++ b/src/caches/low_order_rk_caches.jl @@ -536,7 +536,7 @@ end alg_cache(alg::Anas5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = Anas5ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits)) -@cache struct LDDRK46Cache{uType,rateType,TabType} <: OrdinaryDiffEqMutableCache +@cache struct CFRLDDRK64Cache{uType,rateType,TabType} <: OrdinaryDiffEqMutableCache u::uType uprev::uType k::rateType @@ -545,7 +545,7 @@ alg_cache(alg::Anas5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeN tab::TabType end -struct LDDRK46ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache +struct CFRLDDRK64ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache α1::T α2::T α3::T @@ -563,7 +563,7 @@ struct LDDRK46ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache c5::T2 c6::T2 - function LDDRK46ConstantCache(::Type{T}, ::Type{T2}) where {T,T2} + function CFRLDDRK64ConstantCache(::Type{T}, ::Type{T2}) where {T,T2} α1 = T(0.17985400977138) α2 = T(0.14081893152111) α3 = T(0.08255631629428) @@ -584,14 +584,122 @@ struct LDDRK46ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache end end -function alg_cache(alg::LDDRK46,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}}) +function alg_cache(alg::CFRLDDRK64,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}}) tmp = similar(u) k = zero(rate_prototype) fsalfirst = zero(rate_prototype) - tab = LDDRK46ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits)) - LDDRK46Cache(u,uprev,k,tmp,fsalfirst,tab) + tab = CFRLDDRK64ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits)) + CFRLDDRK64Cache(u,uprev,k,tmp,fsalfirst,tab) end -function alg_cache(alg::LDDRK46,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) - LDDRK46ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits)) +function alg_cache(alg::CFRLDDRK64,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) + CFRLDDRK64ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits)) +end + +@cache struct NDBLSRK144Cache{uType,rateType,TabType} <: OrdinaryDiffEqMutableCache + u::uType + uprev::uType + k::rateType + tmp::uType + fsalfirst::rateType + tab::TabType +end + +struct NDBLSRK144ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache + α2::T + α3::T + α4::T + α5::T + α6::T + α7::T + α8::T + α9::T + α10::T + α11::T + α12::T + α13::T + α14::T + β1::T + β2::T + β3::T + β4::T + β5::T + β6::T + β7::T + β8::T + β9::T + β10::T + β11::T + β12::T + β13::T + β14::T + c2::T2 + c3::T2 + c4::T2 + c5::T2 + c6::T2 + c7::T2 + c8::T2 + c9::T2 + c10::T2 + c11::T2 + c12::T2 + c13::T2 + c14::T2 + + function NDBLSRK144ConstantCache(::Type{T}, ::Type{T2}) where {T,T2} + α2 = T(-0.7188012108672410) + α3 = T(-0.7785331173421570) + α4 = T(-0.0053282796654044) + α5 = T(-0.8552979934029281) + α6 = T(-3.9564138245774565) + α7 = T(-1.5780575380587385) + α8 = T(-2.0837094552574054) + α9 = T(-0.7483334182761610) + α10 = T(-0.7032861106563359) + α11 = T(0.0013917096117681) + α12 = T(-0.0932075369637460) + α13 = T(-0.9514200470875948) + α14 = T(-7.1151571693922548) + β1 = T(0.0367762454319673) + β2 = T(0.3136296607553959) + β3 = T(0.1531848691869027) + β4 = T(0.0030097086818182) + β5 = T(0.3326293790646110) + β6 = T(0.2440251405350864) + β7 = T(0.3718879239592277) + β8 = T(0.6204126221582444) + β9 = T(0.1524043173028741) + β10 = T(0.0760894927419266) + β11 = T(0.0077604214040978) + β12 = T(0.0024647284755382) + β13 = T(0.0780348340049386) + β14 = T(5.5059777270269628) + c2 = T2(0.0367762454319673) + c3 = T2(0.1249685262725025) + c4 = T2(0.2446177702277698) + c5 = T2(0.2476149531070420) + c6 = T2(0.2969311120382472) + c7 = T2(0.3978149645802642) + c8 = T2(0.5270854589440328) + c9 = T2(0.6981269994175695) + c10 = T2(0.8190890835352128) + c11 = T2(0.8527059887098624) + c12 = T2(0.8604711817462826) + c13 = T2(0.8627060376969976) + c14 = T2(0.8734213127600976) + new{T,T2}(α2, α3, α4, α5, α6, α7, α8, α9, α10, α11, α12, α13, α14, β1, β2, β3, β4, β5, β6, β7, β8, β9, β10, β11, β12, β13, β14, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11, c12, c13, c14) + end +end + +function alg_cache(alg::NDBLSRK144,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}}) + tmp = similar(u) + k = zero(rate_prototype) + fsalfirst = zero(rate_prototype) + tab = NDBLSRK144ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits)) + NDBLSRK144Cache(u,uprev,k,tmp,fsalfirst,tab) +end + +function alg_cache(alg::NDBLSRK144,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) + NDBLSRK144ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits)) end diff --git a/src/perform_step/low_order_rk_perform_step.jl b/src/perform_step/low_order_rk_perform_step.jl index 7080a0df7b..f0e414c753 100644 --- a/src/perform_step/low_order_rk_perform_step.jl +++ b/src/perform_step/low_order_rk_perform_step.jl @@ -804,7 +804,7 @@ end end end -function initialize!(integrator,cache::LDDRK46ConstantCache) +function initialize!(integrator,cache::CFRLDDRK64ConstantCache) integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal integrator.kshortsize = 1 integrator.k = typeof(integrator.k)(undef, integrator.kshortsize) @@ -814,7 +814,7 @@ function initialize!(integrator,cache::LDDRK46ConstantCache) integrator.k[1] = integrator.fsalfirst end -@muladd function perform_step!(integrator,cache::LDDRK46ConstantCache,repeat_step=false) +@muladd function perform_step!(integrator,cache::CFRLDDRK64ConstantCache,repeat_step=false) @unpack t,dt,uprev,u,f,p = integrator @unpack α1,α2,α3,α4,α5,β1,β2,β3,β4,β5,β6,c2,c3,c4,c5,c6 = cache @@ -844,7 +844,7 @@ end integrator.u = u end -function initialize!(integrator,cache::LDDRK46Cache) +function initialize!(integrator,cache::CFRLDDRK64Cache) @unpack k,fsalfirst = cache integrator.fsalfirst = fsalfirst integrator.fsallast = k @@ -854,7 +854,7 @@ function initialize!(integrator,cache::LDDRK46Cache) integrator.f(integrator.fsalfirst,integrator.uprev,integrator.p,integrator.t) # FSAL for interpolation end -@muladd function perform_step!(integrator,cache::LDDRK46Cache,repeat_step=false) +@muladd function perform_step!(integrator,cache::CFRLDDRK64Cache,repeat_step=false) @unpack t,dt,uprev,u,f,p = integrator @unpack k,fsalfirst,tmp = cache @unpack α1,α2,α3,α4,α5,β1,β2,β3,β4,β5,β6,c2,c3,c4,c5,c6 = cache.tab @@ -887,3 +887,142 @@ end f( k, u, p, t+dt) end + +function initialize!(integrator,cache::NDBLSRK144ConstantCache) + integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal + integrator.kshortsize = 1 + integrator.k = typeof(integrator.k)(undef, integrator.kshortsize) + + # Avoid undefined entries if k is an array of arrays + integrator.fsallast = zero(integrator.fsalfirst) + integrator.k[1] = integrator.fsalfirst +end + +@muladd function perform_step!(integrator,cache::NDBLSRK144ConstantCache,repeat_step=false) + @unpack t,dt,uprev,u,f,p = integrator + @unpack α2, α3, α4, α5, α6, α7, α8, α9, α10, α11, α12, α13, α14, β1, β2, β3, β4, β5, β6, β7, β8, β9, β10, β11, β12, β13, β14, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11, c12, c13, c14 = cache + + # u0 + u = uprev + #u1 + tmp = dt*integrator.fsalfirst + u = u + β1*tmp + #u2 + tmp = α2*tmp + dt*f(u, p, t+c2*dt) + u = u + β2*tmp + #u3 + tmp = α3*tmp + dt*f(u, p, t+c3*dt) + u = u + β3*tmp + #u4 + tmp = α4*tmp + dt*f(u, p, t+c4*dt) + u = u + β4*tmp + #u5 + tmp = α5*tmp + dt*f(u, p, t+c5*dt) + u = u + β5*tmp + #u6 + tmp = α6*tmp + dt*f(u, p, t+c6*dt) + u = u + β6*tmp + #u7 + tmp = α7*tmp + dt*f(u, p, t+c7*dt) + u = u + β7*tmp + #u8 + tmp = α8*tmp + dt*f(u, p, t+c8*dt) + u = u + β8*tmp + #u9 + tmp = α9*tmp + dt*f(u, p, t+c9*dt) + u = u + β9*tmp + #u10 + tmp = α10*tmp + dt*f(u, p, t+c10*dt) + u = u + β10*tmp + #u11 + tmp = α11*tmp + dt*f(u, p, t+c11*dt) + u = u + β11*tmp + #u12 + tmp = α12*tmp + dt*f(u, p, t+c12*dt) + u = u + β12*tmp + #u13 + tmp = α13*tmp + dt*f(u, p, t+c13*dt) + u = u + β13*tmp + #u14 + tmp = α14*tmp + dt*f(u, p, t+c14*dt) + u = u + β14*tmp + integrator.fsallast = f(u, p, t+dt) # For interpolation, then FSAL'd + integrator.k[1] = integrator.fsalfirst + integrator.u = u +end + +function initialize!(integrator,cache::NDBLSRK144Cache) + @unpack k,fsalfirst = cache + integrator.fsalfirst = fsalfirst + integrator.fsallast = k + integrator.kshortsize = 1 + resize!(integrator.k, integrator.kshortsize) + integrator.k[1] = integrator.fsalfirst + integrator.f(integrator.fsalfirst,integrator.uprev,integrator.p,integrator.t) # FSAL for interpolation +end + +@muladd function perform_step!(integrator,cache::NDBLSRK144Cache,repeat_step=false) + @unpack t,dt,uprev,u,f,p = integrator + @unpack k,fsalfirst,tmp = cache + @unpack α2, α3, α4, α5, α6, α7, α8, α9, α10, α11, α12, α13, α14, β1, β2, β3, β4, β5, β6, β7, β8, β9, β10, β11, β12, β13, β14, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11, c12, c13, c14 = cache.tab + + #u0 + @. u = uprev + #u1 + @. tmp = dt*fsalfirst + @. u = u + β1*tmp + #u2 + f(k, u, p, t+c2*dt) + @. tmp = α2*tmp + dt*k + @. u = u + β2*tmp + #u3 + f(k, u, p, t+c3*dt) + @. tmp = α3*tmp + dt*k + @. u = u + β3*tmp + #u4 + f(k, u, p, t+c4*dt) + @. tmp = α4*tmp + dt*k + @. u = u + β4*tmp + #u5 + f(k, u, p, t+c5*dt) + @. tmp = α5*tmp + dt*k + @. u = u + β5*tmp + #u6 + f(k, u, p, t+c6*dt) + @. tmp = α6*tmp + dt*k + @. u = u + β6*tmp + #u7 + f(k, u, p, t+c7*dt) + @. tmp = α7*tmp + dt*k + @. u = u + β7*tmp + #u8 + f(k, u, p, t+c8*dt) + @. tmp = α8*tmp + dt*k + @. u = u + β8*tmp + #u9 + f(k, u, p, t+c9*dt) + @. tmp = α9*tmp + dt*k + @. u = u + β9*tmp + #u10 + f(k, u, p, t+c10*dt) + @. tmp = α10*tmp + dt*k + @. u = u + β10*tmp + #u11 + f(k, u, p, t+c11*dt) + @. tmp = α11*tmp + dt*k + @. u = u + β11*tmp + #u12 + f(k, u, p, t+c12*dt) + @. tmp = α12*tmp + dt*k + @. u = u + β12*tmp + #u13 + f(k, u, p, t+c13*dt) + @. tmp = α13*tmp + dt*k + @. u = u + β13*tmp + #u14 + f(k, u, p, t+c14*dt) + @. tmp = α14*tmp + dt*k + @. u = u + β14*tmp + + f( k, u, p, t+dt) +end diff --git a/test/ode/ode_ssprk_tests.jl b/test/ode/ode_ssprk_tests.jl index 13a8fb0951..62caa3eabf 100644 --- a/test/ode/ode_ssprk_tests.jl +++ b/test/ode/ode_ssprk_tests.jl @@ -312,7 +312,9 @@ for prob in test_problems_nonlinear @test_broken sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol end -alg = LDDRK46() +# for CFRLDDRK64 to be in asymptotic range +dts = 1 .//2 .^(7:-1:4) +alg = CFRLDDRK64() for prob in test_problems_only_time sim = test_convergence(dts, prob, alg) @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol @@ -325,3 +327,19 @@ for prob in test_problems_nonlinear sim = test_convergence(dts, prob, alg) @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol end +#reverting back to original dts +dts = 1 .//2 .^(8:-1:4) + +alg = NDBLSRK144() +for prob in test_problems_only_time + sim = test_convergence(dts, prob, alg) + @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol +end +for prob in test_problems_linear + sim = test_convergence(dts, prob, alg) + @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol +end +for prob in test_problems_nonlinear + sim = test_convergence(dts, prob, alg) + @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol +end