XBraid: Parallel-in-Time Solver

• XBraid is an open-source library that applies multigrid reduction in time (MGRIT) to introduce parallelism into the temporal dimension of unsteady simulations.
• It uses a non-intrusive wrapper around user-supplied time propagators, enabling simultaneous state, adjoint, and optimization iterations without requiring reformulated discretizations.
• The framework supports discrete adjoint sensitivity analysis and integration into large-scale, multiscale solver architectures for enhanced computational performance.

XBraid is an open-source software library, and a mature, non-intrusive C library, for bringing parallelism to the time dimension of unsteady simulations and unsteady PDE solvers. Its core algorithm is multigrid reduction in time (MGRIT), implemented in FAS form, and its central design choice is to wrap an existing, user-supplied time propagator rather than require a reformulation of the underlying discretization. In this way, XBraid parallelizes not only in space but also in the time domain, computes unsteady solutions parallel in time, extends the same non-intrusive structure to discrete adjoint sensitivity analysis, and can embed primal and adjoint iterations into simultaneous optimization workflows such as the One-shot method (Günther et al., 2017, Günther et al., 2018, Griebel et al., 26 Sep 2025).

1. Space–time formulation and MGRIT structure

For a one-step method, XBraid assumes a user-supplied propagator

$u_{n+1} = \Phi(u_n,p), \qquad n = 0,1,\ldots,N-1,$

or, in the notation used for unsteady PDE solvers,

$u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$

The same progression can be written as a residual equation

$$R_{n+1}(u_{n+1},u_n,p) = u_{n+1} - \Phi(u_n,p) = 0,$$

or, for linear time stepping, as a block lower bidiagonal space–time system. XBraid does not advance this system strictly sequentially. Instead, it views the entire time grid as a single space–time problem and applies an iterative multigrid algorithm in time (Günther et al., 2017, Günther et al., 2018).

MGRIT constructs a hierarchy of temporal grids by coarsening the time grid by a factor $m$ or $k$, classifying time points into C-points and F-points. On each level, relaxation sweeps and coarse-grid correction are applied. F-relaxation updates F-point values by propagating forward from the left C-point using the fine-grid propagator; these intervals are independent and computed concurrently. C-relaxation updates C-points using a coarse operator that approximates the $m$-fold composition of the fine propagator. A common pattern is FCF-relaxation: one sweep over F-points, one over C-points, followed by another F sweep. Intergrid transfer uses injection for restriction and prolongation, with prolongation augmented by relaxation to fill F-points. Standard V- and F-cycles can be used (Günther et al., 2017).

In the nonlinear case, XBraid uses FAS so that coarse-grid equations remain consistent with the current fine-grid iterate. One MGRIT cycle can be written as a nonlinear fixed-point operator

$u_{k+1} = H(u_k,p),$

which converges to the sequential time-stepping solution on the finest grid. The same viewpoint is used in the optimization literature, where the fixed-point equation $u = H(u,\rho)$ becomes the constraint for simultaneous optimization (Günther et al., 2017, Günther et al., 2018).

XBraid’s relation to other parallel-in-time methods is explicit in the literature. With two levels, MGRIT is equivalent to Parareal; with more levels, it achieves better concurrency and communication scaling. PFASST also uses nonlinear FAS ideas, whereas XBraid emphasizes wrapping arbitrary one-step propagators with minimal code changes. Schwarz-type time-parallel preconditioners offer concurrency bounded by subdomain sizes; XBraid’s concurrency is set by the temporal coarsening factor, which can be as small as $2$ (Günther et al., 2017).

2. Non-intrusive software interface and execution model

The library’s non-intrusive character is defined operationally. XBraid requires a small wrapper layer around an existing time-stepping code rather than a rewrite of the solver. In the adjoint paper, this wrapper consists of an application structure App holding time-grid metadata, MPI communicators, parameters, and any running objective; a vector structure storing a solution snapshot at one time point; a my_Step routine that advances the state from one time point to the next; a my_Access routine that evaluates per-time-step outputs; vector utilities such as init, clone, sum, norm, and free; and MPI pack/unpack routines so that states can move across processors and grids (Günther et al., 2017).

