Timesteppers API

Time integration schemes for evolving PDEs in time.

Tarang's timesteppers are IMEX schemes: the linear operator L (everything on the left-hand side of an equation) is advanced implicitly, and the nonlinear right-hand side is advanced explicitly.

TimeStepper Types

Abstract Type

All timesteppers inherit from the abstract type Tarang.TimeStepper. It is not exported — write it qualified if you need to refer to it (e.g. in a type annotation):

using Tarang

supertype(RK222)                  # Tarang.TimeStepper
RK222() isa Tarang.TimeStepper    # true

The concrete schemes below are exported, unless the section says otherwise.

IMEX Runge-Kutta

RK111

1st-order IMEX Runge-Kutta.

RK111()

Properties:

• Order: 1
• Stages: 2 rows (explicit first stage), one implicit solve per step — the Ascher–Ruuth–Spiteri (1,1,1) RK111 form
• Implicit part: Backward Euler for linear terms
• Explicit part: Forward Euler for nonlinear terms, evaluated at the OLD time inside the implicit solve, so both tableaux are stiffly accurate and boundary conditions hold exactly after the step
• Workspace field sets: 3 (stage state, one RHS per stage)

RK222

2nd-order IMEX Runge-Kutta (Ascher, Ruuth & Spiteri 1997, ARS(2,2,2)) with three stored stages: an explicit first stage and two diagonally implicit stages.

RK222()

Properties:

• Order: 2
• Stages: 3 — an ESDIRK form whose first stage is explicit, so only 2 of the 3 stages require an implicit solve. (RK222().stages == 3; the explicit and implicit tableaux share the abscissae c = [0, γ, 1], which is what makes the IMEX pair 2nd order.)
• Implicit part: L-stable ESDIRK with γ = 1 - 1/sqrt(2) on stages 2 and 3
• Explicit part: matched three-row explicit tableau
• Workspace field sets: 4 (one stage state plus one retained RHS per stage)

Recommended for: General purpose problems.

RK443

3rd-order IMEX Runge-Kutta using the Ascher-Ruuth-Spiteri (4,4,3) pair. Five stored rows include the initial stage; four stages require implicit solves. The common stage times are c = [0, 1/2, 2/3, 1/2, 1].

RK443()

Properties:

