WIP: Add BHR 553 methods by nickki-m · Pull Request #142 · NumericalMathematics/Ariadne.jl

nickki-m · 2026-04-28T08:09:38Z

Adding methods from https://doi.org/10.1137/080713562.

ranocha · 2026-05-05T08:05:50Z

julia> alg = Theseus.BHR553G1()
Theseus.BHR553G1()

julia> tab = Theseus.RKTableau(alg, Float64)
Theseus.IMEXButcher{Matrix{Float64}, Vector{Float64}}([0.0 0.0 … 0.0 0.0; 0.8717330430169644 0.0 … 0.0 0.0; … ; -0.8009984530656288 0.0 … 0.0 0.0; 0.35675320777963954 -0.19733989037886918 … -0.04136215879462423 0.0], [0.4128980428124743, 0.0, 0.0, 0.19733989037886918, 0.4358665215084822], [0.0, 0.8717330430169644, 0.8717330430168252, 2.3402125751086804, 1.0], [0.0 0.0 … 0.0 0.0; 0.4358665215084822 0.4358665215084822 … 0.0 0.0; … ; -0.06675868701958376 -6.129659659417483e-13 … 0.4358665215084822 0.0; 0.4128980428124743 0.0 … 0.19733989037886918 0.4358665215084822], [0.4128980428124743, 0.0, 0.19733989037886918, -0.04610445469982568, 0.4358665215084822], [0.0, 0.8717330430169644, 0.8717330430168252, 2.3402125751086804, 1.0])

julia> tab.b_ex' * tab.c_ex
0.8976838145436803

…some tests

nickki-m · 2026-05-07T17:54:22Z

concerning the BHR553 methods: i changed the bs of the explicit parts to the same b as the implicit part, i also changed the A-matrix of the implicit parts so that the methods are stiffly accurate. In the last testset it seems that some of the individual methods are of higher order than the one given in the paper.
concerning the ESDIRKSA2-method: in some tests i changed the newton-tolerances to 1e-12

into BHR_IMEX_method

MarcoArtiano · 2026-05-08T07:47:28Z

LGTM!

ranocha

Thanks

github-actions · 2026-05-08T13:03:35Z

Your PR requires formatting changes to meet the project's style guidelines.
Please consider running Runic (git runic main) to apply these changes.

Click here to view the suggested changes.

diff --git a/libs/Theseus/src/dirk/tableau.jl b/libs/Theseus/src/dirk/tableau.jl
index bcf50eb..2f4e4b2 100644
--- a/libs/Theseus/src/dirk/tableau.jl
+++ b/libs/Theseus/src/dirk/tableau.jl
@@ -140,7 +140,7 @@ function RKTableau(alg::ESDIRK43SA2, RealT)
 
     b = a[6, :]
     #weights for the embedded method:
-    #b_e = zeros(RealT, nstage) 
+    #b_e = zeros(RealT, nstage)
     #b_e[1] = -1007911106287 // 12117826057527
     #b_e[2] = b_e[1]
     #b_e[3] = 17694008993113 // 35931961998873
diff --git a/libs/Theseus/src/imex/tableau.jl b/libs/Theseus/src/imex/tableau.jl
index edf391a..9c30705 100644
--- a/libs/Theseus/src/imex/tableau.jl
+++ b/libs/Theseus/src/imex/tableau.jl
@@ -549,36 +549,36 @@ function RKTableau(alg::BHR553G1, RealT)
     nstage = 5
     gamma = 424782 // 974569
     a = zeros(RealT, nstage, nstage)
-    a[2, 1] = 2*gamma
+    a[2, 1] = 2 * gamma
     a[3, 1] = gamma
     a[3, 2] = gamma
     a[4, 1] = -475883375220285986033264 // 594112726933437845704163
-    a[4, 3] = 1866233449822026827708736//594112726933437845704163
-    a[5, 1] = 62828845818073169585635881686091391737610308247 //176112910684412105319781630311686343715753056000
-    a[5, 2] = -302987763081184622639300143137943089 //1535359944203293318639180129368156500
+    a[4, 3] = 1866233449822026827708736 // 594112726933437845704163
+    a[5, 1] = 62828845818073169585635881686091391737610308247 // 176112910684412105319781630311686343715753056000
+    a[5, 2] = -302987763081184622639300143137943089 // 1535359944203293318639180129368156500
     a[5, 3] = 262315887293043739337088563996093207 // 297427554730376353252081786906492000
     a[5, 4] = -987618231894176581438124717087 // 23877337660202969319526901856000
     b = zeros(RealT, nstage)
