Skip to content

Commit be2b2ad

Browse files
Merge pull request #287 from ErikQQY/qqy/global_error
Add global error control adaptivity
2 parents b99ea13 + 361b8bd commit be2b2ad

File tree

22 files changed

+786
-165
lines changed

22 files changed

+786
-165
lines changed

docs/pages.jl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -4,7 +4,7 @@ pages = ["index.md",
44
"Getting Started with BVP solving in Julia" => "tutorials/getting_started.md",
55
"Tutorials" => Any["tutorials/continuation.md", "tutorials/solve_nlls_bvp.md"],
66
"Basics" => Any["basics/bvp_problem.md", "basics/bvp_functions.md",
7-
"basics/solve.md", "basics/autodiff.md"],
7+
"basics/solve.md", "basics/autodiff.md", "basics/error_control.md"],
88
"Solver Summaries and Recommendations" => Any[
99
"solvers/mirk.md", "solvers/firk.md", "solvers/shooting.md", "solvers/mirkn.md",
1010
"solvers/ascher.md", "solvers/simple_solvers.md", "solvers/wrappers.md"],

docs/src/assets/Project.toml

Lines changed: 40 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,40 @@
1+
[deps]
2+
ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
3+
BoundaryValueDiffEq = "764a87c0-6b3e-53db-9096-fe964310641d"
4+
BoundaryValueDiffEqAscher = "7227322d-7511-4e07-9247-ad6ff830280e"
5+
BoundaryValueDiffEqCore = "56b672f2-a5fe-4263-ab2d-da677488eb3a"
6+
BoundaryValueDiffEqFIRK = "85d9eb09-370e-4000-bb32-543851f73618"
7+
BoundaryValueDiffEqMIRK = "1a22d4ce-7765-49ea-b6f2-13c8438986a6"
8+
BoundaryValueDiffEqMIRKN = "9255f1d6-53bf-473e-b6bd-23f1ff009da4"
9+
BoundaryValueDiffEqShooting = "ed55bfe0-3725-4db6-871e-a1dc9f42a757"
10+
DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
11+
Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
12+
DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"
13+
DocumenterInterLinks = "d12716ef-a0f6-4df4-a9f1-a5a34e75c656"
14+
LineSearch = "87fe0de2-c867-4266-b59a-2f0a94fc965b"
15+
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
16+
LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
17+
OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
18+
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
19+
SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
20+
SimpleBoundaryValueDiffEq = "be0294bd-f90f-4760-ac4e-3421ce2b2da0"
21+
22+
[compat]
23+
ADTypes = "1.9"
24+
BoundaryValueDiffEq = "5.13.0"
25+
BoundaryValueDiffEqAscher = "1.2.0"
26+
BoundaryValueDiffEqCore = "1.3.0"
27+
BoundaryValueDiffEqFIRK = "1.3.0"
28+
BoundaryValueDiffEqMIRK = "1.3.0"
29+
BoundaryValueDiffEqShooting = "1.3.0"
30+
DiffEqBase = "6.158.3"
31+
Documenter = "1"
32+
DocumenterCitations = "1"
33+
DocumenterInterLinks = "1.0.0"
34+
LineSearch = "0.1.4"
35+
LinearAlgebra = "1.10"
36+
LinearSolve = "2.36.2, 3"
37+
OrdinaryDiffEq = "6.90.1"
38+
Plots = "1"
39+
SciMLBase = "2.60.0"
40+
SimpleBoundaryValueDiffEq = "1.1.0"

docs/src/basics/error_control.md