• Order: 3
• Stages: 5 rows, first stage explicit, so 4 implicit solves per step
• Implicit part: L-stable ESDIRK for linear terms
• Explicit part: four explicit evaluations for nonlinear terms (the fifth row's F has zero weight everywhere and is not evaluated)
• Workspace field sets: 6 (one stage state plus one retained RHS per stage)

Recommended for: High accuracy requirements.

IMEX Multistep Methods

CNAB (Crank-Nicolson Adams-Bashforth)

CNAB1()  # 1st order
CNAB2()  # 2nd order

Treatment:

• Linear terms: Crank-Nicolson (implicit)
• Nonlinear terms: Adams-Bashforth (explicit)

Properties:

• CNAB1: Order 1, 1 history level
• CNAB2: Order 2, 2 history levels

SBDF (Semi-implicit BDF)

SBDF1()  # 1st order
SBDF2()  # 2nd order
SBDF3()  # 3rd order
SBDF4()  # 4th order

Treatment:

• Linear terms: BDF (implicit)
• Nonlinear terms: Extrapolation (explicit)

Properties: | Method | Order | History Levels | Stability | |––––|–––-|––––––––|–––––-| | SBDF1 | 1 | 1 | A-stable | | SBDF2 | 2 | 2 | A-stable | | SBDF3 | 3 | 3 | A(α)-stable | | SBDF4 | 4 | 4 | A(α)-stable |

Multistep schemes start from a single state, so the history has to be seeded. CNAB2 and SBDF2 bootstrap with their own 1st-order member (CNAB1 / SBDF1), so their first step is reduced order. SBDF3 / SBDF4 use the lower-order SBDF sequence on every architecture. Starting from a single state can therefore limit overall convergence to order two, despite the higher formal order once history is established. Only SBDF1/SBDF2 are A-stable: the higher-order members trade stability angle for formal order, so they suit smooth, accuracy-limited integration rather than the stiffest linear terms.

MCNAB2 / CNLF2

Additional two-step multistep schemes. Neither is exported — they must be qualified:

Tarang.MCNAB2()   # Three-level modified CNAB2
Tarang.CNLF2()    # Crank-Nicolson Leapfrog (CN implicit + centered leapfrog explicit)

MCNAB2() uses the three-level modified CNAB scheme. Its implicit weights are (9/16, 6/16, 1/16) at constant timestep, with variable-step corrections based on the current/previous timestep ratio. It has no theta parameter.

CNLF2 is 2nd order and 2-step. Its explicit leapfrog part is unstable for a negative real eigenvalue: for dt(u) = -a*u, a > 0, the parasitic amplification factor is -a*dt - sqrt(1 + (a*dt)^2), whose magnitude exceeds one at every positive timestep. Put dissipative terms on the implicit LHS when using CNLF2.

Exponential Time Differencing (ETD)

Exponential integrators that treat the linear term exactly via the matrix exponential, which is well suited to very stiff linear operators.

ETD_RK222()  # 2nd-order exponential Runge-Kutta
ETD_CNAB2()  # 2nd-order exponential Crank-Nicolson Adams-Bashforth
ETD_SBDF2()  # 2nd-order exponential semi-implicit BDF

Properties:

• Order: 2
• Linear term: exact exponential propagation (via φ-functions)
• Nonlinear term: explicit (RK / Adams-Bashforth / BDF extrapolation)
• Best for problems where the linear part dominates the stiffness.

Specialized IMEX-RK

Additive Runge-Kutta (ARK) schemes. RKSMR is exported; RKGFY is not, and must be qualified.

RKSMR

RKSMR()

The Spalart–Moser–Rogers scheme, the workhorse IMEX integrator of incompressible spectral DNS. Its nonlinear/explicit part is third order and its linear/implicit part (Crank-Nicolson on the viscous term) is second order — that asymmetric accuracy profile is the standard SMR one, not an implementation limitation.

It is stored in a four-stage additive-RK (ESDIRK) tableau — a trivial explicit first stage plus the three SMR substeps — so it shares the same generic IMEX runtime paths as RK222 and RK443.

RKGFY

Tarang.RKGFY()       # 2nd-order predictor-corrector, three stored stages

For pure implicit decay dt(u) + a*u = 0, RKGFY has amplification (1 - a*dt/2) / (1 + a*dt/2). It approaches -1 for very stiff decay, so the method is not L-stable and does not rapidly damp unresolved stiff modes.

RK443_IMEX

Tarang.RK443_IMEX is a const alias for RK443 — the same type, not a copy of its tableau. It was a separate type until it was found to diverge: the per-mode diagonal implicit solve is selected by ts isa Union{RK222, RK443, SBDF2}, so the twin lost that path and a single-GPU run with an implicit L that RK443() completed was refused under the other name.

Diagonal IMEX

RK222()  # 2nd order; automatic diagonal implementation on GPU
RK443()  # 3rd order
SBDF2()  # 2nd order multistep

For GPU RK222(), RK443(), and SBDF2(), diagonal implementations are internal; no GPU-specific timestepper constructors are exported. CPU and boundary-value problems retain the appropriate matrix or subproblem solve. An attached spectral linear operator also selects the internal diagonal path on serial CPU.

The implicit solve is diagonal in spectral space — the linear operator is applied per Fourier mode ((I + γ·dt·L̂)⁻¹ per wavenumber) rather than through a global matrix solve. Much cheaper for pure-Fourier / diagonalizable linear terms, and because the solve is an element-wise division there is no sparse factorization to move off-device on GPU. It requires the implicit operator to be diagonal in the chosen basis.

These schemes read the implicit operator from your equation's L and diagonalize it automatically, so ∂t(u) - nu*lap(u) = 0 is treated implicitly with no extra setup. You can still attach an explicit SpectralLinearOperator with set_spectral_linear_operator! on serial CPU or a single GPU to override that — an attached operator takes precedence and is applied to every field. MPI rejects attached operators with ArgumentError; write the linear term on the equation's LHS instead. Scalar division such as dt(u) - lap(u)/2 = 0 is supported by the diagonal parser on all three paths.

Serial attachment example:

using Tarang

coords = CartesianCoordinates("x")
dist   = Distributor(coords; dtype=Float64, device=CPU())
xb     = RealFourier(coords["x"]; size=16, bounds=(0.0, 2π))

u = ScalarField(dist, "u", (xb,), Float64)
x = local_grid(xb, dist, 1)
ensure_layout!(u, :g)
get_grid_data(u) .= cos.(2 .* x)          # a single mode, k = 2

problem = InitialValueProblem([u])
add_equation!(problem, "∂t(u) = 0")       # the linear term lives in L, NOT in the equation

solver = InitialValueSolver(problem, RK222(); dt=0.005)
L = SpectralLinearOperator(dist, (xb,), :laplacian; ν=0.5)   # L̂(k) = ν k²
set_spectral_linear_operator!(solver, L)

run!(solver; stop_iteration=200, progress=false)   # t = 1.0

ensure_layout!(u, :g)
maximum(abs, get_grid_data(u))   # 0.13533 ≈ exp(-ν k² t) = exp(-2) = 0.13534

SpectralLinearOperator(dist, bases, kind; ν, order) supports :laplacian (L̂ = ν k²), :hyperviscosity (L̂ = ν (k²)^order), :biharmonic (L̂ = ν k⁴), and :custom (pass coefficients yourself).

Usage with Solver

using Tarang

coords = CartesianCoordinates("x", "y")
dist   = Distributor(coords; dtype=Float64, device=CPU())
xb     = RealFourier(coords["x"]; size=16, bounds=(0.0, 2π))
yb     = RealFourier(coords["y"]; size=16, bounds=(0.0, 2π))
domain = Domain(dist, (xb, yb))

u = ScalarField(domain, "u")
problem = InitialValueProblem([u])
add_parameters!(problem, nu=0.1)
add_equation!(problem, "∂t(u) - nu*lap(u) = 0")
set!(u, (x, y) -> sin(x) * cos(y))

# Create solver with a timestepper — note the `()`
solver = InitialValueSolver(problem, RK222(); dt=1e-3)

# Drive it with run!
run!(solver; stop_iteration=50, progress=false)

ensure_layout!(u, :g)
maximum(abs, get_grid_data(u))   # 0.99004983, matching exp(-2νt) to 1.6e-11

Or drive the loop yourself:

solver = InitialValueSolver(problem, SBDF2(); dt=1e-3)
while solver.sim_time < 0.05
    step!(solver)
end

Timestepper State

TimestepperState

Tarang.TimestepperState holds everything a scheme needs between steps: the state history, the timestep history (for variable-dt multistep schemes), pre-allocated workspace fields, and any forcing. It is not exported and is managed automatically by the solver — you reach it as solver.timestepper_state and normally never touch it.

fieldnames(Tarang.TimestepperState)
# (:timestepper, :dt, :history, :dt_history, :stage, :timestepper_data,
#  :workspace_fields, :workspace_allocated, :forcing, :current_substep, :forcing_generated)

Stability

Explicit (advective) restriction

The explicit half of every IMEX scheme is limited by the advective CFL condition:

\[\Delta t < C \frac{\Delta x}{u_{max}}\]

The admissible factor C depends on the PDE, the spatial discretization and the scheme; Tarang does not define universal CFL constants for individual steppers. Rather than hard-coding C, use the CFL controller, which recomputes dt from the grid spacing and the current velocity field with a safety factor. It is constructed from the solver (never from a problem), and its velocities must be VectorFields:

v = VectorField(domain, "v")
set!(v.components[1], (x, y) -> 0.5)

prob = InitialValueProblem([u, v])
add_parameters!(prob, nu=0.1)
add_equation!(prob, "∂t(u) - nu*lap(u) = -v⋅∇(u)")
add_equation!(prob, "∂t(v) - nu*lap(v) = 0")
set!(u, (x, y) -> sin(x) * cos(y))

solver = InitialValueSolver(prob, RK222(); dt=1e-3)

cfl = CFL(solver; initial_dt=1e-3, cadence=5, safety=0.4, max_dt=0.01)
add_velocity!(cfl, v)

run!(solver; stop_iteration=10, cfl=cfl, progress=false)
solver.dt   # updated by the controller

Implicit (diffusive) restriction

Because the linear operator is advanced implicitly, the diffusive stability restriction (Δt ≲ Δx²/ν) is removed. Implicit linear treatment removes the explicit diffusive limit, but it does not remove the advective CFL limit above: only advective CFL constrains dt. That is the whole point of the IMEX splitting, and it is why all of Tarang's schemes are IMEX.

Method Selection Guide

Performance

Computational Cost

The RHS is evaluated once per stage; an implicit solve is needed on each stage whose implicit diagonal coefficient A_implicit[s,s] is nonzero (the ESDIRK first stage is free). Each IMEX-RK step then finishes with one mass-matrix solve for the update; the factorizations of (M + a·dt·L) and M are cached across steps and reused as long as dt is unchanged.

The multistep rows are the per-step cost once the history ring is full; the startup steps cost the same but run at reduced order.

For the internal diagonal implementations the stage counts are the same, but each "implicit solve" is an element-wise division by (1 + a·dt·L̂(k)) instead of a sparse factorization.

Memory Requirements

Workspace storage is preallocated per timestepper type and reused: these are scratch buffers, not fresh state copies allocated on every step. Workspace field sets are allocated per state field (each set is one full field):

These are the counts _workspace_count returns per scheme; the multistep and distributed paths keep their own recycled history buffers on top of them.

The pool does not grow on demand: get_workspace_field! returns a freshly allocated ScalarField for any index past the end, and the distributed diagonal-IMEX RK stepper asks for one per (stage, field) on every step. A scheme whose count is below what its path indexes therefore allocates a full field per stage per step — which is what RKGFY and RKSMR did before they were given a row above.

The multistep schemes additionally retain previous state and RHS levels in a history ring (2–5 levels, depending on order), which the RK schemes do not.

See Also

• Solvers: Using timesteppers
• User Guide: Timesteppers: Conceptual guide