Skip to content

Zero Dissipation ERK6(5) #1033

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 7 commits into from
Feb 22, 2020
Merged

Zero Dissipation ERK6(5) #1033

Show file tree
Hide file tree
Changes from 6 commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
304840e
Add FRK65 Butcher Table and basic performstep implementation
utkarsh530 Feb 14, 2020
51d4f37
WIP: Add Adaptive EEST and sim tests
utkarsh530 Feb 14, 2020
a5fa042
Update src/caches/low_order_rk_caches.jl
ChrisRackauckas Feb 15, 2020
17c83f3
Update src/caches/low_order_rk_caches.jl
ChrisRackauckas Feb 15, 2020
c2e67de
Add zero-dissipation fitted coefficients
utkarsh530 Feb 21, 2020
a4e0642
Add def value to omega
utkarsh530 Feb 21, 2020
c15a279
Add parameterisation to FRK65 and update convergence tests
utkarsh530 Feb 21, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion src/OrdinaryDiffEq.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
3 changes: 2 additions & 1 deletion src/alg_utils.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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])
Expand Down Expand Up @@ -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
Expand Down
5 changes: 5 additions & 0 deletions src/algorithms.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -398,6 +398,11 @@ struct Vern9 <: OrdinaryDiffEqAdaptiveAlgorithm
lazy::Bool
Vern9(;lazy=true) = new(lazy)
end
struct FRK65 <: OrdinaryDiffEqAdaptiveAlgorithm
Copy link
Member

@ChrisRackauckas ChrisRackauckas Feb 21, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Parameterize this

Copy link
Member

@ChrisRackauckas ChrisRackauckas Feb 21, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
struct FRK65 <: OrdinaryDiffEqAdaptiveAlgorithm
struct FRK65{T} <: OrdinaryDiffEqAdaptiveAlgorithm

omega::Float64
Copy link
Member

@ChrisRackauckas ChrisRackauckas Feb 21, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
omega::Float64
omega::T

FRK65(omega=0.0) = new(omega)
end


################################################################################

Expand Down
245 changes: 245 additions & 0 deletions src/caches/low_order_rk_caches.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -617,3 +617,248 @@ 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,uNoUnitsType,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
atmp::uNoUnitsType
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

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
α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//2351253316)

β1 = T1(1026331676//33222204855)
#β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(1157843020//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)

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

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)
atmp = similar(u,uEltypeNoUnits)
tmp = similar(u)
FRK65Cache(u, uprev, utilde, k1, k2, k3, k4, k5, k6, k7, k8, k9, tmp, atmp, tab)
end
Loading