Skip to content

Commit 78053c1

Browse files
ChrisRackauckasChrisRackauckas
authored
Merge pull request SciML#2274 from ParamThakkar123/Verner
Added Verner Methods
2 parents 31b1a49 + 512c7df commit 78053c1

22 files changed

+1450
-1321
lines changed

‎lib/OrdinaryDiffEqDefault/Project.toml‎

Lines changed: 21 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,21 @@
1+
name = "OrdinaryDiffEqDefault"
2+
uuid = "50262376-6c5a-4cf5-baba-aaf4f84d72d7"
3+
authors = ["ParamThakkar123 <[email protected]>"]
4+
version = "1.0.0"
5+
6+
[deps]
7+
EnumX = "4e289a0a-7415-4d19-859d-a7e5c4648b56"
8+
LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
9+
OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
10+
11+
[extras]
12+
DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d"
13+
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
14+
SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
15+
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
16+
17+
[compat]
18+
julia = "1.10"
19+
20+
[targets]
21+
test = ["DiffEqDevTools", "Random", "SafeTestsets", "Test"]

‎lib/OrdinaryDiffEqDefault/src/OrdinaryDiffEqDefault.jl‎

Lines changed: 15 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,15 @@
1+
module OrdinaryDiffEqDefault
2+
3+
using OrdinaryDiffEq: Vern7, Vern8, Vern9, Vern6, Tsit5, Rosenbrock23, Rodas5P, FBDF,
4+
alg_stability_size, beta2_default, beta1_default, AutoSwitchCache, ODEIntegrator,
5+
CompositeAlgorithm, OrdinaryDiffEqAlgorithm, OrdinaryDiffEqMutableCache, AutoAlgSwitch
6+
import OrdinaryDiffEq: is_mass_matrix_alg, default_autoswitch
7+
import LinearSolve
8+
using LinearAlgebra: I, isdiag
9+
using EnumX
10+
11+
include("default_alg.jl")
12+
13+
export DefaultODEAlgorithm
14+
15+
end # module OrdinaryDiffEqDefault

‎lib/OrdinaryDiffEqDefault/src/default_alg.jl‎

