Skip to content

Commit 830d51f

Browse files
ErikQQYErikQQY
committed
Add improved MIRK6 method
1 parent 845a5e3 commit 830d51f

File tree

6 files changed

+268
-1
lines changed

6 files changed

+268
-1
lines changed

‎lib/BoundaryValueDiffEqMIRK/src/BoundaryValueDiffEqMIRK.jl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -157,7 +157,7 @@ include("sparse_jacobians.jl")
157157
end
158158
end
159159

160-
export MIRK2, MIRK3, MIRK4, MIRK5, MIRK6
160+
export MIRK2, MIRK3, MIRK4, MIRK5, MIRK6, MIRK6I
161161
export BVPJacobianAlgorithm
162162

163163
end

‎lib/BoundaryValueDiffEqMIRK/src/adaptivity.jl‎

Lines changed: 116 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -377,3 +377,119 @@ for order in (2, 3, 4, 5, 6)
377377
end
378378
end
379379
end
380+
381+
for order in (6,)
382+
alg = Symbol("MIRK$(order)I")
383+
@eval begin
384+
function interp_weights(τ::T, ::$(alg)) where {T}
385+
if $(order == 6)
386+
w = [
387+
-(12233 + 1450 * sqrt(7)) *
388+
(800086000 * τ^5 + 63579600 * sqrt(7) * τ^4 - 2936650584 * τ^4 +
389+
4235152620 * τ^3 - 201404565 * sqrt(7) * τ^3 +
390+
232506630 * sqrt(7) * τ^2 - 3033109390 * τ^2 + 1116511695 * τ -
391+
116253315 * sqrt(7) * τ + 22707000 * sqrt(7) - 191568780) *
392+
τ / 2112984835740,
393+
-(-10799 + 650 * sqrt(7)) *
394+
(24962000 * τ^4 + 473200 * sqrt(7) * τ^3 - 67024328 * τ^3 -
395+
751855 * sqrt(7) * τ^2 + 66629600 * τ^2 - 29507250 * τ +
396+
236210 * sqrt(7) * τ +
397+
5080365 +
398+
50895 * sqrt(7)) *
399+
τ^2 / 29551834260,
400+
7 / 1274940 *
401+
(259 + 50 * sqrt(7)) *
402+
(14000 * τ^4 - 48216 * τ^3 + 1200 * sqrt(7) * τ^3 -
403+
3555 * sqrt(7) * τ^2 +
404+
62790 * τ^2 +
405+
3610 * sqrt(7) * τ - 37450 * τ + 9135 - 1305 * sqrt(7)) *
406+
τ^2,
407+
7 / 1274940 *
408+
(259 + 50 * sqrt(7)) *
409+
(14000 * τ^4 - 48216 * τ^3 + 1200 * sqrt(7) * τ^3 -
410+
3555 * sqrt(7) * τ^2 +
411+
62790 * τ^2 +
412+
3610 * sqrt(7) * τ - 37450 * τ + 9135 - 1305 * sqrt(7)) *
413+
τ^2,
414+
16 / 2231145 *
415+
(259 + 50 * sqrt(7)) *
416+
(14000 * τ^4 - 48216 * τ^3 + 1200 * sqrt(7) * τ^3 -
417+
3555 * sqrt(7) * τ^2 +
418+
62790 * τ^2 +
419+
3610 * sqrt(7) * τ - 37450 * τ + 9135 - 1305 * sqrt(7)) *
420+
τ^2,
421+
4 / 1227278493 *
422+
(740 * sqrt(7) - 6083) *
423+
(1561000 * τ^2 - 2461284 * τ - 109520 * sqrt(7) * τ +
424+
979272 +
425+
86913 * sqrt(7)) *
426+
(τ - 1)^2 *
427+
τ^2,
428+
-49 / 63747 *
429+
sqrt(7) *
430+
(20000 * τ^2 - 20000 * τ + 3393) *
431+
(τ - 1)^2 *
432+
τ^2,
433+
-1250000000 / 889206903 * (28 * τ^2 - 28 * τ + 9) * (τ - 1)^2 * τ^2]
434+
435+
# Derivative polynomials.
436+
437+
wp = [
438+
(1450 * sqrt(7) + 12233) *
439+
(14 * τ - 7 + sqrt(7)) *
440+
(τ - 1) *
441+
(-400043 * τ + 75481 + 2083 * sqrt(7)) *
442+
(100 * τ - 87) *
443+
(2 * τ - 1) / 493029795006,
444+
-(650 * sqrt(7) - 10799) *
445+
(14 * τ - 7 + sqrt(7)) *
446+
(37443 * τ - 13762 - 2083 * sqrt(7)) *
447+
(100 * τ - 87) *
448+
(2 * τ - 1) *
449+
τ / 20686283982,
450+
7 / 42498 *
451+
(259 + 50 * sqrt(7)) *
452+
(14 * τ - 7 + sqrt(7)) *
453+
(τ - 1) *
454+
(100 * τ - 87) *
455+
(2 * τ - 1) *
456+
τ,
457+
7 / 42498 *
458+
(259 + 50 * sqrt(7)) *
459+
(14 * τ - 7 + sqrt(7)) *
460+
(τ - 1) *
461+
(100 * τ - 87) *
462+
(2 * τ - 1) *
463+
τ,
464+
32 / 148743 *
465+
(259 + 50 * sqrt(7)) *
466+
(14 * τ - 7 + sqrt(7)) *
467+
(τ - 1) *
468+
(100 * τ - 87) *
469+
(2 * τ - 1) *
470+
τ,
471+
4 / 1227278493 *
472+
(740 * sqrt(7) - 6083) *
473+
(14 * τ - 7 + sqrt(7)) *
474+
(τ - 1) *
475+
(100 * τ - 87) *
476+
(6690 * τ - 4085 - 869 * sqrt(7)) *
477+
τ,
478+
-98 / 21249 *
479+
sqrt(7) *
480+
(τ - 1) *
481+
(100 * τ - 13) *
482+
(100 * τ - 87) *
483+
(2 * τ - 1) *
484+
τ,
485+
-1250000000 / 2074816107 *
486+
(14 * τ - 7 + sqrt(7)) *
487+
(τ - 1) *
488+
(14 * τ - 7 - sqrt(7)) *
489+
(2 * τ - 1) *
490+
τ]
491+
end
492+
return T.(w), T.(wp)
493+
end
494+
end
495+
end

