juliagni/geometricintegrators.jl

Geometric Numerical Integration in Julia

GeometricIntegrators.jl

Julia library of geometric integrators for differential equations.

GeometricIntegrators.jl is a library of geometric integrators for ordinary differential equations, stochastic differential equations and differential algebraic equations in Julia. Its main aim is the democratization and proliferation of geometric integrators, providing a comprehensive collection of standard and structure-preserving algorithms under a unified interface. Furthermore it serves as testbed for the implementation and verification of novel geometric integrators, in particular their analysis with respect to long-time stability and conservation of geometric structures. GeometricIntegrators.jl can be used either interactively or as computational core in other codes. It provides both, a high-level interface that requires only very few lines of code to solve an actual problem, and a lean low-level interface that allows for straightforward integration into application codes via the exchange of very small data structures. Due to the modular structure and the use of the multiple dispatch paradigm, the library can easily be extended, e.g., towards new algorithms or new types of equations. GeometricIntegrators.jl is designed to minimize overhead and maximize performance in order to be able to perform simulations with millions or even billions of time steps to facilitate the study of the long-time behaviour of both numerical algorithms and dynamical systems.

Installation

GeometricIntegrators.jl and all of its dependencies can be installed via the Julia REPL by typing

]add GeometricIntegrators

Basic usage

In the simplest cases, the use of GeometricIntegrators.jl requires the construction of two objects, an equation and an integrator. For many standard methods, the integrator can be implicitly selected by specifying an equation and a tableau.

Before any use, we need to load GeometricIntegrators,

using GeometricIntegrators

Then we can create an ODE problem for the equation $\dot{x} (t) = x(t)$ with integration time span $(0, 1)$. a time step of $\Delta t = 0.1$, and initial condition $x(0) = 1$,

prob = ODEProblem((ẋ, t, x, params) -> ẋ[1] = x[1], (0.0, 1.0), 0.1, [1.0])

An integrator for this ODE, using the explicit Euler method is obtained by

int = GeometricIntegrator(prob, ExplicitEuler())

With that, the solution is simply computed by

sol = integrate(int)

which returns an object holding the solution for all time steps. With the help of the Plots.jl package we can visualise the result and compare with the exact solution:

using Plots
plot(xlims=[0,1], xlab="t", ylab="x(t)", legend=:bottomright)
plot!(sol.t, sol.q[:,1], label="numeric")
plot!(sol.t, exp.(sol.t), label="exact")

References

If you use GeometricIntegrators.jl in your work, please consider citing it by