-    b[1] = 487698502336740678603511//1181159636928185920260208
-    b[3] = 302987763081184622639300143137943089//1535359944203293318639180129368156500
-    b[4] = -105235928335100616072938218863//2282554452064661756575727198000
+    b[1] = 487698502336740678603511 // 1181159636928185920260208
+    b[3] = 302987763081184622639300143137943089 // 1535359944203293318639180129368156500
+    b[4] = -105235928335100616072938218863 // 2282554452064661756575727198000
     #
     b[5] = gamma
     c = zeros(RealT, nstage)
-    c[2] = 2* gamma
-    c[3] = 902905985686//1035759735069
-    c[4] = 2684624//1147171
+    c[2] = 2 * gamma
+    c[3] = 902905985686 // 1035759735069
+    c[4] = 2684624 // 1147171
     c[5] = 1
     a_im = zeros(RealT, nstage, nstage)
     a_im[2, 1] = gamma
     a_im[2, 2] = gamma
     a_im[3, 1] = gamma
-    a_im[4, 1] = -3012378541084922027361996761794919360516301377809610//45123394056585269977907753045030512597955897345819349
-    a_im[5, 1] = b[1]   
-    a_im[3, 2] = -31733082319927313//455705377221960889379854647102
+    a_im[4, 1] = -3012378541084922027361996761794919360516301377809610 // 45123394056585269977907753045030512597955897345819349
+    a_im[5, 1] = b[1]
+    a_im[3, 2] = -31733082319927313 // 455705377221960889379854647102
     a_im[3, 3] = gamma
-    a_im[4, 2] = -62865589297807153294268//102559673441610672305587327019095047
-    a_im[4, 3] = 418769796920855299603146267001414900945214277000//212454360385257708555954598099874818603217167139
+    a_im[4, 2] = -62865589297807153294268 // 102559673441610672305587327019095047
+    a_im[4, 3] = 418769796920855299603146267001414900945214277000 // 212454360385257708555954598099874818603217167139
     a_im[4, 4] = gamma
     a_im[5, 3] = b[3]
     a_im[5, 4] = b[4]
@@ -612,37 +612,37 @@ struct BHR553G2 <: RKIMEX{5} end
 function RKTableau(alg::BHR553G2, RealT)
     # BHR(5,5,3)_g2 IMEX Runge-Kutta - Third order
     nstage = 5
-    gamma2= 2051948 // 3582211
+    gamma2 = 2051948 // 3582211
     a = zeros(RealT, nstage, nstage)
-    a[2, 1] = 2*gamma2
-    a[3, 1] = 473447115440655855452482357894373//1226306256343706154920072735579148
-    a[4, 1] = 37498105210828143724516848//172642583546398006173766007
-    a[5, 1] =  -3409975860212064612303539855622639333030782744869519//5886704102363745137792385361113084313351870216475136
-    a[3, 2] = 129298766034131882323069978722019//1226306256343706154920072735579148
-    a[5, 2] =  -237416352433826978856941795734073//554681702576878342891447163499456
-    a[4, 3] = 76283359742561480140804416//172642583546398006173766007
-    a[5, 3] = 4298159710546228783638212411650783228275//2165398513352098924587211488610407046208
-    a[5, 4] = 6101865615855760853571922289749//272863973025878249803640374568448
+    a[2, 1] = 2 * gamma2
+    a[3, 1] = 473447115440655855452482357894373 // 1226306256343706154920072735579148
+    a[4, 1] = 37498105210828143724516848 // 172642583546398006173766007
+    a[5, 1] = -3409975860212064612303539855622639333030782744869519 // 5886704102363745137792385361113084313351870216475136
+    a[3, 2] = 129298766034131882323069978722019 // 1226306256343706154920072735579148
+    a[5, 2] = -237416352433826978856941795734073 // 554681702576878342891447163499456
+    a[4, 3] = 76283359742561480140804416 // 172642583546398006173766007
+    a[5, 3] = 4298159710546228783638212411650783228275 // 2165398513352098924587211488610407046208
+    a[5, 4] = 6101865615855760853571922289749 // 272863973025878249803640374568448
     b = zeros(RealT, nstage)
-    b[1] =  -2032971420760927701493589//38017147656515384190997416
-    b[3] = 2197602776651676983265261109643897073447//945067123279139583549933947379097184164
-    b[4] = -128147215194260398070666826235339//69468482710687503388562952626424
+    b[1] = -2032971420760927701493589 // 38017147656515384190997416
+    b[3] = 2197602776651676983265261109643897073447 // 945067123279139583549933947379097184164
+    b[4] = -128147215194260398070666826235339 // 69468482710687503388562952626424
     b[5] = gamma2
     c = zeros(RealT, nstage)