‎lib/BoundaryValueDiffEqMIRK/src/alg_utils.jl‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -4,4 +4,10 @@ for order in (2, 3, 4, 5, 6)
44
@eval alg_stage(::$(alg)) = $(order - 1)
55
end
66

7+
for order in (6,)
8+
alg = Symbol("MIRK$(order)I")
9+
@eval alg_order(::$(alg)) = $order
10+
@eval alg_stage(::$(alg)) = $(order - 1)
11+
end
12+
713
SciMLBase.isadaptive(alg::AbstractMIRK) = true

‎lib/BoundaryValueDiffEqMIRK/src/algorithms.jl‎

Lines changed: 53 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -53,3 +53,56 @@ for order in (2, 3, 4, 5, 6)
5353
end
5454
end
5555
end
56+
57+
for order in (6)
58+
alg = Symbol("MIRK$(order)I")
59+
60+
@eval begin
61+
"""
62+
$($alg)(; nlsolve = NewtonRaphson(), jac_alg = BVPJacobianAlgorithm(),
63+
defect_threshold = 0.1, max_num_subintervals = 3000)
64+
65+
$($order)th order Monotonic Implicit Runge Kutta method.
66+
67+
## Keyword Arguments
68+
69+
- `nlsolve`: Internal Nonlinear solver. Any solver which conforms to the SciML
70+
`NonlinearProblem` interface can be used. Note that any autodiff argument for
71+
the solver will be ignored and a custom jacobian algorithm will be used.
72+
- `jac_alg`: Jacobian Algorithm used for the nonlinear solver. Defaults to
73+
`BVPJacobianAlgorithm()`, which automatically decides the best algorithm to
74+
use based on the input types and problem type.
75+
- For `TwoPointBVProblem`, only `diffmode` is used (defaults to
76+
`AutoSparse(AutoForwardDiff())` if possible else `AutoSparse(AutoFiniteDiff())`).
77+
- For `BVProblem`, `bc_diffmode` and `nonbc_diffmode` are used. For
78+
`nonbc_diffmode` defaults to `AutoSparse(AutoForwardDiff())` if possible else
79+
`AutoSparse(AutoFiniteDiff())`. For `bc_diffmode`, defaults to `AutoForwardDiff` if
80+
possible else `AutoFiniteDiff`.
81+
- `defect_threshold`: Threshold for defect control.
82+
- `max_num_subintervals`: Number of maximal subintervals, default as 3000.
83+
84+
!!! note
85+
For type-stability, the chunksizes for ForwardDiff ADTypes in
86+
`BVPJacobianAlgorithm` must be provided.
87+
88+
## References
89+
90+
```bibtex
91+
@article{Enright1996RungeKuttaSW,
92+
title={Runge-Kutta Software with Defect Control for Boundary Value ODEs},
93+
author={Wayne H. Enright and Paul H. Muir},
94+
journal={SIAM J. Sci. Comput.},
95+
year={1996},
96+
volume={17},
97+
pages={479-497}
98+
}
99+
```
100+
"""
101+
@kwdef struct $(alg){N, J <: BVPJacobianAlgorithm, T} <: AbstractMIRK
102+
nlsolve::N = nothing
103+
jac_alg::J = BVPJacobianAlgorithm()
104+
defect_threshold::T = 0.1
105+
max_num_subintervals::Int = 3000
106+
end
107+
end
108+
end