The simultaneous optimization paper describes the same pattern in multigrid language. Users supply callbacks for state stepping, state creation/cloning and destruction, basic linear algebra actions, and communication hooks to move data between levels and time partitions. For nonlinear problems, a residual or “apply” routine consistent with FAS is also supplied. The functional iterations that solve an implicit Crank–Nicolson step can remain untouched inside $\Phi^i$; XBraid orchestrates multilevel relaxation, coarse-grid correction, and inter-level communication around them (Günther et al., 2018).

The 2025 combination-method paper adds implementation details that are important for large-scale use. XBraid is described there as supporting variable/adaptive time-step sizes, distributing time steps across MPI processes in blocks that include each C point and the subsequent F points, and increasing per-process memory by a factor $u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$0 relative to sequential stepping. By default, XBraid can store only C points and reconstruct F points by F-relaxation when needed; this was the mode used in that paper. A plausible implication is that XBraid’s non-intrusiveness is not limited to a narrow class of fixed-step integrators, but extends to broader temporal configurations so long as the propagator interface is respected (Griebel et al., 26 Sep 2025).

This interface design is the main reason XBraid is described as non-intrusive across the papers. The user exposes the existing stepper and object-management routines through XBraid’s interface; XBraid then reformulates the serial time march into a global iterative solve over a hierarchy of time grids (Günther et al., 2017, Günther et al., 2018).

3. Discrete adjoint formulation and parallel-in-time sensitivities

XBraid’s adjoint extension was introduced as a consistent discrete adjoint solver for the MGRIT iteration. For objectives of the form

$u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$1

the standard discrete adjoint recursion for a one-step method is

$u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$2

with gradient

$u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$3

In XBraid, however, the primal solver computes $u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$4 by iterating $u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$5, so the adjoint is formulated for the fixed-point equation:

$u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$6

The primal and adjoint can also be advanced piggyback:

$u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$7

The papers report that this simultaneous iteration converges at the same asymptotic rate, with the adjoint typically lagging slightly behind the primal (Günther et al., 2017, Günther et al., 2018).

The implementation mechanism is reverse-mode AD over an “action tape” of the primal XBraid iteration. The tape records the sequence of XBraid calls to user routines such as my_Step, my_Access, clone, sum, and MPI pack/unpack. Adjoint replay traverses this tape in reverse and applies differentiated versions of those actions. Schematically, my_Step_adjoint updates the previous-time adjoint and accumulates parameter-gradient contributions through Jacobian-transpose times vector products, while my_Access_adjoint adds the derivatives of the instantaneous objective with respect to state and parameters. Linear utility routines propagate adjoints according to their linear operations, and MPI pack/unpack reverse as send↔recv of adjoint buffers. The prototype implementation used reverse-mode AD with CoDiPack, but the interface also admits hand-coded adjoints or algorithmic Jacobian-vector products (Günther et al., 2017).

The validation case in the 2017 adjoint paper couples a far-wake advection-diffusion PDE to a near-wake nonlinear oscillator. The PDE is

$u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$8

with

$u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$9

and the oscillator is

$$R_{n+1}(u_{n+1},u_n,p) = u_{n+1} - \Phi(u_n,p) = 0,$$0

The parameters were $$R_{n+1}(u_{n+1},u_n,p) = u_{n+1} - \Phi(u_n,p) = 0,$$1, $$R_{n+1}(u_{n+1},u_n,p) = u_{n+1} - \Phi(u_n,p) = 0,$$2, and $$R_{n+1}(u_{n+1},u_n,p) = u_{n+1} - \Phi(u_n,p) = 0,$$3, with $$R_{n+1}(u_{n+1},u_n,p) = u_{n+1} - \Phi(u_n,p) = 0,$$4 as the design variable. Time discretization used implicit Crank–Nicolson, advection used second-order linear upwind, diffusion used central differences, and nonlinear solves at each time step used functional iterations wrapped in my_Step. The objective was

