diff --git a/docs/src/globalerrorcontrol/GlobalDiffEq.md b/docs/src/globalerrorcontrol/GlobalDiffEq.md index 0999e40f80..f2ae509fd0 100644 --- a/docs/src/globalerrorcontrol/GlobalDiffEq.md +++ b/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 ``` diff --git a/lib/GlobalDiffEq/Project.toml b/lib/GlobalDiffEq/Project.toml index a20f4a3f5d..3087070381 100644 --- a/lib/GlobalDiffEq/Project.toml +++ b/lib/GlobalDiffEq/Project.toml @@ -1,30 +1,59 @@ name = "GlobalDiffEq" uuid = "1d72d19b-84cc-4cb7-a099-7cbdb9ccc67c" authors = ["Chris Rackauckas 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"] diff --git a/lib/GlobalDiffEq/README.md b/lib/GlobalDiffEq/README.md index 7a881a22ff..9ff142ddaf 100644 --- a/lib/GlobalDiffEq/README.md +++ b/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. 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..9260563706 100644 --- a/lib/GlobalDiffEq/src/GlobalDiffEq.jl +++ b/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 diff --git a/lib/GlobalDiffEq/src/adjoint.jl b/lib/GlobalDiffEq/src/adjoint.jl new file mode 100644 index 0000000000..d2d17bca89 --- /dev/null +++ b/lib/GlobalDiffEq/src/adjoint.jl @@ -0,0 +1,394 @@ +""" + 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 + +""" + 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), + 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/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/correction.jl b/lib/GlobalDiffEq/src/correction.jl new file mode 100644 index 0000000000..904ab77d6b --- /dev/null +++ b/lib/GlobalDiffEq/src/correction.jl @@ -0,0 +1,129 @@ +""" + GlobalDefectCorrection(alg; gtol=nothing, correction_alg=alg, maxiters=6, + safety=0.8, companion_abstol, companion_reltol) + +Wrap an ODE algorithm with global error estimation and control based on +solving for the correction (defect correction). After a dense forward solve +producing the interpolant `P(t)`, the global error `ε(t) ≈ y(t) - P(t)` is +estimated by integrating the nonlinear companion ODE + +```math +ε'(t) = f(P(t) + ε(t), p, t) - P'(t), \\qquad ε(t_0) = 0, +``` + +which requires no Jacobian, only extra evaluations of `f` on perturbed +arguments. This is the "solving for the correction" technique of Zadunaisky +(1976) and Dormand, Duckers and Prince (1984), in the continuous per-step +dense-output form of Dormand, Lockyer, McGorrigan and Prince (1989); its +validity under standard local-error step control was proven by Calvo, Higham, +Montijano and Rández (1996). + +Set the requested absolute endpoint error with `gtol`: + +```julia +solve(prob, GlobalDefectCorrection(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). +`correction_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 + + - P. E. Zadunaisky, On the estimation of errors propagated in the numerical + integration of ordinary differential equations, Numerische Mathematik 27 + (1976). + - J. R. Dormand, R. R. Duckers and P. J. Prince, Global error estimation with + Runge-Kutta methods, IMA Journal of Numerical Analysis 4 (1984). + - J. R. Dormand, M. A. Lockyer, N. E. McGorrigan and P. J. Prince, Global + error estimation with Runge-Kutta triples, Computers & Mathematics with + Applications 18 (1989). +""" +struct GlobalDefectCorrection{A, CA, G, V} <: GlobalDiffEqAlgorithm + alg::A + correction_alg::CA + gtol::G + options::V +end + +function GlobalDefectCorrection( + alg; + correction_alg = alg, + 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 GlobalDefectCorrection(alg, correction_alg, gtol, options) +end + +SciMLBase.allows_arbitrary_number_types(::GlobalDefectCorrection) = false +SciMLBase.allowscomplex(::GlobalDefectCorrection) = false +SciMLBase.isautodifferentiable(::GlobalDefectCorrection) = false + +function _correction_rhs(sol) + prob = sol.prob + foop = _oop_rhs(prob) + return let sol = sol, foop = foop, p = prob.p + function (ε, _, t) + u = sol(t, continuity = :right) + du = sol(t, Val{1}, continuity = :right) + return foop(u + ε, p, t) - du + end + end +end + +""" + global_error_estimate(prob, alg::GlobalDefectCorrection; 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 solving for the correction of the dense numerical solution. `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::GlobalDefectCorrection, 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( + _correction_rhs, "GlobalDefectCorrection", + prob, alg.alg, alg.correction_alg, args...; + abstol, reltol, companion_abstol, companion_reltol, kwargs... + ) +end + +function SciMLBase.__solve( + prob::SciMLBase.AbstractODEProblem, alg::GlobalDefectCorrection, args...; + abstol = something(alg.gtol, 1.0e-6), + reltol = something(alg.gtol, 1.0e-3), + kwargs... + ) + alg.gtol === nothing && throw( + ArgumentError("GlobalDefectCorrection requires a positive `gtol` constructor keyword") + ) + _validate_tolerances(abstol, reltol, "local") + haskey(kwargs, :callback) && + throw(ArgumentError("GlobalDefectCorrection 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/src/glee/algorithms.jl b/lib/GlobalDiffEq/src/glee/algorithms.jl new file mode 100644 index 0000000000..eb6145a98a --- /dev/null +++ b/lib/GlobalDiffEq/src/glee/algorithms.jl @@ -0,0 +1,90 @@ +# 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 + +""" + 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) + +_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/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..c6abd3e1f1 --- /dev/null +++ b/lib/GlobalDiffEq/src/glee/tableaus.jl @@ -0,0 +1,189 @@ +# 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 + +""" + 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/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/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/companion_tests.jl b/lib/GlobalDiffEq/test/companion_tests.jl new file mode 100644 index 0000000000..b3b58ead39 --- /dev/null +++ b/lib/GlobalDiffEq/test/companion_tests.jl @@ -0,0 +1,68 @@ +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()), GlobalDefectCorrection(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()), GlobalDefectCorrection(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), + GlobalDefectCorrection(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, GlobalDefectCorrection) + 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/glee_tests.jl b/lib/GlobalDiffEq/test/glee_tests.jl new file mode 100644 index 0000000000..1ce77d6eb4 --- /dev/null +++ b/lib/GlobalDiffEq/test/glee_tests.jl @@ -0,0 +1,111 @@ +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 + 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(), 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] + @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(), 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] + 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 "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( + (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 + @test SciMLBase.alg_order(MM5GEE()) == 5 +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..5e8338bdda 100644 --- a/lib/GlobalDiffEq/test/runtests.jl +++ b/lib/GlobalDiffEq/test/runtests.jl @@ -11,6 +11,15 @@ 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 + @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