Developer documentation

Warning

This page specifies developer extension APIs for packages implementing operator-splitting solvers. They are versioned for the OrdinaryDiffEq ecosystem, but are not supported end-user APIs. Application code should use the APIs in the API Reference instead.

Synchronizers API

A key part of operator splitting algorithms is the synchronization logic. Parameters of one subproblem might need to be kept in sync with the solution of other subproblems and vice versa. To handle this efficiently OrdinaryDiffEqOperatorSplitting.jl provides a small set of utils.

OrdinaryDiffEqOperatorSplitting.NoExternalSynchronization — Type

NoExternalSynchronization()

Indicator that no synchronization between parameters and solution vectors is necessary.

source

OrdinaryDiffEqOperatorSplitting.forward_sync_subintegrator! — Function

forward_sync_subintegrator!(parent_integrator::AnySplitIntegrator, inner_integrator::DEIntegrator, solution_indices, sync)

Synchronize one child with its parent immediately before advancing that child.

Arguments

parent_integrator: Splitting node that owns the current full solution.
inner_integrator: Direct leaf child to receive its local state.
solution_indices: Indices selecting the child state in the parent's solution.
sync: Synchronizer object for any parameter or cross-child synchronization.

Interface requirements

The built-in method copies the parent's selected state into inner_integrator and marks its derivative cache stale. To synchronize parameters or coupled child state, implement forward_sync_external!(parent_integrator, inner_integrator, sync) for the concrete sync type.

Returns

nothing.

source

OrdinaryDiffEqOperatorSplitting.backward_sync_subintegrator! — Function

backward_sync_subintegrator!(parent_integrator::AnySplitIntegrator, inner_integrator::DEIntegrator, solution_indices, sync)

Synchronize one advanced child back into its parent.

Arguments

parent_integrator: Splitting node that owns the full solution.
inner_integrator: Direct leaf child whose state is copied back.
solution_indices: Indices selecting the child state in the parent's solution.
sync: Synchronizer object for parameter or cross-child synchronization.

Interface requirements

The built-in method copies the child state into the parent's selected state. To synchronize parameters or coupled child state, implement backward_sync_external!(parent_integrator, inner_integrator, sync) for the concrete sync type.

Returns

nothing.

source

OrdinaryDiffEqOperatorSplitting.forward_sync_external! — Function

forward_sync_external!(parent_integrator, inner_integrator, sync)

Synchronize external state into a child immediately before that child advances.

Arguments

parent_integrator: Splitting parent that owns the full state.
inner_integrator: Direct child receiving the synchronized state.
sync: Synchronizer object selected in the GenericSplitFunction tree.

Interface requirements

Implement this method for a concrete sync type when parameters or coupled child state must be refreshed before the child solves. The default NoExternalSynchronization implementation is a no-op.

Returns

nothing.

This is a developer extension API, not a supported end-user API.

source

OrdinaryDiffEqOperatorSplitting.backward_sync_external! — Function

backward_sync_external!(parent_integrator, inner_integrator, sync)

Synchronize external state after a child has advanced.

Arguments

parent_integrator: Splitting parent that owns the full state.
inner_integrator: Direct child whose result may update coupled state.
sync: Synchronizer object selected in the GenericSplitFunction tree.

Interface requirements

Implement this method for a concrete sync type when parameters or coupled child state must be refreshed after the child solves. The default NoExternalSynchronization implementation is a no-op.

Returns

nothing.

This is a developer extension API, not a supported end-user API.

source

OrdinaryDiffEqOperatorSplitting.need_sync — Function

need_sync(a, b)

Return whether copying solution information from b into a is necessary.

Arguments

a: Destination vector or view.
b: Source vector or view.

Returns

false when both arguments provably share the same backing storage and copying would be redundant; true otherwise (including whenever that cannot be decided cheaply). SubArray arguments are compared by parent identity. Plain DenseVectors – which includes a GPU array package's own contiguous view, since those are commonly returned as a new dense wrapper over the same buffer rather than a SubArray – are compared by base pointer and length: equal on both means the same memory. Extend this function for custom array wrappers whose aliasing relationship cannot be determined by the built-in methods.

source

OrdinaryDiffEqOperatorSplitting.sync_vectors! — Function

sync_vectors!(a, b)

Copy solution information from b into a when need_sync determines that their storage does not alias.

Arguments

a: Destination vector or view.
b: Source vector or view.