Lines changed: 121 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,121 @@
1+
EnumX.@enumx DefaultSolverChoice begin
2+
Tsit5 = 1
3+
Vern7 = 2
4+
Rosenbrock23 = 3
5+
Rodas5P = 4
6+
FBDF = 5
7+
KrylovFBDF = 6
8+
end
9+
10+
const NUM_NONSTIFF = 2
11+
const NUM_STIFF = 4
12+
const LOW_TOL = 1e-6
13+
const MED_TOL = 1e-2
14+
const EXTREME_TOL = 1e-9
15+
const SMALLSIZE = 50
16+
const MEDIUMSIZE = 500
17+
const STABILITY_SIZES = (alg_stability_size(Tsit5()), alg_stability_size(Vern7()))
18+
const DEFAULTBETA2S = (
19+
beta2_default(Tsit5()), beta2_default(Vern7()), beta2_default(Rosenbrock23()),
20+
beta2_default(Rodas5P()), beta2_default(FBDF()), beta2_default(FBDF()))
21+
const DEFAULTBETA1S = (
22+
beta1_default(Tsit5(), DEFAULTBETA2S[1]), beta1_default(Vern7(), DEFAULTBETA2S[2]),
23+
beta1_default(Rosenbrock23(), DEFAULTBETA2S[3]), beta1_default(
24+
Rodas5P(), DEFAULTBETA2S[4]),
25+
beta1_default(FBDF(), DEFAULTBETA2S[5]), beta1_default(FBDF(), DEFAULTBETA2S[6]))
26+
27+
callbacks_exists(integrator) = !isempty(integrator.opts.callbacks)
28+
current_nonstiff(current) = ifelse(current <= NUM_NONSTIFF, current, current - NUM_STIFF)
29+
30+
function DefaultODEAlgorithm(; lazy = true, stiffalgfirst = false, kwargs...)
31+
nonstiff = (Tsit5(), Vern7(lazy = lazy))
32+
stiff = (Rosenbrock23(; kwargs...), Rodas5P(; kwargs...), FBDF(; kwargs...),
33+
FBDF(; linsolve = LinearSolve.KrylovJL_GMRES(), kwargs...))
34+
AutoAlgSwitch(nonstiff, stiff; stiffalgfirst)
35+
end
36+
37+
function is_stiff(integrator, alg, ntol, stol, is_stiffalg, current)
38+
eigen_est, dt = integrator.eigen_est, integrator.dt
39+
stiffness = abs(eigen_est * dt /
40+
STABILITY_SIZES[nonstiffchoice(integrator.opts.reltol)]) # `abs` here is just for safety
41+
tol = is_stiffalg ? stol : ntol
42+
os = oneunit(stiffness)
43+
bool = stiffness > os * tol
44+
45+
if !bool
46+
integrator.alg.choice_function.successive_switches += 1
47+
else
48+
integrator.alg.choice_function.successive_switches = 0
49+
end
50+
51+
integrator.do_error_check = (integrator.alg.choice_function.successive_switches >
52+
integrator.alg.choice_function.switch_max || !bool) ||
53+
is_stiffalg
54+
bool
55+
end
56+
57+
function nonstiffchoice(reltol)
58+
x = if reltol < LOW_TOL
59+
DefaultSolverChoice.Vern7
60+
else
61+
DefaultSolverChoice.Tsit5
62+
end
63+
Int(x)
64+
end
65+
66+
function stiffchoice(reltol, len, mass_matrix)
67+
x = if len > MEDIUMSIZE
68+
DefaultSolverChoice.KrylovFBDF
69+
elseif len > SMALLSIZE
70+
DefaultSolverChoice.FBDF
71+
else
72+
if reltol < LOW_TOL || !isdiag(mass_matrix)
73+
DefaultSolverChoice.Rodas5P
74+
else
75+
DefaultSolverChoice.Rosenbrock23
76+
end
77+
end
78+
Int(x)
79+
end
80+
81+
function default_autoswitch(AS::AutoSwitchCache, integrator)
82+
len = length(integrator.u)
83+
reltol = integrator.opts.reltol
84+
85+
# Choose the starting method
86+
if AS.current == 0
87+
choice = if AS.stiffalgfirst || integrator.f.mass_matrix != I
88+
stiffchoice(reltol, len, integrator.f.mass_matrix)
89+
else
90+
nonstiffchoice(reltol)
91+
end
92+
AS.current = choice
93+
return AS.current
94+
end
95+
96+
dt = integrator.dt
97+
# Successive stiffness test positives are counted by a positive integer,
98+
# and successive stiffness test negatives are counted by a negative integer
99+
AS.count = is_stiff(integrator, AS.nonstiffalg, AS.nonstifftol, AS.stifftol,
100+
AS.is_stiffalg, AS.current) ?
101+
AS.count < 0 ? 1 : AS.count + 1 :
102+
AS.count > 0 ? -1 : AS.count - 1
103+
if integrator.f.mass_matrix != I
104+
#don't change anything
105+
elseif (!AS.is_stiffalg && AS.count > AS.maxstiffstep)
106+
integrator.dt = dt * AS.dtfac
107+
AS.is_stiffalg = true
108+
AS.current = stiffchoice(reltol, len, integrator.f.mass_matrix)
109+
elseif (AS.is_stiffalg && AS.count < -AS.maxnonstiffstep)
110+
integrator.dt = dt / AS.dtfac
111+
AS.is_stiffalg = false
112+
AS.current = nonstiffchoice(reltol)
113+
end
114+
return AS.current
115+
end
116+
117+
# hack for the default alg
118+
function is_mass_matrix_alg(alg::CompositeAlgorithm{
119+
<:Any, <:Tuple{Tsit5, Vern7, Rosenbrock23, Rodas5P, FBDF, FBDF}})
120+
true
121+
end

