Merge branch 'nystrom'

src/OrdinaryDiffEq.jl

@@ -133,5 +133,5 @@ module OrdinaryDiffEq
 
   export SplitEuler
 
-  export Nystrom4, Nystrom4VelocityIndependent
+  export Nystrom4, Nystrom4VelocityIndependent, Nystrom5VelocityIndependent
 end # module

src/alg_utils.jl

@@ -69,6 +69,7 @@ isfsal(alg::SofSpa10) = true
 
 isfsal(alg::Nystrom4) = true
 isfsal(alg::Nystrom4VelocityIndependent) = true
+isfsal(alg::Nystrom5VelocityIndependent) = true
 
 fsal_typeof(alg::OrdinaryDiffEqAlgorithm,rate_prototype) = typeof(rate_prototype)
 #fsal_typeof(alg::LawsonEuler,rate_prototype) = Vector{typeof(rate_prototype)}
@@ -171,6 +172,7 @@ alg_order(alg::SofSpa10) = 10
 
 alg_order(alg::Nystrom4) = 4
 alg_order(alg::Nystrom4VelocityIndependent) = 4
+alg_order(alg::Nystrom5VelocityIndependent) = 5
 
 alg_order(alg::Midpoint) = 2
 alg_order(alg::IIF1) = 1

src/algorithms.jl

@@ -65,6 +65,7 @@ struct SofSpa10 <: OrdinaryDiffEqAlgorithm end
 
 struct Nystrom4 <: OrdinaryDiffEqAlgorithm end
 struct Nystrom4VelocityIndependent <: OrdinaryDiffEqAlgorithm end
+struct Nystrom5VelocityIndependent <: OrdinaryDiffEqAlgorithm end
 
 ################################################################################
 

src/caches/rkn_caches.jl