$$R_{n+1}(u_{n+1},u_n,p) = u_{n+1} - \Phi(u_n,p) = 0,$$5

For $$R_{n+1}(u_{n+1},u_n,p) = u_{n+1} - \Phi(u_n,p) = 0,$$6, $$R_{n+1}(u_{n+1},u_n,p) = u_{n+1} - \Phi(u_n,p) = 0,$$7, and $$R_{n+1}(u_{n+1},u_n,p) = u_{n+1} - \Phi(u_n,p) = 0,$$8, with temporal coarsening $$R_{n+1}(u_{n+1},u_n,p) = u_{n+1} - \Phi(u_n,p) = 0,$$9 and up to three levels, the adjoint solver’s gradient matched finite differences within roughly $m$0–$m$1, and strong scaling curves for the adjoint closely followed the primal (Günther et al., 2017).

The same study reported weak-scaling speedups over serial of approximately $m$2, $m$3, and $m$4 for the primal at $m$5, $m$6, and $m$7, and approximately $m$8, $m$9, and $k$0 for the adjoint at the same configurations. The adjoint runtimes were roughly $k$1 the primal, close to a $k$2 ratio also observed in serial runs; the extra approximately $k$3 overhead was attributed to AD-based recording/replay and adjoint orchestration in XBraid (Günther et al., 2017).

4. Embedding XBraid in simultaneous optimization

The 2018 optimization paper extends XBraid from simulation and sensitivity analysis to simultaneous optimization with unsteady PDEs. The method embeds the primal and adjoint XBraid iterations into the One-shot framework, which updates design variables after each state and adjoint update rather than waiting for full convergence of the PDE and adjoint solves. The optimization problem is written as

$k$4

and the coupled iteration is

$k$5

$k$6

$k$7

In the reported tests, a simple constant preconditioner $k$8 with $k$9 yielded robust descent (Günther et al., 2018).

The paper also introduces an augmented Lagrangian merit function,

$m$0

and states that if $m$1, where $m$2 is the MGRIT contractivity and $m$3 quantifies adjoint time lag, then $m$4 is an exact penalty and descent of $m$5 implies convergence to an optimal solution. The paper reports convergence even with $m$6 in the test problem, that is, when using the reduced gradient only (Günther et al., 2018).

The validation problem is again advection-dominated flow control based on the Van-der-Pol oscillator and the 1D advection–diffusion equation, but with a tracking-type cost functional and Tikhonov regularization:

$m$7

with $m$8 and $m$9 computed from $u_{k+1} = H(u_k,p),$0. The discretization used $u_{k+1} = H(u_k,p),$1, $u_{k+1} = H(u_k,p),$2, $u_{k+1} = H(u_k,p),$3, $u_{k+1} = H(u_k,p),$4 spatial points, $u_{k+1} = H(u_k,p),$5, implicit Crank–Nicolson in time, second-order linear upwind for advection, central differences for diffusion, MGRIT coarsening $u_{k+1} = H(u_k,p),$6, and three temporal levels. CoDiPack generated the derivatives required for adjoint computations (Günther et al., 2018).

At fixed design $u_{k+1} = H(u_k,p),$7, both primal and adjoint residuals dropped simultaneously, with the adjoint showing the expected time lag. Strong scaling showed that speedup over the time-serial forward/backward loops begins beyond approximately $u_{k+1} = H(u_k,p),$8 processors; at $u_{k+1} = H(u_k,p),$9 processors, speedups of approximately $u = H(u,\rho)$0 for the state and approximately $u = H(u,\rho)$1 for the adjoint were reported. In optimization, the One-shot method reached the stopping criterion $u = H(u,\rho)$2 in $u = H(u,\rho)$3 iterations, and the objective leveled off at the regularization scale. With $u = H(u,\rho)$4 time processors, the reported timings were $u = H(u,\rho)$5 for time-serial reduced-space, $u = H(u,\rho)$6 for time-parallel reduced-space, and $u = H(u,\rho)$7 for time-parallel One-shot, corresponding to a $u = H(u,\rho)$8 speedup of the time-parallel One-shot method over the classical time-serial reduced-space method. Optimization overhead relative to a pure time-serial simulation of $u = H(u,\rho)$9 was reduced from $2$0 to $2$1 (Günther et al., 2018).

