API

AbstractPararealBackend

How the fine sweep of Parareal executes. The algorithm lives ONCE in parareal/parareal.jl; a backend supplies only these primitives:

• is_root(backend) — whether this process runs the serial parts (coarse init, correction, convergence check). true everywhere except non-root MPI ranks.
• sync_starts!(backend, cache, k) — make the chunk starting values U[n] available wherever chunk n will be solved. No-op off MPI.
• fine_map!(backend, cache, solution, problem, parareal, k; …) — solve chunks k:N with the fine solver; deposit boundary values into cache.F[n+1] on the root and chunk solutions into solution.lastiterate (MPI ranks keep their chunk local until collect_chunks!).
• snapshot!(backend, solution, k; …) — record iterate k when saveiterates is on. In memory by default; on MPI the ranks have already written their chunk to the run directory inside fine_map!.
• sync_converged(backend, flag) — agree the convergence decision.
• collect_chunks!(backend, solution, cache; …) — bring final chunks (and, under MPI with saveiterates, the per-iterate history) to the root.

An abstract type for parameters of Parareals.

An abstract type for caching intermediate computations in AbstractTimeParallelSolvers.

An abstract type for iterates in AbstractTimeParallelSolutions.

An abstract type for parameters in AbstractTimeParallelSolvers.

An abstract type for time-parallel solutions of NSDEBase.AbstractInitialValueProblems.

An abstract type for time-parallel solvers of NSDEBase.AbstractInitialValueProblems.

An abstract type for tolerance parameters in AbstractTimeParallelSolvers.

An abstract type for proximity-function weights used by AbstractTimeParallelSolvers.

Fine chunks via pmap over the default worker pool (Distributed.jl).

SPMD over MPI ranks: rank r owns chunk r + 1. Requires MPI.Init().

Parareal <: AbstractTimeParallelSolver

A composite type for the Parareal algorithm.

Constructors

Parareal(finesolver, coarsesolver, parameters, tolerance)
Parareal(finesolver, coarsesolver; parameters=PararealParameters(), tolerance=Tolerance())

Arguments

• finesolver :: AbstractInitialValueSolver : fine solver (accurate but expensive).
• coarsesolver :: AbstractInitialValueSolver : coarse solver (rough but quick).
• parameters :: AbstractPararealParameters : parameters for the correction step.
• tolerance :: AbstractTolerance : tolerance and error mechanism.

Methods

(parareal::Parareal)(solution::PararealSolution, problem::AbstractInitialValueProblem)
(parareal::Parareal)(problem::AbstractInitialValueProblem)

returns the solution of a problem using parareal.

PararealCache <: AbstractTimeParallelCache

Pre-allocated state for one Parareal run: the chunk boundary values U, the fine and coarse boundary results F and G, the previous-iterate snapshot U_ for the error function, the chunk time grid T, the chunk problems (built once, their u0 slots referencing U so the correction updates them in place), and reusable solver caches.

PararealIterate <: AbstractTimeParallelIterate

A composite type for a single iterate in a PararealSolution: one fine chunk solution per time chunk. The chunk storage is CONCRETE — chunks are built eagerly from the fine solver, so eltype(iterate.chunks) is the fine solver's concrete solution type, not an abstract box. This matters: these objects cross thread and process boundaries in the hot path.

Constructors

PararealIterate(chunks::AbstractVector{𝕊}) where 𝕊<:AbstractInitialValueSolution
PararealIterate(problem::AbstractInitialValueProblem, parareal::Parareal)

Functions

• firstindex/getindex/lastindex/setindex! : chunk access
• length/numchunks : number of chunks

Methods

(iterate::PararealIterate)(t::Real)

returns the value of iterate at t via interpolation of the owning chunk.

PararealParameters <: AbstractPararealParameters

A composite type for the basic parameters of Parareal.

Constructors

PararealParameters(N, K)
PararealParameters(; N=10, K=N)

Arguments

