Multigrid-in-Time Method

• Multigrid-in-Time methods are advanced algorithms that recast sequential time discretization as a hierarchical space-time system to enable parallel and efficient solution of temporal problems.
• They utilize temporal coarsening, various relaxation schemes (F-, C-, and FCF-relaxation), and transfer operators to correct low-frequency errors and accelerate convergence.
• Applications span discontinuous Galerkin discretizations, MGRIT frameworks, and preconditioning for KKT systems, demonstrating significant speedups and scalability in large-scale simulations.

A multigrid-in-time method is a multilevel algorithm that applies multigrid principles in the temporal dimension of a time-dependent problem. After time discretization, sequential stepping can be recast as a block lower bidiagonal or more general space-time system, and the method replaces a purely serial forward solve by relaxation on fine time points and coarse-grid correction on successively coarser temporal grids. In the literature, the term covers classical time multigrid for discontinuous Galerkin discretizations, multigrid-reduction-in-time (MGRIT), PFASST interpreted as a two-level multigrid method on composite collocation systems, waveform-relaxation variants for fractional evolution equations, and multigrid-in-time preconditioners for coupled optimality systems (Gander et al., 2014, Bolten et al., 2016, Hessenthaler et al., 2018, Gaspar et al., 2016, Vuchkov et al., 2024).

1. Foundational formulations

A standard starting point is the one-step evolution law

$$\mathbf{u}_{n+1}=\Phi_n \mathbf{u}_n+\mathbf{g}_{n+1},$$

or, for time-independent stepping,

$$A\mathbf{u}=\mathbf{g}, \qquad A= \begin{bmatrix} I \ -\Phi & I \ & \ddots & \ddots \end{bmatrix}.$$

In this form, conventional time marching is simply a block forward solve. Multigrid-in-time methods instead treat the entire temporal history as an algebraic object and build a hierarchy of coarser time grids on which low-frequency temporal error can be corrected (Hessenthaler et al., 2018).

The earliest formulations in the provided corpus already display the modern structure. For discontinuous Galerkin discretizations in time, the global system is block lower bidiagonal, the smoother is a damped block Jacobi iteration on time-slab blocks, and a two-grid cycle is analyzed by local Fourier analysis. In that setting, the method is designed as the time component of a space-time multigrid solver for parabolic problems and is reported to have excellent strong and weak scaling properties (Gander et al., 2014).

A distinct but closely related formulation appears in PFASST. Under linear-autonomous, two-level assumptions, PFASST can be written rigorously as a multigrid-in-time method applied to a composite collocation problem over the whole time interval. In that interpretation, spectral deferred correction acts as the smoother, the fine level corresponds to a parallel block Jacobi iteration across time subintervals, and the coarse level corresponds to a serial block Gauss-Seidel-like iteration with FAS coarse-grid correction (Bolten et al., 2016).

Within current practice, MGRIT is the most common nonintrusive realization. It works directly with the time-stepper as a black box, constructs a hierarchy of temporal grids, and solves the all-at-once space-time system iteratively rather than by sequential stepping. This nonintrusive character is emphasized in applications such as isogeometric analysis, where MGRIT is coupled to an external spatial solver without changing the underlying spatial discretization (Tielen et al., 2021).

2. Core algorithmic structure

The defining multigrid ingredients are temporal coarsening, relaxation, transfer operators, and recursive coarse-grid correction. With coarsening factor $m$, every $m$-th time point is a coarse-grid point or $C$-point, and the remaining time points are fine-only or $F$-points. The ideal coarse operator is the exact coarse propagator $\Phi^m$, but practical methods replace it by a cheaper approximation $\Psi$ or $\Phi_\Delta$ built by rediscretization or by a specialized coarse-grid construction (Friedhoff et al., 2019, Sterck et al., 2022).