Returns

nothing. If no copy is required, a is left unchanged.

source

Adding Synchronizers

Warning

The API is not stable yet and subject to breaking changes.

You need to provide dispatches for

OrdinaryDiffEqOperatorSplitting.forward_sync_subintegrator! — Function

forward_sync_subintegrator!(parent_integrator::AnySplitIntegrator, inner_integrator::DEIntegrator, solution_indices, sync)

Synchronize one child with its parent immediately before advancing that child.

Arguments

parent_integrator: Splitting node that owns the current full solution.
inner_integrator: Direct leaf child to receive its local state.
solution_indices: Indices selecting the child state in the parent's solution.
sync: Synchronizer object for any parameter or cross-child synchronization.

Interface requirements

Returns

nothing.

source

OrdinaryDiffEqOperatorSplitting.backward_sync_subintegrator! — Function

backward_sync_subintegrator!(parent_integrator::AnySplitIntegrator, inner_integrator::DEIntegrator, solution_indices, sync)

Synchronize one advanced child back into its parent.

Arguments

parent_integrator: Splitting node that owns the full solution.
inner_integrator: Direct leaf child whose state is copied back.
solution_indices: Indices selecting the child state in the parent's solution.
sync: Synchronizer object for parameter or cross-child synchronization.

Interface requirements

Returns

nothing.

source

with your custom synchronizer object and add it to the split function construction as follows:

f1, f2 = generate_individual_functions() # assuming 3 unknowns each
i1, i2 = generate_solution_indices()     # e.g. ([1,2,3], Int[])
synchronizer_tree = generate_my_synchronizer_tree() # e.g. (MySynchronizer([1,2,3]), NoExternalSynchronization())
f = GenericSplitFunction((f1, f2), (i1, i2), synchronizer_tree)
u0 = [-1.0, 1.0, 0.0]
tspan = (0.0, 1.0)
prob = OperatorSplittingProblem(f, u0, tspan)

Solver extension API

OrdinaryDiffEqOperatorSplitting.AbstractOperatorSplittingAlgorithm — Type

AbstractOperatorSplittingAlgorithm

Developer-only abstract supertype for algorithms that advance an OperatorSplittingProblem.

Interface requirements

Store an inner_algs tuple whose shape mirrors the associated GenericSplitFunction.
Implement init_cache to construct a cache for every node.
Implement _perform_step! to advance that node.

This interface is versioned for OrdinaryDiffEq solver developers. It is not a supported end-user API.

source

OrdinaryDiffEqOperatorSplitting.AbstractOperatorSplittingCache — Type

AbstractOperatorSplittingCache

Developer-only abstract supertype for an algorithm's per-node cache. A concrete subtype holds references to the node's u and uprev buffers and any additional temporary storage required by the splitting scheme. Construct it from init_cache.

This interface is versioned for OrdinaryDiffEq solver developers. It is not a supported end-user API.

source

OrdinaryDiffEqOperatorSplitting.init_cache — Function

init_cache(f::GenericSplitFunction, alg::AbstractOperatorSplittingAlgorithm; uprev, u)

Construct the per-node cache for a developer-defined operator-splitting algorithm.

Arguments

f: Operator tree at the node being initialized.
alg: Algorithm used at that node.

Keyword Arguments

uprev: Mutable state buffer for the preceding accepted state.
u: Mutable state buffer for the state currently being advanced.

Returns

A concrete AbstractOperatorSplittingCache. The cache must retain the provided buffers by reference; the integrator owns their allocation and restores them after rejected steps.

This is a developer extension API, not a supported end-user API.

source

OrdinaryDiffEqOperatorSplitting._perform_step! — Function

_perform_step!(parent, children::Tuple, cache::AbstractOperatorSplittingCache, dt)

Advance one node of an operator-splitting tree by dt.

Arguments

parent: Integrator node that owns the current solution and rollback buffers.
children: Direct child integrators, either DEIntegrators or nested splitting nodes.
cache: Cache returned by init_cache for parent's algorithm.
dt: Signed duration of the splitting step.

Interface requirements

Synchronize a child before and after advancing it.
Leave parent.uprev unchanged; rejection restores the node from that buffer.
Set parent.force_stepfail = true and return immediately when a child fails.
For adaptive algorithms, pass the tolerance-scaled local error estimate to OrdinaryDiffEqCore.set_EEst! when parent.controller_cache !== nothing.