These results suggest that XBraid is not only a solver wrapper, but also a fixed-point framework in which state, adjoint, and design variables can be evolved simultaneously when the application admits a discrete adjoint and a reduced-gradient update.

5. Convergence theory, Runge–Kutta choices, and practical tuning

A substantial part of XBraid-related practice is selecting time integrators, coarsening factors, and relaxation schemes that lead to robust MGRIT convergence. For the linear model problem

$2$2

with $2$3 symmetric positive definite or skew-symmetric, the 2019 analysis studies two-level Parareal and MGRIT with Runge–Kutta time integration. With fine-grid and coarse-grid propagators derived from RK stability functions,

$2$4

the two-level convergence bounds are expressed by

$2$5

For SPD $2$6, these are the main a priori predictors of the two-level convergence factor (Friedhoff et al., 2019).

The principal conclusion is that not all RK schemes are equal from the perspective of parallel-in-time. For SPD problems, if both fine and coarse RK methods are L-stable, then

$2$7

and the method attains $2$8- and $2$9-independent convergence across the full spectrum including very stiff modes. Backward Euler, SDIRK22, SDIRK33, and ESDIRK-32 are listed as L-stable examples. By contrast, A-stable but not L-stable schemes such as trapezoidal/Crank–Nicolson, implicit midpoint, Gauss–Legendre, and ESDIRK-33 can lose this property; two-level convergence is then restricted to particular ranges of $\Phi^i$0, and $\Phi^i$1 changes in $\Phi^i$2 or the coarsening factor can change the worst-case factor from $\Phi^i$3 to divergence (Friedhoff et al., 2019).

For XBraid configuration, the paper’s guidance is explicit. If both fine and coarse schemes are L-stable, then coarsening factors $\Phi^i$4–$\Phi^i$5 typically yield robust, $\Phi^i$6-independent two-level convergence; FCF still improves rates, but F-relaxation may suffice for small to moderate $\Phi^i$7. If an A-stable fine scheme must be used, the coarse scheme should be L-stable, FCF relaxation should be turned on, and larger $\Phi^i$8, often in $\Phi^i$9, is much more robust; odd $u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$00 may help for certain A-stable schemes. For skew-symmetric operators, $u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$01-independent convergence cannot be obtained in general, and the recommendation is backward Euler with small $u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$02 and limited coarse levels (Friedhoff et al., 2019).

The same study connects these results to practical XBraid diagnostics. If per-iteration space–time residual reduction stagnates or $u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$03, the proposed remedies are to check coarse-grid stability, reduce $u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$04, switch the coarse integrator to an L-stable method, increase $u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$05 and switch on FCF, or reduce $u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$06. A Mathematica notebook for a priori two-grid analysis is available at https://github.com/XBraid/xbraid-convergence-est (Friedhoff et al., 2019).

Within the broader XBraid literature, these theoretical results clarify earlier empirical observations. The 2017 adjoint paper reports that larger coarsening led to instability in the nonlinear solver at coarse levels in an advection-dominated test, and the 2018 optimization paper notes that hyperbolic or advection-dominated problems often need careful parameter tuning and converge more slowly than parabolic ones. The 2019 analysis provides a mechanism for those observations in terms of RK stability, spectral location, and coarsening sensitivity (Günther et al., 2017, Günther et al., 2018, Friedhoff et al., 2019).

6. Application domains, large-scale couplings, and limitations