‎lib/OrdinaryDiffEqVerner/Project.toml‎

Lines changed: 21 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,21 @@
1+
name = "OrdinaryDiffEqVerner"
2+
uuid = "79d7bb75-1356-48c1-b8c0-6832512096c2"
3+
authors = ["ParamThakkar123 <[email protected]>"]
4+
version = "1.0.0"
5+
6+
[deps]
7+
OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
8+
TruncatedStacktraces = "781d530d-4396-4725-bb49-402e4bee1e77"
9+
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
10+
11+
[extras]
12+
DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d"
13+
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
14+
SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
15+
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
16+
17+
[compat]
18+
julia = "1.10"
19+
20+
[targets]
21+
test = ["DiffEqDevTools", "Random", "SafeTestsets", "Test"]

‎lib/OrdinaryDiffEqVerner/src/OrdinaryDiffEqVerner.jl‎

Lines changed: 35 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,35 @@
1+
module OrdinaryDiffEqVerner
2+
3+
import OrdinaryDiffEq: alg_order, calculate_residuals!,
4+
initialize!, perform_step!, @unpack, unwrap_alg,
5+
calculate_residuals, alg_stability_size,
6+
OrdinaryDiffEqAlgorithm,
7+
CompositeAlgorithm, AbstractController, PIDController,
8+
OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache,
9+
OrdinaryDiffEqAdaptiveAlgorithm, CompiledFloats, uses_uprev,
10+
alg_cache, _vec, _reshape, @cache, isfsal, full_cache,
11+
constvalue, _unwrap_val, du_alias_or_new,
12+
explicit_rk_docstring, trivial_limiter!, _ode_interpolant,
13+
_ode_interpolant!, _ode_addsteps!, @fold,
14+
@OnDemandTableauExtract, AutoAlgSwitch,
15+
DerivativeOrderNotPossibleError
16+
using DiffEqBase, FastBroadcast, Polyester, MuladdMacro, RecursiveArrayTools
17+
using DiffEqBase: @def, @tight_loop_macros
18+
using Static: False
19+
using TruncatedStacktraces
20+
using LinearAlgebra: norm
21+
22+
include("algorithms.jl")
23+
include("alg_utils.jl")
24+
include("verner_tableaus.jl")
25+
include("verner_caches.jl")
26+
include("verner_addsteps.jl")
27+
include("interp_func.jl")
28+
include("interpolants.jl")
29+
include("controllers.jl")
30+
include("verner_rk_perform_step.jl")
31+
32+
export Vern6, Vern7, Vern8, Vern9
33+
export AutoVern6, AutoVern7, AutoVern8, AutoVern9
34+
35+
end

‎lib/OrdinaryDiffEqVerner/src/alg_utils.jl‎

Lines changed: 13 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,13 @@
1+
isfsal(alg::Vern7) = false
2+
isfsal(alg::Vern8) = false
3+
isfsal(alg::Vern9) = false
4+
5+
alg_order(alg::Vern6) = 6
6+
alg_order(alg::Vern7) = 7
7+
alg_order(alg::Vern8) = 8
8+
alg_order(alg::Vern9) = 9
9+
10+
alg_stability_size(alg::Vern6) = 4.8553
11+
alg_stability_size(alg::Vern7) = 4.6400
12+
alg_stability_size(alg::Vern8) = 5.8641
13+
alg_stability_size(alg::Vern9) = 4.4762

‎lib/OrdinaryDiffEqVerner/src/algorithms.jl‎