This is a developer extension API, not a supported end-user API.

source

OrdinaryDiffEqOperatorSplitting.advance_solution_by! — Function

advance_solution_by!(parent, child, dt)

Advance one direct child of an operator-splitting node by the signed duration dt.

Arguments

parent: The node that owns child and coordinates the splitting step.
child: A leaf DEIntegrator or nested splitting node.
dt: Signed duration to apply to the child.

Returns

nothing. This function only advances the child; a custom _perform_step! must synchronize the child before and after the call and use child_failed before accepting its state.

This is a developer extension API, not a supported end-user API.

source

OrdinaryDiffEqOperatorSplitting.child_failed — Function

child_failed(child)

Return whether a direct child integrator has failed while a developer-defined splitting step is executing.

Arguments

child: A leaf DEIntegrator or nested operator-splitting integrator.

Returns

true when the child cannot be used for the remainder of the current splitting step. In that case, _perform_step! must set parent.force_stepfail = true and return without synchronizing the failed state back to its parent.

This is a developer extension API, not a supported end-user API.

source

OrdinaryDiffEqOperatorSplitting.alg_adaptive_order — Function

alg_adaptive_order(alg::AbstractOperatorSplittingAlgorithm)

Order of the error estimator of an adaptive operator splitting algorithm; every algorithm with SciMLBase.isadaptive(alg) == true has to implement this.

source

OrdinaryDiffEqOperatorSplitting.splitting_interpolant — Function

splitting_interpolant(integrator, cache, Θ, dt, y₀, y₁, idxs, ::Type{Val{D}})
splitting_interpolant!(out, integrator, cache, Θ, dt, y₀, y₁, idxs, ::Type{Val{D}})

Dense output of a single operator splitting step, evaluated in the step-local coordinate Θ = (t - tprev) / dt, where y₀/y₁ are the state at tprev/t and dt is the step that was taken. D is derivative order 0 or 1, and idxs selects components (nothing for all of them).

This is the extension point for higher order dense output: dispatch on the algorithm's cache type. The fallback for any AbstractOperatorSplittingCache is linear interpolation.

source

OrdinaryDiffEqOperatorSplitting.splitting_interpolant! — Function

splitting_interpolant(integrator, cache, Θ, dt, y₀, y₁, idxs, ::Type{Val{D}})
splitting_interpolant!(out, integrator, cache, Θ, dt, y₀, y₁, idxs, ::Type{Val{D}})

This is the extension point for higher order dense output: dispatch on the algorithm's cache type. The fallback for any AbstractOperatorSplittingCache is linear interpolation.

source

Adding solvers

Warning

This is a developer extension API. It may change independently of the end-user API and should not be used in application code.

To add a new solver, define two structs – one describing the algorithm, one for its cache – and dispatch the developer extension functions on them:

init_cache(f, alg; uprev, u) builds the cache for one node of the splitting tree.
_perform_step!(parent, children, cache, dt) advances that node by dt.
alg_adaptive_order(alg), only if the algorithm is adaptive (see below).

The algorithm struct has to carry the inner algorithms of the problem sequence in a field named inner_algs, because that is how the tree of integrators is built alongside the tree of split functions. The cache is where a scheme keeps the buffers it needs beyond u/uprev.

using SciMLBase, OrdinaryDiffEqOperatorSplitting
import OrdinaryDiffEqOperatorSplitting as OS

struct MySimpleFirstOrderAlgorithm{InnerAlgorithmTypes} <:
    OS.AbstractOperatorSplittingAlgorithm
    inner_algs::InnerAlgorithmTypes # Tuple of solvers for the problem sequence
end

struct MySimpleFirstOrderCache{uType, uprevType} <: OS.AbstractOperatorSplittingCache
    u::uType
    uprev::uprevType
end

function OS.init_cache(
        f::GenericSplitFunction, alg::MySimpleFirstOrderAlgorithm;
        uprev::AbstractArray, u::AbstractVector
    )
    # `u` and `uprev` are the buffers of *this* node; the integrator owns them and
    # the cache only keeps references. Allocate additional buffers with
    # `similar(u)` if the scheme needs them.
    return MySimpleFirstOrderCache(u, uprev)
end

The stepping function receives the node whose step is being performed (parent) and the tuple of its child integrators, which may be leaf DEIntegrators or nested SplitSubIntegrators – the same code handles both. Everything else a step needs is reachable from parent: parent.child_solution_indices[i] are the indices of the i-th child in this node's solution vector, and parent.child_synchronizers[i] is its synchronizer.

