From dd0e0e137dba308f33a1bf0f952dbc36a1b92073 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Mon, 17 Jan 2022 13:11:16 +0100 Subject: [PATCH] time integration methods and multirate infinitesimal split methods (#43) * add RootedTrees API reference * use RungeKuttaMethod, AdditiveRungeKuttaMethod * update CI * WIP: MultirateInfinitesimalSplitMethod * fix some type instabilities * another type stability fix * polynomial multiplication is fast * reexport Polynomial and add first tests * test for MIS methods Co-authored-by: OsKnoth * update CI to use parallel coveralls * fix typo * fix typo * increase test coverage * mode=PKGMODE_MANIFEST for Pkg.status in the docs * fix typo * fix typo * add TODO note Co-authored-by: OsKnoth --- .github/workflows/CI.yml | 75 ++++++- Project.toml | 4 +- docs/src/api_reference.md | 13 +- docs/src/benchmarks.md | 3 +- docs/src/index.md | 26 ++- src/BSeries.jl | 263 +++++++++++++++++++++--- test/runtests.jl | 415 +++++++++++++++++++++++++++++++++++++- 7 files changed, 761 insertions(+), 38 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 90141061..48f79fd3 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -11,11 +11,18 @@ on: jobs: test: runs-on: ubuntu-latest + strategy: + matrix: + group: + - Core + version: + - '1' + - '1.6' steps: - uses: actions/checkout@v2 - uses: julia-actions/setup-julia@v1 with: - version: 1 + version: ${{ matrix.version }} - uses: julia-actions/cache@v1 - uses: julia-actions/julia-buildpkg@v1 env: @@ -49,7 +56,73 @@ jobs: - uses: codecov/codecov-action@v1 with: file: lcov.info + # The standard setup of Coveralls is just annoying for parallel builds, see, e.g., + # https://github.com/trixi-framework/Trixi.jl/issues/691 + # https://github.com/coverallsapp/github-action/issues/47 + # https://github.com/coverallsapp/github-action/issues/67 + # This standard setup is reproduced below for completeness. + # - uses: coverallsapp/github-action@master + # with: + # github-token: ${{ secrets.GITHUB_TOKEN }} + # flag-name: run-${{ matrix.version }}-${{ github.run_id }} + # parallel: true + # path-to-lcov: ./lcov.info + # Instead, we use a more tedious approach: + # - Store all individual coverage files as artifacts (directly below) + # - Download and merge individual coverage reports in another step + # - Upload only the merged coverage report to Coveralls + - shell: bash + run: | + cp ./lcov.info ./lcov-${{ matrix.version }}-${{ github.run_id }}.info + - uses: actions/upload-artifact@v2 + with: + name: lcov-${{ matrix.version }}-${{ github.run_id }} + path: ./lcov-${{ matrix.version }}-${{ github.run_id }}.info + + finish: + needs: test + runs-on: ubuntu-latest + steps: + # The standard setup of Coveralls is just annoying for parallel builds, see, e.g., + # https://github.com/trixi-framework/Trixi.jl/issues/691 + # https://github.com/coverallsapp/github-action/issues/47 + # https://github.com/coverallsapp/github-action/issues/67 + # This standard setup is reproduced below for completeness. + # - name: Coveralls Finished + # uses: coverallsapp/github-action@master + # with: + # github-token: ${{ secrets.GITHUB_TOKEN }} + # parallel-finished: true + # Instead, we use the more tedious approach described above. + # At first, we check out the repository and download all artifacts + # (and list files for debugging). + - uses: actions/checkout@v2 + - uses: actions/download-artifact@v2 + - run: ls -R + # Next, we merge the individual coverage files and upload + # the combined results to Coveralls. + - name: Merge lcov files using Coverage.jl + shell: julia --color=yes {0} + run: | + using Pkg + Pkg.activate(temp=true) + Pkg.add("Coverage") + using Coverage + coverage = LCOV.readfolder(".") + for cov in coverage + cov.filename = replace(cov.filename, "\\" => "/") + end + coverage = merge_coverage_counts(coverage) + @show covered_lines, total_lines = get_summary(coverage) + LCOV.writefile("./lcov.info", coverage) - uses: coverallsapp/github-action@master with: github-token: ${{ secrets.GITHUB_TOKEN }} path-to-lcov: ./lcov.info + # Upload merged coverage data as artifact for debugging + - uses: actions/upload-artifact@v2 + with: + name: lcov + path: ./lcov.info + # That's it + - run: echo "Finished testing BSeries" diff --git a/Project.toml b/Project.toml index f31ef21c..fbe2f2b6 100644 --- a/Project.toml +++ b/Project.toml @@ -7,6 +7,7 @@ version = "0.1.17-pre" Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" +Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" Requires = "ae029012-a4dd-5104-9daa-d747884805df" RootedTrees = "47965b36-3f3e-11e9-0dcf-4570dfd42a8c" @@ -14,7 +15,8 @@ RootedTrees = "47965b36-3f3e-11e9-0dcf-4570dfd42a8c" [compat] Latexify = "0.15" OrderedCollections = "1" +Polynomials = "2.0.23" Reexport = "1" Requires = "1" -RootedTrees = "2.6.3" +RootedTrees = "2.8.3" julia = "1.6" diff --git a/docs/src/api_reference.md b/docs/src/api_reference.md index 3a80c1ec..5dd0c08d 100644 --- a/docs/src/api_reference.md +++ b/docs/src/api_reference.md @@ -1,4 +1,6 @@ -# BSeries.jl API +# API reference + +## BSeries.jl API ```@meta CurrentModule = BSeries @@ -8,3 +10,12 @@ CurrentModule = BSeries Modules = [BSeries] ``` +## RootedTrees.jl API + +```@meta +CurrentModule = RootedTrees +``` + +```@autodocs +Modules = [RootedTrees] +``` diff --git a/docs/src/benchmarks.md b/docs/src/benchmarks.md index f688bbb5..72230254 100644 --- a/docs/src/benchmarks.md +++ b/docs/src/benchmarks.md @@ -104,7 +104,8 @@ using InteractiveUtils versioninfo() using Pkg -Pkg.status(["BSeries", "RootedTrees", "SymEngine", "SymPy", "Symbolics"]) +Pkg.status(["BSeries", "RootedTrees", "SymEngine", "SymPy", "Symbolics"], + mode=PKGMODE_MANIFEST) nothing # hide ``` diff --git a/docs/src/index.md b/docs/src/index.md index d84d3861..29dffec5 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -1,10 +1,14 @@ # BSeries.jl -The Julia library [BSeries.jl](https://github.com/ranocha/BSeries.jl) -is work in progress. Nevertheless, we follow semantic versioning. Thus, you can -safely use the package right now. Extended documentation will be provided in the -future. +is a collection of functionality around +[B-series](https://en.wikipedia.org/wiki/Butcher_group) +in [Julia](https://julialang.org/). See for example + +- Philippe Chartier, Ernst Hairer, Gilles Vilmart (2010) + Algebraic Structures of B-series. + Foundations of Computational Mathematics + [DOI: 10.1007/s10208-010-9065-1](https://doi.org/10.1007/s10208-010-9065-1) BSeries.jl re-exports everything from [RootedTrees.jl](https://github.com/SciML/RootedTrees.jl). @@ -13,6 +17,20 @@ you should also include it explicitly in your project dependencies to track breaking changes, since the version numbers of RootedTrees.jl and BSeries.jl are not necessarily synchronized. +The main API of BSeries.jl consists of the following components. + +- B-series behave like `AbstractDict`s mapping + (abstract) `RootedTree`s to coefficients. +- The B-series of time integration methods such as Runge-Kutta methods + can be constructed by the function [`bseries`](@ref). +- The algebraic structures of the composition law and the substitution law are + implemented via [`compose`](@ref) and [`substitute`](@ref). +- Backward error analysis can be performed via + [`modified_equation`](@ref)s and [`modifying_integrator`](@ref)s. + +Further information is provided in the following tutorials and +API documentation. + ## Installation diff --git a/src/BSeries.jl b/src/BSeries.jl index dbb07a43..cfb47c89 100644 --- a/src/BSeries.jl +++ b/src/BSeries.jl @@ -6,7 +6,7 @@ module BSeries using Reexport: @reexport @reexport using RootedTrees -using RootedTrees: RootedTree +using RootedTrees: AbstractRootedTree @reexport using OrderedCollections: OrderedDict @@ -14,6 +14,8 @@ using Requires: @require using Latexify: Latexify, LaTeXString +@reexport using Polynomials: Polynomials, Polynomial + export TruncatedBSeries, ExactSolution @@ -23,6 +25,8 @@ export modified_equation, modifying_integrator export elementary_differentials +export MultirateInfinitesimalSplitMethod + # Types used for traits @@ -55,7 +59,7 @@ Generally, this kind of `struct` should be constructed via [`bseries`](@ref) or one of the other functions returning a B-series, e.g., [`modified_equation`](@ref) or [`modifying_integrator`](@ref). """ -struct TruncatedBSeries{T<:RootedTree, V} <: AbstractDict{T, V} +struct TruncatedBSeries{T<:AbstractRootedTree, V} <: AbstractDict{T, V} coef::OrderedDict{T, V} end @@ -72,18 +76,18 @@ TruncatedBSeries{T, V}() where {T, V} = TruncatedBSeries{T, V}(OrderedDict{T, V} @inline Base.length(series::TruncatedBSeries) = length(series.coef) -@inline Base.getindex(series::TruncatedBSeries, t::RootedTree) = getindex(series.coef, t) -@inline Base.setindex!(series::TruncatedBSeries, val, t::RootedTree) = setindex!(series.coef, val, t) +@inline Base.getindex(series::TruncatedBSeries, t::AbstractRootedTree) = getindex(series.coef, t) +@inline Base.setindex!(series::TruncatedBSeries, val, t::AbstractRootedTree) = setindex!(series.coef, val, t) -@inline Base.get(series::TruncatedBSeries, t::RootedTree, default) = get(series.coef, t, default) -@inline Base.get(f::F, series::TruncatedBSeries, t::RootedTree) where {F<:Function} = get(f, series.coef, t) +@inline Base.get(series::TruncatedBSeries, t::AbstractRootedTree, default) = get(series.coef, t, default) +@inline Base.get(f::F, series::TruncatedBSeries, t::AbstractRootedTree) where {F<:Function} = get(f, series.coef, t) -@inline Base.getkey(series::TruncatedBSeries, t::RootedTree, default) = getkey(series.coef, t, default) +@inline Base.getkey(series::TruncatedBSeries, t::AbstractRootedTree, default) = getkey(series.coef, t, default) -@inline Base.delete!(series::TruncatedBSeries, t::RootedTree) = (delete!(series.coef, t); series) +@inline Base.delete!(series::TruncatedBSeries, t::AbstractRootedTree) = (delete!(series.coef, t); series) -@inline Base.pop!(series::TruncatedBSeries, t::RootedTree) = pop!(series.coef, t) -@inline Base.pop!(series::TruncatedBSeries, t::RootedTree, default) = pop!(series.coef, t, default) +@inline Base.pop!(series::TruncatedBSeries, t::AbstractRootedTree) = pop!(series.coef, t) +@inline Base.pop!(series::TruncatedBSeries, t::AbstractRootedTree, default) = pop!(series.coef, t, default) @inline Base.sizehint!(series::TruncatedBSeries, n) = sizehint!(series.coef, n) @@ -118,7 +122,7 @@ differential equation using coefficients of type at least as representative as """ struct ExactSolution{V} end -Base.getindex(::ExactSolution{V}, t::RootedTree) where {V} = convert(V, 1//1) / γ(t) +Base.getindex(::ExactSolution{V}, t::AbstractRootedTree) where {V} = convert(V, 1//1) / γ(t) # general interface methods of iterators for `ExactSolution` Base.IteratorSize(::Type{<:ExactSolution}) = Base.SizeUnknown() @@ -187,9 +191,10 @@ end """ + bseries(rk::RungeKuttaMethod, order) bseries(A::AbstractMatrix, b::AbstractVector, c::AbstractVector, order) -Compute the B-series of the Runge-Kutta method with Butcher coefficients +Compute the B-series of the Runge-Kutta method `rk` with Butcher coefficients `A, b, c` up to a prescribed integer `order`. !!! note "Normalization by elementary differentials" @@ -199,9 +204,8 @@ Compute the B-series of the Runge-Kutta method with Butcher coefficients of the input vector field ``f``. See also [`evaluate`](@ref). """ -function bseries(A::AbstractMatrix, b::AbstractVector, c::AbstractVector, - order) - V_tmp = promote_type(eltype(A), eltype(b), eltype(c)) +function bseries(rk::RungeKuttaMethod, order) + V_tmp = eltype(rk) if V_tmp <: Integer # If people use integer coefficients, they will likely want to have results # as exact as possible. However, general terms are not integers. Thus, we @@ -215,14 +219,191 @@ function bseries(A::AbstractMatrix, b::AbstractVector, c::AbstractVector, series[rootedtree(Int[])] = one(V) for o in 1:order for t in RootedTreeIterator(o) - series[copy(t)] = elementary_weight(t, A, b, c) + series[copy(t)] = elementary_weight(t, rk) + end + end + + return series +end + +function bseries(A::AbstractMatrix, b::AbstractVector, c::AbstractVector, + order) + rk = RungeKuttaMethod(A, b, c) + bseries(rk, order) +end + +# TODO: bseries(rk::RungeKuttaMethod) +# should create a lazy version, optionally a memoized one + + +""" + bseries(ark::AdditiveRungeKuttaMethod, order) + +Compute the B-series of the additive Runge-Kutta method `ark` up to a prescribed +integer `order`. + +!!! note "Normalization by elementary differentials" + The coefficients of the B-series returned by this method need to be + multiplied by a power of the time step divided by the `symmetry` of the + colored rooted tree and multiplied by the corresponding elementary + differential of the input vector fields ``f^\\nu``. + See also [`evaluate`](@ref). +""" +function bseries(ark::AdditiveRungeKuttaMethod, order) + if length(ark.rks) != 2 + throw(ArgumentError("Only AdditiveRungeKuttaMethod with a dual splitting are supported. Got an ARK with $(length(ark.rks)) methods.")) + end + + V_tmp = eltype(ark) + if V_tmp <: Integer + # If people use integer coefficients, they will likely want to have results + # as exact as possible. However, general terms are not integers. Thus, we + # use rationals instead. + V = Rational{V_tmp} + else + V = V_tmp + end + series = TruncatedBSeries{BicoloredRootedTree{Int, Vector{Int}, Vector{Bool}}, V}() + + series[rootedtree(Int[], Bool[])] = one(V) + for o in 1:order + for t in BicoloredRootedTreeIterator(o) + series[copy(t)] = elementary_weight(t, ark) + end + end + + return series +end + +# TODO: bseries(ark::AdditiveRungeKuttaMethod) +# should create a lazy version, optionally a memoized one + + +# TODO: Documentation, Base.show, export etc. +""" + MultirateInfinitesimalSplitMethod(A, D, G, c) + +# References + +- Knoth, Oswald, and Joerg Wensch. + "Generalized split-explicit Runge-Kutta methods for the + compressible Euler equations". + Monthly Weather Review 142, no. 5 (2014): 2067-2081. + [DOI: 10.1175/MWR-D-13-00068.1](https://doi.org/10.1175/MWR-D-13-00068.1) + +!!! warning "Experimental code" + This code is considered to be experimental at the moment + and can change any time. +""" +struct MultirateInfinitesimalSplitMethod{T, PolyMatT<:AbstractMatrix{<:Polynomial{T}}, MatT<:AbstractMatrix{T}, VecT<:AbstractVector{T}} <: RootedTrees.AbstractTimeIntegrationMethod + A::PolyMatT + D::MatT + G::MatT + c::VecT +end + +# TODO: Deduce `c` from other parameters? +function MultirateInfinitesimalSplitMethod(A::AbstractMatrix{<:Polynomial}, + D::AbstractMatrix, + G::AbstractMatrix, + c::AbstractVector) + T = promote_type(eltype(eltype(A)), eltype(D), eltype(G), eltype(c)) + PolyT = typeof(zero(first(A)) + zero(T)) + _A = PolyT.(A) + _D = T.(D) + _G = T.(G) + _c = T.(c) + return MultirateInfinitesimalSplitMethod(_A, _D, _G, _c) +end + +Base.eltype(mis::MultirateInfinitesimalSplitMethod{T}) where {T} = T + +""" + bseries(mis::MultirateInfinitesimalSplitMethod, order) + +Compute the B-series of the multirate infinitesimal split method `mis` +up to a prescribed integer `order`. + +!!! note "Normalization by elementary differentials" + The coefficients of the B-series returned by this method need to be + multiplied by a power of the time step divided by the `symmetry` of the + colored rooted tree and multiplied by the corresponding elementary + differential of the input vector fields ``f^\\nu``. + See also [`evaluate`](@ref). +""" +function bseries(mis::MultirateInfinitesimalSplitMethod, order) + V_tmp = eltype(mis) + if V_tmp <: Integer + # If people use integer coefficients, they will likely want to have results + # as exact as possible. However, general terms are not integers. Thus, we + # use rationals instead. + V = Rational{V_tmp} + else + V = V_tmp + end + + prototype_scalar = TruncatedBSeries{BicoloredRootedTree{Int, Vector{Int}, Vector{Bool}}, V}() + prototype_scalar[rootedtree(Int[], Bool[])] = one(V) + + poly_one = Polynomial{V, :x}([one(V)]) + poly_zero = zero(poly_one) + prototype_polynomial = TruncatedBSeries{BicoloredRootedTree{Int, Vector{Int}, Vector{Bool}}, typeof(poly_one)}() + prototype_polynomial[rootedtree(Int[], Bool[])] = poly_one + + A = mis.A + D = mis.D + G = mis.G + ns = size(A, 2) + Z = Vector{typeof(prototype_polynomial)}(undef, ns + 1) + η = Vector{typeof(prototype_scalar)}(undef, ns + 1) + for i in 1:ns+1 + Z[i] = copy(prototype_polynomial) + η[i] = copy(prototype_scalar) + end + + for o in 1:order + for t_ in BicoloredRootedTreeIterator(o) + t = copy(t_) + for i in 1:ns+1 + phi = poly_zero + for j in 1:i-1 + r = poly_one + for subtree in RootedTrees.SubtreeIterator(t) + if subtree.color_sequence[1] + r = r * Z[i][subtree] + else + r = r * η[j][subtree] + end + end + + v = [one(V)] + for k in 0:length(A[i, j])-1 + phi = phi + Polynomial{V, :x}(v) * A[i, j][k] * r + v = [0; v] + end + end + + for j in 1:i-1 + phi = phi + G[i,j] * (η[j][t] - η[1][t]) + end + + phi = Polynomials.integrate(phi) + phi[0] = 0 + for j in 1:i-1 + phi[0] = phi[0] + D[i,j] * (η[j][t] - η[1][t]) + end + + Z[i][t] = phi + η[i][t] = phi(1) + end end end + series = η[end] return series end -# TODO: bseries(A::AbstractMatrix, b::AbstractVector, c::AbstractVector) +# TODO: bseries(mis::MultirateInfinitesimalSplitMethod) # should create a lazy version, optionally a memoized one @@ -503,11 +684,12 @@ function modified_equation(series_integrator, ::EagerEvaluation) end """ + modified_equation(rk::RungeKuttaMethod, order) modified_equation(A::AbstractMatrix, b::AbstractVector, c::AbstractVector, order) Compute the B-series of the [`modified_equation`](@ref) of the Runge-Kutta -method with Butcher coefficients `A, b, c` up to the prescribed `order`. +method `rk` with Butcher coefficients `A, b, c` up to the prescribed `order`. !!! note "Normalization by elementary differentials" The coefficients of the B-series returned by this method need to be @@ -516,13 +698,18 @@ method with Butcher coefficients `A, b, c` up to the prescribed `order`. of the input vector field ``f``. See also [`evaluate`](@ref). """ -function modified_equation(A::AbstractMatrix, b::AbstractVector, - c::AbstractVector, order) +function modified_equation(rk::RungeKuttaMethod, order) # B-series of the Runge-Kutta method - series = bseries(A, b, c, order) + series = bseries(rk, order) modified_equation(series) end +function modified_equation(A::AbstractMatrix, b::AbstractVector, + c::AbstractVector, order) + rk = RungeKuttaMethod(A, b, c) + modified_equation(rk, order) +end + """ modified_equation(f, u, dt, series_integrator) @@ -550,12 +737,13 @@ function modified_equation(f, u, dt, series_integrator) end """ + modified_equation(f, u, dt, rk::RungeKuttaMethod, order) modified_equation(f, u, dt, A::AbstractMatrix, b::AbstractVector, c::AbstractVector, order) Compute the B-series of the [`modified_equation`](@ref) of the Runge-Kutta -method with Butcher coefficients `A, b, c` up to the prescribed `order` with +method `rk` with Butcher coefficients `A, b, c` up to the prescribed `order` with respect to the ordinary differential equation ``u'(t) = f(u(t))`` with vector field `f` and dependent variables `u` for a time step size `dt`. @@ -569,11 +757,16 @@ from are supported. """ +function modified_equation(f, u, dt, rk::RungeKuttaMethod, order) + series_integrator = bseries(rk, order) + modified_equation(f, u, dt, series_integrator) +end + function modified_equation(f, u, dt, A::AbstractMatrix, b::AbstractVector, c::AbstractVector, order) - series_integrator = bseries(A, b, c, order) - modified_equation(f, u, dt, series_integrator) + rk = RungeKuttaMethod(A, b, c) + modified_equation(f, u, dt, rk, order) end @@ -657,6 +850,7 @@ function modifying_integrator(series_integrator, ::EagerEvaluation) end """ + modifying_integrator(rk::RungeKuttaMethod, order) modifying_integrator(A::AbstractMatrix, b::AbstractVector, c::AbstractVector, order) @@ -671,13 +865,18 @@ Runge-Kutta method with Butcher coefficients `A, b, c` up to the prescribed of the input vector field ``f``. See also [`evaluate`](@ref). """ -function modifying_integrator(A::AbstractMatrix, b::AbstractVector, - c::AbstractVector, order) +function modifying_integrator(rk::RungeKuttaMethod, order) # B-series of the Runge-Kutta method - series = bseries(A, b, c, order) + series = bseries(rk, order) modifying_integrator(series) end +function modifying_integrator(A::AbstractMatrix, b::AbstractVector, + c::AbstractVector, order) + rk = RungeKuttaMethod(A, b, c) + modifying_integrator(rk, order) +end + """ modifying_integrator(f, u, dt, series_integrator) @@ -705,6 +904,7 @@ function modifying_integrator(f, u, dt, series_integrator) end """ + modifying_integrator(f, u, dt, rk::RungeKuttaMethod, order) modifying_integrator(f, u, dt, A::AbstractMatrix, b::AbstractVector, c::AbstractVector, order) @@ -724,11 +924,16 @@ from are supported. """ +function modifying_integrator(f, u, dt, rk::RungeKuttaMethod, order) + series_integrator = bseries(rk, order) + modifying_integrator(f, u, dt, series_integrator) +end + function modifying_integrator(f, u, dt, A::AbstractMatrix, b::AbstractVector, c::AbstractVector, order) - series_integrator = bseries(A, b, c, order) - modifying_integrator(f, u, dt, series_integrator) + rk = RungeKuttaMethod(A, b, c) + modifying_integrator(f, u, dt, rk, order) end diff --git a/test/runtests.jl b/test/runtests.jl index be133c36..5c44aabd 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,7 +3,7 @@ using BSeries using BSeries.Latexify: latexify -using StaticArrays: @SArray +using StaticArrays: @SArray, @SMatrix, @SVector using SymEngine: SymEngine using SymPy: SymPy @@ -732,4 +732,417 @@ end end +@testset "Runge-Kutta methods interface" begin + # classical RK4 + A = [0 0 0 0 + 1//2 0 0 0 + 0 1//2 0 0 + 0 0 1 0] + b = [1//6, 1//3, 1//3, 1//6] + rk = RungeKuttaMethod(A, b) + + # fourth-order accurate + series_integrator = @inferred bseries(rk, 4) + series_exact = @inferred ExactSolution(series_integrator) + @test mapreduce(==, &, values(series_integrator), values(series_exact)) + + # not fifth-order accurate + series_integrator = @inferred bseries(rk, 5) + series_exact = @inferred ExactSolution(series_integrator) + @test mapreduce(==, &, values(series_integrator), values(series_exact)) == false +end + + +@testset "additive Runge-Kutta methods interface" begin + @testset "IMEX Euler" begin + ex_euler = RungeKuttaMethod( + @SMatrix([0]), @SVector [1] + ) + im_euler = RungeKuttaMethod( + @SMatrix([1]), @SVector [1] + ) + ark = AdditiveRungeKuttaMethod([ex_euler, im_euler]) + + # first-order accurate + series_integrator = @inferred bseries(ark, 1) + series_exact = @inferred ExactSolution(series_integrator) + @test mapreduce(==, &, values(series_integrator), values(series_exact)) + + # not second-order accurate + series_integrator = @inferred bseries(ark, 2) + series_exact = @inferred ExactSolution(series_integrator) + @test mapreduce(==, &, values(series_integrator), values(series_exact)) == false + end + + @testset "Störmer-Verlet" begin + # Hairer, Lubich, Wanner (2002) + # Geometric numerical integration + # Table II.2.1 + As = [ + [0 0; 1//2 1//2], + [1//2 0; 1//2 0] + ] + bs = [ + [1//2, 1//2], + [1//2, 1//2] + ] + ark = AdditiveRungeKuttaMethod(As, bs) + + # second-order accurate + series_integrator = @inferred bseries(ark, 2) + series_exact = @inferred ExactSolution(series_integrator) + @test mapreduce(==, &, values(series_integrator), values(series_exact)) + + # not third-order accurate + series_integrator = @inferred bseries(ark, 3) + series_exact = @inferred ExactSolution(series_integrator) + @test mapreduce(==, &, values(series_integrator), values(series_exact)) == false + end +end + + +@testset "multirate infinitesimal split methods interface" begin + # NOTE: This is still considered experimental and might change at any time! + + @testset "KW43" begin + # Oswald Knoth, Ralf Wolke + # Implicit-explicit Runge-Kutta methods for computing atmospheric reactive flows + # Applied Numerical Mathematics, Volume 28, (1998) Issue 2-4 + # https://doi.org/10.1016/S0168-9274(98)00051-8 + ns = 4 + A = zeros(Polynomial{Rational{Int}, :x}, ns+1, ns) + D = zeros(Int, ns+1, ns) + G = zeros(Rational{Int}, ns+1, ns) + c = zeros(Rational{Int}, ns+1) + dts = zeros(Rational{Int}, ns+1) + + A[2, 1] = 1 // 2 + A[3, 1] = -1 // 6 + A[3, 2] = 2 // 3 + A[4, 1] = 1 // 3 + A[4, 2] = -1 // 3 + A[4, 3] = 1 + A[5, 1] = 1 // 6 + A[5, 2] = 1 // 3 + A[5, 3] = 1 // 3 + A[5, 4] = 1 // 6 + A[5, :] = A[5, :] - A[4, :] + A[4, :] = A[4, :] - A[3, :] + A[3, :] = A[3, :] - A[2, :] + + c[2] = 1 // 2 + c[3] = 1 // 2 + c[4] = 1 + c[5] = 1 + + D[2,1] = 1 + D[3,2] = 1 + D[4,3] = 1 + D[5,4] = 1 + + dts[1] = 0 + dts[2] = 1 // 2 + dts[3] = 0 + dts[4] = 1 // 2 + dts[5] = 0 + + mis = @inferred MultirateInfinitesimalSplitMethod(A, D, G, c) + + # third-order accurate + series_integrator = @inferred bseries(mis, 3) + series_exact = @inferred ExactSolution(series_integrator) + @test mapreduce(==, &, values(series_integrator), values(series_exact)) + + # not fourth-order accurate + series_integrator = @inferred bseries(mis, 4) + series_exact = @inferred ExactSolution(series_integrator) + @test mapreduce(==, &, values(series_integrator), values(series_exact)) == false + end + + @testset "MRI-GARK-ERK33" begin + # Sandu, Adrian. + # "A class of multirate infinitesimal GARK methods". + # SIAM Journal on Numerical Analysis 57, no. 5 (2019): 2300-2327. + # https://doi.org/10.1137/18M1205492 + ns = 3 + A = zeros(Polynomial{Float64, :x}, ns+1, ns) + D = zeros(ns+1, ns) + G = zeros(ns+1, ns) + c = zeros(ns+1) + dts = zeros(ns+1) + + delta = -0.5 + + # Γ^i in Sandu (2019) is the i-th term + A[2, 1] = Polynomial([1 / 3]) + A[3, 1] = Polynomial([(-6 * delta - 7) / 12, (2 * delta + 1) / 2]) + A[3, 2] = Polynomial([(6 * delta + 11) / 12, -(2 * delta + 1) / 2]) + A[4, 1] = Polynomial([0.0, 0.5]) + A[4, 2] = Polynomial([(6 * delta - 5) / 12, -(2 * delta + 1) / 2]) + A[4, 3] = Polynomial([(3 - 2 * delta) / 4, delta]) + + # These should be the usual `c` coefficients → c[3] = 2 / 3? + c[2] = 1 / 3 + c[3] = 2 / 3 + c[4] = 1 + + # All one??? + D[2, 1] = 1 + D[3, 2] = 1 + D[4, 3] = 1 + + # Differences of `c`s? + dts[1] = 0 + dts[2] = 1 / 3 + dts[3] = 1 / 3 + dts[4] = 1 / 3 + + mis = @inferred MultirateInfinitesimalSplitMethod(A, D, G, c) + + # third-order accurate + series_integrator = @inferred bseries(mis, 3) + series_exact = @inferred ExactSolution(series_integrator) + @test mapreduce(≈, &, values(series_integrator), values(series_exact)) + + # not fourth-order accurate + series_integrator = @inferred bseries(mis, 4) + series_exact = @inferred ExactSolution(series_integrator) + @test mapreduce(≈, &, values(series_integrator), values(series_exact)) == false + end + + @testset "MRI-GARK-ERK45a" begin + # Sandu, Adrian. + # "A class of multirate infinitesimal GARK methods". + # SIAM Journal on Numerical Analysis 57, no. 5 (2019): 2300-2327. + # https://doi.org/10.1137/18M1205492 + ns = 5 + A = zeros(Polynomial{Float64, :x}, ns+1, ns) + D = zeros(ns+1, ns) + G = zeros(ns+1, ns) + c = zeros(ns+1) + dts = zeros(ns+1) + + # Γ^i in Sandu (2019) is the i-th term + A[2,1] = Polynomial([1 / 5]) + + A[3,1] = Polynomial([-53 / 16, 503 / 80]) + A[3,2] = Polynomial([281 / 80, -503 / 80]) + + A[4,1] = Polynomial([-36562993.0/71394880.0, -1365537.0/35697440.0]) + A[4,2] = Polynomial([34903117.0/17848720.0, 4963773.0/7139488.0]) + A[4,3] = Polynomial([-88770499.0/71394880.0, -1465833.0/2231090.0]) + + A[5,1] = Polynomial([-7631593.0/71394880.0, 66974357.0/35697440.0]) + A[5,2] = Polynomial([-166232021.0/35697440.0, 21445367.0/7139488.0]) + A[5,3] = Polynomial([6068517.0/1519040.0, -3.0]) + A[5,4] = Polynomial([8644289.0/8924360.0, -8388609.0/4462180.0]) + + A[6,1] = Polynomial([277061.0/303808.0, -18227.0/7520.0]) + A[6,2] = Polynomial([-209323.0/1139280.0, 2.0]) + A[6,3] = Polynomial([-1360217.0/1139280.0, 1.0]) + A[6,4] = Polynomial([-148789.0/56964.0, 5.0]) + A[6,5] = Polynomial([147889.0/45120.0, -41933.0/7520.0]) + + # These should be the usual `c` coefficients + c[2] = 1 / 5 + c[3] = 2 / 5 + c[4] = 3 / 5 + c[5] = 4 / 5 + c[6] = 1 + + D[2,1] = 0 + D[3,2] = 1 + D[4,3] = 1 + D[5,4] = 1 + D[6,5] = 1 + + # Differences of `c`s? + dts[1] = 0.0 + dts[2] = 1 / 5 + dts[3] = 1 / 5 + dts[4] = 1 / 5 + dts[5] = 1 / 5 + dts[6] = 1 / 5 + + mis = @inferred MultirateInfinitesimalSplitMethod(A, D, G, c) + + # fourth-order accurate + series_integrator = @inferred bseries(mis, 4) + series_exact = @inferred ExactSolution(series_integrator) + @test mapreduce(≈, &, values(series_integrator), values(series_exact)) + + # not fifth-order accurate + series_integrator = @inferred bseries(mis, 5) + series_exact = @inferred ExactSolution(series_integrator) + @test mapreduce(≈, &, values(series_integrator), values(series_exact)) == false + end + + @testset "EB4" begin + # Krogstad, Stein + # "Generalized integrating factor methods for stiff PDEs". + # Journal of Computational Physics 203, no. 1 (2005): 72-88. + # https://doi.org/10.1016/j.jcp.2004.08.006 + ns = 4 + A = zeros(Polynomial{Float64, :x}, ns+1, ns) + D = zeros(ns+1, ns) + G = zeros(ns+1, ns) + c = zeros(ns+1) + dts = zeros(ns+1) + + fac1 = 1.0 + fac2 = 2.0 + + A[2,1] = Polynomial([0.5]) + + A[3,1] = Polynomial([0.5, -1 / fac1]) + A[3,2] = Polynomial([0.0, 1 / fac1]) + + A[4,1] = Polynomial([1.0, -2 / fac1]) + A[4,3] = Polynomial([0.0, 2 / fac1]) + + A[5,1] = Polynomial([1.0, -3 / fac1, 4 / fac2]) + A[5,2] = Polynomial([0.0, 2 / fac1, -4 / fac2]) + A[5,3] = Polynomial([0.0, 2 / fac1, -4 / fac2]) + A[5,4] = Polynomial([0.0, -1 / fac1, 4 / fac2]) + + c[1] = 0 + c[2] = 0.5 + c[3] = 0.5 + c[4] = 1 + c[5] = 1 + + mis = @inferred MultirateInfinitesimalSplitMethod(A, D, G, c) + + # fourth-order accurate + series_integrator = @inferred bseries(mis, 4) + series_exact = @inferred ExactSolution(series_integrator) + @test mapreduce(≈, &, values(series_integrator), values(series_exact)) + + # not fifth-order accurate + series_integrator = @inferred bseries(mis, 5) + series_exact = @inferred ExactSolution(series_integrator) + @test mapreduce(≈, &, values(series_integrator), values(series_exact)) == false + end + + @testset "OwrenRK4" begin + # Owren, Brynjulf. + # "Order conditions for commutator-free Lie group methods". + # Journal of Physics A: Mathematical and General 39, no. 19 (2006): 5585. + # https://doi.org/10.1088/0305-4470/39/19/S15 + ns = 5 + A = zeros(Polynomial{Float64, :x}, ns+1, ns) + D = zeros(ns+1, ns) + G = zeros(ns+1, ns) + c = zeros(ns+1) + dts = zeros(ns+1) + + A[2,1] = Polynomial([1 / 2]) + + A[3,2] = Polynomial([1 / 2]) + + A[4,1] = Polynomial([-1 / 2]) + A[4,3] = Polynomial([1.0]) + + A[5,1] = Polynomial([1 / 4]) + A[5,2] = Polynomial([1 / 6]) + A[5,3] = Polynomial([1 / 6]) + A[5,4] = Polynomial([-1 / 12]) + + A[6,1] = Polynomial([-1 / 12]) + A[6,2] = Polynomial([1 / 6]) + A[6,3] = Polynomial([1 / 6]) + A[6,4] = Polynomial([1 / 4]) + + D[4,2] = 1 + D[6,5] = 1 + + dts[1] = 0 + dts[2] = 1 / 2 + dts[3] = 1 / 2 + dts[4] = 1 / 2 + dts[5] = 1 / 2 + dts[6] = 1 / 2 + + mis = @inferred MultirateInfinitesimalSplitMethod(A, D, G, c) + + # fourth-order accurate + series_integrator = @inferred bseries(mis, 4) + series_exact = @inferred ExactSolution(series_integrator) + @test mapreduce(≈, &, values(series_integrator), values(series_exact)) + + # not fifth-order accurate + series_integrator = @inferred bseries(mis, 5) + series_exact = @inferred ExactSolution(series_integrator) + @test mapreduce(≈, &, values(series_integrator), values(series_exact)) == false + end + + @testset "MISBI4_4" begin + # Knoth, Oswald, and Joerg Wensch. + # "Generalized split-explicit Runge–Kutta methods for the compressible Euler equations". + # Monthly Weather Review 142, no. 5 (2014): 2067-2081. + # https://doi.org/10.1175/MWR-D-13-00068.1 + # Table 4, method MIS4a + ns = 4 + A = zeros(Polynomial{Float64, :x}, ns+1, ns) + D = zeros(ns+1, ns) + G = zeros(ns+1, ns) + c = zeros(ns+1) + dts = zeros(ns+1) + + # β in Knoth, Wensch (2014) + A[2,1] = Polynomial([1 / 2]) + A[2,1] = Polynomial([0.38758444641450318]) + A[3,1] = Polynomial([-2.5318448354142823e-002]) + A[3,2] = Polynomial([0.38668943087310403]) + A[4,1] = Polynomial([0.20899983523553325]) + A[4,2] = Polynomial([-0.45856648476371231]) + A[4,3] = Polynomial([0.43423187573425748]) + A[5,1] = Polynomial([-0.10048822195663100]) + A[5,2] = Polynomial([-0.46186171956333327]) + A[5,3] = Polynomial([0.83045062122462809]) + A[5,4] = Polynomial([0.27014914900250392]) + + # seems to be unused? + c[2] = 0.38758444641450318 + c[3] = 0.61521685655017821 + c[4] = 0.23254717315441453 + c[5] = 1.0000000000000002 + + # α in Knoth, Wensch (2014) + D[3,2] = 0.52349249922385610 + D[4,2] = 1.1683374366893629 + D[4,3] = -0.75762080241712637 + D[5,2] = -3.6477233846797109e-002 + D[5,3] = 0.56936148730740477 + D[5,4] = 0.47746263002599681 + + # γ in Knoth, Wensch (2014) + G[3,2] = 0.13145089796226542 + G[4,2] = -0.36855857648747881 + G[4,3] = 0.33159232636600550 + G[5,2] = -6.5767130537473045e-002 + G[5,3] = 4.0591093109036858e-002 + G[5,4] = 6.4902111640806712e-002 + + for i in 1:ns+1 + for j in 1:i-1 + dts[i] = dts[i] + A[i,j][0] + end + end + + mis = @inferred MultirateInfinitesimalSplitMethod(A, D, G, c) + + # third-order accurate + series_integrator = @inferred bseries(mis, 3) + series_exact = @inferred ExactSolution(series_integrator) + @test_broken mapreduce(≈, &, values(series_integrator), values(series_exact)) + + # not fourth-order accurate + series_integrator = @inferred bseries(mis, 4) + series_exact = @inferred ExactSolution(series_integrator) + @test mapreduce(≈, &, values(series_integrator), values(series_exact)) == false + end +end + + end # @testset "BSeries"