The basic relaxation schemes are standardized. F-relaxation propagates from each $C$-point through the intervening $$A\mathbf{u}=\mathbf{g}, \qquad A= \begin{bmatrix} I \ -\Phi & I \ & \ddots & \ddots \end{bmatrix}.$$0-points; C-relaxation updates the $$A\mathbf{u}=\mathbf{g}, \qquad A= \begin{bmatrix} I \ -\Phi & I \ & \ddots & \ddots \end{bmatrix}.$$1-points from the preceding $$A\mathbf{u}=\mathbf{g}, \qquad A= \begin{bmatrix} I \ -\Phi & I \ & \ddots & \ddots \end{bmatrix}.$$2-points; and FCF-relaxation applies F-relaxation, then C-relaxation, then F-relaxation again. In the eddy-current formulation of MGRIT, restriction is injection to $$A\mathbf{u}=\mathbf{g}, \qquad A= \begin{bmatrix} I \ -\Phi & I \ & \ddots & \ddots \end{bmatrix}.$$3-points, and ideal interpolation is defined as injection followed by an F-relaxation. Multilevel $$A\mathbf{u}=\mathbf{g}, \qquad A= \begin{bmatrix} I \ -\Phi & I \ & \ddots & \ddots \end{bmatrix}.$$4- and $$A\mathbf{u}=\mathbf{g}, \qquad A= \begin{bmatrix} I \ -\Phi & I \ & \ddots & \ddots \end{bmatrix}.$$5-cycles are obtained by recursively applying the two-level correction process on deeper temporal hierarchies (Friedhoff et al., 2019).

For nonlinear problems, multigrid-in-time commonly uses FAS. In nonlinear hyperbolic conservation laws, the coarse-grid right-hand side is corrected by the nonlinear residual, and interpolation is implemented as injection followed by fine-point propagation with the fine time-stepper. In chaotic initial value problems, the same FAS framework is modified further by adding a nonlinear $$A\mathbf{u}=\mathbf{g}, \qquad A= \begin{bmatrix} I \ -\Phi & I \ & \ddots & \ddots \end{bmatrix}.$$6-correction that matches fine and coarse tangent dynamics more closely (Danieli et al., 2021, Vargas et al., 2022).

The same structural ideas appear outside classical initial-value problems. For KKT systems arising in time-dependent optimal control, virtual state variables and continuity constraints reorder the augmented system into a block-tridiagonal form in time. A block-Jacobi iteration then serves as a trivially parallel smoother, while geometric multigrid-style restriction and prolongation act along the time dimension and a coarse-grid GMRES iteration preconditioned by symmetric block Gauss-Seidel supplies the coarse solve (Vuchkov et al., 2024).

3. Convergence theory

Convergence theory for multigrid-in-time is strongly spectral. A central assumption in much of the MGRIT literature is simultaneous diagonalization of the time-stepping operators, so that bounds can be reduced to mode-wise scalar expressions involving eigenvalues of fine and coarse propagators. On this basis, multilevel frameworks derive residual and error propagation operators for $$A\mathbf{u}=\mathbf{g}, \qquad A= \begin{bmatrix} I \ -\Phi & I \ & \ddots & \ddots \end{bmatrix}.$$7- and $$A\mathbf{u}=\mathbf{g}, \qquad A= \begin{bmatrix} I \ -\Phi & I \ & \ddots & \ddots \end{bmatrix}.$$8-cycles, for F-relaxation and FCF-relaxation, and for varying numbers of levels (Hessenthaler et al., 2018).

For DG-in-time discretizations, the smoother can be analyzed almost completely. The damping behavior is governed by the DG/RADAU IA stability function $$A\mathbf{u}=\mathbf{g}, \qquad A= \begin{bmatrix} I \ -\Phi & I \ & \ddots & \ddots \end{bmatrix}.$$9, asymptotically optimal relaxation parameters can be derived, and for the lowest-order case $m$0 the two-grid convergence factor is

$m$1

which is uniformly bounded by $m$2 and improves as $m$3 grows (Gander et al., 2014).

The multilevel MGRIT theory generalizes this picture. A priori upper bounds are available for residual and error convergence, and the theory captures several empirical observations: two-level Parareal and multilevel MGRIT with F-relaxation do not yield scalable algorithms; stronger relaxation is beneficial; diffusion problems can achieve nearly constant convergence factors; and hyperbolic problems deteriorate with increasing numbers of levels. The same framework also indicates that L-stable Runge-Kutta schemes are more amendable to multilevel parallel-in-time integration with MGRIT than A-stable Runge-Kutta schemes (Hessenthaler et al., 2018).

A notable refinement appears for the unsteady fractional Laplacian. There, the fine and coarse propagators are time-dependent because the analysis allows nonuniform time stepping, and the paper introduces a temporal eigenvalue approximation property (TEAP) that measures how accurately coarse-grid eigenvalues approximate products of fine-grid eigenvalues. This removes the previous unitary diagonalization assumption and yields a generalized two-level convergence theory with a sharp upper bound that is validated numerically (Yue et al., 2019).