@@ -47,3 +47,27 @@ function alg_cache(alg::Nystrom4VelocityIndependent,u,rate_prototype,uEltypeNoUn
   tmp = similar(u)
   Nystrom4VelocityIndependentCache(u,uprev,k₁,k₂,k₃,k,tmp)
 end
+
+struct Nystrom5VelocityIndependentCache{uType,rateType} <: OrdinaryDiffEqMutableCache
+  u::uType
+  uprev::uType
+  fsalfirst::rateType
+  k₂::rateType
+  k₃::rateType
+  k₄::rateType
+  k::rateType
+  tmp::uType
+end
+
+u_cache(c::Nystrom5VelocityIndependentCache) = ()
+du_cache(c::Nystrom5VelocityIndependentCache) = (c.fsalfirst,c.k₂,c.k₃,c.k₄,c.k)
+
+function alg_cache(alg::Nystrom5VelocityIndependent,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,reltol,::Type{Val{true}})
+  k₁ = zeros(rate_prototype)
+  k₂ = zeros(rate_prototype)
+  k₃ = zeros(rate_prototype)
+  k₄ = zeros(rate_prototype)
+  k  = zeros(rate_prototype)
+  tmp = similar(u)
+  Nystrom5VelocityIndependentCache(u,uprev,k₁,k₂,k₃,k₄,k,tmp)
+end

src/integrators/rkn_integrators.jl

@@ -90,6 +90,7 @@ end
   half_dtsq = dtsq/2
   ttmp = t+halfdt
 
+  f[2](ttmp,uprev,duprev,k₁.x[2])
   @tight_loop_macros for i in uidx
     ## y₁ = y₀ + hy'₀ + h²∑b̄ᵢk'ᵢ
     @inbounds ku[i] = @muladd uprev[i] + halfdt*duprev[i] + eighth_dtsq*k₁.x[2][i]
@@ -108,3 +109,52 @@ end
   f[1](t+dt,ku,du,k.x[1])
   f[2](t+dt,ku,du,k.x[2])
 end
+
+
+@inline function initialize!(integrator,cache::Nystrom5VelocityIndependentCache,f=integrator.f)
+  @unpack tmp,fsalfirst,k₂,k = cache
+  uprev,duprev = integrator.uprev.x
+
+  integrator.fsalfirst = fsalfirst
+  integrator.fsallast = k
+  integrator.kshortsize = 2
+  integrator.k = eltype(integrator.sol.k)(integrator.kshortsize)
+  integrator.k[1] = integrator.fsalfirst
+  integrator.k[2] = integrator.fsallast
+  f[1](integrator.t,uprev,duprev,integrator.k[2].x[1])
+  f[2](integrator.t,uprev,duprev,integrator.k[2].x[2])
+end
+
+@inline function perform_step!(integrator,cache::Nystrom5VelocityIndependentCache,f=integrator.f)
+  @unpack t,dt,k = integrator
+  u,du = integrator.u.x
+  uprev,duprev = integrator.uprev.x
+  uidx = eachindex(integrator.uprev.x[1])
+  @unpack tmp,fsalfirst,k₂,k₃,k₄,k = cache
+  ku, kdu = integrator.cache.tmp.x[1], integrator.cache.tmp.x[2]
+  k₁ = fsalfirst
+  dtsq = dt^2
+
+  f[2](t+1//5*dt,uprev,duprev,k₁.x[2])
+  @tight_loop_macros for i in uidx
+    @inbounds ku[i] = @muladd uprev[i] + (1//5*dt)*duprev[i] + (1//50*dtsq)*k₁.x[2][i]
+  end
+
+  f[2](t+1//5*dt,ku,du,k₂.x[2])
+  @tight_loop_macros for i in uidx
+    @inbounds ku[i] = @muladd uprev[i] + (2//3*dt)*duprev[i] + (-1//27*dtsq)*k₁.x[2][i] + (7//27*dtsq)*k₂.x[2][i]
+  end
+
+  f[2](t+2//3*dt,ku,du,k₃.x[2])
+  @tight_loop_macros for i in uidx
+    @inbounds ku[i] = @muladd uprev[i] + dt*duprev[i] + (3//10*dtsq)*k₁.x[2][i] + (-2//35*dtsq)*k₂.x[2][i] + (9//35*dtsq)*k₃.x[2][i]
+  end
+
+  f[2](t+dt,ku,du,k₄.x[2])
+  @tight_loop_macros for i in uidx
+    @inbounds u[i]  = @muladd uprev[i] + dt*duprev[i] + (14//336*dtsq)*k₁.x[2][i] + (100//336*dtsq)*k₂.x[2][i] + (54//336*dtsq)*k₃.x[2][i]
+    @inbounds du[i] = @muladd duprev[i] + (14//336*dt)*k₁.x[2][i] + (125//336*dt)*k₂.x[2][i] + (162//336*dt)*k₃.x[2][i] + (35//336*dt)*k₄.x[2][i]
+  end
+  f[1](t+dt,ku,du,k.x[1])
+  f[2](t+dt,ku,du,k.x[2])
+end

test/partitioned_methods_tests.jl

@@ -110,8 +110,12 @@ sim = test_convergence(dts,prob,Nystrom4(),dense_errors=true)
 @test sim.𝒪est[:l2] ≈ 4 rtol = 1e-1
 @test sim.𝒪est[:L2] ≈ 4 rtol = 1e-1
 sim = test_convergence(dts,prob,Nystrom4VelocityIndependent(),dense_errors=true)
-@test_broken sim.𝒪est[:l2] ≈ 4 rtol = 1e-1
-@test_broken sim.𝒪est[:L2] ≈ 4 rtol = 1e-1
+@test sim.𝒪est[:l2] ≈ 4 rtol = 1e-1
+@test sim.𝒪est[:L2] ≈ 4 rtol = 1e-1
+dts = 1.0./2.0.^(5:-1:0)
+sim = test_convergence(dts,prob,Nystrom5VelocityIndependent(),dense_errors=true)
+@test sim.𝒪est[:l2] ≈ 5 rtol = 1e-1
+@test sim.𝒪est[:L2] ≈ 5 rtol = 1e-1
 
 dts = 1.0./2.0.^(2:-1:-2)
 sim = test_convergence(dts,prob,SofSpa10(),dense_errors=true)