Lines changed: 35 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,35 @@
1+
# Error Control Adaptivity
2+
3+
Adaptivity helps ensure the quality of the our numerical solution, and when our solution exhibits significant estimating errors, adaptivity automatically refine the mesh based on the error distribution, and providing a final satisfying solution.
4+
5+
When comes to solving ill-conditioned BVP, for example the singular pertubation problem where the small parameters become extremally small leading to the layers phonemona, the error control adaptivity becomes even more critical, because the minor pertubations can lead to large deviation in the solution. In such cases, adaptivity autimatically figure out where to use refined mesh and where to use coarse mesh to achieve the balance of computational efficiency and accuracy.
6+
7+
BoundaryValuDiffEq.jl support error control adaptivity for collocation methods, and the adaptivity is default as defect control adaptivity when using adaptive collocation solvers:
8+
9+
```julia
10+
sol = solve(prob, MIRK4(), dt = 0.01, adaptive = true)
11+
```
12+
13+
Actually, BoundaryValueDiffEq.jl supports both defect and global error control adaptivity(while the defect control is the default controller) [boisvert2013runge](@Citet), to specify different error control metods, we simply need to specify the `controller` keyword in `solve`:
14+
15+
```julia
16+
sol = solve(prob, MIRK4(), dt = 0.01, controller = GlobalErrorControl()) # Use global error control
17+
sol = solve(prob, MIRK4(), dt = 0.01, controller = SequentialErrorControl()) # Use Sequential error control
18+
sol = solve(prob, MIRK4(), dt = 0.01, controller = HybridErrorControl()) # Use Hybrid error control
19+
```
20+
21+
## Error control methods
22+
23+
```@docs
24+
DefectControl
25+
GlobalErrorControl
26+
SequentialErrorControl
27+
HybridErrorControl
28+
```
29+
30+
While we can achieve global error control in different ways, we can use different methods to estimate the global error:
31+
32+
```@docs
33+
HOErrorControl
34+
REErrorControl
35+
```

