From e6a9f6b7804715d6a93eda9b5eeec75d23ecbf1b Mon Sep 17 00:00:00 2001 From: deeepeshthakur Date: Mon, 21 Jan 2019 01:21:46 +0530 Subject: [PATCH] Implemented NDBLSRK144 scheme --- src/OrdinaryDiffEq.jl | 2 +- src/alg_utils.jl | 1 + src/algorithms.jl | 1 + src/caches/low_order_rk_caches.jl | 108 ++++++++++++++ src/perform_step/low_order_rk_perform_step.jl | 139 ++++++++++++++++++ test/ode/ode_ssprk_tests.jl | 18 +++ 6 files changed, 268 insertions(+), 1 deletion(-) diff --git a/src/OrdinaryDiffEq.jl b/src/OrdinaryDiffEq.jl index fafdbd3867..044e9c8adb 100644 --- a/src/OrdinaryDiffEq.jl +++ b/src/OrdinaryDiffEq.jl @@ -154,7 +154,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, CFRLDDRK64, 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 3c0f55a878..69896ded63 100644 --- a/src/alg_utils.jl +++ b/src/alg_utils.jl @@ -186,6 +186,7 @@ alg_order(alg::OwrenZen4) = 4 alg_order(alg::OwrenZen5) = 5 alg_order(alg::LDDRK64) = 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 fbfc460f4a..a4bc9233ec 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -47,6 +47,7 @@ struct OwrenZen4 <: OrdinaryDiffEqAdaptiveAlgorithm end struct OwrenZen5 <: OrdinaryDiffEqAdaptiveAlgorithm end struct LDDRK64 <: 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 7295d8c8e9..c7d7f60be6 100644 --- a/src/caches/low_order_rk_caches.jl +++ b/src/caches/low_order_rk_caches.jl @@ -595,3 +595,111 @@ end 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 fa1ae9b92f..f0e414c753 100644 --- a/src/perform_step/low_order_rk_perform_step.jl +++ b/src/perform_step/low_order_rk_perform_step.jl @@ -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 17281bee48..62caa3eabf 100644 --- a/test/ode/ode_ssprk_tests.jl +++ b/test/ode/ode_ssprk_tests.jl @@ -312,6 +312,8 @@ for prob in test_problems_nonlinear @test_broken sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol end +# 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) @@ -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