GlobalDiffEq: add GlobalDefectCorrection global error estimation and control

docs/src/globalerrorcontrol/GlobalDiffEq.md

@@ -14,12 +14,63 @@ To use these methods:
 using GlobalDiffEq
 ```
 
-[`GlobalRichardson`](@ref) wraps any fixed-step method in global Richardson
-extrapolation over whole solves, interpreting `abstol` and `reltol` as global
-tolerances. It is the most robust and most expensive option.
+## Choosing a method
+
+  - For a solution accompanied by a running, asymptotically correct estimate of
+    its global error at every time point, use the global-error-estimating
+    solvers [`GLEE24`](@ref), [`GLEE35`](@ref) (Constantinescu 2016), or the
+    Dormand-Prince-based [`MM5GEE`](@ref) (Makazaga and Murua 2003). These cost
+    only a few extra stages per step over a plain method of the same order and
+    require nothing beyond the right-hand side `f`.
+  - To *control* the endpoint global error to a tolerance `gtol`, wrap any
+    adaptive solver in [`GlobalErrorTransport`](@ref) (linearized
+    error-transport equation, Jacobian-vector products via automatic
+    differentiation), [`GlobalDefectCorrection`](@ref) (solving for the
+    correction; no Jacobian needed), or [`GlobalAdjoint`](@ref) (adjoint-based,
+    for endpoint functionals; requires SciMLSensitivity and QuadGK to be
+    loaded). Each solves the problem, estimates the endpoint global error, and
+    tightens the local tolerances until the requested global tolerance is met.
+  - [`GlobalRichardson`](@ref) wraps any fixed-step method in global Richardson
+    extrapolation over whole solves, interpreting `abstol` and `reltol` as
+    global tolerances. It is the most robust and most expensive option.
+
+For example, solving with a controlled endpoint global error of `1e-8`:
+
+```julia
+using GlobalDiffEq, OrdinaryDiffEqTsit5
+
+function lorenz!(du, u, p, t)
+    du[1] = 10.0(u[2] - u[1])
+    du[2] = u[1] * (28.0 - u[3]) - u[2]
+    du[3] = u[1] * u[2] - (8 / 3) * u[3]
+end
+prob = ODEProblem(lorenz!, [1.0; 0.0; 0.0], (0.0, 10.0))
+sol = solve(prob, GlobalDefectCorrection(Tsit5(); gtol = 1.0e-8))
+```
+
+or solving while tracking the global error along the trajectory:
+
+```julia
+sol = solve(prob, GLEE35(); abstol = 1.0e-8, reltol = 1.0e-8)
+errs = global_error_estimate(sol)  # global error estimate at every sol.t
+```
+
+## Global-error-estimating solvers
+
+```@docs
+GLEE23
+GLEE24
+GLEE35
+MM5GEE
+global_error_estimate
+```
 
 ## Global error controlling wrappers
 
 ```@docs
 GlobalRichardson