Weighted relaxation extends the same line of analysis. By weighting the C-relaxation step, one obtains new two-level convergence bounds involving the fine-grid eigenvalues $m$4, coarse-grid eigenvalues $m$5, and the relaxation weight $m$6. The numerical evidence shows that non-unitary weights consistently improve convergence, with iteration savings of about $m$7-$m$8 in most cases, and that under-relaxation can restore convergence for some A-stable integration schemes for which unweighted relaxation is not convergent (Sugiyama et al., 2021).

4. Specialized variants

Multigrid-in-time is not a single algorithmic template. One specialized branch is multigrid waveform relaxation for the time-fractional heat equation. Because Caputo derivatives produce dense lower-triangular Toeplitz-like time coupling, the method updates the entire time waveform at each spatial point, coarsens only in space, and exploits FFT-based Toeplitz operations. The resulting computational cost is $m$9 with storage $m$0, and the analysis is performed by semi-algebraic mode analysis rather than by standard local Fourier analysis (Gaspar et al., 2016).

Another variant is time-periodic MGRIT for periodic steady states. Here the fine-grid step function is modified so that the solution at $m$1 is sent back and used to update the initial value at $m$2. The algorithm therefore solves the fixed-point condition $m$3 directly over one period instead of simulating many full cycles to reach periodicity. The reported results show that the periodic steady state obtained is the same as that from sequential time-stepping (Hessenthaler et al., 2021).

AT-MGRIT replaces the global coarsest time grid by truncated local coarse grids, each covering only a short temporal subinterval. These local coarse problems can be solved completely independently, which reduces the sequential part of the algorithm. The two-level theory shows that truncation introduces an additional term in the convergence bound, but the numerical results indicate that this perturbation is typically mild and temporary relative to the reduction in coarse-grid serial work (Hahne et al., 2021).

Chaotic dynamics require more drastic modification. Standard MGRIT performs poorly because coarse-grid errors are amplified exponentially along unstable directions and naive time coarsening changes the Lyapunov spectrum. A modified algorithm therefore combines a $m$4-method coarse propagator with a nonlinear $m$5-correction that matches the linearized fine and coarse dynamics; for PDEs, a low-rank variant targets only the unstable subspace and simultaneously produces estimates of unstable Lyapunov vectors (Vargas et al., 2022).

Multigrid-in-time also extends naturally to coupled optimality systems. For augmented KKT systems in time-dependent optimization, the solver keeps the full primal-dual structure intact, introduces virtual state variables to relax continuity across time intervals, and builds a geometric multigrid method in the time dimension with a block-Jacobi smoother and FGMRES as the outer solver because the coarse solve is inexact (Vuchkov et al., 2024).

5. Representative applications

Electromagnetic simulation provides a clear example of the standard MGRIT workflow. For a voltage-driven eddy current problem on a 2D $m$6-$m$7 cross section of a coaxial cable, spatial discretization with edge shape functions produces a nonlinear index-1 DAE, backward Euler yields a nonlinear system $m$8 at each time step, and MGRIT solves the resulting all-at-once space-time system. On a fine uniform grid with $m$9 intervals on $C$0, 2269 spatial degrees of freedom, and a coarsest grid containing 4 time points, all tested $C$1- and $C$2-cycle variants converged in 7 to 9 iterations to residual tolerance $C$3, and the best speedup over sequential time stepping was about $C$4 (Friedhoff et al., 2019).

Fractional diffusion has motivated several distinct multigrid-in-time designs. For the unsteady fractional Laplacian, the fully discrete finite element formulation leads to time-dependent propagators $C$5, spatial solves are handled by the Chen–Nochetto–Otárola–Salgado multigrid with line smoother, and MGRIT iteration counts remain around roughly 11–14 over a range of fractional orders and time steps. For the 2-D space-fractional diffusion equation, a fully finite element scheme coupled to nonintrusive MGRIT achieves a best reported speedup of about $C$6, compared with about $C$7 for parareal, while the observed convergence factors are very close to the two-level theoretical bounds (Yue et al., 2019, Yue et al., 2018).

In isogeometric analysis, MGRIT has been combined with a $C$8-multigrid spatial solver for transient diffusion. The time-parallel iteration counts are typically about 5–10 and remain independent of mesh width $C$9, polynomial degree $F$0, and number of time steps $F$1. Replacing conjugate gradient by $F$2-multigrid greatly reduces the $F$3-dependence of the runtime, and strong scaling was demonstrated up to 2048 cores, with reported speed-up at 128 processors around $F$4 to $F$5 depending on $F$6 (Tielen et al., 2021).

