Skip to content

Commit 154e1b5

Browse files
Add GLEE23/GLEE24/GLEE35 global-error-estimating GLM solvers
Implements the explicit general linear methods with built-in global error estimation from Constantinescu (2016), SIAM J. Numer. Anal. 54(6) (arXiv:1503.05166; PETSc's TSGLEE23/TSGLEE24/TSGLEE35), addressing SciML/GlobalDiffEq.jl#6 and SciML/GlobalDiffEq.jl#22, as proper OrdinaryDiffEqCore algorithms. The methods propagate the partitioned state (y, ε) where ε is an asymptotically correct global error estimate; solve/init on plain ODEProblems transparently extend to the ArrayPartition form and global_error_estimate(sol) extracts the estimates. The per-step ε increment is an asymptotically correct local error estimate driving standard step-size adaptivity. The paper's method A4 is defective as printed (violates b·c = 1/2, converges at order 1, absent from PETSc) and is deliberately not implemented. Verified locally on the unstable Prince42 problem: solution orders 2/2/3, estimate accuracy orders 3/3/4, est/true-error ratio 1.00. Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
1 parent ead709c commit 154e1b5

11 files changed

Lines changed: 671 additions & 13 deletions

File tree

docs/src/globalerrorcontrol/GlobalDiffEq.md

Lines changed: 39 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -14,15 +14,45 @@ To use these methods:
1414
using GlobalDiffEq
1515
```
1616

17-
[`GlobalRichardson`](@ref) wraps any fixed-step method in global Richardson
18-
extrapolation over whole solves, interpreting `abstol` and `reltol` as global
19-
tolerances. It is the most robust and most expensive option.
20-
21-
To *control* the endpoint global error to a tolerance `gtol`, wrap any
22-
adaptive solver in [`GlobalAdjoint`](@ref) (adjoint-based, for endpoint
23-
functionals; requires SciMLSensitivity and QuadGK to be loaded). It solves the
24-
problem, estimates the endpoint global error, and tightens the local
25-
tolerances until the requested global tolerance is met.
17+
## Choosing a method
18+
19+
- For a solution accompanied by a running, asymptotically correct estimate of
20+
its global error at every time point, use the global-error-estimating
21+
solvers [`GLEE24`](@ref) or [`GLEE35`](@ref) (Constantinescu 2016). These
22+
cost only a few extra stages per step over a plain method of the same order
23+
and require nothing beyond the right-hand side `f`.
24+
- To *control* the endpoint global error to a tolerance `gtol`, wrap any
25+
adaptive solver in [`GlobalAdjoint`](@ref) (adjoint-based, for endpoint
26+
functionals; requires SciMLSensitivity and QuadGK to be loaded). It solves
27+
the problem, estimates the endpoint global error, and tightens the local
28+
tolerances until the requested global tolerance is met.
29+
- [`GlobalRichardson`](@ref) wraps any fixed-step method in global Richardson
30+
extrapolation over whole solves, interpreting `abstol` and `reltol` as
31+
global tolerances. It is the most robust and most expensive option.
32+
33+
For example, solving while tracking the global error along the trajectory:
34+
35+
```julia
36+
using GlobalDiffEq
37+
38+
function lorenz!(du, u, p, t)
39+
du[1] = 10.0(u[2] - u[1])
40+
du[2] = u[1] * (28.0 - u[3]) - u[2]
41+
du[3] = u[1] * u[2] - (8 / 3) * u[3]
42+
end
43+
prob = ODEProblem(lorenz!, [1.0; 0.0; 0.0], (0.0, 10.0))
44+
sol = solve(prob, GLEE35(); abstol = 1.0e-8, reltol = 1.0e-8)
45+
errs = global_error_estimate(sol) # global error estimate at every sol.t
46+
```
47+
48+
## Global-error-estimating solvers
49+
50+
```@docs
51+
GLEE23
52+
GLEE24
53+
GLEE35
54+
global_error_estimate
55+
```
2656

2757
## Global error controlling wrappers
2858

lib/GlobalDiffEq/Project.toml

Lines changed: 12 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,14 +1,18 @@
11
name = "GlobalDiffEq"
22
uuid = "1d72d19b-84cc-4cb7-a099-7cbdb9ccc67c"
33
authors = ["Chris Rackauckas <accounts@chrisrackauckas.com> and contributors"]
4-
version = "1.3.0"
4+
version = "1.4.0"
55

66
[deps]
77
DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
8+
FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898"
89
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
10+
MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
11+
OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8"
912
OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
1013
PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
1114
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
15+
RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
1216
Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
1317
Richardson = "708f8203-808e-40c0-ba2d-98a6953ed40d"
1418
SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
@@ -23,16 +27,21 @@ GlobalDiffEqSciMLSensitivityExt = ["SciMLSensitivity", "QuadGK"]
2327

2428
[sources]
2529
DiffEqBase = {path = "../DiffEqBase"}
30+
OrdinaryDiffEqCore = {path = "../OrdinaryDiffEqCore"}
2631
OrdinaryDiffEqTsit5 = {path = "../OrdinaryDiffEqTsit5"}
2732

2833
[compat]
2934
DiffEqBase = "7"
35+
FastBroadcast = "1.3"
3036
LinearAlgebra = "1"
37+
MuladdMacro = "0.2.1"
38+
OrdinaryDiffEqCore = "4"
3139
OrdinaryDiffEqSSPRK = "2"
3240
OrdinaryDiffEqTsit5 = "2"
3341
PrecompileTools = "1.1"
3442
QuadGK = "2.11"
3543
Random = "1"
44+
RecursiveArrayTools = "4.2.0"
3645
Reexport = "0.2, 1.0"
3746
Richardson = "1.2"
3847
SafeTestsets = "0.0.1, 0.1"
@@ -49,10 +58,11 @@ OrdinaryDiffEqSSPRK = "669c94d9-1f4b-4b64-b377-1aa079aa2388"
4958
OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
5059
QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
5160
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
61+
RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
5262
SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
5363
SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1"
5464
SciMLTesting = "09d9d899-5365-40a9-917a-5f67fddea283"
5565
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
5666

5767
[targets]
58-
test = ["LinearAlgebra", "OrdinaryDiffEqSSPRK", "OrdinaryDiffEqTsit5", "QuadGK", "Random", "SafeTestsets", "SciMLSensitivity", "SciMLTesting", "Test"]
68+
test = ["LinearAlgebra", "OrdinaryDiffEqSSPRK", "OrdinaryDiffEqTsit5", "QuadGK", "Random", "RecursiveArrayTools", "SafeTestsets", "SciMLSensitivity", "SciMLTesting", "Test"]

lib/GlobalDiffEq/README.md

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -12,6 +12,9 @@ provides:
1212

1313
- `GlobalRichardson`: global Richardson extrapolation of whole solves of any
1414
fixed-step method, interpreting `abstol`/`reltol` as global tolerances.
15+
- `GLEE23`, `GLEE24`, `GLEE35`: explicit general linear methods with built-in
16+
global error estimation (Constantinescu 2016), which propagate the solution
17+
together with an asymptotically correct estimate of its global error.
1518
- `GlobalAdjoint`: adjoint-based a posteriori endpoint error estimation and
1619
control (Cao and Petzold 2004), available as a package extension when
1720
SciMLSensitivity and QuadGK are loaded.

lib/GlobalDiffEq/src/GlobalDiffEq.jl

Lines changed: 12 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -3,16 +3,26 @@ module GlobalDiffEq
33
using Reexport: @reexport
44
@reexport using DiffEqBase
55

6-
import LinearAlgebra, OrdinaryDiffEqTsit5, Random, Richardson, SciMLBase,
7-
SciMLStructures
6+
import LinearAlgebra, OrdinaryDiffEqCore, OrdinaryDiffEqTsit5, Random,
7+
RecursiveArrayTools, Richardson, SciMLBase, SciMLStructures
8+
import DiffEqBase: initialize!, calculate_residuals, calculate_residuals!
9+
import OrdinaryDiffEqCore: perform_step!, @cache
10+
using FastBroadcast: FastBroadcast, @..
11+
using MuladdMacro: MuladdMacro, @muladd
812
using PrecompileTools: @setup_workload, @compile_workload
913

1014
abstract type GlobalDiffEqAlgorithm <: SciMLBase.AbstractODEAlgorithm end
1115

1216
include("richardson.jl")
1317
include("adjoint.jl")
18+
include("glee/tableaus.jl")
19+
include("glee/algorithms.jl")
20+
include("glee/solve.jl")
21+
include("glee/caches.jl")
22+
include("glee/perform_step.jl")
1423

1524
export GlobalAdjoint, GlobalRichardson, adjoint_error_estimate
25+
export GLEE23, GLEE24, GLEE35, global_error_estimate
1626

1727
@setup_workload begin
1828
# Simple test ODE: exponential decay du/dt = -u
Lines changed: 70 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,70 @@
1+
# Explicit general linear methods with built-in global error estimation
2+
# (GLEE methods) from Constantinescu (2016), doi:10.1137/15M1014633.
3+
abstract type AbstractGLEEAlgorithm <:
4+
OrdinaryDiffEqCore.OrdinaryDiffEqAdaptiveAlgorithm end
5+
6+
const _GLEE_DOCS_SHARED = """
7+
GLEE methods propagate the solution `y` together with an asymptotically
8+
correct estimate `ε` of its global (accumulated) error, at the cost of a few
9+
extra stages per step. Solving any `ODEProblem` with a GLEE method produces a
10+
solution whose states are `ArrayPartition`s: `sol.u[i].x[1]` is the solution
11+
and `sol.u[i].x[2]` is the global error estimate at `sol.t[i]` (see
12+
[`global_error_estimate`](@ref)). The per-step increment of `ε` is an
13+
asymptotically correct local error estimate, which drives standard step-size
14+
adaptivity, so local tolerances behave exactly as for ordinary adaptive
15+
Runge-Kutta methods while the global error is estimated for free.
16+
17+
Only explicit, mass-matrix-free ODEs are supported. The reference for the
18+
methods and their theory is:
19+
20+
Emil M. Constantinescu, *Estimating Global Errors in Time Stepping*, SIAM
21+
Journal on Numerical Analysis 54(6), 2016. [arXiv:1503.05166](https://arxiv.org/abs/1503.05166)
22+
"""
23+
24+
"""
25+
GLEE23()
26+
27+
3-stage, second-order explicit general linear method with global error
28+
estimation (Constantinescu 2016, eq. (4.6); PETSc's `TSGLEE23`). The cheapest
29+
GLEE method: only the first decoupling condition holds, so prefer
30+
[`GLEE24`](@ref) for long-time integration or mildly stiff problems.
31+
32+
$(_GLEE_DOCS_SHARED)
33+
"""
34+
struct GLEE23 <: AbstractGLEEAlgorithm end
35+
36+
"""
37+
GLEE24()
38+
39+
4-stage, second-order explicit general linear method with global error
40+
estimation (Constantinescu 2016, eq. (A.3); PETSc's `TSGLEE24`). Satisfies
41+
both decoupling conditions (`B·U` and `B·A·U` diagonal), which keeps the error
42+
estimate faithful in long-time integration; this is the recommended
43+
second-order GLEE method.
44+
45+
$(_GLEE_DOCS_SHARED)
46+
"""
47+
struct GLEE24 <: AbstractGLEEAlgorithm end
48+
49+
"""
50+
GLEE35()
51+
52+
5-stage, third-order explicit general linear method with global error
53+
estimation (Constantinescu 2016, eq. (4.9); PETSc's `TSGLEE35`). Satisfies
54+
both decoupling conditions and has a large negative-real-axis stability
55+
region; this is the recommended third-order GLEE method.
56+
57+
$(_GLEE_DOCS_SHARED)
58+
"""
59+
struct GLEE35 <: AbstractGLEEAlgorithm end
60+
61+
SciMLBase.alg_order(::GLEE23) = 2
62+
SciMLBase.alg_order(::GLEE24) = 2
63+
SciMLBase.alg_order(::GLEE35) = 3
64+
65+
OrdinaryDiffEqCore.alg_adaptive_order(alg::AbstractGLEEAlgorithm) =
66+
SciMLBase.alg_order(alg)
67+
68+
_glee_tableau_for(::GLEE23, ::Type{T}, ::Type{T2}) where {T, T2} = GLEE23Tableau(T, T2)
69+
_glee_tableau_for(::GLEE24, ::Type{T}, ::Type{T2}) where {T, T2} = GLEE24Tableau(T, T2)
70+
_glee_tableau_for(::GLEE35, ::Type{T}, ::Type{T2}) where {T, T2} = GLEE35Tableau(T, T2)
Lines changed: 52 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,52 @@
1+
@cache struct GLEECache{uType, yType, rateType, yRateType, uNoUnitsType, TabType} <:
2+
OrdinaryDiffEqCore.OrdinaryDiffEqMutableCache
3+
u::uType
4+
uprev::uType
5+
tmp::uType
6+
fsalfirst::rateType
7+
fsallast::rateType
8+
ks::Vector{yRateType}
9+
ytmp::yType
10+
εloc::yType
11+
atmp::uNoUnitsType
12+
tab::TabType
13+
end
14+
15+
OrdinaryDiffEqCore.get_fsalfirstlast(cache::GLEECache, u) =
16+
(cache.fsalfirst, cache.fsallast)
17+
18+
struct GLEEConstantCache{TabType} <: OrdinaryDiffEqCore.OrdinaryDiffEqConstantCache
19+
tab::TabType
20+
end
21+
22+
function OrdinaryDiffEqCore.alg_cache(
23+
alg::AbstractGLEEAlgorithm, u, rate_prototype, ::Type{uEltypeNoUnits},
24+
::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
25+
dt, reltol, p, calck,
26+
::Val{true}, verbose
27+
) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
28+
tab = _glee_tableau_for(
29+
alg, OrdinaryDiffEqCore.constvalue(uBottomEltypeNoUnits),
30+
OrdinaryDiffEqCore.constvalue(tTypeNoUnits)
31+
)
32+
y = u.x[1]
33+
yrate = rate_prototype.x[1]
34+
ks = [zero(yrate) for _ in 1:nstages(tab)]
35+
return GLEECache(
36+
u, uprev, zero(u), zero(rate_prototype), zero(rate_prototype), ks,
37+
zero(y), zero(y), fill!(similar(y, uEltypeNoUnits), 0), tab
38+
)
39+
end
40+
41+
function OrdinaryDiffEqCore.alg_cache(
42+
alg::AbstractGLEEAlgorithm, u, rate_prototype, ::Type{uEltypeNoUnits},
43+
::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
44+
dt, reltol, p, calck,
45+
::Val{false}, verbose
46+
) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
47+
tab = _glee_tableau_for(
48+
alg, OrdinaryDiffEqCore.constvalue(uBottomEltypeNoUnits),
49+
OrdinaryDiffEqCore.constvalue(tTypeNoUnits)
50+
)
51+
return GLEEConstantCache(tab)
52+
end

0 commit comments

Comments
 (0)