add adaptivity

Shreyas-Ekanathan · Sep 18, 2024 · fed9fc5 · fed9fc5
1 parent 533f9a8
commit fed9fc5
Show file tree

Hide file tree

Showing 2 changed files with 63 additions and 49 deletions.
diff --git a/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl b/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl
@@ -1354,7 +1354,7 @@ end
 @muladd function perform_step!(integrator, cache::AdaptiveRadauConstantCache,
     repeat_step = false)
     @unpack t, dt, uprev, u, f, p = integrator
-    @unpack T, TI, γ, α, β, c, #=e,=# num_stages = cache.tab
+    @unpack T, TI, γ, α, β, c, e, num_stages = cache.tab
     @unpack κ, cont = cache
     @unpack internalnorm, abstol, reltol, adaptive = integrator.opts
     alg = unwrap_alg(integrator, true)
@@ -1393,9 +1393,9 @@ end
     if integrator.iter == 1 || integrator.u_modified || alg.extrapolant == :constant
         cache.dtprev = one(cache.dtprev)
         for i in 1 : num_stages
-            z[i] = map(zero, u)
-            w[i] = map(zero, u)
-            cache.cont[i] = map(zero, u)
+            z[i] = @.. map(zero, u)
+            w[i] = @.. map(zero, u)
+            cache.cont[i] = @.. map(zero, u)
         end
     else
         c_prime = Vector{typeof(u)}(undef, num_stages) #time stepping
@@ -1416,7 +1416,7 @@ end
         for i in 1:num_stages
             w[i] = @.. zero(u)
             for j in 1:num_stages
-                w[i] =@.. w[i] + TI[i,j] * z[j]
+                w[i] = @.. w[i] + TI[i,j] * z[j]
             end
         end
     end
@@ -1537,7 +1537,7 @@ end
         edt = e ./ dt
         tmp = dot(edt, z)
         mass_matrix != I && (tmp = mass_matrix * tmp)
-        utilde = @.. broadcast=false integrator.fsalfirst+tmp
+        utilde = @.. broadcast=false 1 / γ * dt * integrator.fsalfirst+tmp
         if alg.smooth_est
             utilde = _reshape(LU1 \ _vec(utilde), axes(u))
             integrator.stats.nsolve += 1
@@ -1549,7 +1549,7 @@ end
         integrator.u_modified
             f0 = f(uprev .+ utilde, p, t)
             integrator.stats.nf += 1
-            utilde = @.. broadcast=false f0+tmp
+            utilde = @.. broadcast=false 1 / γ * dt * f0 + tmp
             if alg.smooth_est
                 utilde = _reshape(LU1 \ _vec(utilde), axes(u))
                 integrator.stats.nsolve += 1
@@ -1574,7 +1574,7 @@ end
                 end
             end
             for i in 1 : num_stages
-                cache.cont[i] = derivatives[i, num_stages]
+                cache.cont[i] = @.. derivatives[i, num_stages]
             end
         end
     end
@@ -1589,7 +1589,7 @@ 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 T, TI, γ, α, β, c, e, num_stages = cache.tab
     @unpack κ, cont, derivatives, z, w, c_prime = cache
     @unpack dw1, ubuff, dw2, cubuff, dw = cache
     @unpack ks, k, fw, J, W1, W2 = cache
@@ -1738,9 +1738,9 @@ end
             end
         end
 
-        w[1] = @.. w[1] -  dw1
+        @.. w[1] = w[1] -  dw1
         for i in 2 : num_stages
-            w[i] = @.. w[i] - dw[i - 1]
+            @.. w[i] = w[i] - dw[i - 1]
         end
 
         # transform `w` to `z`
@@ -1784,7 +1784,7 @@ end
     if adaptive
         utilde = w2
         edt = e./dt
-        @.. tmp= dot(edt, z)
+        @.. tmp = dot(edt, z) + 1 / γ * dt * fsalfirst
         mass_matrix != I && (mul!(w1, mass_matrix, tmp); copyto!(tmp, w1))
         @.. ubuff=integrator.fsalfirst + tmp
 

diff --git a/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl b/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl
@@ -118,7 +118,7 @@ struct RadauIIATableau{T1, T2}
     γ::T1
     α::Vector{T1}
     β::Vector{T1}
