From 2f2d4e51c543501f997b7e015d780de44344478a Mon Sep 17 00:00:00 2001 From: David Widmann Date: Sun, 10 Sep 2017 14:40:58 +0200 Subject: [PATCH 1/5] Fix Vern6 interpolant --- src/dense/interpolants.jl | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/dense/interpolants.jl b/src/dense/interpolants.jl index 246a61d94a..447081be03 100644 --- a/src/dense/interpolants.jl +++ b/src/dense/interpolants.jl @@ -576,7 +576,13 @@ end b12Θ = @evalpoly(Θ, 0, 0, r122, r123, r124, r125, r126) if out == nothing - return y₀[idxs] + dt*(k[1][idxs]*b1Θ + k[4][idxs]*b4Θ + k[5][idxs]*b5Θ + k[6][idxs]*b6Θ + k[7][idxs]*b7Θ + k[8][idxs]*b8Θ + k[9][idxs]*b9Θ + k[10][idxs]*b10Θ + k[11][idxs]*b11Θ + k[12][idxs]*b12Θ) + if idxs == nothing + # return @. y₀ + dt*(k[1]*b1Θ + k[4]*b4Θ + k[5]*b5Θ + k[6]*b6Θ + k[7]*b7Θ + k[8]*b8Θ + k[9]*b9Θ + k[10]*b10Θ + k[11]*b11Θ + k[12]*b12Θ) + return y₀ + dt*(k[1]*b1Θ + k[4]*b4Θ + k[5]*b5Θ + k[6]*b6Θ + k[7]*b7Θ + k[8]*b8Θ + k[9]*b9Θ + k[10]*b10Θ + k[11]*b11Θ + k[12]*b12Θ) + else + # return @. y₀[idxs] + dt*(k[1][idxs]*b1Θ + k[4][idxs]*b4Θ + k[5][idxs]*b5Θ + k[6][idxs]*b6Θ + k[7][idxs]*b7Θ + k[8][idxs]*b8Θ + k[9][idxs]*b9Θ + k[10][idxs]*b10Θ + k[11][idxs]*b11Θ + k[12][idxs]*b12Θ) + return y₀[idxs] + dt*(k[1][idxs]*b1Θ + k[4][idxs]*b4Θ + k[5][idxs]*b5Θ + k[6][idxs]*b6Θ + k[7][idxs]*b7Θ + k[8][idxs]*b8Θ + k[9][idxs]*b9Θ + k[10][idxs]*b10Θ + k[11][idxs]*b11Θ + k[12][idxs]*b12Θ) + end elseif idxs == nothing #@. out = y₀ + dt*(k[1]*b1Θ + k[4]*b4Θ + k[5]*b5Θ + k[6]*b6Θ + k[7]*b7Θ + k[8]*b8Θ + k[9]*b9Θ + k[10]*b10Θ + k[11]*b11Θ + k[12]*b12Θ) @inbounds for i in eachindex(out) From 2480a9b0bea108ea9abfa682a975a40227f7c3f1 Mon Sep 17 00:00:00 2001 From: YingboMa Date: Tue, 12 Sep 2017 13:04:18 -0400 Subject: [PATCH 2/5] Add ERKN5 tableau --- src/tableaus/rkn_tableaus.jl | 82 ++++++++++++++++++++++++++++++++++++ 1 file changed, 82 insertions(+) diff --git a/src/tableaus/rkn_tableaus.jl b/src/tableaus/rkn_tableaus.jl index 57a9caa6f4..b346936d7b 100644 --- a/src/tableaus/rkn_tableaus.jl +++ b/src/tableaus/rkn_tableaus.jl @@ -227,6 +227,88 @@ ERKN4ConstantCache( T(0.0016835016835016834)) end +struct ERKN5ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache + c1::T2 + c2::T2 + c3::T2 + a21::T + a31::T + a32::T + a41::T + a42::T + a43::T + b1::T + b2::T + b3::T + b4::T + bp1::T # bp denotes bprime + bp2::T + bp3::T + bp4::T + btilde1::T + btilde2::T + btilde3::T + btilde4::T + # bptilde1::T + # bptilde2::T + # bptilde3::T + # bptilde4::T +end + +function ERKN5ConstantCache(T::Type,T2::Type) + c1 = T2(1//2) + c2 = T2(19//70) + c3 = T2(44//51) + a21 = T(1//8) + a31 = T(2907//343000) + a32 = T(1216//42875) + a41 = T(6624772//128538819) + a42 = T(6273905//54121608) + a43 = T(210498365//1028310552) + b1 = T(479//5016) + b2 = T(235//1776) + b3 = T(145775//641744) + b4 = T(309519//6873416) + btilde1 = T(479//5016 - 184883//2021250) + btilde2 = T(235//1776 - 411163//3399375) + btilde3 = T(145775//641744 - 6//25) + btilde4 = T(309519//6873416 - 593028//12464375) + bp1 = b1 + bp2 = T(235//888) + bp3 = T(300125//962616) + bp4 = T(2255067//6873416) + # bptilde1 = T(0) + # bptilde2 = T(0) + # bptilde3 = T(0) + # bptilde4 = T(0) + ERKN5ConstantCache(c1, c2, c3, a21, a31, a32, a41, a42, a43, b1, b2, b3, b4, bp1, bp2, bp3, bp4, btilde1, btilde2, btilde3, btilde4) +end + +Base.@pure function ERKN5ConstantCache{T<:CompiledFloats,T2<:CompiledFloats}(::Type{T},::Type{T2}) + ERKN5ConstantCache( + T2(0.5), + T2(0.2714285714285714), + T2(0.8627450980392157), + T(0.125), + T(0.008475218658892128), + T(0.028361516034985424), + T(0.051539076300366506), + T(0.11592236875149756), + T(0.20470310704348388), + T(0.09549441786283891), + T(0.13231981981981983), + T(0.22715444164651324), + T(0.04503132067082801), + T(0.09549441786283891), + T(0.26463963963963966), + T(0.3117806061814888), + T(0.32808533631603265), + T(0.004024782736060931), + T(0.011367291781577495), + T(-0.012845558353486749), + T(-0.0025465161641516788)) +end + struct DPRKN6ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache c1::T2 c2::T2 From aca4f762fd8b76dd0a160701fca9d7b86170afcb Mon Sep 17 00:00:00 2001 From: YingboMa Date: Tue, 12 Sep 2017 13:05:09 -0400 Subject: [PATCH 3/5] Add ERKN5 utilities --- src/OrdinaryDiffEq.jl | 2 +- src/alg_utils.jl | 2 ++ src/algorithms.jl | 1 + src/caches/rkn_caches.jl | 31 +++++++++++++++++++++++++++++++ 4 files changed, 35 insertions(+), 1 deletion(-) diff --git a/src/OrdinaryDiffEq.jl b/src/OrdinaryDiffEq.jl index b1a9f54d85..32c6dcd998 100644 --- a/src/OrdinaryDiffEq.jl +++ b/src/OrdinaryDiffEq.jl @@ -152,5 +152,5 @@ module OrdinaryDiffEq export SplitEuler export Nystrom4, Nystrom4VelocityIndependent, Nystrom5VelocityIndependent, - IRKN3, IRKN4, DPRKN6, DPRKN8, DPRKN12, ERKN4 + IRKN3, IRKN4, DPRKN6, DPRKN8, DPRKN12, ERKN4, ERKN5 end # module diff --git a/src/alg_utils.jl b/src/alg_utils.jl index c27627d49d..1ebe26e614 100644 --- a/src/alg_utils.jl +++ b/src/alg_utils.jl @@ -102,6 +102,7 @@ isfsal(alg::DPRKN6) = true isfsal(alg::DPRKN8) = true isfsal(alg::DPRKN12) = true isfsal(alg::ERKN4) = true +isfsal(alg::ERKN5) = true fsal_typeof(alg::OrdinaryDiffEqAlgorithm,rate_prototype) = typeof(rate_prototype) #fsal_typeof(alg::LawsonEuler,rate_prototype) = Vector{typeof(rate_prototype)} @@ -269,6 +270,7 @@ alg_order(alg::DPRKN6) = 6 alg_order(alg::DPRKN8) = 8 alg_order(alg::DPRKN12) = 12 alg_order(alg::ERKN4) = 4 +alg_order(alg::ERKN5) = 5 alg_order(alg::Midpoint) = 2 alg_order(alg::GenericIIF1) = 1 diff --git a/src/algorithms.jl b/src/algorithms.jl index c71ca9f994..d84cad4219 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -133,6 +133,7 @@ struct DPRKN6 <: OrdinaryDiffEqAdaptiveAlgorithm end struct DPRKN8 <: OrdinaryDiffEqAdaptiveAlgorithm end struct DPRKN12 <: OrdinaryDiffEqAdaptiveAlgorithm end struct ERKN4 <: OrdinaryDiffEqAdaptiveAlgorithm end +struct ERKN5 <: OrdinaryDiffEqAdaptiveAlgorithm end ################################################################################ diff --git a/src/caches/rkn_caches.jl b/src/caches/rkn_caches.jl index 94d479484f..00342bcf6f 100644 --- a/src/caches/rkn_caches.jl +++ b/src/caches/rkn_caches.jl @@ -294,3 +294,34 @@ function alg_cache(alg::ERKN4,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev tmp = similar(u) ERKN4Cache(u,uprev,k1,k2,k3,k4,k,utilde,tmp,atmp,tab) end + +struct ERKN5Cache{uType,uArrayType,rateType,reducedRateType,uEltypeNoUnits,TabType} <: OrdinaryDiffEqMutableCache + u::uType + uprev::uType + fsalfirst::rateType + k2::reducedRateType + k3::reducedRateType + k4::reducedRateType + k::rateType + utilde::uArrayType + tmp::uType + atmp::uEltypeNoUnits + tab::TabType +end + +u_cache(c::ERKN5Cache) = (c.atmp,c.utilde) +du_cache(c::ERKN5Cache) = (c.fsalfirst,c.k2,c.k3,c.k4,c.k) + +function alg_cache(alg::ERKN5,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,::Type{Val{true}}) + reduced_rate_prototype = rate_prototype.x[2] + tab = ERKN5ConstantCache(real(uEltypeNoUnits),real(tTypeNoUnits)) + k1 = zeros(rate_prototype) + k2 = zeros(reduced_rate_prototype) + k3 = zeros(reduced_rate_prototype) + k4 = zeros(reduced_rate_prototype) + k = zeros(rate_prototype) + utilde = similar(u,indices(u)) + atmp = similar(u,uEltypeNoUnits) + tmp = similar(u) + ERKN5Cache(u,uprev,k1,k2,k3,k4,k,utilde,tmp,atmp,tab) +end From c166d05f8649abebc44d7c193a7a7f6ba5ad2f0f Mon Sep 17 00:00:00 2001 From: YingboMa Date: Tue, 12 Sep 2017 13:05:27 -0400 Subject: [PATCH 4/5] Add ERKN5 algorithm --- src/integrators/rkn_integrators.jl | 41 ++++++++++++++++++++++++++++-- test/partitioned_methods_tests.jl | 9 +++++-- 2 files changed, 46 insertions(+), 4 deletions(-) diff --git a/src/integrators/rkn_integrators.jl b/src/integrators/rkn_integrators.jl index ff7b22cc9e..16f544d9b7 100644 --- a/src/integrators/rkn_integrators.jl +++ b/src/integrators/rkn_integrators.jl @@ -7,9 +7,9 @@ const NystromDefaultInitialization = Union{Nystrom4Cache, Nystrom4VelocityIndependentCache, Nystrom5VelocityIndependentCache, - IRKN3Cache, IRKN4Cache, + IRKN3Cache, IRKN4Cache, DPRKN8Cache, DPRKN12Cache, - ERKN4Cache} + ERKN4Cache, ERKN5Cache} function initialize!(integrator,cache::NystromDefaultInitialization) @unpack fsalfirst,k = cache @@ -492,3 +492,40 @@ end integrator.EEst = integrator.opts.internalnorm(atmp) end end + +@muladd function perform_step!(integrator,cache::ERKN5Cache,repeat_step=false) + @unpack t,dt,f = integrator + u,du = integrator.u.x + uprev,duprev = integrator.uprev.x + @unpack tmp,atmp,fsalfirst,k2,k3,k4,k,utilde = cache + @unpack c1, c2, c3, a21, a31, a32, a41, a42, a43, b1, b2, b3, b4, bp1, bp2, bp3, bp4, btilde1, btilde2, btilde3, btilde4 = cache.tab + ku, kdu = integrator.cache.tmp.x[1], integrator.cache.tmp.x[2] + uidx = eachindex(integrator.uprev.x[2]) + k1 = fsalfirst.x[2] + + @. ku = uprev + dt*(c1*duprev + dt*a21*k1) + + f.f2(t+dt*c1,ku,du,k2) + @. ku = uprev + dt*(c2*duprev + dt*(a31*k1 + a32*k2)) + + f.f2(t+dt*c2,ku,du,k3) + @. ku = uprev + dt*(c3*duprev + dt*(a41*k1 + a42*k2 + a43*k3)) + + f.f2(t+dt*c3,ku,du,k4) + @tight_loop_macros for i in uidx + @inbounds u[i] = uprev[i] + dt*(duprev[i] + dt*(b1 *k1[i] + b2 *k2[i] + b3 *k3[i] + b4 *k4[i])) + @inbounds du[i] = duprev[i] + dt*(bp1*k1[i] + bp2*k2[i] + bp3*k3[i] + bp4*k4[i]) + end + + f.f1(t+dt,u,du,k.x[1]) + f.f2(t+dt,u,du,k.x[2]) + if integrator.opts.adaptive + uhat, duhat = utilde.x + dtsq = dt^2 + @tight_loop_macros for i in uidx + @inbounds uhat[i] = dtsq*(btilde1*k1[i] + btilde2*k2[i] + btilde3*k3[i] + btilde4*k4[i]) + end + calculate_residuals!(atmp, uhat, integrator.uprev, integrator.u, integrator.opts.abstol, integrator.opts.reltol) + integrator.EEst = integrator.opts.internalnorm(atmp) + end +end diff --git a/test/partitioned_methods_tests.jl b/test/partitioned_methods_tests.jl index 48b0b1b538..7b65341e61 100644 --- a/test/partitioned_methods_tests.jl +++ b/test/partitioned_methods_tests.jl @@ -144,6 +144,9 @@ sim = test_convergence(dts,prob_big,DPRKN12(),dense_errors=true) sim = test_convergence(dts,prob_big,ERKN4(),dense_errors=true) @test sim.𝒪est[:l2] ≈ 5 rtol = 1e-1 @test sim.𝒪est[:L2] ≈ 4 rtol = 1e-1 +sim = test_convergence(dts,prob_big,ERKN5(),dense_errors=true) +@test sim.𝒪est[:l2] ≈ 5 rtol = 1e-1 +@test sim.𝒪est[:L2] ≈ 4 rtol = 1e-1 # Adaptive methods regression test sol = solve(prob, DPRKN6()) @@ -152,8 +155,10 @@ sol = solve(prob, DPRKN8()) @test length(sol.u) < 13 sol = solve(prob, DPRKN12()) @test length(sol.u) < 9 -sol = solve(prob, ERKN4()) -@test length(sol.u) < 15 +sol = solve(prob, ERKN4(),reltol=1e-8) +@test length(sol.u) < 38 +sol = solve(prob, ERKN5(),reltol=1e-8) +@test length(sol.u) < 29 f = function (t,u,du) du.x[1] .= u.x[2] From 1abeb37e5ab82054c9f1152f66f0bd6a9956d447 Mon Sep 17 00:00:00 2001 From: YingboMa Date: Tue, 12 Sep 2017 15:04:32 -0400 Subject: [PATCH 5/5] Mark Int64 for long integers --- src/tableaus/rkn_tableaus.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/tableaus/rkn_tableaus.jl b/src/tableaus/rkn_tableaus.jl index b346936d7b..381566c50c 100644 --- a/src/tableaus/rkn_tableaus.jl +++ b/src/tableaus/rkn_tableaus.jl @@ -262,9 +262,9 @@ function ERKN5ConstantCache(T::Type,T2::Type) a21 = T(1//8) a31 = T(2907//343000) a32 = T(1216//42875) - a41 = T(6624772//128538819) - a42 = T(6273905//54121608) - a43 = T(210498365//1028310552) + a41 = T(6624772//Int64(128538819)) + a42 = T(6273905//Int64(54121608)) + a43 = T(Int64(210498365)//Int64(1028310552)) b1 = T(479//5016) b2 = T(235//1776) b3 = T(145775//641744) @@ -272,7 +272,7 @@ function ERKN5ConstantCache(T::Type,T2::Type) btilde1 = T(479//5016 - 184883//2021250) btilde2 = T(235//1776 - 411163//3399375) btilde3 = T(145775//641744 - 6//25) - btilde4 = T(309519//6873416 - 593028//12464375) + btilde4 = T(309519//6873416 - 593028//Int64(12464375)) bp1 = b1 bp2 = T(235//888) bp3 = T(300125//962616)