Lines changed: 119 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,119 @@
1+
@doc explicit_rk_docstring(
2+
"Verner's “Most Efficient” 6/5 Runge-Kutta method. (lazy 6th order interpolant).",
3+
"Vern6",
4+
references = "@article{verner2010numerically,
5+
title={Numerically optimal Runge--Kutta pairs with interpolants},
6+
author={Verner, James H},
7+
journal={Numerical Algorithms},
8+
volume={53},
9+
number={2-3},
10+
pages={383--396},
11+
year={2010},
12+
publisher={Springer}
13+
}",
14+
extra_keyword_description = """- `lazy`: determines if the lazy interpolant is used.
15+
""",
16+
extra_keyword_default = "lazy = true")
17+
Base.@kwdef struct Vern6{StageLimiter, StepLimiter, Thread} <:
18+
OrdinaryDiffEqAdaptiveAlgorithm
19+
stage_limiter!::StageLimiter = trivial_limiter!
20+
step_limiter!::StepLimiter = trivial_limiter!
21+
thread::Thread = False()
22+
lazy::Bool = true
23+
end
24+
TruncatedStacktraces.@truncate_stacktrace Vern6 3
25+
# for backwards compatibility
26+
function Vern6(stage_limiter!, step_limiter! = trivial_limiter!; lazy = true)
27+
Vern6(stage_limiter!, step_limiter!, False(), lazy)
28+
end
29+
30+
@doc explicit_rk_docstring(
31+
"Verner's “Most Efficient” 7/6 Runge-Kutta method. (lazy 7th order interpolant).",
32+
"Vern7",
33+
references = "@article{verner2010numerically,
34+
title={Numerically optimal Runge--Kutta pairs with interpolants},
35+
author={Verner, James H},
36+
journal={Numerical Algorithms},
37+
volume={53},
38+
number={2-3},
39+
pages={383--396},
40+
year={2010},
41+
publisher={Springer}
42+
}",
43+
extra_keyword_description = """- `lazy`: determines if the lazy interpolant is used.
44+
""",
45+
extra_keyword_default = "lazy = true")
46+
Base.@kwdef struct Vern7{StageLimiter, StepLimiter, Thread} <:
47+
OrdinaryDiffEqAdaptiveAlgorithm
48+
stage_limiter!::StageLimiter = trivial_limiter!
49+
step_limiter!::StepLimiter = trivial_limiter!
50+
thread::Thread = False()
51+
lazy::Bool = true
52+
end
53+
TruncatedStacktraces.@truncate_stacktrace Vern7 3
54+
# for backwards compatibility
55+
function Vern7(stage_limiter!, step_limiter! = trivial_limiter!; lazy = true)
56+
Vern7(stage_limiter!, step_limiter!, False(), lazy)
57+
end
58+
59+
@doc explicit_rk_docstring(
60+
"Verner's “Most Efficient” 8/7 Runge-Kutta method. (lazy 8th order interpolant).",
61+
"Vern8",
62+
references = "@article{verner2010numerically,
63+
title={Numerically optimal Runge--Kutta pairs with interpolants},
64+
author={Verner, James H},
65+
journal={Numerical Algorithms},
66+
volume={53},
67+
number={2-3},
68+
pages={383--396},
69+
year={2010},
70+
publisher={Springer}
71+
}",
72+
extra_keyword_description = """- `lazy`: determines if the lazy interpolant is used.
73+
""",
74+
extra_keyword_default = "lazy = true")
75+
Base.@kwdef struct Vern8{StageLimiter, StepLimiter, Thread} <:
76+
OrdinaryDiffEqAdaptiveAlgorithm
77+
stage_limiter!::StageLimiter = trivial_limiter!
78+
step_limiter!::StepLimiter = trivial_limiter!
79+
thread::Thread = False()
80+
lazy::Bool = true
81+
end
82+
TruncatedStacktraces.@truncate_stacktrace Vern8 3
83+
# for backwards compatibility
84+
function Vern8(stage_limiter!, step_limiter! = trivial_limiter!; lazy = true)
85+
Vern8(stage_limiter!, step_limiter!, False(), lazy)
86+
end
87+
88+
@doc explicit_rk_docstring(
89+
"Verner's “Most Efficient” 9/8 Runge-Kutta method. (lazy9th order interpolant).",
90+
"Vern9",
91+
references = "@article{verner2010numerically,
92+
title={Numerically optimal Runge--Kutta pairs with interpolants},
93+
author={Verner, James H},
94+
journal={Numerical Algorithms},
95+
volume={53},
96+
number={2-3},
97+
pages={383--396},
98+
year={2010},
99+
publisher={Springer}
100+
}",
101+
extra_keyword_description = """- `lazy`: determines if the lazy interpolant is used.
102+
""", extra_keyword_default = "lazy = true")
103+
Base.@kwdef struct Vern9{StageLimiter, StepLimiter, Thread} <:
104+
OrdinaryDiffEqAdaptiveAlgorithm
105+
stage_limiter!::StageLimiter = trivial_limiter!
106+
step_limiter!::StepLimiter = trivial_limiter!
107+
thread::Thread = False()
108+
lazy::Bool = true
109+
end
110+
TruncatedStacktraces.@truncate_stacktrace Vern9 3
111+
# for backwards compatibility
112+
function Vern9(stage_limiter!, step_limiter! = trivial_limiter!; lazy = true)
113+
Vern9(stage_limiter!, step_limiter!, False(), lazy)
114+
end
115+
116+
AutoVern6(alg; lazy = true, kwargs...) = AutoAlgSwitch(Vern6(lazy = lazy), alg; kwargs...)
117+
AutoVern7(alg; lazy = true, kwargs...) = AutoAlgSwitch(Vern7(lazy = lazy), alg; kwargs...)
118+
AutoVern8(alg; lazy = true, kwargs...) = AutoAlgSwitch(Vern8(lazy = lazy), alg; kwargs...)
119+
AutoVern9(alg; lazy = true, kwargs...) = AutoAlgSwitch(Vern9(lazy = lazy), alg; kwargs...)