-    c[2] = 2* gamma2
-    c[3] = 12015769930846//24446477850549
-    c[4] = 3532944//5360597
+    c[2] = 2 * gamma2
+    c[3] = 12015769930846 // 24446477850549
+    c[4] = 3532944 // 5360597
     c[5] = 1
     a_im = zeros(RealT, nstage, nstage)
     a_im[2, 1] = gamma2
     a_im[2, 2] = gamma2
-    a_im[3, 1] = 259252258169672523902708425780469319755//4392887760843243968922388674191715336228
-    a_im[4, 1] = 1103202061574553405285863729195740268785131739395559693754//9879457735937277070641522414590493459028264677925767305837
-    a_im[5, 1] = b[1]   
-    a_im[3, 2] = -172074174703261986564706189586177//1226306256343706154920072735579148
+    a_im[3, 1] = 259252258169672523902708425780469319755 // 4392887760843243968922388674191715336228
+    a_im[4, 1] = 1103202061574553405285863729195740268785131739395559693754 // 9879457735937277070641522414590493459028264677925767305837
+    a_im[5, 1] = b[1]
+    a_im[3, 2] = -172074174703261986564706189586177 // 1226306256343706154920072735579148
     a_im[3, 3] = gamma2
-    a_im[4, 2] =  -103754520567058969566542556296087324094//459050363888246734833121482275319954529
-    a_im[4, 3] = 3863207083069979654596872190377240608602701071947128//19258690251287609765240683320611425745736762681950551
+    a_im[4, 2] = -103754520567058969566542556296087324094 // 459050363888246734833121482275319954529
+    a_im[4, 3] = 3863207083069979654596872190377240608602701071947128 // 19258690251287609765240683320611425745736762681950551
     a_im[4, 4] = gamma2
     a_im[5, 3] = b[3]
     a_im[5, 4] = b[4]
diff --git a/libs/Theseus/test/convergence.jl b/libs/Theseus/test/convergence.jl
index fd39364..5a1cd86 100644
--- a/libs/Theseus/test/convergence.jl
+++ b/libs/Theseus/test/convergence.jl
@@ -6,7 +6,7 @@ function compute_eoc(dts, errors)
     eocs = similar(errors)
     eocs[begin] = 0 # no EOC defined for the first grid
     for idx in Iterators.drop(eachindex(errors, dts, eocs), 1)
-        eocs[idx] = log(errors[idx] / errors[idx - 1]) / 
+        eocs[idx] = log(errors[idx] / errors[idx - 1]) /
             log(dts[idx] / dts[idx - 1])
     end
     return mean(eocs[(begin + 1):end])
@@ -16,7 +16,7 @@ function compute_errors(prob, u_ana, alg, dts; kwargs...)
     errors = similar(dts)
     for (i, dt) in enumerate(dts)
         sol = @inferred solve(
-            prob, alg; 
+            prob, alg;
             dt, adaptive = false, kwargs...
         )
         u_num = sol.u[end]
@@ -216,7 +216,7 @@ end
             eoc = compute_eoc(dts, errors)
             @test isapprox(eoc, order; atol = 0.1)
         end
-  
+
         @testset "KenCarpARK437" begin
             alg = Theseus.KenCarpARK437()
             order = 4
@@ -234,7 +234,7 @@ end
             eoc = compute_eoc(dts, errors)
             @test isapprox(eoc, order; atol = 0.1)
         end
-  
+
         @testset "KenCarpARK548" begin
             alg = Theseus.KenCarpARK548()
             order = 5
@@ -745,7 +745,7 @@ end
         @test isapprox(eoc, order; atol = 0.1)
     end
 
-     @testset "BHR553G1" begin
+    @testset "BHR553G1" begin
         alg = Theseus.BHR553G1()
         order = 3
         dts = 2.0 .^ (-3:-1:-7)
@@ -756,9 +756,9 @@ end
         )
         eoc = compute_eoc(dts, errors)
         @test isapprox(eoc, order; atol = 0.1)
-     end
+    end
 
-     @testset "BHR553G2" begin
+    @testset "BHR553G2" begin
         alg = Theseus.BHR553G2()
         order = 3
         dts = 2.0 .^ (-3:-1:-7)
