From d1c20535af49f052d65eaea742b6e2e532f7ec3c Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Sat, 18 Jul 2026 00:01:02 -0400 Subject: [PATCH 1/4] Add adjoint-based global error control as a SciMLSensitivity extension Carries SciML/GlobalDiffEq.jl#56 (GlobalAdjoint, adjoint_error_estimate; Cao and Petzold 2004, closes SciML/GlobalDiffEq.jl#4) into the monorepo, restructured as the package extension GlobalDiffEqSciMLSensitivityExt on SciMLSensitivity + QuadGK weakdeps. sensealg defaults to nothing and resolves to QuadratureAdjoint(autojacvec = true) when the extension loads; solving without the extension raises an instructive error. Public API and semantics otherwise match the original PR. Co-Authored-By: Chris Rackauckas --- docs/src/globalerrorcontrol/GlobalDiffEq.md | 8 + lib/GlobalDiffEq/Project.toml | 21 +- lib/GlobalDiffEq/README.md | 3 + .../ext/GlobalDiffEqSciMLSensitivityExt.jl | 65 +++ lib/GlobalDiffEq/src/GlobalDiffEq.jl | 6 +- lib/GlobalDiffEq/src/adjoint.jl | 383 ++++++++++++++++++ lib/GlobalDiffEq/test/adjoint_tests.jl | 61 +++ lib/GlobalDiffEq/test/qa/Project.toml | 4 + lib/GlobalDiffEq/test/qa/qa.jl | 17 +- lib/GlobalDiffEq/test/runtests.jl | 3 + 10 files changed, 563 insertions(+), 8 deletions(-) create mode 100644 lib/GlobalDiffEq/ext/GlobalDiffEqSciMLSensitivityExt.jl create mode 100644 lib/GlobalDiffEq/src/adjoint.jl create mode 100644 lib/GlobalDiffEq/test/adjoint_tests.jl diff --git a/docs/src/globalerrorcontrol/GlobalDiffEq.md b/docs/src/globalerrorcontrol/GlobalDiffEq.md index 0999e40f80..f055dc6205 100644 --- a/docs/src/globalerrorcontrol/GlobalDiffEq.md +++ b/docs/src/globalerrorcontrol/GlobalDiffEq.md @@ -18,8 +18,16 @@ using GlobalDiffEq extrapolation over whole solves, interpreting `abstol` and `reltol` as global tolerances. It is the most robust and most expensive option. +To *control* the endpoint global error to a tolerance `gtol`, wrap any +adaptive solver in [`GlobalAdjoint`](@ref) (adjoint-based, for endpoint +functionals; requires SciMLSensitivity and QuadGK to be loaded). It solves the +problem, estimates the endpoint global error, and tightens the local +tolerances until the requested global tolerance is met. + ## Global error controlling wrappers ```@docs GlobalRichardson +GlobalAdjoint +adjoint_error_estimate ``` diff --git a/lib/GlobalDiffEq/Project.toml b/lib/GlobalDiffEq/Project.toml index a20f4a3f5d..25a57c1734 100644 --- a/lib/GlobalDiffEq/Project.toml +++ b/lib/GlobalDiffEq/Project.toml @@ -1,15 +1,25 @@ name = "GlobalDiffEq" uuid = "1d72d19b-84cc-4cb7-a099-7cbdb9ccc67c" authors = ["Chris Rackauckas and contributors"] -version = "1.2.1" +version = "1.3.0" [deps] DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a" PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" 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"} @@ -21,10 +31,14 @@ LinearAlgebra = "1" OrdinaryDiffEqSSPRK = "2" OrdinaryDiffEqTsit5 = "2" PrecompileTools = "1.1" +QuadGK = "2.11" +Random = "1" 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 +47,12 @@ 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" 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", "SafeTestsets", "SciMLSensitivity", "SciMLTesting", "Test"] diff --git a/lib/GlobalDiffEq/README.md b/lib/GlobalDiffEq/README.md index 7a881a22ff..30c8dba47d 100644 --- a/lib/GlobalDiffEq/README.md +++ b/lib/GlobalDiffEq/README.md @@ -12,6 +12,9 @@ provides: - `GlobalRichardson`: global Richardson extrapolation of whole solves of any fixed-step method, interpreting `abstol`/`reltol` as global tolerances. + - `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. diff --git a/lib/GlobalDiffEq/ext/GlobalDiffEqSciMLSensitivityExt.jl b/lib/GlobalDiffEq/ext/GlobalDiffEqSciMLSensitivityExt.jl new file mode 100644 index 0000000000..9e09e74830 --- /dev/null +++ b/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 diff --git a/lib/GlobalDiffEq/src/GlobalDiffEq.jl b/lib/GlobalDiffEq/src/GlobalDiffEq.jl index 77fd942f54..73e7f7c712 100644 --- a/lib/GlobalDiffEq/src/GlobalDiffEq.jl +++ b/lib/GlobalDiffEq/src/GlobalDiffEq.jl @@ -3,14 +3,16 @@ module GlobalDiffEq using Reexport: @reexport @reexport using DiffEqBase -import OrdinaryDiffEqTsit5, Richardson, SciMLBase +import LinearAlgebra, OrdinaryDiffEqTsit5, Random, Richardson, SciMLBase, + SciMLStructures using PrecompileTools: @setup_workload, @compile_workload abstract type GlobalDiffEqAlgorithm <: SciMLBase.AbstractODEAlgorithm end include("richardson.jl") +include("adjoint.jl") -export GlobalRichardson +export GlobalAdjoint, GlobalRichardson, adjoint_error_estimate @setup_workload begin # Simple test ODE: exponential decay du/dt = -u diff --git a/lib/GlobalDiffEq/src/adjoint.jl b/lib/GlobalDiffEq/src/adjoint.jl new file mode 100644 index 0000000000..32ffddd7d4 --- /dev/null +++ b/lib/GlobalDiffEq/src/adjoint.jl @@ -0,0 +1,383 @@ +""" + GlobalAdjoint(alg; gtol=nothing, adjoint_alg=alg, sensealg=nothing, samples=2, + rng=Random.default_rng(), maxiters=6, safety=0.8) + +Wrap an ODE algorithm with adjoint-based global error estimation and control. +Set the requested absolute endpoint error with the `gtol` constructor keyword: + +```julia +using SciMLSensitivity, QuadGK +solve(prob, GlobalAdjoint(Tsit5(); gtol = 1.0e-6)) +``` + +The algorithm solves the ODE with local tolerances, estimates the 2-norm of the +endpoint error using `samples` orthogonal random projections, and tightens the +local tolerances until the estimate is at most `gtol`. The estimate is +probabilistic for systems with more than one state; the default of two samples +is accurate within a factor of ten with probability greater than 99% under the +small-sample estimate in Cao and Petzold (2004). + +Omit `gtol` when using the algorithm only with [`adjoint_error_estimate`](@ref). +`maxiters` limits the number of forward-solve refinements, while `safety` +controls each local-tolerance reduction. The `adjoint_abstol`, +`adjoint_reltol`, `quadrature_abstol`, and `quadrature_reltol` constructor +keywords control the numerical adjoints and defect quadrature. + +The adjoint machinery lives in a package extension: both SciMLSensitivity and +QuadGK must be loaded to solve with `GlobalAdjoint`. `sensealg` must be a +`SciMLSensitivity.QuadratureAdjoint`; the default (`nothing`) resolves to +`SciMLSensitivity.QuadratureAdjoint(autojacvec = true)`, which uses +SciMLSensitivity's ForwardDiff-based state Jacobian construction. Parameters +are presented to SciMLSensitivity as fixed data because this estimator does not +differentiate with respect to them. `adjoint_alg` selects the solver for the +reverse-time adjoint problems. + +This implementation supports forward-time, standard-mass-matrix ODEs with +real vector states and no callbacks. The wrapped solver must provide a dense +first-derivative interpolation. +""" +struct GlobalAdjoint{A, AA, S, R, G, V} <: GlobalDiffEqAlgorithm + alg::A + adjoint_alg::AA + sensealg::S + samples::Int + rng::R + gtol::G + options::V +end + +function GlobalAdjoint( + alg; + adjoint_alg = alg, + sensealg = nothing, + samples = 2, + rng = Random.default_rng(), + gtol = nothing, + maxiters = 6, + safety = 0.8, + adjoint_abstol = gtol === nothing ? 1.0e-10 : min(gtol / 100, 1.0e-10), + adjoint_reltol = gtol === nothing ? 1.0e-8 : min(gtol / 100, 1.0e-8), + quadrature_abstol = gtol === nothing ? 0.0 : zero(gtol), + quadrature_reltol = 1.0e-6 + ) + samples isa Integer || throw(ArgumentError("samples must be an integer")) + samples > 0 || throw(ArgumentError("samples must be positive")) + sensealg === nothing || _is_quadrature_adjoint(sensealg) || + throw(ArgumentError(_ADJOINT_SENSEALG_MESSAGE)) + (gtol === nothing || _positive_finite_real(gtol)) || + throw(ArgumentError("gtol must be a positive finite real number")) + maxiters isa Integer && maxiters > 0 || + throw(ArgumentError("maxiters must be a positive integer")) + safety isa Real && isfinite(safety) && 0 < safety < 1 || + throw(ArgumentError("safety must be between zero and one")) + _validate_tolerances(adjoint_abstol, adjoint_reltol, "adjoint") + _validate_tolerances(quadrature_abstol, quadrature_reltol, "quadrature") + options = (; + maxiters = Int(maxiters), safety, + adjoint_abstol, adjoint_reltol, + quadrature_abstol, quadrature_reltol, + ) + return GlobalAdjoint(alg, adjoint_alg, sensealg, Int(samples), rng, gtol, options) +end + +# Extension hooks: GlobalDiffEqSciMLSensitivityExt adds methods for these when +# both SciMLSensitivity and QuadGK are loaded. +function _adjoint_solution end +function _defect_projection end +function _default_quadrature_sensealg end +_is_quadrature_adjoint(sensealg) = false + +const _ADJOINT_SENSEALG_MESSAGE = "sensealg must be a SciMLSensitivity.QuadratureAdjoint, \ + and both SciMLSensitivity and QuadGK must be loaded for GlobalAdjoint \ + (run `using SciMLSensitivity, QuadGK`)" + +_adjoint_ext_loaded() = !isempty(methods(_default_quadrature_sensealg)) + +function _require_adjoint_ext() + _adjoint_ext_loaded() || throw( + ArgumentError( + "GlobalAdjoint requires the SciMLSensitivity extension: run " * + "`using SciMLSensitivity, QuadGK` before solving" + ) + ) + return nothing +end + +function _resolve_sensealg(alg::GlobalAdjoint) + _require_adjoint_ext() + sensealg = alg.sensealg === nothing ? _default_quadrature_sensealg() : alg.sensealg + _is_quadrature_adjoint(sensealg) || throw(ArgumentError(_ADJOINT_SENSEALG_MESSAGE)) + return sensealg +end + +_positive_finite_real(value) = value isa Real && isfinite(value) && value > 0 + +function _validate_tolerances(abstol, reltol, name) + abstol isa Real && isfinite(abstol) && abstol >= 0 || + throw(ArgumentError("$(name)_abstol must be a nonnegative finite real number")) + reltol isa Real && isfinite(reltol) && reltol >= 0 || + throw(ArgumentError("$(name)_reltol must be a nonnegative finite real number")) + iszero(abstol) && iszero(reltol) && + throw(ArgumentError("$(name)_abstol and $(name)_reltol cannot both be zero")) + return nothing +end + +SciMLBase.allows_arbitrary_number_types(::GlobalAdjoint) = false +SciMLBase.allowscomplex(::GlobalAdjoint) = false +SciMLBase.isautodifferentiable(::GlobalAdjoint) = false + +struct _FixedParameter{P, T} + value::P + tunables::Vector{T} +end + +SciMLStructures.isscimlstructure(::_FixedParameter) = true +SciMLStructures.ismutablescimlstructure(::_FixedParameter) = false +SciMLStructures.hasportion(::SciMLStructures.AbstractPortion, ::_FixedParameter) = false +SciMLStructures.hasportion(::SciMLStructures.Tunable, ::_FixedParameter) = true + +function SciMLStructures.canonicalize(::SciMLStructures.Tunable, p::_FixedParameter) + repack = _ -> p + return p.tunables, repack, true +end + +function SciMLStructures.canonicalize( + ::SciMLStructures.AbstractPortion, ::_FixedParameter + ) + return nothing, nothing, nothing +end + +function SciMLStructures.replace( + ::SciMLStructures.Tunable, p::_FixedParameter, tunables + ) + isempty(tunables) || throw(DimensionMismatch("fixed parameters have no tunable values")) + return p +end + +function _fixed_parameter_problem(prob) + _validate_adjoint_problem(prob) + fixed_parameter = _FixedParameter(prob.p, eltype(prob.u0)[]) + original_f = SciMLBase.unwrapped_f(prob.f) + wrapped_f = if SciMLBase.isinplace(prob) + let f = original_f + (du, u, p, t) -> f(du, u, p.value, t) + end + else + let f = original_f + (u, p, t) -> f(u, p.value, t) + end + end + full_specialized_f = SciMLBase.ODEFunction{ + SciMLBase.isinplace(prob), SciMLBase.FullSpecialize, + }(wrapped_f) + return SciMLBase.remake(prob; f = full_specialized_f, p = fixed_parameter) +end + +function _validate_adjoint_problem(prob) + prob.u0 isa AbstractVector{<:AbstractFloat} || + throw(ArgumentError("GlobalAdjoint requires a real floating-point vector state")) + isempty(prob.u0) && throw(ArgumentError("GlobalAdjoint requires a nonempty state")) + prob.tspan[1] < prob.tspan[2] || + throw(ArgumentError("GlobalAdjoint requires a forward-time ODEProblem")) + prob.f.mass_matrix == LinearAlgebra.I || + throw(ArgumentError("GlobalAdjoint currently requires the standard mass matrix")) + problem_kwargs = values(prob.kwargs) + if haskey(problem_kwargs, :callback) && problem_kwargs.callback !== nothing + throw(ArgumentError("GlobalAdjoint does not currently support callbacks")) + end + return nothing +end + +function _validate_adjoint_solution(sol) + SciMLBase.successful_retcode(sol) || + throw(ErrorException("the forward solve failed with retcode $(sol.retcode)")) + _validate_adjoint_problem(sol.prob) + length(sol.t) >= 2 || + throw(ArgumentError("the forward solution must save at least its two endpoints")) + return nothing +end + +function _orthogonal_directions(u0, requested_samples, rng) + sample_count = min(requested_samples, length(u0)) + T = eltype(u0) + directions = Vector{Vector{T}}() + sizehint!(directions, sample_count) + + for _ in 1:sample_count + accepted = false + for _ in 1:10 + direction = T[Base.randn(rng) for _ in eachindex(u0)] + for _ in 1:2, previous in directions + direction .-= LinearAlgebra.dot(previous, direction) .* previous + end + direction_norm = LinearAlgebra.norm(direction) + if direction_norm > sqrt(eps(T)) + push!(directions, direction ./ direction_norm) + accepted = true + break + end + end + accepted || error("failed to generate orthogonal random directions") + end + + return directions +end + +function _sphere_projection_expectation(dimension, ::Type{T}) where {T} + dimension > 0 || throw(ArgumentError("dimension must be positive")) + if dimension == 1 + return one(T) + end + + expectation = isodd(dimension) ? one(T) : T(2) / T(pi) + first_numerator = isodd(dimension) ? 1 : 2 + for numerator in first_numerator:2:(dimension - 2) + expectation *= T(numerator) / T(numerator + 1) + end + return expectation +end + +function _adjoint_error_estimate( + sol, alg, directions; + adjoint_abstol, adjoint_reltol, quadrature_abstol, quadrature_reltol + ) + sensealg = _resolve_sensealg(alg) + projections = map(directions) do direction + adjoint_sol = _adjoint_solution( + sol, sensealg, alg.adjoint_alg, direction; + abstol = adjoint_abstol, reltol = adjoint_reltol + ) + _defect_projection( + sol, adjoint_sol; + abstol = quadrature_abstol, reltol = quadrature_reltol + ) + end + T = eltype(sol.prob.u0) + sample_count = length(directions) + factor = _sphere_projection_expectation(sample_count, T) / + _sphere_projection_expectation(length(sol.prob.u0), T) + return factor * LinearAlgebra.norm(projections) +end + +""" + adjoint_error_estimate(prob, alg::GlobalAdjoint; abstol=1e-6, reltol=1e-3, kwargs...) + +Solve an ODE problem and estimate the 2-norm of its global error at the final +time. The estimate evaluates the dense interpolation defect and projects it +through reverse-time adjoints constructed by SciMLSensitivity. Both +SciMLSensitivity and QuadGK must be loaded to enable the implementation. + +`abstol` and `reltol` control the forward solve. Keyword arguments +`adjoint_abstol`, `adjoint_reltol`, `quadrature_abstol`, and `quadrature_reltol` +control the numerical adjoint solves and defect quadrature. The random +directions come from `alg.rng`; their number is `alg.samples`, capped at the +state dimension. +""" +function adjoint_error_estimate( + prob::SciMLBase.AbstractODEProblem, alg::GlobalAdjoint, args...; + abstol = 1.0e-6, + reltol = 1.0e-3, + adjoint_abstol = alg.options.adjoint_abstol, + adjoint_reltol = alg.options.adjoint_reltol, + quadrature_abstol = alg.options.quadrature_abstol, + quadrature_reltol = alg.options.quadrature_reltol, + kwargs... + ) + _require_adjoint_ext() + haskey(kwargs, :callback) && + throw(ArgumentError("adjoint_error_estimate does not currently support callbacks")) + estimation_prob = _fixed_parameter_problem(prob) + solve_kwargs = merge( + (; kwargs...), + (; + dense = true, + save_everystep = true, + save_start = true, + save_end = true, + saveat = (), + save_idxs = nothing, + ) + ) + sol = SciMLBase.solve( + estimation_prob, alg.alg, args...; abstol, reltol, solve_kwargs... + ) + _validate_adjoint_solution(sol) + directions = _orthogonal_directions(prob.u0, alg.samples, alg.rng) + return _adjoint_error_estimate( + sol, alg, directions; + adjoint_abstol, + adjoint_reltol, + quadrature_abstol, + quadrature_reltol + ) +end + +function SciMLBase.__solve( + prob::SciMLBase.AbstractODEProblem, alg::GlobalAdjoint, args...; + abstol = something(alg.gtol, 1.0e-6), + reltol = something(alg.gtol, 1.0e-3), + kwargs... + ) + alg.gtol === nothing && + throw(ArgumentError("GlobalAdjoint requires a positive `gtol` constructor keyword")) + _require_adjoint_ext() + _validate_tolerances(abstol, reltol, "local") + haskey(kwargs, :callback) && + throw(ArgumentError("GlobalAdjoint does not currently support callbacks")) + + gtol = alg.gtol + options = alg.options + estimation_prob = _fixed_parameter_problem(prob) + directions = _orthogonal_directions(prob.u0, alg.samples, alg.rng) + local_abstol = abstol + local_reltol = reltol + trial_kwargs = merge( + (; kwargs...), + (; + dense = true, + save_everystep = true, + save_start = true, + save_end = true, + saveat = (), + save_idxs = nothing, + ) + ) + last_estimate = oftype(float(gtol), Inf) + + for _ in 1:options.maxiters + trial_sol = SciMLBase.solve( + estimation_prob, alg.alg, args...; + abstol = local_abstol, reltol = local_reltol, trial_kwargs... + ) + _validate_adjoint_solution(trial_sol) + last_estimate = _adjoint_error_estimate( + trial_sol, alg, directions; + adjoint_abstol = options.adjoint_abstol, + adjoint_reltol = options.adjoint_reltol, + quadrature_abstol = options.quadrature_abstol, + quadrature_reltol = options.quadrature_reltol + ) + isfinite(last_estimate) || + throw(ErrorException("the adjoint global error estimate is not finite")) + + if last_estimate <= gtol + return SciMLBase.solve( + prob, alg.alg, args...; + abstol = local_abstol, reltol = local_reltol, kwargs... + ) + end + + tolerance_scale = min(0.5, options.safety * gtol / last_estimate) + local_abstol *= tolerance_scale + local_reltol *= tolerance_scale + iszero(local_abstol) && iszero(local_reltol) && + throw(ErrorException("local tolerances underflowed during global error control")) + end + + throw( + ErrorException( + "failed to meet gtol=$(gtol) after $(options.maxiters) iterations; " * + "last estimated global error was $(last_estimate)" + ) + ) +end diff --git a/lib/GlobalDiffEq/test/adjoint_tests.jl b/lib/GlobalDiffEq/test/adjoint_tests.jl new file mode 100644 index 0000000000..0ea49ecac0 --- /dev/null +++ b/lib/GlobalDiffEq/test/adjoint_tests.jl @@ -0,0 +1,61 @@ +using GlobalDiffEq, OrdinaryDiffEqTsit5, Random +using SciMLSensitivity, QuadGK +using Test +import SciMLBase + +@testset "Adjoint error estimation and control" begin + function linear!(du, u, p, t) + du[1] = p * u[1] + return nothing + end + + rate = 2.0 + tspan = (0.0, 2.0) + linear_prob = ODEProblem(linear!, [1.0], tspan, rate) + estimator_alg = GlobalAdjoint(Tsit5(); samples = 1, rng = Xoshiro(123)) + + coarse_sol = solve( + linear_prob, Tsit5(); + abstol = 1.0e-4, reltol = 1.0e-4, dense = true, save_everystep = true + ) + estimated_error = adjoint_error_estimate( + linear_prob, estimator_alg; abstol = 1.0e-4, reltol = 1.0e-4 + ) + actual_error = abs(coarse_sol.u[end][1] - exp(rate * tspan[2])) + + @test estimated_error > 0 + @test estimated_error / actual_error ≈ 1 rtol = 0.15 + + gtol = 1.0e-6 + controlled_sol = solve( + linear_prob, + GlobalAdjoint(Tsit5(); gtol, samples = 1, rng = Xoshiro(321)); + abstol = 1.0e-3, reltol = 1.0e-3 + ) + @test abs(controlled_sol.u[end][1] - exp(rate * tspan[2])) <= gtol + @test controlled_sol.prob.p === rate + + forcing = t -> 0.1sin(t) + forced(u, p, t) = [-u[1] + p(t)] + + forced_prob = ODEProblem(forced, [0.0], (0.0, 1.0), forcing) + callable_parameter_estimate = adjoint_error_estimate( + forced_prob, GlobalAdjoint(Tsit5(); samples = 1, rng = Xoshiro(456)); + abstol = 1.0e-4, reltol = 1.0e-4 + ) + @test isfinite(callable_parameter_estimate) + @test callable_parameter_estimate >= 0 + + adjoint_alg = GlobalAdjoint(Tsit5()) + @test !SciMLBase.allows_arbitrary_number_types(adjoint_alg) + @test !SciMLBase.allowscomplex(adjoint_alg) + @test !SciMLBase.isautodifferentiable(adjoint_alg) + @test_throws ArgumentError GlobalAdjoint(Tsit5(); samples = 0) + @test_throws ArgumentError GlobalAdjoint(Tsit5(); gtol = -1.0) + @test_throws ArgumentError GlobalAdjoint(Tsit5(); sensealg = :not_a_sensealg) + @test_throws ArgumentError solve(linear_prob, adjoint_alg) + + callback = DiscreteCallback((u, t, integrator) -> false, integrator -> nothing) + callback_prob = ODEProblem(linear!, [1.0], tspan, rate; callback) + @test_throws ArgumentError adjoint_error_estimate(callback_prob, adjoint_alg) +end diff --git a/lib/GlobalDiffEq/test/qa/Project.toml b/lib/GlobalDiffEq/test/qa/Project.toml index 5245dd4b41..195725c8f6 100644 --- a/lib/GlobalDiffEq/test/qa/Project.toml +++ b/lib/GlobalDiffEq/test/qa/Project.toml @@ -4,6 +4,8 @@ GlobalDiffEq = "1d72d19b-84cc-4cb7-a099-7cbdb9ccc67c" JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b" OrdinaryDiffEqSSPRK = "669c94d9-1f4b-4b64-b377-1aa079aa2388" OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a" +QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" +SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1" SciMLTesting = "09d9d899-5365-40a9-917a-5f67fddea283" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" @@ -17,5 +19,7 @@ GlobalDiffEq = "1" JET = "0.9, 0.10, 0.11" OrdinaryDiffEqSSPRK = "2" OrdinaryDiffEqTsit5 = "2" +QuadGK = "2.11" +SciMLSensitivity = "7.116" SciMLTesting = "2.1" Test = "1" diff --git a/lib/GlobalDiffEq/test/qa/qa.jl b/lib/GlobalDiffEq/test/qa/qa.jl index c86904e621..eb1a27404e 100644 --- a/lib/GlobalDiffEq/test/qa/qa.jl +++ b/lib/GlobalDiffEq/test/qa/qa.jl @@ -1,6 +1,7 @@ using SciMLTesting, GlobalDiffEq, Test using JET using OrdinaryDiffEqTsit5, OrdinaryDiffEqSSPRK +using SciMLSensitivity, QuadGK # `@reexport using DiffEqBase` republishes DiffEqBase's API; those names are # documented and rendered at their owning packages, not in the OrdinaryDiffEq @@ -16,6 +17,7 @@ reexported_names = Tuple( run_qa( GlobalDiffEq; + reexports_allow = union(public_api_names(DiffEqBase), (:DiffEqBase,)), explicit_imports = true, # GlobalDiffEq's rendered documentation lives in the monorepo docs, two # directories up from the sublibrary root. @@ -26,10 +28,17 @@ run_qa( ), ei_kwargs = (; all_qualified_accesses_are_public = (; - # `SciMLBase.__solve` is SciMLBase's internal solve entry point (not - # part of the public API); GlobalDiffEq overloads it via its owner - # SciMLBase. - ignore = (:__solve,), + ignore = ( + # `SciMLBase.__solve` is SciMLBase's internal solve entry point (not + # part of the public API); GlobalDiffEq overloads it via its owner + # SciMLBase. + :__solve, + # Extension hooks owned by GlobalDiffEq itself: the + # SciMLSensitivity extension adds methods to these internal + # functions of its parent package. + :_adjoint_solution, :_defect_projection, + :_default_quadrature_sensealg, :_is_quadrature_adjoint, + ), ), ), # `@reexport using DiffEqBase` deliberately reexports DiffEqBase's API, so diff --git a/lib/GlobalDiffEq/test/runtests.jl b/lib/GlobalDiffEq/test/runtests.jl index 1d8bdfd6a6..e9f910a69e 100644 --- a/lib/GlobalDiffEq/test/runtests.jl +++ b/lib/GlobalDiffEq/test/runtests.jl @@ -11,6 +11,9 @@ run_tests(; @safetestset "Algorithm traits forwarding" begin include("algorithm_traits_tests.jl") end + @safetestset "Adjoint error estimation and control" begin + include("adjoint_tests.jl") + end return @safetestset "BigFloat support" begin include("bigfloat_tests.jl") end From 6b4606cdbad14136063cb7b74d269cb6a576979d Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Sat, 18 Jul 2026 00:12:21 -0400 Subject: [PATCH 2/4] Add GLEE23/GLEE24/GLEE35 global-error-estimating GLM solvers MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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 --- docs/src/globalerrorcontrol/GlobalDiffEq.md | 48 +++++-- lib/GlobalDiffEq/Project.toml | 14 +- lib/GlobalDiffEq/README.md | 3 + lib/GlobalDiffEq/src/GlobalDiffEq.jl | 14 +- lib/GlobalDiffEq/src/glee/algorithms.jl | 70 ++++++++++ lib/GlobalDiffEq/src/glee/caches.jl | 52 +++++++ lib/GlobalDiffEq/src/glee/perform_step.jl | 131 ++++++++++++++++++ lib/GlobalDiffEq/src/glee/solve.jl | 106 ++++++++++++++ lib/GlobalDiffEq/src/glee/tableaus.jl | 144 ++++++++++++++++++++ lib/GlobalDiffEq/test/glee_tests.jl | 99 ++++++++++++++ lib/GlobalDiffEq/test/runtests.jl | 3 + 11 files changed, 671 insertions(+), 13 deletions(-) create mode 100644 lib/GlobalDiffEq/src/glee/algorithms.jl create mode 100644 lib/GlobalDiffEq/src/glee/caches.jl create mode 100644 lib/GlobalDiffEq/src/glee/perform_step.jl create mode 100644 lib/GlobalDiffEq/src/glee/solve.jl create mode 100644 lib/GlobalDiffEq/src/glee/tableaus.jl create mode 100644 lib/GlobalDiffEq/test/glee_tests.jl diff --git a/docs/src/globalerrorcontrol/GlobalDiffEq.md b/docs/src/globalerrorcontrol/GlobalDiffEq.md index f055dc6205..2c7794ff70 100644 --- a/docs/src/globalerrorcontrol/GlobalDiffEq.md +++ b/docs/src/globalerrorcontrol/GlobalDiffEq.md @@ -14,15 +14,45 @@ 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. - -To *control* the endpoint global error to a tolerance `gtol`, wrap any -adaptive solver in [`GlobalAdjoint`](@ref) (adjoint-based, for endpoint -functionals; requires SciMLSensitivity and QuadGK to be loaded). It solves the -problem, estimates the endpoint global error, and tightens the local -tolerances until the requested global tolerance is met. +## 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) or [`GLEE35`](@ref) (Constantinescu 2016). 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 [`GlobalAdjoint`](@ref) (adjoint-based, for endpoint + functionals; requires SciMLSensitivity and QuadGK to be loaded). It 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 while tracking the global error along the trajectory: + +```julia +using GlobalDiffEq + +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, 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 +global_error_estimate +``` ## Global error controlling wrappers diff --git a/lib/GlobalDiffEq/Project.toml b/lib/GlobalDiffEq/Project.toml index 25a57c1734..15149c2682 100644 --- a/lib/GlobalDiffEq/Project.toml +++ b/lib/GlobalDiffEq/Project.toml @@ -1,14 +1,18 @@ name = "GlobalDiffEq" uuid = "1d72d19b-84cc-4cb7-a099-7cbdb9ccc67c" authors = ["Chris Rackauckas and contributors"] -version = "1.3.0" +version = "1.4.0" [deps] DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" +FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898" 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" @@ -23,16 +27,21 @@ GlobalDiffEqSciMLSensitivityExt = ["SciMLSensitivity", "QuadGK"] [sources] DiffEqBase = {path = "../DiffEqBase"} +OrdinaryDiffEqCore = {path = "../OrdinaryDiffEqCore"} OrdinaryDiffEqTsit5 = {path = "../OrdinaryDiffEqTsit5"} [compat] DiffEqBase = "7" +FastBroadcast = "1.3" 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" @@ -49,10 +58,11 @@ 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", "QuadGK", "Random", "SafeTestsets", "SciMLSensitivity", "SciMLTesting", "Test"] +test = ["LinearAlgebra", "OrdinaryDiffEqSSPRK", "OrdinaryDiffEqTsit5", "QuadGK", "Random", "RecursiveArrayTools", "SafeTestsets", "SciMLSensitivity", "SciMLTesting", "Test"] diff --git a/lib/GlobalDiffEq/README.md b/lib/GlobalDiffEq/README.md index 30c8dba47d..72e63b51d5 100644 --- a/lib/GlobalDiffEq/README.md +++ b/lib/GlobalDiffEq/README.md @@ -12,6 +12,9 @@ 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. - `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. diff --git a/lib/GlobalDiffEq/src/GlobalDiffEq.jl b/lib/GlobalDiffEq/src/GlobalDiffEq.jl index 73e7f7c712..a3d9934bb3 100644 --- a/lib/GlobalDiffEq/src/GlobalDiffEq.jl +++ b/lib/GlobalDiffEq/src/GlobalDiffEq.jl @@ -3,16 +3,26 @@ module GlobalDiffEq using Reexport: @reexport @reexport using DiffEqBase -import LinearAlgebra, OrdinaryDiffEqTsit5, Random, Richardson, SciMLBase, - SciMLStructures +import 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") include("adjoint.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 GLEE23, GLEE24, GLEE35, global_error_estimate @setup_workload begin # Simple test ODE: exponential decay du/dt = -u diff --git a/lib/GlobalDiffEq/src/glee/algorithms.jl b/lib/GlobalDiffEq/src/glee/algorithms.jl new file mode 100644 index 0000000000..960119d737 --- /dev/null +++ b/lib/GlobalDiffEq/src/glee/algorithms.jl @@ -0,0 +1,70 @@ +# Explicit general linear methods with built-in global error estimation +# (GLEE methods) from Constantinescu (2016), doi:10.1137/15M1014633. +abstract type AbstractGLEEAlgorithm <: +OrdinaryDiffEqCore.OrdinaryDiffEqAdaptiveAlgorithm end + +const _GLEE_DOCS_SHARED = """ +GLEE methods propagate the solution `y` together with an asymptotically +correct estimate `ε` of its global (accumulated) error, at the cost of a few +extra stages per step. Solving any `ODEProblem` with a GLEE method produces a +solution whose states are `ArrayPartition`s: `sol.u[i].x[1]` is the solution +and `sol.u[i].x[2]` is the global error estimate at `sol.t[i]` (see +[`global_error_estimate`](@ref)). The per-step increment of `ε` is an +asymptotically correct local error estimate, which drives standard step-size +adaptivity, so local tolerances behave exactly as for ordinary adaptive +Runge-Kutta methods while the global error is estimated for free. + +Only explicit, mass-matrix-free ODEs are supported. The reference for the +methods and their theory is: + +Emil M. Constantinescu, *Estimating Global Errors in Time Stepping*, SIAM +Journal on Numerical Analysis 54(6), 2016. [arXiv:1503.05166](https://arxiv.org/abs/1503.05166) +""" + +""" + GLEE23() + +3-stage, second-order explicit general linear method with global error +estimation (Constantinescu 2016, eq. (4.6); PETSc's `TSGLEE23`). The cheapest +GLEE method: only the first decoupling condition holds, so prefer +[`GLEE24`](@ref) for long-time integration or mildly stiff problems. + +$(_GLEE_DOCS_SHARED) +""" +struct GLEE23 <: AbstractGLEEAlgorithm end + +""" + GLEE24() + +4-stage, second-order explicit general linear method with global error +estimation (Constantinescu 2016, eq. (A.3); PETSc's `TSGLEE24`). Satisfies +both decoupling conditions (`B·U` and `B·A·U` diagonal), which keeps the error +estimate faithful in long-time integration; this is the recommended +second-order GLEE method. + +$(_GLEE_DOCS_SHARED) +""" +struct GLEE24 <: AbstractGLEEAlgorithm end + +""" + GLEE35() + +5-stage, third-order explicit general linear method with global error +estimation (Constantinescu 2016, eq. (4.9); PETSc's `TSGLEE35`). Satisfies +both decoupling conditions and has a large negative-real-axis stability +region; this is the recommended third-order GLEE method. + +$(_GLEE_DOCS_SHARED) +""" +struct GLEE35 <: AbstractGLEEAlgorithm end + +SciMLBase.alg_order(::GLEE23) = 2 +SciMLBase.alg_order(::GLEE24) = 2 +SciMLBase.alg_order(::GLEE35) = 3 + +OrdinaryDiffEqCore.alg_adaptive_order(alg::AbstractGLEEAlgorithm) = + SciMLBase.alg_order(alg) + +_glee_tableau_for(::GLEE23, ::Type{T}, ::Type{T2}) where {T, T2} = GLEE23Tableau(T, T2) +_glee_tableau_for(::GLEE24, ::Type{T}, ::Type{T2}) where {T, T2} = GLEE24Tableau(T, T2) +_glee_tableau_for(::GLEE35, ::Type{T}, ::Type{T2}) where {T, T2} = GLEE35Tableau(T, T2) diff --git a/lib/GlobalDiffEq/src/glee/caches.jl b/lib/GlobalDiffEq/src/glee/caches.jl new file mode 100644 index 0000000000..6619a40768 --- /dev/null +++ b/lib/GlobalDiffEq/src/glee/caches.jl @@ -0,0 +1,52 @@ +@cache struct GLEECache{uType, yType, rateType, yRateType, uNoUnitsType, TabType} <: + OrdinaryDiffEqCore.OrdinaryDiffEqMutableCache + u::uType + uprev::uType + tmp::uType + fsalfirst::rateType + fsallast::rateType + ks::Vector{yRateType} + ytmp::yType + εloc::yType + atmp::uNoUnitsType + tab::TabType +end + +OrdinaryDiffEqCore.get_fsalfirstlast(cache::GLEECache, u) = + (cache.fsalfirst, cache.fsallast) + +struct GLEEConstantCache{TabType} <: OrdinaryDiffEqCore.OrdinaryDiffEqConstantCache + tab::TabType +end + +function OrdinaryDiffEqCore.alg_cache( + alg::AbstractGLEEAlgorithm, u, rate_prototype, ::Type{uEltypeNoUnits}, + ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, + dt, reltol, p, calck, + ::Val{true}, verbose + ) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} + tab = _glee_tableau_for( + alg, OrdinaryDiffEqCore.constvalue(uBottomEltypeNoUnits), + OrdinaryDiffEqCore.constvalue(tTypeNoUnits) + ) + y = u.x[1] + yrate = rate_prototype.x[1] + ks = [zero(yrate) for _ in 1:nstages(tab)] + return GLEECache( + u, uprev, zero(u), zero(rate_prototype), zero(rate_prototype), ks, + zero(y), zero(y), fill!(similar(y, uEltypeNoUnits), 0), tab + ) +end + +function OrdinaryDiffEqCore.alg_cache( + alg::AbstractGLEEAlgorithm, u, rate_prototype, ::Type{uEltypeNoUnits}, + ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, + dt, reltol, p, calck, + ::Val{false}, verbose + ) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} + tab = _glee_tableau_for( + alg, OrdinaryDiffEqCore.constvalue(uBottomEltypeNoUnits), + OrdinaryDiffEqCore.constvalue(tTypeNoUnits) + ) + return GLEEConstantCache(tab) +end diff --git a/lib/GlobalDiffEq/src/glee/perform_step.jl b/lib/GlobalDiffEq/src/glee/perform_step.jl new file mode 100644 index 0000000000..d5e70fae07 --- /dev/null +++ b/lib/GlobalDiffEq/src/glee/perform_step.jl @@ -0,0 +1,131 @@ +function initialize!(integrator, cache::GLEECache) + integrator.kshortsize = 2 + resize!(integrator.k, integrator.kshortsize) + integrator.k[1] = cache.fsalfirst + integrator.k[2] = cache.fsallast + integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t) + OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1) + return nothing +end + +@muladd function perform_step!(integrator, cache::GLEECache, repeat_step = false) + (; t, dt, uprev, u, f, p) = integrator + (; ks, ytmp, εloc, atmp, tab) = cache + inner = _glee_inner(f) + yprev = uprev.x[1] + εprev = uprev.x[2] + s = nstages(tab) + + for i in 1:s + if i == 1 && tab.stage1_fsal + copyto!(ks[1], integrator.fsalfirst.x[1]) + continue + end + @.. ytmp = tab.U[i, 1] * yprev + tab.U[i, 2] * εprev + for j in 1:(i - 1) + iszero(tab.A[i, j]) && continue + @.. ytmp = ytmp + dt * tab.A[i, j] * ks[j] + end + inner(ks[i], ytmp, p, t + tab.c[i] * dt) + OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1) + end + + ynew = u.x[1] + εnew = u.x[2] + fill!(εloc, zero(eltype(εloc))) + for j in 1:s + iszero(tab.B[2, j]) && continue + @.. εloc = εloc + dt * tab.B[2, j] * ks[j] + end + copyto!(ynew, yprev) + for j in 1:s + iszero(tab.B[1, j]) && continue + @.. ynew = ynew + dt * tab.B[1, j] * ks[j] + end + @.. εnew = εprev + εloc + + if tab.solution_stage == 0 + f(integrator.fsallast, u, p, t + dt) + OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1) + else + copyto!(integrator.fsallast.x[1], ks[tab.solution_stage]) + end + εslope = integrator.fsallast.x[2] + @.. εslope = εloc / dt + + if integrator.opts.adaptive + calculate_residuals!( + atmp, εloc, yprev, ynew, integrator.opts.abstol, + integrator.opts.reltol, integrator.opts.internalnorm, t + ) + OrdinaryDiffEqCore.set_EEst!(integrator, integrator.opts.internalnorm(atmp, t)) + end + return nothing +end + +function initialize!(integrator, cache::GLEEConstantCache) + integrator.kshortsize = 2 + integrator.k = typeof(integrator.k)(undef, integrator.kshortsize) + integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) + OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1) + integrator.fsallast = zero(integrator.fsalfirst) + integrator.k[1] = integrator.fsalfirst + integrator.k[2] = integrator.fsallast + return nothing +end + +@muladd function perform_step!(integrator, cache::GLEEConstantCache, repeat_step = false) + (; t, dt, uprev, u, f, p) = integrator + tab = cache.tab + inner = _glee_inner(f) + yprev = uprev.x[1] + εprev = uprev.x[2] + s = nstages(tab) + + k1 = if tab.stage1_fsal + integrator.fsalfirst.x[1] + else + Y1 = tab.U[1, 1] * yprev + tab.U[1, 2] * εprev + OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1) + inner(Y1, p, t + tab.c[1] * dt) + end + ks = Vector{typeof(k1)}(undef, s) + ks[1] = k1 + for i in 2:s + Y = tab.U[i, 1] * yprev + tab.U[i, 2] * εprev + for j in 1:(i - 1) + iszero(tab.A[i, j]) && continue + Y = Y + dt * tab.A[i, j] * ks[j] + end + ks[i] = inner(Y, p, t + tab.c[i] * dt) + OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1) + end + + ynew = yprev + εloc = zero(yprev) + for j in 1:s + ynew = ynew + dt * tab.B[1, j] * ks[j] + εloc = εloc + dt * tab.B[2, j] * ks[j] + end + εnew = εprev + εloc + integrator.u = RecursiveArrayTools.ArrayPartition(ynew, εnew) + + fy = if tab.solution_stage == 0 + OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1) + inner(ynew, p, t + dt) + else + ks[tab.solution_stage] + end + integrator.fsallast = RecursiveArrayTools.ArrayPartition(fy, εloc / dt) + integrator.k[1] = integrator.fsalfirst + integrator.k[2] = integrator.fsallast + + if integrator.opts.adaptive + atmp = calculate_residuals( + εloc, yprev, ynew, integrator.opts.abstol, + integrator.opts.reltol, integrator.opts.internalnorm, t + ) + OrdinaryDiffEqCore.set_EEst!(integrator, integrator.opts.internalnorm(atmp, t)) + end + return nothing +end diff --git a/lib/GlobalDiffEq/src/glee/solve.jl b/lib/GlobalDiffEq/src/glee/solve.jl new file mode 100644 index 0000000000..8e2d16e945 --- /dev/null +++ b/lib/GlobalDiffEq/src/glee/solve.jl @@ -0,0 +1,106 @@ +# The GLEE methods integrate a partitioned state (y, ε). Users pass a plain +# ODEProblem; __solve extends it to ArrayPartition form with this RHS wrapper, +# which applies the user's f to the solution partition. The ε partition has no +# ODE of its own (the GL method updates it directly), so its rate is zero +# wherever the wrapper is evaluated (initialization, dense output derivatives). +struct GLEEExtendedRHS{iip, F} + f::F +end + +function (rhs::GLEEExtendedRHS{true})(du, u, p, t) + rhs.f(du.x[1], u.x[1], p, t) + fill!(du.x[2], zero(eltype(du.x[2]))) + return nothing +end + +function (rhs::GLEEExtendedRHS{false})(u, p, t) + return RecursiveArrayTools.ArrayPartition(rhs.f(u.x[1], p, t), zero(u.x[1])) +end + +function _glee_inner(f) + rhs = f isa SciMLBase.AbstractODEFunction ? f.f : f + rhs isa GLEEExtendedRHS || throw( + ArgumentError( + "GLEE methods integrate plain ODEProblems through their own " * + "partitioned-state extension; do not pass a manually partitioned problem" + ) + ) + return rhs.f +end + +function _glee_extended_problem(prob) + prob.u0 isa RecursiveArrayTools.ArrayPartition && throw( + ArgumentError( + "GLEE methods construct their own partitioned (y, ε) state; " * + "pass the plain ODEProblem instead of an ArrayPartition state" + ) + ) + prob.u0 isa AbstractArray || + throw(ArgumentError("GLEE methods require an array state")) + prob.f.mass_matrix == LinearAlgebra.I || + throw(ArgumentError("GLEE methods require the standard mass matrix")) + iip = SciMLBase.isinplace(prob) + rhs = GLEEExtendedRHS{iip, typeof(prob.f)}(prob.f) + extended_f = SciMLBase.ODEFunction{iip, SciMLBase.FullSpecialize}(rhs) + u0 = RecursiveArrayTools.ArrayPartition(copy(prob.u0), zero(prob.u0)) + return SciMLBase.remake(prob; f = extended_f, u0 = u0) +end + +_is_glee_extended(prob) = prob.f.f isa GLEEExtendedRHS + +# init/solve on a plain ODEProblem transparently extend it to the partitioned +# (y, ε) state; the invoke dispatches into OrdinaryDiffEqCore's generic __init +# (its exact five-argument form, so the DiffEqBase default-algorithm catch-all +# cannot be selected). +function SciMLBase.__init( + prob::SciMLBase.AbstractODEProblem, alg::AbstractGLEEAlgorithm, + timeseries_init = (), ts_init = (), ks_init = (); + kwargs... + ) + extended_prob = _is_glee_extended(prob) ? prob : _glee_extended_problem(prob) + return invoke( + SciMLBase.__init, + Tuple{ + SciMLBase.AbstractODEProblem, + OrdinaryDiffEqCore.OrdinaryDiffEqAlgorithm, + Any, Any, Any, + }, + extended_prob, alg, timeseries_init, ts_init, ks_init; kwargs... + ) +end + +function SciMLBase.__solve( + prob::SciMLBase.AbstractODEProblem, alg::AbstractGLEEAlgorithm, args...; + kwargs... + ) + integrator = SciMLBase.__init(prob, alg, args...; kwargs...) + SciMLBase.solve!(integrator) + return integrator.sol +end + +""" + global_error_estimate(sol) + global_error_estimate(sol, i) + +Extract the global error estimate from a solution computed with a GLEE method +([`GLEE23`](@ref), [`GLEE24`](@ref), [`GLEE35`](@ref)). + +GLEE solutions carry the partitioned state `(y, ε)`: `sol.u[i].x[1]` is the +solution value and `sol.u[i].x[2]` the estimate of its global error at +`sol.t[i]`. `global_error_estimate(sol, i)` returns the error estimate at the +`i`-th saved time point and `global_error_estimate(sol)` returns the vector of +estimates at every saved time point. +""" +function global_error_estimate(sol::SciMLBase.AbstractODESolution) + return [global_error_estimate(sol, i) for i in eachindex(sol.u)] +end + +function global_error_estimate(sol::SciMLBase.AbstractODESolution, i::Integer) + u = sol.u[i] + u isa RecursiveArrayTools.ArrayPartition && length(u.x) == 2 || throw( + ArgumentError( + "global_error_estimate expects a solution produced by a GLEE method" + ) + ) + return u.x[2] +end diff --git a/lib/GlobalDiffEq/src/glee/tableaus.jl b/lib/GlobalDiffEq/src/glee/tableaus.jl new file mode 100644 index 0000000000..e0334d6e5d --- /dev/null +++ b/lib/GlobalDiffEq/src/glee/tableaus.jl @@ -0,0 +1,144 @@ +# GLEE (global-error-estimating general linear method) tableaus in the y-ε form +# of Constantinescu (2016), doi:10.1137/15M1014633 (arXiv:1503.05166): +# +# Y_i = U[i,1] y⁽ⁿ⁻¹⁾ + U[i,2] ε⁽ⁿ⁻¹⁾ + Δt Σ_j A[i,j] f(Y_j) +# y⁽ⁿ⁾ = y⁽ⁿ⁻¹⁾ + Δt Σ_j B[1,j] f(Y_j) +# ε⁽ⁿ⁾ = ε⁽ⁿ⁻¹⁾ + Δt Σ_j B[2,j] f(Y_j) +# +# ε⁽ⁿ⁾ is an asymptotically correct estimate of the global error at t_n, and the +# per-step increment ε⁽ⁿ⁾ − ε⁽ⁿ⁻¹⁾ is an asymptotically correct local error +# estimate for the order-p solution y⁽ⁿ⁾. +struct GLEETableau{MType, M2Type, BType, cType, T} + A::MType + U::M2Type + B::BType + c::cType + γ::T + order::Int + # Stage 1 equals the propagated solution register (U row [1, 0], zero A row), + # so its f evaluation can reuse fsalfirst. + stage1_fsal::Bool + # Index of a stage that equals the step's end solution y⁽ⁿ⁾ (FSAL-style), so + # its f evaluation doubles as fsallast; 0 when no such stage exists. + solution_stage::Int +end + +nstages(tab::GLEETableau) = length(tab.c) + +function _detect_solution_stage(A, U, B, c) + for i in 2:length(c) + isone(U[i, 1]) && iszero(U[i, 2]) || continue + isone(c[i]) || continue + all(j -> A[i, j] == B[1, j], 1:(i - 1)) || continue + all(iszero, @view B[1, i:end]) || continue + return i + end + return 0 +end + +# Convert a tableau from the paper's y-ỹ form to the y-ε form via Lemma 4.1 with +# T = [1 0; 1 1-γ]: U_yε[:,1] = U_yỹ*1, U_yε[:,2] = (1-γ) U_yỹ[:,2], +# B_yε = T⁻¹ B_yỹ. +function _yytilde_to_yeps(U::AbstractMatrix, B::AbstractMatrix, γ) + Uyε = hcat(U[:, 1] .+ U[:, 2], (1 - γ) .* U[:, 2]) + Byε = vcat(B[1:1, :], (B[2:2, :] .- B[1:1, :]) ./ (1 - γ)) + return Uyε, Byε +end + +function _glee_tableau(::Type{T}, ::Type{T2}, A, U, B, γ, order) where {T, T2} + c_exact = vec(sum(A, dims = 2)) + stage1_fsal = all(iszero, A[1, :]) && isone(U[1, 1]) && iszero(U[1, 2]) + solution_stage = _detect_solution_stage(A, U, B, c_exact) + return GLEETableau( + map(T, A), map(T, U), map(T, B), map(T2, c_exact), T(γ), order, + stage1_fsal, solution_stage + ) +end + +""" + GLEE23Tableau(T, T2) + +Tableau of the 3-stage, order-2 explicit GLEE method of Constantinescu (2016), +eq. (4.6) (PETSc's `TSGLEE23`). `γ = 0`, `B·U` diagonal. +""" +function GLEE23Tableau(::Type{T}, ::Type{T2}) where {T, T2} + A = [ + 0 0 0 + 1 0 0 + 1 // 4 1 // 4 0 + ] + U = [ + 1 0 + 1 10 + 1 -1 + ] + B = [ + 1 // 12 1 // 12 5 // 6 + 1 // 12 1 // 12 -1 // 6 + ] + return _glee_tableau(T, T2, A, U, B, 0, 2) +end + +""" + GLEE24Tableau(T, T2) + +Tableau of the 4-stage, order-2 explicit GLEE method of Constantinescu (2016), +eq. (A.3) (PETSc's `TSGLEE24`). `γ = 0`; both decoupling conditions (`B·U` and +`B·A·U` diagonal) hold, making it the robust second-order choice for long-time +integration. +""" +function GLEE24Tableau(::Type{T}, ::Type{T2}) where {T, T2} + A = [ + 0 0 0 0 + 3 // 4 0 0 0 + 1 // 4 29 // 60 0 0 + -21 // 44 145 // 44 -20 // 11 0 + ] + # y-ỹ form as printed in the paper + U = [ + 0 1 + 75 // 58 -17 // 58 + 0 1 + 0 1 + ] + B = [ + 109 // 275 58 // 75 -37 // 110 1 // 6 + 3 // 11 0 75 // 88 -1 // 8 + ] + Uyε, Byε = _yytilde_to_yeps(U, B, 0 // 1) + return _glee_tableau(T, T2, A, Uyε, Byε, 0, 2) +end + +""" + GLEE35Tableau(T, T2) + +Tableau of the 5-stage, order-3 explicit GLEE method of Constantinescu (2016), +eq. (4.9) (PETSc's `TSGLEE35`). `γ = 0`; both decoupling conditions hold. The +coefficients are exact rationals whose numerators exceed `Int64`, so the +tableau is constructed with `BigInt` rationals before conversion. +""" +function GLEE35Tableau(::Type{T}, ::Type{T2}) where {T, T2} + R(num, den) = big(num) // big(den) + Z = R(0, 1) + A = [ + Z Z Z Z Z; + R(-2169604947363702313, 24313474998937147335) Z Z Z Z; + R(46526746497697123895, 94116917485856474137) R(-10297879244026594958, 49199457603717988219) Z Z Z; + R(23364788935845982499, 87425311444725389446) R(-79205144337496116638, 148994349441340815519) R(40051189859317443782, 36487615018004984309) Z Z; + R(42089522664062539205, 124911313006412840286) R(-15074384760342762939, 137927286865289746282) R(-62274678522253371016, 125918573676298591413) R(13755475729852471739, 79257927066651693390) Z + ] + # y-ỹ form as printed in the paper + U = [ + R(70820309139834661559, 80863923579509469826) R(10043614439674808267, 80863923579509469826); + R(161694774978034105510, 106187653640211060371) R(-55507121337823045139, 106187653640211060371); + R(78486094644566264568, 88171030896733822981) R(9684936252167558413, 88171030896733822981); + R(65394922146334854435, 84570853840405479554) R(19175931694070625119, 84570853840405479554); + R(8607282770183754108, 108658046436496925911) R(100050763666313171803, 108658046436496925911) + ] + B = [ + R(61546696837458703723, 56982519523786160813) R(-55810892792806293355, 206957624151308356511) R(24061048952676379087, 158739347956038723465) R(3577972206874351339, 7599733370677197135) R(-59449832954780563947, 137360038685338563670); + R(-9738262186984159168, 99299082461487742983) R(-32797097931948613195, 61521565616362163366) R(42895514606418420631, 71714201188501437336) R(22608567633166065068, 55371917805607957003) R(94655809487476459565, 151517167160302729021) + ] + Uyε, Byε = _yytilde_to_yeps(U, B, Z) + return _glee_tableau(T, T2, A, Uyε, Byε, 0, 3) +end diff --git a/lib/GlobalDiffEq/test/glee_tests.jl b/lib/GlobalDiffEq/test/glee_tests.jl new file mode 100644 index 0000000000..f4ffc2e2c4 --- /dev/null +++ b/lib/GlobalDiffEq/test/glee_tests.jl @@ -0,0 +1,99 @@ +using GlobalDiffEq, OrdinaryDiffEqTsit5, LinearAlgebra +using RecursiveArrayTools: ArrayPartition +using Test +import SciMLBase + +# Unstable Prince42 problem: y' = y - sin(t) + cos(t), y(0) = 0, exact y = sin(t). +# Perturbations grow like e^t, making it the classic global-error-estimation test. +f_prince!(du, u, p, t) = (du[1] = u[1] - sin(t) + cos(t); nothing) +f_prince(u, p, t) = [u[1] - sin(t) + cos(t)] +const prince_tspan = (0.0, 2.0) +prince_exact(t) = sin(t) + +function glee_convergence(prob, alg, dts) + errs = Float64[] + est_errs = Float64[] + for dt in dts + sol = solve(prob, alg; dt = dt, adaptive = false) + err = prince_exact(prince_tspan[2]) - sol.u[end].x[1][1] + est = global_error_estimate(sol, length(sol.u))[1] + push!(errs, abs(err)) + push!(est_errs, abs(est - err)) + end + sol_order = log2(errs[end - 1] / errs[end]) + est_order = log2(est_errs[end - 1] / est_errs[end]) + return sol_order, est_order +end + +@testset "GLEE convergence and estimate accuracy" begin + dts = 2.0 .^ -(5:8) + for iip in (true, false) + prob = if iip + ODEProblem(f_prince!, [0.0], prince_tspan) + else + ODEProblem(f_prince, [0.0], prince_tspan) + end + for (alg, p) in ( + (GLEE23(), 2), (GLEE24(), 2), (GLEE35(), 3), + ) + sol_order, est_order = glee_convergence(prob, alg, dts) + @test sol_order ≈ p atol = 0.15 + # ε is an asymptotically correct global error estimate: it converges + # to the true error one order faster than the solution converges + @test est_order ≈ p + 1 atol = 0.2 + end + end +end + +@testset "GLEE estimate tracks the true error" begin + prob = ODEProblem(f_prince!, [0.0], prince_tspan) + for alg in (GLEE23(), GLEE24(), GLEE35()) + sol = solve(prob, alg; dt = 1 / 128, adaptive = false) + err = prince_exact(prince_tspan[2]) - sol.u[end].x[1][1] + est = global_error_estimate(sol)[end][1] + @test est / err ≈ 1 rtol = 0.1 + end +end + +@testset "GLEE adaptive stepping" begin + prob = ODEProblem(f_prince!, [0.0], prince_tspan) + for alg in (GLEE24(), GLEE35()) + sol = solve(prob, alg; abstol = 1.0e-8, reltol = 1.0e-8) + @test SciMLBase.successful_retcode(sol) + err = prince_exact(prince_tspan[2]) - sol.u[end].x[1][1] + est = global_error_estimate(sol)[end][1] + @test isfinite(est) + @test est / err ≈ 1 rtol = 0.6 + end +end + +@testset "GLEE integrator interface" begin + prob = ODEProblem(f_prince!, [0.0], prince_tspan) + integrator = init(prob, GLEE24(); dt = 1 / 64, adaptive = false) + @test integrator.u isa ArrayPartition + solve!(integrator) + sol = integrator.sol + @test SciMLBase.successful_retcode(sol) + err = prince_exact(prince_tspan[2]) - sol.u[end].x[1][1] + @test global_error_estimate(sol)[end][1] / err ≈ 1 rtol = 0.1 +end + +@testset "GLEE argument validation" begin + prob = ODEProblem(f_prince!, [0.0], prince_tspan) + partitioned_prob = ODEProblem( + (du, u, p, t) -> nothing, ArrayPartition([0.0], [0.0]), prince_tspan + ) + @test_throws ArgumentError solve(partitioned_prob, GLEE24(); dt = 0.1) + + mm_prob = ODEProblem( + ODEFunction(f_prince!; mass_matrix = [2.0;;]), [0.0], prince_tspan + ) + @test_throws ArgumentError solve(mm_prob, GLEE24(); dt = 0.1) + + sol = solve(prob, Tsit5()) + @test_throws ArgumentError global_error_estimate(sol, 1) + + @test SciMLBase.alg_order(GLEE23()) == 2 + @test SciMLBase.alg_order(GLEE24()) == 2 + @test SciMLBase.alg_order(GLEE35()) == 3 +end diff --git a/lib/GlobalDiffEq/test/runtests.jl b/lib/GlobalDiffEq/test/runtests.jl index e9f910a69e..df4cbe5018 100644 --- a/lib/GlobalDiffEq/test/runtests.jl +++ b/lib/GlobalDiffEq/test/runtests.jl @@ -14,6 +14,9 @@ run_tests(; @safetestset "Adjoint error estimation and control" begin include("adjoint_tests.jl") end + @safetestset "GLEE solvers" begin + include("glee_tests.jl") + end return @safetestset "BigFloat support" begin include("bigfloat_tests.jl") end From e362cb543b2ac8b2165c31e8b37cbf1d84bd1662 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Sat, 18 Jul 2026 00:19:49 -0400 Subject: [PATCH 3/4] Add the Makazaga-Murua MM5GEE global-error-estimating solver Adds the DOPRI5-based order-5 scheme with cheap global error estimation of Makazaga and Murua (BIT 43, 2003) as another tableau in the GLEE general linear framework (SciML/GlobalDiffEq.jl#6, the coefficient- complete published representative of the Dormand-Duckers-Prince RK triple family). Stages 1-7 are the standard DOPRI5 stages, three extra stages propagate an order-6 companion solution, and a generic solution-stage FSAL detection reuses stage 7's evaluation for fsallast, giving exactly 9 function evaluations per step (verified by an nf accounting test). The tableau was verified with exact rational arithmetic against the order and independency conditions. Verified locally on the unstable Prince42 problem: solution order 5, estimate accuracy order 6, est/true-error ratio 1.04, 9.01 f-evals/step. Co-Authored-By: Chris Rackauckas --- docs/src/globalerrorcontrol/GlobalDiffEq.md | 8 ++-- lib/GlobalDiffEq/Project.toml | 2 +- lib/GlobalDiffEq/README.md | 2 + lib/GlobalDiffEq/src/GlobalDiffEq.jl | 2 +- lib/GlobalDiffEq/src/glee/algorithms.jl | 20 +++++++++ lib/GlobalDiffEq/src/glee/tableaus.jl | 45 +++++++++++++++++++++ lib/GlobalDiffEq/test/glee_tests.jl | 16 +++++++- 7 files changed, 88 insertions(+), 7 deletions(-) diff --git a/docs/src/globalerrorcontrol/GlobalDiffEq.md b/docs/src/globalerrorcontrol/GlobalDiffEq.md index 2c7794ff70..d4dc8ec673 100644 --- a/docs/src/globalerrorcontrol/GlobalDiffEq.md +++ b/docs/src/globalerrorcontrol/GlobalDiffEq.md @@ -18,9 +18,10 @@ using GlobalDiffEq - 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) or [`GLEE35`](@ref) (Constantinescu 2016). 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`. + 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 [`GlobalAdjoint`](@ref) (adjoint-based, for endpoint functionals; requires SciMLSensitivity and QuadGK to be loaded). It solves @@ -51,6 +52,7 @@ errs = global_error_estimate(sol) # global error estimate at every sol.t GLEE23 GLEE24 GLEE35 +MM5GEE global_error_estimate ``` diff --git a/lib/GlobalDiffEq/Project.toml b/lib/GlobalDiffEq/Project.toml index 15149c2682..fa4f80330d 100644 --- a/lib/GlobalDiffEq/Project.toml +++ b/lib/GlobalDiffEq/Project.toml @@ -1,7 +1,7 @@ name = "GlobalDiffEq" uuid = "1d72d19b-84cc-4cb7-a099-7cbdb9ccc67c" authors = ["Chris Rackauckas and contributors"] -version = "1.4.0" +version = "1.5.0" [deps] DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" diff --git a/lib/GlobalDiffEq/README.md b/lib/GlobalDiffEq/README.md index 72e63b51d5..d90331a751 100644 --- a/lib/GlobalDiffEq/README.md +++ b/lib/GlobalDiffEq/README.md @@ -15,6 +15,8 @@ provides: - `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. - `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. diff --git a/lib/GlobalDiffEq/src/GlobalDiffEq.jl b/lib/GlobalDiffEq/src/GlobalDiffEq.jl index a3d9934bb3..04fa5172a8 100644 --- a/lib/GlobalDiffEq/src/GlobalDiffEq.jl +++ b/lib/GlobalDiffEq/src/GlobalDiffEq.jl @@ -22,7 +22,7 @@ include("glee/caches.jl") include("glee/perform_step.jl") export GlobalAdjoint, GlobalRichardson, adjoint_error_estimate -export GLEE23, GLEE24, GLEE35, global_error_estimate +export GLEE23, GLEE24, GLEE35, MM5GEE, global_error_estimate @setup_workload begin # Simple test ODE: exponential decay du/dt = -u diff --git a/lib/GlobalDiffEq/src/glee/algorithms.jl b/lib/GlobalDiffEq/src/glee/algorithms.jl index 960119d737..eb6145a98a 100644 --- a/lib/GlobalDiffEq/src/glee/algorithms.jl +++ b/lib/GlobalDiffEq/src/glee/algorithms.jl @@ -58,9 +58,28 @@ $(_GLEE_DOCS_SHARED) """ struct GLEE35 <: AbstractGLEEAlgorithm end +""" + MM5GEE() + +Fifth-order Dormand-Prince-based scheme with cheap global error estimation +(Makazaga and Murua, BIT Numerical Mathematics 43, 2003). Propagates the +standard DOPRI5 solution together with an order-6 companion solution through +three extra stages (9 function evaluations per step in total), whose +difference is an asymptotically correct global error estimate. This is the +recommended method of the Runge-Kutta triple family of Dormand, Duckers and +Prince, as it is the coefficient-complete published scheme of that type. + +$(_GLEE_DOCS_SHARED) + + - J. Makazaga and A. Murua, New Runge-Kutta based schemes for ODEs with + cheap global error estimation, BIT Numerical Mathematics 43 (2003). +""" +struct MM5GEE <: AbstractGLEEAlgorithm end + SciMLBase.alg_order(::GLEE23) = 2 SciMLBase.alg_order(::GLEE24) = 2 SciMLBase.alg_order(::GLEE35) = 3 +SciMLBase.alg_order(::MM5GEE) = 5 OrdinaryDiffEqCore.alg_adaptive_order(alg::AbstractGLEEAlgorithm) = SciMLBase.alg_order(alg) @@ -68,3 +87,4 @@ OrdinaryDiffEqCore.alg_adaptive_order(alg::AbstractGLEEAlgorithm) = _glee_tableau_for(::GLEE23, ::Type{T}, ::Type{T2}) where {T, T2} = GLEE23Tableau(T, T2) _glee_tableau_for(::GLEE24, ::Type{T}, ::Type{T2}) where {T, T2} = GLEE24Tableau(T, T2) _glee_tableau_for(::GLEE35, ::Type{T}, ::Type{T2}) where {T, T2} = GLEE35Tableau(T, T2) +_glee_tableau_for(::MM5GEE, ::Type{T}, ::Type{T2}) where {T, T2} = MM5GEETableau(T, T2) diff --git a/lib/GlobalDiffEq/src/glee/tableaus.jl b/lib/GlobalDiffEq/src/glee/tableaus.jl index e0334d6e5d..c6abd3e1f1 100644 --- a/lib/GlobalDiffEq/src/glee/tableaus.jl +++ b/lib/GlobalDiffEq/src/glee/tableaus.jl @@ -142,3 +142,48 @@ function GLEE35Tableau(::Type{T}, ::Type{T2}) where {T, T2} Uyε, Byε = _yytilde_to_yeps(U, B, Z) return _glee_tableau(T, T2, A, Uyε, Byε, 0, 3) end + +""" + MM5GEETableau(T, T2) + +Tableau of the Makazaga-Murua global-error-estimating scheme built on the +Dormand-Prince 5(4) method (Makazaga and Murua, BIT Numerical Mathematics 43, +2003, Table 4.1), rewritten in the y-ε general linear form: stages 1-7 are the +standard DOPRI5 stages (stage 7 the FSAL stage), stages 8-10 propagate the +order-6 companion solution `ȳ = y + ε` through the mixing coefficients +`U[i,2] = 1 - μ_i`, and the second output row is `b̄ - b`. +""" +function MM5GEETableau(::Type{T}, ::Type{T2}) where {T, T2} + R(num, den) = big(num) // big(den) + Z = R(0, 1) + b = [ + R(35, 384), Z, R(500, 1113), R(125, 192), R(-2187, 6784), R(11, 84), Z, + Z, Z, Z, + ] + A = [ + Z Z Z Z Z Z Z Z Z Z; + R(1, 5) Z Z Z Z Z Z Z Z Z; + R(3, 40) R(9, 40) Z Z Z Z Z Z Z Z; + R(44, 45) R(-56, 15) R(32, 9) Z Z Z Z Z Z Z; + R(19372, 6561) R(-25360, 2187) R(64448, 6561) R(-212, 729) Z Z Z Z Z Z; + R(9017, 3168) R(-355, 33) R(46732, 5247) R(49, 176) R(-5103, 18656) Z Z Z Z Z; + R(35, 384) Z R(500, 1113) R(125, 192) R(-2187, 6784) R(11, 84) Z Z Z Z; + R(26251126, 75292183) R(-30511879, 68834945) R(11490887, 155205387) R(700737845, 174891007) R(-5336, 941) R(5735, 1214) R(-2507, 898) Z Z Z; + R(-126276029, 115017392) R(153409379, 49308629) R(-107711621, 48274693) R(-675136779, 64711289) R(559269939, 36928210) R(-669687859, 52442748) R(193952703, 25738526) R(169021117, 130072535) Z Z; + R(89178409, 82486612) R(-275044175, 99029299) R(115406143, 68971088) R(140298385, 24130572) R(-344040692, 42025591) R(121333564, 17575013) R(-190380249, 47005513) R(-12078143, 165601005) R(56747365, 92317949) Z + ] + # U[:, 2] = 1 - μ_i: the weight of the companion register ȳ = y + ε in each + # stage; μ_i = 1 for the DOPRI5 stages 1-7. + one_minus_μ = [ + Z, Z, Z, Z, Z, Z, Z, + R(140719960, 143529893), R(941, 896), R(92493035, 95359057), + ] + U = hcat(fill(R(1, 1), 10), one_minus_μ) + b̄ = [ + R(56696811, 789712427), Z, R(-47431484, 279691831), R(72791025, 357831874), + R(17490085, 349505178), R(-66245097, 563676842), R(-24, 611), + R(40757463, 82884629), R(33159666, 111811519), R(42422453, 199331202), + ] + B = permutedims(hcat(b, b̄ .- b)) + return _glee_tableau(T, T2, A, U, B, 0, 5) +end diff --git a/lib/GlobalDiffEq/test/glee_tests.jl b/lib/GlobalDiffEq/test/glee_tests.jl index f4ffc2e2c4..1ce77d6eb4 100644 --- a/lib/GlobalDiffEq/test/glee_tests.jl +++ b/lib/GlobalDiffEq/test/glee_tests.jl @@ -42,12 +42,15 @@ end # to the true error one order faster than the solution converges @test est_order ≈ p + 1 atol = 0.2 end + sol_order, est_order = glee_convergence(prob, MM5GEE(), 2.0 .^ -(2:5)) + @test sol_order ≈ 5 atol = 0.35 + @test est_order > 5.5 end end @testset "GLEE estimate tracks the true error" begin prob = ODEProblem(f_prince!, [0.0], prince_tspan) - for alg in (GLEE23(), GLEE24(), GLEE35()) + for alg in (GLEE23(), GLEE24(), GLEE35(), MM5GEE()) sol = solve(prob, alg; dt = 1 / 128, adaptive = false) err = prince_exact(prince_tspan[2]) - sol.u[end].x[1][1] est = global_error_estimate(sol)[end][1] @@ -57,7 +60,7 @@ end @testset "GLEE adaptive stepping" begin prob = ODEProblem(f_prince!, [0.0], prince_tspan) - for alg in (GLEE24(), GLEE35()) + for alg in (GLEE24(), MM5GEE()) sol = solve(prob, alg; abstol = 1.0e-8, reltol = 1.0e-8) @test SciMLBase.successful_retcode(sol) err = prince_exact(prince_tspan[2]) - sol.u[end].x[1][1] @@ -78,6 +81,14 @@ end @test global_error_estimate(sol)[end][1] / err ≈ 1 rtol = 0.1 end +@testset "MM5GEE function evaluation count" begin + prob = ODEProblem(f_prince!, [0.0], prince_tspan) + sol = solve(prob, MM5GEE(); dt = 1 / 32, adaptive = false) + # 9 f-evals per step (FSAL first stage + solution stage reuse) plus the + # initialization evaluation + @test sol.stats.nf == 9 * (length(sol.t) - 1) + 1 +end + @testset "GLEE argument validation" begin prob = ODEProblem(f_prince!, [0.0], prince_tspan) partitioned_prob = ODEProblem( @@ -96,4 +107,5 @@ end @test SciMLBase.alg_order(GLEE23()) == 2 @test SciMLBase.alg_order(GLEE24()) == 2 @test SciMLBase.alg_order(GLEE35()) == 3 + @test SciMLBase.alg_order(MM5GEE()) == 5 end From 94843bc4c564d07404fd8a44c1dca36140662bce Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Sat, 18 Jul 2026 00:28:42 -0400 Subject: [PATCH 4/4] Add GlobalErrorTransport global error estimation and control MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Implements the linearized error-transport estimator (SciML/GlobalDiffEq.jl#5; Shampine 1986 / Berzins 1988 / Lang-Verwer 2007): after a dense forward solve, the companion linear ODE ε' = Jε + d(t) driven by the dense-output defect d(t) = f(P(t)) - P'(t) is integrated for the endpoint global error, with matrix-free Jacobian-vector products through DifferentiationInterface (autodiff accepts any ADTypes backend, default AutoForwardDiff). gtol-controlled solves tighten local tolerances via the shared refinement loop in companion.jl, mirroring GlobalAdjoint. global_error_estimate(prob, alg) is the standalone estimator, also added for GlobalAdjoint as an alias of adjoint_error_estimate. Verified locally on Lotka-Volterra vs a 1e-13 reference: estimate to true-error ratios 0.995-0.999 at tolerances 1e-4/1e-6, and gtol=1e-7 control yields 9.8e-9 final error. Co-Authored-By: Chris Rackauckas --- docs/src/globalerrorcontrol/GlobalDiffEq.md | 11 +- lib/GlobalDiffEq/Project.toml | 8 +- lib/GlobalDiffEq/README.md | 3 + lib/GlobalDiffEq/src/GlobalDiffEq.jl | 8 +- lib/GlobalDiffEq/src/adjoint.jl | 11 ++ lib/GlobalDiffEq/src/companion.jl | 143 ++++++++++++++++++++ lib/GlobalDiffEq/src/transport.jl | 136 +++++++++++++++++++ lib/GlobalDiffEq/test/companion_tests.jl | 67 +++++++++ lib/GlobalDiffEq/test/runtests.jl | 3 + 9 files changed, 383 insertions(+), 7 deletions(-) create mode 100644 lib/GlobalDiffEq/src/companion.jl create mode 100644 lib/GlobalDiffEq/src/transport.jl create mode 100644 lib/GlobalDiffEq/test/companion_tests.jl diff --git a/docs/src/globalerrorcontrol/GlobalDiffEq.md b/docs/src/globalerrorcontrol/GlobalDiffEq.md index d4dc8ec673..6f3f945940 100644 --- a/docs/src/globalerrorcontrol/GlobalDiffEq.md +++ b/docs/src/globalerrorcontrol/GlobalDiffEq.md @@ -23,10 +23,12 @@ using GlobalDiffEq 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 [`GlobalAdjoint`](@ref) (adjoint-based, for endpoint - functionals; requires SciMLSensitivity and QuadGK to be loaded). It solves - the problem, estimates the endpoint global error, and tightens the local - tolerances until the requested global tolerance is met. + adaptive solver in [`GlobalErrorTransport`](@ref) (linearized + error-transport equation, Jacobian-vector products via automatic + differentiation) 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. @@ -60,6 +62,7 @@ global_error_estimate ```@docs GlobalRichardson +GlobalErrorTransport GlobalAdjoint adjoint_error_estimate ``` diff --git a/lib/GlobalDiffEq/Project.toml b/lib/GlobalDiffEq/Project.toml index fa4f80330d..bc4e0c5651 100644 --- a/lib/GlobalDiffEq/Project.toml +++ b/lib/GlobalDiffEq/Project.toml @@ -1,11 +1,14 @@ name = "GlobalDiffEq" uuid = "1d72d19b-84cc-4cb7-a099-7cbdb9ccc67c" authors = ["Chris Rackauckas and contributors"] -version = "1.5.0" +version = "1.6.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" @@ -31,8 +34,11 @@ 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" diff --git a/lib/GlobalDiffEq/README.md b/lib/GlobalDiffEq/README.md index d90331a751..f478cf8fcd 100644 --- a/lib/GlobalDiffEq/README.md +++ b/lib/GlobalDiffEq/README.md @@ -17,6 +17,9 @@ provides: 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). - `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. diff --git a/lib/GlobalDiffEq/src/GlobalDiffEq.jl b/lib/GlobalDiffEq/src/GlobalDiffEq.jl index 04fa5172a8..9aa957f7bf 100644 --- a/lib/GlobalDiffEq/src/GlobalDiffEq.jl +++ b/lib/GlobalDiffEq/src/GlobalDiffEq.jl @@ -3,8 +3,9 @@ module GlobalDiffEq using Reexport: @reexport @reexport using DiffEqBase -import LinearAlgebra, OrdinaryDiffEqCore, OrdinaryDiffEqTsit5, Random, - RecursiveArrayTools, Richardson, SciMLBase, SciMLStructures +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, @.. @@ -15,6 +16,8 @@ abstract type GlobalDiffEqAlgorithm <: SciMLBase.AbstractODEAlgorithm end include("richardson.jl") include("adjoint.jl") +include("companion.jl") +include("transport.jl") include("glee/tableaus.jl") include("glee/algorithms.jl") include("glee/solve.jl") @@ -22,6 +25,7 @@ include("glee/caches.jl") include("glee/perform_step.jl") export GlobalAdjoint, GlobalRichardson, adjoint_error_estimate +export GlobalErrorTransport export GLEE23, GLEE24, GLEE35, MM5GEE, global_error_estimate @setup_workload begin diff --git a/lib/GlobalDiffEq/src/adjoint.jl b/lib/GlobalDiffEq/src/adjoint.jl index 32ffddd7d4..d2d17bca89 100644 --- a/lib/GlobalDiffEq/src/adjoint.jl +++ b/lib/GlobalDiffEq/src/adjoint.jl @@ -312,6 +312,17 @@ function adjoint_error_estimate( ) end +""" + global_error_estimate(prob, alg::GlobalAdjoint; kwargs...) + +Alias for [`adjoint_error_estimate`](@ref). +""" +function global_error_estimate( + prob::SciMLBase.AbstractODEProblem, alg::GlobalAdjoint, args...; kwargs... + ) + return adjoint_error_estimate(prob, alg, args...; kwargs...) +end + function SciMLBase.__solve( prob::SciMLBase.AbstractODEProblem, alg::GlobalAdjoint, args...; abstol = something(alg.gtol, 1.0e-6), diff --git a/lib/GlobalDiffEq/src/companion.jl b/lib/GlobalDiffEq/src/companion.jl new file mode 100644 index 0000000000..1b3e50b5ec --- /dev/null +++ b/lib/GlobalDiffEq/src/companion.jl @@ -0,0 +1,143 @@ +# Shared machinery for the defect-companion global error estimators +# (GlobalErrorTransport and GlobalDefectCorrection): both solve an auxiliary +# ODE for the global error driven by the defect of the dense numerical +# solution, d(t) = f(P(t)) - P'(t), where P is the solution interpolant. + +function _validate_estimation_problem(prob, name) + prob.u0 isa AbstractVector{<:AbstractFloat} || + throw(ArgumentError("$name requires a real floating-point vector state")) + isempty(prob.u0) && throw(ArgumentError("$name requires a nonempty state")) + prob.tspan[1] < prob.tspan[2] || + throw(ArgumentError("$name requires a forward-time ODEProblem")) + prob.f.mass_matrix == LinearAlgebra.I || + throw(ArgumentError("$name currently requires the standard mass matrix")) + problem_kwargs = values(prob.kwargs) + if haskey(problem_kwargs, :callback) && problem_kwargs.callback !== nothing + throw(ArgumentError("$name does not currently support callbacks")) + end + return nothing +end + +function _validate_estimation_solution(sol, name) + SciMLBase.successful_retcode(sol) || + throw(ErrorException("the forward solve failed with retcode $(sol.retcode)")) + length(sol.t) >= 2 || + throw(ArgumentError("the forward solution must save at least its two endpoints")) + return nothing +end + +# Out-of-place view of the (possibly in-place) user RHS, safe for dual numbers. +function _oop_rhs(prob) + f = SciMLBase.unwrapped_f(prob.f) + if SciMLBase.isinplace(prob) + return let f = f + function (u, p, t) + du = similar(u) + f(du, u, p, t) + return du + end + end + else + return f + end +end + +const _DENSE_SOLVE_KWARGS = (; + dense = true, + save_everystep = true, + save_start = true, + save_end = true, + saveat = (), + save_idxs = nothing, +) + +# Endpoint-error refinement loop shared by the companion estimators: solve, +# estimate the endpoint global error, and tighten the local tolerances until +# the estimate is at most gtol; then redo the solve with the user's original +# saving options. +function _refine_to_gtol( + estimator, prob, inner_alg, gtol, options, args...; + abstol, reltol, kwargs... + ) + local_abstol = abstol + local_reltol = reltol + last_estimate = oftype(float(gtol), Inf) + + for _ in 1:options.maxiters + last_estimate = estimator(local_abstol, local_reltol) + isfinite(last_estimate) || + throw(ErrorException("the global error estimate is not finite")) + + if last_estimate <= gtol + return SciMLBase.solve( + prob, inner_alg, args...; + abstol = local_abstol, reltol = local_reltol, kwargs... + ) + end + + tolerance_scale = min(0.5, options.safety * gtol / last_estimate) + local_abstol *= tolerance_scale + local_reltol *= tolerance_scale + iszero(local_abstol) && iszero(local_reltol) && + throw(ErrorException("local tolerances underflowed during global error control")) + end + + throw( + ErrorException( + "failed to meet gtol=$(gtol) after $(options.maxiters) iterations; " * + "last estimated global error was $(last_estimate)" + ) + ) +end + +function _companion_options(gtol, maxiters, safety, companion_abstol, companion_reltol) + maxiters isa Integer && maxiters > 0 || + throw(ArgumentError("maxiters must be a positive integer")) + safety isa Real && isfinite(safety) && 0 < safety < 1 || + throw(ArgumentError("safety must be between zero and one")) + (gtol === nothing || _positive_finite_real(gtol)) || + throw(ArgumentError("gtol must be a positive finite real number")) + _validate_tolerances(companion_abstol, companion_reltol, "companion") + return (; + maxiters = Int(maxiters), safety, + companion_abstol, companion_reltol, + ) +end + +# Dense forward solve + companion error solve, returning the endpoint global +# error 2-norm estimate. rhs_for(sol) builds the companion RHS from the dense +# forward solution. +function _companion_error_estimate( + rhs_for, name, prob, inner_alg, companion_alg, args...; + abstol, reltol, companion_abstol, companion_reltol, kwargs... + ) + haskey(kwargs, :callback) && + throw(ArgumentError("$name does not currently support callbacks")) + _validate_estimation_problem(prob, name) + solve_kwargs = merge((; kwargs...), _DENSE_SOLVE_KWARGS) + sol = SciMLBase.solve(prob, inner_alg, args...; abstol, reltol, solve_kwargs...) + _validate_estimation_solution(sol, name) + endpoint_error = _companion_endpoint( + rhs_for(sol), sol, companion_alg; + abstol = companion_abstol, reltol = companion_reltol + ) + return LinearAlgebra.norm(endpoint_error) +end + +# Solve the companion error ODE ε' = rhs(ε, t), ε(t0) = 0, over the forward +# solution's time span and return the endpoint error estimate ε(T). +function _companion_endpoint(rhs, sol, companion_alg; abstol, reltol) + ε0 = zero(sol.prob.u0) + companion_prob = SciMLBase.ODEProblem{false}(rhs, ε0, sol.prob.tspan) + companion_sol = SciMLBase.solve( + companion_prob, companion_alg; + abstol, reltol, dense = false, save_everystep = false, + save_start = false, save_end = true + ) + SciMLBase.successful_retcode(companion_sol) || throw( + ErrorException( + "the companion error solve failed with retcode $(companion_sol.retcode)" + ) + ) + return companion_sol.u[end] +end diff --git a/lib/GlobalDiffEq/src/transport.jl b/lib/GlobalDiffEq/src/transport.jl new file mode 100644 index 0000000000..bc3b4f97e2 --- /dev/null +++ b/lib/GlobalDiffEq/src/transport.jl @@ -0,0 +1,136 @@ +""" + GlobalErrorTransport(alg; gtol=nothing, transport_alg=alg, + autodiff=ADTypes.AutoForwardDiff(), maxiters=6, + safety=0.8, companion_abstol, companion_reltol) + +Wrap an ODE algorithm with global error estimation and control based on the +linearized error-transport (first variational) equation. After a dense forward +solve producing the interpolant `P(t)`, the global error `ε(t) ≈ y(t) - P(t)` +is estimated by integrating the linear companion ODE + +```math +ε'(t) = J(P(t), p, t) \\, ε(t) + d(t), \\qquad ε(t_0) = 0, +``` + +where `d(t) = f(P(t), p, t) - P'(t)` is the defect of the dense output and `J` +is the Jacobian of `f`, applied matrix-free through Jacobian-vector products +computed with `autodiff` (any `ADTypes` backend supported by +DifferentiationInterface). This is the error-transport approach of Shampine +(1986) and Berzins (1988), in the continuous defect-driven form summarized by +Lang and Verwer (2007). + +Set the requested absolute endpoint error with `gtol`: + +```julia +solve(prob, GlobalErrorTransport(Tsit5(); gtol = 1.0e-6)) +``` + +The solver then tightens local tolerances until the estimated endpoint global +error 2-norm is at most `gtol`, like [`GlobalAdjoint`](@ref). Omit `gtol` when +using the algorithm only with [`global_error_estimate`](@ref). `transport_alg` +selects the solver for the companion equation, and `companion_abstol` / +`companion_reltol` control its tolerances. + +This implementation supports forward-time, standard-mass-matrix ODEs with real +vector states and no callbacks. The wrapped solver must provide a dense +first-derivative interpolation. + +## References + + - L. F. Shampine, Global error estimation with one-step methods, Computers & + Mathematics with Applications 12A (1986). + - M. Berzins, Global error estimation in the method of lines for parabolic + equations, SIAM Journal on Scientific and Statistical Computing 9 (1988). + - J. Lang and J. Verwer, On global error estimation and control for initial + value problems, SIAM Journal on Scientific Computing 29 (2007). +""" +struct GlobalErrorTransport{A, TA, AD, G, V} <: GlobalDiffEqAlgorithm + alg::A + transport_alg::TA + autodiff::AD + gtol::G + options::V +end + +function GlobalErrorTransport( + alg; + transport_alg = alg, + autodiff = ADTypes.AutoForwardDiff(), + gtol = nothing, + maxiters = 6, + safety = 0.8, + companion_abstol = gtol === nothing ? 1.0e-10 : min(gtol / 100, 1.0e-10), + companion_reltol = gtol === nothing ? 1.0e-8 : min(gtol / 100, 1.0e-8) + ) + options = _companion_options( + gtol, maxiters, safety, companion_abstol, companion_reltol + ) + return GlobalErrorTransport(alg, transport_alg, autodiff, gtol, options) +end + +SciMLBase.allows_arbitrary_number_types(::GlobalErrorTransport) = false +SciMLBase.allowscomplex(::GlobalErrorTransport) = false +SciMLBase.isautodifferentiable(::GlobalErrorTransport) = false + +function _transport_rhs(sol, autodiff) + prob = sol.prob + foop = _oop_rhs(prob) + return let sol = sol, foop = foop, p = prob.p, backend = autodiff + function (ε, _, t) + u = sol(t, continuity = :right) + du = sol(t, Val{1}, continuity = :right) + defect = foop(u, p, t) - du + jv = only( + DifferentiationInterface.pushforward( + x -> foop(x, p, t), backend, u, (ε,) + ) + ) + return jv + defect + end + end +end + +""" + global_error_estimate(prob, alg::GlobalErrorTransport; abstol=1e-6, reltol=1e-3, kwargs...) + +Solve an ODE problem and estimate the 2-norm of its global error at the final +time by integrating the linearized error-transport equation driven by the +dense-output defect. `abstol` and `reltol` control the forward solve; +`companion_abstol` and `companion_reltol` control the companion error solve. +""" +function global_error_estimate( + prob::SciMLBase.AbstractODEProblem, alg::GlobalErrorTransport, args...; + abstol = 1.0e-6, + reltol = 1.0e-3, + companion_abstol = alg.options.companion_abstol, + companion_reltol = alg.options.companion_reltol, + kwargs... + ) + return _companion_error_estimate( + sol -> _transport_rhs(sol, alg.autodiff), "GlobalErrorTransport", + prob, alg.alg, alg.transport_alg, args...; + abstol, reltol, companion_abstol, companion_reltol, kwargs... + ) +end + +function SciMLBase.__solve( + prob::SciMLBase.AbstractODEProblem, alg::GlobalErrorTransport, args...; + abstol = something(alg.gtol, 1.0e-6), + reltol = something(alg.gtol, 1.0e-3), + kwargs... + ) + alg.gtol === nothing && throw( + ArgumentError("GlobalErrorTransport requires a positive `gtol` constructor keyword") + ) + _validate_tolerances(abstol, reltol, "local") + haskey(kwargs, :callback) && + throw(ArgumentError("GlobalErrorTransport does not currently support callbacks")) + estimator = (local_abstol, local_reltol) -> global_error_estimate( + prob, alg, args...; + abstol = local_abstol, reltol = local_reltol, kwargs... + ) + return _refine_to_gtol( + estimator, prob, alg.alg, alg.gtol, alg.options, args...; + abstol, reltol, kwargs... + ) +end diff --git a/lib/GlobalDiffEq/test/companion_tests.jl b/lib/GlobalDiffEq/test/companion_tests.jl new file mode 100644 index 0000000000..2021ff166c --- /dev/null +++ b/lib/GlobalDiffEq/test/companion_tests.jl @@ -0,0 +1,67 @@ +using GlobalDiffEq, OrdinaryDiffEqTsit5, LinearAlgebra +using Test +import SciMLBase + +lv!(du, u, p, t) = (du[1] = 1.5u[1] - u[1] * u[2]; du[2] = -3.0u[2] + u[1] * u[2]; nothing) +lv(u, p, t) = [1.5u[1] - u[1] * u[2], -3.0u[2] + u[1] * u[2]] +const lv_tspan = (0.0, 10.0) + +@testset "Companion global error estimators" begin + prob = ODEProblem(lv!, [1.0, 1.0], lv_tspan) + ref = solve(prob, Tsit5(); abstol = 1.0e-13, reltol = 1.0e-13) + + for alg in (GlobalErrorTransport(Tsit5()),) + for tol in (1.0e-4, 1.0e-6) + est = global_error_estimate(prob, alg; abstol = tol, reltol = tol) + sol = solve(prob, Tsit5(); abstol = tol, reltol = tol) + true_err = norm(sol.u[end] - ref.u[end]) + @test est / true_err ≈ 1 rtol = 0.1 + end + end + + # out-of-place problems + prob_oop = ODEProblem(lv, [1.0, 1.0], lv_tspan) + for alg in (GlobalErrorTransport(Tsit5()),) + est = global_error_estimate(prob_oop, alg; abstol = 1.0e-6, reltol = 1.0e-6) + sol = solve(prob_oop, Tsit5(); abstol = 1.0e-6, reltol = 1.0e-6) + true_err = norm(sol.u[end] - ref.u[end]) + @test est / true_err ≈ 1 rtol = 0.1 + end +end + +@testset "Companion global error control" begin + prob = ODEProblem(lv!, [1.0, 1.0], lv_tspan) + ref = solve(prob, Tsit5(); abstol = 1.0e-13, reltol = 1.0e-13) + gtol = 1.0e-7 + + for alg in ( + GlobalErrorTransport(Tsit5(); gtol), + ) + sol = solve(prob, alg; abstol = 1.0e-3, reltol = 1.0e-3) + @test SciMLBase.successful_retcode(sol) + @test norm(sol.u[end] - ref.u[end]) <= gtol + end +end + +@testset "Companion estimator argument validation" begin + prob = ODEProblem(lv!, [1.0, 1.0], lv_tspan) + + for Alg in (GlobalErrorTransport,) + alg = Alg(Tsit5()) + @test !SciMLBase.allows_arbitrary_number_types(alg) + @test !SciMLBase.allowscomplex(alg) + @test !SciMLBase.isautodifferentiable(alg) + @test_throws ArgumentError Alg(Tsit5(); gtol = -1.0) + @test_throws ArgumentError Alg(Tsit5(); maxiters = 0) + @test_throws ArgumentError Alg(Tsit5(); safety = 2.0) + # gtol is required to solve with the wrapper + @test_throws ArgumentError solve(prob, alg) + + callback = DiscreteCallback((u, t, integrator) -> false, integrator -> nothing) + callback_prob = ODEProblem(lv!, [1.0, 1.0], lv_tspan; callback) + @test_throws ArgumentError global_error_estimate(callback_prob, alg) + + backwards_prob = ODEProblem(lv!, [1.0, 1.0], (10.0, 0.0)) + @test_throws ArgumentError global_error_estimate(backwards_prob, alg) + end +end diff --git a/lib/GlobalDiffEq/test/runtests.jl b/lib/GlobalDiffEq/test/runtests.jl index df4cbe5018..5e8338bdda 100644 --- a/lib/GlobalDiffEq/test/runtests.jl +++ b/lib/GlobalDiffEq/test/runtests.jl @@ -17,6 +17,9 @@ run_tests(; @safetestset "GLEE solvers" begin include("glee_tests.jl") end + @safetestset "Companion error estimators" begin + include("companion_tests.jl") + end return @safetestset "BigFloat support" begin include("bigfloat_tests.jl") end