Skip to content

Commit

Permalink
caches
Browse files Browse the repository at this point in the history
  • Loading branch information
@Shreyas-Ekanathan
Shreyas-Ekanathan committed Sep 12, 2024
1 parent 3318e37 commit 9c15269
Show file tree
Hide file tree
Showing 3 changed files with 45 additions and 21 deletions.
7 changes: 6 additions & 1 deletion lib/OrdinaryDiffEqFIRK/src/firk_caches.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -527,10 +527,12 @@ mutable struct AdaptiveRadauCache{uType, cuType, uNoUnitsType, rateType, JType,
uprev::uType
z::Vector{uType}
w::Vector{uType}
c_prime::Vector{BigFloat}
dw1::uType
ubuff::uType
dw2::Vector{cuType}
cubuff::Vector{cuType}
dw::Vector{uType}
cont::Vector{uType}
derivatives:: Matrix{uType}
du1::rateType
Expand Down Expand Up @@ -587,12 +589,15 @@ function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits}
z[i] = w[i] = zero(u)
end

c_prime = Vector{BigFloat}(undef, num_stages) #time stepping

dw1 = zero(u)
ubuff = zero(u)
dw2 = [similar(u, Complex{eltype(u)}) for _ in 1 : (num_stages - 1) ÷ 2]
recursivefill!.(dw2, false)
cubuff = [similar(u, Complex{eltype(u)}) for _ in 1 : (num_stages - 1) ÷ 2]
recursivefill!.(cubuff, false)
dw = Vector{typeof(u)}(undef, num_stages - 1)

cont = Vector{typeof(u)}(undef, num_stages)
for i in 1 : num_stages
Expand Down Expand Up @@ -641,7 +646,7 @@ function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits}
atol = reltol isa Number ? reltol : zero(reltol)