-    #e::Vector{T1}
+    e::Vector{T1}
     num_stages::Int
 end
 
@@ -135,9 +135,9 @@ function BigRadauIIA5Tableau(T1, T2)
     c[3] = big"1"
 
     e = Vector{T1}(undef, 3)
-    e[1] = big"-0.804701356815835379608495496358640916569322134539215617920280276511680200030933806355291481868922518805459899199875734619185214695254668403298825163805293365"
-    e[2] = big"-0.267446751803505087778945794929857825182629030352446373700645445786652166171126710221557624849436998601914181921833023811045195644454184323577861370786198096"
-    e[3] = big"-0.202740720976336900360387312310916037542642249112497662477129527815832849767696924955091220025846425370353707227521045881781871836521093030105009933962412198"
+    e[1] = big"-0.428298294115368098113417591057340987284723986878723769598436582629558514867595819404541110575367601847354683220540647741034052880983125451920841851066713815"
+    e[2] = big"0.245039074384916438547779982042524963814308131639809054920681599475247366387530452311400517264433163448821319862243292031101333835037542080594323438762615605"
+    e[3] = big"-0.0916296098652259003910911359018870983677439465358993542342383692664256286871842771687905889462502859518430848371143253189423645209875225522045944109790661043"
 
     TI = Matrix{T1}(undef, 3, 3)
     TI[1, 1] = big"4.32557989006315535102435095295614882731995158490590784287320458848019483341979047442263696495019938973156007686663488090615420049217658854859024016717169837"
@@ -161,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, γ, α, β, #=e,=# 3)
+    c, γ, α, β, e, 3)
 end
 
 function BigRadauIIA9Tableau(T1, T2)
@@ -181,11 +181,11 @@ function BigRadauIIA9Tableau(T1, T2)
     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"
+    e[1] = big"-0.239909571163200476817707991076962793618683493830916562279975225042872448414819259070978815977101189851237591634144816820425592764432710089981892023441549743"
+    e[2] = big"0.125293484229223300606887443525929336197638450194973243323835481816625682669684634271160851996902445722310139152840852195000603355301430153561918341655137084"
+    e[3] = big"-0.0688288849083782089370741375422279772873399871312158026536514369967825732836040693366396751988019623495452730460577336089848105495733304937016477629621990433"
+    e[4] = big"0.0372433600301293198267284480468961585078575499935477902539140092558239369583143189611737274891328175087350382605569395937945872776839987254373869550943049195"
+    e[5] = big"-0.012863950751139890646895902730137465239579845437088427263654957605488133042638147254426913683751171160691603649073170415735165315443538347036196755548109703"
 
     TI = Matrix{T1}(undef, 5, 5)
     TI[1, 1] = big"30.0415677215444016277146611632467970747634862837368422955138470463852339244593400023985957753164599415374157317627305099177616927640413043608408838747985125"
@@ -242,7 +242,7 @@ function BigRadauIIA9Tableau(T1, T2)
     T[5, 5] = big"0.0"
 
     RadauIIATableau{T1, T2}(T, TI,  
-    c, γ, α, β, #=e,=# 5)
+    c, γ, α, β, e, 5)
 end
 
 
@@ -415,13 +415,13 @@ function BigRadauIIA13Tableau(T1, T2)
     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"
+    e[1] = big"-0.171003707892600662399316289094713451418682228802554671094543075316220821099884263713032871578607576486535539488632407200766379971279557791263375813632287147"
+    e[2] = big"0.0934967172358652400317534533028674569308657324394316331629203486361371292312231403973668315582870547753526899857449840409175610009437530537219068836721686211"
+    e[3] = big"-0.0538908303114758775848180855518003793385120454908028879947132475960997563222416509683480350817114056356343433378213334826358824351525243402758389616051172681"
+    e[4] = big"0.03036786965048439581219923157368590250090409822952169396779792168990510618756404452728392288892998298088147691907128776816545685599760715439221674418662785"
+    e[5] = big"-0.0169792974425458224044481617230998766694942329759644144506911809530904808476835739189122151558601434810772520378036293579816345384682687602578758514350075723"
+    e[6] = big"0.00942688256820236884916415666439281573527695349338346787075550606528435808025071461023926432308932314282041885090975780812273515256740094297311541275151861292"
+    e[7] = big"-0.00331409873565629283448601269346047459594635696619041493081994712789731442072563377354737487903843138987115421498455722938021358621090485566426506726181005806"
 
     TI = Matrix{T1}(undef, 7, 7)
     TI[1, 1] = big"258.131926319982229276108947425184471333411128774462923076434633414645220927977539758484670571338176678808837829326061674950321562391576244286310404028770676"