@misc{Kraus:2020:GeometricIntegrators,
  title={GeometricIntegrators.jl: Geometric Numerical Integration in Julia},
  author={Kraus, Michael},
  year={2020},
  howpublished={\url{https://github.com/JuliaGNI/GeometricIntegrators.jl}},
  doi={10.5281/zenodo.3648325}
}

Development

We are using git hooks, e.g., to enforce that all tests pass before pushing. In order to activate these hooks, the following command must be executed once:

git config core.hooksPath .githooks

geometricintegrators.jl's Issues

Possibly generalize the indices?

I see lots of codes like:

for k in 1:sol.nd

Note that this will only loop though Array in Julia. If someone uses an OffsetArray or many other things in AbstractArray, this will not get the indices right, which is why it's normally nice to use eachindex(u) where u is the array if you're updating each element.

Do you plan on supporting non-traditional arrays? Note that this is something we'll be compiling a compatibility list of, so it's definitely okay to not go out of your way to do this. I think only OrdinaryDiffEq.jl is compatible with most of the odd types like this, while things like Sundials.jl, LSODA.jl, and ODEInterface.jl are on the other end and need Vector{Float64}.

Pull the Solvers Out to Another Package?

It seems you have a good set of linear and nonlinear solvers you're developing. What about pulling these out to separate packages, so that way one could use the linear solvers directly, or use the nonlinear solvers directly? I think those would be two very nice packages in the Julia ecosystem.

Add documentation for projection methods

Test GeometricIntegrators fails

[0234f1f7] HDF5_jll v1.12.2+2
[1d5cc7b8] IntelOpenMP_jll v2018.0.3+2
[1d63c593] LLVMOpenMP_jll v15.0.4+0
[856f044c] MKL_jll v2022.2.0+0
[458c3c95] OpenSSL_jll v1.1.19+0
[efe28fd5] OpenSpecFun_jll v0.5.5+0
[f50d1b31] Rmath_jll v0.3.0+0
[0dad84c5] ArgTools v1.1.1 @stdlib/ArgTools
[56f22d72] Artifacts @stdlib/Artifacts
[2a0f44e3] Base64 @stdlib/Base64
[ade2ca70] Dates @stdlib/Dates
[8bb1440f] DelimitedFiles @stdlib/DelimitedFiles
[8ba89e20] Distributed @stdlib/Distributed
[f43a241f] Downloads v1.6.0 @stdlib/Downloads
[7b1f6079] FileWatching @stdlib/FileWatching
[9fa8497b] Future @stdlib/Future
[b77e0a4c] InteractiveUtils @stdlib/InteractiveUtils
[4af54fe1] LazyArtifacts @stdlib/LazyArtifacts
[b27032c2] LibCURL v0.6.3 @stdlib/LibCURL
[76f85450] LibGit2 @stdlib/LibGit2
[8f399da3] Libdl @stdlib/Libdl
[37e2e46d] LinearAlgebra @stdlib/LinearAlgebra
[56ddb016] Logging @stdlib/Logging
[d6f4376e] Markdown @stdlib/Markdown
[a63ad114] Mmap @stdlib/Mmap
[ca575930] NetworkOptions v1.2.0 @stdlib/NetworkOptions
[44cfe95a] Pkg v1.8.0 @stdlib/Pkg
[de0858da] Printf @stdlib/Printf
[3fa0cd96] REPL @stdlib/REPL
[9a3f8284] Random @stdlib/Random
[ea8e919c] SHA v0.7.0 @stdlib/SHA
[9e88b42a] Serialization @stdlib/Serialization
[1a1011a3] SharedArrays @stdlib/SharedArrays
[6462fe0b] Sockets @stdlib/Sockets
[2f01184e] SparseArrays @stdlib/SparseArrays
[10745b16] Statistics @stdlib/Statistics
[4607b0f0] SuiteSparse @stdlib/SuiteSparse
[fa267f1f] TOML v1.0.0 @stdlib/TOML
[a4e569a6] Tar v1.10.1 @stdlib/Tar
[8dfed614] Test @stdlib/Test
[cf7118a7] UUIDs @stdlib/UUIDs
[4ec0a83e] Unicode @stdlib/Unicode
[e66e0078] CompilerSupportLibraries_jll v1.0.1+0 @stdlib/CompilerSupportLibraries_jll
[781609d7] GMP_jll v6.2.1+2 @stdlib/GMP_jll
[deac9b47] LibCURL_jll v7.84.0+0 @stdlib/LibCURL_jll
[29816b5a] LibSSH2_jll v1.10.2+0 @stdlib/LibSSH2_jll
[3a97d323] MPFR_jll v4.1.1+1 @stdlib/MPFR_jll
[c8ffd9c3] MbedTLS_jll v2.28.0+0 @stdlib/MbedTLS_jll
[14a3606d] MozillaCACerts_jll v2022.2.1 @stdlib/MozillaCACerts_jll
[4536629a] OpenBLAS_jll v0.3.20+0 @stdlib/OpenBLAS_jll
[05823500] OpenLibm_jll v0.8.1+0 @stdlib/OpenLibm_jll
[83775a58] Zlib_jll v1.2.12+3 @stdlib/Zlib_jll
[8e850b90] libblastrampoline_jll v5.1.1+0 @stdlib/libblastrampoline_jll
[8e850ede] nghttp2_jll v1.48.0+0 @stdlib/nghttp2_jll
[3f19e933] p7zip_jll v17.4.0+0 @stdlib/p7zip_jll
Info Packages marked with ⌅ have new versions available but compatibility constraints restrict them from upgrading.
Precompiling project...
✗ GeometricIntegrators
9 dependencies successfully precompiled in 70 seconds. 194 already precompiled.
1 dependency errored. To see a full report either run import Pkg; Pkg.precompile() or load the package
Testing Running tests...
ERROR: LoadError: TypeError: in Type{...} expression, expected UnionAll, got Type{SimpleSolvers.NonlinearSolver}
Stacktrace:
[1] top-level scope
@ C:\Users\Denis.julia\packages\GeometricIntegrators\mEZqa\src\integrators\rk\integrators_firk.jl:65
[2] include(mod::Module, _path::String)
@ Base .\Base.jl:419
[3] include(x::String)
@ GeometricIntegrators.Integrators C:\Users\Denis.julia\packages\GeometricIntegrators\mEZqa\src\Integrators.jl:1
[4] top-level scope
@ C:\Users\Denis.julia\packages\GeometricIntegrators\mEZqa\src\Integrators.jl:99
[5] include(mod::Module, _path::String)
@ Base .\Base.jl:419
[6] include(x::String)
@ GeometricIntegrators C:\Users\Denis.julia\packages\GeometricIntegrators\mEZqa\src\GeometricIntegrators.jl:1
[7] top-level scope
@ C:\Users\Denis.julia\packages\GeometricIntegrators\mEZqa\src\GeometricIntegrators.jl:15
[8] include
@ .\Base.jl:419 [inlined]
[9] include_package_for_output(pkg::Base.PkgId, input::String, depot_path::Vector{String}, dl_load_path::Vector{String}, load_path::Vector{String}, concrete_deps::Vector{Pair{Base.PkgId, UInt64}}, source::String)
@ Base .\loading.jl:1554
[10] top-level scope
@ stdin:1
in expression starting at C:\Users\Denis.julia\packages\GeometricIntegrators\mEZqa\src\integrators\rk\integrators_firk.jl:65
in expression starting at C:\Users\Denis.julia\packages\GeometricIntegrators\mEZqa\src\Integrators.jl:1
in expression starting at C:\Users\Denis.julia\packages\GeometricIntegrators\mEZqa\src\GeometricIntegrators.jl:1
in expression starting at stdin:1
Solution Step : Error During Test at C:\Users\Denis.julia\packages\SafeTestsets\A83XK\src\SafeTestsets.jl:25
Got exception outside of a @test
LoadError: Failed to precompile GeometricIntegrators [dcce2d33-59f6-5b8d-9047-0defad88ae06] to C:\Users\Denis.julia\compiled\v1.8\GeometricIntegrators\jl_3AD2.tmp.
Stacktrace:
[1] error(s::String)
@ Base .\error.jl:35
[2] compilecache(pkg::Base.PkgId, path::String, internal_stderr::IO, internal_stdout::IO, keep_loaded_modules::Bool)
@ Base .\loading.jl:1707
[3] compilecache
@ .\loading.jl:1651 [inlined]
[4] _require(pkg::Base.PkgId)
@ Base .\loading.jl:1337
[5] _require_prelocked(uuidkey::Base.PkgId)
@ Base .\loading.jl:1200
[6] macro expansion
@ .\loading.jl:1180 [inlined]
[7] macro expansion
@ .\lock.jl:223 [inlined]
[8] require(into::Module, mod::Symbol)
@ Base .\loading.jl:1144
[9] include(mod::Module, _path::String)
@ Base .\Base.jl:419
[10] include(x::String)
@ Main.var"##312".var"##313" C:\Users\Denis.julia\packages\SafeTestsets\A83XK\src\SafeTestsets.jl:23
[11] macro expansion
@ C:\Users\Denis.julia\packages\GeometricIntegrators\mEZqa\test\solutions\solutions_tests.jl:5 [inlined]
[12] macro expansion
@ C:\Users\Denis.julia\juliaup\julia-1.8.4+0.x64.w64.mingw32\share\julia\stdlib\v1.8\Test\src\Test.jl:1363 [inlined]
[13] top-level scope
@ C:\Users\Denis.julia\packages\GeometricIntegrators\mEZqa\test\solutions\solutions_tests.jl:5
[14] eval(m::Module, e::Any)
@ Core .\boot.jl:368
[15] top-level scope
@ C:\Users\Denis.julia\packages\SafeTestsets\A83XK\src\SafeTestsets.jl:23
[16] include(mod::Module, _path::String)
@ Base .\Base.jl:419
[17] include(x::String)
@ Main.var"##312" C:\Users\Denis.julia\packages\SafeTestsets\A83XK\src\SafeTestsets.jl:23
[18] macro expansion
@ C:\Users\Denis.julia\packages\GeometricIntegrators\mEZqa\test\runtests.jl:4 [inlined]
[19] macro expansion
@ C:\Users\Denis.julia\juliaup\julia-1.8.4+0.x64.w64.mingw32\share\julia\stdlib\v1.8\Test\src\Test.jl:1363 [inlined]
[20] top-level scope
@ C:\Users\Denis.julia\packages\GeometricIntegrators\mEZqa\test\runtests.jl:4
[21] eval(m::Module, e::Any)
@ Core .\boot.jl:368
[22] top-level scope
@ C:\Users\Denis.julia\packages\SafeTestsets\A83XK\src\SafeTestsets.jl:23
[23] include(fname::String)
@ Base.MainInclude .\client.jl:476
[24] top-level scope
@ none:6
[25] eval
@ .\boot.jl:368 [inlined]
[26] exec_options(opts::Base.JLOptions)
@ Base .\client.jl:276
[27] _start()
@ Base .\client.jl:522
in expression starting at C:\Users\Denis.julia\packages\GeometricIntegrators\mEZqa\test\solutions\solution_step_tests.jl:3
ERROR: LoadError: TypeError: in Type{...} expression, expected UnionAll, got Type{SimpleSolvers.NonlinearSolver}
Stacktrace:
[1] top-level scope
@ C:\Users\Denis.julia\packages\GeometricIntegrators\mEZqa\src\integrators\rk\integrators_firk.jl:65
[2] include(mod::Module, _path::String)
@ Base .\Base.jl:419
[3] include(x::String)
@ GeometricIntegrators.Integrators C:\Users\Denis.julia\packages\GeometricIntegrators\mEZqa\src\Integrators.jl:1
[4] top-level scope
@ C:\Users\Denis.julia\packages\GeometricIntegrators\mEZqa\src\Integrators.jl:99
[5] include(mod::Module, _path::String)
@ Base .\Base.jl:419
[6] include(x::String)
@ GeometricIntegrators C:\Users\Denis.julia\packages\GeometricIntegrators\mEZqa\src\GeometricIntegrators.jl:1
[7] top-level scope
@ C:\Users\Denis.julia\packages\GeometricIntegrators\mEZqa\src\GeometricIntegrators.jl:15
[8] include
@ .\Base.jl:419 [inlined]
[9] include_package_for_output(pkg::Base.PkgId, input::String, depot_path::Vector{String}, dl_load_path::Vector{String}, load_path::Vector{String}, concrete_deps::Vector{Pair{Base.PkgId, UInt64}}, source::String)
@ Base .\loading.jl:1554
[10] top-level scope
@ stdin:1
in expression starting at C:\Users\Denis.julia\packages\GeometricIntegrators\mEZqa\src\integrators\rk\integrators_firk.jl:65
in expression starting at C:\Users\Denis.julia\packages\GeometricIntegrators\mEZqa\src\Integrators.jl:1
in expression starting at C:\Users\Denis.julia\packages\GeometricIntegrators\mEZqa\src\GeometricIntegrators.jl:1
in expression starting at stdin:1
Deterministic Solutions : Error During Test at C:\Users\Denis.julia\packages\SafeTestsets\A83XK\src\SafeTestsets.jl:25
Got exception outside of a @test
LoadError: Failed to precompile GeometricIntegrators [dcce2d33-59f6-5b8d-9047-0defad88ae06] to C:\Users\Denis.julia\compiled\v1.8\GeometricIntegrators\jl_7BD4.tmp.
Stacktrace:
[1] error(s::String)
@ Base .\error.jl:35
[2] compilecache(pkg::Base.PkgId, path::String, internal_stderr::IO, internal_stdout::IO, keep_loaded_modules::Bool)
@ Base .\loading.jl:1707
[3] compilecache
@ .\loading.jl:1651 [inlined]
[4] _require(pkg::Base.PkgId)
@ Base .\loading.jl:1337
[5] _require_prelocked(uuidkey::Base.PkgId)
@ Base .\loading.jl:1200
[6] macro expansion
@ .\loading.jl:1180 [inlined]
[7] macro expansion
@ .\lock.jl:223 [inlined]
[8] require(into::Module, mod::Symbol)
@ Base .\loading.jl:1144
[9] include(mod::Module, _path::String)
@ Base .\Base.jl:419
[10] include(x::String)
@ Main.var"##312".var"##314" C:\Users\Denis.julia\packages\SafeTestsets\A83XK\src\SafeTestsets.jl:23
[11] macro expansion
@ C:\Users\Denis.julia\packages\GeometricIntegrators\mEZqa\test\solutions\solutions_tests.jl:8 [inlined]
[12] macro expansion
@ C:\Users\Denis.julia\juliaup\julia-1.8.4+0.x64.w64.mingw32\share\julia\stdlib\v1.8\Test\src\Test.jl:1363 [inlined]
[13] top-level scope
@ C:\Users\Denis.julia\packages\GeometricIntegrators\mEZqa\test\solutions\solutions_tests.jl:8
[14] eval(m::Module, e::Any)
@ Core .\boot.jl:368
[15] top-level scope
@ C:\Users\Denis.julia\packages\SafeTestsets\A83XK\src\SafeTestsets.jl:23
[16] include(mod::Module, _path::String)
@ Base .\Base.jl:419
[17] include(x::String)
@ Main.var"##312" C:\Users\Denis.julia\packages\SafeTestsets\A83XK\src\SafeTestsets.jl:23
[18] macro expansion
@ C:\Users\Denis.julia\packages\GeometricIntegrators\mEZqa\test\runtests.jl:4 [inlined]
[19] macro expansion
@ C:\Users\Denis.julia\juliaup\julia-1.8.4+0.x64.w64.mingw32\share\julia\stdlib\v1.8\Test\src\Test.jl:1363 [inlined]
[20] top-level scope
@ C:\Users\Denis.julia\packages\GeometricIntegrators\mEZqa\test\runtests.jl:4
[21] eval(m::Module, e::Any)
@ Core .\boot.jl:368
[22] top-level scope
@ C:\Users\Denis.julia\packages\SafeTestsets\A83XK\src\SafeTestsets.jl:23
[23] include(fname::String)
@ Base.MainInclude .\client.jl:476
[24] top-level scope
@ none:6
[25] eval
@ .\boot.jl:368 [inlined]
[26] exec_options(opts::Base.JLOptions)
@ Base .\client.jl:276
[27] _start()
@ Base .\client.jl:522
in expression starting at C:\Users\Denis.julia\packages\GeometricIntegrators\mEZqa\test\solutions\deterministic_solutions_tests.jl:2
Test Summary: | Error Total Time
Solution Tests | 2 2 30.3s
Solution Step | 1 1 15.3s
Deterministic Solutions | 1 1 15.0s
ERROR: LoadError: Some tests did not pass: 0 passed, 0 failed, 2 errored, 0 broken.
in expression starting at C:\Users\Denis.julia\packages\GeometricIntegrators\mEZqa\test\runtests.jl:4
ERROR: Package GeometricIntegrators errored during testing

(@v1.8) pkg>

Register v0.4.1

@JuliaRegistrator register

Example in the README.md fails

I get the following error using Julia v1.5.3.

pkg> st
Status `Project.toml`
  [dcce2d33] GeometricIntegrators v0.7.0
  [91a5bcdd] Plots v1.10.1

julia> using GeometricIntegrators

julia> ode = ODE((t, x, ẋ) -> ẋ[1] = x[1], [1.0]);

julia> int = Integrator(ode, TableauExplicitEuler(), 0.1);

julia> sol = integrate(ode, int, 10);

julia> using Plots

julia> plot(xlims=[0,1], xlab="t", ylab="x(t)", legend=:bottomright)

julia> plot!(sol.t, sol.q[1,:], label="numeric")
Error showing value of type Plots.Plot{Plots.GRBackend}:
ERROR: BoundsError: attempt to access 11-element OffsetArray(::Array{Float64,1}, 0:10) with eltype Float64 with indices 0:10 at index [0:11]
Stacktrace:
 [1] throw_boundserror(::OffsetArrays.OffsetArray{Float64,1,Array{Float64,1}}, ::Tuple{UnitRange{Int64}}) at ./abstractarray.jl:541
 [2] checkbounds at ./abstractarray.jl:506 [inlined]
 [3] _getindex at ./multidimensional.jl:742 [inlined]
 [4] getindex(::OffsetArrays.OffsetArray{Float64,1,Array{Float64,1}}, ::UnitRange{Int64}) at ./abstractarray.jl:1060
 [5] gr_draw_segments(::Plots.Series, ::OffsetArrays.OffsetArray{Float64,1,Array{Float64,1}}, ::Array{Float64,1}, ::Nothing, ::Tuple{Float64,Float64}) at ~/.julia/packages/Plots/lmp2A/src/backends/gr.jl:1645
 [6] gr_add_series(::Plots.Subplot{Plots.GRBackend}, ::Plots.Series) at ~/.julia/packages/Plots/lmp2A/src/backends/gr.jl:1560
 [7] gr_display(::Plots.Subplot{Plots.GRBackend}, ::Measures.Length{:mm,Float64}, ::Measures.Length{:mm,Float64}, ::Array{Float64,1}) at ~/.julia/packages/Plots/lmp2A/src/backends/gr.jl:935
 [8] gr_display(::Plots.Plot{Plots.GRBackend}, ::String) at ~/.julia/packages/Plots/lmp2A/src/backends/gr.jl:626
 [9] gr_display at ~/.julia/packages/Plots/lmp2A/src/backends/gr.jl:593 [inlined]
 [10] _display(::Plots.Plot{Plots.GRBackend}) at ~/.julia/packages/Plots/lmp2A/src/backends/gr.jl:1845
 [11] display(::Plots.PlotsDisplay, ::Plots.Plot{Plots.GRBackend}) at ~/.julia/packages/Plots/lmp2A/src/output.jl:150
 [12] display(::Any) at ./multimedia.jl:328
 [13] #invokelatest#1 at ./essentials.jl:710 [inlined]
 [14] invokelatest at ./essentials.jl:709 [inlined]
 [15] print_response(::IO, ::Any, ::Bool, ::Bool, ::Any) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.5/REPL/src/REPL.jl:238
 [16] print_response(::REPL.AbstractREPL, ::Any, ::Bool, ::Bool) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.5/REPL/src/REPL.jl:223
 [17] (::REPL.var"#do_respond#54"{Bool,Bool,REPL.var"#64#73"{REPL.LineEditREPL,REPL.REPLHistoryProvider},REPL.LineEditREPL,REPL.LineEdit.Prompt})(::Any, ::Any, ::Any) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.5/REPL/src/REPL.jl:822
 [18] (::REPL.var"#69#78")(::Any, ::Any, ::Vararg{Any,N} where N) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.5/REPL/src/REPL.jl:1086
 [19] #invokelatest#1 at ./essentials.jl:710 [inlined]
 [20] invokelatest at ./essentials.jl:709 [inlined]
 [21] (::REPL.LineEdit.var"#22#23"{REPL.var"#69#78",String})(::Any, ::Any) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.5/REPL/src/LineEdit.jl:1364
 [22] prompt!(::REPL.Terminals.TextTerminal, ::REPL.LineEdit.ModalInterface, ::REPL.LineEdit.MIState) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.5/REPL/src/LineEdit.jl:2447
 [23] run_interface(::REPL.Terminals.TextTerminal, ::REPL.LineEdit.ModalInterface, ::REPL.LineEdit.MIState) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.5/REPL/src/LineEdit.jl:2350
 [24] run_frontend(::REPL.LineEditREPL, ::REPL.REPLBackendRef) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.5/REPL/src/REPL.jl:1144
 [25] (::REPL.var"#38#42"{REPL.LineEditREPL,REPL.REPLBackendRef})() at ./task.jl:356

julia> plot!(sol.t, exp.(sol.t), label="exact")
Error showing value of type Plots.Plot{Plots.GRBackend}:
ERROR: BoundsError: attempt to access 11-element OffsetArray(::Array{Float64,1}, 0:10) with eltype Float64 with indices 0:10 at index [0:11]
Stacktrace:
 [1] throw_boundserror(::OffsetArrays.OffsetArray{Float64,1,Array{Float64,1}}, ::Tuple{UnitRange{Int64}}) at ./abstractarray.jl:541
 [2] checkbounds at ./abstractarray.jl:506 [inlined]
 [3] _getindex at ./multidimensional.jl:742 [inlined]
 [4] getindex(::OffsetArrays.OffsetArray{Float64,1,Array{Float64,1}}, ::UnitRange{Int64}) at ./abstractarray.jl:1060
 [5] gr_draw_segments(::Plots.Series, ::OffsetArrays.OffsetArray{Float64,1,Array{Float64,1}}, ::Array{Float64,1}, ::Nothing, ::Tuple{Float64,Float64}) at ~/.julia/packages/Plots/lmp2A/src/backends/gr.jl:1645
 [6] gr_add_series(::Plots.Subplot{Plots.GRBackend}, ::Plots.Series) at ~/.julia/packages/Plots/lmp2A/src/backends/gr.jl:1560
 [7] gr_display(::Plots.Subplot{Plots.GRBackend}, ::Measures.Length{:mm,Float64}, ::Measures.Length{:mm,Float64}, ::Array{Float64,1}) at ~/.julia/packages/Plots/lmp2A/src/backends/gr.jl:935
 [8] gr_display(::Plots.Plot{Plots.GRBackend}, ::String) at ~/.julia/packages/Plots/lmp2A/src/backends/gr.jl:626
 [9] gr_display at ~/.julia/packages/Plots/lmp2A/src/backends/gr.jl:593 [inlined]
 [10] _display(::Plots.Plot{Plots.GRBackend}) at ~/.julia/packages/Plots/lmp2A/src/backends/gr.jl:1845
 [11] display(::Plots.PlotsDisplay, ::Plots.Plot{Plots.GRBackend}) at ~/.julia/packages/Plots/lmp2A/src/output.jl:150
 [12] display(::Any) at ./multimedia.jl:328
 [13] #invokelatest#1 at ./essentials.jl:710 [inlined]
 [14] invokelatest at ./essentials.jl:709 [inlined]
 [15] print_response(::IO, ::Any, ::Bool, ::Bool, ::Any) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.5/REPL/src/REPL.jl:238
 [16] print_response(::REPL.AbstractREPL, ::Any, ::Bool, ::Bool) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.5/REPL/src/REPL.jl:223
 [17] (::REPL.var"#do_respond#54"{Bool,Bool,REPL.var"#64#73"{REPL.LineEditREPL,REPL.REPLHistoryProvider},REPL.LineEditREPL,REPL.LineEdit.Prompt})(::Any, ::Any, ::Any) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.5/REPL/src/REPL.jl:822
 [18] #invokelatest#1 at ./essentials.jl:710 [inlined]
 [19] invokelatest at ./essentials.jl:709 [inlined]
 [20] run_interface(::REPL.Terminals.TextTerminal, ::REPL.LineEdit.ModalInterface, ::REPL.LineEdit.MIState) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.5/REPL/src/LineEdit.jl:2355
 [21] run_frontend(::REPL.LineEditREPL, ::REPL.REPLBackendRef) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.5/REPL/src/REPL.jl:1144
 [22] (::REPL.var"#38#42"{REPL.LineEditREPL,REPL.REPLBackendRef})() at ./task.jl:356

Get references for algorithms

Currently, I do not know how to obtain a reference for a given algorithm that I should cite when using it in a research paper. For example, none of

julia> using GeometricIntegrators

julia> TableauQinZhang
TableauQinZhang (generic function with 1 method)

julia> TableauQinZhang()
TableauDIRK{Float64}(:qinzhang, 2, 2, Runge-Kutta Coefficients qinzhang with 2 stages and order 2  a = [0.25 0.0; 0.5 0.25]  b = [0.5, 0.5]  c = [0.25, 0.75])

help?> TableauQinZhang()
  Tableau of Qin and Zhang's symplectic two-stage, 2nd order method

help?> TableauQinZhang
search: TableauQinZhang

  Tableau of Qin and Zhang's symplectic two-stage, 2nd order method

works for this and I can't also find no reference in the source code.

Package Name

@musm brought up a good point in the Gitter chats that the package may need to be re-named if you want to eventually register it in the package ecosystem. I thought I might as well bring this up because you seem to be promoting it on your personal site, and it would be a bummer to get too deep and then have to rename it!

In case you're not aware, names which are acronyms and abbreviations are not allowed. In fact, this was said by the gatekeeper of the package management system @tkelman:

abbreviations are usually okay, but ODE and SDE are bad names that probably would be spelled out more if registered today. people not in the field will want to know what the package is for and possibly use it as well. longer names are clearer.

SciML/DifferentialEquations.jl#59 (comment)

Given that, if you'd want to register this package you'd very likely have to change the name so that "people not in the field will want to know what the package is for", i.e. don't abbreviate DAE.

That said, I think something like GeometricIntegrators.jl would work. Just something to think about.

Add documentation and references for degenerate variational integrators

Register v0.4.0.

@JuliaRegistrator register

Register v0.3.1

@JuliaRegistrator register

Improve docstrings

Some docstrings are not really helpful, e.g.

julia> using GeometricIntegrators
help?> Integrator
search: Integrator Integrators IntegratorRK IntegratorERK IntegratorIPRK IntegratorFLRK IntegratorFIRK

  Print error for integrators not implemented, yet.

  ─────────────────────────────────────────────────────────────────────────────────────────────────────────────

  Create integrator for explicit Runge-Kutta tableau.

  ─────────────────────────────────────────────────────────────────────────────────────────────────────────────

  Create integrator for diagonally implicit Runge-Kutta tableau.

  ─────────────────────────────────────────────────────────────────────────────────────────────────────────────

  Create integrator for fully implicit Runge-Kutta tableau.

  ─────────────────────────────────────────────────────────────────────────────────────────────────────────────

  Create integrator for explicit partitioned Runge-Kutta tableau.

  ─────────────────────────────────────────────────────────────────────────────────────────────────────────────

  Create integrator for implicit partitioned Runge-Kutta tableau.

  ─────────────────────────────────────────────────────────────────────────────────────────────────────────────

  Create integrator for variational partitioned Runge-Kutta tableau.

  ─────────────────────────────────────────────────────────────────────────────────────────────────────────────

  Create integrator for formal Lagrangian Runge-Kutta tableau.

  ─────────────────────────────────────────────────────────────────────────────────────────────────────────────

  Create integrator for Projected Gauss-Legendre Runge-Kutta tableau.

  ─────────────────────────────────────────────────────────────────────────────────────────────────────────────

  Create integrator for variational partitioned additive Runge-Kutta tableau.

  ─────────────────────────────────────────────────────────────────────────────────────────────────────────────

  Create integrator for special partitioned additive Runge-Kutta tableau.

  ─────────────────────────────────────────────────────────────────────────────────────────────────────────────

  Create integrator for variational special partitioned additive Runge-Kutta tableau.

  ─────────────────────────────────────────────────────────────────────────────────────────────────────────────

  Create integrator for variational special partitioned additive Runge-Kutta tableau with projection on primary
  constraint.

  ─────────────────────────────────────────────────────────────────────────────────────────────────────────────

  Create integrator for variational special partitioned additive Runge-Kutta tableau with projection on
  secondary constraint.

  ─────────────────────────────────────────────────────────────────────────────────────────────────────────────

  Create integrator for Hamiltonian partitioned additive Runge-Kutta tableau.

  ─────────────────────────────────────────────────────────────────────────────────────────────────────────────

  Create integrator for Hamiltonian special partitioned additive Runge-Kutta tableau.

  ─────────────────────────────────────────────────────────────────────────────────────────────────────────────

  Create integrator for Hamiltonian special partitioned additive Runge-Kutta tableau with projection on primary
  constraint.

  ─────────────────────────────────────────────────────────────────────────────────────────────────────────────

  Create integrator for splitting tableau.

  ─────────────────────────────────────────────────────────────────────────────────────────────────────────────

  Create integrator for stochastic explicit Runge-Kutta tableau.

  ─────────────────────────────────────────────────────────────────────────────────────────────────────────────

  Create integrator for weak explicit Runge-Kutta tableau.

  ─────────────────────────────────────────────────────────────────────────────────────────────────────────────

  Create integrator for stochastic fully implicit Runge-Kutta tableau.

  ─────────────────────────────────────────────────────────────────────────────────────────────────────────────

  Create integrator for stochastic fully implicit partitioned Runge-Kutta tableau.

  ─────────────────────────────────────────────────────────────────────────────────────────────────────────────

  Create integrator for stochastic fully implicit split partitioned Runge-Kutta tableau.

  ─────────────────────────────────────────────────────────────────────────────────────────────────────────────

  Create integrator for weak fully implicit Runge-Kutta tableau.

help?> integrate
search: integrate integrate! integrate_step! IntegratorERK IntegratorEPRK IntegratorExactODE

  Apply integrator for ntime time steps and return solution.

  ─────────────────────────────────────────────────────────────────────────────────────────────────────────────

  Integrate given equation with given tableau for ntime time steps and return solution.

  ─────────────────────────────────────────────────────────────────────────────────────────────────────────────

  Integrate ODE specified by vector field and initial condition with given tableau for ntime time steps and
  return solution.

  ─────────────────────────────────────────────────────────────────────────────────────────────────────────────

  Integrate PODE specified by two vector fields and initial conditions with given tableau for ntime time steps
  and return solution.

help?> TableauExplicitEuler
search: TableauExplicitEuler TableauExplicitMidpoint

  Tableau for explicit Euler method

I still don't know how to call the corresponding methods.

TestProblems not defined

Running the example from https://juliagni.github.io/GeometricIntegrators.jl/stable/integrators/

julia> using GeometricIntegrators

julia> ode = TestProblems.HarmonicOscillatorProblem.harmonic_oscillator_ode()
ERROR: UndefVarError: TestProblems not defined
Stacktrace:
 [1] top-level scope at REPL[2]:1

Thread-safety

Cache of implicit integrators is not thread-safe. e.g.

cache = NonlinearFunctionCacheSIPRK{ST}(D, M, S)

Possible solution: pre-generate cache Base.Threads.nthreads times and use according to threadid.

Register v0.2.0 to v0.3.0

@JuliaRegistrator register

Register v0.5.0

@JuliaRegistrator register

Add the JuliaDiffEq Common Interface

Hey,
I was wondering if you would like to include the code to plug into the JuliaDiffEq "common interface". This common interface is a type-based interface which makes all of the differential equations solvers have a common solver signature. The full discussion is here:

SciML/Roadmap#5

Since the discussion is long, let me give a quick summary. There is a problem type which defines the Cauchy problem to solve:

type ODEProblem{uType,tType,isinplace,F} <: AbstractODEProblem{uType,tType,isinplace,F}
  f::F
  u0::uType
  tspan::Vector{tType}
end

https://github.com/JuliaDiffEq/DiffEqBase.jl/blob/common_interface/src/problems.jl#L26

Each algorithm is a subtype of AbstractODEAlgorithm and thus the command solve(prob,alg;kwargs....) dispatches to the appropriate solvers from the appropriate packages. A similar setup exists for DAE problems (as exemplified by Sundials.jl and now DASKR.jl).

There are many advantages to adding this interface. For one, people already know and use this API, meaning that it will be easy for users to pick up your integrators as use them as a backend simply by changing the algorithm choice. Secondly, packages can write generic algorithms on the common interface, which then allows any backend solver to be used. For example, DiffEqParamEstim.jl performs parameter estimation using an ODE solver, and any ODE solver which implements the common interface can be used due to dispatch magic (you just pass in the algorithm type you wish to use). Another example is DiffEqSensitivity.jl which uses whichever backend to perform sensitivity analysis. Other tools include DiffEqUncertainty.jl which will be able to perform uncertainty quantification on numerical methods (once common callback interfaces are sorted out), and DiffEqDevTools which provides methods to test / benchmark any solver which implements the common interface (which is used to make these benchmarks).

You can see the full set of algorithms we already have collected here (this is missing DASKR.jl and LSODE.jl which have recently implemented this same API as well):

https://juliadiffeq.github.io/DiffEqDocs.jl/latest/solvers/ode_solve.html

It would be really nice to have geometric integrators! Let me know if you need starter code.

TagBot trigger issue

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

Partitoned Methods not working for HODEs

It seems like there is some typing issue with converting an HODEProblem to a PODEProblem such that partitioning methods don't working for Hamiltonian problems. Here's some code along with an error I get

function H(t, q, p, params)
    return 1/2*p[1]^2 - q[1]^2;
end

function f1(f1, t, q, p, params)
    f1[1] = -2*q[1]
end

function v(v, t, q, p, params);
    v[1] = p[1]
end

hprob = HODEProblem(v, f1, H, (0.0, 10.0), 0.1, (q = [0.0], p = [1.0]));
integrator = GeometricIntegrator(hprob, LobattoIIIAIIIB(2));
integrate(integrator);

and I get this error

MethodError: no method matching compute_stage_q!(::SolutionStepPODE{Float64, Float64, Vector{Float64}, Vector{Float64}, Vector{Float64}, NamedTuple{(:t, :q, :p, :v, :f), Tuple{OffsetArrays.OffsetVector{Float64, Vector{Float64}}, OffsetArrays.OffsetVector{Vector{Float64}, Vector{Vector{Float64}}}, OffsetArrays.OffsetVector{Vector{Float64}, Vector{Vector{Float64}}}, OffsetArrays.OffsetVector{Vector{Float64}, Vector{Vector{Float64}}}, OffsetArrays.OffsetVector{Vector{Float64}, Vector{Vector{Float64}}}}}, NamedTuple{(), Tuple{}}, NullParameters, 2}, ::HODEProblem{Float64, Float64, Vector{Float64}, HODE{typeof(v), typeof(f1), typeof(H), NullInvariants, NullParameters, NullPeriodicity}, NamedTuple{(:v, :f, :h), Tuple{typeof(v), typeof(f1), typeof(H)}}, NamedTuple{(), Tuple{}}, NamedTuple{(:q, :p), Tuple{Vector{Float64}, Vector{Float64}}}, NullParameters}, ::EPRK{PartitionedTableau{Float64, 2, Int64, Int64, Int64, 4}}, ::GeometricIntegrators.Integrators.EPRKCache{Float64, 1, 2}, ::Int64, ::Int64, ::Float64)
Stacktrace:
 [1] integrate_step!(int::GeometricIntegrator{EPRK{PartitionedTableau{Float64, 2, Int64, Int64, Int64, 4}}, HODEProblem{Float64, Float64, Vector{Float64}, HODE{typeof(v), typeof(f1), typeof(H), NullInvariants, NullParameters, NullPeriodicity}, NamedTuple{(:v, :f, :h), Tuple{typeof(v), typeof(f1), typeof(H)}}, NamedTuple{(), Tuple{}}, NamedTuple{(:q, :p), Tuple{Vector{Float64}, Vector{Float64}}}, NullParameters}, GeometricIntegrators.Integrators.CacheDict{HODEProblem{Float64, Float64, Vector{Float64}, HODE{typeof(v), typeof(f1), typeof(H), NullInvariants, NullParameters, NullPeriodicity}, NamedTuple{(:v, :f, :h), Tuple{typeof(v), typeof(f1), typeof(H)}}, NamedTuple{(), Tuple{}}, NamedTuple{(:q, :p), Tuple{Vector{Float64}, Vector{Float64}}}, NullParameters}, EPRK{PartitionedTableau{Float64, 2, Int64, Int64, Int64, 4}}}, NoSolver, NoInitialGuess, SolutionStepPODE{Float64, Float64, Vector{Float64}, Vector{Float64}, Vector{Float64}, NamedTuple{(:t, :q, :p, :v, :f), Tuple{OffsetArrays.OffsetVector{Float64, Vector{Float64}}, OffsetArrays.OffsetVector{Vector{Float64}, Vector{Vector{Float64}}}, OffsetArrays.OffsetVector{Vector{Float64}, Vector{Vector{Float64}}}, OffsetArrays.OffsetVector{Vector{Float64}, Vector{Vector{Float64}}}, OffsetArrays.OffsetVector{Vector{Float64}, Vector{Vector{Float64}}}}}, NamedTuple{(), Tuple{}}, NullParameters, 2}})
   @ GeometricIntegrators.Integrators C:\Users\ndefi\.julia\packages\GeometricIntegrators\uE6qJ\src\integrators\rk\integrators_eprk.jl:0
 [2] integrate!
   @ C:\Users\ndefi\.julia\packages\GeometricIntegrators\uE6qJ\src\integrators\integrator.jl:57 [inlined]
 [3] integrate!(sol::GeometricSolution{Float64, Float64, NamedTuple{(:q, :p), Tuple{DataSeries{Float64, Vector{Float64}}, DataSeries{Float64, Vector{Float64}}}}, HODEProblem{Float64, Float64, Vector{Float64}, HODE{typeof(v), typeof(f1), typeof(H), NullInvariants, NullParameters, NullPeriodicity}, NamedTuple{(:v, :f, :h), Tuple{typeof(v), typeof(f1), typeof(H)}}, NamedTuple{(), Tuple{}}, NamedTuple{(:q, :p), Tuple{Vector{Float64}, Vector{Float64}}}, NullParameters}, NamedTuple{(:q, :p), Tuple{NullPeriodicity, NullPeriodicity}}}, int::GeometricIntegrator{EPRK{PartitionedTableau{Float64, 2, Int64, Int64, Int64, 4}}, HODEProblem{Float64, Float64, Vector{Float64}, HODE{typeof(v), typeof(f1), typeof(H), NullInvariants, NullParameters, NullPeriodicity}, NamedTuple{(:v, :f, :h), Tuple{typeof(v), typeof(f1), typeof(H)}}, NamedTuple{(), Tuple{}}, NamedTuple{(:q, :p), Tuple{Vector{Float64}, Vector{Float64}}}, NullParameters}, GeometricIntegrators.Integrators.CacheDict{HODEProblem{Float64, Float64, Vector{Float64}, HODE{typeof(v), typeof(f1), typeof(H), NullInvariants, NullParameters, NullPeriodicity}, NamedTuple{(:v, :f, :h), Tuple{typeof(v), typeof(f1), typeof(H)}}, NamedTuple{(), Tuple{}}, NamedTuple{(:q, :p), Tuple{Vector{Float64}, Vector{Float64}}}, NullParameters}, EPRK{PartitionedTableau{Float64, 2, Int64, Int64, Int64, 4}}}, NoSolver, NoInitialGuess, SolutionStepPODE{Float64, Float64, Vector{Float64}, Vector{Float64}, Vector{Float64}, NamedTuple{(:t, :q, :p, :v, :f), Tuple{OffsetArrays.OffsetVector{Float64, Vector{Float64}}, OffsetArrays.OffsetVector{Vector{Float64}, Vector{Vector{Float64}}}, OffsetArrays.OffsetVector{Vector{Float64}, Vector{Vector{Float64}}}, OffsetArrays.OffsetVector{Vector{Float64}, Vector{Vector{Float64}}}, OffsetArrays.OffsetVector{Vector{Float64}, Vector{Vector{Float64}}}}}, NamedTuple{(), Tuple{}}, NullParameters, 2}}, n₁::Int64, n₂::Int64)
   @ GeometricIntegrators.Integrators C:\Users\ndefi\.julia\packages\GeometricIntegrators\uE6qJ\src\integrators\integrator.jl:89
 [4] integrate!(sol::GeometricSolution{Float64, Float64, NamedTuple{(:q, :p), Tuple{DataSeries{Float64, Vector{Float64}}, DataSeries{Float64, Vector{Float64}}}}, HODEProblem{Float64, Float64, Vector{Float64}, HODE{typeof(v), typeof(f1), typeof(H), NullInvariants, NullParameters, NullPeriodicity}, NamedTuple{(:v, :f, :h), Tuple{typeof(v), typeof(f1), typeof(H)}}, NamedTuple{(), Tuple{}}, NamedTuple{(:q, :p), Tuple{Vector{Float64}, Vector{Float64}}}, NullParameters}, NamedTuple{(:q, :p), Tuple{NullPeriodicity, NullPeriodicity}}}, int::GeometricIntegrator{EPRK{PartitionedTableau{Float64, 2, Int64, Int64, Int64, 4}}, HODEProblem{Float64, Float64, Vector{Float64}, HODE{typeof(v), typeof(f1), typeof(H), NullInvariants, NullParameters, NullPeriodicity}, NamedTuple{(:v, :f, :h), Tuple{typeof(v), typeof(f1), typeof(H)}}, NamedTuple{(), Tuple{}}, NamedTuple{(:q, :p), Tuple{Vector{Float64}, Vector{Float64}}}, NullParameters}, GeometricIntegrators.Integrators.CacheDict{HODEProblem{Float64, Float64, Vector{Float64}, HODE{typeof(v), typeof(f1), typeof(H), NullInvariants, NullParameters, NullPeriodicity}, NamedTuple{(:v, :f, :h), Tuple{typeof(v), typeof(f1), typeof(H)}}, NamedTuple{(), Tuple{}}, NamedTuple{(:q, :p), Tuple{Vector{Float64}, Vector{Float64}}}, NullParameters}, EPRK{PartitionedTableau{Float64, 2, Int64, Int64, Int64, 4}}}, NoSolver, NoInitialGuess, SolutionStepPODE{Float64, Float64, Vector{Float64}, Vector{Float64}, Vector{Float64}, NamedTuple{(:t, :q, :p, :v, :f), Tuple{OffsetArrays.OffsetVector{Float64, Vector{Float64}}, OffsetArrays.OffsetVector{Vector{Float64}, Vector{Vector{Float64}}}, OffsetArrays.OffsetVector{Vector{Float64}, Vector{Vector{Float64}}}, OffsetArrays.OffsetVector{Vector{Float64}, Vector{Vector{Float64}}}, OffsetArrays.OffsetVector{Vector{Float64}, Vector{Vector{Float64}}}}}, NamedTuple{(), Tuple{}}, NullParameters, 2}})
   @ GeometricIntegrators.Integrators C:\Users\ndefi\.julia\packages\GeometricIntegrators\uE6qJ\src\integrators\integrator.jl:124
 [5] integrate(integrator::GeometricIntegrator{EPRK{PartitionedTableau{Float64, 2, Int64, Int64, Int64, 4}}, HODEProblem{Float64, Float64, Vector{Float64}, HODE{typeof(v), typeof(f1), typeof(H), NullInvariants, NullParameters, NullPeriodicity}, NamedTuple{(:v, :f, :h), Tuple{typeof(v), typeof(f1), typeof(H)}}, NamedTuple{(), Tuple{}}, NamedTuple{(:q, :p), Tuple{Vector{Float64}, Vector{Float64}}}, NullParameters}, GeometricIntegrators.Integrators.CacheDict{HODEProblem{Float64, Float64, Vector{Float64}, HODE{typeof(v), typeof(f1), typeof(H), NullInvariants, NullParameters, NullPeriodicity}, NamedTuple{(:v, :f, :h), Tuple{typeof(v), typeof(f1), typeof(H)}}, NamedTuple{(), Tuple{}}, NamedTuple{(:q, :p), Tuple{Vector{Float64}, Vector{Float64}}}, NullParameters}, EPRK{PartitionedTableau{Float64, 2, Int64, Int64, Int64, 4}}}, NoSolver, NoInitialGuess, SolutionStepPODE{Float64, Float64, Vector{Float64}, Vector{Float64}, Vector{Float64}, NamedTuple{(:t, :q, :p, :v, :f), Tuple{OffsetArrays.OffsetVector{Float64, Vector{Float64}}, OffsetArrays.OffsetVector{Vector{Float64}, Vector{Vector{Float64}}}, OffsetArrays.OffsetVector{Vector{Float64}, Vector{Vector{Float64}}}, OffsetArrays.OffsetVector{Vector{Float64}, Vector{Vector{Float64}}}, OffsetArrays.OffsetVector{Vector{Float64}, Vector{Vector{Float64}}}}}, NamedTuple{(), Tuple{}}, NullParameters, 2}}; kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ GeometricIntegrators.Integrators C:\Users\ndefi\.julia\packages\GeometricIntegrators\uE6qJ\src\integrators\integrator.jl:25
 [6] integrate(integrator::GeometricIntegrator{EPRK{PartitionedTableau{Float64, 2, Int64, Int64, Int64, 4}}, HODEProblem{Float64, Float64, Vector{Float64}, HODE{typeof(v), typeof(f1), typeof(H), NullInvariants, NullParameters, NullPeriodicity}, NamedTuple{(:v, :f, :h), Tuple{typeof(v), typeof(f1), typeof(H)}}, NamedTuple{(), Tuple{}}, NamedTuple{(:q, :p), Tuple{Vector{Float64}, Vector{Float64}}}, NullParameters}, GeometricIntegrators.Integrators.CacheDict{HODEProblem{Float64, Float64, Vector{Float64}, HODE{typeof(v), typeof(f1), typeof(H), NullInvariants, NullParameters, NullPeriodicity}, NamedTuple{(:v, :f, :h), Tuple{typeof(v), typeof(f1), typeof(H)}}, NamedTuple{(), Tuple{}}, NamedTuple{(:q, :p), Tuple{Vector{Float64}, Vector{Float64}}}, NullParameters}, EPRK{PartitionedTableau{Float64, 2, Int64, Int64, Int64, 4}}}, NoSolver, NoInitialGuess, SolutionStepPODE{Float64, Float64, Vector{Float64}, Vector{Float64}, Vector{Float64}, NamedTuple{(:t, :q, :p, :v, :f), Tuple{OffsetArrays.OffsetVector{Float64, Vector{Float64}}, OffsetArrays.OffsetVector{Vector{Float64}, Vector{Vector{Float64}}}, OffsetArrays.OffsetVector{Vector{Float64}, Vector{Vector{Float64}}}, OffsetArrays.OffsetVector{Vector{Float64}, Vector{Vector{Float64}}}, OffsetArrays.OffsetVector{Vector{Float64}, Vector{Vector{Float64}}}}}, NamedTuple{(), Tuple{}}, NullParameters, 2}})
   @ GeometricIntegrators.Integrators C:\Users\ndefi\.julia\packages\GeometricIntegrators\uE6qJ\src\integrators\integrator.jl:24

Check timings on forced SIMD

Have you checked timings when using the forced SIMD functions like the ones here?

https://github.com/DDMGNI/GeomDAE.jl/blob/master/src/utils/matrix_utils.jl#L22

Note that in Julia v0.5, LLVM will automatically apply SIMD to a loop like that when it recognizes it would increase performance (with -O3). So with that I've only seen @simd slow down codes.

Register v0.5.1

@JuliaRegistrator register

Parameter `t` in the function `v` of ODEProblem is not the current time

Maybe I misunderstand something crucial, but I tried to integrate a time varying ODE and the result is just wrong. When I looked at the value t that I get, it seems to be always the timestep taken, not the total time of the current integration step.

Example:

using GeometricIntegrators

function eo(xdot, t, x, params)
    println(t)
    xdot[1] = x[1]
end

prob = ODEProblem(eo, (0.0, 1.0), 0.1, [3.0])
int = Integrator(prob, TableauRK4())

integrate(prob, int)

prints:

0.0
0.05
0.05
0.1
0.0
0.05
0.05
0.1
0.0
0.05
0.05
0.1
0.0
0.05
0.05
0.1
0.0
0.05
0.05
0.1
0.0
0.05
0.05
0.1
0.0
0.05
0.05
0.1
0.0
0.05
0.05
0.1
0.0
0.05
0.05
0.1
0.0
0.05
0.05
0.1

Register v0.3.2

@JuliaRegistrator register

Adaptive timestepping for partitioned symplectic methods

Hi,

I see in the docs that this repository implements explicit and implicit partitionned symplectic RK methods:
https://docs.juliahub.com/GeometricIntegrators/fhpp1/0.11.4/integrators/rk/

I am too creating a (much smaller) python code aiming at doing this, and I'm looking for some advice concerning adaptive timestepping.
Hence the questions: do you have adaptive timestepping capabilities that preserve symplecticity? If yes, could you please explain the concept or point me to a reference?

Thank you

Add developer documentation