@@ -769,5 +769,5 @@ end
         )
         eoc = compute_eoc(dts, errors)
         @test isapprox(eoc, order; atol = 0.1)
-     end
+    end
 end

codecov-commenter · 2026-05-08T13:03:56Z

Codecov Report

✅ All modified and coverable lines are covered by tests.

📢 Thoughts on this report? Let us know!

ranocha · 2026-05-08T13:06:05Z

Could you please check the failing convergence tests?

ranocha · 2026-05-12T07:49:18Z

+            errors = compute_errors(ode_split, u_ana, alg, dts; newton_tol_abs = 1.0e-9, newton_tol_rel = 1.0e-9)
            eoc = compute_eoc(dts, errors)
            @test isapprox(eoc, order; atol = 0.1)
        end


Please add convergence tests in the new test sets at the end of the file.

ranocha · 2026-05-12T07:49:46Z

+    @testset "BHR553G1" begin
+        alg = Theseus.BHR553G1()
+        order = 3
+        dts = 2.0 .^ (-2:-1:-6)
+        for ode_split in (ode_split_1, ode_split_2)
+            errors = compute_errors(ode_split, u_ana, alg, dts; newton_tol_abs = 1.0e-9, newton_tol_rel = 1.0e-9)
+            eoc = compute_eoc(dts, errors)
+            @test isapprox(eoc, order; atol = 0.1)
+        end
+    end


Explanation:

julia> alg = Theseus.BHR553G1() Theseus.BHR553G1() julia> tab = Theseus.RKTableau(alg, Float64) Theseus.IMEXButcher{Matrix{Float64}, Vector{Float64}}([0.0 0.0 … 0.0 0.0; 0.8717330430169644 0.0 … 0.0 0.0; … ; -0.8009984530656288 0.0 … 0.0 0.0; 0.35675320777963954 -0.19733989037886918 … -0.04136215879462423 0.0], [0.4128980428124743, 0.0, 0.19733989037886918, -0.04610445469982568, 0.4358665215084822], [0.0, 0.8717330430169644, 0.8717330430168252, 2.3402125751086804, 1.0], [0.0 0.0 … 0.0 0.0; 0.4358665215084822 0.4358665215084822 … 0.0 0.0; … ; -0.06675868701958376 -6.129659659417483e-13 … 0.4358665215084822 0.0; 0.4128980428124743 0.0 … -0.04610445469982568 0.4358665215084822], [0.4128980428124743, 0.0, 0.19733989037886918, -0.04610445469982568, 0.4358665215084822], [0.0, 0.8717330430169644, 0.8717330430168252, 2.3402125751086804, 1.0]) julia> tab.b_ex' * tab.c_ex 0.5 julia> tab.b_ex' * tab.a_ex * tab.c_ex 0.16666666666666669 julia> tab.b_ex' * tab.a_ex^2 * tab.c_ex 0.04166666666669573 julia> 1/24 0.041666666666666664

Suggested change

@testset "BHR553G1" begin

alg = Theseus.BHR553G1()

order = 3

dts = 2.0 .^ (-2:-1:-6)

for ode_split in (ode_split_1, ode_split_2)

errors = compute_errors(ode_split, u_ana, alg, dts; newton_tol_abs = 1.0e-9, newton_tol_rel = 1.0e-9)

eoc = compute_eoc(dts, errors)

@test isapprox(eoc, order; atol = 0.1)

end

end

@testset "BHR553G1" begin

alg = Theseus.BHR553G1()

order = 3

dts = 2.0 .^ (-2:-1:-6)

let ode_split = ode_split_1

# The explicit method is fourth-order accurate for linear problems

errors = compute_errors(ode_split, u_ana, alg, dts; newton_tol_abs = 1.0e-9, newton_tol_rel = 1.0e-9)

eoc = compute_eoc(dts, errors)

@test isapprox(eoc, order + 1; atol = 0.1)

end

let ode_split = ode_split_2

errors = compute_errors(ode_split, u_ana, alg, dts; newton_tol_abs = 1.0e-9, newton_tol_rel = 1.0e-9)

eoc = compute_eoc(dts, errors)

@test isapprox(eoc, order; atol = 0.1)

end

end

ranocha · 2026-05-12T07:52:07Z

+
+    @testset "BHR553G2" begin
+        alg = Theseus.BHR553G2()
+        order = 3


Suggested change

order = 3

# Both individual methods are fourth-order accurate for linear problems.

order = 3 + 1

Check:

