From 9c15269eeb5bfc75f69769ec5926f0d771191116 Mon Sep 17 00:00:00 2001 From: Shreyas Ekanathan Date: Wed, 11 Sep 2024 21:01:17 -0400 Subject: [PATCH] caches --- lib/OrdinaryDiffEqFIRK/src/firk_caches.jl | 7 +++- .../src/firk_perform_step.jl | 18 ++++---- lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl | 41 ++++++++++++++----- 3 files changed, 45 insertions(+), 21 deletions(-) diff --git a/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl b/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl index 642ebe9b67..f282953419 100644 --- a/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl +++ b/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl @@ -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 @@ -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 @@ -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, diff --git a/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl b/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl index 6d666a738b..24e8c080ab 100644 --- a/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl +++ b/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl @@ -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] @@ -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 @@ -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 @@ -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] @@ -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] @@ -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 @@ -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 diff --git a/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl b/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl index c491c8db72..59382feeaa 100644 --- a/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl +++ b/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl @@ -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" @@ -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) @@ -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" @@ -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 @@ -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" @@ -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 @@ -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 @@ -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 @@ -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 -