‎lib/OrdinaryDiffEqVerner/src/controllers.jl‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,8 @@
1+
# checks whether the controller should accept a step based on the error estimate
2+
@inline function accept_step_controller(integrator, controller::AbstractController)
3+
return integrator.EEst <= 1
4+
end
5+
6+
@inline function accept_step_controller(integrator, controller::PIDController)
7+
return integrator.qold >= controller.accept_safety
8+
end

‎lib/OrdinaryDiffEqVerner/src/interp_func.jl‎

Lines changed: 31 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,31 @@
1+
function DiffEqBase.interp_summary(::Type{cacheType},
2+
dense::Bool) where {
3+
cacheType <:
4+
Union{Vern6Cache, Vern6ConstantCache
5+
}}
6+
dense ? "specialized 6th order lazy interpolation" : "1st order linear"
7+
end
8+
9+
function DiffEqBase.interp_summary(::Type{cacheType},
10+
dense::Bool) where {
11+
cacheType <:
12+
Union{Vern7Cache, Vern7ConstantCache
13+
}}
14+
dense ? "specialized 7th order lazy interpolation" : "1st order linear"
15+
end
16+
17+
function DiffEqBase.interp_summary(::Type{cacheType},
18+
dense::Bool) where {
19+
cacheType <:
20+
Union{Vern8Cache, Vern8ConstantCache
21+
}}
22+
dense ? "specialized 8th order lazy interpolation" : "1st order linear"
23+
end
24+
25+
function DiffEqBase.interp_summary(::Type{cacheType},
26+
dense::Bool) where {
27+
cacheType <:
28+
Union{Vern9Cache, Vern9ConstantCache
29+
}}
30+
dense ? "specialized 9th order lazy interpolation" : "1st order linear"
31+
end

0 commit comments

Comments
 (0)