docs/src/refs.bib

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -27,3 +27,14 @@ @book{ascher1995numerical
2727
URL = {https://epubs.siam.org/doi/abs/10.1137/1.9781611971231},
2828
eprint = {https://epubs.siam.org/doi/pdf/10.1137/1.9781611971231}
2929
}
30+
31+
@article{boisvert2013runge,
32+
title={A Runge-Kutta BVODE solver with global error and defect control},
33+
author={Boisvert, Jason J and Muir, Paul H and Spiteri, Raymond J},
34+
journal={ACM Transactions on Mathematical Software (TOMS)},
35+
volume={39},
36+
number={2},
37+
pages={1--22},
38+
year={2013},
39+
publisher={ACM New York, NY, USA}
40+
}

docs/src/solvers/firk.md

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -20,7 +20,7 @@ solve(prob::TwoPointBVProblem, alg, dt; kwargs...)
2020

2121
### Radau IIA methods
2222

23-
- `RadauIIa1`: 1 stage Radau IIA method, with defect control adaptivity
23+
- `RadauIIa1`: 1 stage Radau IIA method, without defect control adaptivity
2424
- `RadauIIa2`: 2 stage Radau IIA method, with defect control adaptivity.
2525
- `RadauIIa3`: 3 stage Radau IIA method, with defect control adaptivity.
2626
- `RadauIIa5`: 5 stage Radau IIA method, with defect control adaptivity.
@@ -35,14 +35,14 @@ solve(prob::TwoPointBVProblem, alg, dt; kwargs...)
3535

3636
### Lobatto IIIB methods
3737

38-
- `LobattoIIIb2`: 2 stage Lobatto IIIb method, with defect control adaptivity.
38+
- `LobattoIIIb2`: 2 stage Lobatto IIIb method, without defect control adaptivity.
3939
- `LobattoIIIb3`: 3 stage Lobatto IIIb method, with defect control adaptivity.
4040
- `LobattoIIIb4`: 4 stage Lobatto IIIb method, with defect control adaptivity.
4141
- `LobattoIIIb5`: 5 stage Lobatto IIIb method, with defect control adaptivity.
4242

4343
### Lobatto IIIC methods
4444

45-
- `LobattoIIIc2`: 2 stage Lobatto IIIc method, with defect control adaptivity.
45+
- `LobattoIIIc2`: 2 stage Lobatto IIIc method, without defect control adaptivity.
4646
- `LobattoIIIc3`: 3 stage Lobatto IIIc method, with defect control adaptivity.
4747
- `LobattoIIIc4`: 4 stage Lobatto IIIc method, with defect control adaptivity.
4848
- `LobattoIIIc5`: 5 stage Lobatto IIIc method, with defect control adaptivity.

lib/BoundaryValueDiffEqCore/src/BoundaryValueDiffEqCore.jl

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -29,6 +29,7 @@ include("algorithms.jl")
2929
include("alg_utils.jl")
3030
include("default_nlsolve.jl")
3131
include("misc_utils.jl")
32+
include("calc_errors.jl")
3233

3334
function SciMLBase.__solve(
3435
prob::AbstractBVProblem, alg::BoundaryValueDiffEqAlgorithm, args...; kwargs...)
@@ -37,5 +38,7 @@ function SciMLBase.__solve(
3738
end
3839

3940
export BVPJacobianAlgorithm
41+
export DefectControl, GlobalErrorControl, SequentialErrorControl, HybridErrorControl
42+
export HOErrorControl, REErrorControl
4043

4144
end
Lines changed: 111 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,111 @@
1+
"""
2+
AbstractErrorControl
3+
4+
Abstract type for different error control methods.
5+
"""
6+
abstract type AbstractErrorControl end
7+
8+
"""
9+
GlobalErrorControlMethod
10+
11+
Abstract type for different global error control methods, and according to the different global error estimation methods, there are
12+
13+
- `HOErrorControl`: Higher order global error estimation method
14+
- `REErrorControl`: Richardson extrapolation global error estimation method
15+
"""
16+
abstract type GlobalErrorControlMethod end
17+
18+
"""
19+
DefectControl(; defect_threshold = 0.1)
20+
21+
Defect estimation method with defect defined as
22+
23+
```math
24+
defect = \\max\\frac{S'(x) - f(x,S(x))}{1 + |f(x,S(x))|}
25+
```
26+
27+
Defect controller, with the maximum `defect_threshold` as 0.1, when the estimating defect is greater than the `defect_threshold`, the mesh will be refined.
28+
"""
29+
struct DefectControl{T} <: AbstractErrorControl
30+
defect_threshold::T
31+
32+
function DefectControl(; defect_threshold = 0.1)
33+
return new{typeof(defect_threshold)}(defect_threshold)
34+
end
35+
end
36+
37+
"""
38+
GlobalErrorControl(; method = HOErorControl())
39+
40+
Global error controller, use high order global error estimation method `HOErrorControl` as default.
41+
"""
42+
struct GlobalErrorControl <: AbstractErrorControl
43+
method::GlobalErrorControlMethod
44+
45+
function GlobalErrorControl(; method = HOErrorControl())
46+
return new(method)
47+
end
48+
end
49+
50+
"""
51+
SequentialErrorControl(; defect = DefectControl(), global_error = GlobalErrorControl())
52+
53+
First use the defect controller, if the defect is satisfying, then use global error controller.
54+
"""
55+
struct SequentialErrorControl <: AbstractErrorControl
56+
defect::DefectControl
57+
global_error::GlobalErrorControl
58+
59+
function SequentialErrorControl(;
60+
defect = DefectControl(), global_error = GlobalErrorControl())
61+
return new(defect, global_error)
62+
end
63+
end
64+
65+
"""
66+
HybridErrorControl(; DE = 1.0, GE = 1.0, defect = DefectControl(), global_error = GlobalErrorControl())
67+
68+
Control both of the defect and global error, where the error norm is the linear combination of the defect and global error.
69+
"""
70+
struct HybridErrorControl{T1, T2} <: AbstractErrorControl
71+
DE::T1
72+
GE::T2
73+
defect::DefectControl
74+
global_error::GlobalErrorControl
75+
76+
function HybridErrorControl(; DE = 1.0, GE = 1.0, defect = DefectControl(),
77+
global_error = GlobalErrorControl())
78+
return new{typeof(DE), typeof(GE)}(DE, GE, defect, global_error)
79+
end
80+
end
81+
82+
"""
83+
HOErrorControl()
84+
85+
Higher order global error estimation method
86+
87+
Uses a solution from order+2 method on the original mesh and calculate the error with
88+
89+
```math
90+
error = \\max\\frac{u_p - u_{p+2}}{1 + |u_p|}
91+
```
92+
"""
93+
struct HOErrorControl <: GlobalErrorControlMethod end
94+
95+
"""
96+
REErrorControl()
97+
98+
Richardson extrapolation global error estimation method
99+
100+
Use Richardson extrapolation to calculate the error on the doubled mesh with
101+
102+
```math
103+
error = \\frac{2^p}{2^p-1} * \\max\\frac{u_h - u_{h/2}}{1 + |u_h|}
104+
```
105+
"""
106+
struct REErrorControl <: GlobalErrorControlMethod end
107+
108+
# Some utils for error control adaptivity
109+
# If error control use both defect and global error or not
110+
@inline __use_both_error_control(::HybridErrorControl) = true
111+
@inline __use_both_error_control(_) = false

lib/BoundaryValueDiffEqCore/src/utils.jl

Lines changed: 23 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -144,33 +144,42 @@ function eval_bc_residual!(
144144
bcb!(resid[2], dub, ub, p)
145145
end
146146

147-
__append_similar!(::Nothing, n, _) = nothing
147+
"""
148+
__resize!(x, n, M)
149+
150+
Resizes the input `x` to length `n` and returns the resized array. If `n` is less than the
151+
length of `x`, it truncates the array. If `n` is greater than the length of `x`, it appends
152+
zeros to the array.
148153
149-
# NOTE: We use `last` since the `first` might not conform to the same structure. For eg,
150-
# in the case of residuals
151-
function __append_similar!(x::AbstractVector{<:AbstractArray}, n, _)
154+
!!! note
155+
156+
We use `last` since the `first` might not conform to the same structure. For example, in the case of residuals
157+
"""
158+
function __resize!(x::AbstractVector{<:AbstractArray}, n, M)
152159
N = n - length(x)
153160
N == 0 && return x
154-
N < 0 && throw(ArgumentError("Cannot append a negative number of elements"))
155-
append!(x, [zero(last(x)) for _ in 1:N])
161+
N > 0 ? append!(x, [zero(last(x)) for _ in 1:N]) : resize!(x, n)
156162
return x
157163
end
158164

159-
function __append_similar!(x::AbstractVector{<:MaybeDiffCache}, n, M)
165+
__resize!(::Nothing, n, _) = nothing
166+
167+
function __resize!(x::AbstractVector{<:MaybeDiffCache}, n, M)
160168
N = n - length(x)
161169
N == 0 && return x
162-
N < 0 && throw(ArgumentError("Cannot append a negative number of elements"))
163-
chunksize = pickchunksize(M * (N + length(x)))
164-
append!(x, [__maybe_allocate_diffcache(last(x), chunksize) for _ in 1:N])
170+
if N > 0
171+
chunksize = pickchunksize(M * (N + length(x)))
172+
append!(x, [__maybe_allocate_diffcache(last(x), chunksize) for _ in 1:N])
173+
else
174+
resize!(x, n)
175+
end
165176
return x
166177
end
167178

168-
__append_similar!(::Nothing, n, _, _) = nothing
169-
function __append_similar!(x::AbstractVectorOfArray, n, _)
179+
function __resize!(x::AbstractVectorOfArray, n, M)
170180
N = n - length(x)
171181
N == 0 && return x
172-
N < 0 && throw(ArgumentError("Cannot append a negative number of elements"))
173-
append!(x, VectorOfArray([similar(last(x)) for _ in 1:N]))
182+
N > 0 ? append!(x, VectorOfArray([similar(last(x)) for _ in 1:N])) : resize!(x, n)
174183
return x
175184
end
176185

lib/BoundaryValueDiffEqFIRK/src/BoundaryValueDiffEqFIRK.jl

Lines changed: 7 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -9,14 +9,13 @@ using BoundaryValueDiffEqCore: BoundaryValueDiffEqAlgorithm, BVPJacobianAlgorith
99
__FastShortcutBVPCompatibleNonlinearPolyalg,
1010
__FastShortcutBVPCompatibleNLLSPolyalg, EvalSol,
1111
concrete_jacobian_algorithm, eval_bc_residual,
12-
eval_bc_residual!, get_tmp, __maybe_matmul!,
13-
__append_similar!, __extract_problem_details,
14-
__initial_guess, __maybe_allocate_diffcache,
15-
__restructure_sol, __get_bcresid_prototype, __similar, __vec,
16-
__vec_f, __vec_f!, __vec_bc, __vec_bc!,
17-
recursive_flatten_twopoint!, __internal_nlsolve_problem,
18-
MaybeDiffCache, __extract_mesh, __extract_u0,
19-
__has_initial_guess, __initial_guess_length,
12+
eval_bc_residual!, get_tmp, __maybe_matmul!, __resize!,
13+
__extract_problem_details, __initial_guess,
14+
__maybe_allocate_diffcache, __restructure_sol,
15+
__get_bcresid_prototype, __similar, __vec, __vec_f, __vec_f!,
16+
__vec_bc, __vec_bc!, recursive_flatten_twopoint!,
17+
__internal_nlsolve_problem, MaybeDiffCache, __extract_mesh,
18+
__extract_u0, __has_initial_guess, __initial_guess_length,
2019
__initial_guess_on_mesh, __flatten_initial_guess,
2120
__build_solution, __Fix3, _sparse_like, get_dense_ad
2221

lib/BoundaryValueDiffEqFIRK/src/adaptivity.jl

Lines changed: 7 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -182,9 +182,10 @@ Generate new mesh based on the defect.
182182
n_ = T(0.1) * n
183183
n_predict = ifelse(abs((n_predict - n)) < n_, round(Int, n + n_), n_predict)
184184

185-
if r₁ ρ * r
185+
if r₁ ρ * r
186186
Nsub_star = 2 * (N - 1)
187-
if Nsub_star > cache.alg.max_num_subintervals # Need to determine the too large threshold
187+
# Need to determine the too large threshold
188+
if Nsub_star > cache.alg.max_num_subintervals
188189
info = ReturnCode.Failure
189190
meshₒ = mesh
190191
mesh_dt₀ = mesh_dt
@@ -218,14 +219,14 @@ Generate a new mesh based on the `ŝ`.
218219
"""
219220
function redistribute!(cache::Union{FIRKCacheExpand{iip, T}, FIRKCacheNested{iip, T}},
220221
Nsub_star, ŝ, mesh, mesh_dt) where {iip, T}
221-
N = length(mesh)
222+
N = length(mesh) - 1
222223
ζ = sum(ŝ .* mesh_dt) / Nsub_star
223224
k, i = 1, 0
224-
append!(cache.mesh, Nsub_star + 1 - N)
225+
resize!(cache.mesh, Nsub_star + 1)
225226
cache.mesh[1] = mesh[1]
226227
t = mesh[1]
227228
integral = T(0)
228-
while k N - 1
229+
while k N
229230
next_piece = ŝ[k] * (mesh[k + 1] - t)
230231
_int_next = integral + next_piece
231232
if _int_next > ζ
@@ -240,7 +241,7 @@ function redistribute!(cache::Union{FIRKCacheExpand{iip, T}, FIRKCacheNested{iip
240241
end
241242
end
242243
cache.mesh[end] = mesh[end]
243-
append!(cache.mesh_dt, Nsub_star - N)
244+
resize!(cache.mesh_dt, Nsub_star)
244245
diff!(cache.mesh_dt, cache.mesh)
245246
return cache
246247
end

0 commit comments

Comments
 (0)