+GlobalErrorTransport
+GlobalDefectCorrection
+GlobalAdjoint
+adjoint_error_estimate
 ```

lib/GlobalDiffEq/Project.toml

@@ -1,30 +1,59 @@
 name = "GlobalDiffEq"
 uuid = "1d72d19b-84cc-4cb7-a099-7cbdb9ccc67c"
 authors = ["Chris Rackauckas <accounts@chrisrackauckas.com> and contributors"]
-version = "1.2.1"
+version = "1.7.0"
 
 [deps]
+ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
+DifferentiationInterface = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63"
+FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
+OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8"
 OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Richardson = "708f8203-808e-40c0-ba2d-98a6953ed40d"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
+SciMLStructures = "53ae85a6-f571-4167-b2af-e1d143709226"
+
+[weakdeps]
+QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
+SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1"
+
+[extensions]
+GlobalDiffEqSciMLSensitivityExt = ["SciMLSensitivity", "QuadGK"]
 
 [sources]
 DiffEqBase = {path = "../DiffEqBase"}
+OrdinaryDiffEqCore = {path = "../OrdinaryDiffEqCore"}
 OrdinaryDiffEqTsit5 = {path = "../OrdinaryDiffEqTsit5"}
 
 [compat]
+ADTypes = "1"
 DiffEqBase = "7"
+DifferentiationInterface = "0.6, 0.7"
+FastBroadcast = "1.3"
+ForwardDiff = "0.10, 1"
 LinearAlgebra = "1"
+MuladdMacro = "0.2.1"
+OrdinaryDiffEqCore = "4"
 OrdinaryDiffEqSSPRK = "2"
 OrdinaryDiffEqTsit5 = "2"
 PrecompileTools = "1.1"
+QuadGK = "2.11"
+Random = "1"
+RecursiveArrayTools = "4.2.0"
 Reexport = "0.2, 1.0"
 Richardson = "1.2"
 SafeTestsets = "0.0.1, 0.1"
 SciMLBase = "3"
+SciMLSensitivity = "7.116"
+SciMLStructures = "1.10"
 SciMLTesting = "2.1"
 Test = "1"
 julia = "1.10"
@@ -33,9 +62,13 @@ julia = "1.10"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 OrdinaryDiffEqSSPRK = "669c94d9-1f4b-4b64-b377-1aa079aa2388"
 OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
+QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
+SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1"
 SciMLTesting = "09d9d899-5365-40a9-917a-5f67fddea283"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["LinearAlgebra", "OrdinaryDiffEqSSPRK", "OrdinaryDiffEqTsit5", "SafeTestsets", "SciMLTesting", "Test"]
+test = ["LinearAlgebra", "OrdinaryDiffEqSSPRK", "OrdinaryDiffEqTsit5", "QuadGK", "Random", "RecursiveArrayTools", "SafeTestsets", "SciMLSensitivity", "SciMLTesting", "Test"]

lib/GlobalDiffEq/README.md

@@ -12,6 +12,20 @@ provides:
 
   - `GlobalRichardson`: global Richardson extrapolation of whole solves of any
     fixed-step method, interpreting `abstol`/`reltol` as global tolerances.
+  - `GLEE23`, `GLEE24`, `GLEE35`: explicit general linear methods with built-in
+    global error estimation (Constantinescu 2016), which propagate the solution
+    together with an asymptotically correct estimate of its global error.
+  - `MM5GEE`: the Makazaga-Murua (2003) Dormand-Prince-based order-5 scheme
+    with a cheap built-in global error estimate.
+  - `GlobalErrorTransport`: global error estimation and `gtol`-based control by
+    integrating the linearized error-transport equation driven by the
+    dense-output defect (Shampine 1986, Berzins 1988, Lang-Verwer 2007).
+  - `GlobalDefectCorrection`: global error estimation and control by solving
+    for the correction (Zadunaisky 1976, Dormand-Duckers-Prince 1984/1989),
+    requiring no Jacobian.
+  - `GlobalAdjoint`: adjoint-based a posteriori endpoint error estimation and
+    control (Cao and Petzold 2004), available as a package extension when
+    SciMLSensitivity and QuadGK are loaded.
 
 See the [OrdinaryDiffEq.jl documentation](https://docs.sciml.ai/OrdinaryDiffEq/stable/)
 for the full API reference.

lib/GlobalDiffEq/ext/GlobalDiffEqSciMLSensitivityExt.jl

@@ -0,0 +1,65 @@
+module GlobalDiffEqSciMLSensitivityExt
+
+import GlobalDiffEq, LinearAlgebra, QuadGK, SciMLBase, SciMLSensitivity
+
+function GlobalDiffEq._default_quadrature_sensealg()
+    return SciMLSensitivity.QuadratureAdjoint(autojacvec = true)
+end
+
+GlobalDiffEq._is_quadrature_adjoint(::SciMLSensitivity.QuadratureAdjoint) = true
+
+function GlobalDiffEq._adjoint_solution(
+        sol, sensealg::SciMLSensitivity.QuadratureAdjoint, adjoint_alg, direction;
+        abstol, reltol
+    )
+    terminal_gradient! = let direction = direction
+        function (out, u, p, t, i)
+            copyto!(out, direction)
+            return nothing
+        end
+    end
+    terminal_time = sol.prob.tspan[2]
+    adjoint_prob = SciMLSensitivity.ODEAdjointProblem(
+        sol, sensealg, adjoint_alg, [terminal_time], terminal_gradient!
+    )
+    adjoint_sol = SciMLBase.solve(
+        adjoint_prob, adjoint_alg;
+        abstol, reltol, dense = true, save_everystep = true
+    )
+    SciMLBase.successful_retcode(adjoint_sol) ||
+        throw(ErrorException("the adjoint solve failed with retcode $(adjoint_sol.retcode)"))
+    return adjoint_sol
+end
+
+function GlobalDiffEq._defect_projection(sol, adjoint_sol; abstol, reltol)
+    prob = sol.prob
+    isinplace = SciMLBase.isinplace(prob)
+
+    function integrand(t)
+        u = sol(t, continuity = :right)
+        du = sol(t, Val{1}, continuity = :right)
+        rhs = if isinplace
+            value = similar(u)
+            prob.f(value, u, prob.p, t)
+            value
+        else
+            prob.f(u, prob.p, t)
+        end
+        lambda = adjoint_sol(t, continuity = :right)
+        return LinearAlgebra.dot(lambda, du - rhs)
+    end
+
+    projection = zero(eltype(prob.u0))
+    for i in 1:(length(sol.t) - 1)
+        left = sol.t[i]
+        right = sol.t[i + 1]
+        left == right && continue
+        interval_projection, _ = QuadGK.quadgk(
+            integrand, left, right; atol = abstol, rtol = reltol
+        )
+        projection += interval_projection
+    end
+    return projection
+end
+
+end

lib/GlobalDiffEq/src/GlobalDiffEq.jl

@@ -3,14 +3,31 @@ module GlobalDiffEq
 using Reexport: @reexport
 @reexport using DiffEqBase
 
-import OrdinaryDiffEqTsit5, Richardson, SciMLBase
+import ADTypes, DifferentiationInterface, ForwardDiff, LinearAlgebra,
+    OrdinaryDiffEqCore, OrdinaryDiffEqTsit5, Random, RecursiveArrayTools,
+    Richardson, SciMLBase, SciMLStructures
+import DiffEqBase: initialize!, calculate_residuals, calculate_residuals!
+import OrdinaryDiffEqCore: perform_step!, @cache
+using FastBroadcast: FastBroadcast, @..
+using MuladdMacro: MuladdMacro, @muladd
 using PrecompileTools: @setup_workload, @compile_workload
 
 abstract type GlobalDiffEqAlgorithm <: SciMLBase.AbstractODEAlgorithm end
 
 include("richardson.jl")
-
-export GlobalRichardson
+include("adjoint.jl")
+include("companion.jl")
+include("transport.jl")
+include("correction.jl")
+include("glee/tableaus.jl")
+include("glee/algorithms.jl")
+include("glee/solve.jl")
+include("glee/caches.jl")
+include("glee/perform_step.jl")
+
+export GlobalAdjoint, GlobalRichardson, adjoint_error_estimate
+export GlobalDefectCorrection, GlobalErrorTransport
+export GLEE23, GLEE24, GLEE35, MM5GEE, global_error_estimate
 
 @setup_workload begin
     # Simple test ODE: exponential decay du/dt = -u