@@ -526,7 +526,7 @@ function BigRadauIIA13Tableau(T1, T2)
     T[7, 7] = big"0.0"
 
     RadauIIATableau{T1, T2}(T, TI,  
-    c, γ, α, β, #=e,=# 7)
+    c, γ, α, β, e, 7)
 end
 
 using Polynomials, LinearAlgebra, GenericSchur, RootedTrees, Symbolics
@@ -597,23 +597,37 @@ function adaptiveRadauTableau(T1, T2, num_stages::Int)
         end
     end
     TI = inv(T)
-
-    p = num_stages
-    eb = variables(:b, 1:num_stages + 1)
-    @variables y
-    zz = zeros(size(a, 1) + 1)
-    zz2 = zeros(size(a, 1))
-    eA = [zz'
-          zz2 a]
-    ec = [0; c]
-    constraints = map(Iterators.flatten(RootedTreeIterator(i) for i in 1:p)) do t
-        residual_order_condition(t, RungeKuttaMethod(eA, eb, ec))
+
+    if (num_stages == 9)
+        e = Vector{BigFloat}(undef, 9)
+        e[1] = big"-0.133101731359431287515066981129913748644705107621439651956220186897253838380345034218235538734967567153163030284540660584311040323114847240173627907922903296"
+        e[2] = big"0.0754476228408557299650196603226967248368445025181771896522057250989188754588885465998346476425502117889420021664297319179240040109156780754680742172762707621"
+        e[3] = big"-0.0458369394236156144604575482137179697005739995740615341890112217655441769701945378217626766299683076189687755618065050383493055018324395934911567207485032988"
+        e[4] = big"0.0271430329153098694457979735602502142083095152399102869109830450899844979409229538982100527256348792152825816553434603418662939944133319974874915933773657075"
+        e[5] = big"-0.0156126300301219212217568535995825232086423550686814635293876744035364259647929167763641353639085929285192729729570945658304937255929114458885296622493040224"
+        e[6] = big"0.00890598154557403928205152521539967562877335780940124672915181111908317890891659158654221736499522823959933517986673010006749138291836676520080172845444352328"
+        e[7] = big"-0.00514824122639241252178399021479378841872099572255461304439292434131750195489022869965968028106854978547414579491205935930595041763060069987112580994637398395"
+        e[8] = big"0.00296533914055503317169967748114188676589522458557982039693426239853498956125735811263087631479968309978854200615027412311940897061471388689986239742919640848"
+        e[9] = big"-0.0010634368308888065260482548541946175520274736959410047497431569257848032902381738362547705844630238841535652230832162703806430112125115777122361837311714267"
+    else
+        p = num_stages
+        eb = variables(:b, 1:num_stages + 1)
+        @variables y
+        zz = zeros(size(a, 1) + 1)
+        zz2 = zeros(size(a, 1))
+        eA = [zz'
+            zz2 a]
+        ec = [0; c]
+        constraints = map(Iterators.flatten(RootedTreeIterator(i) for i in 1:2*p-3)) do t
+            residual_order_condition(t, RungeKuttaMethod(eA, eb, ec))
+        end 
+        AA, bb, islinear = Symbolics.linear_expansion(Symbolics.substitute.(constraints, (eb[1]=>1/γ,)), eb[2:end])
+        AA = Float64.(map(unwrap, AA))
+        idxs = qr(AA', ColumnNorm()).p[1:num_stages]
+        @assert rank(AA[idxs, :]) == num_stages
+        @assert islinear
+        b_hat = Symbolics.expand.((AA \ -bb))
+        e = [Symbolics.symbolic_to_float(b_hat[i] - b[i]) for i in 1 : num_stages]    
     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 islinear
-    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)
+    RadauIIATableau{T1, T2}(T, TI, c, γ, α, β, e, num_stages)
 end