Advancing one child means synchronizing into it, stepping it, and synchronizing back:

function advance_one_child!(parent, child, i, dt)
    idxs = parent.child_solution_indices[i]
    sync = parent.child_synchronizers[i]

    OS.forward_sync_subintegrator!(parent, child, idxs, sync)
    OS.advance_solution_by!(parent, child, dt)
    if OS.child_failed(child)
        # A failed child must stop the remaining stages of this step.
        parent.force_stepfail = true
        return
    end
    OS.backward_sync_subintegrator!(parent, child, idxs, sync)
    return
end

function OS._perform_step!(
        parent, children::Tuple, cache::MySimpleFirstOrderCache, dt
    )
    advance_one_child!(parent, children[1], 1, dt)
    parent.force_stepfail && return

    advance_one_child!(parent, children[2], 2, dt)
    parent.force_stepfail && return

    # Done :) The solution of the step is in `parent.u`; `parent.uprev` holds the
    # state at the beginning of the step and must be left untouched, as rollback
    # after a rejected step restores from it.
    return
end

This example is written for exactly two operators. A scheme that works for any number of them loops over children with Unrolled.@unroll, as the built-in algorithms in src/solver.jl do; a plain for loop over the heterogeneously typed tuple would be type unstable.

Adaptive algorithms

A splitting node runs a step size controller only if its algorithm both declares itself adaptive and produces an error estimate. Such an algorithm has to

define SciMLBase.isadaptive(::MyAlgorithm) = true,
define OrdinaryDiffEqOperatorSplitting.alg_adaptive_order, the order of its error estimator, and
hand the tolerance-scaled error estimate to OrdinaryDiffEqCore.set_EEst! at the end of _perform_step!, following the OrdinaryDiffEq convention that a step is accepted when the estimate is <= 1.

Only do the last step when parent.controller_cache !== nothing: the node may have been configured non-adaptive, in which case no controller consumes the estimate. The tolerances and the norm to scale with live in parent.opts:

if parent.controller_cache !== nothing
    (; abstol, reltol, internalnorm) = parent.opts
    @. residual = error_of_the_step / (abstol + max(abs(parent.u), abs(parent.uprev)) * reltol)
    OrdinaryDiffEqCore.set_EEst!(parent, internalnorm(residual, parent.t + dt))
end

set_EEst! and its counterpart OrdinaryDiffEqCore.get_EEst are the public OrdinaryDiffEqCore interface for the estimate; go through them rather than touching the EEst field, whose location on the integrator is an implementation detail.

See PalindromicPairLieTrotterGodunov in src/solver.jl for a complete example, and Adaptive time stepping for how the two layers of adaptivity interact.

Dense output

Saving, saveat and continuous callback root-finding all go through a single hook, so an algorithm only has to describe how to interpolate within one of its steps:

function OrdinaryDiffEqOperatorSplitting.splitting_interpolant(
        integrator, cache::MySimpleFirstOrderCache, Θ, dt, y₀, y₁, idxs, ::Type{Val{D}}
) where {D}
    # Θ = (t - tprev) / dt is the step-local coordinate, y₀ = u(tprev), y₁ = u(t).
    # ...
end

function OrdinaryDiffEqOperatorSplitting.splitting_interpolant!(
        out, integrator, cache::MySimpleFirstOrderCache, Θ, dt, y₀, y₁, idxs,
        ::Type{Val{D}}
) where {D}
    # In-place variant; must return `out`.
end

Both fall back to linear interpolation for any AbstractOperatorSplittingCache, so implementing them is optional – but with the fallback, saved output is first order even for a second-order scheme, and continuous callback event times are the exact roots of a straight chord across the step. One method improves both.

The fallback is linear for a structural reason, and it is the same reason saving and callbacks live on the outer integrator alone: a splitting step advances its children sequentially over staggered subintervals, so a child's own interpolant describes a different sub-problem over a different interval, and they do not compose into an approximation of the split solution. Only the step endpoints, which the outer integrator owns, are states of the full split system – an inner split is a stage, not a step.

When implementing this, note that Θ is derived from integrator.t - integrator.tprev, not from integrator.dt: once a step is accepted, step_accept_controller! has already replaced dt with the step size proposed for the next step.