• N :: Integer : number of time chunks/processors.
• K :: Integer : maximum number of iterations.

PararealSolution <: AbstractTimeParallelSolution

A composite type for an AbstractTimeParallelSolution obtained using Parareal.

Constructors

PararealSolution(lastiterate, errors)
PararealSolution(problem::AbstractInitialValueProblem, parareal::Parareal)

Arguments

• lastiterate :: PararealIterate
• errors :: AbstractVector{ℝ} where ℝ<:Real : iteration errors.

Functions

• firstindex : first index.
• getindex : get iterate.
• lastindex : last index.
• numiterates : number of iterates.
• numchunks : number of chunks of last iterate.
• setindex! : set iterate.

Methods

(solution::PararealSolution)(t::Real)

returns the value of solution at t via interpolation.

PipelinedMPIBackend <: AbstractPararealBackend

Pipelined (task-scheduled) Parareal over MPI ranks, one rank per chunk — the scheduling of Aubanel (thesis §2.7, eq:aubanel_parallel_efficiency). There is no root: the serial coarse chain is passed rank to rank, and each rank runs its fine solve for the next sweep BEFORE blocking on the corrected start from its left neighbour, so the fine work hides the serial coarse cost. Judge measured speed-ups against theoretical_speedup(K, N, ζ; estimate = :aubanel).

This is a SEMANTIC VARIANT of Parareal, not a drop-in adapter. Differences from the primitive-based backends:

• Frontier convergence, two flavours. Convergence sweeps left to right as a frontier; a chunk stops once every chunk to its left has accepted AND its own per-chunk test passes, or it hits the finite-termination diagonal (k = n). The per-chunk test depends on tolerance.ψ:

  • ψ = ψ∞ — the CERTIFIED thesis criterion (§3.4 subsec:local_proximity): the test is the weighted CURRENT defect of the previous iterate at this rank's outgoing boundary, w^(T[1] − T[n+1]) ‖uold − F(start)‖ ≤ ϵ, where uold is the boundary value sent last sweep and F(start) this sweep's fine solve from the very start that produced it — the thesis's one-stage lag (rem:pipelined_protocol). The acceptance flag is the prefix conjunction A_n = A_{n−1} ∧ tₙ, riding the existing messages; when the run accepts, the assembled vector satisfies ψ∞ ≤ ϵ, hence ‖U − U*‖_{W,∞} ≤ s(θ_F) ϵ (thm:equivalence_error_norm_local) — set ϵ = r/s(θ_F) to certify radius r. Requires a PRESET scalar weight: weights.updatew = true is rejected (no root sees all boundaries at a common sweep). Acceptance is suppressed for one sweep whenever the incoming start changed in the same message that raised the flag (an upstream diagonal exit): the test must be re-run against the frozen start, or the certificate would mix iterates — the thesis's freezing cascade, at most one extra sweep per chunk.
  • anything else — relative STAGNATION of the outgoing boundary (update ≤ tolerance.ϵ), the pre-thesis heuristic; the global ψ and its weights are ignored, with a warning unless ψ = ψ₁. Kept as the default for backward compatibility.

• solution.errors[k] records, among chunks still active at sweep k, the largest weighted defect (ψ∞ mode: the running restriction of ψ∞ to active chunks) or the largest relative update (stagnation mode) — a frontier diagnostic either way.

• Chunks may stop at different sweeps, so numiterates(solution) counts the sweeps of the LAST chunk to converge.

• saveiterates is not supported: recording every sweep would force the per-sweep synchronisation the pipeline exists to avoid.