julia> alg = Theseus.BHR553G2() Theseus.BHR553G2() julia> tab = Theseus.RKTableau(alg, Float64); julia> tab.b_ex' * tab.c_ex 0.4999999999999999 julia> tab.b_ex' * tab.a_ex * tab.c_ex 0.16666666666666657 julia> tab.b_ex' * tab.a_ex^2 * tab.c_ex 0.041666666666666664 julia> tab.b_im' * tab.a_im * tab.c_im 0.16666666666666685 julia> tab.b_im' * tab.a_im^2 * tab.c_im 0.041666666666644536 julia> 1/24 0.041666666666666664

ranocha · 2026-05-12T07:53:11Z

+    @testset "BHR553G2" begin
+        alg = Theseus.BHR553G2()
+        order = 3
+        dts = 2.0 .^ (-2:-1:-6)


Maybe

Suggested change

dts = 2.0 .^ (-2:-1:-6)

dts = 2.0 .^ (-3:-1:-7)

or increase the tolarance or make the newton tolerances stricter.

…rgence tests

into BHR_IMEX_method

add Tableau BHR553

7e9eb01

ranocha reviewed May 1, 2026

View reviewed changes

Comment thread libs/Theseus/src/imex/tableau.jl

MarcoArtiano linked an issue May 2, 2026 that may be closed by this pull request

Add more IMEX and DIRK methods #124

Open

13 tasks

i added the second BHR553 method and some tests

ae08d7f

nickki-m and others added 2 commits May 7, 2026 19:40

I added the ESDIRKSA2 method, corrected the BHR553 methods and added …

934f30d

…some tests

Merge branch 'main' into BHR_IMEX_method

cb41538

ranocha reviewed May 7, 2026

View reviewed changes

Comment thread libs/Theseus/src/dirk/tableau.jl Outdated

Comment thread libs/Theseus/src/imex/tableau.jl Outdated

Comment thread libs/Theseus/src/imex/tableau.jl Outdated

nickki-m added 2 commits May 7, 2026 20:47

added stability properties of the BHR553 methods

c686d93

Merge branch 'BHR_IMEX_method' of https://github.com/nickki-m/Ariadne.jl

de7892f

into BHR_IMEX_method

ranocha reviewed May 8, 2026

View reviewed changes

Comment thread libs/Theseus/src/imex/tableau.jl Outdated

Apply suggestion from @ranocha

5259aed

ranocha reviewed May 8, 2026

View reviewed changes

Comment thread libs/Theseus/src/imex/tableau.jl Outdated

Apply suggestion from @ranocha

fd11526

ranocha approved these changes May 8, 2026

View reviewed changes

ranocha removed a link to an issue May 8, 2026

Add more IMEX and DIRK methods #124

Open

13 tasks

ranocha linked an issue May 8, 2026 that may be closed by this pull request

Add DIRK methods of Kennedy and Carpenter (2019) #66

Open

ranocha mentioned this pull request May 8, 2026

Add more IMEX and DIRK methods #124

Open

13 tasks

ranocha marked this pull request as ready for review May 8, 2026 13:02

ranocha requested changes May 12, 2026

View reviewed changes

nickki-m added 2 commits May 12, 2026 18:09

added BHR -emothds to the latest testsets and fixed the failing conve…

3a46cf8

…rgence tests

Merge branch 'BHR_IMEX_method' of https://github.com/nickki-m/Ariadne.jl

e7d5601

into BHR_IMEX_method

	order = 3
	# Both individual methods are fourth-order accurate for linear problems.
	order = 3 + 1

Conversation

nickki-m commented Apr 28, 2026

Uh oh!

Uh oh!

ranocha commented May 5, 2026

Uh oh!

nickki-m commented May 7, 2026

Uh oh!

Uh oh!

Uh oh!

Uh oh!

MarcoArtiano commented May 8, 2026

Uh oh!

Uh oh!

Uh oh!

ranocha left a comment

Choose a reason for hiding this comment

Uh oh!

github-actions Bot commented May 8, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

codecov-commenter commented May 8, 2026

Codecov Report

Uh oh!

ranocha commented May 8, 2026

Uh oh!

ranocha May 12, 2026

Choose a reason for hiding this comment

Uh oh!

ranocha May 12, 2026

Choose a reason for hiding this comment

Uh oh!

ranocha May 12, 2026

Choose a reason for hiding this comment

Uh oh!

ranocha May 12, 2026

Choose a reason for hiding this comment

Uh oh!

ranocha May 12, 2026

Choose a reason for hiding this comment

Uh oh!

ranocha May 12, 2026

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

github-actions Bot commented May 8, 2026 •

edited

Loading