Time-periodic and coupled multiphysics problems have also been treated successfully. In a simplified 2D Stokes model of the ascending aorta, time-periodic MGRIT over one cycle converged to the same periodic solution as sequential time-stepping over multiple cycles, with up to about $F$7 speedup on TheoSim with 256 processors for F-relaxation. In a 2D linear fluid-structure interaction benchmark, the method reduced the runtime from $F$8 for sequential time-stepping over 10 cycles to $F$9 on 32 processors, a speedup of $\Phi^m$0, and it was also applied to a nonlinear 3D FSI problem and to patient-specific cardiac flow (Hessenthaler et al., 2021).

For chaotic systems, the modified MGRIT algorithm extends the application range substantially. In the Lorenz system with $\Phi^m$1 and 8192 points, the combined $\Phi^m$2 method converged in 4 iterations where standard $\Phi^m$3 stalled or showed poor convergence. For the chaotic Kuramoto–Sivashinsky equation, weak scaling reached about $\Phi^m$4 speedup over sequential time stepping for the largest tested problem with rank-9 $\Phi^m$5-correction and $\Phi^m$6 coarse grids; strong scaling reached about $\Phi^m$7 speedup for $\Phi^m$8, and a longer run with $\Phi^m$9 achieved about $\Psi$0 speedup (Vargas et al., 2022).

Beyond initial-value problems, multigrid-in-time has been used as a preconditioner for KKT systems in viscous Burgers’ and van der Pol optimal control, where iteration counts remain stable under time refinement and the method is algorithmically scalable, and as the time component of a combined high-dimensional parabolic solver based on XBraid, sparse-grid combination, and domain decomposition. In the latter setting, the local-in-time combination variant was reported on a 2D linear oscillator with average runtime $\Psi$1 s on 1760 processes, versus $\Psi$2 s for full-grid sequential stepping and $\Psi$3 s for full-grid XBraid (Vuchkov et al., 2024, Griebel et al., 26 Sep 2025).

6. Limitations, misconceptions, and scope

A common misconception is that temporal coarsening alone is uniformly effective. The literature instead shows a sharp problem dependence. Diffusion-dominated and strongly damped problems are the most favorable setting; in contrast, hyperbolic, advection-dominated, and chaotic problems are much more delicate because the coarse-grid operator must capture phase propagation, CFL effects, or unstable tangent dynamics with much higher fidelity (Hessenthaler et al., 2018, Sterck et al., 2022, Vargas et al., 2022).

Advection is the clearest negative case for standard rediscretization. For method-of-lines discretizations of linear advection, a convergence analysis proves that robust MGRIT convergence with respect to CFL number and coarsening factor is not possible if the coarse grid is obtained by direct rediscretization. The failure is traced to inadequate correction of smooth characteristic components, and fast convergence is recovered only after replacing the coarse operator by a semi-Lagrangian discretization with an implicitly treated truncation-error correction (Sterck et al., 2022). A closely related conclusion holds for nonlinear hyperbolic conservation laws, where the principal determining factor for convergence degradation is the coarse-level CFL limit; the reported evidence suggests that effective multigrid-in-time for explicit hyperbolic problems likely needs spatial coarsening together with temporal coarsening (Danieli et al., 2021).

Another practical limitation is the serial bottleneck of the coarse solve. In classical Parareal and MGRIT, the coarsest level spans the whole time interval and remains inherently serial. AT-MGRIT addresses this by using truncated local coarse grids, while standard MGRIT studies on eddy-current problems show that deeper $\Psi$4-cycles may reduce iteration counts slightly but can become slower at scale because of additional communication (Hahne et al., 2021, Friedhoff et al., 2019).

The term also should not be restricted to a single data structure or transfer strategy. Some methods solve block lower bidiagonal all-at-once systems on temporal hierarchies; some, such as waveform-relaxation methods for fractional diffusion, iterate on whole time histories while coarsening only in space; some act on collocation systems, as in PFASST; and some preserve a fully coupled primal-dual structure, as in KKT preconditioners (Gaspar et al., 2016, Bolten et al., 2016, Vuchkov et al., 2024).

Finally, multigrid-in-time is often inseparable from spatial solver design. The practical performance in isogeometric analysis depended strongly on replacing CG by $\Psi$5-multigrid inside each time-step solve, while the robust generalized Stokes multigrid developed for implicit time-dependent Stokes subproblems is explicitly presented as a suitable building block for multigrid-in-time and all-at-once methods because it is robust with respect to both the spatial grid size and the time-step parameter $\Psi$6 (Tielen et al., 2021, Takacs, 2015). This suggests that, in large-scale applications, temporal multigrid and spatial multigrid are most effective when designed as a coupled solver architecture rather than as independent layers.