Accepted chunks remain fine solves from their accepted starts (one final fine solve is issued whenever a start arrived after the last sweep's), so flatten, seams and every downstream statistic read exactly as for the other backends. At ϵ = 0, K = N, finite termination makes the result the chunked fine solve, bitwise — the same fixed point as SerialBackend, reached by a different route.

Requires MPI.Init() and at least N ranks; ranks beyond N idle through the chain and join the final gather. Select with mode = "PIPELINED" or by passing the backend object; pair with Tolerance(ϵ = r/s, ψ = ψ∞, weights = Weights(w = ŵ)) for the certified criterion (ŵ preset — the updater is rejected here).

Runs everything in the calling task. The reference implementation.

Fine chunks via Threads.@threads (shared memory). No schedule argument is given on purpose: from Julia 1.8 the default is :dynamic, which is what we want, and on 1.6–1.7 the argument does not exist — there the loop is scheduled statically, which is correct and merely balances less well when chunks differ in cost.

Tolerance <: AbstractTolerance

A composite type for the tolerance mechanism of an time-parallel solver.

Constructors

Tolerance(ϵ, ψ, weights)
Tolerance(; ϵ=1e-12, ψ=ψ₁, weights=Weights())

Arguments

• ϵ :: Real : tolerance.
• ψ :: Function : error function.
• weights :: Weights : weights for ψ.

Weights <: AbstractWeights

A composite type for the weights of Tolerance.

Constructors

Weights(; w=1.0, updatew=false, δ=1.0)

Arguments

• w :: Real : weighting factor for ψ, a PER-UNIT-TIME base (w = exp(Λ) for per-time rate Λ). A finite, positive scalar. (Earlier versions admitted a vector here; nothing downstream — ψ₂, ψ∞, update!, the MoWi zoom — ever defined what a per-chunk vector of bases meant, and each threw a MethodError on one. A vector is now refused at construction.)
• updatew :: Bool : flags when to update! w using (an approximation of) the Lipschitz function of the fine solver. Bulk-synchronous backends only: under PipelinedMPIBackend with ψ∞ there is no root that sees all boundaries at a common sweep, so the updater cannot run — that path REJECTS updatew = true; preset w (a known rate, or a frozen probe measurement — see the criterion docs).
• δ :: Real : safety divisor applied to the MEASURED rate only (δ < 1 inflates the measured base — the thesis's safety factor w = C·Λ̂ with C = 1/δ). It never touches the running w, so with updatew = false update! is a strict no-op regardless of δ.

Functions

• update! : updates w using (an approximation of) the Lipschitz function of the fine solver.

firstindex(iterate::PararealIterate)

returns the first index of iterate.

firstindex(solution::PararealSolution)

returns the first index of solution.

getindex(iterate::PararealIterate, n::Integer)

returns the n-th chunk of iterate.

getindex(solution::PararealSolution, n::Integer)

returns the n-th chunk of the last iteration of a PararealSolution.

lastindex(iterate::PararealIterate)

returns the last index of iterate.

lastindex(solution::PararealSolution)

returns the last index of solution.

length(iterate::PararealIterate)

returns the number of chunks of iterate.

length(solution::PararealSolution)

returns the number of chunks of solution.

setindex!(iterate::PararealIterate, value::AbstractInitialValueSolution, n::Integer)

stores value into the n-th chunk of iterate.

setindex!(solution::PararealSolution, chunk::AbstractInitialValueSolution, n::Integer)

stores an AbstractInitialValueSolution as the n-th chunk of the last iteration of a PararealSolution.

NSDEBase.initialize_cache(problem, parareal::Parareal) :: PararealCache

returns a reusable PararealCache for problem, honouring the same contract NSDERungeKutta honours for its solvers. Build once, then call parareal(cache, solution, problem; ...) repeatedly — e.g. in timing loops, where rebuilding the cache per solve measures the allocator, not the algorithm.

NSDEBase.initialize_solution(problem, parareal::Parareal; saveiterates=false) :: PararealSolution

returns a PararealSolution for problem. Reusable across solves.

solve!(solution::AbstractTimeParallelSolution, problem, solver::AbstractTimeParallelSolver; kwargs...) :: AbstractTimeParallelSolution

computes the solution of problem using solver.

solve(problem, solver::AbstractTimeParallelSolver; kwargs...) :: AbstractTimeParallelSolution

computes the solution of problem using solver.

TimeParallelSolution(problem::AbstractInitialValueProblem, parareal::Parareal)

returns a PararealSolution constructor for the solution of problem with parareal.

Wnorm(iterate::PararealIterate, reference::AbstractInitialValueSolution, w::Number)

weighted distance between iterate and a reference solution at the chunk boundaries.

boundarytimes(iterate::PararealIterate)

returns the N + 1 chunk-boundary times of iterate: each chunk's start, plus the final chunk's end. These are the points Parareal iterates on — the grid every convergence measure (Wnorm, per-iterate error studies) is taken over.

boundaryvalues(iterate::PararealIterate)

returns the N + 1 chunk-boundary values of iterate: each chunk's first state, plus the final chunk's last state — the U vector of the Parareal iteration, read off the chunk solutions. The returned states ALIAS the chunk storage; copy them before mutating.

coarseinit!(cache, parareal)

runs the serial coarse pass over the chunk grid, filling G and — where makeGs allows — seeding the chunk starting values U. Chunks whose makeGs[n] is false keep whatever U[n] already holds (injected guesses, e.g. from a moving-window driver).

collect_iterates!(solution::PararealSolution; directory)

reads the per-iterate chunk files written by the MPI backend under saveiterates back into solution.iterates. Chunks below the diagonal (n < k) are final from earlier iterations and are copied forward.

contractionrate(errors::AbstractVector{<:Real}) :: Float64
contractionrate(solution::AbstractTimeParallelSolution) :: Float64

the observed per-iteration contraction rate β of an error trace: the least-squares slope of log(errors[k]) against k, exponentiated, so that errors[k] ≈ C βᵏ. Non-finite entries and EXACT ZEROS are excluded before fitting — finite termination drives the last error to bitwise zero, which is a property of the algorithm's fixed point, not of its rate. Returns NaN when fewer than two usable points remain: a run that converges in one iteration exhibits no observable rate, and downstream bounds (e.g. the §4.3 outer radius R = r/βᴷ) must treat that as unbounded rather than fabricate a number.

correct!(cache, parareal, k)

the serial Parareal correction sweep after fine batch k: re-run the coarse solver from the freshest starts and update U[n+1] = F[n+1] + (Gnew − Gold), writing IN PLACE (ping-pong against the U_ snapshot taken by the error function — the old code rebound a fresh array here every step to protect that history).

Two details make finite termination exact in floating point, not just in real arithmetic:

• Order of operations. F + (Gnew − Gold) returns F bitwise whenever the two coarse values cancel bitwise; (Gnew + F) − Gold does not — it rounds F into Gnew first and loses low bits (with Gnew = Gold = 1e16 and F = 1 it returns 0). Do not "simplify" this expression.
• The newly exact interface is assigned, not computed. After batch k the start U[k] is already exact, so F[k+1] is the exact value of U[k+1] and U[k] has not moved since G[k+1] was last computed: the coarse re-run at n = k would reproduce G[k+1] bitwise. Skip it and copy F[k+1] across — one coarse solve saved per sweep, and no arithmetic on the exact prefix at all.

costratio(problem, parareal; repeats=3) :: Real

measures ζ, the coarse/fine cost ratio over ONE chunk of problem under parareal's chunking — the input theoretical_speedup needs. Both solvers are timed on the first chunk with @elapsed (best of repeats, after a compile warm-up). Measure ζ rather than guessing it from step-size ratios: those ignore per-step cost differences between the two solvers.

flatten(iterate::PararealIterate)

concatenates the chunk solutions into one pair of vectors (u, t), dropping each chunk's FIRST point after chunk 1: boundaries are stored twice (chunk n's last point and chunk n+1's corrected start), and keeping both would double-weight every boundary in any downstream statistic. The kept value is the fine endpoint; the dropped corrected start agrees with it only to the solve tolerance. t is monotone.

flatten(solution::PararealSolution)

flattens the last iterate: one (u, t) pair for the whole span, chunk boundaries deduplicated.

iteration_budget(S::Real, N::Integer, ζ::Real; estimate=:improved) :: Int

the largest per-window iteration count K ∈ 1:N for which Parareal over N chunks with cost ratio ζ still meets the target speed-up S, from the parallel-efficiency estimates behind theoretical_speedup; 0 when the target is infeasible — no K ≥ 1 reaches S, not even a single sweep. The budget is capped at N because the algorithm terminates after N sweeps: a target so modest that every sweep count meets it yields N, not a number beyond the algorithm's range.

• estimate = :improved (default): the work-credited bound, S(K) = N / (K (1 + ζN)(1 − (K − 1)/2N)), which is strictly decreasing on 1:N. It is inverted by integer bisection against theoretical_speedup itself, so the budget agrees with the estimate it is derived from bit for bit. (The closed-form root of the underlying quadratic was used before; its discriminant goes negative exactly when the target is BELOW the K = N speed-up — every K feasible — and the old code read that as "infeasible" and returned 0. It also cancelled badly at large N.)
• estimate = :naive: K = ⌊N/(S(1 + ζN))⌋, likewise clamped to 0:N.

Measure ζ with costratio; a step-size guess distorts the budget exactly where it matters (large ζN).

maxseam(iterate::PararealIterate) :: Real

returns the largest chunk-boundary seam of iterate (see seams), or 0.0 for a single-chunk iterate.

maxseam(solution::PararealSolution) :: Real

maxseam of the last iterate.

numchunks(iterate::PararealIterate)

returns the number of chunks of iterate.

numchunks(solution::PararealSolution)

returns the number of chunks of solution.

numiterates(solution::PararealSolution)

returns the number of iterates of solution.

parareal!(cache, solution, problem, parareal, backend; saveiterates, directory, nocollect)

runs Parareal with the given AbstractPararealBackend. Under MPI this is SPMD: every rank calls it; is_root guards the serial parts.

seams(iterate::PararealIterate) :: Vector{<:Real}

returns the N − 1 chunk-boundary seam sizes of iterate: at each interior boundary, ‖(start of chunk n+1) − (end of chunk n)‖ — the distance between the corrected starting value U[n+1] and the fine endpoint F(U[n]) it should agree with at convergence. This is exactly the discontinuity that flatten hides by dropping the corrected start, so any statistical claim made on flattened output should report maxseam alongside it: under the weighted criterion ψ₂, late-window seams are unconstrained BY DESIGN (the discount forgives them), and only the seam sizes say which regime a run actually exercised. Empty for a single-chunk iterate.

seams(solution::PararealSolution) :: Vector{<:Real}

seams of the last iterate.

shiftwindow!(cache::PararealCache, τ0, τN)

re-targets a cache at the window [τ0, τN]: recomputes the chunk grid T and updates every chunk problem's tspan in place. This is the official seam for window drivers (e.g. NSDEMovingWindow) — chunk problems are built once at cache construction, so their time spans must be moved through this function, never by poking cache.T alone.

theoretical_speedup(K::Integer, N::Integer, ζ::Real; wallclock=false, estimate=:work) :: Real

Parareal's own speed-up ceiling over the serial fine solve, for K iterations over N chunks with cost ratio ζ = (coarse solve time over one chunk) / (fine solve time over one chunk). Measure ζ — costratio does it on one chunk of each — rather than guessing it. Measured speed-ups must be judged against this ceiling, not against the serial fine solve wishfully.

Three estimates (thesis §2.7), selected with estimate:

• :work (default) — the classical WORK bound, whose (1 - (K - 1)/2N) factor credits the shrinking k:N sweeps (eq:improved_parallel_efficiency). That credit is saved work, not saved wall time: with workers ≥ chunks a parallel sweep costs one chunk-time however few chunks remain, so wall-clock measurements at small N (and any K = N point) sit BELOW this bound by construction.
• :wallclock — the work bound with the credit dropped, N / (K (1 + ζN)): the right comparator for a timing benchmark with a full worker pool. wallclock = true is an equivalent spelling, kept for backward compatibility; an explicit estimate wins if both are given.
• :aubanel — the task-scheduled (pipelined) ceiling (eq:aubanel_parallel_efficiency), N / (K (1 - (K - 1)/2N)(1 + ζ) + ζ(N - 1)): the serial coarse chain is passed rank to rank and hidden behind fine work, so the ζN term no longer multiplies K. This is the comparator for PipelinedMPIBackend. At ζ = 0 it coincides with :work — with a free coarse solver the pipeline has nothing to hide.

update!(weights::Weights, U, F, T)

updates weights.w from the chunk-boundary values — bulk-synchronous backends only (the pipelined ψ∞ path has no global view and rejects updatew = true): the largest measured per-chunk amplification ‖F[i+1] − F[i]‖ / ‖U[i] − U[i−1]‖, converted to a PER-UNIT-TIME base via r^(1/(T[i+1] − T[i])). The conversion is not optional: ψ₂ consumes w as w^(T[1] − T[n]) — an exponent in TIME — so w must be exp(Λ) with Λ a per-time rate. Storing the raw per-chunk ratio (the old behaviour) is right only when chunks are exactly one time unit long, as in the Lorenz thesis set-up; with 10-unit chunks the discount runs at ten times the honest rate, and ψ₂ goes blind past the first boundary.

ψ₁(cache, k, weights)

standard relative error function: the mean over chunk boundaries of ‖U[n] − V[n]‖ / ‖U[n]‖, where V is the previous iterate (the coarse prediction F at k = 1). Ignores weights.

ψ₂(cache, k, weights)

weighted error function: the mean over chunk boundaries of ‖w^(T[1] − T[n]) (U[n] − F[n])‖. With w = exp(Λ) for a problem with Lyapunov exponent Λ, this DISCOUNTS each boundary error at the rate the dynamics amplifies it — the moving-window criterion of the thesis. On long chaotic spans it deliberately forgives late-window errors, so convergence in ψ₂ means "early boundaries settled", not "uniformly within ϵ"; check the terminal error against a serial fine solve when that distinction matters (e.g. in speed-up benchmarks, where K sets the Amdahl ceiling).

ψ∞(cache, k, weights)

local (pipelined) proximity function — the ℓ∞ member of the weighted family (thesis §3.4, subsec:local_proximity): the MAXIMUM over chunk boundaries of ‖w^(T[1] − T[n]) (U[n] − F[n])‖, against ψ₂'s mean. Same defect, same per-unit-time base w, so ψ₂ ≤ ψ∞ ≤ N ψ₂ on any state.

Why it exists — two properties ψ₂ cannot have:

• Sharper certification. ψ∞(U) ≤ ϵ certifies the weighted-max error max_n w^(T[1]−T[n]) ‖U[n] − U*[n]‖ ≤ s(θ_F) ϵ with s(θ_F) = Σ_{j<N} θ_F^j, θ_F = Λ_F/w — no factor N (thesis thm:equivalence_error_norm_local; the ℓ¹ constant is N·s(θ_F)). To certify a target radius r, set ϵ = r / s(θ_F).
• Prefix-decomposable acceptance. ψ∞ ≤ ϵ is the conjunction of per-chunk tests, decidable left to right with one boolean riding the messages a pipeline already sends. PipelinedMPIBackend therefore HONOURS ψ∞ — the one criterion it can enforce without re-serialising the pipeline — while ψ₁/ψ₂ remain bulk-synchronous only.

Like ψ₂, evaluated over the stored chunk STARTS (the terminal boundary is omitted — see the convention footnote in the docs) and clamping w ≥ 1. In the bulk-synchronous loop the n = 1 term is identically zero (F[1] mirrors U[1]).