‎lib/BoundaryValueDiffEqMIRK/src/mirk_tableaus.jl‎

Lines changed: 44 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,12 @@ for order in (2, 3, 4, 5, 6)
44
@eval constructMIRK(::$(alg), ::Type{T}) where {T} = $(f)(T)
55
end
66

7+
for order in (6,)
8+
alg = Symbol("MIRK$(order)I")
9+
f = Symbol("constructMIRK$(order)I")
10+
@eval constructMIRK(::$(alg), ::Type{T}) where {T} = $(f)(T)
11+
end
12+
713
function constructMIRK2(::Type{T}) where {T}
814
# RK coefficients tableau
915
s = 1
@@ -117,3 +123,41 @@ function constructMIRK6(::Type{T}) where {T}
117123
ITU = MIRKInterpTableau(s_star, T.(c_star), T.(v_star), T.(x_star), T(τ_star))
118124
return TU, ITU
119125
end
126+
127+
function constructMIRK6I(::Type{T}) where {T}
128+
# RK coefficients tableau
129+
s = 5
130+
c = [0, 1, 1 // 2 - sqrt(21) / 14, 1 // 2 + sqrt(21) / 14, 1 // 2]
131+
v = [0, 1, 1 // 2 - 9 * sqrt(21) / 98, 1 // 2 + 9 * sqrt(21) / 98, 1 // 2]
132+
b = [7 // 90, 7 // 90, 16 // 45, 16 // 45, 2 // 15, 0, 0, 0, 0]
133+
b = [1 // 20, 1 // 20, 49 // 180, 49 // 180, 16 // 45]
134+
x = [0 0 0 0 0
135+
0 0 0 0 0
136+
9//64 -3//64 0 0 0
137+
3//64 -9//64 0 0 0
138+
-5//24 5//24 2//3 -2//3 0]
139+
140+
x = [0 0 0 0 0
141+
0 0 0 0 0
142+
1 // 14+sqrt(21) / 98 -1 // 14+sqrt(21) / 98 0 0 0
143+
1 // 14-sqrt(21) / 98 -1 // 14-sqrt(21) / 98 0 0 0
144+
-5//128 5//128 7 * sqrt(21)/128 -7 * sqrt(21)/128 0]
145+
146+
# Interpolant tableau
147+
s_star = 8
148+
c_star = [1 // 2, 1 // 2 - sqrt(7) / 14, 87 // 100]
149+
v_star = [1 // 2, 1 // 2 - sqrt(7) / 14, 87 // 100]
150+
x_star = [1547//32768 -1225//32768 749//4096 -287//2048 -861//16384 0 0 0 0
151+
83//1536 -13//384 283//1536 -167//1536 -49//512 0 0 0 0
152+
1225//32768 -1547//32768 287//2048 -749//4096 861//16384 0 0 0 0
153+
233//3456 -19//1152 0 0 0 -5//72 7//72 -17//216 0]
154+
x_star = [1//64 -1//64 7 / 192*sqrt(21) -7 / 192*sqrt(21) 0 0 0 0
155+
3 // 112+9 / 1960 * sqrt(7) -3 // 112+9 / 1960 * sqrt(7) 11 / 840 * sqrt(7)+3 / 112 * sqrt(7) * sqrt(3) 11 / 840 * sqrt(7)-3 / 112 * sqrt(7) * sqrt(3) 88 / 5145*sqrt(7) -18 / 343*sqrt(7) 0 0
156+
2707592511 // 1000000000000-1006699707 / 1000000000000 * sqrt(7) -51527976591 // 1000000000000-1006699707 / 1000000000000 * sqrt(7) -610366393//75000000000+7046897949 / 1000000000000*sqrt(7)+14508670449/1000000000000*sqrt(7)*sqrt(3) -610366393 // 75000000000 + 7046897949 / 1000000000000 * sqrt(7)-14508670449 / 1000000000000 * sqrt(7) * sqrt(3) -12456457 // 1171875000+1006699707 / 109375000000 * sqrt(7) 3020099121 / 437500000000 * sqrt(7)+47328957 // 625000000 -7046897949 / 250000000000*sqrt(7) 0]
157+
158+
τ_star = 0.4
159+
160+
TU = MIRKTableau(s, T.(c), T.(v), T.(b), T.(x))
161+
ITU = MIRKInterpTableau(s_star, T.(c_star), T.(v_star), T.(x_star), T(τ_star))
162+
return TU, ITU
163+
end

‎lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl‎

Lines changed: 48 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -93,6 +93,11 @@ end
9393
sol = solve(prob, mirk_solver(Val(order)); dt = 0.2)
9494
@test norm(diff(first.(sol.u)) .+ 0.2, Inf) + abs(sol.u[1][1] - 5) < affineTol
9595
end
96+
97+
@testset "MIRK$(order)I" for order in (6,)
98+
sol = solve(prob, MIRK6I(); dt = 0.2)
99+
@test norm(diff(first.(sol.u)) .+ 0.2, Inf) + abs(sol.u[1][1] - 5) < affineTol
100+
end
96101
end
97102
end
98103

@@ -109,6 +114,15 @@ end
109114
@test_call target_modules=(BoundaryValueDiffEqMIRK,) solve(
110115
prob, solver; dt = 0.2)
111116
end
117+
118+
@testset "MIRK$(order)I" for order in (6,)
119+
solver = MIRK6I(; nlsolve = NewtonRaphson(),
120+
jac_alg = BVPJacobianAlgorithm(AutoForwardDiff(; chunksize = 2)))
121+
@test_opt target_modules=(BoundaryValueDiffEqMIRK,) solve(
122+
prob, solver; dt = 0.2)
123+
@test_call target_modules=(BoundaryValueDiffEqMIRK,) solve(
124+
prob, solver; dt = 0.2)
125+
end
112126
end
113127
end
114128

@@ -122,6 +136,11 @@ end
122136
dts, prob, mirk_solver(Val(order)); abstol = 1e-8, reltol = 1e-8)
123137
@test sim.𝒪est[:final]≈order atol=testTol
124138
end
139+
140+
@testset "MIRK$(order)I" for (_, order) in enumerate((6,))
141+
sim = test_convergence(dts, prob, MIRK6I(); abstol = 1e-8, reltol = 1e-8)
142+
@test sim.𝒪est[:final]≈order atol=testTol
143+
end
125144
end
126145
end
127146

@@ -153,6 +172,7 @@ end
153172
@test_nowarn solve(bvp1, MIRK4(; jac_alg); dt = 0.05)
154173
@test_nowarn solve(bvp1, MIRK5(; jac_alg); dt = 0.05)
155174
@test_nowarn solve(bvp1, MIRK6(; jac_alg); dt = 0.05)
175+
@test_nowarn solve(bvp1, MIRK6I(; jac_alg); dt = 0.05)
156176
end
157177

158178
@testitem "Interpolation" begin
@@ -219,6 +239,34 @@ end
219239
@test sol(0.04, Val{0})≈sol_analytic atol=testTol
220240
@test sol(0.04, Val{1})≈dsol_analytic atol=testTol
221241
end
242+
243+
@testset "Interpolation for adaptive MIRK$(order)I" for order in (6,)
244+
sol = solve(prob_bvp_linear, MIRK6I(); dt = 0.001)
245+
sol_analytic = prob_bvp_linear_analytic(nothing, λ, 0.001)
246+
247+
@test sol(0.001)≈sol_analytic atol=testTol
248+
@test sol(0.001; idxs = [1, 2])≈sol_analytic atol=testTol
249+
@test sol(0.001; idxs = 1)≈sol_analytic[1] atol=testTol
250+
@test sol(0.001; idxs = 2)≈sol_analytic[2] atol=testTol
251+
end
252+
253+
@testset "Interpolation for non-adaptive MIRK$(order)I" for order in (6,)
254+
sol = solve(prob_bvp_linear, MIRK6I(); dt = 0.001, adaptive = false)
255+
256+
@test_nowarn sol(0.01)
257+
@test_nowarn sol(0.01; idxs = [1, 2])
258+
@test_nowarn sol(0.01; idxs = 1)
259+
@test_nowarn sol(0.01; idxs = 2)
260+
end
261+
262+
@testset "Interpolation for solution derivative" for order in (6,)
263+
sol = solve(prob_bvp_linear, MIRK6I(); dt = 0.001)
264+
sol_analytic = prob_bvp_linear_analytic(nothing, λ, 0.04)
265+
dsol_analytic = prob_bvp_linear_analytic_derivative(nothing, λ, 0.04)
266+
267+
@test sol(0.04, Val{0})≈sol_analytic atol=testTol
268+
@test sol(0.04, Val{1})≈dsol_analytic atol=testTol
269+
end
222270
end
223271

224272
@testitem "Swirling Flow III" begin

0 commit comments

Comments
 (0)