diff --git a/.github/workflows/changelog.yml b/.github/workflows/changelog.yml new file mode 100644 index 0000000..826527b --- /dev/null +++ b/.github/workflows/changelog.yml @@ -0,0 +1,12 @@ +name: Check Changelog +on: + pull_request: + +jobs: + Check-Changelog: + name: Check Changelog Action + runs-on: ubuntu-latest + steps: + - uses: tarides/changelog-check-action@v3 + with: + changelog: NEWS.md \ No newline at end of file diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 2633b84..6351584 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -11,7 +11,7 @@ jobs: runs-on: ${{ matrix.os }} strategy: matrix: - julia-version: ["1.6", "1.10", "~1.11.0-0"] + julia-version: ["1.10", "~1.11.0-0"] os: [ubuntu-latest, macOS-latest, windows-latest] steps: - uses: actions/checkout@v4 diff --git a/.github/workflows/documenter.yml b/.github/workflows/documenter.yml index 4ecdb1e..6cab182 100644 --- a/.github/workflows/documenter.yml +++ b/.github/workflows/documenter.yml @@ -14,7 +14,7 @@ jobs: - uses: actions/checkout@v4 - uses: julia-actions/setup-julia@latest with: - version: 1.6 + version: "1.10" - uses: julia-actions/julia-docdeploy@v1 env: PYTHON: "" diff --git a/Project.toml b/Project.toml index 4b9c67f..d2f475a 100644 --- a/Project.toml +++ b/Project.toml @@ -1,29 +1,35 @@ name = "ManifoldDiffEq" uuid = "1143c485-9b25-4e23-a65f-701df382ec90" authors = ["Seth Axen ", "Mateusz Baran ", "Ronny Bergmann ", "Antoine Levitt "] -version = "0.1.9" +version = "0.2.0" [deps] +Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697" +ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9" DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Manifolds = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e" ManifoldsBase = "3362f125-f0bb-47a3-aa74-596ffd7ef2fb" Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a" -OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" +OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8" RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" +SimpleUnPack = "ce78b400-467f-4804-87d8-8f486da07d0a" [compat] +Accessors = "0.1" +ConstructionBase = "1" DiffEqBase = "6" DoubleFloats = ">= 0.9.2" LinearAlgebra = "1.6" Manifolds = "0.10" ManifoldsBase = "0.15" Markdown = "1.6" -OrdinaryDiffEq = "6" +OrdinaryDiffEqCore = "1" RecursiveArrayTools = "2, 3" -SciMLBase = "1.0, 2 - 2.39" -julia = "1.6" +SciMLBase = "1, 2" +SimpleUnPack = "1" +julia = "1.10" [extras] DoubleFloats = "497a8b3b-efae-58df-a0af-a86822472b78" diff --git a/docs/make.jl b/docs/make.jl index 76debdb..e0fab47 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -21,6 +21,7 @@ makedocs( "Error estimation" => "error_estimation.md", "Notation" => "notation.md", "Examples" => "examples.md", + "Internals" => "internals.md", "References" => "references.md", ], plugins = [bib], diff --git a/docs/src/examples.md b/docs/src/examples.md index 3d363e5..0b0ee8d 100644 --- a/docs/src/examples.md +++ b/docs/src/examples.md @@ -2,7 +2,7 @@ We take a look at the simple example from -In the following code an ODE on a sphere is solved the introductionary example from the +In the following code an ODE on a sphere is solved the introductory example from the [lecture notes](https://www.unige.ch/~hairer/poly-sde-mani.pdf) by [E. Hairer](http://www.unige.ch/~hairer/). We solve the ODE system on the sphere ``\mathbb S^2`` given by @@ -27,7 +27,7 @@ We solve the ODE system on the sphere ``\mathbb S^2`` given by ``` ```julia -using ManifoldDiffEq, OrdinaryDiffEq, Manifolds +using ManifoldDiffEq, Manifolds using GLMakie, LinearAlgebra, Colors n = 25 @@ -72,15 +72,15 @@ S2 = Manifolds.Sphere(2) u0 = [0.0, sqrt(9/10), sqrt(1/10)] tspan = (0, 20.0) -A_lie = ManifoldDiffEq.LieManifoldDiffEqOperator{Float64}() do u, p, t +A_lie = LieManifoldDiffEqOperator{Float64}() do u, p, t return hat(SpecialOrthogonal(3), Matrix(I(3)), cross(u, f2(u...))) end -prob_lie = ManifoldDiffEq.ManifoldODEProblem(A_lie, u0, tspan, S2) +prob_lie = ManifoldODEProblem(A_lie, u0, tspan, S2) -A_frozen = ManifoldDiffEq.FrozenManifoldDiffEqOperator{Float64}() do u, p, t +A_frozen = FrozenManifoldDiffEqOperator{Float64}() do u, p, t return f2(u...) end -prob_frozen = ManifoldDiffEq.ManifoldODEProblem(A_frozen, u0, tspan, S2) +prob_frozen = ManifoldODEProblem(A_frozen, u0, tspan, S2) action = RotationAction(Euclidean(3), SpecialOrthogonal(3)) alg_lie_euler = ManifoldDiffEq.ManifoldLieEuler(S2, ExponentialRetraction(), action) diff --git a/docs/src/internals.md b/docs/src/internals.md new file mode 100644 index 0000000..ffd53af --- /dev/null +++ b/docs/src/internals.md @@ -0,0 +1,12 @@ +# Internal documentation + +This page documents the internal types and methods of `ManifoldDiffEq.jl`. + +## Functions + +```@docs +ManifoldDiffEq.AbstractManifoldDiffEqAlgorithm +ManifoldDiffEq.ManifoldInterpolationData +ManifoldDiffEq.AbstractManifoldDiffEqAdaptiveAlgorithm +ManifoldDiffEq.ManifoldODESolution +``` diff --git a/src/ManifoldDiffEq.jl b/src/ManifoldDiffEq.jl index d37cac9..ebad339 100644 --- a/src/ManifoldDiffEq.jl +++ b/src/ManifoldDiffEq.jl @@ -2,34 +2,683 @@ module ManifoldDiffEq using ManifoldsBase using Manifolds + +using LinearAlgebra + +import ConstructionBase: constructorof + +using Accessors: @set + +using SciMLBase: isinplace, promote_tspan, solve! + using SciMLBase: - SciMLBase, AbstractODEProblem, AbstractODEFunction, NullParameters, promote_tspan -import SciMLBase: build_solution -using OrdinaryDiffEq + AbstractDiffEqOperator, + AbstractODEFunction, + AbstractODEProblem, + CallbackSet, + NullParameters, + ODEFunction, + ReturnCode, + SciMLBase -using OrdinaryDiffEq: - InterpolationData, +import SciMLBase: alg_order, isadaptive, solution_new_retcode, solve + +using DiffEqBase: + DiffEqBase, + ODE_DEFAULT_NORM, + ODE_DEFAULT_ISOUTOFDOMAIN, + ODE_DEFAULT_UNSTABLE_CHECK, + ODE_DEFAULT_PROG_MESSAGE + +using OrdinaryDiffEqCore: + default_controller, + fsal_typeof, + gamma_default, + get_differential_vars, + handle_dt!, + initialize_callbacks!, + initialize_d_discontinuities, + initialize_saveat, + initialize_tstops, + isdtchangeable, + qmax_default, + qmin_default, + qsteady_max_default, + qsteady_min_default, + uses_uprev + +import OrdinaryDiffEqCore: perform_step! + +using OrdinaryDiffEqCore: + DEOptions, + ODEIntegrator, OrdinaryDiffEqAlgorithm, - OrdinaryDiffEqAdaptiveAlgorithm, - OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache, - trivial_limiter!, - constvalue, - @muladd, - @unpack, - @cache, - @.. - -import OrdinaryDiffEq: alg_cache, alg_order, initialize!, perform_step! + OrdinaryDiffEqCore, + OrdinaryDiffEqMutableCache using RecursiveArrayTools +using SimpleUnPack + +""" + abstract type AbstractManifoldDiffEqAlgorithm end + +A subtype of `OrdinaryDiffEqAlgorithm` for manifold-aware algorithms. +""" +abstract type AbstractManifoldDiffEqAlgorithm <: OrdinaryDiffEqAlgorithm end + +isadaptive(::AbstractManifoldDiffEqAlgorithm) = false + +""" + AbstractManifoldDiffEqAdaptiveAlgorithm <: AbstractManifoldDiffEqAlgorithm + +An abstract subtype of `AbstractManifoldDiffEqAlgorithm` for adaptive algorithms. +This is the manifold-aware analogue of `OrdinaryDiffEqAdaptiveAlgorithm`. + +""" +abstract type AbstractManifoldDiffEqAdaptiveAlgorithm <: AbstractManifoldDiffEqAlgorithm end + +isadaptive(::AbstractManifoldDiffEqAdaptiveAlgorithm) = true + +""" + struct ManifoldODESolution{T} end + +Counterpart of `SciMLBase.ODESolution`. It doesn't use the `N` parameter (because it +is not a generic manifold concept) and fields `u_analytic`, `errors`, `alg_choice`, +`original`, `tslocation` and `resid` (because we don't use them currently in +`ManifoldDiffEq.jl`). + +Type parameter `T` denotes scalar floating point type of the solution + +Fields: +* `u`: the representation of the ODE solution. Uses a nested power manifold representation. +* `t`: time points at which values in `u` were calculated. +* `k`: the representation of the `f` function evaluations at time points `k`. Uses a nested + power manifold representation. +* `prob`: original problem that was solved. +* `alg`: [`AbstractManifoldDiffEqAlgorithm`](@ref) used to obtain the solution. +* `interp` [`ManifoldInterpolationData`](@ref). It is used for calculating solution values + at times `t` other then the ones at which it was saved. +* `dense`: `true` if ODE solution is saved at every step and `false` otherwise. +* `stats`: [`DEStats`](https://docs.sciml.ai/DiffEqDocs/stable/basics/solution/#SciMLBase.DEStats) of the solver +* `retcode`: [`ReturnCode`](https://docs.sciml.ai/SciMLBase/stable/interfaces/Solutions/#retcodes) of the solution. +""" +struct ManifoldODESolution{ + T<:Real, + uType, + tType, + rateType, + P, + A<:AbstractManifoldDiffEqAlgorithm, + IType, + S, +} + u::uType + t::tType + k::rateType + prob::P + alg::A + interp::IType + dense::Bool + stats::S + retcode::ReturnCode.T +end + +function ManifoldODESolution{T}( + u, + t, + k, + prob, + alg, + interp, + dense, + stats, + retcode, +) where {T<:Real} + return ManifoldODESolution{ + T, + typeof(u), + typeof(t), + typeof(k), + typeof(prob), + typeof(alg), + typeof(interp), + typeof(stats), + }( + u, + t, + k, + prob, + alg, + interp, + dense, + stats, + retcode, + ) +end + +constructorof(::Type{<:ManifoldODESolution{T}}) where {T<:Real} = ManifoldODESolution{T} + +function solution_new_retcode(sol::ManifoldODESolution, retcode) + return @set sol.retcode = retcode +end + +function (sol::ManifoldODESolution)( + t, + ::Type{deriv} = Val{0}; + idxs = nothing, + continuity = :left, +) where {deriv} + return sol(t, deriv, idxs, continuity) +end + +function (sol::ManifoldODESolution)( + t::Number, + ::Type{deriv}, + idxs::Nothing, + continuity, +) where {deriv} + return sol.interp(t, idxs, deriv, sol.prob.p, continuity) +end + include("utils.jl") include("error_estimation.jl") include("operators.jl") -include("problems.jl") include("interpolation.jl") +include("problems.jl") + +alg_extrapolates(::AbstractManifoldDiffEqAlgorithm) = false + +struct DefaultInit end + + +# Adapted from OrdinaryDiffEq.jl: +# https://github.com/SciML/OrdinaryDiffEq.jl/blob/1eef9db17600766bb71e7dce0cb105ae5f99b2a5/lib/OrdinaryDiffEqCore/src/solve.jl#L11 +function SciMLBase.__init( + prob::ManifoldODEProblem, + alg::AbstractManifoldDiffEqAlgorithm, + timeseries_init = (), + ts_init = (), + ks_init = (), + recompile::Type{Val{recompile_flag}} = Val{true}; + saveat = (), + tstops = (), + d_discontinuities = (), + save_everystep = isempty(saveat), + save_on = true, + save_start = save_everystep || + isempty(saveat) || + saveat isa Number || + prob.tspan[1] in saveat, + save_end = nothing, + callback = nothing, + dense::Bool = save_everystep && isempty(saveat), + calck = (callback !== nothing && callback !== CallbackSet()) || + (dense) || + !isempty(saveat), # and no dense output + dt = eltype(prob.tspan)(0), + dtmin = eltype(prob.tspan)(0), + dtmax = eltype(prob.tspan)((prob.tspan[end] - prob.tspan[1])), + force_dtmin = false, + adaptive = isadaptive(alg), + gamma = gamma_default(alg), + abstol::Union{Nothing,Real} = nothing, + reltol::Union{Nothing,Real} = nothing, + qmin = qmin_default(alg), + qmax = qmax_default(alg), + qsteady_min = qsteady_min_default(alg), + qsteady_max = qsteady_max_default(alg), + qoldinit = isadaptive(alg) ? 1 // 10^4 : 0, + controller = nothing, + failfactor = 2, + maxiters::Int = isadaptive(alg) ? 1000000 : typemax(Int), + internalopnorm = opnorm, + isoutofdomain = ODE_DEFAULT_ISOUTOFDOMAIN, + unstable_check = ODE_DEFAULT_UNSTABLE_CHECK, + verbose = true, + timeseries_errors = true, + dense_errors = false, + advance_to_tstop = false, + stop_at_next_tstop = false, + initialize_save = true, + progress = false, + progress_steps = 1000, + progress_name = "ODE", + progress_message = ODE_DEFAULT_PROG_MESSAGE, + progress_id = :OrdinaryDiffEq, + userdata = nothing, + allow_extrapolation = alg_extrapolates(alg), + initialize_integrator = true, + initializealg = DefaultInit(), + kwargs..., +) where {recompile_flag} + isdae = false + + if !isempty(saveat) && dense + @warn( + "Dense output is incompatible with saveat. Please use the SavingCallback from the Callback Library to mix the two behaviors." + ) + end + + tType = eltype(prob.tspan) + tspan = prob.tspan + tdir = sign(tspan[end] - tspan[1]) + + t = tspan[1] + + _alg = alg + f = prob.f + p = prob.p + + # Get the control variables + + M = prob.manifold + + u = copy(M, prob.u0) + + du = nothing + duprev = nothing + + uType = typeof(u) + uBottomEltype = recursive_bottom_eltype(u) + uBottomEltypeNoUnits = recursive_unitless_bottom_eltype(u) + + uEltypeNoUnits = recursive_unitless_eltype(u) + tTypeNoUnits = typeof(one(tType)) + + if abstol === nothing + abstol_internal = real(convert(uBottomEltype, oneunit(uBottomEltype) * 1 // 10^6)) + else + abstol_internal = real(abstol) + end + + if reltol === nothing + reltol_internal = real(convert(uBottomEltype, oneunit(uBottomEltype) * 1 // 10^3)) + else + reltol_internal = real(reltol) + end + + dtmax > zero(dtmax) && tdir < 0 && (dtmax *= tdir) # Allow positive dtmax, but auto-convert + # dtmin is all abs => does not care about sign already. + + if isinplace(prob) && + u isa AbstractArray && + eltype(u) <: Number && + uBottomEltypeNoUnits == uBottomEltype && + tType == tTypeNoUnits # Could this be more efficient for other arrays? + rate_prototype = copy(M, u) + else + if (uBottomEltypeNoUnits == uBottomEltype && tType == tTypeNoUnits) || + eltype(u) <: Enum + rate_prototype = u + else # has units! + rate_prototype = u / oneunit(tType) + end + end + rateType = typeof(rate_prototype) ## Can be different if united + + tstops_internal = initialize_tstops(tType, tstops, d_discontinuities, tspan) + saveat_internal = initialize_saveat(tType, saveat, tspan) + d_discontinuities_internal = + initialize_d_discontinuities(tType, d_discontinuities, tspan) + + callbacks_internal = CallbackSet(callback) + + max_len_cb = DiffEqBase.max_vector_callback_length_int(callbacks_internal) + if max_len_cb !== nothing + uBottomEltypeReal = real(uBottomEltype) + if isinplace(prob) + callback_cache = DiffEqBase.CallbackCache( + u, + max_len_cb, + uBottomEltypeReal, + uBottomEltypeReal, + ) + else + callback_cache = + DiffEqBase.CallbackCache(max_len_cb, uBottomEltypeReal, uBottomEltypeReal) + end + else + callback_cache = nothing + end + + ### Algorithm-specific defaults ### + ksEltype = Vector{rateType} + + # Have to convert in case passed in wrong. + timeseries = timeseries_init === () ? uType[] : convert(Vector{uType}, timeseries_init) + + ts = ts_init === () ? tType[] : convert(Vector{tType}, ts_init) + ks = ks_init === () ? ksEltype[] : convert(Vector{ksEltype}, ks_init) + + if (!adaptive || !isadaptive(_alg)) && save_everystep && tspan[2] - tspan[1] != Inf + if dt == 0 + steps = length(tstops) + else + # For fixed dt, the only time dtmin makes sense is if it's smaller than eps(). + # Therefore user specified dtmin doesn't matter, but we need to ensure dt>=eps() + # to prevent infinite loops. + abs(dt) < dtmin && throw(ArgumentError("Supplied dt is smaller than dtmin")) + steps = ceil(Int, abs((tspan[2] - tspan[1]) / dt)) + end + sizehint!(timeseries, steps + 1) + sizehint!(ts, steps + 1) + sizehint!(ks, steps + 1) + elseif save_everystep + sizehint!(timeseries, 50) + sizehint!(ts, 50) + sizehint!(ks, 50) + elseif !isempty(saveat_internal) + savelength = length(saveat_internal) + 1 + if save_start == false + savelength -= 1 + end + if save_end == false && prob.tspan[2] in saveat_internal.valtree + savelength -= 1 + end + sizehint!(timeseries, savelength) + sizehint!(ts, savelength) + sizehint!(ks, savelength) + else + sizehint!(timeseries, 2) + sizehint!(ts, 2) + sizehint!(ks, 2) + end + + QT = number_eltype(u) + EEstT = number_eltype(u) + + k = rateType[] + + if uses_uprev(_alg, adaptive) || calck + uprev = copy(M, u) + else + # Some algorithms do not use `uprev` explicitly. In that case, we can save + # some memory by aliasing `uprev = u`, e.g. for "2N" low storage methods. + uprev = u + end + if allow_extrapolation + uprev2 = copy(M, u) + else + uprev2 = uprev + end + + cache = alg_cache( + _alg, + u, + rate_prototype, + uEltypeNoUnits, + uBottomEltypeNoUnits, + tTypeNoUnits, + uprev, + uprev2, + f, + t, + dt, + reltol_internal, + p, + calck, + Val(isinplace(prob)), + ) + + # Setting up the step size controller + if controller === nothing + controller = default_controller(_alg, cache, qoldinit, nothing, nothing) + end + + save_end_user = save_end + save_end = + save_end === nothing ? + save_everystep || isempty(saveat) || saveat isa Number || prob.tspan[2] in saveat : + save_end + + opts = DEOptions{ + typeof(abstol_internal), + typeof(reltol_internal), + QT, + tType, + typeof(controller), + typeof(ODE_DEFAULT_NORM), + typeof(internalopnorm), + typeof(save_end_user), + typeof(callbacks_internal), + typeof(isoutofdomain), + typeof(progress_message), + typeof(unstable_check), + typeof(tstops_internal), + typeof(d_discontinuities_internal), + typeof(userdata), + typeof(nothing), + typeof(maxiters), + typeof(tstops), + typeof(saveat), + typeof(d_discontinuities), + }( + maxiters, + save_everystep, + adaptive, + abstol_internal, + reltol_internal, + QT(gamma), + QT(qmax), + QT(qmin), + QT(qsteady_max), + QT(qsteady_min), + QT(qoldinit), + QT(failfactor), + tType(dtmax), + tType(dtmin), + controller, + ODE_DEFAULT_NORM, + internalopnorm, + nothing, # save_idxs + tstops_internal, + saveat_internal, + d_discontinuities_internal, + tstops, + saveat, + d_discontinuities, + userdata, + progress, + progress_steps, + progress_name, + progress_message, + progress_id, + timeseries_errors, + dense_errors, + dense, + save_on, + save_start, + save_end, + save_end_user, + callbacks_internal, + isoutofdomain, + unstable_check, + verbose, + calck, + force_dtmin, + advance_to_tstop, + stop_at_next_tstop, + ) + + stats = SciMLBase.DEStats(0) + differential_vars = get_differential_vars(f, u) + + + manifold_interp = ManifoldInterpolationData(f, timeseries, ts, ks, dense, cache, M) + sol = build_solution( + prob, + _alg, + ts, + timeseries, + dense = dense, + k = ks, + manifold_interp = manifold_interp, + stats = stats, + ) + + if recompile_flag == true + FType = typeof(f) + SolType = typeof(sol) + cacheType = typeof(cache) + else + FType = Function + SolType = DiffEqBase.AbstractDAESolution + cacheType = DAECache + end + + # rate/state = (state/time)/state = 1/t units, internalnorm drops units + # we don't want to differentiate through eigenvalue estimation + eigen_est = inv(one(tType)) + tprev = t + dtcache = tType(dt) + dtpropose = tType(dt) + iter = 0 + kshortsize = 0 + reeval_fsal = false + u_modified = false + EEst = EEstT(1) + just_hit_tstop = false + isout = false + accept_step = false + force_stepfail = false + last_stepfail = false + do_error_check = true + event_last_time = 0 + vector_event_last_time = 1 + last_event_error = zero(uBottomEltypeNoUnits) + dtchangeable = isdtchangeable(_alg) + q11 = QT(1) + success_iter = 0 + erracc = QT(1) + dtacc = tType(1) + reinitiailize = true + saveiter = 0 # Starts at 0 so first save is at 1 + saveiter_dense = 0 + #fsalfirst, fsallast = get_fsalfirstlast(cache, rate_prototype) + fsalfirst, fsallast = allocate(rate_prototype), allocate(rate_prototype) + + integrator = ODEIntegrator{ + typeof(_alg), + isinplace(prob), + uType, + typeof(du), + tType, + typeof(p), + typeof(eigen_est), + typeof(EEst), + QT, + typeof(tdir), + typeof(k), + SolType, + FType, + cacheType, + typeof(opts), + typeof(fsalfirst), + typeof(last_event_error), + typeof(callback_cache), + typeof(initializealg), + typeof(differential_vars), + }( + sol, + u, + du, + k, + t, + tType(dt), + f, + p, + uprev, + uprev2, + duprev, + tprev, + _alg, + dtcache, + dtchangeable, + dtpropose, + tdir, + eigen_est, + EEst, + QT(qoldinit), + q11, + erracc, + dtacc, + success_iter, + iter, + saveiter, + saveiter_dense, + cache, + callback_cache, + kshortsize, + force_stepfail, + last_stepfail, + just_hit_tstop, + do_error_check, + event_last_time, + vector_event_last_time, + last_event_error, + accept_step, + isout, + reeval_fsal, + u_modified, + reinitiailize, + isdae, + opts, + stats, + initializealg, + differential_vars, + fsalfirst, + fsallast, + ) + + if initialize_integrator + if SciMLBase.has_initializeprob(prob.f) + update_uprev!(integrator) + end + + if save_start + integrator.saveiter += 1 # Starts at 1 so first save is at 2 + integrator.saveiter_dense += 1 + copyat_or_push!(ts, 1, t) + # N.B.: integrator.u can be modified by initialized_dae! + copyat_or_push!(timeseries, 1, integrator.u) + copyat_or_push!(ks, 1, [rate_prototype]) + else + integrator.saveiter = 0 # Starts at 0 so first save is at 1 + integrator.saveiter_dense = 0 + end + + initialize_callbacks!(integrator, initialize_save) + initialize!(integrator, integrator.cache) + + end + + handle_dt!(integrator) + return integrator +end + +function solve( + prob::ManifoldODEProblem, + alg::AbstractManifoldDiffEqAlgorithm, + args...; + u0 = nothing, + p = nothing, + kwargs..., +) + u0 = u0 !== nothing ? u0 : prob.u0 + p = p !== nothing ? p : prob.p + + integrator = SciMLBase.__init(prob, alg, args...; kwargs...) + solve!(integrator) + return integrator.sol +end + include("frozen_solvers.jl") include("lie_solvers.jl") +export alg_order, solve + +export FrozenManifoldDiffEqOperator, LieManifoldDiffEqOperator, ManifoldODEProblem + end # module diff --git a/src/NEWS.md b/src/NEWS.md new file mode 100644 index 0000000..7a92565 --- /dev/null +++ b/src/NEWS.md @@ -0,0 +1,14 @@ +# Changelog + +All notable changes to this project will be documented in this file. + +The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), +and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). + +## [0.2.0] – 2024-09-03 + +### Changed + +* ODE solutions are no longer stored in `SciMLBase.ODESolution` but in `ManifoldODESolution` to avoid having to define `ndims` for all points. +* `SciMLBase.jl` newer than 2.39 is supported again. +* Julia 1.10 is now required due to SciML dependencies. diff --git a/src/frozen_solvers.jl b/src/frozen_solvers.jl index d0c05f8..8f95670 100644 --- a/src/frozen_solvers.jl +++ b/src/frozen_solvers.jl @@ -6,7 +6,7 @@ The manifold Euler algorithm for problems in the [`ExplicitManifoldODEProblemTyp formulation. """ struct ManifoldEuler{TM<:AbstractManifold,TR<:AbstractRetractionMethod} <: - OrdinaryDiffEqAlgorithm + AbstractManifoldDiffEqAlgorithm manifold::TM retraction_method::TR end @@ -16,7 +16,7 @@ alg_order(::ManifoldEuler) = 1 """ ManifoldEulerCache -Cache for [`ManifoldEuler`](@ref). +Mutable cache for [`ManifoldEuler`](@ref). """ struct ManifoldEulerCache <: OrdinaryDiffEqMutableCache end @@ -48,10 +48,12 @@ function alg_cache( end function perform_step!(integrator, ::ManifoldEulerCache, repeat_step = false) - @unpack t, dt, uprev, u, f, p, alg = integrator + u = integrator.u + t = integrator.t + alg = integrator.alg - k = f(u, p, t) - retract!(alg.manifold, u, u, k, dt, alg.retraction_method) + integrator.k[1] = integrator.f(u, integrator.p, t) + retract!(alg.manifold, u, u, integrator.k[1], integrator.dt, alg.retraction_method) return integrator.stats.nf += 1 end @@ -84,7 +86,8 @@ The Butcher tableau is identical to the Euclidean RK2: \end{array} ``` """ -struct CG2{TM<:AbstractManifold,TR<:AbstractRetractionMethod} <: OrdinaryDiffEqAlgorithm +struct CG2{TM<:AbstractManifold,TR<:AbstractRetractionMethod} <: + AbstractManifoldDiffEqAlgorithm manifold::TM retraction_method::TR end @@ -94,7 +97,7 @@ alg_order(::CG2) = 2 """ CG2Cache -Cache for [`CG2`](@ref). +Mutable cache for [`CG2`](@ref). """ struct CG2Cache{TX,TK2u} <: OrdinaryDiffEqMutableCache X1::TX @@ -171,7 +174,7 @@ The Butcher tableau reads (see tableau (5) of [EngøMarthinsen:1998](@cite)): The last row is used for error estimation. """ struct CG2_3{TM<:AbstractManifold,TR<:AbstractRetractionMethod} <: - OrdinaryDiffEqAdaptiveAlgorithm + AbstractManifoldDiffEqAdaptiveAlgorithm manifold::TM retraction_method::TR end @@ -219,6 +222,10 @@ function alg_cache( ) end +function DiffEqBase.get_tmp_cache(integrator, ::CG2_3, cache::CG2_3Cache) + return (cache.X1,) +end + function initialize!(integrator, cache::CG2_3Cache) integrator.kshortsize = 2 integrator.k = typeof(integrator.k)(undef, integrator.kshortsize) @@ -306,7 +313,8 @@ A Crouch-Grossmann algorithm of second order for problems in the ``` """ -struct CG3{TM<:AbstractManifold,TR<:AbstractRetractionMethod} <: OrdinaryDiffEqAlgorithm +struct CG3{TM<:AbstractManifold,TR<:AbstractRetractionMethod} <: + AbstractManifoldDiffEqAlgorithm manifold::TM retraction_method::TR end @@ -316,7 +324,7 @@ alg_order(::CG3) = 3 """ CG3Cache -Cache for [`CG3`](@ref). +Mutable cache for [`CG3`](@ref). """ struct CG3Cache{TX,TP} <: OrdinaryDiffEqMutableCache X1::TX @@ -404,7 +412,8 @@ A Crouch-Grossmann algorithm of second order for problems in the [`ExplicitManifoldODEProblemType`](@ref) formulation. See coefficients from Example 1 of [JackiewiczMarthinsenOwren:2000](@cite). """ -struct CG4a{TM<:AbstractManifold,TR<:AbstractRetractionMethod} <: OrdinaryDiffEqAlgorithm +struct CG4a{TM<:AbstractManifold,TR<:AbstractRetractionMethod} <: + AbstractManifoldDiffEqAlgorithm manifold::TM retraction_method::TR end @@ -414,7 +423,7 @@ alg_order(::CG4a) = 4 """ CG4aCache -Cache for [`CG4a`](@ref). +Mutable cache for [`CG4a`](@ref). """ struct CG4aCache{TX,TP} <: OrdinaryDiffEqMutableCache X1::TX diff --git a/src/interpolation.jl b/src/interpolation.jl index c7e5642..ec7ddad 100644 --- a/src/interpolation.jl +++ b/src/interpolation.jl @@ -1,7 +1,12 @@ -struct ManifoldInterpolationData{F,uType,tType,kType,cacheType,TM} <: - OrdinaryDiffEq.OrdinaryDiffEqInterpolation{cacheType} +""" + struct ManifoldInterpolationData end + +Inspired by [`OrdinaryDiffEq.InterpolationData`](https://github.com/SciML/OrdinaryDiffEq.jl/blob/41333beef24655d43d370af19b37efd9888179f6/lib/OrdinaryDiffEqCore/src/interp_func.jl#L4). +The main difference is using on-manifold interpolation instead of the Euclidean one. +""" +struct ManifoldInterpolationData{F,uType,tType,kType,cacheType,TM} f::F timeseries::uType ts::tType @@ -54,12 +59,13 @@ function ode_interpolation( if continuity === :left # we have i₋ = i₊ = 1 if tval = ts[1], i₊ = i₋ + 1 = lastindex(ts) if tval > ts[end], # and otherwise i₋ and i₊ satisfy ts[i₋] < tval ≤ ts[i₊] - i₊ = min(lastindex(ts), OrdinaryDiffEq._searchsortedfirst(ts, tval, 2, tdir > 0)) + i₊ = + min(lastindex(ts), OrdinaryDiffEqCore._searchsortedfirst(ts, tval, 2, tdir > 0)) i₋ = i₊ > 1 ? i₊ - 1 : i₊ else # we have i₋ = i₊ - 1 = 1 if tval < ts[1], i₊ = i₋ = lastindex(ts) if tval = ts[end], # and otherwise i₋ and i₊ satisfy ts[i₋] ≤ tval < ts[i₊] - i₋ = max(1, OrdinaryDiffEq._searchsortedlast(ts, tval, 1, tdir > 0)) + i₋ = max(1, OrdinaryDiffEqCore._searchsortedlast(ts, tval, 1, tdir > 0)) i₊ = i₋ < lastindex(ts) ? i₋ + 1 : i₋ end diff --git a/src/lie_solvers.jl b/src/lie_solvers.jl index 01e56ff..d0bb45b 100644 --- a/src/lie_solvers.jl +++ b/src/lie_solvers.jl @@ -9,7 +9,7 @@ struct ManifoldLieEuler{ TM<:AbstractManifold, TR<:AbstractRetractionMethod, TG<:AbstractGroupAction, -} <: OrdinaryDiffEqAlgorithm +} <: AbstractManifoldDiffEqAlgorithm manifold::TM retraction_method::TR action::TG @@ -20,7 +20,7 @@ alg_order(::ManifoldLieEuler) = 1 """ ManifoldLieEulerCache -Cache for [`ManifoldLieEuler`](@ref). +Mutable cache for [`ManifoldLieEuler`](@ref). """ struct ManifoldLieEulerCache{TID<:Identity} <: OrdinaryDiffEqMutableCache id::TID @@ -102,7 +102,7 @@ The Butcher tableau is: For more details see [MuntheKaasOwren:1999](@cite). """ struct RKMK4{TM<:AbstractManifold,TR<:AbstractRetractionMethod,TG<:AbstractGroupAction} <: - OrdinaryDiffEqAlgorithm + AbstractManifoldDiffEqAlgorithm manifold::TM retraction_method::TR action::TG @@ -113,7 +113,7 @@ alg_order(::RKMK4) = 4 """ RKMK4Cache -Cache for [`RKMK4`](@ref). +Mutable cache for [`RKMK4`](@ref). """ struct RKMK4Cache{TID<:Identity} <: OrdinaryDiffEqMutableCache id::TID diff --git a/src/operators.jl b/src/operators.jl index 68e531a..fefad13 100644 --- a/src/operators.jl +++ b/src/operators.jl @@ -32,12 +32,12 @@ function (vto::DefaultVectorTransportOperator)( end """ - FrozenManifoldDiffEqOperator{T<:Number,TM<:AbstractManifold,TF,TVT} <: SciMLBase.AbstractDiffEqOperator{T} + FrozenManifoldDiffEqOperator{T<:Number,TM<:AbstractManifold,TF,TVT} <: AbstractDiffEqOperator{T} DiffEq operator on manifolds in the frozen vector field formulation. """ struct FrozenManifoldDiffEqOperator{T<:Number,TF,TVT<:AbstractVectorTransportOperator} <: - SciMLBase.AbstractDiffEqOperator{T} + AbstractDiffEqOperator{T} func::TF operator_vector_transport::TVT end @@ -65,7 +65,7 @@ end DiffEq operator on manifolds in the Lie group action formulation. """ -struct LieManifoldDiffEqOperator{T<:Number,TF} <: SciMLBase.AbstractDiffEqOperator{T} +struct LieManifoldDiffEqOperator{T<:Number,TF} <: AbstractDiffEqOperator{T} func::TF end diff --git a/src/problems.jl b/src/problems.jl index f8e0205..b164625 100644 --- a/src/problems.jl +++ b/src/problems.jl @@ -97,7 +97,7 @@ struct ManifoldODEProblem{uType,tType,isinplace,P,F,K,PT,TM} <: end """ - ManifoldODEProblem{isinplace}(f,u0,tspan,p=NullParameters(),callback=CallbackSet()) + ManifoldODEProblem{isinplace}(f, u0, tspan, p=NullParameters(), callback=CallbackSet()) Define an ODE problem with the specified function. `isinplace` optionally sets whether the function is inplace or not. @@ -120,39 +120,6 @@ struct ManifoldODEProblem{uType,tType,isinplace,P,F,K,PT,TM} <: kwargs..., ) end - - function ManifoldODEProblem{iip,recompile}( - f, - u0, - tspan, - manifold::AbstractManifold, - p = NullParameters(); - kwargs..., - ) where {iip,recompile} - if !recompile - if iip - ManifoldODEProblem{iip}( - wrapfun_iip(f, (u0, u0, p, tspan[1])), - u0, - tspan, - manifold, - p; - kwargs..., - ) - else - ManifoldODEProblem{iip}( - wrapfun_oop(f, (u0, p, tspan[1])), - u0, - tspan, - manifold, - p; - kwargs..., - ) - end - else - ManifoldODEProblem{iip}(f, u0, tspan, p, manifold; kwargs...) - end - end end """ @@ -182,7 +149,7 @@ function ManifoldODEProblem( return ManifoldODEProblem(ODEFunction(f), u0, tspan, manifold, p; kwargs...) end -function OrdinaryDiffEq.ode_determine_initdt( +function ode_determine_initdt( u0, t, tdir, @@ -203,40 +170,22 @@ function build_solution( alg, t, u; - timeseries_errors = length(u) > 2, dense = false, - dense_errors = dense, - calculate_error = true, k = nothing, - interp::InterpolationData, + manifold_interp::ManifoldInterpolationData, retcode = ReturnCode.Default, stats = nothing, - kwargs..., ) T = eltype(eltype(u)) - - manifold_interp = ManifoldInterpolationData( - interp.f, - interp.timeseries, - interp.ts, - interp.ks, - interp.dense, - interp.cache, - prob.manifold, - ) - return ODESolution{T,1}( + return ManifoldODESolution{T}( u, - nothing, - nothing, t, k, prob, alg, manifold_interp, dense, - 0, stats, - nothing, retcode, ) end diff --git a/test/runtests.jl b/test/runtests.jl index 428b4f7..3f17296 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,10 +1,10 @@ using Test using ManifoldDiffEq using Manifolds -using OrdinaryDiffEq: OrdinaryDiffEq, alg_order using LinearAlgebra -using DiffEqBase using RecursiveArrayTools +using DiffEqBase +using OrdinaryDiffEq function test_solver_frozen(manifold_to_alg; expected_order = nothing, adaptive = false) expected_order !== nothing && @testset "alg_order" begin @@ -12,14 +12,14 @@ function test_solver_frozen(manifold_to_alg; expected_order = nothing, adaptive @test alg_order(alg) == expected_order end - @testset "Sphere" begin - M = Sphere(2) - A = ManifoldDiffEq.FrozenManifoldDiffEqOperator{Float64}() do u, p, t + M = Sphere(2) + alg = manifold_to_alg(M) + @testset "$alg on sphere" begin + A = FrozenManifoldDiffEqOperator{Float64}() do u, p, t return cross(u, [1.0, 0.0, 0.0]) end u0 = [0.0, 1.0, 0.0] - alg = manifold_to_alg(M) - prob = ManifoldDiffEq.ManifoldODEProblem(A, u0, (0, 2.0), M) + prob = ManifoldODEProblem(A, u0, (0, 2.0), M) sol1 = if adaptive solve(prob, alg) else @@ -30,14 +30,15 @@ function test_solver_frozen(manifold_to_alg; expected_order = nothing, adaptive @test is_point(M, sol1(1.0)) end - @testset "Product manifold" begin - M = ProductManifold(Sphere(2), Euclidean(3)) - A = ManifoldDiffEq.FrozenManifoldDiffEqOperator{Float64}() do u, p, t + M = ProductManifold(Sphere(2), Euclidean(3)) + alg = manifold_to_alg(M) + @testset "$alg on product manifold" begin + + A = FrozenManifoldDiffEqOperator{Float64}() do u, p, t return ArrayPartition(cross(u.x[1], [1.0, 0.0, 0.0]), u.x[2]) end u0 = ArrayPartition([0.0, 1.0, 0.0], [1.0, 0.0, 0.0]) - alg = manifold_to_alg(M) - prob = ManifoldDiffEq.ManifoldODEProblem(A, u0, (0, 2.0), M) + prob = ManifoldODEProblem(A, u0, (0, 2.0), M) sol1 = if adaptive solve(prob, alg) else @@ -61,12 +62,12 @@ function test_solver_lie(manifold_to_alg; expected_order = nothing) M = Sphere(2) action = RotationAction(Euclidean(3), SpecialOrthogonal(3)) - A = ManifoldDiffEq.LieManifoldDiffEqOperator{Float64}() do u, p, t + A = LieManifoldDiffEqOperator{Float64}() do u, p, t return hat(SpecialOrthogonal(3), Matrix(I(3)), cross(u, [1.0, 0.0, 0.0])) end u0 = [0.0, 1.0, 0.0] alg = manifold_to_alg(M, action) - prob = ManifoldDiffEq.ManifoldODEProblem(A, u0, (0, 2.0), M) + prob = ManifoldODEProblem(A, u0, (0, 2.0), M) sol1 = solve(prob, alg, dt = 1 / 8) @test sol1(0.0) ≈ u0 @@ -170,12 +171,12 @@ function compare_with_diffeq_frozen(manifold_to_alg, tableau) k = -0.25 c = 1 f(u, p, t) = [-c * u[1] - k * u[2], u[1]] - A = ManifoldDiffEq.FrozenManifoldDiffEqOperator{Float64}(f) + A = FrozenManifoldDiffEqOperator{Float64}(f) u0 = [-1.0, 1.0] alg = manifold_to_alg(M) tspan = (0, 2.0) dt = 1 / 8 - prob_frozen = ManifoldDiffEq.ManifoldODEProblem(A, u0, tspan, M) + prob_frozen = ManifoldODEProblem(A, u0, tspan, M) sol_frozen = solve(prob_frozen, alg, dt = dt) alg_diffeq = OrdinaryDiffEq.ExplicitRK(tableau) @@ -193,12 +194,12 @@ function compare_with_diffeq_lie(manifold_to_alg, tableau) f(u, p, t) = [-c * u[1] - k * u[2], u[1]] action = TranslationAction(Euclidean(2), TranslationGroup(2)) - A = ManifoldDiffEq.LieManifoldDiffEqOperator{Float64}(f) + A = LieManifoldDiffEqOperator{Float64}(f) u0 = [-1.0, 1.0] alg = manifold_to_alg(M, action) tspan = (0, 2.0) dt = 1 / 8 - prob_lie = ManifoldDiffEq.ManifoldODEProblem(A, u0, tspan, M) + prob_lie = ManifoldODEProblem(A, u0, tspan, M) sol_lie = solve(prob_lie, alg, dt = dt) alg_diffeq = OrdinaryDiffEq.ExplicitRK(tableau) @@ -240,4 +241,28 @@ end (M, action) -> ManifoldDiffEq.RKMK4(M, ExponentialRetraction(), action) test_solver_lie(manifold_to_alg_rkmk4; expected_order = 4) compare_with_diffeq_lie(manifold_to_alg_rkmk4, constructRKMK4()) + + @testset "abstol and reltol" begin + M = Sphere(2) + alg = ManifoldDiffEq.CG2_3(M, ExponentialRetraction()) + A = FrozenManifoldDiffEqOperator{Float64}() do u, p, t + return cross(u, [1.0, 0.0, 0.0]) + end + u0 = [0.0, 1.0, 0.0] + prob = ManifoldODEProblem(A, u0, (0, 2.0), M) + sol = solve(prob, alg; abstol = 0.1, reltol = 0.3) + end + + @testset "saveat" begin + M = Sphere(2) + alg = ManifoldDiffEq.CG2_3(M, ExponentialRetraction()) + A = FrozenManifoldDiffEqOperator{Float64}() do u, p, t + return cross(u, [1.0, 0.0, 0.0]) + end + u0 = [0.0, 1.0, 0.0] + prob = ManifoldODEProblem(A, u0, (0, 2.0), M) + saveat = [0.0, 0.2, 2.0] + sol = solve(prob, alg; saveat = saveat) + @test sol.t ≈ saveat + end end