From de52b9d122c2480ee803a1733114936a7648841b Mon Sep 17 00:00:00 2001 From: Henry Langner <128313572+HenryLangner@users.noreply.github.com> Date: Fri, 4 Aug 2023 09:27:48 +0000 Subject: [PATCH 1/2] Dense Output for FineRKN5, with high order hermite polynomial --- src/caches/rkn_caches.jl | 29 ++-- src/dense/interpolants.jl | 142 ++++++++++++++++++ src/interp_func.jl | 7 + src/perform_step/rkn_perform_step.jl | 105 +++++++++++-- src/tableaus/rkn_tableaus.jl | 211 ++++++++++++++++++++++++++- 5 files changed, 462 insertions(+), 32 deletions(-) diff --git a/src/caches/rkn_caches.jl b/src/caches/rkn_caches.jl index 7ffa67bf23..301dfb8060 100644 --- a/src/caches/rkn_caches.jl +++ b/src/caches/rkn_caches.jl @@ -78,24 +78,6 @@ function alg_cache(alg::FineRKN4, u, rate_prototype, ::Type{uEltypeNoUnits}, FineRKN4ConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits)) end -@cache struct FineRKN5Cache{uType, rateType, reducedRateType, uNoUnitsType, TabType} <: - OrdinaryDiffEqMutableCache - u::uType - uprev::uType - fsalfirst::rateType - k2::reducedRateType - k3::reducedRateType - k4::reducedRateType - k5::reducedRateType - k6::reducedRateType - k7::reducedRateType - k::rateType - utilde::uType - tmp::uType - atmp::uNoUnitsType - tab::TabType -end - function alg_cache(alg::FineRKN5, u, rate_prototype, ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, dt, reltol, p, calck, @@ -109,12 +91,21 @@ function alg_cache(alg::FineRKN5, u, rate_prototype, ::Type{uEltypeNoUnits}, k5 = zero(reduced_rate_prototype) k6 = zero(reduced_rate_prototype) k7 = zero(reduced_rate_prototype) + k8 = zero(reduced_rate_prototype) + k9 = zero(reduced_rate_prototype) k = zero(rate_prototype) utilde = zero(u) atmp = similar(u, uEltypeNoUnits) recursivefill!(atmp, false) tmp = zero(u) - FineRKN5Cache(u, uprev, k1, k2, k3, k4, k5, k6, k7, k, utilde, tmp, atmp, tab) + FineRKN5Cache(u, uprev, k1, k2, k3, k4, k5, k6, k7, k8, k9, k, utilde, tmp, atmp, tab) +end + +function alg_cache(alg::FineRKN5, u, rate_prototype, ::Type{uEltypeNoUnits}, + ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, + dt, reltol, p, calck, + ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} + FineRKN5ConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits)) end function alg_cache(alg::FineRKN5, u, rate_prototype, ::Type{uEltypeNoUnits}, diff --git a/src/dense/interpolants.jl b/src/dense/interpolants.jl index 50067ce777..8e9ee8a747 100644 --- a/src/dense/interpolants.jl +++ b/src/dense/interpolants.jl @@ -2844,6 +2844,148 @@ end """ +""" +@def finerkn5unpack begin + if typeof(cache) <: OrdinaryDiffEqMutableCache + @unpack r12, r13, r14, r15, r16, r33, r34, + r35, r36, r43, + r44, r45, r46, r53, r54, r55, r56, r63, r64, r65, r66, r73, r74, r75, r76, r83, r84, + r85, r86, r93, r94, r95, r96, rp11, rp12, rp13, rp14, rp15, rp32, rp33, rp34, rp35, + rp42, rp43, rp44, rp45, rp52, rp53, rp54, rp55, rp62, rp63, rp64, rp65, rp72, rp73, + rp74, rp75, rp82, rp83, rp84, rp85, rp92, rp93, rp94, rp95 = cache.tab + else + @unpack r12, r13, r14, r15, r16, r33, r34, + r35, r36, r43, + r44, r45, r46, r53, r54, r55, r56, r63, r64, r65, r66, r73, r74, r75, r76, r83, r84, + r85, r86, r93, r94, r95, r96, rp11, rp12, rp13, rp14, rp15, rp32, rp33, rp34, rp35, + rp42, rp43, rp44, rp45, rp52, rp53, rp54, rp55, rp62, rp63, rp64, rp65, rp72, rp73, + rp74, rp75, rp82, rp83, rp84, rp85, rp92, rp93, rp94, rp95 = cache + end +end + +@def finerkn5pre0 begin + @finerkn5unpack + b1Θ = Θ^2 * @evalpoly(Θ, r12, r13, r14, r15, r16) + b3Θ = Θ^3 * @evalpoly(Θ, r33, r34, r35, r36) + b4Θ = Θ^3 * @evalpoly(Θ, r43, r44, r45, r46) + b5Θ = Θ^3 * @evalpoly(Θ, r53, r54, r55, r56) + b6Θ = Θ^3 * @evalpoly(Θ, r63, r64, r65, r66) + b7Θ = Θ^3 * @evalpoly(Θ, r73, r74, r75, r76) + b8Θ = Θ^3 * @evalpoly(Θ, r83, r84, r85, r86) + b9Θ = Θ^3 * @evalpoly(Θ, r93, r94, r95, r96) + + bp1Θ = Θ * @evalpoly(Θ, rp11, rp12, rp13, rp14, rp15) + bp3Θ = Θ^2 * @evalpoly(Θ, rp32, rp33, rp34, rp35) + bp4Θ = Θ^2 * @evalpoly(Θ, rp42, rp43, rp44, rp45) + bp5Θ = Θ^2 * @evalpoly(Θ, rp52, rp53, rp54, rp55) + bp6Θ = Θ^2 * @evalpoly(Θ, rp62, rp63, rp64, rp65) + bp7Θ = Θ^2 * @evalpoly(Θ, rp72, rp73, rp74, rp75) + bp8Θ = Θ^2 * @evalpoly(Θ, rp82, rp83, rp84, rp85) + bp9Θ = Θ^2 * @evalpoly(Θ, rp92, rp93, rp94, rp95) + + kk1, kk2, kk3, kk4 = k + k1, k3 = kk1.x + k4, k5 = kk2.x + k6, k7 = kk3.x + k8, k9 = kk4.x + + duprev, uprev = y₀.x + dtsq = dt^2 +end + +@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, + cache::Union{FineRKN5ConstantCache, FineRKN5Cache}, + idxs::Nothing, T::Type{Val{0}}) + @finerkn5pre0 + return ArrayPartition(duprev + + dt * Θ * + (bp1Θ * k1 + bp3Θ * k3 + + bp4Θ * k4 + bp5Θ * k5 + bp6Θ * k6 + bp7Θ * k7 + bp8Θ * k8 + + bp9Θ * k9), + uprev + + dt * Θ * + (duprev + + dt * Θ * + (b1Θ * k1 + b3Θ * k3 + + b4Θ * k4 + b5Θ * k5 + b6Θ * k6 + b7Θ * k7 + b8Θ * k8 + b9Θ * k9))) +end + +@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, + cache::Union{FineRKN5ConstantCache, FineRKN5Cache}, idxs, + T::Type{Val{0}}) + @finerkn5pre0 + return ArrayPartition(duprev[idxs] + + dt * Θ * + (bp1Θ * k1[idxs] + bp3Θ * k3[idxs] + + bp4Θ * k4[idxs] + bp5Θ * k5[idxs] + bp6Θ * k6[idxs] + + bp7Θ * k7[idxs] + bp8Θ * k8[idxs] + bp9Θ * k9[idxs]), + uprev[idxs] + + dt * Θ * + (duprev[idxs] + + dt * Θ * + (b1Θ * k1[idxs] + + b3Θ * k3[idxs] + + b4Θ * k4[idxs] + b5Θ * k5[idxs] + b6Θ * k6[idxs] + b7Θ * k7[idxs] + + b8Θ * k8[idxs] + b9Θ * k9[idxs]))) +end + +@muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, + cache::Union{FineRKN5ConstantCache, FineRKN5Cache}, + idxs::Nothing, T::Type{Val{0}}) + @finerkn5pre0 + @inbounds @.. broadcast=false out.x[2]=uprev + + dt * Θ * + (duprev + + dt * Θ * + (b1Θ * k1 + + b3Θ * k3 + + b4Θ * k4 + b5Θ * k5 + b6Θ * k6 + b7Θ * k7 + + b8Θ * k8 + b9Θ * k9)) + @inbounds @.. broadcast=false out.x[1]=duprev + + dt * Θ * + (bp1Θ * k1 + bp3Θ * k3 + + bp4Θ * k4 + bp5Θ * k5 + bp6Θ * k6 + bp7Θ * k7 + + bp8Θ * k8 + bp9Θ * k9) + #for i in eachindex(out.x[1]) + # out.x[2][i] = uprev[i] + dt*Θ*(duprev[i] + dt*Θ*(b1Θ*k1[i] + + # b3Θ*k3[i] + + # b4Θ*k4[i] + b5Θ*k5[i] + b6Θ*k6[i])) + # out.x[1][i] = duprev[i] + dt*Θ*(bp1Θ*k1[i] + bp3Θ*k3[i] + + # bp4Θ*k4[i] + bp5Θ*k5[i] + bp6Θ*k6[i]) + #end + out +end + +@muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, + cache::Union{FineRKN5ConstantCache, FineRKN5Cache}, idxs, + T::Type{Val{0}}) + @finerkn5pre0 + @inbounds @.. broadcast=false out.x[2]=uprev[idxs] + + dt * Θ * + (duprev[idxs] + + dt * Θ * + (b1Θ * k1[idxs] + + b3Θ * k3[idxs] + + b4Θ * k4[idxs] + b5Θ * k5[idxs] + + b6Θ * k6[idxs] + b7Θ * k7[idxs] + + b8Θ * k8[idxs] + b9Θ * k9[idxs])) + @inbounds @.. broadcast=false out.x[1]=duprev[idxs] + + dt * Θ * + (bp1Θ * k1[idxs] + bp3Θ * k3[idxs] + + bp4Θ * k4[idxs] + bp5Θ * k5[idxs] + + bp6Θ * k6[idxs] + bp7Θ * k7[idxs] + + bp8Θ * k8[idxs] + bp9Θ * k9[idxs]) + #for (j,i) in enumerate(idxs) + # out.x[2][j] = uprev[i] + dt*Θ*(duprev[i] + dt*Θ*(b1Θ*k1[i] + + # b3Θ*k3[i] + + # b4Θ*k4[i] + b5Θ*k5[i] + b6Θ*k6[i])) + # out.x[1][j] = duprev[i] + dt*Θ*(bp1Θ*k1[i] + bp3Θ*k3[i] + + # bp4Θ*k4[i] + bp5Θ*k5[i] + bp6Θ*k6[i]) + #end + out +end +""" + """ @def dprkn6unpack begin if typeof(cache) <: OrdinaryDiffEqMutableCache diff --git a/src/interp_func.jl b/src/interp_func.jl index 0f43051a87..063d3cf3e2 100644 --- a/src/interp_func.jl +++ b/src/interp_func.jl @@ -156,6 +156,13 @@ function DiffEqBase.interp_summary(::Type{cacheType}, caches = fieldtype(cacheType, :caches) join([DiffEqBase.interp_summary(ct, dense) for ct in fieldtypes(caches)], ", ") end +function DiffEqBase.interp_summary(::Type{cacheType}, + dense::Bool) where { + cacheType <: + Union{FineRKN5ConstantCache, + FineRKN5Cache}} + dense ? "specialized 6th order interpolation" : "1st order linear" +end function DiffEqBase.interp_summary(::Type{cacheType}, dense::Bool) where { cacheType <: diff --git a/src/perform_step/rkn_perform_step.jl b/src/perform_step/rkn_perform_step.jl index 4d93f3c9c6..ffe5b56e82 100644 --- a/src/perform_step/rkn_perform_step.jl +++ b/src/perform_step/rkn_perform_step.jl @@ -5,7 +5,6 @@ ## y'₁ = y'₀ + h∑bᵢk'ᵢ const NystromCCDefaultInitialization = Union{Nystrom4ConstantCache, FineRKN4ConstantCache, - FineRKN5ConstantCache, Nystrom4VelocityIndependentConstantCache, Nystrom5VelocityIndependentConstantCache, IRKN3ConstantCache, IRKN4ConstantCache, @@ -26,7 +25,7 @@ function initialize!(integrator, cache::NystromCCDefaultInitialization) integrator.fsalfirst = ArrayPartition((kdu, ku)) end -const NystromDefaultInitialization = Union{Nystrom4Cache, FineRKN4Cache, FineRKN5Cache, +const NystromDefaultInitialization = Union{Nystrom4Cache, FineRKN4Cache, Nystrom4VelocityIndependentCache, Nystrom5VelocityIndependentCache, IRKN3Cache, IRKN4Cache, @@ -232,11 +231,36 @@ end end end +function initialize!(integrator, cache::FineRKN5ConstantCache) + duprev, uprev = integrator.uprev.x + integrator.kshortsize = 4 + integrator.k = typeof(integrator.k)(undef, integrator.kshortsize) + + kdu = integrator.f.f1(duprev, uprev, integrator.p, integrator.t) + ku = integrator.f.f2(duprev, uprev, integrator.p, integrator.t) + integrator.stats.nf += 1 + integrator.stats.nf2 += 1 + integrator.fsalfirst = ArrayPartition((kdu, ku)) + integrator.fsallast = zero(integrator.fsalfirst) + + integrator.k[1] = integrator.fsalfirst + @inbounds for i in 2:(integrator.kshortsize - 1) + integrator.k[i] = zero(integrator.fsalfirst) + end + integrator.k[integrator.kshortsize] = integrator.fsallast +end + @muladd function perform_step!(integrator, cache::FineRKN5ConstantCache, repeat_step = false) @unpack t, dt, f, p = integrator duprev, uprev = integrator.uprev.x - @unpack c2, c3, c4, c5, c6, c7, a21, a31, a32, a41, a43, a51, a52, a53, a54, a61, a62, a63, a64, a71, a73, a74, a75, abar21, abar31, abar32, abar41, abar42, abar43, abar51, abar52, abar53, abar54, abar61, abar62, abar63, abar64, abar65, abar71, abar73, abar74, abar75, abar76, b1, b3, b4, b5, bbar1, bbar3, bbar4, bbar5, bbar6, btilde1, btilde3, btilde4, btilde5, bptilde1, bptilde3, bptilde4, bptilde5, bptilde6, bptilde7 = cache + @unpack c2, c3, c4, c5, c6, c7, c8, c9, a21, a31, a32, a41, a43, a51, + a52, a53, a54, a61, a62, a63, a64, a71, a73, a74, a75, a81, a82, a83, a84, a85, a87, a91, a92, a93, a94, a95, a97, + abar21, abar31, abar32, abar41, abar42, abar43, abar51, + abar52, abar53, abar54, abar61, abar62, abar63, abar64, abar65, + abar71, abar73, abar74, abar75, abar76, abar81, abar82, abar83, abar84, abar85, abar87, abar91, abar93, abar94, abar95, abar97, b1, b3, b4, + b5, bbar1, bbar3, bbar4, bbar5, bbar6, btilde1, btilde3, btilde4, btilde5, bptilde1, + bptilde3, bptilde4, bptilde5, bptilde6, bptilde7 = cache k1 = integrator.fsalfirst.x[1] ku = uprev + dt * (c2 * duprev + dt * (a21 * k1)) @@ -272,12 +296,31 @@ end u = uprev + dt * (duprev + dt * (b1 * k1 + b3 * k3 + b4 * k4 + b5 * k5)) # no b6, b7 du = duprev + dt * (bbar1 * k1 + bbar3 * k3 + bbar4 * k4 + bbar5 * k5 + bbar6 * k6) # no b2, b7 + ku = uprev + + dt * (c8 * duprev + + dt * (a81 * k1 + a82 * k2 + a83 * k3 + a84 * k4 + a85 * k5 + a87 * k7)) # a86 = 0 + kdu = duprev + + dt * (abar81 * k1 + abar82 * k2 + abar83 * k3 + abar84 * k4 + abar85 * k5 + + abar87 * k7) # abar86 = 0 + + k8 = f.f1(kdu, ku, p, t + dt * c8) + ku = uprev + + dt * (c9 * duprev + + dt * (a91 * k1 + a92 * k2 + a93 * k3 + a94 * k4 + a95 * k5 + a97 * k7)) # a96 = a98 = 0 + kdu = duprev + + dt * (abar91 * k1 + abar93 * k3 + abar94 * k4 + abar95 * k5 + + abar97 * k7) # abar92 = abar96 = abar98 = 0 + + k9 = f.f1(kdu, ku, p, t + dt * c9) + integrator.u = ArrayPartition((du, u)) integrator.fsallast = ArrayPartition((f.f1(du, u, p, t + dt), f.f2(du, u, p, t + dt))) - integrator.stats.nf += 7 + integrator.stats.nf += 9 integrator.stats.nf2 += 1 - integrator.k[1] = integrator.fsalfirst - integrator.k[2] = integrator.fsallast + integrator.k[1] = ArrayPartition(integrator.fsalfirst.x[1], k3) + integrator.k[2] = ArrayPartition(k4, k5) + integrator.k[3] = ArrayPartition(k6, k7) + integrator.k[4] = ArrayPartition(k8, k9) if integrator.opts.adaptive dtsq = dt^2 @@ -292,12 +335,36 @@ end end end +function initialize!(integrator, cache::FineRKN5Cache) + @unpack fsalfirst, k = cache + duprev, uprev = integrator.uprev.x + + integrator.fsalfirst = fsalfirst + integrator.fsallast = k + integrator.kshortsize = 4 + resize!(integrator.k, integrator.kshortsize) + integrator.k[1] = ArrayPartition(cache.fsalfirst.x[1], cache.k3) + integrator.k[2] = ArrayPartition(cache.k4, cache.k5) + integrator.k[3] = ArrayPartition(cache.k6, cache.k7) + integrator.k[4] = ArrayPartition(cache.k8, cache.k9) + integrator.f.f1(integrator.fsallast.x[1], duprev, uprev, integrator.p, integrator.t) + integrator.f.f2(integrator.fsallast.x[2], duprev, uprev, integrator.p, integrator.t) + integrator.stats.nf += 1 + integrator.stats.nf2 += 1 +end + @muladd function perform_step!(integrator, cache::FineRKN5Cache, repeat_step = false) @unpack t, dt, f, p = integrator du, u = integrator.u.x duprev, uprev = integrator.uprev.x - @unpack tmp, atmp, fsalfirst, k2, k3, k4, k5, k6, k7, k, utilde = cache - @unpack c1, c2, c3, c4, c5, c6, c7, a21, a31, a32, a41, a43, a51, a52, a53, a54, a61, a62, a63, a64, a71, a73, a74, a75, abar21, abar31, abar32, abar41, abar42, abar43, abar51, abar52, abar53, abar54, abar61, abar62, abar63, abar64, abar65, abar71, abar73, abar74, abar75, abar76, b1, b3, b4, b5, bbar1, bbar3, bbar4, bbar5, bbar6, btilde1, btilde3, btilde4, btilde5, bptilde1, bptilde3, bptilde4, bptilde5, bptilde6, bptilde7 = cache.tab + @unpack tmp, atmp, fsalfirst, k2, k3, k4, k5, k6, k7, k8, k9, k, utilde = cache + @unpack c2, c3, c4, c5, c6, c7, c8, c9, a21, a31, a32, a41, a43, a51, + a52, a53, a54, a61, a62, a63, a64, a71, a73, a74, a75, a81, a82, a83, a84, a85, a87, a91, a92, a93, a94, a95, a97, + abar21, abar31, abar32, abar41, abar42, abar43, abar51, + abar52, abar53, abar54, abar61, abar62, abar63, abar64, abar65, + abar71, abar73, abar74, abar75, abar76, abar81, abar82, abar83, abar84, abar85, abar87, abar91, abar93, abar94, abar95, abar97, b1, b3, b4, + b5, bbar1, bbar3, bbar4, bbar5, bbar6, btilde1, btilde3, btilde4, btilde5, bptilde1, + bptilde3, bptilde4, bptilde5, bptilde6, bptilde7 = cache.tab kdu, ku = integrator.cache.tmp.x[1], integrator.cache.tmp.x[2] uidx = eachindex(integrator.uprev.x[2]) k1 = integrator.fsalfirst.x[1] @@ -344,9 +411,29 @@ end dt * (bbar1 * k1 + bbar3 * k3 + bbar4 * k4 + bbar5 * k5 + bbar6 * k6) + @.. broadcast=false ku=uprev + + dt * (c8 * duprev + + dt * (a81 * k1 + a82 * k2 + a83 * k3 + a84 * k4 + a85 * k5 + + a87 * k7)) # a86 = 0 + @.. broadcast=false kdu=duprev + + dt * (abar81 * k1 + abar82 * k2 + abar83 * k3 + abar84 * k4 + + abar85 * k5 + + abar87 * k7) # abar86 = 0 + + f.f1(k8, kdu, ku, p, t + dt * c8) + @.. broadcast=false ku=uprev + + dt * (c9 * duprev + + dt * (a91 * k1 + a92 * k2 + a93 * k3 + a94 * k4 + a95 * k5 + + a97 * k7)) # a96 = a98 = 0 + @.. broadcast=false kdu=duprev + + dt * (abar91 * k1 + abar93 * k3 + abar94 * k4 + abar95 * k5 + + abar97 * k7) # abar92 = abar96 = abar98 = 0 + + f.f1(k9, kdu, ku, p, t + dt * c9) + f.f1(k.x[1], du, u, p, t + dt) f.f2(k.x[2], du, u, p, t + dt) - integrator.stats.nf += 7 + integrator.stats.nf += 9 integrator.stats.nf2 += 1 if integrator.opts.adaptive duhat, uhat = utilde.x diff --git a/src/tableaus/rkn_tableaus.jl b/src/tableaus/rkn_tableaus.jl index d4e89d2a15..3463424e65 100644 --- a/src/tableaus/rkn_tableaus.jl +++ b/src/tableaus/rkn_tableaus.jl @@ -107,6 +107,8 @@ struct FineRKN5ConstantCache{T, T2} <: OrdinaryDiffEqConstantCache c5::T2 c6::T2 c7::T2 + c8::T2 + c9::T2 a21::T a31::T a32::T @@ -128,6 +130,21 @@ struct FineRKN5ConstantCache{T, T2} <: OrdinaryDiffEqConstantCache a74::T a75::T #a76::T + a81::T + a82::T + a83::T + a84::T + a85::T + #a86::T + a87::T + a91::T + a92::T + a93::T + a94::T + a95::T + #a96::T + a97::T + #a98::T abar21::T abar31::T abar32::T @@ -149,6 +166,21 @@ struct FineRKN5ConstantCache{T, T2} <: OrdinaryDiffEqConstantCache abar74::T abar75::T abar76::T + abar81::T + abar82::T + abar83::T + abar84::T + abar85::T + #abar86::T + abar87::T + abar91::T + #abar92::T + abar93::T + abar94::T + abar95::T + #abar96::T + abar97::T + #abar98::T b1::T #b2::T b3::T @@ -177,6 +209,72 @@ struct FineRKN5ConstantCache{T, T2} <: OrdinaryDiffEqConstantCache bptilde5::T bptilde6::T bptilde7::T + r12::T + r13::T + r14::T + r15::T + r16::T + r33::T + r34::T + r35::T + r36::T + r43::T + r44::T + r45::T + r46::T + r53::T + r54::T + r55::T + r56::T + r63::T + r64::T + r65::T + r66::T + r73::T + r74::T + r75::T + r76::T + r83::T + r84::T + r85::T + r86::T + r93::T + r94::T + r95::T + r96::T + rp11::T + rp12::T + rp13::T + rp14::T + rp15::T + rp32::T + rp33::T + rp34::T + rp35::T + rp42::T + rp43::T + rp44::T + rp45::T + rp52::T + rp53::T + rp54::T + rp55::T + rp62::T + rp63::T + rp64::T + rp65::T + rp72::T + rp73::T + rp74::T + rp75::T + rp82::T + rp83::T + rp84::T + rp85::T + rp92::T + rp93::T + rp94::T + rp95::T end function FineRKN5ConstantCache(T::Type, T2::Type) @@ -187,6 +285,8 @@ function FineRKN5ConstantCache(T::Type, T2::Type) c5 = convert(T2, 43 // 47) c6 = convert(T2, 1 // 1) # 36463 // 36464 c7 = convert(T2, 1 // 1) + c8 = convert(T2, 2 // 5) + c9 = convert(T2, 1 // 5) a21 = convert(T, 32 // 1521) a31 = convert(T, 4 // 169) a32 = convert(T, 4 // 169) @@ -208,6 +308,21 @@ function FineRKN5ConstantCache(T::Type, T2::Type) a74 = convert(T, 3276 // 23575) a75 = convert(T, -1142053 // 22015140) #a76 = convert(T, 0 // 1) + a81 = convert(T, 3166724675977 // 89400687626250) + a82 = convert(T, 12182 // 175275) + a83 = convert(T, -2308196389073 // 59333908961250) + a84 = convert(T, 113223739712 // 2656609580625) + a85 = convert(T, -4985173058548 // 281912811350625) + #a86 = convert(T, 0 // 1) + a87 = convert(T, -108322 // 9884875) + a91 = convert(T, 13703589067379 // 1021293449694000) + a92 = convert(T, -16588 // 178955775) + a93 = convert(T, 393366167467741 // 44058124399590000) + a94 = convert(T, -129051960428 // 18967820851875) + a95 = convert(T, 286445641484101 // 88563993965842500) + #a96 = convert(T, 0 // 1) + a97 = convert(T, 185739 // 141153250) + #a98 = convert(T, 0 // 1) abar21 = convert(T, 8 // 39) abar31 = convert(T, 1 // 13) abar32 = convert(T, 3 // 13) @@ -229,6 +344,21 @@ function FineRKN5ConstantCache(T::Type, T2::Type) abar74 = convert(T, 19656 // 23575) abar75 = convert(T, -53676491 // 88060560) abar76 = convert(T, 53 // 240) + abar81 = convert(T, -153602563 // 3630543750) + abar82 = convert(T, 4 // 5) + abar83 = convert(T, -15809974379 // 36693821250) + abar84 = convert(T, 25076304 // 131584375) + abar85 = convert(T, -20756397983 // 250144464375) + #abar86 = convert(T, 0 // 1) + abar87 = convert(T, -251893 // 7317375) + abar91 = convert(T, 549292232911 // 4942380225000) + #abar92 = convert(T, 0 // 1) + abar93 = convert(T, 38673447228481 // 349667653965000) + abar94 = convert(T, -14447155986 // 237691990625) + abar95 = convert(T, 239970566676929 // 8434666091985000) + #abar96 = convert(T, 0 // 1) + abar97 = convert(T, 48691361 // 4597563000) + #abar98 = convert(T, 0 // 1) b1 = convert(T, 4817 // 51600) #b2 = convert(T, 0 // 1) b3 = convert(T, 388869 // 1216880) @@ -257,13 +387,86 @@ function FineRKN5ConstantCache(T::Type, T2::Type) bptilde5 = convert(T, -846261273 // 4494757750) bptilde6 = convert(T, 4187 // 36750) bptilde7 = convert(T, -1 // 25) - FineRKN5ConstantCache(c1, c2, c3, c4, c5, c6, c7, a21, a31, a32, a41, a43, a51, - a52, a53, a54, a61, a62, a63, a64, a71, a73, a74, a75, + r12 = convert(T, 1 // 2) + r13 = convert(T, -65425102193 // 34422618000) + r14 = convert(T, 246336178993 // 68845236000) + r15 = convert(T, -17271401477 // 5737103000) + r16 = convert(T, 2646330829 // 2868551500) + r33 = convert(T, -103408733716249 // 21918241774800) + r34 = convert(T, 36216248499769 // 2087451597600) + r35 = convert(T, -34348334365943 // 1826520147900) + r36 = convert(T, 11946681497647 // 1826520147900) + r43 = convert(T, 10008729576 // 15727000375) + r44 = convert(T, -45252884088 // 15727000375) + r45 = convert(T, 67501517184 // 15727000375) + r46 = convert(T, -28901603736 // 15727000375) + r53 = convert(T, -42869319978551 // 58745639878800) + r54 = convert(T, 26279956317109 // 8392234268400) + r55 = convert(T, -5270387298308 // 1223867497475) + r56 = convert(T, 103642379853661 // 58745639878800) + r63 = convert(T, 10265285443 // 27377989200) + r64 = convert(T, -14253109853 // 9125996400) + r65 = convert(T, 1166320544 // 570374775) + r66 = convert(T, -2474620297 // 3041998800) + r73 = convert(T, 201702341 // 608399760) + r74 = convert(T, -765268783 // 608399760) + r75 = convert(T, 861365839 // 608399760) + r76 = convert(T, -925430543 // 1825199280) + r83 = convert(T, 39325 // 31704) + r84 = convert(T, -96525 // 21136) + r85 = convert(T, 53625 // 10568) + r86 = convert(T, -7150 // 3963) + r93 = convert(T, 25525 // 4848) + r94 = convert(T, -25525 // 1616) + r95 = convert(T, 25525 // 1616) + r96 = convert(T, -25525 // 4848) + rp11 = convert(T, 1) + rp12 = convert(T, -65425102193 // 11474206000) + rp13 = convert(T, 246336178993 // 17211309000) + rp14 = convert(T, -17271401477 // 1147420600) + rp15 = convert(T, 7938992487 // 1434275750) + rp32 = convert(T, -103408733716249 // 7306080591600) + rp33 = convert(T, 36216248499769 // 521862899400) + rp34 = convert(T, -34348334365943 // 365304029580) + rp35 = convert(T, 11946681497647 // 304420024650) + rp42 = convert(T, 30026188728 // 15727000375) + rp43 = convert(T, -181011536352 // 15727000375) + rp44 = convert(T, 67501517184 // 3145400075) + rp45 = convert(T, -173409622416 // 15727000375) + rp52 = convert(T, -42869319978551 // 19581879959600) + rp53 = convert(T, 26279956317109 // 2098058567100) + rp54 = convert(T, -5270387298308 // 244773499495) + rp55 = convert(T, 103642379853661 // 9790939979800) + rp62 = convert(T, 10265285443 // 9125996400) + rp63 = convert(T, -14253109853 // 2281499100) + rp64 = convert(T, 1166320544 // 114074955) + rp65 = convert(T, -2474620297 // 506999800) + rp72 = convert(T, 201702341 // 202799920) + rp73 = convert(T, -765268783 // 152099940) + rp74 = convert(T, 861365839 // 121679952) + rp75 = convert(T, -925430543 // 304199880) + rp82 = convert(T, 39325 // 10568) + rp83 = convert(T, -96525 // 5284) + rp84 = convert(T, 268125 // 10568) + rp85 = convert(T, -14300 // 1321) + rp92 = convert(T, 25525 // 1616) + rp93 = convert(T, -25525 // 404) + rp94 = convert(T, 127625 // 1616) + rp95 = convert(T, -25525 // 808) + FineRKN5ConstantCache(c1, c2, c3, c4, c5, c6, c7, c8, c9, a21, a31, a32, a41, a43, a51, + a52, a53, a54, a61, a62, a63, a64, a71, a73, a74, a75, a81, a82, a83, a84, a85, a87, + a91, a92, a93, a94, a95, a97, abar21, abar31, abar32, abar41, abar42, abar43, abar51, abar52, abar53, abar54, abar61, abar62, abar63, abar64, abar65, - abar71, abar73, abar74, abar75, abar76, b1, b3, b4, + abar71, abar73, abar74, abar75, abar76, abar81, abar82, abar83, abar84, abar85, + abar87, abar91, abar93, abar94, abar95, abar97, b1, b3, b4, b5, bbar1, bbar3, bbar4, bbar5, bbar6, btilde1, btilde3, btilde4, btilde5, bptilde1, - bptilde3, bptilde4, bptilde5, bptilde6, bptilde7) + bptilde3, bptilde4, bptilde5, bptilde6, bptilde7, r12, r13, r14, r15, r16, r33, r34, + r35, r36, r43, + r44, r45, r46, r53, r54, r55, r56, r63, r64, r65, r66, r73, r74, r75, r76, r83, r84, + r85, r86, r93, r94, r95, r96, rp11, rp12, rp13, rp14, rp15, rp32, rp33, rp34, rp35, + rp42, rp43, rp44, rp45, rp52, rp53, rp54, rp55, rp62, rp63, rp64, rp65, rp72, rp73, + rp74, rp75, rp82, rp83, rp84, rp85, rp92, rp93, rp94, rp95) end struct IRKN3ConstantCache{T, T2} <: OrdinaryDiffEqConstantCache From 4f08178b826867ea3e736a9ed04c734e730ba811 Mon Sep 17 00:00:00 2001 From: Henry Langner <128313572+HenryLangner@users.noreply.github.com> Date: Fri, 10 Nov 2023 11:10:00 +0000 Subject: [PATCH 2/2] Changed the form of the specialized Interpolent from RK-Step to Hermit. --- src/algorithms.jl | 1 + src/caches/rkn_caches.jl | 33 +++- src/dense/interpolants.jl | 136 ++++--------- src/perform_step/rkn_perform_step.jl | 24 ++- src/tableaus/rkn_tableaus.jl | 281 ++++++++++++++++++--------- 5 files changed, 275 insertions(+), 200 deletions(-) diff --git a/src/algorithms.jl b/src/algorithms.jl index 0b588e8271..ddadd1cb88 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -5278,6 +5278,7 @@ struct FineRKN4 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end A 5th order explicit Runge-Kutta-Nyström method which can be applied directly to second order ODEs. In particular, this method allows the acceleration equation to depend on the velocity. +Algorithm uses a specialized 6th order interpolant. ## References ``` diff --git a/src/caches/rkn_caches.jl b/src/caches/rkn_caches.jl index 301dfb8060..b8779ac79f 100644 --- a/src/caches/rkn_caches.jl +++ b/src/caches/rkn_caches.jl @@ -78,6 +78,28 @@ function alg_cache(alg::FineRKN4, u, rate_prototype, ::Type{uEltypeNoUnits}, FineRKN4ConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits)) end +@cache struct FineRKN5Cache{uType, rateType, reducedRateType, uNoUnitsType, TabType} <: + OrdinaryDiffEqMutableCache + u::uType + uprev::uType + fsalfirst::rateType + k2::reducedRateType + k3::reducedRateType + k4::reducedRateType + k5::reducedRateType + k6::reducedRateType + k7::reducedRateType + k8::reducedRateType + k9::reducedRateType + k10::reducedRateType + k11::reducedRateType + k::rateType + utilde::uType + tmp::uType + atmp::uNoUnitsType + tab::TabType +end + function alg_cache(alg::FineRKN5, u, rate_prototype, ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, dt, reltol, p, calck, @@ -93,19 +115,14 @@ function alg_cache(alg::FineRKN5, u, rate_prototype, ::Type{uEltypeNoUnits}, k7 = zero(reduced_rate_prototype) k8 = zero(reduced_rate_prototype) k9 = zero(reduced_rate_prototype) + k10 = zero(reduced_rate_prototype) + k11 = zero(reduced_rate_prototype) k = zero(rate_prototype) utilde = zero(u) atmp = similar(u, uEltypeNoUnits) recursivefill!(atmp, false) tmp = zero(u) - FineRKN5Cache(u, uprev, k1, k2, k3, k4, k5, k6, k7, k8, k9, k, utilde, tmp, atmp, tab) -end - -function alg_cache(alg::FineRKN5, u, rate_prototype, ::Type{uEltypeNoUnits}, - ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, - dt, reltol, p, calck, - ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} - FineRKN5ConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits)) + FineRKN5Cache(u, uprev, k1, k2, k3, k4, k5, k6, k7, k8, k9, k10, k11, k, utilde, tmp, atmp, tab) end function alg_cache(alg::FineRKN5, u, rate_prototype, ::Type{uEltypeNoUnits}, diff --git a/src/dense/interpolants.jl b/src/dense/interpolants.jl index 8e9ee8a747..b6b8b90958 100644 --- a/src/dense/interpolants.jl +++ b/src/dense/interpolants.jl @@ -2847,49 +2847,52 @@ end """ @def finerkn5unpack begin if typeof(cache) <: OrdinaryDiffEqMutableCache - @unpack r12, r13, r14, r15, r16, r33, r34, - r35, r36, r43, - r44, r45, r46, r53, r54, r55, r56, r63, r64, r65, r66, r73, r74, r75, r76, r83, r84, - r85, r86, r93, r94, r95, r96, rp11, rp12, rp13, rp14, rp15, rp32, rp33, rp34, rp35, - rp42, rp43, rp44, rp45, rp52, rp53, rp54, rp55, rp62, rp63, rp64, rp65, rp72, rp73, - rp74, rp75, rp82, rp83, rp84, rp85, rp92, rp93, rp94, rp95 = cache.tab + @unpack r10, r11, + r12, r13, r14, r15, r16, r20, r21, r22, r23, r24, r25, r26, r30, r31, r32, r33, r34, + r35, r36, r40, r41, r42, r43, r44, r45, r46, r50, r51, r52, r53, r54, r55, r56, r60, + r61, r62, r63, r64, r65, r66, r70, r71, r72, r73, r74, r75, r76, rp10, rp11, rp12, + rp13, rp14, rp15, rp16, rp20, rp21, rp22, rp23, rp24, rp25, rp26, rp30, rp31, rp32, + rp33, rp34, rp35, rp36, rp40, rp41, rp42, rp43, rp44, rp45, rp46, rp50, rp51, rp52, + rp53, rp54, rp55, rp56, rp60, rp61, rp62, rp63, rp64, rp65, rp66, rp70, rp71, rp72, + rp73, rp74, rp75, rp76 = cache.tab else - @unpack r12, r13, r14, r15, r16, r33, r34, - r35, r36, r43, - r44, r45, r46, r53, r54, r55, r56, r63, r64, r65, r66, r73, r74, r75, r76, r83, r84, - r85, r86, r93, r94, r95, r96, rp11, rp12, rp13, rp14, rp15, rp32, rp33, rp34, rp35, - rp42, rp43, rp44, rp45, rp52, rp53, rp54, rp55, rp62, rp63, rp64, rp65, rp72, rp73, - rp74, rp75, rp82, rp83, rp84, rp85, rp92, rp93, rp94, rp95 = cache + @unpack r10, r11, r12, r13, r14, r15, r16, r20, r21, r22, r23, r24, r25, r26, r30, r31, r32, r33, r34, + r35, r36, r40, r41, r42, r43, r44, r45, r46, r50, r51, r52, r53, r54, r55, r56, r60, + r61, r62, r63, r64, r65, r66, r70, r71, r72, r73, r74, r75, r76, rp10, rp11, rp12, + rp13, rp14, rp15, rp16, rp20, rp21, rp22, rp23, rp24, rp25, rp26, rp30, rp31, rp32, + rp33, rp34, rp35, rp36, rp40, rp41, rp42, rp43, rp44, rp45, rp46, rp50, rp51, rp52, + rp53, rp54, rp55, rp56, rp60, rp61, rp62, rp63, rp64, rp65, rp66, rp70, rp71, rp72, + rp73, rp74, + rp75, rp76 = cache end end @def finerkn5pre0 begin @finerkn5unpack - b1Θ = Θ^2 * @evalpoly(Θ, r12, r13, r14, r15, r16) - b3Θ = Θ^3 * @evalpoly(Θ, r33, r34, r35, r36) - b4Θ = Θ^3 * @evalpoly(Θ, r43, r44, r45, r46) - b5Θ = Θ^3 * @evalpoly(Θ, r53, r54, r55, r56) - b6Θ = Θ^3 * @evalpoly(Θ, r63, r64, r65, r66) - b7Θ = Θ^3 * @evalpoly(Θ, r73, r74, r75, r76) - b8Θ = Θ^3 * @evalpoly(Θ, r83, r84, r85, r86) - b9Θ = Θ^3 * @evalpoly(Θ, r93, r94, r95, r96) - - bp1Θ = Θ * @evalpoly(Θ, rp11, rp12, rp13, rp14, rp15) - bp3Θ = Θ^2 * @evalpoly(Θ, rp32, rp33, rp34, rp35) - bp4Θ = Θ^2 * @evalpoly(Θ, rp42, rp43, rp44, rp45) - bp5Θ = Θ^2 * @evalpoly(Θ, rp52, rp53, rp54, rp55) - bp6Θ = Θ^2 * @evalpoly(Θ, rp62, rp63, rp64, rp65) - bp7Θ = Θ^2 * @evalpoly(Θ, rp72, rp73, rp74, rp75) - bp8Θ = Θ^2 * @evalpoly(Θ, rp82, rp83, rp84, rp85) - bp9Θ = Θ^2 * @evalpoly(Θ, rp92, rp93, rp94, rp95) + b1Θ = @evalpoly(Θ, r10, r11, r12, r13, r14, r15, r16) + b2Θ = @evalpoly(Θ, r20, r21, r22, r23, r24, r25, r26) + b3Θ = @evalpoly(Θ, r30, r31, r32, r33, r34, r35, r36) + b4Θ = @evalpoly(Θ, r40, r41, r42, r43, r44, r45, r46) + b5Θ = @evalpoly(Θ, r50, r51, r52, r53, r54, r55, r56) + b6Θ = @evalpoly(Θ, r60, r61, r62, r63, r64, r65, r66) + b7Θ = @evalpoly(Θ, r70, r71, r72, r73, r74, r75, r76) + + bp1Θ = @evalpoly(Θ, rp10, rp11, rp12, rp13, rp14, rp15) + bp2Θ = @evalpoly(Θ, rp20, rp21, rp22, rp23, rp24, rp25) + bp3Θ = @evalpoly(Θ, rp30, rp31, rp32, rp33, rp34, rp35) + bp4Θ = @evalpoly(Θ, rp40, rp41, rp42, rp43, rp44, rp45) + bp5Θ = @evalpoly(Θ, rp50, rp51, rp52, rp53, rp54, rp55) + bp6Θ = @evalpoly(Θ, rp60, rp61, rp62, rp63, rp64, rp65) + bp7Θ = @evalpoly(Θ, rp70, rp71, rp72, rp73, rp74, rp75) kk1, kk2, kk3, kk4 = k k1, k3 = kk1.x k4, k5 = kk2.x k6, k7 = kk3.x - k8, k9 = kk4.x + u1, u2 = kk4.x duprev, uprev = y₀.x + dunew, unew = y₁.x dtsq = dt^2 end @@ -2897,62 +2900,23 @@ end cache::Union{FineRKN5ConstantCache, FineRKN5Cache}, idxs::Nothing, T::Type{Val{0}}) @finerkn5pre0 - return ArrayPartition(duprev + - dt * Θ * - (bp1Θ * k1 + bp3Θ * k3 + - bp4Θ * k4 + bp5Θ * k5 + bp6Θ * k6 + bp7Θ * k7 + bp8Θ * k8 + - bp9Θ * k9), - uprev + - dt * Θ * - (duprev + - dt * Θ * - (b1Θ * k1 + b3Θ * k3 + - b4Θ * k4 + b5Θ * k5 + b6Θ * k6 + b7Θ * k7 + b8Θ * k8 + b9Θ * k9))) + return ArrayPartition((bp1Θ*uprev)/dt + bp2Θ*duprev + bp3Θ*dt*k1 + (bp4Θ/dt)*u2 + bp5Θ*dt*k7 + bp6Θ*dunew + (bp7Θ*u1)/dt, b1Θ*uprev + b2Θ*dt*duprev + b3Θ*dtsq*k1 + b4Θ*u2 + b5Θ*dtsq*k7 + b6Θ*dt*dunew + b7Θ*u1) end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{FineRKN5ConstantCache, FineRKN5Cache}, idxs, T::Type{Val{0}}) @finerkn5pre0 - return ArrayPartition(duprev[idxs] + - dt * Θ * - (bp1Θ * k1[idxs] + bp3Θ * k3[idxs] + - bp4Θ * k4[idxs] + bp5Θ * k5[idxs] + bp6Θ * k6[idxs] + - bp7Θ * k7[idxs] + bp8Θ * k8[idxs] + bp9Θ * k9[idxs]), - uprev[idxs] + - dt * Θ * - (duprev[idxs] + - dt * Θ * - (b1Θ * k1[idxs] + - b3Θ * k3[idxs] + - b4Θ * k4[idxs] + b5Θ * k5[idxs] + b6Θ * k6[idxs] + b7Θ * k7[idxs] + - b8Θ * k8[idxs] + b9Θ * k9[idxs]))) + return ArrayPartition((bp1Θ*uprev[idxs])/dt + bp2Θ*duprev[idxs] + bp3Θ*dt*k1[idxs] + (bp4Θ/dt)*u2[idxs] + bp5Θ*dt*k7[idxs] + bp6Θ*dunew[idxs] + (bp7Θ*u1)/dt[idxs], + b1Θ*uprev[idxs] + b2Θ*dt*duprev[idxs] + b3Θ*dtsq*k1[idxs] + b4Θ*u2[idxs] + b5Θ*dtsq*k7[idxs] + b6Θ*dt*dunew[idxs] + b7Θ*u1[idxs]) end @muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, cache::Union{FineRKN5ConstantCache, FineRKN5Cache}, idxs::Nothing, T::Type{Val{0}}) @finerkn5pre0 - @inbounds @.. broadcast=false out.x[2]=uprev + - dt * Θ * - (duprev + - dt * Θ * - (b1Θ * k1 + - b3Θ * k3 + - b4Θ * k4 + b5Θ * k5 + b6Θ * k6 + b7Θ * k7 + - b8Θ * k8 + b9Θ * k9)) - @inbounds @.. broadcast=false out.x[1]=duprev + - dt * Θ * - (bp1Θ * k1 + bp3Θ * k3 + - bp4Θ * k4 + bp5Θ * k5 + bp6Θ * k6 + bp7Θ * k7 + - bp8Θ * k8 + bp9Θ * k9) - #for i in eachindex(out.x[1]) - # out.x[2][i] = uprev[i] + dt*Θ*(duprev[i] + dt*Θ*(b1Θ*k1[i] + - # b3Θ*k3[i] + - # b4Θ*k4[i] + b5Θ*k5[i] + b6Θ*k6[i])) - # out.x[1][i] = duprev[i] + dt*Θ*(bp1Θ*k1[i] + bp3Θ*k3[i] + - # bp4Θ*k4[i] + bp5Θ*k5[i] + bp6Θ*k6[i]) - #end + @inbounds @.. broadcast=false out.x[1]=(bp1Θ*uprev)/dt + bp2Θ*duprev + bp3Θ*dt*k1 + (bp4Θ/dt)*u2 + bp5Θ*dt*k7 + bp6Θ*dunew + (bp7Θ*u1)/dt + @inbounds @.. broadcast=false out.x[2]=b1Θ*uprev + b2Θ*dt*duprev + b3Θ*dtsq*k1 + b4Θ*u2 + b5Θ*dtsq*k7 + b6Θ*dt*dunew + b7Θ*u1 out end @@ -2960,28 +2924,8 @@ end cache::Union{FineRKN5ConstantCache, FineRKN5Cache}, idxs, T::Type{Val{0}}) @finerkn5pre0 - @inbounds @.. broadcast=false out.x[2]=uprev[idxs] + - dt * Θ * - (duprev[idxs] + - dt * Θ * - (b1Θ * k1[idxs] + - b3Θ * k3[idxs] + - b4Θ * k4[idxs] + b5Θ * k5[idxs] + - b6Θ * k6[idxs] + b7Θ * k7[idxs] + - b8Θ * k8[idxs] + b9Θ * k9[idxs])) - @inbounds @.. broadcast=false out.x[1]=duprev[idxs] + - dt * Θ * - (bp1Θ * k1[idxs] + bp3Θ * k3[idxs] + - bp4Θ * k4[idxs] + bp5Θ * k5[idxs] + - bp6Θ * k6[idxs] + bp7Θ * k7[idxs] + - bp8Θ * k8[idxs] + bp9Θ * k9[idxs]) - #for (j,i) in enumerate(idxs) - # out.x[2][j] = uprev[i] + dt*Θ*(duprev[i] + dt*Θ*(b1Θ*k1[i] + - # b3Θ*k3[i] + - # b4Θ*k4[i] + b5Θ*k5[i] + b6Θ*k6[i])) - # out.x[1][j] = duprev[i] + dt*Θ*(bp1Θ*k1[i] + bp3Θ*k3[i] + - # bp4Θ*k4[i] + bp5Θ*k5[i] + bp6Θ*k6[i]) - #end + @inbounds @.. broadcast=false out.x[1]=b1Θ*uprev[idxs] + b2Θ*dt*duprev[idxs] + b3Θ*dtsq*k1[idxs] + b4Θ*u2[idxs] + b5Θ*dtsq*k7[idxs] + b6Θ*dt*dunew[idxs] + b7Θ*u1[idxs] + @inbounds @.. broadcast=false out.x[2]=(bp1Θ*uprev[idxs])/dt + bp2Θ*duprev[idxs] + bp3Θ*dt*k1[idxs] + (bp4Θ/dt)*u2[idxs] + bp5Θ*dt*k7[idxs] + bp6Θ*dunew[idxs] + (bp7Θ*u1)/dt[idxs] out end """ diff --git a/src/perform_step/rkn_perform_step.jl b/src/perform_step/rkn_perform_step.jl index ffe5b56e82..3b7715642a 100644 --- a/src/perform_step/rkn_perform_step.jl +++ b/src/perform_step/rkn_perform_step.jl @@ -260,7 +260,8 @@ end abar52, abar53, abar54, abar61, abar62, abar63, abar64, abar65, abar71, abar73, abar74, abar75, abar76, abar81, abar82, abar83, abar84, abar85, abar87, abar91, abar93, abar94, abar95, abar97, b1, b3, b4, b5, bbar1, bbar3, bbar4, bbar5, bbar6, btilde1, btilde3, btilde4, btilde5, bptilde1, - bptilde3, bptilde4, bptilde5, bptilde6, bptilde7 = cache + bptilde3, bptilde4, bptilde5, bptilde6, bptilde7, btilde1_1, + btilde2_1, btilde3_1, btilde4_1, btilde5_1, btilde6_1, btilde7_1, btilde8_1, btilde9_1, btilde1_12, btilde2_12, btilde3_12, btilde4_12, btilde5_12, btilde6_12, btilde7_12, btilde8_12, btilde9_12 = cache k1 = integrator.fsalfirst.x[1] ku = uprev + dt * (c2 * duprev + dt * (a21 * k1)) @@ -312,15 +313,17 @@ end abar97 * k7) # abar92 = abar96 = abar98 = 0 k9 = f.f1(kdu, ku, p, t + dt * c9) - + k10 = uprev + dt * (duprev + dt * (btilde1_1 * k1 + btilde2_1 * k2 + btilde3_1 * k3 + btilde4_1 * k4 + btilde5_1 * k5 + btilde6_1 * k6 + btilde7_1 * k7 + btilde8_1 * k8 + btilde9_1 * k9)) # higher order Evaluations for ytilde_{n+1} (as donated in corresponding Paper) + k11 = uprev + dt * (duprev * 0.5 + dt * (btilde1_12 * k1 + btilde2_12 * k2 + btilde3_12 * k3 + btilde4_12 * k4 + btilde5_12 * k5 + btilde6_12 * k6 + btilde7_12 * k7 + btilde8_12 * k8 + btilde9_12 * k9)) # higher order Evaluations for ytilde_{n+1/2} (as donated in corresponding Paper) + integrator.u = ArrayPartition((du, u)) integrator.fsallast = ArrayPartition((f.f1(du, u, p, t + dt), f.f2(du, u, p, t + dt))) - integrator.stats.nf += 9 + integrator.stats.nf += 11 integrator.stats.nf2 += 1 integrator.k[1] = ArrayPartition(integrator.fsalfirst.x[1], k3) integrator.k[2] = ArrayPartition(k4, k5) integrator.k[3] = ArrayPartition(k6, k7) - integrator.k[4] = ArrayPartition(k8, k9) + integrator.k[4] = ArrayPartition(k10, k11) if integrator.opts.adaptive dtsq = dt^2 @@ -346,7 +349,7 @@ function initialize!(integrator, cache::FineRKN5Cache) integrator.k[1] = ArrayPartition(cache.fsalfirst.x[1], cache.k3) integrator.k[2] = ArrayPartition(cache.k4, cache.k5) integrator.k[3] = ArrayPartition(cache.k6, cache.k7) - integrator.k[4] = ArrayPartition(cache.k8, cache.k9) + integrator.k[4] = ArrayPartition(cache.k10, cache.k11) integrator.f.f1(integrator.fsallast.x[1], duprev, uprev, integrator.p, integrator.t) integrator.f.f2(integrator.fsallast.x[2], duprev, uprev, integrator.p, integrator.t) integrator.stats.nf += 1 @@ -357,14 +360,15 @@ end @unpack t, dt, f, p = integrator du, u = integrator.u.x duprev, uprev = integrator.uprev.x - @unpack tmp, atmp, fsalfirst, k2, k3, k4, k5, k6, k7, k8, k9, k, utilde = cache + @unpack tmp, atmp, fsalfirst, k2, k3, k4, k5, k6, k7, k8, k9, k10, k11, k, utilde = cache @unpack c2, c3, c4, c5, c6, c7, c8, c9, a21, a31, a32, a41, a43, a51, a52, a53, a54, a61, a62, a63, a64, a71, a73, a74, a75, a81, a82, a83, a84, a85, a87, a91, a92, a93, a94, a95, a97, abar21, abar31, abar32, abar41, abar42, abar43, abar51, abar52, abar53, abar54, abar61, abar62, abar63, abar64, abar65, abar71, abar73, abar74, abar75, abar76, abar81, abar82, abar83, abar84, abar85, abar87, abar91, abar93, abar94, abar95, abar97, b1, b3, b4, b5, bbar1, bbar3, bbar4, bbar5, bbar6, btilde1, btilde3, btilde4, btilde5, bptilde1, - bptilde3, bptilde4, bptilde5, bptilde6, bptilde7 = cache.tab + bptilde3, bptilde4, bptilde5, bptilde6, bptilde7, btilde1_1, + btilde2_1, btilde3_1, btilde4_1, btilde5_1, btilde6_1, btilde7_1, btilde8_1, btilde9_1, btilde1_12, btilde2_12, btilde3_12, btilde4_12, btilde5_12, btilde6_12, btilde7_12, btilde8_12, btilde9_12 = cache.tab kdu, ku = integrator.cache.tmp.x[1], integrator.cache.tmp.x[2] uidx = eachindex(integrator.uprev.x[2]) k1 = integrator.fsalfirst.x[1] @@ -433,7 +437,11 @@ end f.f1(k.x[1], du, u, p, t + dt) f.f2(k.x[2], du, u, p, t + dt) - integrator.stats.nf += 9 + + @.. broadcast=false k10= uprev + dt * (duprev + dt * (btilde1_1 * k1 + btilde2_1 * k2 + btilde3_1 * k3 + btilde4_1 * k4 + btilde5_1 * k5 + btilde6_1 * k6 + btilde7_1 * k7 + btilde8_1 * k8 + btilde9_1 * k9)) # higher order Evaluations for ytilde_{n+1} (as donated in corresponding Paper) + @.. broadcast=false k11= uprev + dt * (duprev * 0.5 + dt * (btilde1_12 * k1 + btilde2_12 * k2 + btilde3_12 * k3 + btilde4_12 * k4 + btilde5_12 * k5 + btilde6_12 * k6 + btilde7_12 * k7 + btilde8_12 * k8 + btilde9_12 * k9)) # higher order Evaluations for ytilde_{n+1/2} (as donated in corresponding Paper) + + integrator.stats.nf += 11 integrator.stats.nf2 += 1 if integrator.opts.adaptive duhat, uhat = utilde.x diff --git a/src/tableaus/rkn_tableaus.jl b/src/tableaus/rkn_tableaus.jl index 3463424e65..de4cc02005 100644 --- a/src/tableaus/rkn_tableaus.jl +++ b/src/tableaus/rkn_tableaus.jl @@ -209,72 +209,122 @@ struct FineRKN5ConstantCache{T, T2} <: OrdinaryDiffEqConstantCache bptilde5::T bptilde6::T bptilde7::T + btilde1_1::T + btilde2_1::T + btilde3_1::T + btilde4_1::T + btilde5_1::T + btilde6_1::T + btilde7_1::T + btilde8_1::T + btilde9_1::T + btilde1_12::T + btilde2_12::T + btilde3_12::T + btilde4_12::T + btilde5_12::T + btilde6_12::T + btilde7_12::T + btilde8_12::T + btilde9_12::T + r10::T + r11::T r12::T r13::T r14::T r15::T r16::T + r20::T + r21::T + r22::T + r23::T + r24::T + r25::T + r26::T + r30::T + r31::T + r32::T r33::T r34::T r35::T r36::T + r40::T + r41::T + r42::T r43::T r44::T r45::T r46::T + r50::T + r51::T + r52::T r53::T r54::T r55::T r56::T + r60::T + r61::T + r62::T r63::T r64::T r65::T r66::T + r70::T + r71::T + r72::T r73::T r74::T r75::T r76::T - r83::T - r84::T - r85::T - r86::T - r93::T - r94::T - r95::T - r96::T + rp10::T rp11::T rp12::T rp13::T rp14::T rp15::T + rp16::T + rp20::T + rp21::T + rp22::T + rp23::T + rp24::T + rp25::T + rp26::T + rp30::T + rp31::T rp32::T rp33::T rp34::T rp35::T + rp36::T + rp40::T + rp41::T rp42::T rp43::T rp44::T rp45::T + rp46::T + rp50::T + rp51::T rp52::T rp53::T rp54::T rp55::T + rp56::T + rp60::T + rp61::T rp62::T rp63::T rp64::T rp65::T + rp66::T + rp70::T + rp71::T rp72::T rp73::T rp74::T rp75::T - rp82::T - rp83::T - rp84::T - rp85::T - rp92::T - rp93::T - rp94::T - rp95::T + rp76::T end function FineRKN5ConstantCache(T::Type, T2::Type) @@ -387,72 +437,122 @@ function FineRKN5ConstantCache(T::Type, T2::Type) bptilde5 = convert(T, -846261273 // 4494757750) bptilde6 = convert(T, 4187 // 36750) bptilde7 = convert(T, -1 // 25) - r12 = convert(T, 1 // 2) - r13 = convert(T, -65425102193 // 34422618000) - r14 = convert(T, 246336178993 // 68845236000) - r15 = convert(T, -17271401477 // 5737103000) - r16 = convert(T, 2646330829 // 2868551500) - r33 = convert(T, -103408733716249 // 21918241774800) - r34 = convert(T, 36216248499769 // 2087451597600) - r35 = convert(T, -34348334365943 // 1826520147900) - r36 = convert(T, 11946681497647 // 1826520147900) - r43 = convert(T, 10008729576 // 15727000375) - r44 = convert(T, -45252884088 // 15727000375) - r45 = convert(T, 67501517184 // 15727000375) - r46 = convert(T, -28901603736 // 15727000375) - r53 = convert(T, -42869319978551 // 58745639878800) - r54 = convert(T, 26279956317109 // 8392234268400) - r55 = convert(T, -5270387298308 // 1223867497475) - r56 = convert(T, 103642379853661 // 58745639878800) - r63 = convert(T, 10265285443 // 27377989200) - r64 = convert(T, -14253109853 // 9125996400) - r65 = convert(T, 1166320544 // 570374775) - r66 = convert(T, -2474620297 // 3041998800) - r73 = convert(T, 201702341 // 608399760) - r74 = convert(T, -765268783 // 608399760) - r75 = convert(T, 861365839 // 608399760) - r76 = convert(T, -925430543 // 1825199280) - r83 = convert(T, 39325 // 31704) - r84 = convert(T, -96525 // 21136) - r85 = convert(T, 53625 // 10568) - r86 = convert(T, -7150 // 3963) - r93 = convert(T, 25525 // 4848) - r94 = convert(T, -25525 // 1616) - r95 = convert(T, 25525 // 1616) - r96 = convert(T, -25525 // 4848) - rp11 = convert(T, 1) - rp12 = convert(T, -65425102193 // 11474206000) - rp13 = convert(T, 246336178993 // 17211309000) - rp14 = convert(T, -17271401477 // 1147420600) - rp15 = convert(T, 7938992487 // 1434275750) - rp32 = convert(T, -103408733716249 // 7306080591600) - rp33 = convert(T, 36216248499769 // 521862899400) - rp34 = convert(T, -34348334365943 // 365304029580) - rp35 = convert(T, 11946681497647 // 304420024650) - rp42 = convert(T, 30026188728 // 15727000375) - rp43 = convert(T, -181011536352 // 15727000375) - rp44 = convert(T, 67501517184 // 3145400075) - rp45 = convert(T, -173409622416 // 15727000375) - rp52 = convert(T, -42869319978551 // 19581879959600) - rp53 = convert(T, 26279956317109 // 2098058567100) - rp54 = convert(T, -5270387298308 // 244773499495) - rp55 = convert(T, 103642379853661 // 9790939979800) - rp62 = convert(T, 10265285443 // 9125996400) - rp63 = convert(T, -14253109853 // 2281499100) - rp64 = convert(T, 1166320544 // 114074955) - rp65 = convert(T, -2474620297 // 506999800) - rp72 = convert(T, 201702341 // 202799920) - rp73 = convert(T, -765268783 // 152099940) - rp74 = convert(T, 861365839 // 121679952) - rp75 = convert(T, -925430543 // 304199880) - rp82 = convert(T, 39325 // 10568) - rp83 = convert(T, -96525 // 5284) - rp84 = convert(T, 268125 // 10568) - rp85 = convert(T, -14300 // 1321) - rp92 = convert(T, 25525 // 1616) - rp93 = convert(T, -25525 // 404) - rp94 = convert(T, 127625 // 1616) - rp95 = convert(T, -25525 // 808) + btilde1_1 = convert(T, 20342293 // 227212000) + btilde2_1 = convert(T, 0 // 1) + btilde3_1 = convert(T, 159248338847 // 434024589600) + btilde4_1 = convert(T, 33225336 // 155712875) + btilde5_1 = convert(T, -27213980937 // 193879999600) + btilde6_1 = convert(T, 12057023 // 271069200) + btilde7_1 = convert(T, -19822 // 1129455) + btilde8_1 = convert(T, -3575 // 63408) + btilde9_1 = convert(T, 0 // 1) + btilde1_12 = convert(T, 26173973 // 833856000) + btilde2_12 = convert(T, 0 // 1) + btilde3_12 = convert(T, 692804177 // 75849868800) + btilde4_12 = convert(T, 488241 // 95243000) + btilde5_12 = convert(T, -1019853329 // 406588185600) + btilde6_12 = convert(T, 590579 // 1326412800) + btilde7_12 = convert(T, -75401 // 88427520) + btilde8_12 = convert(T, 0 // 1) + btilde9_12 = convert(T, 25525 // 310272) + r10 = convert(T, 1 // 1) + r11 = convert(T, 0 // 1) + r12 = convert(T, 0 // 1) + r13 = convert(T, -42 // 1) + r14 = convert(T, 111 // 1) + r15 = convert(T, -102 // 1) + r16 = convert(T, 32 // 1) + r20 = convert(T, 0 // 1) + r21 = convert(T, 1 // 1) + r22 = convert(T, 0 // 1) + r23 = convert(T, -16 // 1) + r24 = convert(T, 38 // 1) + r25 = convert(T, -33 // 1) + r26 = convert(T, 10 // 1) + r30 = convert(T, 0 // 2) + r31 = convert(T, 0 // 2) + r32 = convert(T, 1 // 2) + r33 = convert(T, -5 // 2) + r34 = convert(T, 9 // 2) + r35 = convert(T, -7 // 2) + r36 = convert(T, 1 // 1) + r40 = convert(T, 0 // 1) + r41 = convert(T, 0 // 1) + r42 = convert(T, 0 // 1) + r43 = convert(T, 64 // 1) + r44 = convert(T, -192 // 1) + r45 = convert(T, 192 // 1) + r46 = convert(T, -64 // 1) + r50 = convert(T, 0 // 1) + r51 = convert(T, 0 // 1) + r52 = convert(T, 0 // 1) + r53 = convert(T, -1 // 2) + r54 = convert(T, 2 // 1) + r55 = convert(T, -5 // 2) + r56 = convert(T, 1 // 1) + r60 = convert(T, 0 // 1) + r61 = convert(T, 0 // 1) + r62 = convert(T, 0 // 1) + r63 = convert(T, 6 // 1) + r64 = convert(T, -23 // 1) + r65 = convert(T, 27 // 1) + r66 = convert(T, -10 // 1) + r70 = convert(T, 0 // 1) + r71 = convert(T, 0 // 1) + r72 = convert(T, 0 // 1) + r73 = convert(T, -22 // 1) + r74 = convert(T, 81 // 1) + r75 = convert(T, -90 // 1) + r76 = convert(T, 32 // 1) + rp10 = convert(T, 0 // 1) + rp11 = convert(T, 0 // 1) + rp12 = convert(T, -126 // 1) + rp13 = convert(T, 444 // 1) + rp14 = convert(T, -510 // 1) + rp15 = convert(T, 192 // 1) + rp16 = convert(T, 0 // 1) + rp20 = convert(T, 1 // 1) + rp21 = convert(T, 0 // 1) + rp22 = convert(T, -48 // 1) + rp23 = convert(T, 152 // 1) + rp24 = convert(T, -165 // 1) + rp25 = convert(T, 60 // 1) + rp26 = convert(T, 0 // 1) + rp30 = convert(T, 0 // 1) + rp31 = convert(T, 1 // 1) + rp32 = convert(T, -15 // 2) + rp33 = convert(T, 18 // 1) + rp34 = convert(T, -35 // 2) + rp35 = convert(T, 6 // 1) + rp36 = convert(T, 0 // 1) + rp40 = convert(T, 0 // 1) + rp41 = convert(T, 0 // 1) + rp42 = convert(T, 192 // 1) + rp43 = convert(T, -768 // 1) + rp44 = convert(T, 960 // 1) + rp45 = convert(T, -384 // 1) + rp46 = convert(T, 0 // 1) + rp50 = convert(T, 0 // 1) + rp51 = convert(T, 0 // 1) + rp52 = convert(T, -3 // 2) + rp53 = convert(T, 8 // 1) + rp54 = convert(T, -25 // 2) + rp55 = convert(T, 6 // 1) + rp56 = convert(T, 0 // 1) + rp60 = convert(T, 0 // 1) + rp61 = convert(T, 0 // 1) + rp62 = convert(T, 18 // 1) + rp63 = convert(T, -92 // 1) + rp64 = convert(T, 135 // 1) + rp65 = convert(T, -60 // 1) + rp66 = convert(T, 0 // 1) + rp70 = convert(T, 0 // 1) + rp71 = convert(T, 0 // 1) + rp72 = convert(T, -66 // 1) + rp73 = convert(T, 324 // 1) + rp74 = convert(T, -450 // 1) + rp75 = convert(T, 192 // 1) + rp76 = convert(T, 0 // 1) FineRKN5ConstantCache(c1, c2, c3, c4, c5, c6, c7, c8, c9, a21, a31, a32, a41, a43, a51, a52, a53, a54, a61, a62, a63, a64, a71, a73, a74, a75, a81, a82, a83, a84, a85, a87, a91, a92, a93, a94, a95, a97, @@ -461,12 +561,17 @@ function FineRKN5ConstantCache(T::Type, T2::Type) abar71, abar73, abar74, abar75, abar76, abar81, abar82, abar83, abar84, abar85, abar87, abar91, abar93, abar94, abar95, abar97, b1, b3, b4, b5, bbar1, bbar3, bbar4, bbar5, bbar6, btilde1, btilde3, btilde4, btilde5, bptilde1, - bptilde3, bptilde4, bptilde5, bptilde6, bptilde7, r12, r13, r14, r15, r16, r33, r34, - r35, r36, r43, - r44, r45, r46, r53, r54, r55, r56, r63, r64, r65, r66, r73, r74, r75, r76, r83, r84, - r85, r86, r93, r94, r95, r96, rp11, rp12, rp13, rp14, rp15, rp32, rp33, rp34, rp35, - rp42, rp43, rp44, rp45, rp52, rp53, rp54, rp55, rp62, rp63, rp64, rp65, rp72, rp73, - rp74, rp75, rp82, rp83, rp84, rp85, rp92, rp93, rp94, rp95) + bptilde3, bptilde4, bptilde5, bptilde6, bptilde7, btilde1_1, + btilde2_1, btilde3_1, btilde4_1, btilde5_1, btilde6_1, btilde7_1, btilde8_1, btilde9_1, btilde1_12, btilde2_12, btilde3_12, btilde4_12, btilde5_12, btilde6_12, btilde7_12, btilde8_12, btilde9_12, r10, + r11, + r12, r13, r14, r15, r16, r20, r21, r22, r23, r24, r25, r26, r30, r31, r32, r33, r34, + r35, r36, r40, r41, r42, r43, r44, r45, r46, r50, r51, r52, r53, r54, r55, r56, r60, + r61, r62, r63, r64, r65, r66, r70, r71, r72, r73, r74, r75, r76, rp10, rp11, rp12, + rp13, rp14, rp15, rp16, rp20, rp21, rp22, rp23, rp24, rp25, rp26, rp30, rp31, rp32, + rp33, rp34, rp35, rp36, rp40, rp41, rp42, rp43, rp44, rp45, rp46, rp50, rp51, rp52, + rp53, rp54, rp55, rp56, rp60, rp61, rp62, rp63, rp64, rp65, rp66, rp70, rp71, rp72, + rp73, rp74, + rp75, rp76) end struct IRKN3ConstantCache{T, T2} <: OrdinaryDiffEqConstantCache