Skip to content

Commit d565aed

Browse files
ChrisRackauckasChrisRackauckas
authored
Merge pull request #1047 from utkarsh530/PFRK87
Phase Fitted RK8(7) with zero dissipation
2 parents 75ec5dd + 1f870c8 commit d565aed

File tree

7 files changed

+383
-1
lines changed

7 files changed

+383
-1
lines changed

‎src/OrdinaryDiffEq.jl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -176,7 +176,7 @@ module OrdinaryDiffEq
176176

177177
export FunctionMap, Euler, Heun, Ralston, Midpoint, RK4, ExplicitRK, OwrenZen3, OwrenZen4, OwrenZen5,
178178
BS3, BS5, RK46NL, DP5, Tsit5, DP8, Vern6, Vern7, Vern8, TanYam7, TsitPap8,
179-
Vern9, Feagin10, Feagin12, Feagin14, CompositeAlgorithm, Anas5, RKO65, FRK65
179+
Vern9, Feagin10, Feagin12, Feagin14, CompositeAlgorithm, Anas5, RKO65, FRK65, PFRK87
180180

181181
export SSPRK22, SSPRK33, SSPRK53, SSPRK53_2N1, SSPRK53_2N2, SSPRK63, SSPRK73, SSPRK83, SSPRK432,
182182
SSPRKMSVS32, SSPRKMSVS43, SSPRK932, SSPRK54, SSPRK104

‎src/alg_utils.jl‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -290,6 +290,7 @@ alg_order(alg::Hairer42) = 4
290290
alg_order(alg::Feagin10) = 10
291291
alg_order(alg::Feagin12) = 12
292292
alg_order(alg::Feagin14) = 14
293+
alg_order(alg::PFRK87) = 8
293294

294295
alg_order(alg::Rosenbrock23) = 2
295296
alg_order(alg::Rosenbrock32) = 3

‎src/algorithms.jl‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -403,6 +403,11 @@ struct FRK65{T} <: OrdinaryDiffEqAdaptiveAlgorithm
403403
FRK65(omega=0.0) = new{typeof(omega)}(omega)
404404
end
405405

406+
struct PFRK87{T} <: OrdinaryDiffEqAdaptiveAlgorithm
407+
omega::T
408+
PFRK87(omega=0.0) = new{typeof(omega)}(omega)
409+
end
410+
406411

407412
################################################################################
408413

‎src/caches/high_order_rk_caches.jl‎

Lines changed: 47 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -123,3 +123,50 @@ function alg_cache(alg::TsitPap8,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNo
123123
end
124124

125125
alg_cache(alg::TsitPap8,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false}) = TsitPap8ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))
126+
127+
@cache struct PFRK87Cache{uType,rateType,uNoUnitsType,TabType} <: OrdinaryDiffEqMutableCache
128+
u::uType
129+
uprev::uType
130+
fsalfirst::rateType
131+
k2::rateType
132+
k3::rateType
133+
k4::rateType
134+
k5::rateType
135+
k6::rateType
136+
k7::rateType
137+
k8::rateType
138+
k9::rateType
139+
k10::rateType
140+
k11::rateType
141+
k12::rateType
142+
k13::rateType
143+
utilde::uType
144+
tmp::uType
145+
atmp::uNoUnitsType
146+
k::rateType
147+
tab::TabType
148+
end
149+
150+
function alg_cache(alg::PFRK87,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true})
151+
tab = PFRK87ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))
152+
k1 = zero(rate_prototype)
153+
k2 = zero(rate_prototype)
154+
k3 = zero(rate_prototype)
155+
k4 = zero(rate_prototype)
156+
k5 = zero(rate_prototype)
157+
k6 = zero(rate_prototype)
158+
k7 = zero(rate_prototype)
159+
k8 = zero(rate_prototype)
160+
k9 = zero(rate_prototype)
161+
k10 = zero(rate_prototype)
162+
k11 = zero(rate_prototype)
163+
k12 = zero(rate_prototype)
164+
k13 = zero(rate_prototype)
165+
utilde = similar(u)
166+
k = zero(rate_prototype)
167+
tmp = similar(u)
168+
atmp = similar(u,uEltypeNoUnits)
169+
PFRK87Cache(u,uprev,k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,utilde,tmp,atmp,k,tab)
170+
end
171+
172+
alg_cache(alg::PFRK87,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false}) = PFRK87ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))

‎src/perform_step/high_order_rk_perform_step.jl‎