AdaptiveRadauCache(u, uprev,
z, w, dw1, ubuff, dw2, cubuff, cont, derivatives,
z, w, c_prime, dw1, ubuff, dw2, cubuff, dw, cont, derivatives,
du1, fsalfirst, ks, k, fw,
J, W1, W2,
uf, tab, κ, one(uToltype), 10000, tmp,
Expand Down
18 changes: 9 additions & 9 deletions lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1398,7 +1398,7 @@ end
cache.cont[i] = map(zero, u)
end
else
c_prime = Vector{eltype(u)}(undef, num_stages) #time stepping
c_prime = Vector{typeof(u)}(undef, num_stages) #time stepping
c_prime[num_stages] = dt / cache.dtprev
for i in 1 : num_stages - 1
c_prime[i] = c[i] * c_prime[num_stages]
Expand Down Expand Up @@ -1491,8 +1491,10 @@ end
break
end
end

w = @.. w - dw

for i in 1 : num_stages
w[i] -= dw[i]
end

# transform `w` to `z`
#z = T * w
Expand Down Expand Up @@ -1574,8 +1576,8 @@ end
@muladd function perform_step!(integrator, cache::AdaptiveRadauCache, repeat_step = false)
@unpack t, dt, uprev, u, f, p, fsallast, fsalfirst = integrator
@unpack T, TI, γ, α, β, c, #=e,=# num_stages = cache.tab
@unpack κ, cont, derivatives, z, w = cache
@unpack dw1, ubuff, dw2, cubuff = cache
@unpack κ, cont, derivatives, z, w, c_prime = cache
@unpack dw1, ubuff, dw2, cubuff, dw = cache
@unpack ks, k, fw, J, W1, W2 = cache
@unpack tmp, atmp, jac_config, linsolve1, linsolve2, rtol, atol, step_limiter! = cache
@unpack internalnorm, abstol, reltol, adaptive = integrator.opts
Expand Down Expand Up @@ -1604,7 +1606,6 @@ end
z[i] = w[i] = cache.cont[i] = map(zero, u)
end
else
c_prime = Vector{eltype(u)}(undef, num_stages) #time stepping
c_prime[num_stages] = dt / cache.dtprev
for i in 1 : num_stages - 1
c_prime[i] = c[i] * c_prime[num_stages]
Expand Down Expand Up @@ -1692,7 +1693,6 @@ end
end

integrator.stats.nsolve += (num_stages + 1) / 2
dw = Vector{typeof(u)}(undef, num_stages - 1)

for i in 1 : (num_stages - 1) ÷ 2
dw[2 * i - 1] = z[2 * i - 1]
Expand Down Expand Up @@ -1757,7 +1757,7 @@ end
@.. broadcast=false u=uprev + z[num_stages]

step_limiter!(u, integrator, p, t + dt)
#=

if adaptive
utilde = w2
edt = e./dt
Expand Down Expand Up @@ -1795,7 +1795,7 @@ end
integrator.EEst = internalnorm(atmp, t)
end
end
=#

if integrator.EEst <= oneunit(integrator.EEst)
cache.dtprev = dt
if alg.extrapolant != :constant
Expand Down
41 changes: 30 additions & 11 deletions lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -134,6 +134,11 @@ function BigRadauIIA5Tableau(T1, T2)
c[2] = big"0.644948974278317809819728407470589139196594748065667012843269256725096037745731502653985943310464023481859460122661418912485886545983775734162578395123729143"
c[3] = big"1"

e = Vector{T1}(undef, 3)
e[1] = big"-0.804701356815835379608495496358640916569322134539215617920280276511680200030933806355291481868922518805459899199875734619185214695254668403298825163805293365"
e[2] = big"-0.267446751803505087778945794929857825182629030352446373700645445786652166171126710221557624849436998601914181921833023811045195644454184323577861370786198096"
e[3] = big"-0.202740720976336900360387312310916037542642249112497662477129527815832849767696924955091220025846425370353707227521045881781871836521093030105009933962412198"

TI = Matrix{T1}(undef, 3, 3)
TI[1, 1] = big"4.32557989006315535102435095295614882731995158490590784287320458848019483341979047442263696495019938973156007686663488090615420049217658854859024016717169837"
TI[1, 2] = big"0.339199251815809869542824974053410987511771566126056902312311333553438988409693737874718833892037643701271502187763370262948704203562215007824701228014200056"
Expand All @@ -156,7 +161,7 @@ function BigRadauIIA5Tableau(T1, T2)
T[3, 2] = big"1.0"
T[3, 3] = big"0.0"
RadauIIATableau{T1, T2}(T, TI,
c, γ, α, β, 3)
c, γ, α, β, #=e,=# 3)
end

function BigRadauIIA9Tableau(T1, T2)
Expand All @@ -175,6 +180,13 @@ function BigRadauIIA9Tableau(T1, T2)
c[4] = big"0.8602401356562194478479129188751197667383780225872255049242335941839742579301655644134901549264276106897445531811874851737136468026848125542506920602484255"
c[5] = big"1.0"

e = Vector{T1}(undef, 5)
e[1] = big"-0.396056873040772391443753928838733350903268235649241407109949157055321706077305169410373530093363946563049059516774269126208180048957098799522070282580085012"
e[2] = big"-0.120998893046492111917470082824942310714828787321581605224083897201222931079742440239750023176022706685332409182564676958242591761944047403771685667130014796"
e[3] = big"-0.428099657316704068620981167438991676522506275261025355075012671900974663138721382194232832263754491637468895406090947322821492290199844487667268567641865902"
e[4] = big"-0.14209725213800672440012694437858306290960212386817941668712296105013559223718245243280058605734579343272550597862802439475182620322497804597976958854332963"
e[5] = big"-0.0718131688854938240955830308703126002625513555250069460240421718019137231332378610692892428976217345796439675210144794505060225760814921842351985264738071248"

TI = Matrix{T1}(undef, 5, 5)
TI[1, 1] = big"30.0415677215444016277146611632467970747634862837368422955138470463852339244593400023985957753164599415374157317627305099177616927640413043608408838747985125"
TI[1, 2] = big"13.8651078562714131651762946846279728486098595017962436746405940971751244384714668104145151259298432908422191238542910724677205181071665482818120092330632702"
Expand Down Expand Up @@ -230,7 +242,7 @@ function BigRadauIIA9Tableau(T1, T2)
T[5, 5] = big"0.0"

RadauIIATableau{T1, T2}(T, TI,
c, γ, α, β, 5)
c, γ, α, β, #=e,=# 5)
end


Expand Down Expand Up @@ -402,6 +414,15 @@ function BigRadauIIA13Tableau(T1, T2)
c[6] = big"0.926945671319741114851873965819682011056172419542283252724467079656645202452528243814339480013587391545656707320049986592771178724621938506933715568048004783"
c[7] = big"1.0"

e = Vector{T1}(undef, 7)
e[1] = big"-0.252864387074253458829206302755008427271632071911523081751767012887135523702827548481788009254605639434091010898676803950224099333883157851379911155069978279"
e[2] = big"-0.0431154147693765556786027719094601654492921814596937312134031119112531513440588153185568190763617489229387640388525470783702009674718741042530929470636829615"
e[3] = big"-0.301177734439979327067358132631475024908856612462075071661429433362113545599206315692850620370071719845567311214260856007516359780524163248261047915189592181"
e[4] = big"-0.152771711192720814532185222660572424298879489997094076620761921659113290348236035957988677484677250113716116249595661396377139314351659646376567165252091199"
e[5] = big"-0.246155816769719505509940312338612337785792811109326830573144303824690959223435503683715228237563474852776043407904088603594853209284614752559595389755490267"
e[6] = big"-0.0794180284601028524502179548802248076178607044829584655895971993642306464734826800923021367047713399664978080601493771667070406823265242401710433650128072783"
e[7] = big"-0.0363933725938825618683946400054160074125285023799690190921600209776133289723506736807240580451513859987883185020494128433220917384498561166906858467063724684"

TI = Matrix{T1}(undef, 7, 7)
TI[1, 1] = big"258.131926319982229276108947425184471333411128774462923076434633414645220927977539758484670571338176678808837829326061674950321562391576244286310404028770676"
TI[1, 2] = big"189.073763081398508951976143411165126555759459745371576264125287430947886413126866952443113984840310549596923934762141954737541643761162558070450614795561734"
Expand Down Expand Up @@ -505,7 +526,7 @@ function BigRadauIIA13Tableau(T1, T2)
T[7, 7] = big"0.0"

RadauIIATableau{T1, T2}(T, TI,
c, γ, α, β, 7)
c, γ, α, β, #=e,=# 7)
end

using Polynomials, LinearAlgebra, GenericSchur, RootedTrees, Symbolics
Expand All @@ -524,7 +545,7 @@ function adaptiveRadauTableau(T1, T2, num_stages::Int)
for i in 1:(num_stages - 1)
radau_p = derivative(radau_p)
end
c = roots(radau_p)
c = real(roots(radau_p))
c[num_stages] = 1
c_powers = Matrix{BigFloat}(undef, num_stages, num_stages)
for i in 1 : num_stages
Expand Down Expand Up @@ -576,7 +597,7 @@ function adaptiveRadauTableau(T1, T2, num_stages::Int)
end
end
TI = inv(T)
#=

p = num_stages
eb = variables(:b, 1:num_stages + 1)
@variables y
Expand All @@ -588,13 +609,11 @@ function adaptiveRadauTableau(T1, T2, num_stages::Int)
constraints = map(Iterators.flatten(RootedTreeIterator(i) for i in 1:p)) do t
residual_order_condition(t, RungeKuttaMethod(eA, eb, ec))
end
AA, bb, islinear = Symbolics.linear_expansion(substitute.(constraints, (eb[1]=>y,)), eb[2:end])
AA, bb, islinear = Symbolics.linear_expansion(substitute.(constraints, (eb[1]=>1/γ,)), eb[2:end])
AA = BigFloat.(map(unwrap, AA))
idxs = qr(AA', ColumnNorm()).p[1:num_stages]
@assert rank(AA[idxs, :]) == num_stages
@assert islinear
Symbolics.expand.((AA[idxs, :] \ -bb[idxs]) - b)=#
#e = b_hat - b
RadauIIATableau{T1, T2}(T, TI, c, γ, α, β, num_stages)
b_hat = Symbolics.expand.((AA[idxs, :] \ -bb[idxs]) - b)
#e = symbolic_to_float(b_hat - b)
RadauIIATableau{T1, T2}(T, TI, c, γ, α, β, #=e,=# num_stages)
end

0 comments on commit 9c15269

Please sign in to comment.