XBraid has been used in at least three distinct roles: as a time-parallel primal solver, as a time-parallel adjoint and optimization engine, and as the temporal component of larger multiscale or multidimensional solver stacks. In the 2025 parabolic-problem paper, XBraid is the time-parallel component of a three-scale solver consisting of MGRIT in time, the sparse grid combination method in space, and a space-filling-curve domain decomposition method for each spatial subproblem. The authors chose XBraid for its maturity, non-intrusive C implementation, and support for variable/adaptive time-step sizes (Griebel et al., 26 Sep 2025).

That paper formulates the semi-discrete problem for each spatial subproblem $u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$07 as

$u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$08

with backward Euler typically used,

$u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$09

Time coarsening uses a factor $u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$10, with $u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$11 in the experiments, and XBraid implements F-relaxation, C-relaxation, and combinations like CF and FCF, with prolongation given by injection followed by F-relaxation and restriction discarding F-points. Two integration strategies are described: $u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$12, which applies MGRIT on the sparse grid with one global time partition, and $u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$13, which applies MGRIT on each subproblem with subproblem-specific time partitions and recombination only at selected times (Griebel et al., 26 Sep 2025).

The reported applications include the heat equation up to $u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$14, chemical master equation models of 2D and 3D genetic toggle switches, and Fokker–Planck equations associated with linear SDEs. For the heat equation, the median MGRIT iterations per space-time subproblem were essentially independent of dimension $u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$15, global level $u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$16, and anisotropy. For $u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$17 with $u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$18, $u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$19, and $u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$20, the numbers of available processes were reported as $u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$21 for $u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$22, $u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$23 for $u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$24, $u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$25 for $u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$26, $u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$27 for $u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$28, and $u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$29 for $u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$30 (Griebel et al., 26 Sep 2025).

For the 2D oscillator Fokker–Planck example on $u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$31, with backward Euler, $u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$32, coarse factor $u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$33, and $u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$34, the reported average runtime of $u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$35 was $u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$36 on $u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$37 processes, compared with $u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$38 for the full-grid sequential method on $u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$39 processes and $u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$40 for full-grid XBraid on $u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$41 processes. The maximum error over time at the origin was $u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$42 for $u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$43, $u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$44 for the sequential full-grid method, and $u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$45 for full-grid XBraid. When the sequential full-grid method was forced to the same $u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$46 processes, runtime degraded to $u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$47, which the paper interprets as communication overhead dominating at high spatial concurrency; XBraid’s time concurrency enabled speedups beyond pure spatial parallelization (Griebel et al., 26 Sep 2025).

The literature also states XBraid’s limitations with comparable precision. Coarse-level instability for stiff or highly advective problems can limit the number of levels or require tailored coarse operators and relaxation tuning. Both serial time stepping and MGRIT are $u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$48 optimal, but MGRIT has a larger constant per iteration, so a sufficient number of time processors is required to beat the serial solve. The multigrid constant-factor overhead is described as often $u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$49–$u^i = \Phi^i(u^{i-1},\rho), \qquad i=1,\ldots,N.$50 compared to optimal sequential stepping. Adjoint-based optimization further requires consistent derivatives of the discrete integrator and objective, whether supplied by automatic differentiation, hand-coded derivatives, or Jacobian-vector products. In large multilevel settings, memory and communication trade-offs remain significant, although storing only C points mitigates the memory footprint (Günther et al., 2017, Günther et al., 2018, Griebel et al., 26 Sep 2025).

Across these studies, XBraid is positioned as a multilevel, non-intrusive infrastructure for temporal parallelism whose practical identity is defined by three recurring properties: it wraps existing one-step propagators with minimal code changes, it extends the same abstraction to adjoints and gradients, and it can serve as the time-parallel backbone inside more elaborate solver architectures. The reported evidence spans unsteady flow, adjoint sensitivity analysis, simultaneous optimization, and parabolic problems with up to six space dimensions (Günther et al., 2017, Günther et al., 2018, Griebel et al., 26 Sep 2025).