From 304840eb4474d19031d35c5b6eae1abf22d3f90f Mon Sep 17 00:00:00 2001 From: Utkarsh Date: Sat, 15 Feb 2020 00:52:08 +0530 Subject: [PATCH 1/7] Add FRK65 Butcher Table and basic performstep implementation --- src/OrdinaryDiffEq.jl | 2 +- src/alg_utils.jl | 3 +- src/algorithms.jl | 4 + src/caches/low_order_rk_caches.jl | 168 ++++++++++++++++++ src/perform_step/low_order_rk_perform_step.jl | 98 ++++++++++ 5 files changed, 273 insertions(+), 2 deletions(-) diff --git a/src/OrdinaryDiffEq.jl b/src/OrdinaryDiffEq.jl index 695cca1d26..1ba230e341 100644 --- a/src/OrdinaryDiffEq.jl +++ b/src/OrdinaryDiffEq.jl @@ -176,7 +176,7 @@ module OrdinaryDiffEq export FunctionMap, Euler, Heun, Ralston, Midpoint, RK4, ExplicitRK, OwrenZen3, OwrenZen4, OwrenZen5, BS3, BS5, RK46NL, DP5, Tsit5, DP8, Vern6, Vern7, Vern8, TanYam7, TsitPap8, - Vern9, Feagin10, Feagin12, Feagin14, CompositeAlgorithm, Anas5, RKO65 + Vern9, Feagin10, Feagin12, Feagin14, CompositeAlgorithm, Anas5, RKO65, FRK65 export SSPRK22, SSPRK33, SSPRK53, SSPRK53_2N1, SSPRK53_2N2, SSPRK63, SSPRK73, SSPRK83, SSPRK432, SSPRKMSVS32, SSPRKMSVS43, SSPRK932, SSPRK54, SSPRK104 diff --git a/src/alg_utils.jl b/src/alg_utils.jl index 4c637166cd..52da7dfb7e 100644 --- a/src/alg_utils.jl +++ b/src/alg_utils.jl @@ -25,6 +25,7 @@ isfsal(alg::NDBLSRK144) = false isfsal(alg::PDIRK44) = false isfsal(alg::DImplicitEuler) = false isfsal(alg::RKO65) = false +isfsal(alg::FRK65) = true get_current_isfsal(alg, cache) = isfsal(alg) get_current_isfsal(alg::CompositeAlgorithm, cache) = isfsal(alg.algs[cache.current]) @@ -164,7 +165,7 @@ alg_order(alg::Anas5) = 5 alg_order(alg::RK46NL) = 4 alg_order(alg::KuttaPRK2p5) = 5 alg_order(alg::RKO65) = 5 - +alg_order(alg::FRK65) = 6 alg_order(alg::SymplecticEuler) = 1 alg_order(alg::VelocityVerlet) = 2 diff --git a/src/algorithms.jl b/src/algorithms.jl index 272e907126..7a72c8f77a 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -398,6 +398,10 @@ struct Vern9 <: OrdinaryDiffEqAdaptiveAlgorithm lazy::Bool Vern9(;lazy=true) = new(lazy) end +struct FRK65 <: OrdinaryDiffEqAdaptiveAlgorithm + omega::Float64 +end + ################################################################################ diff --git a/src/caches/low_order_rk_caches.jl b/src/caches/low_order_rk_caches.jl index 5795654bf7..45fe41b2d4 100644 --- a/src/caches/low_order_rk_caches.jl +++ b/src/caches/low_order_rk_caches.jl @@ -617,3 +617,171 @@ function alg_cache(alg::RKO65,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUni tab = RKO65ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits)) RKO65Cache(u,uprev,k, k1,k2,k3,k4,k5,k6, tmp, fsalfirst, tab) end + +@cache struct FRK65Cache{uType,rateType,TabType} <: OrdinaryDiffEqMutableCache + u::uType + uprev::uType + utilde::uType + k1::rateType + k2::rateType + k3::rateType + k4::rateType + k5::rateType + k6::rateType + k7::rateType + k8::rateType + k9::rateType + tmp::uType + fsalfirst::rateType + tab::TabType +end + +struct FRK65ConstantCache{T1,T2} <: OrdinaryDiffEqConstantCache + + α21::T1 + α31::T1 + α41::T1 + α51::T1 + α61::T1 + α71::T1 + α81::T1 + α91::T1 + + α32::T1 + + α43::T1 + α53::T1 + α63::T1 + α73::T1 + α83::T1 + + α54::T1 + α64::T1 + α74::T1 + α84::T1 + α94::T1 + + α65::T1 + α75::T1 + α85::T1 + α95::T1 + + α76::T1 + α86::T1 + α96::T1 + + α87::T1 + α97::T1 + + α98::T1 + + β1::T1 + β4::T1 + β5::T1 + β6::T1 + β7::T1 + β8::T1 + + β1tilde::T1 + β4tilde::T1 + β5tilde::T1 + β6tilde::T1 + β7tilde::T1 + β8tilde::T1 + β9tilde::T1 + + c2::T2 + c3::T2 + c4::T2 + c5::T2 + c6::T2 + c7::T2 + c8::T2 + c9::T2 + + function FRK65ConstantCache(T1,T2) + + # elements of Butcher Table + α21 = T1(1//89) + α31 = T1(-38624//142129) + α41 = T1(51//1508) + α51 = T1(3259284578//3517556363) + α61 = T1(-108363632681//45875676369) + α71 = T1(7137368591//11299833148) + α81 = T1(8898824396//9828950919) + α91 = T1(1026331676//33222204855) + + α32 = T1(51442//142129) + + α43 = T1(153//1508) + α53 = T1(-69727055112//19553806387) + α63 = T1(80902506271//8700424616) + α73 = T1(-33088067061//10572251159) + α83 = T1(25673454973//11497947835) + + + α54 = T1(36230363390//11788838981) + α64 = T1(-120088218786//17139312481) + α74 = T1(11481363823//3650030081) + α84 = T1(-74239028301//15737704666) + α94 = T1(1450675392//5936579813) + + α65 = T1(4533285649//6676940598) + α75 = T1(-4096673444//7349814937) + α85 = T1(222688842816//44196813415) + α95 = T1(4617877550//16762182457) + + α76 = T1(9911918171//12847192605) + α86 = T1(-105204445705//30575217706) + α96 = T1(1144867463//6520294355) + + α87 = T1(8799291910//8966990271) + α97 = T1(1822809703//7599996644) + + α98 = T1(79524953//2361253316) + + β1 = T1(1026331676//33222204855) + β4 = T1(1450675392//5936579813) + β5 = T1(4617877550//16762182457) + β6 = T1(114867463//6520294355) + β7 = T1(1822809703//7599996644) + β8 = T1(79524953//235125316) + + β1tilde = T1(413034411//13925408836) + β4tilde = T1(1865954212//7538591735) + β5tilde = T1(4451980162//16576017119) + β6tilde = T1(115784320//6320223511) + β7tilde = T1(802708729//3404369569) + β8tilde = T1(-251398161//17050111121) + β9tilde = T1(1//20) + + c2 = T2(1//89) + c3 = T2(34//377) + c4 = T2(51//377) + c5 = T2(14497158//33407747) + c6 = T2(9744566553//16002998914) + c7 = T2(330//383) + c8 = T2(1//1) + c9 = T2(1//1) + new{T1,T2}(α21, α31, α41, α51, α61, α71, α81, α91, α32, α43, α53, α63, α73, α83, α54, α64, α74, α84, α94, α65, α75, α85, α95, α76, α86, α96, α87, α97, α98, β1, β4, β5, β6, β7, β8, β1tilde, β4tilde, β5tilde, β6tilde, β7tilde, β8tilde, β9tilde, c2, c3, c4, c5, c6, c7, c8, c9) + end +end + +alg_cache(alg::FRK65,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false}) = FRK65ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits)) + +function alg_cache(alg::FRK65,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true}) + tab = FRK65ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits)) + k1 = zero(rate_prototype) + k2 = zero(rate_prototype) + k3 = zero(rate_prototype) + k4 = zero(rate_prototype) + k5 = zero(rate_prototype) + k6 = zero(rate_prototype) + k7 = zero(rate_prototype) + k8 = zero(rate_prototype) + k9 = zero(rate_prototype) + utilde = similar(u) + tmp = similar(u) + fsalfirst = zero(rate_prototype) + FRK65Cache(u, uprev, utilde, k1, k2, k3, k4, k5, k6, k7, k8, k9, tmp, fsalfirst, tab) +end \ No newline at end of file diff --git a/src/perform_step/low_order_rk_perform_step.jl b/src/perform_step/low_order_rk_perform_step.jl index 581df94130..f8a9213a02 100644 --- a/src/perform_step/low_order_rk_perform_step.jl +++ b/src/perform_step/low_order_rk_perform_step.jl @@ -1029,3 +1029,101 @@ end #return nothing end + +function initialize!(integrator, cache::FRK65ConstantCache) + integrator.kshortsize = 9 + integrator.k = typeof(integrator.k)(undef, integrator.kshortsize) + integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal + integrator.destats.nf += 1 + + # Avoid undefined entries if k is an array of arrays + 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::FRK65ConstantCache, repeat_step=false) + @unpack t, dt, uprev, u, f, p = integrator + @unpack α21, α31, α41, α51, α61, α71, α81, α91, α32, α43, α53, α63, α73, α83, α54, α64, α74, α84, α94, α65, α75, α85, α95, α76, α86, α96, α87, α97, α98, β1, β4, β5, β6, β7, β8, β1tilde, β4tilde, β5tilde, β6tilde, β7tilde, β8tilde, β9tilde, c2, c3, c4, c5, c6, c7, c8, c9 = cache + alg = unwrap_alg(integrator, true) + omega = alg.omega + + k1 = integrator.fsalfirst + k2 = f(uprev+α21*dt*k1, p, t+c2*dt) + k3 = f(uprev+α31*dt*k1+α32*dt*k2, p, t+c3*dt) + k4 = f(uprev+α41*dt*k1+α43*dt*k3, p, t+c4*dt) + k5 = f(uprev+α51*dt*k1+α53*dt*k3+α54*dt*k4, p, t+c5*dt) + k6 = f(uprev+α61*dt*k1+α63*dt*k3+α64*dt*k4+α65*dt*k5, p, t+c6*dt) + k7 = f(uprev+α71*dt*k1+α73*dt*k3+α74*dt*k4+α75*dt*k5+α76*dt*k6, p, t+c7*dt) + k8 = f(uprev+α81*dt*k1+α83*dt*k3+α84*dt*k4+α85*dt*k5+α86*dt*k6+α87*dt*k7, p, t+c8*dt) + u = uprev+dt*(β1*k1+β4*k4+β5*k5+β6*k6+β7*k7+β8*k8) + integrator.fsallast = f(u, p, t+dt) + k9 = integrator.fsallast + + integrator.destats.nf += 8 + integrator.k[1]=k1 + integrator.k[2]=k2 + integrator.k[3]=k3 + integrator.k[4]=k4 + integrator.k[5]=k5 + integrator.k[6]=k6 + integrator.k[7]=k7 + integrator.k[8]=k8 + integrator.k[9]=k9 + integrator.u = u +end + + +function initialize!(integrator, cache::FRK65Cache) + integrator.kshortsize = 9 + integrator.fsalfirst = cache.k1 + integrator.fsallast = cache.k9 + + resize!(integrator.k, integrator.kshortsize) + + integrator.k[1]=cache.k1 + integrator.k[2]=cache.k2 + integrator.k[3]=cache.k3 + integrator.k[4]=cache.k4 + integrator.k[5]=cache.k5 + integrator.k[6]=cache.k6 + integrator.k[7]=cache.k7 + integrator.k[8]=cache.k8 + integrator.k[9]=cache.k9 + + integrator.f(integrator.fsalfirst,integrator.uprev,integrator.p,integrator.t) # Pre-start fsal + integrator.destats.nf += 1 +end + + + +@muladd function perform_step!(integrator, cache::FRK65Cache, repeat_step=false) + @unpack t, dt, uprev, u, f, p = integrator + @unpack tmp, k1, k2, k3, k4, k5, k6, k7, k8, k9, utilde = cache + @unpack α21, α31, α41, α51, α61, α71, α81, α91, α32, α43, α53, α63, α73, α83, α54, α64, α74, α84, α94, α65, α75, α85, α95, α76, α86, α96, α87, α97, α98, β1, β4, β5, β6, β7, β8, β1tilde, β4tilde, β5tilde, β6tilde, β7tilde, β8tilde, β9tilde, c2, c3, c4, c5, c6, c7, c8, c9 = cache.tab + + @.. tmp = uprev+α21*dt*k1 + f(k2, tmp, p, t+c2*dt) + @.. tmp = uprev+α31*dt*k1+α32*dt*k2 + f(k3, tmp, p, t+c3*dt) + @.. tmp = uprev+α41*dt*k1+α43*dt*k3 + f(k4, tmp, p, t+c4*dt) + @.. tmp = uprev+α51*dt*k1+α53*dt*k3+α54*dt*k4 + f(k5, tmp, p, t+c5*dt) + @.. tmp = uprev+α61*dt*k1+α63*dt*k3+α64*dt*k4+α65*dt*k5 + f(k6, tmp, p, t+c6*dt) + @.. tmp = uprev+α71*dt*k1+α73*dt*k3+α74*dt*k4+α75*dt*k5+α76*dt*k6 + f(k7, tmp, p, t+c7*dt) + @.. tmp = uprev+α81*dt*k1+α83*dt*k3+α84*dt*k4+α85*dt*k5+α86*dt*k6+α87*dt*k7 + f(k8, tmp, p, t+c5*dt) + @.. u = uprev+dt*(β1*k1+β4*k4+β5*k5+β6*k6+β7*k7+β8*k8) + f(k9, u, p, t+dt) + + integrator.destats.nf += 8 + + #return nothing +end + From 51d4f375275ae5937c2734d956ca1c0777a31362 Mon Sep 17 00:00:00 2001 From: Utkarsh Date: Sat, 15 Feb 2020 04:57:03 +0530 Subject: [PATCH 2/7] WIP: Add Adaptive EEST and sim tests --- src/caches/low_order_rk_caches.jl | 8 ++++---- src/perform_step/low_order_rk_perform_step.jl | 20 ++++++++++++++----- test/algconvergence/ode_convergence_tests.jl | 3 +++ 3 files changed, 22 insertions(+), 9 deletions(-) diff --git a/src/caches/low_order_rk_caches.jl b/src/caches/low_order_rk_caches.jl index 45fe41b2d4..43787e1581 100644 --- a/src/caches/low_order_rk_caches.jl +++ b/src/caches/low_order_rk_caches.jl @@ -618,7 +618,7 @@ function alg_cache(alg::RKO65,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUni RKO65Cache(u,uprev,k, k1,k2,k3,k4,k5,k6, tmp, fsalfirst, tab) end -@cache struct FRK65Cache{uType,rateType,TabType} <: OrdinaryDiffEqMutableCache +@cache struct FRK65Cache{uType,rateType,uNoUnitsType,TabType} <: OrdinaryDiffEqMutableCache u::uType uprev::uType utilde::uType @@ -632,7 +632,7 @@ end k8::rateType k9::rateType tmp::uType - fsalfirst::rateType + atmp::uNoUnitsType tab::TabType end @@ -781,7 +781,7 @@ function alg_cache(alg::FRK65,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUni k8 = zero(rate_prototype) k9 = zero(rate_prototype) utilde = similar(u) + atmp = similar(u,uEltypeNoUnits) tmp = similar(u) - fsalfirst = zero(rate_prototype) - FRK65Cache(u, uprev, utilde, k1, k2, k3, k4, k5, k6, k7, k8, k9, tmp, fsalfirst, tab) + FRK65Cache(u, uprev, utilde, k1, k2, k3, k4, k5, k6, k7, k8, k9, tmp, atmp, tab) end \ No newline at end of file diff --git a/src/perform_step/low_order_rk_perform_step.jl b/src/perform_step/low_order_rk_perform_step.jl index f8a9213a02..d8822d68f1 100644 --- a/src/perform_step/low_order_rk_perform_step.jl +++ b/src/perform_step/low_order_rk_perform_step.jl @@ -1062,7 +1062,11 @@ end u = uprev+dt*(β1*k1+β4*k4+β5*k5+β6*k6+β7*k7+β8*k8) integrator.fsallast = f(u, p, t+dt) k9 = integrator.fsallast - + if integrator.opts.adaptive + utilde = dt*(β1tilde*k1+β4tilde*k4+β5tilde*k5+β6tilde*k6+β7tilde*k7+β8tilde*k8+β9tilde*k9) + atmp = calculate_residuals(utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t) + integrator.EEst = integrator.opts.internalnorm(atmp,t) + end integrator.destats.nf += 8 integrator.k[1]=k1 integrator.k[2]=k2 @@ -1102,7 +1106,7 @@ end @muladd function perform_step!(integrator, cache::FRK65Cache, repeat_step=false) @unpack t, dt, uprev, u, f, p = integrator - @unpack tmp, k1, k2, k3, k4, k5, k6, k7, k8, k9, utilde = cache + @unpack tmp, k1, k2, k3, k4, k5, k6, k7, k8, k9, utilde, atmp = cache @unpack α21, α31, α41, α51, α61, α71, α81, α91, α32, α43, α53, α63, α73, α83, α54, α64, α74, α84, α94, α65, α75, α85, α95, α76, α86, α96, α87, α97, α98, β1, β4, β5, β6, β7, β8, β1tilde, β4tilde, β5tilde, β6tilde, β7tilde, β8tilde, β9tilde, c2, c3, c4, c5, c6, c7, c8, c9 = cache.tab @.. tmp = uprev+α21*dt*k1 @@ -1118,12 +1122,18 @@ end @.. tmp = uprev+α71*dt*k1+α73*dt*k3+α74*dt*k4+α75*dt*k5+α76*dt*k6 f(k7, tmp, p, t+c7*dt) @.. tmp = uprev+α81*dt*k1+α83*dt*k3+α84*dt*k4+α85*dt*k5+α86*dt*k6+α87*dt*k7 - f(k8, tmp, p, t+c5*dt) + f(k8, tmp, p, t+c8*dt) @.. u = uprev+dt*(β1*k1+β4*k4+β5*k5+β6*k6+β7*k7+β8*k8) f(k9, u, p, t+dt) integrator.destats.nf += 8 - - #return nothing + if integrator.opts.adaptive + @.. utilde = dt*(β1tilde*k1+β4tilde*k4+β5tilde*k5+β6tilde*k6+β7tilde*k7+β8tilde*k8+β9tilde*k9) + calculate_residuals(atmp, utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t) + integrator.EEst = integrator.opts.internalnorm(atmp,t) + end + return nothing end + + diff --git a/test/algconvergence/ode_convergence_tests.jl b/test/algconvergence/ode_convergence_tests.jl index 3041740b95..d5dd5a8339 100644 --- a/test/algconvergence/ode_convergence_tests.jl +++ b/test/algconvergence/ode_convergence_tests.jl @@ -36,6 +36,9 @@ testTol = 0.2 sim3 = test_convergence(dts2,prob,RKO65()) @test sim3.𝒪est[:l∞] ≈ 5 atol=testTol + sim3 = test_convergence(dts2,prob,FRK65(10)) + @test sim3.𝒪est[:l∞] ≈ 6 atol=testTol + sim4 = test_convergence(dts,prob,BS3()) @test sim4.𝒪est[:l2] ≈ 3 atol=testTol sim5 = test_convergence(dts, prob, AB3()) From a5fa04275ed96608475057744b63e1b58ec2c89e Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sat, 15 Feb 2020 16:42:10 -0500 Subject: [PATCH 3/7] Update src/caches/low_order_rk_caches.jl Co-Authored-By: Kanav Gupta <33966400+kanav99@users.noreply.github.com> --- src/caches/low_order_rk_caches.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/caches/low_order_rk_caches.jl b/src/caches/low_order_rk_caches.jl index 43787e1581..f7f3faef57 100644 --- a/src/caches/low_order_rk_caches.jl +++ b/src/caches/low_order_rk_caches.jl @@ -743,7 +743,7 @@ struct FRK65ConstantCache{T1,T2} <: OrdinaryDiffEqConstantCache β1 = T1(1026331676//33222204855) β4 = T1(1450675392//5936579813) β5 = T1(4617877550//16762182457) - β6 = T1(114867463//6520294355) + β6 = T1(1144867463//6520294355) β7 = T1(1822809703//7599996644) β8 = T1(79524953//235125316) @@ -784,4 +784,4 @@ function alg_cache(alg::FRK65,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUni atmp = similar(u,uEltypeNoUnits) tmp = similar(u) FRK65Cache(u, uprev, utilde, k1, k2, k3, k4, k5, k6, k7, k8, k9, tmp, atmp, tab) -end \ No newline at end of file +end From 17c83f3a1a03f0b8553dc7e95cdfb6d8bb8b4272 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sat, 15 Feb 2020 16:51:20 -0500 Subject: [PATCH 4/7] Update src/caches/low_order_rk_caches.jl Co-Authored-By: Kanav Gupta <33966400+kanav99@users.noreply.github.com> --- src/caches/low_order_rk_caches.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/caches/low_order_rk_caches.jl b/src/caches/low_order_rk_caches.jl index f7f3faef57..4d16fb5a0b 100644 --- a/src/caches/low_order_rk_caches.jl +++ b/src/caches/low_order_rk_caches.jl @@ -745,7 +745,7 @@ struct FRK65ConstantCache{T1,T2} <: OrdinaryDiffEqConstantCache β5 = T1(4617877550//16762182457) β6 = T1(1144867463//6520294355) β7 = T1(1822809703//7599996644) - β8 = T1(79524953//235125316) + β8 = T1(79524953//2351253316) β1tilde = T1(413034411//13925408836) β4tilde = T1(1865954212//7538591735) From c2e67de0e09c121570a80c9821650012c790cb42 Mon Sep 17 00:00:00 2001 From: Utkarsh Date: Fri, 21 Feb 2020 12:13:23 +0530 Subject: [PATCH 5/7] Add zero-dissipation fitted coefficients --- src/algorithms.jl | 2 +- src/caches/low_order_rk_caches.jl | 95 +++++++++++++++++-- src/perform_step/low_order_rk_perform_step.jl | 25 +++-- test/algconvergence/ode_convergence_tests.jl | 3 +- 4 files changed, 107 insertions(+), 18 deletions(-) diff --git a/src/algorithms.jl b/src/algorithms.jl index 7a72c8f77a..dcaeff6042 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -399,7 +399,7 @@ struct Vern9 <: OrdinaryDiffEqAdaptiveAlgorithm Vern9(;lazy=true) = new(lazy) end struct FRK65 <: OrdinaryDiffEqAdaptiveAlgorithm - omega::Float64 + omega::Float64 #if not specified set default 0(Basic RK6(5))? end diff --git a/src/caches/low_order_rk_caches.jl b/src/caches/low_order_rk_caches.jl index 4d16fb5a0b..fa2ec96dc4 100644 --- a/src/caches/low_order_rk_caches.jl +++ b/src/caches/low_order_rk_caches.jl @@ -676,9 +676,9 @@ struct FRK65ConstantCache{T1,T2} <: OrdinaryDiffEqConstantCache α98::T1 β1::T1 - β4::T1 - β5::T1 - β6::T1 + # β4::T1 + # β5::T1 + # β6::T1 β7::T1 β8::T1 @@ -699,6 +699,44 @@ struct FRK65ConstantCache{T1,T2} <: OrdinaryDiffEqConstantCache c8::T2 c9::T2 + d1::T1 + d2::T1 + d3::T1 + d4::T1 + d5::T1 + d6::T1 + d7::T1 + d8::T1 + d9::T1 + d10::T1 + d11::T1 + d12::T1 + d13::T1 + + e1::T1 + e2::T1 + e3::T1 + e4::T1 + e5::T1 + e6::T1 + e7::T1 + e8::T1 + e9::T1 + e10::T1 + e11::T1 + + f1::T1 + f2::T1 + f3::T1 + f4::T1 + f5::T1 + f6::T1 + f7::T1 + f8::T1 + f9::T1 + f10::T1 + f11::T1 + function FRK65ConstantCache(T1,T2) # elements of Butcher Table @@ -738,19 +776,19 @@ struct FRK65ConstantCache{T1,T2} <: OrdinaryDiffEqConstantCache α87 = T1(8799291910//8966990271) α97 = T1(1822809703//7599996644) - α98 = T1(79524953//2361253316) + α98 = T1(79524953//2351253316) β1 = T1(1026331676//33222204855) - β4 = T1(1450675392//5936579813) - β5 = T1(4617877550//16762182457) - β6 = T1(1144867463//6520294355) + #β4 = T1(1450675392//5936579813) + #β5 = T1(4617877550//16762182457) + #β6 = T1(1144867463//6520294355) β7 = T1(1822809703//7599996644) β8 = T1(79524953//2351253316) β1tilde = T1(413034411//13925408836) β4tilde = T1(1865954212//7538591735) β5tilde = T1(4451980162//16576017119) - β6tilde = T1(115784320//6320223511) + β6tilde = T1(1157843020//6320223511) β7tilde = T1(802708729//3404369569) β8tilde = T1(-251398161//17050111121) β9tilde = T1(1//20) @@ -763,7 +801,46 @@ struct FRK65ConstantCache{T1,T2} <: OrdinaryDiffEqConstantCache c7 = T2(330//383) c8 = T2(1//1) c9 = T2(1//1) - new{T1,T2}(α21, α31, α41, α51, α61, α71, α81, α91, α32, α43, α53, α63, α73, α83, α54, α64, α74, α84, α94, α65, α75, α85, α95, α76, α86, α96, α87, α97, α98, β1, β4, β5, β6, β7, β8, β1tilde, β4tilde, β5tilde, β6tilde, β7tilde, β8tilde, β9tilde, c2, c3, c4, c5, c6, c7, c8, c9) + + d1 = T1(140209127//573775965) + d2 = T1(-8530039//263747097) + d3 = T1(-308551//104235790) + d4 = T1(233511//333733259) + d5 = T1(9126//184950985) + d6 = T1(22//50434083) + d7 = T1(19//424427471) + d8 = T1(-28711//583216934059) + d9 = T1(-3831531//316297807) + d10 = T1(551767//187698280) + d11 = T1(9205//210998423) + d12 = T1(-250//519462673) + d13 = T1(67//327513887) + + e1=T1(437217689//1587032700) + e2=T1(-15824413//592362279) + e3=T1(-1563775//341846569) + e4=T1(270497//369611210) + e5=T1(-26623//453099487) + e6=T1(-616297487849) + e7=T1(-47682337//491732789) + e8=T1(-4778275//287766311) + e9=T1(641177//265376522) + e10=T1(44633//291742143) + e11=T1(611//223639880) + + f1=T1(44861261//255495624) + f2=T1(-11270940//352635157) + f3=T1(-182222//232874507) + f4=T1(164263//307215200) + f5=T1(32184//652060417) + f6=T1(-352//171021903) + f7=T1(-18395427//101056291) + f8=T1(-621686//139501937) + f9=T1(2030024//612171255) + f10=T1(-711049//7105160932) + f11=T1(267//333462710) + + new{T1,T2}(α21, α31, α41, α51, α61, α71, α81, α91, α32, α43, α53, α63, α73, α83, α54, α64, α74, α84, α94, α65, α75, α85, α95, α76, α86, α96, α87, α97, α98, β1, β7, β8, β1tilde, β4tilde, β5tilde, β6tilde, β7tilde, β8tilde, β9tilde, c2, c3, c4, c5, c6, c7, c8, c9, d1, d2, d3, d4, d5, d6, d7, d8, d9, d10, d11, d12, d13, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, e11, f1, f2, f3, f4, f5, f6, f7, f8, f9, f10, f11) end end diff --git a/src/perform_step/low_order_rk_perform_step.jl b/src/perform_step/low_order_rk_perform_step.jl index d8822d68f1..b16668ab4f 100644 --- a/src/perform_step/low_order_rk_perform_step.jl +++ b/src/perform_step/low_order_rk_perform_step.jl @@ -1047,9 +1047,13 @@ end @muladd function perform_step!(integrator, cache::FRK65ConstantCache, repeat_step=false) @unpack t, dt, uprev, u, f, p = integrator - @unpack α21, α31, α41, α51, α61, α71, α81, α91, α32, α43, α53, α63, α73, α83, α54, α64, α74, α84, α94, α65, α75, α85, α95, α76, α86, α96, α87, α97, α98, β1, β4, β5, β6, β7, β8, β1tilde, β4tilde, β5tilde, β6tilde, β7tilde, β8tilde, β9tilde, c2, c3, c4, c5, c6, c7, c8, c9 = cache + @unpack α21, α31, α41, α51, α61, α71, α81, α91, α32, α43, α53, α63, α73, α83, α54, α64, α74, α84, α94, α65, α75, α85, α95, α76, α86, α96, α87, α97, α98, β1, β7, β8, β1tilde, β4tilde, β5tilde, β6tilde, β7tilde, β8tilde, β9tilde, c2, c3, c4, c5, c6, c7, c8, c9, d1, d2, d3, d4, d5, d6, d7, d8, d9, d10, d11, d12, d13, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, e11, f1, f2, f3, f4, f5, f6, f7, f8, f9, f10, f11 = cache alg = unwrap_alg(integrator, true) - omega = alg.omega + ν = alg.omega*dt + νsq = ν^2 + β4 = (d1 + νsq*(d2 + νsq*(d3 + νsq*(d4 + νsq*(d5 + νsq*(d6 + +νsq*d7))))))/(1 + νsq*(d8 + νsq*(d9 + νsq*(d10 + νsq*(d11 + νsq*(d12 + +νsq*d13)))))) + β5 = (e1 + νsq*(e2 + νsq*(e3 + νsq*(e4 + νsq*(e5 + νsq*e6)))))/(1 + νsq*(e8 + νsq*(e9 + νsq*(e10 + νsq*e11)))) + β6 = (f1 + νsq*(f2 + νsq*(f3 + νsq*(f4 + νsq*(f5 + νsq*f6)))))/(1 + νsq*(f8 + νsq*(f9 + νsq*(f10 + νsq*f11)))) k1 = integrator.fsalfirst k2 = f(uprev+α21*dt*k1, p, t+c2*dt) @@ -1063,7 +1067,7 @@ end integrator.fsallast = f(u, p, t+dt) k9 = integrator.fsallast if integrator.opts.adaptive - utilde = dt*(β1tilde*k1+β4tilde*k4+β5tilde*k5+β6tilde*k6+β7tilde*k7+β8tilde*k8+β9tilde*k9) + utilde = dt*(β1tilde*k1 + β4tilde*k4 + β5tilde*k5 + β6tilde*k6 + β7tilde*k7 + β8tilde*k8 + β9tilde*k9) atmp = calculate_residuals(utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t) integrator.EEst = integrator.opts.internalnorm(atmp,t) end @@ -1107,8 +1111,15 @@ end @muladd function perform_step!(integrator, cache::FRK65Cache, repeat_step=false) @unpack t, dt, uprev, u, f, p = integrator @unpack tmp, k1, k2, k3, k4, k5, k6, k7, k8, k9, utilde, atmp = cache - @unpack α21, α31, α41, α51, α61, α71, α81, α91, α32, α43, α53, α63, α73, α83, α54, α64, α74, α84, α94, α65, α75, α85, α95, α76, α86, α96, α87, α97, α98, β1, β4, β5, β6, β7, β8, β1tilde, β4tilde, β5tilde, β6tilde, β7tilde, β8tilde, β9tilde, c2, c3, c4, c5, c6, c7, c8, c9 = cache.tab - + @unpack α21, α31, α41, α51, α61, α71, α81, α91, α32, α43, α53, α63, α73, α83, α54, α64, α74, α84, α94, α65, α75, α85, α95, α76, α86, α96, α87, α97, α98, β1, β7, β8, β1tilde, β4tilde, β5tilde, β6tilde, β7tilde, β8tilde, β9tilde, c2, c3, c4, c5, c6, c7, c8, c9, d1, d2, d3, d4, d5, d6, d7, d8, d9, d10, d11, d12, d13, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, e11, f1, f2, f3, f4, f5, f6, f7, f8, f9, f10, f11 = cache.tab + alg = unwrap_alg(integrator, true) + + ν = alg.omega*dt + νsq = ν^2 + β4 = (d1 + νsq*(d2 + νsq*(d3 + νsq*(d4 + νsq*(d5 + νsq*(d6 + +νsq*d7))))))/(1 + νsq*(d8 + νsq*(d9 + νsq*(d10 + νsq*(d11 + νsq*(d12 + +νsq*d13)))))) + β5 = (e1 + νsq*(e2 + νsq*(e3 + νsq*(e4 + νsq*(e5 + νsq*e6)))))/(1 + νsq*(e8 + νsq*(e9 + νsq*(e10 + νsq*e11)))) + β6 = (f1 + νsq*(f2 + νsq*(f3 + νsq*(f4 + νsq*(f5 + νsq*f6)))))/(1 + νsq*(f8 + νsq*(f9 + νsq*(f10 + νsq*f11)))) + @.. tmp = uprev+α21*dt*k1 f(k2, tmp, p, t+c2*dt) @.. tmp = uprev+α31*dt*k1+α32*dt*k2 @@ -1128,8 +1139,8 @@ end integrator.destats.nf += 8 if integrator.opts.adaptive - @.. utilde = dt*(β1tilde*k1+β4tilde*k4+β5tilde*k5+β6tilde*k6+β7tilde*k7+β8tilde*k8+β9tilde*k9) - calculate_residuals(atmp, utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t) + @.. utilde = dt*(β1tilde*k1 + β4tilde*k4 + β5tilde*k5 + β6tilde*k6 + β7tilde*k7 + β8tilde*k8 + β9tilde*k9) + calculate_residuals!(atmp, utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t) integrator.EEst = integrator.opts.internalnorm(atmp,t) end return nothing diff --git a/test/algconvergence/ode_convergence_tests.jl b/test/algconvergence/ode_convergence_tests.jl index d5dd5a8339..f70b5a9326 100644 --- a/test/algconvergence/ode_convergence_tests.jl +++ b/test/algconvergence/ode_convergence_tests.jl @@ -10,6 +10,7 @@ ODEProblemLibrary.importodeproblems() dts1 = 1 .//2 .^(9:-1:5) dts2 = 1 .//2 .^(7:-1:3) dts3 = 1 .//2 .^(12:-1:7) +dts4 = 1 .//2 .^(5:-1:3) testTol = 0.2 @testset "Explicit Solver Convergence Tests ($(["out-of-place", "in-place"][i]))" for i in 1:2 @@ -36,7 +37,7 @@ testTol = 0.2 sim3 = test_convergence(dts2,prob,RKO65()) @test sim3.𝒪est[:l∞] ≈ 5 atol=testTol - sim3 = test_convergence(dts2,prob,FRK65(10)) + sim3 = test_convergence(dts4,prob,FRK65(0)) #standard value for non zero-dissipation tests @test sim3.𝒪est[:l∞] ≈ 6 atol=testTol sim4 = test_convergence(dts,prob,BS3()) From a4e064258857bb18ffa1d0042b66f7536def55ed Mon Sep 17 00:00:00 2001 From: Utkarsh Date: Fri, 21 Feb 2020 16:07:40 +0530 Subject: [PATCH 6/7] Add def value to omega --- src/algorithms.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/algorithms.jl b/src/algorithms.jl index dcaeff6042..15d277aade 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -399,7 +399,8 @@ struct Vern9 <: OrdinaryDiffEqAdaptiveAlgorithm Vern9(;lazy=true) = new(lazy) end struct FRK65 <: OrdinaryDiffEqAdaptiveAlgorithm - omega::Float64 #if not specified set default 0(Basic RK6(5))? + omega::Float64 + FRK65(omega=0.0) = new(omega) end From c15a2792bf3aa532a5e484a0baa6370962605c1c Mon Sep 17 00:00:00 2001 From: Utkarsh Date: Fri, 21 Feb 2020 22:31:46 +0530 Subject: [PATCH 7/7] Add parameterisation to FRK65 and update convergence tests --- src/algorithms.jl | 6 +++--- test/algconvergence/ode_convergence_tests.jl | 4 ++-- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/algorithms.jl b/src/algorithms.jl index 15d277aade..64f958c0d4 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -398,9 +398,9 @@ struct Vern9 <: OrdinaryDiffEqAdaptiveAlgorithm lazy::Bool Vern9(;lazy=true) = new(lazy) end -struct FRK65 <: OrdinaryDiffEqAdaptiveAlgorithm - omega::Float64 - FRK65(omega=0.0) = new(omega) +struct FRK65{T} <: OrdinaryDiffEqAdaptiveAlgorithm + omega::T + FRK65(omega=0.0) = new{typeof(omega)}(omega) end diff --git a/test/algconvergence/ode_convergence_tests.jl b/test/algconvergence/ode_convergence_tests.jl index f70b5a9326..2a18b50f73 100644 --- a/test/algconvergence/ode_convergence_tests.jl +++ b/test/algconvergence/ode_convergence_tests.jl @@ -37,8 +37,8 @@ testTol = 0.2 sim3 = test_convergence(dts2,prob,RKO65()) @test sim3.𝒪est[:l∞] ≈ 5 atol=testTol - sim3 = test_convergence(dts4,prob,FRK65(0)) #standard value for non zero-dissipation tests - @test sim3.𝒪est[:l∞] ≈ 6 atol=testTol + sim3 = test_convergence(dts4,prob,FRK65()) + @test sim3.𝒪est[:l∞] ≈ 6 atol=0.6 sim4 = test_convergence(dts,prob,BS3()) @test sim4.𝒪est[:l2] ≈ 3 atol=testTol