Lines changed: 173 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -601,3 +601,176 @@ end
601601
integrator.destats.nf += 1
602602
end
603603
=#
604+
605+
function initialize!(integrator, cache::PFRK87ConstantCache)
606+
integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
607+
integrator.destats.nf += 1
608+
integrator.kshortsize = 2
609+
integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
610+
611+
# Avoid undefined entries if k is an array of arrays
612+
integrator.fsallast = zero(integrator.fsalfirst)
613+
integrator.k[1] = integrator.fsalfirst
614+
integrator.k[2] = integrator.fsallast
615+
end
616+
617+
@muladd function perform_step!(integrator, cache::PFRK87ConstantCache, repeat_step=false)
618+
@unpack t,dt,uprev,u,f,p = integrator
619+
@unpack α0201, α0301, α0401, α0501, α0601, α0701, α0302, α0403, α0503, α0504, α0604, α0704, α0605, α0705, α0706, α0908, α1008, α1108, α1208, α1308, α1009, α1109, α1209, α1309, α1110, α1210, α1310, α1211, α1311, β1, β6, β7, β8, β9, β10, β11, β12, β13, β1tilde, β6tilde, β7tilde, β8tilde, β9tilde, β10tilde, β11tilde, β12tilde, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11, c12, c13 = cache
620+
alg = unwrap_alg(integrator, true)
621+
ν = alg.omega*dt
622+
νsq = ν^2
623+
624+
c_ν = -0.19781108078634084 - (0.164050909125528499*νsq) + (0.042578310088756321*(νsq^2)) - (0.002300513610963998*(νsq^3)) + (0.000033467244551879287*(νsq^4)) - (7.8661142036921924*(1/100000000)*(νsq^5))
625+
d_ν = 1-(0.296457092123567400*νsq) + (0.0015793885907465726*(νsq^2)) - (0.00018913011771688527*(νsq^3)) + (0.000017089234650765179*(νsq^4)) - (1.2705211682518626*(1/10000000)*(νsq^5))
626+
627+
α0807 = c_ν/d_ν
628+
α0801 = 0.026876256 + 0.0576576*α0807
629+
α0804 = 0.22464336 - 0.944944*α0807
630+
α0805 = 0.000369024 - 0.2061696*α0807
631+
α0806 = 0.21311136 + 0.093456*α0807
632+
α0901 = 0.07239997637512857 + 0.01913119863380767*α0807
633+
α0904 = -0.688400520601143 - 0.3135390887207368*α0807
634+
α0905 = -0.17301267570583073 - 0.06840852844816077*α0807
635+
α0906 = 0.1440060555560846 + 0.031009360422930017*α0807
636+
α0907 = 0.9982362892760762 + 0.33180705811215994*α0807
637+
α1001 = 0.16261514523236525 - 0.12125171966747463*α0807
638+
α1004 = -2.1255544052061124 + 1.9871809612169453*α0807
639+
α1005 = -0.216403903283323 + 0.43356675517460624*α0807
640+
α1006 = -0.060417230254934076 - 0.1965343807796979*α0807
641+
α1007 = 2.4846281621788395 - 2.102961615944379*α0807
642+
α1101 = -1.0320124180911034 + 1.061943768952537*α0807
643+
α1104 = 13.666683232895137 - 17.40407843561103*α0807
644+
α1105 = 0.25990355211486116 - 3.797253476860588*α0807
645+
α1106 = -5.759316475814002 + 1.7212824826428488*α0807
646+
α1107 = -12.822511612651839 + 18.41810566087623*α0807
647+
α1201 = 0.2478349764611783 - 0.06383934946543009*α0807
648+
α1204 = -4.593782880309185 + 1.046256005127882*α0807
649+
α1205 = -0.39566692537411896 + 0.22827403748244698*α0807
650+
α1206 = -3.0673550479691665 - 0.10347586863902129*α0807
651+
α1207 = 5.386688702227177 - 1.1072148245058775*α0807
652+
α1301 = 0.7332242174431163 - 0.5164807626867616*α0807
653+
α1304 = -10.196728938160977 + 8.464545832921925*α0807
654+
α1305 = -0.43865244706547707 + 1.846809999910238*α0807
655+
α1306 = 0.5693856884667226 - 0.8371528845746959*α0807
656+
α1307 = 10.52865228002416 - 8.957722185570706*α0807
657+
658+
k1 = integrator.fsalfirst
659+
k2 = f(uprev+dt*α0201*k1, p, t + c2*dt)
660+
k3 = f(uprev+dt*(α0301*k1+α0302*k2), p, t + c3*dt)
661+
k4 = f(uprev+dt*(α0401*k1 +α0403*k3), p, t + c4*dt)
662+
k5 = f(uprev+dt*(α0501*k1 +α0503*k3+α0504*k4), p, t + c5*dt)
663+
k6 = f(uprev+dt*(α0601*k1 +α0604*k4+α0605*k5), p, t + c6*dt)
664+
k7 = f(uprev+dt*(α0701*k1 +α0704*k4+α0705*k5+α0706*k6), p, t + c7*dt)
665+
k8 = f(uprev+dt*(α0801*k1 +α0804*k4+α0805*k5+α0806*k6+α0807*k7), p, t + c8*dt)
666+
k9 = f(uprev+dt*(α0901*k1 +α0904*k4+α0905*k5+α0906*k6+α0907*k7+α0908*k8), p, t + c9*dt)
667+
k10 =f(uprev+dt*(α1001*k1 +α1004*k4+α1005*k5+α1006*k6+α1007*k7+α1008*k8+α1009*k9), p, t + c10*dt)
668+
k11= f(uprev+dt*(α1101*k1 +α1104*k4+α1105*k5+α1106*k6+α1107*k7+α1108*k8+α1109*k9+α1110*k10), p, t + c11*dt)
669+
k12= f(uprev+dt*(α1201*k1 +α1204*k4+α1205*k5+α1206*k6+α1207*k7+α1208*k8+α1209*k9+α1210*k10+α1211*k11), p, t+c12*dt)
670+
k13= f(uprev+dt*(α1301*k1 +α1304*k4+α1305*k5+α1306*k6+α1307*k7+α1308*k8+α1309*k9+α1310*k10+α1311*k11), p, t+c13*dt)
671+
integrator.destats.nf += 12
672+
u = uprev + dt*(β1*k1+β6*k6+β7*k7+β8*k8+β9*k9+β10*k10+β11*k11+β12*k12+β13*k13)
673+
if integrator.opts.adaptive
674+
utilde = dt*(β1tilde*k1 + β6tilde*k6 + β7tilde*k7 + β8tilde*k8 + β9tilde*k9 + β10tilde*k10 + β11tilde*k11 + β12tilde*k12)
675+
atmp = calculate_residuals(utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
676+
integrator.EEst = integrator.opts.internalnorm(atmp,t)
677+
end
678+
integrator.fsallast = f(u, p, t+dt)
679+
integrator.destats.nf += 1
680+
integrator.k[1] = integrator.fsalfirst
681+
integrator.k[2] = integrator.fsallast
682+
integrator.u = u
683+
end
684+
685+
686+
function initialize!(integrator, cache::PFRK87Cache)
687+
integrator.fsalfirst = cache.fsalfirst
688+
integrator.fsallast = cache.k
689+
integrator.kshortsize = 2
690+
resize!(integrator.k, integrator.kshortsize)
691+
integrator.k[1] = integrator.fsalfirst
692+
integrator.k[2] = integrator.fsallast
693+
integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
694+
integrator.destats.nf += 1
695+
end
696+
697+
@muladd function perform_step!(integrator, cache::PFRK87Cache, repeat_step=false)
698+
@unpack t,dt,uprev,u,f,p = integrator
699+
@unpack α0201, α0301, α0401, α0501, α0601, α0701, α0302, α0403, α0503, α0504, α0604, α0704, α0605, α0705, α0706, α0908, α1008, α1108, α1208, α1308, α1009, α1109, α1209, α1309, α1110, α1210, α1310, α1211, α1311, β1, β6, β7, β8, β9, β10, β11, β12, β13, β1tilde, β6tilde, β7tilde, β8tilde, β9tilde, β10tilde, β11tilde, β12tilde, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11, c12, c13 = cache.tab
700+
@unpack k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,utilde,tmp,atmp,k = cache
701+
702+
alg = unwrap_alg(integrator, true)
703+
ν = alg.omega*dt
704+
νsq = ν^2
705+
706+
c_ν = -0.19781108078634084 - (0.164050909125528499*νsq) + (0.042578310088756321*(νsq^2)) - (0.002300513610963998*(νsq^3)) + (0.000033467244551879287*(νsq^4)) - (7.8661142036921924*(10^(-8))*(νsq^5))
707+
d_ν = 1-(0.296457092123567400*νsq) + (0.0015793885907465726*(νsq^2)) - (0.00018913011771688527*(νsq^3)) + (0.000017089234650765179*(νsq^4)) - (1.2705211682518626*(10^(-7))*(νsq^5))
708+
709+
α0807 = c_ν/d_ν
710+
α0801 = 0.026876256 + 0.0576576*α0807
711+
α0804 = 0.22464336 - 0.944944*α0807
712+
α0805 = 0.000369024 - 0.2061696*α0807
713+
α0806 = 0.21311136 + 0.093456*α0807
714+
α0901 = 0.07239997637512857 + 0.01913119863380767*α0807
715+
α0904 = -0.688400520601143 - 0.3135390887207368*α0807
716+
α0905 = -0.17301267570583073 - 0.06840852844816077*α0807
717+
α0906 = 0.1440060555560846 + 0.031009360422930017*α0807
718+
α0907 = 0.9982362892760762 + 0.33180705811215994*α0807
719+
α1001 = 0.16261514523236525 - 0.12125171966747463*α0807
720+
α1004 = -2.1255544052061124 + 1.9871809612169453*α0807
721+
α1005 = -0.216403903283323 + 0.43356675517460624*α0807
722+
α1006 = -0.060417230254934076 - 0.1965343807796979*α0807
723+
α1007 = 2.4846281621788395 - 2.102961615944379*α0807
724+
α1101 = -1.0320124180911034 + 1.061943768952537*α0807
725+
α1104 = 13.666683232895137 - 17.40407843561103*α0807
726+
α1105 = 0.25990355211486116 - 3.797253476860588*α0807
727+
α1106 = -5.759316475814002 + 1.7212824826428488*α0807
728+
α1107 = -12.822511612651839 + 18.41810566087623*α0807
729+
α1201 = 0.2478349764611783 - 0.06383934946543009*α0807
730+
α1204 = -4.593782880309185 + 1.046256005127882*α0807
731+
α1205 = -0.39566692537411896 + 0.22827403748244698*α0807
732+
α1206 = -3.0673550479691665 - 0.10347586863902129*α0807
733+
α1207 = 5.386688702227177 - 1.1072148245058775*α0807
734+
α1301 = 0.7332242174431163 - 0.5164807626867616*α0807
735+
α1304 = -10.196728938160977 + 8.464545832921925*α0807
736+
α1305 = -0.43865244706547707 + 1.846809999910238*α0807
737+
α1306 = 0.5693856884667226 - 0.8371528845746959*α0807
738+
α1307 = 10.52865228002416 - 8.957722185570706*α0807
739+
740+
k1 = cache.fsalfirst
741+
f(k1, uprev, p, t)
742+
@.. tmp = uprev+dt*α0201*k1
743+
f(k2, tmp, p, t + c2*dt)
744+
@.. tmp = uprev+dt*(α0301*k1+α0302*k2)
745+
f(k3, tmp, p, t + c3*dt)
746+
@.. tmp = uprev+dt*(α0401*k1+α0403*k3)
747+
f(k4, tmp, p, t + c4*dt)
748+
@.. tmp = uprev+dt*(α0501*k1+α0503*k3+α0504*k4)
749+
f(k5, tmp, p, t + c5*dt)
750+
@.. tmp = uprev+dt*(α0601*k1+α0604*k4+α0605*k5)
751+
f(k6, tmp, p, t + c6*dt)
752+
@.. tmp = uprev+dt*(α0701*k1+α0704*k4+α0705*k5+α0706*k6)
753+
f(k7, tmp, p, t + c7*dt)
754+
@.. tmp = uprev+dt*(α0801*k1+α0804*k4+α0805*k5+α0806*k6+α0807*k7)
755+
f(k8, tmp, p, t + c8*dt)
756+
@.. tmp = uprev+dt*(α0901*k1+α0904*k4+α0905*k5+α0906*k6+α0907*k7+α0908*k8)
757+
f(k9, tmp, p, t + c9*dt)
758+
@.. tmp = uprev+dt*(α1001*k1+α1004*k4+α1005*k5+α1006*k6+α1007*k7+α1008*k8+α1009*k9)
759+
f(k10, tmp, p, t + c10*dt)
760+
@.. tmp = uprev+dt*(α1101*k1+α1104*k4+α1105*k5+α1106*k6+α1107*k7+α1108*k8+α1109*k9+α1110*k10)
761+
f(k11, tmp, p, t + c11*dt)
762+
@.. tmp = uprev+dt*(α1201*k1+α1204*k4+α1205*k5+α1206*k6+α1207*k7+α1208*k8+α1209*k9+α1210*k10+α1211*k11)
763+
f(k12, tmp, p, t+c12*dt)
764+
@.. tmp = uprev+dt*(α1301*k1+α1304*k4+α1305*k5+α1306*k6+α1307*k7+α1308*k8+α1309*k9+α1310*k10+α1311*k11)
765+
f(k13, tmp, p, t+c13*dt)
766+
@.. u = uprev + dt*(β1*k1+β6*k6+β7*k7+β8*k8+β9*k9+β10*k10+β11*k11+β12*k12+β13*k13)
767+
integrator.destats.nf += 13
768+
if integrator.opts.adaptive
769+
@.. utilde = dt*(β1tilde*k1 + β6tilde*k6 + β7tilde*k7 + β8tilde*k8 + β9tilde*k9 + β10tilde*k10 + β11tilde*k11 + β12tilde*k12)
770+
calculate_residuals!(atmp, utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
771+
integrator.EEst = integrator.opts.internalnorm(atmp,t)
772+
end
773+
f(k, u, p, t+dt)
774+
integrator.destats.nf += 1
775+
return nothing
776+
end

0 commit comments

Comments
 (0)