Convergence analysis for parallel-in-time solution of hyperbolic systems

TL;DR

This paper analyzes the convergence of multigrid methods for parallel-in-time solutions of hyperbolic PDEs, introducing a combined analysis approach that improves understanding and potential stability of these methods.

Contribution

It develops a combined analysis framework integrating reduction analysis and semi-algebraic mode analysis for hyperbolic problems, enhancing convergence insights.

Findings

Combined analysis improves stability predictions for hyperbolic problems.

Numerical validation confirms the effectiveness of the combined analysis.

Insights apply to linear advection and elasticity models.

Abstract

Parallel-in-time algorithms have been successfully employed for reducing time-to-solution of a variety of partial differential equations, especially for diffusive (parabolic-type) equations. A major failing of parallel-in-time approaches to date, however, is that most methods show instabilities or poor convergence for hyperbolic problems. This paper focuses on the analysis of the convergence behavior of multigrid methods for the parallel-in-time solution of hyperbolic problems. Three analysis tools are considered that differ, in particular, in the treatment of the time dimension: (1) space-time local Fourier analysis, using a Fourier ansatz in space and time, (2) semi-algebraic mode analysis, coupling standard local Fourier analysis approaches in space with algebraic computation in time, and (3) a two-level reduction analysis, considering error propagation only on the coarse time grid. In this paper, we show how insights from reduction analysis can be used to improve feasibility of the semi-algebraic mode analysis, resulting in a tool that offers the best features of both analysis techniques. Following validating numerical results, we investigate what insights the combined analysis framework can offer for two model hyperbolic problems, the linear advection equation in one space dimension and linear elasticity in two space dimensions.

Figures

Equations

$$u_t+c{u}_x=0\text{ in }\Omega \times [0,T],$$

$$\left(1+\frac{c\Delta t}{\Delta x}\right)u_{j,i}-\frac{c\Delta t}{\Delta x}{u}_{j-1,i}=u_{j,i-1},i=1,\dots ,N_t,j=0,1,\dots ,N_x.$$

$$\rho u_{tt}-{\nabla }_X\cdot \sigma$$

$$\det F$$

$$\rho v_t$$

$$u_t$$

$$\nabla \cdot v$$

$$\rho v_i-\Delta t\mu {\nabla }^2{u}_i+\Delta t\nabla p_i=\rho v_{i-1},$$

$$u_i-\Delta t{v}_i=u_{i-1}.$$

$$⟨\rho v_i-\Delta t\mu {\nabla }^2{u}_i+\Delta t\nabla p_i,\chi ⟩$$

$$\iff \rho ⟨v_i,\chi ⟩+\Delta t^2\mu ⟨\nabla v_i,\nabla \chi ⟩-\Delta t⟨p_i,\nabla \cdot \chi ⟩$$

$$⟨u_i,\varphi ⟩-\Delta t⟨v_i,\varphi ⟩$$

$$-⟨\nabla \cdot v_i,\psi ⟩$$

$$\rho M{v}_i+\Delta t^2\mu K{v}_i+\Delta tB{p}_i$$

$$M{u}_i-\Delta tM{v}_i$$

$$\Delta t{B}^T{v}_i$$

$$\left[\begin{array}{ccc}\rho M+\Delta t^2\mu K& 0& \Delta tB\\ -\Delta tM& M& 0\\ \Delta t{B}^T& 0& 0\end{array}\right]\left[\begin{array}{c}{v}_i\\ u_i\\ p_i\end{array}\right]=\left[\begin{array}{c}\rho M{v}_{i-1}-\Delta t\mu K{u}_{i-1}\\ M{u}_{i-1}\\ 0\end{array}\right].$$

$$\left[\begin{array}{ccc}\rho M+\Delta t^2\mu K& 0& \Delta tB\\ -\Delta tM& M& 0\\ \Delta t{B}^T& 0& 0\end{array}\right]=\left[\begin{array}{ccc}\rho M+\Delta t^2\mu K& 0& 0\\ -\Delta tM& M& 0\\ \Delta t{B}^T& 0& -\Delta t^{2}S\end{array}\right]\left[\begin{array}{ccc}I& 0& \Delta t(\rho M+\Delta t^2\mu K{)}^{-1}B\\ 0& I& \Delta t^2(\rho M+\Delta t^2\mu K{)}^{-1}B\\ 0& 0& I\end{array}\right],$$

$$-\Delta t^{2}S{p}_i=-\Delta t{B}^T(\rho M+\Delta t^2\mu K{)}^{-1}\left(\rho M{v}_{i-1}-\Delta t\mu K{u}_{i-1}\right),$$

$$\left[\begin{array}{cc}\rho M+\Delta t^2\mu K& 0\\ -\Delta tM& M\end{array}\right]\left[\begin{array}{c}{v}_i\\ u_i\end{array}\right]$$

$$=\left[\begin{array}{cc}\rho \left(I-B{S}^{-1}{B}^T(\rho M+\Delta t^2\mu K{)}^{-1}\right)M& -\Delta t\mu \left(I-B{S}^{-1}{B}^T(\rho M+\Delta t^2\mu K{)}^{-1}\right)K\\ 0& M\end{array}\right]\left[\begin{array}{c}{v}_{i-1}\\ u_{i-1}\end{array}\right].$$

$$\left[\begin{array}{c}{v}_i\\ u_i\end{array}\right]$$

$$=\left[\begin{array}{cc}\rho P{\left(\rho M+\Delta t^2\mu K\right)}^{-1}M& -\Delta t\mu P{\left(\rho M+\Delta t^2\mu K\right)}^{-1}K\\ \rho \Delta tP{\left(\rho M+\Delta t^2\mu K\right)}^{-1}M& -(\Delta t{)}^2\mu P{\left(\rho M+\Delta t^2\mu K\right)}^{-1}K+I\end{array}\right]\left[\begin{array}{c}{v}_{i-1}\\ u_{i-1}\end{array}\right].$$

$$u^{\prime }(t)=f(t,u(t)),u(0)=g_0,t\in [0,T],$$

$$u_i=\Phi u_{i-1}+g_i,i=1,\dots ,N_t,$$

$$Au\equiv \left[\begin{array}{cccc}I& & & \\ -\Phi & I& & \\ & \ddots & \ddots & \\ & & -\Phi & I\end{array}\right]\left[\begin{array}{c}{u}_0\\ u_1\\ ⋮\\ u_{N_t}\end{array}\right]=\left[\begin{array}{c}{g}_0\\ g_1\\ ⋮\\ g_{N_t}\end{array}\right]\equiv g.$$

$$u_{\Delta ,j}={\Phi }_c{u}_{\Delta ,j-1}+g_{\Delta ,j},j=1,\dots ,N_T,$$

$$E^F=(I-P_{\Phi }{A}_c^{-1}{R}_{I}A)P_{\Phi }{R}_I=P_{\Phi }(I-A_c^{-1}{A}_S)R_I,$$

$$A_S=\left[\begin{array}{cccc}I& & & \\ -{\Phi }^m& I& & \\ & \ddots & \ddots & \\ & & -{\Phi }^m& I\end{array}\right].$$

$$E^{FCF}=(I-P_{\Phi }{A}_c^{-1}{R}_{I}A)P_{\Phi }(I-A_S)R_I=P_{\Phi }(I-A_c^{-1}{A}_S)(I-A_S)R_I.$$

$$E^{V(1,3)-FCF}=\left(I-P_{\Phi }\left(I_c-E^{(2,3)-FCF}\right)A_c^{-1}{R}_{I}A\right)P_{\Phi }(I-A_S)R_I,$$

$$E^{F(1,3)-FCF}=\left(I-P_{\Phi }\left(I_c-E^{(2,3)-FCF}{E}^{(2,3)-FCF}\right)A_c^{-1}{R}_{I}A\right)P_{\Phi }(I-A_S)R_I,$$

$$G_{\Delta x}=\left\{x_j=j\Delta x:j\in \mathbb{Z}\right\}.$$

$$\phi (\theta ,x_j)$$

$${\tilde{\Phi }}_{\theta }={\left[1+\frac{c\Delta t}{\Delta x}\left(1-e^{-ı\theta }\right)\right]}^{-1}.$$

Full text

††articletype: Article Type

1 Introduction

While spatial parallelism is a well-known tool in scientific computing, hardware trends and scaling limits have served to renew interest in algorithms that also allow “space-time” parallelism. These techniques consider the solution of time-dependent systems of partial differential equations (PDEs), and aim to compute the solution in an “all-at-once” manner, breaking the sequential nature of traditional time-stepping approaches. While this is not a new idea1, 2, recent years have seen significant effort devoted to space-time and time-parallel approaches for the solution of time-dependent PDEs3, 4, 5, 6, 7, 8, 9, 10. With such intense interest in the development of new schemes, there is a pressing need for complementary analysis tools, to provide understanding of the relative performance of related schemes and to inform the optimization of algorithmic parameters as schemes are adapted to new problems. The central aim of this paper is to compare and contrast three such analysis schemes for parallel-in-time algorithms of the Parareal11, 12 and multigrid-reduction-in-time3 (MGRIT) methodologies.

From the multigrid perspective, mode analysis is an attractive tool for analyzing convergence of these methods, using either eigenvectors or Fourier modes (and related techniques) to predict convergence rates. In the context of space-time discretizations, three approaches to mode analysis have been discussed in the recent literature: space-time local Fourier analysis13, 14, 4 (LFA), semi-algebraic mode analysis15, 9, 16 (SAMA), and reduction-based theory17, 8, 18, 19. A fourth approach, Fourier-Laplace mode analysis20, 21, 22, can be seen as an intermediate between space-time LFA and the other two approaches, and will be discussed only in that context. While not considered here, the recent analysis of Gander et al.23 extends the earlier result of Gander and Hairer 24 from the Parareal context to MGRIT. These approaches bound the convergence of the algorithm based on Lipschitz continuity and other properties of the time propagators at the coarse and fine levels; while this offers insight into the general convergence properties of the algorithms, it is difficult to compare fairly to the mode analysis tools that are the focus here.

In broad terms, the advantages and disadvantages of the three mode analysis tools are as follows. LFA is a well-known and widely used tool for analysis of spatial multigrid methods25, 26, for which it allows quantitative predictions of two- and multi-level convergence. However, numerical experience13, 14 shows that it is not predictive for space-time methods in many situations, particularly for meshes with few points in time or discretizing “short” temporal domains. The extension of LFA to semi-algebraic mode analysis15 overcomes this limitation on space-time LFA; however, as implemented by Friedhoff and MacLachlan, the quantitative estimates of SAMA required the computation of Euclidean operator norms of matrices with size equal to the number of time steps multiplied by the dimension of the spatial LFA symbol. When done exactly, this can be prohibitively expensive for complicated problems over long time intervals. Two-level reduction analysis17, 8 overcomes this computational expense, under certain conditions, but is based on solving eigenvalue problems of size equal to the number of degrees of freedom in the spatial discretization. In the results below, we explore each of these methods for two specific hyperbolic model problems, linear advection in one dimension and incompressible linear elasticity in two dimensions, and compare and contrast their predictions in these settings.

Several important observations come from this comparison. First, we see that the bounds used to make the reduction analysis computationally tractable can be directly applied to SAMA. This results in a tool that combines the best of both approaches, using the flexibility of spatial LFA in settings where the spatial eigenproblem is intractable, but the ease of computing a bound for large numbers of time steps that comes from reduction analysis. Secondly, we expose the complications of assuming unitary diagonalizability, as is done in the reduction analysis, and detail how this can be avoided in the SAMA bound. Overall, this allows us to get more insight into MGRIT convergence, particularly as we vary both spatial and temporal problem sizes. Such insight, in turn, can enable the development of improved MGRIT algorithms for difficult problems, as has recently been done for high-order explicit discretizations of the linear advection equation 27. A final benefit is the ability to critically compare space-time LFA with the improved SAMA analysis, and gain better insight into problems and parameter regimes for which space-time LFA may be a feasible and sufficiently accurate option.

The remainder of this paper is organized as follows. Details of the model problems, linear advection in one spatial dimension and incompressible linear elasticity in two spatial dimensions, are given in Section 2. The Parareal and MGRIT algorithms are reviewed in Section 3. Section 4 describes the three mode analysis tools in detail, with discussion on how to combine aspects of SAMA and reduction analysis in Section 4.4. Numerical results comparing the analysis methods as they exist in the literature are given in Section 5.1. Sections 5.2 and 5.3 present numerical results giving new insight into convergence of the Parareal and MGRIT algorithms for linear advection and elasticity, respectively, based on SAMA improved by insights from reduction analysis. Concluding remarks are given in Section 6.

2 Model problems

In this section, we discuss the details of the two model problems for which the analysis tools will be compared in Section 5. For linear advection, we consider a simple upwinding-in-space and implicit-Euler-in-time discretization, corresponding to that used by De Sterck et al. 10 For linear elasticity, we use a mixed-finite-element-in-space and implicit-Euler-in-time discretization corresponding to that used by Hessenthaler et al. 8 We note that “better” (primarily less diffusive) discretizations of both equations certainly exist, but (to our knowledge) successful application of Parareal and MGRIT to those discretizations either requires a modification of the algorithm (e. g., the PITA framework28) or is not well documented in the published literature, aside from some forthcoming results27. The analysis presented herein can be applied to such discretizations in the same manner as used here; we leave such details to future work, to reflect the (hopefully) continued success of developing effective parallel-in-time algorithms for a wider class of discretizations and problems.

2.1 Linear advection

We consider the advection of a scalar quantity, $u(x,t)$, subject to a known non-zero constant flow speed, $c$, in the domain $\Omega\times[0,T]$, where $\Omega$ is a finite interval, $[a,b]$, and $T$ denotes the final time. For example, this models advection of nonreactive particles by an incompressible fluid, where particle density, $u(x,t)$, depends only on advection of particles by the fluid. The governing equation is given by

[TABLE]

where we assume $c>0$ for the subsequent discussion. We prescribe the initial condition $u(x,0)=u_{0}(x)$ and impose the periodic boundary condition $u(a,t)=u(b,t)$ in all that follows. We discretize (1) on a uniform space-time grid, with $N_{x}$ spatial intervals of width ${\Delta x}=(b-a)/N_{x}$ and $N_{t}$ temporal intervals of time step $\Delta t=T/N_{t}$, using a first-order implicit upwind scheme:

[TABLE]

2.2 Linear elasticity

We consider the dynamic and linear elastic response of an incompressible solid structure in the domain $\Omega\times[0,T]$, where $\Omega$ is an open domain in $\mathbb{R}^{2}$, and $T$ denotes the final time. Denoting the current and reference position of a material point by $x$ and $X$, respectively, and the respective Eulerian and Lagrangian gradient operators by $\nabla$ and $\nabla_{X}$, we define the deformation gradient, $F$, by $F=\nabla_{X}x=I+\nabla_{X}u$, where $x(X,t)=X+u(X,t)$ defines the displacement, $u$, of the material point, $x$, in the current configuration at time $t$ with respect to its position in the reference configuration, $X$, and where $I$ denotes the identity matrix of corresponding size. Then, the governing equations are given by

[TABLE]

where $\rho$ denotes material density, and where $\sigma=\sigma(u,p)=\mu(F-I)-pI$ is the Cauchy stress tensor for an incompressible linear elastic material with stiffness parameter, $\mu$, and hydrostatic pressure, $p$. We prescribe $u$ and $u_{t}$ at $t=0$, $u(X,0)=0$ and $u_{t}(X,0)=\widehat{v}_{0}$, and impose homogeneous Dirichlet boundary conditions, $u(X,t)=0$ for $X\in\Gamma^{D}$ and $t\in[0,T]$, where $\Gamma^{D}$ denotes the Dirichlet boundary of the domain $\Omega$. Using $\sigma=\mu\nabla_{X}u-pI$, neglecting higher-order effects of the deforming domain, and transforming Equation (3) to a system of first-order equations, we obtain

[TABLE]

for displacement, $u$, velocity, $v$, and hydrostatic pressure, $p$.

We discretize Equations (5)-(7) on an equidistant time grid consisting of $N_{t}\in\mathbb{N}$ time intervals using a time-step size $\Delta t=T/N_{t}$. Motivated by existing results8, we use implicit Euler for the time discretization. Denoting displacement, velocity, and pressure at time $t_{i}=i\Delta t$, $i=0,\ldots,N_{t}$, by $u_{i},v_{i}$, and $p_{i}$, respectively, for $i=1,\ldots,N_{t}$, Equations (5) and (6) are discretized as

[TABLE]

and

[TABLE]

The weak form of Equations (8), (9), and (7) is found by multiplying by test functions, $\chi\in(H^{1}(\Omega))^{2}$, $\phi\in(L_{2}(\Omega))^{2}$, and $\psi\in L^{2}(\Omega)$, respectively, giving

[TABLE]

We discretize the spatial domain, $\Omega$, on a uniform quadrilateral grid with mesh size ${\Delta x}$ using Taylor-Hood $(Q2-Q1)$ elements 29, 30 for velocity, $v$, and pressure, $p$, and $Q2$ elements for displacement, $u$, ensuring uniqueness of the solution. Denoting the mass and stiffness matrices of the discretized vector Laplacian by $M$ and $K$, respectively, and the (negative) discrete divergence operator by $B^{T}$, Equations (10)-(12) are discretized to

[TABLE]

which are equivalent to the following linear system of equations

[TABLE]

It is tempting to simply remove the rows and columns corresponding to the pressure variable from this system, since $p_{i-1}$ does not appear in the equations at time step $t_{i}$. Indeed, this approach was taken in the analysis by Hessenthaler et al.8. However, to do so ignores the important role that $p_{i}$ plays as a Lagrange multiplier, enforcing the incompressiblity constraint for the solution at time step $t_{i}$ (but was properly accounted for in simulations). Instead, we eliminate the contribution from this block by considering the block factorization of the matrix in (16) and the resulting Schur complement. Factoring this matrix gives

[TABLE]

for $S=B^{T}(\rho M+\Delta t^{2}\mu K)^{-1}B$. Thus, $p_{i}$ satisfies

[TABLE]

and can be directly eliminated from the equation for $v_{i}$ by subtracting $\Delta tBp_{i}$ from the right-hand side of (13). Noting that both velocity and displacement are two-dimensional vector fields of $Q2$ degrees of freedom, this yields the reduced system with four scalar functions,

[TABLE]

It is this form of the propagator that we analyse below. We note that although the only nonzero operator acting on the pressure acts on $p_{i}$, its effect is to change the dependence of $v_{i}$ and $u_{i}$ on $v_{i-1}$ and $u_{i-1}$. This can easily be seen to take the form of an orthogonal projection operator acting on the data from the previous time-step; it is also easy to check that this formulation guarantees that $B^{T}v_{i}=0$, as required by the incompressibility constraint.

As is common in MGRIT, we primarily analyse this propagator in “$\Phi$-form”, i. e., in the form $[v_{i}~{}u_{i}]^{T}=\Phi[v_{i-1}~{}u_{i-1}]^{T}$, writing the projection operator, $P=I-(\rho M+\Delta t^{2}\mu K)^{-1}BS^{-1}B^{T}$, to give

[TABLE]

3 Parallel-in-time methods

3.1 Parareal

The Parareal algorithm 12 is a parallel method for solving systems of ordinary differential equations of the form

[TABLE]

arising, for example, when solving a system of PDEs using a method-of-lines approximation to discretize the spatial domain. Parareal can be interpreted in a variety of ways, including as multiple shooting, domain decomposition, and multigrid methods 11, 3. Here, we describe Parareal as a two-level time-multigrid scheme. For ease of presentation, we only describe the method in a linear setting, i. e., in the case that $f$ is a linear function of $u(t)$; the full approximation storage (FAS) approach 31 can be applied in the same manner to accomodate nonlinear problems.

Parareal combines time stepping on the discretized temporal domain, the fine grid, with time stepping on a coarser temporal grid that uses a larger time step. More precisely, consider a fine temporal grid with points $t_{i}=i\Delta t$, $i=0,1,\ldots,N_{t}$, with constant time step $\Delta t=T/N_{t}>0$, and let $u_{i}$ be an approximation to $u(t_{i})$ for $i=1,\ldots,N_{t}$, with $u_{0}=u(0)$. Further, consider a one-step time integration method,

[TABLE]

with time-stepping operator, $\Phi$, that takes a solution at time $t_{i-1}$ to that at time $t_{i}$, along with a time-dependent forcing term, $g_{i}$. (Note that both the assumption of constant time step and of a time-independent single-step time-stepping operator are for notational convenience only, and can easily be relaxed.) The discrete approximation to the solution of (18) can be represented as a forward solve of the block-structured linear system

[TABLE]

Note that, in the time dimension, this forward solve is completely sequential.

Parareal enables parallelism in the solution process by introducing a coarse temporal grid, or (using multigrid terminology) a set of C-points, derived from the original (fine) temporal grid by considering only every $m$-th temporal point, where $m>1$ is an integer called the coarsening factor. Thus, the coarse temporal grid consists of $N_{T}=N_{t}/m$ intervals, with points $T_{j}=j\Delta T$, $j=0,1,\ldots,N_{T}$, with spacing $\Delta T=m\Delta t$; the remaining temporal points define the set of F-points, as visualized in Figure 1.

The coarse-grid problem, $A_{c}u_{\Delta}=g_{\Delta}$, is defined by considering a one-step method with time-stepping operator, $\Phi_{c}$, using time step $\Delta T$, given by

[TABLE]

where $u_{\Delta,j}$ denotes an approximation to $u(T_{j})$ for $j=1,\ldots,N_{T}$, and $u_{\Delta,0}=u(0)$.

Rather than sequentially solving for each $u_{i}$, Parareal alternates between a “relaxation scheme” on the F-points and a sequential solve on the C-points. The first of these two processes, known as F-relaxation, updates the unknowns at F-points by propagating the values from each C-point, $T_{j}$, across a coarse-scale time interval, $(T_{j},T_{j+1})$, $j=0,1,\ldots,N_{T}-1$ using the fine-grid time-stepper, $\Phi$. Note that within each coarse-scale time interval, these updates are sequential, but that there is no dependency across coarse time intervals, enabling parallelism in the relaxation process. After F-relaxation, the residual is evaluated only at the C-points and “injected” to the coarse temporal grid. After a sequential solve of the coarse-grid equations (for which typical choices in Parareal use high-order integration schemes over longer time intervals), a correction is interpolated to the fine grid, using “ideal interpolation” (which is “ideal” as it corresponds to using the Schur complement on the coarse grid). This interpolation operator is defined by taking the corrected approximate solution at at each C-point, $T_{j}$, and time-stepping across the coarse-scale time interval, $(T_{j},T_{j+1})$, $j=0,1,\ldots,N_{T}-1$, again using the fine-grid time-stepper, $\Phi$. The injection and ideal interpolation operators are block operators of size $(N_{T}+1)\times(N_{t}+1)$ or $(N_{t}+1)\times(N_{T}+1)$, respectively, given by

R_{I}=\setcounter{MaxMatrixCols}{14}\begin{bmatrix}I&0&\cdots&0\\ &&&&I&0&\cdots&0\\ &&&&&&&&\ddots\\ &&&&&&&&&I&0&\cdots&0&\\ &&&&&&&&&&&&&I\end{bmatrix}\quad\text{and}\quad P_{\Phi}=\begin{bmatrix}Z^{(\Phi)}\\ &Z^{(\Phi)}\\ &&\ddots\\ &&&Z^{(\Phi)}\\ &&&&I\end{bmatrix}~{}\text{with }Z^{(\Phi)}=\begin{bmatrix}I\\ \Phi\\ \vdots\\ \Phi^{m-1}\end{bmatrix}.$$m blocks

As the error propagator of F-relaxation can be written as $P_{\Phi}R_{I}$, the Parareal algorithm may be represented by the two-level iteration matrix,

[TABLE]

where equivalence holds since $R_{I}AP_{\Phi}$ defines the Schur complement coarse-grid operator,

[TABLE]

3.2 Multigrid-reduction-in-time

The major sequential bottleneck for Parareal is the sequential solve of the coarse-grid system, $A_{c}u_{\Delta}=g_{\Delta}$. Directly applying multigrid principles and using Parareal recursively to solve this system is observed to often have significantly degraded convergence properties. This motivates the multigrid-reduction-in-time (MGRIT) algorithm 3, which is based on applying multigrid reduction techniques 32, 33 to the time dimension. Like Parareal, MGRIT uses injection and ideal interpolation for intergrid transfer operations, and it uses the same coarse-grid operator, $A_{c}$, as can be used in the two-level Parareal method. However, in order to obtain a scalable multilevel algorithm, MGRIT augments F-relaxation, using relaxation also at the C-points. This C-relaxation is defined analogously to F-relaxation, by leaving the values of F-points unchanged, and updating the unknowns at C-points using the time-stepper, $\Phi$, and values at neighboring F-points. Relaxation in MGRIT consisting of combined sweeps of F- and C-relaxation; typically, FCF-relaxation, F-relaxation, followed by C-relaxation, and again F-relaxation, is employed. This scheme can be shown to be equivalent to Richardson relaxation on the coarse time grid with the Schur complement coarse-grid operator, $A_{S}$, followed by ideal interpolation3.

Writing the error propagator of FCF-relaxation as $P_{\Phi}(I-A_{S})R_{I}$, the two-level MGRIT algorithm may be represented by the two-level iteration matrix,

[TABLE]

Instead of inverting $A_{c}$ directly as in the two-level MGRIT algorithm (21), in the three-level method, the system on the first coarse grid is approximated by a two-grid cycle with zero initial guess, that is, the term $A_{c}^{-1}$ in the two-level iteration matrix (21) is replaced by $(I_{c}-E^{(2,3)-FCF})A_{c}^{-1}$ (where $E^{(2,3)-FCF}$ is defined analogously to (21)) to obtain the three-level V-cycle iteration matrix,

[TABLE]

or, it is replaced by $(I_{c}-E^{(2,3)-FCF}E^{(2,3)-FCF})A_{c}^{-1}$ to obtain the three-level F-cycle (or, equivalently, W-cycle) iteration matrix

[TABLE]

respectively. Iteration matrices of three-level V- and F-cycle Parareal methods can be derived analogously.

While there are clear and important differences, both historical and in practice, between Parareal and MGRIT, in this paper we view them as being two methods within the same broader family. As such, we will consider in detail only the case where the coarse-grid operator, $A_{c}$, is given by rediscretization of the fine-grid operator with the coarse time-step. Since Parareal methods often make use of different coarse- and fine-propagators, we will instead refer to the methods where this condition is imposed as MGRIT with F- and FCF-relaxation, instead of as Parareal and MGRIT, respectively, but use the more distinctive terminology for statements that are true regardless of the choice of coarse propagator. An interesting property of Parareal and MGRIT is the following exactness property: Assuming exact arithmetic, one iteration of F-relaxation computes the exact solution at the first $m-1$ time steps, i. e., at all F-points in the first coarse-scale time interval. Furthermore, one iteration of FCF-relaxation computes the exact solution at the first $2m-1$ time steps. Therefore, MGRIT with F-relaxation solves for the exact solution in $N_{t}/m$ iterations, while MGRIT with FCF-relaxation solves for the exact solution in $N_{t}/(2m)$ iterations, corresponding to half the number of points on the coarse grid.

4 Mode analysis tools for parallel-in-time methods

While both Parareal and MGRIT can be analysed based on abstracted properties of the coarse and fine time integrators24, 23, these tools do not provide deep insight into the success or failure of algorithmic choices within the methodology. Instead, here we present three tools based on mode analysis principles applied to the iteration matrices, $E^{F}$ and $E^{FCF}$, given in Equations (20) and (21), respectively, and to the iteration matrices of the respective three-level methods.

4.1 Space-time local Fourier analysis

Local Fourier analysis is a classical tool for analysing convergence of multigrid methods31, 26, that has been applied to space-time problems with mixed results in terms of accuracy in predictions 13, 4. Here we present a brief review of the calculation of spatial Fourier symbols (as will be used in all of the mode analysis tools considered in this section), along with their coupling with local Fourier analysis in the time direction.

4.1.1 Spatial Fourier symbols

The local Fourier analysis of the time-stepping operator, $\Phi$, makes use of its Fourier representation, also called the Fourier symbol. For the scalar advection equation, this follows from standard analysis; however, LFA for systems of PDEs requires more advanced tools34, 35. Note that we do not consider any coarsening of the spatial domain here to limit the scope of the paper and, thus, we only need the spatial Fourier symbol of the time-stepping operator.

Standard LFA for a problem in one spatial dimension considers a scalar discrete problem posed on an infinite uniform grid with mesh size ${\Delta x}>0$,

[TABLE]

Extending a constant-coefficient discrete operator from a finite spatial domain, such as that on the left of Equation (2), to the infinite grid $G_{\Delta x}$ leads to a Toeplitz operator acting on the space of square-summable sequences ($\ell_{2}$) that is formally diagonalizable by standard Fourier modes,

[TABLE]

Noting that the operator $\Phi$ is the inverse of the lower bidiagonal matrix on the left of (2), we find that its Fourier symbol, denoted $\widetilde{\Phi}_{\theta}$, is the scalar

[TABLE]

Furthermore, the Fourier symbols for the time integrators $\Phi_{c}$ and $\Phi_{cc}$ on the first and second coarse grids can be obtained by replacing $\Delta t$ by $m\Delta t$ in (24) to account for factor-$m$ temporal coarsening. We note that for a spatial problem with periodic boundary conditions, the eigenvalues of $\Phi$ are given precisely by sampling $\widetilde{\Phi}_{\theta}$ at evenly spaced points in $(-\pi,\pi]$; for other boundary conditions, a good estimate of the discrete eigenvalues of $\Phi$ is given by the range of its Fourier symbol36.

For the elasticity operator in (17), we must generalize this analysis to consider a two-dimensional infinite uniform grid, along with adaptations to handle the block structure imposed by the mixed finite-element discretization. Accordingly, we consider the two-dimensional infinite uniform grid, decomposed as

[TABLE]

with

[TABLE]

This decomposition arises by considering the variation in basis functions for the standard $Q2$ approximation space, with nodal basis functions defined for mesh points at corners of ${\Delta x}\times{\Delta x}$ quadrilateral elements, $\mathbf{x}^{N}\in G_{{\Delta x},N}$, mesh points located on the $x$- and $y$-edges of elements, $\mathbf{x}^{XE}\in G_{{\Delta x},XE}$ or $\mathbf{x}^{YE}\in G_{{\Delta x},YE}$, and cell-centered mesh points, $\mathbf{x}^{C}\in G_{{\Delta x},C}$.

While the finite-element discretization of a scalar PDE on a structured two-dimensional infinite uniform grid using nodal $(Q1)$ finite elements leads to an operator that is block Toeplitz with Toeplitz blocks (BTTB) that has a scalar symbol, the same is not true for higher-order elements. However, by using the same basis on each element and permuting the operator into block form ordered by “type” of basis function, we readily obtain a block operator with blocks that are themselves BTTB operators and can be used to generate a matrix-based Fourier symbol37, 35. On $G_{\Delta x}$, this arises by considering a discrete operator, $L_{\Delta x}$, defined in terms of its stencil for each “type” of degree of freedom. For example, for the nodal degrees of freedom, we have

[TABLE]

with constant coefficients, $s_{\boldsymbol{\kappa}}\in\mathbb{C}$, and $\boldsymbol{\kappa}=(\kappa_{1},\kappa_{2})$ taken from a finite index set, $V:=V_{N}\cup V_{C}\cup V_{XE}\cup V_{YE}$, with $V_{N}\subset\mathbb{Z}^{2}$, $V_{XE}\subset\{(z_{1}+1/2,z_{2}):(z_{1},z_{2})\in\mathbb{Z}^{2}\}$, $V_{YE}\subset\{(z_{1},z_{2}+1/2):(z_{1},z_{2})\in\mathbb{Z}^{2}\}$, and $V_{C}\subset\{(z_{1}+1/2,z_{2}+1/2):(z_{1},z_{2})\in\mathbb{Z}^{2}\}$. Similar definitions are used for the discrete operator acting on degrees of freedom corresponding to mesh edges and cell centers. Note that the decomposition of the index set, $V$, naturally defines $L_{\Delta x}^{N}$ as a block operator, with each block corresponding to one type of mesh point. For the elasticity operator in (17), we obtain a block operator of size $16\times 16$, since there are four scalar functions in the system, each discretized using Q2 elements with four types of degrees of freedom.

With this, the Fourier representation (symbol), of an operator, $L_{\Delta x}$, denoted by $\widetilde{L}_{\Delta x}(\boldsymbol{\theta})$, with $\boldsymbol{\theta}\in(-\pi,\pi]^{2}$, is a block matrix that can be computed using a Fourier basis that accounts for the different types of mesh points of $G_{\Delta x}$. This Fourier basis is given by

[TABLE]

for $\boldsymbol{\theta}\in(-\pi,\pi]^{2}$, with

[TABLE]

Considering the discretization of a system of PDEs with $q$ types of degrees of freedom (corresponding to both different functions in the PDE system and different basis functions used in higher-order discretizations of a single function, e. g., $q=16$ for the elasticity system) on a quadrilateral grid, the discretized operator can be written as a block operator,

[TABLE]

with scalar differential operators, $L_{\Delta x}^{i,j}$, $i,j=1,\ldots,q$. The Fourier symbol, $\widetilde{\mathcal{L}}_{\Delta x}(\boldsymbol{\theta})$, of $\mathcal{L}_{\Delta x}$ is computed from the Fourier symbols of each block,

[TABLE]

noting that $L_{\Delta x}^{i,j}$ is a map from a function associated with degree of freedom $j$ to that associated with degree of freedom $i$, which may be defined on different types of meshpoints. Formally, we think of $\widetilde{\mathcal{L}}_{\Delta x}(\theta)$ as describing the action of $\mathcal{L}_{\Delta x}$ on the $q$-dimensional space given by linear combinations of the Fourier modes associated with each type of degree of freedom in the system. Written as a $q\times q$ matrix, $\widetilde{\mathcal{L}}_{\Delta x}(\theta)$ maps the coefficients of a vector-valued function, $\vec{u}_{\Delta x}$, in this basis to the coefficients of the function described by $\mathcal{L}_{\Delta x}\vec{u}_{\Delta x}$.

For the linear elasticity problem, the time-stepping operator, $\Phi$, is defined in Equation (17),

[TABLE]

with projection operator, $P=I-(\rho M+\Delta t^{2}\mu K)^{-1}BS^{-1}B^{T}$, and $S=B^{T}(\rho M+\Delta t^{2}\mu K)^{-1}B$. Noting that both the velocity and displacement degrees of freedom are two-dimensional vector fields of $Q2$ degrees of freedom, we have four scalar functions in the system, whose basis is naturally partitioned into four types, with nodal, edge-based, and cell-centered degrees of freedom. Thus, the Fourier symbol, $\widetilde{\Phi}$, of $\Phi$ is a block matrix of size $16\times 16$, computed from the Fourier symbols of its component parts,

[TABLE]

with $8\times 8$ symbols $\widetilde{P}=I-(\rho\widetilde{M}+\Delta t^{2}\mu\widetilde{K})^{-1}\widetilde{B}\widetilde{S}^{-1}\widetilde{B}^{T}$, $\widetilde{S}=\widetilde{B}^{T}(\rho\widetilde{M}+\Delta t^{2}\mu\widetilde{K})^{-1}\widetilde{B}$,

[TABLE]

The Fourier symbols of $M$, $K$, and $B$ are derived from standard calculations38, 39, and given for reference in Appendix A.

4.1.2 LFA in time

We consider the infinite uniform fine and coarse temporal grids,

[TABLE]

where the grid $G_{T}$ is derived from $G_{t}$ by multiplying the mesh size, $\Delta t$, by a positive coarsening factor, $m$. The fundamental quantities in LFA are the Fourier modes, given by the grid-functions,

[TABLE]

with frequency, $\omega$, taken to vary continuously in the interval $\left(-\frac{\pi}{m},2\pi-\frac{\pi}{m}\right]$, although any interval of length $2\pi$ could be used instead. Fourier modes on the coarse grid are defined analogously on the interval $(-\pi,\pi]$ as explained below. Considering a coarsening factor of $m$, a constant-stencil restriction operator maps $m$ fine-grid functions, the Fourier harmonics, into one coarse-grid function. More precisely, these $m$ functions consist of the mode associated with some base index $\omega^{(0)}\in\left(-\frac{\pi}{m},\frac{\pi}{m}\right]$ and those associated with the frequencies

[TABLE]

The Fourier harmonics associated with base frequency, $\omega^{(0)}$, define subspaces,

[TABLE]

that are left invariant by the iteration operator. As a consequence, using the matrix of Fourier modes, ordered by Fourier harmonics, we can block-diagonalize the infinite grid operators with blocks corresponding to the spaces of harmonic modes. Each block describes the action on the coefficients $\{\alpha_{p}\}_{p=0,\ldots,m-1}$ of the expansion,

[TABLE]

of a function $e_{l}\in\mathcal{E}_{t}^{\omega^{(0)}}$.

Instead of analyzing the error propagation on this basis for the space of harmonics, $\mathcal{E}_{t}^{\omega^{(0)}}$, we consider the transformed basis

[TABLE]

where $\widehat{E}_{r}$ denotes the vector with entries equal to one at time points $t_{Lm+r}$, for $L\in\mathbb{N}$, and zero otherwise. This transformed basis is motivated by the following: Consider a function, $e_{l}\in\mathcal{E}_{t}^{\omega^{(0)}}$, with $l=Lm+r$ for $L\in\mathbb{N}$ and $r\in\{0,1,\ldots,m-1\}$. Then,

[TABLE]

Thus, any function in $\mathcal{E}_{t}^{\omega^{(0)}}$ can be expressed in terms of the basis defining $\widehat{\mathcal{E}}_{t}^{\omega^{(0)}}$, with coefficients $\{\widehat{\alpha}_{r}\}_{r=0,\ldots,m-1}$ that depend only on the “offset” of fine-grid point $l=Lm+r$ from the coarse-grid indices where $r=0$. Moreover, implicit in this expression is an associated coarse-grid frequency of $m\omega^{(0)}$, so that $\operatorname{e}^{\imath\omega^{(0)}\cdot t_{Lm+r}/\Delta t}=\operatorname{e}^{\imath(m\omega^{(0)})\cdot t_{Lm+r}/\Delta T}$. In many senses, this is a more natural basis for the space of Fourier harmonics, since it directly associates a repeating pattern of coefficients in the basis of $\widehat{\mathcal{E}}_{t}^{\omega^{(0)}}$ with the temporal mesh, as illustrated in Figure 2.

For the three-level analysis, we consider a hierarchy of three time grids, with grid spacings $\Delta t$, $m\Delta t$, and $m_{2}m\Delta t$, respectively. Correspondingly, we consider base frequencies, $\omega^{(0)}\in\left(-\frac{\pi}{m_{2}m},\frac{\pi}{m_{2}m}\right]$, and time-series coefficients $\{\widehat{\alpha}_{r}\}_{r=0,\ldots,m_{2}m-1}$, $\{\widehat{\alpha}_{r}\}_{r=0,m,\ldots,(m_{2}-1)m}$, and $\widehat{\alpha}_{0}$, respectively. Figure 3 shows the time-series coefficients on the time grids for the case $m=4$ and $m_{2}=2$.

4.1.3 Space-time LFA

For simplicity, we only describe two-grid space-time LFA; the three-level analysis can be done using inductive arguments. We consider two infinite uniform grids, $\Omega_{h}:=G_{\Delta x}\otimes G_{t}$, and $\Omega_{H}:=G_{\Delta x}\otimes G_{T}$ with $G_{\Delta x},G_{t}$, and $G_{T}$ defined in Equations (25) and (27). The subscript $h$ represents the pair $({\Delta x},\Delta t)$, where ${\Delta x}$ denotes the mesh size in both spatial dimensions, and $\Delta t$ is the time-step size. The grid $\Omega_{H}$ is derived from $\Omega_{h}$ by multiplying the mesh size in the time dimension only, i. e., $H$ represents the pair $(\Delta x,m\Delta t)$. On $\Omega_{h}$, we consider an infinite-grid (multilevel) Toeplitz operator, $A_{h}$, corresponding to the discretization of a time-dependent PDE, or a system of PDEs, in 2D space. Since the operator, $A_{h}$, is (multilevel) Toeplitz, it can be block-diagonalized by a set of continuous space-time Fourier modes,

[TABLE]

Considering semicoarsening in time, we analyze the error propagation on the space $\mathcal{E}^{(\theta,\omega^{(0)})}:=\operatorname{span}\{\operatorname{e}^{\imath\boldsymbol{\theta}\cdot\mathbf{x}_{j,k}/\Delta x}\}\otimes\widehat{\mathcal{E}}_{t}^{\omega^{(0)}}$, for spatial frequencies $\boldsymbol{\theta}\in(-\pi,\pi]^{2}$ and temporal base frequencies $\omega^{(0)}\in\left(-\frac{\pi}{m},\frac{\pi}{m}\right]$, i. e., we compute the action of the operators in the two-grid cycle on the coefficients $\{\widehat{\alpha}_{r}\}_{r=0,\ldots,m-1}$ of the expansion

[TABLE]

of a function $e_{j,k,l}\in\mathcal{E}^{(\theta,\omega^{(0)})}$.

Definition 4.1**.**

Let $L_{h}$ be an infinite-grid (multilevel) Toeplitz operator, and for given $\boldsymbol{\theta}\in(-\pi,\pi]^{2},\omega^{(0)}\in\left(-\frac{\pi}{m},\frac{\pi}{m}\right]$, let

[TABLE]

be a function with uniquely defined coefficients $\{\widehat{\alpha}_{r}\}_{r=0,\ldots,m-1}$. The transformation of the coefficients under the operator $L_{h}$,

[TABLE]

defined by

[TABLE]

is called the space-time Fourier symbol, $\widehat{L}_{h}=\widehat{L}_{h}(\boldsymbol{\theta},\omega^{(0)})$, of $L_{h}$.

A complete convergence analysis of the convergence properties of a Parareal or MGRIT algorithm arises from computing the symbols of each individual operator of the corresponding iteration matrix, $E$, to obtain the symbol, $\widehat{E}(\boldsymbol{\theta},\omega^{(0)})$, of the iteration matrix as a whole. The asymptotic convergence behavior can then be predicted by calculating the asymptotic convergence factor

[TABLE]

Note that, in practice, LFA has to be taken from its infinite-grid setting, and we maximize, instead, over a finite set of frequency tuples, $(\boldsymbol{\theta},\omega^{(0)})$.

Similarly, we introduce the error reduction factor

[TABLE]

again maximizing over finite sets of values of $\boldsymbol{\theta}$ and $\omega^{(0)}$ to get a computable prediction. Since for both the Parareal and MGRIT approaches, iteration operators have only a single eigenvalue of zero, only the error reduction factor provides insight into the convergence behavior of a Parareal or MGRIT algorithm. Furthermore, due to the non-normality of the iteration operators, it is crucial to not only consider the iteration matrix itself, but also powers of the iteration matrix to examine the short-term convergence behavior. More precisely, for any initial error, $\boldsymbol{e}_{0}$, $\|E^{k}\boldsymbol{e}_{0}\|\leq\|E^{k}\|\|\boldsymbol{e}_{0}\|\leq\|E\|^{k}\|\boldsymbol{e}_{0}\|$, where $k$ denotes the number of iterations of the multigrid scheme. Thus, we introduce the error reduction factors

[TABLE]

providing the worst-case error reduction, i. e., an upper bound for the error reduction, in iteration $k$ for the time-periodic problem. In practice, we observe that this also provides upper bounds for the worst-case error reduction in the non-periodic case but this, of course, is no longer a rigorous bound.

Space-time Fourier symbols of MGRIT with F- and FCF-relaxation

Let $\widetilde{\Phi}_{\theta}$ and $\widetilde{\Phi}_{c,\theta}$ denote the $q\times q$ ($q\geq 1$) spatial Fourier symbols of the fine- and coarse-scale time integrators, $\Phi$ and $\Phi_{c}$, respectively, and let $I$ denote an identity matrix of size $q\times q$. Then, the space-time Fourier symbol of the fine-grid operator, $A$, is given by the $qm\times qm$ matrix

[TABLE]

The set of modes $\psi(\boldsymbol{\theta},m\omega^{(0)})$ for $\omega^{(0)}\in\left(-\frac{\pi}{m},\frac{\pi}{m}\right]$ is a complete set of space-time Fourier modes on the coarse grid, $\Omega_{H}$. Further, any function in $\mathcal{E}^{(\theta,\omega^{(0)})}$ is aliased on the coarse grid with the function $\widehat{\alpha}_{0}\psi(\boldsymbol{\theta},m\omega^{(0)})$ with a fixed coefficient $\widehat{\alpha}_{0}$. As a consequence, the space-time Fourier symbol of the coarse-grid operator, $A_{c}$, is a block matrix with one block of size $q\times q$, given by

[TABLE]

The space-time Fourier symbols of the interpolation and restriction operators, $P_{\Phi}$, and $R_{I}$, are block matrices with $m\times 1$ blocks or $1\times m$ blocks of size $q\times q$,

[TABLE]

Since the F- and FCF-relaxation operators, $S^{F}$ and $S^{FCF}$, are equal to $P_{\Phi}R_{I}$ and $P_{\Phi}(I-A_{S})R_{I}$, respectively, the space-time Fourier symbols of F- and FCF-relaxation are defined by

[TABLE]

This completes the definition of the space-time Fourier symbols of the operators for F- and FCF-relaxation and, by using the expressions (20) and (21), this defines the space-time Fourier symbols of the two-level methods as a whole.

A similar approach extends to the three-grid case, however, the grid hierarchy affects the block size of the space-time Fourier symbols. Considering temporal semicoarsening by factors of $m$ and $m_{2}$, respectively, to obtain the first and second coarse grids, the space-time Fourier symbol of the fine-grid operator, $A$, is the same as in the two-level case, but with $m_{2}m\times m_{2}m$ blocks of size $q\times q$ instead of $m\times m$ blocks of size $q\times q$. Denoting the spatial Fourier symbols of the time integrators, $\Phi_{c}$ and $\Phi_{cc}$ on the first and second coarse grid by $\widetilde{\Phi}_{c,\theta}$ and $\widetilde{\Phi}_{cc,\theta}$, respectively, the space-time Fourier symbol of the coarse-grid operator, $A_{c}$, on the first coarse grid is a block matrix with $m_{2}\times m_{2}$ blocks of size $q\times q$,

[TABLE]

and the space-time Fourier symbol of the coarse-grid operator, $A_{cc}$, on the second coarse grid is a block matrix with one block of size $q\times q$,

[TABLE]

The space-time Fourier symbols of the Schur complement coarse-grid operators on the first and second coarse grids, $A_{S}$ and $A_{c,S}$, are defined as those of the coarse-grid operators, $A_{c}$ and $A_{cc}$, with $\widetilde{\Phi}_{c,\theta}$ and $\widetilde{\Phi}_{cc,\theta}$ replaced by $\widetilde{\Phi}_{\theta}^{m}$ and $\widetilde{\Phi}_{\theta}^{m_{2}m}$.

The space-time Fourier symbols of the restriction operators, $R_{I}$ and $R_{c,I}$, from the fine grid to the first coarse grid and from the first coarse grid to the second coarse grid are block matrices with $m_{2}\times m_{2}m$ blocks of size $q\times q$ or $1\times m_{2}$ blocks of size $q\times q$, respectively,

\widehat{R}_{I}=\setcounter{MaxMatrixCols}{13}\begin{bmatrix}I&0&\cdots&0\\ &&&&I&0&\cdots&0\\ &&&&&&&&\ddots\\ &&&&&&&&&I&0&\cdots&0\end{bmatrix}\quad\text{and}\quad\widehat{R}_{c,I}=\begin{bmatrix}I&0&\cdots&0\end{bmatrix}.$$m blocks

Similarly, the space-time Fourier symbols of the interpolation operators, $P_{\Phi}$ and $P_{\Phi_{c}}$, from the first coarse grid to the fine grid and from the second coarse grid to the first coarse grid, respectively, are block matrices with $m_{2}m\times m_{2}$ blocks of size $q\times q$ or $m_{2}\times 1$ blocks of size $q\times q$, respectively,

[TABLE]

4.2 Semi-algebraic mode analysis

LFA focuses on the local character of the operators defining the multilevel algorithms. This means that the effects of boundary conditions, including the initial condition in the time dimension, are ignored. Unless long time intervals are considered, LFA can fail to produce its usual quality of predictive results for convection-dominated or parabolic problems 40, 14, 41, 13. Semi-algebraic mode analysis (SAMA) 15 is one applicable approach of mode analysis that enables accurate predictions of the convergence behavior for multigrid and related multilevel methods. The analysis combines spatial LFA with algebraic computations that account for the non-local character of the operators in time. Other applicable approaches, not considered here, include half-space mode analysis 40, 42, 43, 44, 45 that considers convergence on a discrete half-plane instead of the full infinite lattice used in LFA, and Fourier-Laplace analysis 20, 21, 22, 46, 47, based on using discrete Laplace transforms.

For simplicity, we describe SAMA for the two-level MGRIT algorithm with FCF-relaxation, represented by the two-level iteration matrix, $E^{FCF}$, defined in Equation (21),

[TABLE]

the analysis of other variants is done analogously by considering the corresponding iteration matrix. Motivated by the global block-Toeplitz space-time structure of the operators that define the iteration matrix, the idea of SAMA is to use a Fourier ansatz in space to block-diagonalize the spatial blocks of each individual operator. More precisely, we block-diagonalize the blocks of each operator using the matrix, $\Psi$, of discretized spatial Fourier modes,

[TABLE]

with $N_{x}$ Fourier frequency pairs, $\boldsymbol{\theta}\in(-\pi,\pi]^{2}$, sampled on a uniform quadrilateral mesh given by the tensor product of an equally-spaced mesh over the interval of length $2\pi$ with itself, assuming $N_{x}$ degrees of freedom in the spatial discretization. For example, the spatial blocks of a scalar coarse-grid operator, $A_{c}$, are block-diagonalized by computing $\mathcal{F}^{-1}A_{c}\mathcal{F}$, where $\mathcal{F}=I_{N_{T}+1}\otimes\Psi$. The resulting matrices (on the fine temporal grid) are then reordered from $N_{t}+1\times N_{t}+1$ blocks of size $N_{x}\times N_{x}$ to a block matrix with $N_{x}\times N_{x}$ blocks of size $N_{t}+1\times N_{t}+1$ to obtain a block diagonal structure, with each block corresponding to the evolution of one spatial Fourier mode over time. Applying SAMA to the iteration matrix, $E^{FCF}$, we obtain a block diagonal matrix,

[TABLE]

with a discrete choice of Fourier frequencies, $\boldsymbol{\theta}$, and

[TABLE]

where the diagonal blocks of the Fourier-transformed and permuted operators are denoted by $B_{\boldsymbol{\theta}}^{(\cdot)}$ marking the respective operator in the superscript and the spatial Fourier frequency pair, $\boldsymbol{\theta}$, in the subscript. Note that the grid hierarchy and size of the spatial Fourier symbol of the time integrators, $\Phi$ and $\Phi_{c}$, affect the block size of the transformed operators. Considering temporal semicoarsening by a factor of $m$ in the MGRIT approach and assuming that the spatial Fourier symbol of the time-integration operators, $\Phi$ and $\Phi_{c}$, are of size $q\times q$, the blocks $B_{\boldsymbol{\theta}}^{(A_{c})}$ and $B_{\theta}^{(A_{S})}$ are of size $(N_{t}/m+1)q\times(N_{t}/m+1)q$, and the blocks $B_{\boldsymbol{\theta}}^{(P_{\Phi})}$ and $B_{\boldsymbol{\theta}}^{(R_{I})}$ are of size $(N_{t}+1)q\times(N_{t}/m+1)q$ and $(N_{t}/m+1)q\times(N_{t}+1)q$, respectively.

The short-term convergence behavior of MGRIT with iteration matrix, $E$, can then be predicted by calculating the error reduction factor (corresponding to worst-case error reduction in iteration $k$)

[TABLE]

Similarly to LFA, in practice, we maximize over a finite set of frequencies, $\boldsymbol{\theta}$. We note that the prediction given by SAMA is exact in special cases, notably when $\Phi$ has (block-)circulant structure and we sample at $N_{x}$ evenly-spaced Fourier points for the underlying LFA symbols (so that LFA coincides with rigorous Fourier analysis). In that setting, (32) provides an exact worst-case convergence bound for MGRIT performance, corresponding to the worst-case initial space-time error, against which the quality of predictions by other approaches can be measured. Note, however, that for any given initial value for the underlying PDE problem, the frequency content of the initial error may or may not sample the worst-case modes and, so, this prediction can still be an over-estimate of the actual MGRIT convergence factor, consistent with Bal’s analysis for Parareal48.

SAMA block matrices for MGRIT with F- and FCF-relaxation

Consider temporal semicoarsening by factors of $m$ and $m_{2}$, respectively, to obtain the first and second coarse grids, and let $\widetilde{\Phi}_{\theta}$, $\widetilde{\Phi}_{c,\theta}$, and $\widetilde{\Phi}_{cc,\theta}$ again denote the spatial Fourier symbols of the time integrators, $\Phi$, $\Phi_{c}$, and $\Phi_{cc}$ on the fine, first, and second coarse grid, respectively; for a two-level hierarchy, consider the first coarse grid only. Then, the diagonal blocks of the Fourier-transformed and permuted operators, $A$, $A_{c}$, and $A_{cc}$ are block bidiagonal matrices with $(N_{t}+1)\times(N_{t}+1)$, $(N_{t}/m+1)\times(N_{t}/m+1)$, or $(N_{t}/(m_{2}m)+1)\times(N_{t}/(m_{2}m)+1)$ blocks, respectively, of size $q\times q$, given by

[TABLE]

The diagonal blocks of the Fourier-transformed and permuted Schur complement coarse-grid operators on the first and second coarse grids, $A_{S}$ and $A_{c,S}$, are defined analogously, with $-\widetilde{\Phi}_{\theta}^{m}$ and $-\widetilde{\Phi}_{\theta}^{m_{2}m}$ on the first subdiagonal.

The diagonal blocks of the Fourier-transformed and permuted restriction operators, $R_{I}$ and $R_{c,I}$, from the fine grid to the first coarse grid and from the first coarse grid to the second coarse grid are block matrices with $(N_{t}/m+1)\times(N_{t}+1)$ blocks of size $q\times q$ or $(N_{t}/(m_{2}m)+1)\times(N_{t}/m+1)$ blocks of size $q\times q$, respectively,

B_{\boldsymbol{\theta}}^{(R_{I})}=\setcounter{MaxMatrixCols}{14}\begin{bmatrix}I&0&\cdots&0\\ &&&&I&0&\cdots&0\\ &&&&&&&&\ddots\\ &&&&&&&&&I&0&\cdots&0\\ &&&&&&&&&&&&&I\end{bmatrix}\quad\text{and}\quad B_{\boldsymbol{\theta}}^{(R_{c,I})}=\setcounter{MaxMatrixCols}{14}\begin{bmatrix}I&0&\cdots&0\\ &&&&I&0&\cdots&0\\ &&&&&&&&\ddots\\ &&&&&&&&&I&0&\cdots&0\\ &&&&&&&&&&&&&I\end{bmatrix}.$$m blocks$m_{2}$ blocks

Similarly, the diagonal blocks of the Fourier-transformed and permuted interpolation operators, $P_{\Phi}$ and $P_{\Phi_{c}}$, from the first coarse grid to the fine grid and from the second coarse grid to the first coarse grid, respectively, are block matrices with $(N_{t}+1)\times(N_{t}/m+1)$ blocks of size $q\times q$ or $(N_{t}/m+1)\times(N_{t}/(m_{2}m)+1)$ blocks of size $q\times q$, respectively,

[TABLE]

4.3 Two-level reduction analysis

Two-level reduction analysis 17, 19 provides convergence bounds for two-level Parareal and MGRIT algorithms applied to linear time-stepping problems with fine- and coarse-scale time-stepping operators that can be diagonalized by the same set of eigenvectors. In particular, the analysis allows predictions of the convergence behavior for time-stepping problems arising from a method-of-lines approximation of a time-dependent PDE when considering a fixed spatial discretization in the time-grid hierarchy. In Section 4.3.1, we review the ideas behind the two-level reduction analysis for time-stepping problems arising from scalar PDEs and discuss how it can be combined with LFA in space. Section 4.3.2 is devoted to extending the analysis to time-stepping problems arising from systems of PDEs using finite-element discretizations for the discretization in space.

4.3.1 Two-level reduction analysis for scalar PDEs

We consider solving a time-stepping problem (19), arising from a scalar PDE, by two-level MGRIT. As above, let $\Phi$ and $\Phi_{c}$ denote the two time-stepping operators on the fine time grid with $N_{t}$ time intervals, and on the coarse time grid with $N_{T}=N_{t}/m$ time intervals, respectively. Furthermore, assume that $N_{x}$ degrees of freedom are used for the discretization in space, so that $\Phi$ and $\Phi_{c}$ are matrices of size $N_{x}\times N_{x}$. Motivated by the reduction aspect of MGRIT, in contrast to analyzing full iteration matrices as in the SAMA approach, we consider the iteration matrices only on the coarse grid,

[TABLE]

where $A_{c}$ and $A_{S}$ denote the coarse-grid operator and the Schur complement coarse-grid operator, respectively, introduced in Section 3.1.

Two-level reduction analysis is based on the same analysis techniques used in SAMA, but applied to the coarse-grid iteration matrix, $E_{\Delta}$, instead of to the fine-grid iteration matrix, $E$, of the algorithm. Furthermore, instead of using a Fourier ansatz in space, the eigenvectors of the fine-grid time-stepping operator are used to diagonalize the fine- and coarse-scale time integrators (under the assumption that they are simultaneously unitarily diagonalizable) and, thus, the spatial blocks of the coarse-grid iteration matrix, $E_{\Delta}$. More precisely, using the unitary transformation, $X$, that diagonalizes the fine- and coarse-scale time integrators, $\Phi$ and $\Phi_{c}$, the blocks of the coarse-grid iteration matrix, $E_{\Delta}$ are diagonalized by computing $\mathcal{X}^{-1}E_{\Delta}\mathcal{X}$, where $\mathcal{X}=I_{N_{T}+1}\otimes X$. Similarly to SAMA, the resulting matrix is then permuted to obtain a block diagonal structure, with each block corresponding to the evolution of one eigenvector over time. Denoting the eigenvectors of the fine- and coarse-scale time integrators, $\Phi$ and $\Phi_{c}$, by $\left\{\mathbf{x}_{n}\right\}$, and the corresponding eigenvalues by $\left\{\lambda_{n}\right\}$ and $\left\{\mu_{n}\right\}$, respectively, for $n=1,\ldots,N_{x}$, we obtain

[TABLE]

with

[TABLE]

Using the standard norm inequality $\|E_{\Delta}\|_{2}\leq\sqrt{\|E_{\Delta}\|_{1}\|E_{\Delta}\|_{\infty}}$ for an operator $E_{\Delta}$, the convergence behavior of a two-level method is then predicted by calculating the bounds,

[TABLE]

that have since been shown to be accurate to order $\mathcal{O}(1/N_{T})$19. Assuming that $|\mu_{n}|\neq 1$ for all $n=1,\ldots,N_{x}$, we obtain17

[TABLE]

and

[TABLE]

and, thus,

[TABLE]

Remark 4.2**.**

Note that the two-level reduction analysis requires solving an eigenvalue problem of the size of the degrees of freedom of the spatial discretization to compute the eigenvalues of the two time-integrators, which are assumed to be simultaneous unitarily diagonalizable. To avoid the high computational cost of this eigensolve for high spatial resolutions (and complications in the non-normal case), LFA in space can be applied, providing predictions of the spectra of the time-stepping operators. Thus, the eigenvalues $\left\{\lambda_{n}\right\}$ and $\left\{\mu_{n}\right\}$ of $\Phi$ and $\Phi_{c}$ can be replaced by the spatial Fourier symbols, $\left\{\widetilde{\Phi}_{\boldsymbol{\theta}}\right\}$ and $\left\{\widetilde{\Phi}_{c,\boldsymbol{\theta}}\right\}$, respectively, choosing $N_{x}$ frequency values for $\boldsymbol{\theta}$.

Remark 4.3**.**

Considering powers of the matrices $E_{\Delta,n}^{F}$ and $E_{\Delta,n}^{FCF}$, we can analogously derive bounds of the $L_{2}$-norm of powers of the coarse-grid iteration matrices and, thus, capture the exactness property of MGRIT. For $k\geq 2$, we obtain

\left(E_{\Delta,n}^{F}\right)^{k}=(\lambda_{n}^{m}-\mu_{n})^{k}\begin{bmatrix}0\\ \vdots\\ 0\\ 1&0&0&0\\ {k\choose 1}\mu_{n}&1&0&0&0\\[3.60004pt] {k+1\choose 2}\mu_{n}^{2}&{k\choose 1}\mu_{n}&1&0&0&0\\[3.60004pt] {k+2\choose 3}\mu_{n}^{3}&{k+1\choose 2}\mu_{n}^{2}&\hskip 5.0pt{k\choose 1}\mu_{n}&1&0&0&0\\[6.0pt] \vdots&&\ddots&\ddots&\ddots&\ddots&&\ddots\\[6.0pt] {N_{T}-1\choose N_{T}-k}\mu_{n}^{N_{T}-k}&\cdots&&{k+1\choose 2}\mu_{n}^{2}&{k\choose 1}\mu_{n}&1&0&\cdots&0\end{bmatrix}$$k* rows*$k$* columns*

and

\left(E_{\Delta,n}^{FCF}\right)^{k}=(\lambda_{n}^{m}-\mu_{n})^{k}(\lambda_{n}^{m})^{k}\begin{bmatrix}0\\ \vdots\\ 0\\ 1&0&0&0\\ {k\choose 1}\mu_{n}&1&0&0&0\\[3.60004pt] {k+1\choose 2}\mu_{n}^{2}&{k\choose 1}\mu_{n}&1&0&0&0\\[3.60004pt] {k+2\choose 3}\mu_{n}^{3}&{k+1\choose 2}\mu_{n}^{2}&{k\choose 1}\mu_{n}&1&0&0&0\\[6.0pt] \vdots&&\ddots&\ddots&\ddots&\ddots&&\ddots\\[6.0pt] {N_{T}-1-k\choose N_{T}-2k}\mu_{n}^{N_{T}-2k}&\cdots&&{k+1\choose 2}\mu_{n}^{2}&{k\choose 1}\mu_{n}&1&0&\cdots&0\end{bmatrix}$$2k* rows*$2k$* columns*

Thus,

[TABLE]

and

[TABLE]

Using the definition of the binomial coefficient, for $k\geq 2$ we obtain

[TABLE]

and

[TABLE]

4.3.2 Systems of PDEs

The assumption that a time-stepping operator can be diagonalized by a unitary transformation does not necessarily hold true for time-stepping operators arising from systems of PDEs, since such time integrators can easily be non-symmetric. However, when considering a two-level time-grid hierarchy and a fixed spatial discretization on both grid levels, the assumption that the fine- and coarse-scale time integrators, $\Phi$ and $\Phi_{c}$, can be simultaneously diagonalized may still hold true for time-stepping problems arising from systems of PDEs. One possible generalization of the two-level reduction analysis to handle this case is to derive bounds for the coarse-grid iteration matrix in a mass matrix-induced norm instead of in the $L_{2}$-norm17, 19. Here, to enable a comparison of the different analysis techniques for the systems case, we derive bounds of the $L_{2}$-norm of the coarse-grid iteration matrices.

We consider $q\times q$ block time-stepping operators, $\Phi$ and $\Phi_{c}$, arising from a semi-discretization in space of a system of $q$ scalar PDEs with $q$ unknown functions; note that in the case of a PDE system with vector equations, we first have to break down the vectors into their scalar components. Furthermore, let $E_{\Delta}$ be the coarse-grid iteration matrix of a two-level method, given in Equation (33), defined using the block time integrators, $\Phi$ and $\Phi_{c}$. To derive an upper bound on the $L_{2}$-norm of $E_{\Delta}$, we first consider a Fourier ansatz in space. More precisely, we use the Fourier matrix, $F$, of discretized Fourier modes of a basis that accounts for the $q\times q$-coupling within $\Phi$ and $\Phi_{c}$ and, possibly, different nodal coordinates. Note that $F$ is a square matrix of dimension equal to that of $\Phi$ and $\Phi_{c}$. We then reorder the transformed block matrix, $\mathcal{F}^{-1}E_{\Delta}\mathcal{F}$, where $\mathcal{F}=I_{N_{T}+1}\otimes F$, from $N_{T}+1\times N_{T}+1$ blocks of size $qN_{x}\times qN_{x}$, where $qN_{x}=\dim(\Phi)$, to a block-diagonal matrix with $N_{x}\times N_{x}$ blocks of size $q(N_{T}+1)\times q(N_{T}+1)$. For the two-level algorithm with $F$-relaxation, for example, we obtain

[TABLE]

and where $F^{-1}\Phi F=\operatorname{diag}(\widetilde{\Phi}_{1},\ldots,\widetilde{\Phi}_{N_{x}})$ and $F^{-1}\Phi_{c}F=\operatorname{diag}(\widetilde{\Phi}_{c,1},\ldots,\widetilde{\Phi}_{c,N_{x}})$. Note that the spatial Fourier symbols, $\widetilde{\Phi}_{n}$ and $\widetilde{\Phi}_{c,n}$, $n=1,\ldots,N_{x}$, are dense block matrices of size $q\times q$. Therefore, for each $n=1,\ldots,N_{x}$, we simultaneously diagonalize the Fourier symbols, $\widetilde{\Phi}_{n}$ and $\widetilde{\Phi}_{c,n}$, using the eigenvector matrix, $U_{n}$, of $\widetilde{\Phi}_{n}$,

[TABLE]

Finally, we reorder the transformed block-Toeplitz matrix with diagonal blocks, $\mathcal{U}_{n}^{-1}\widetilde{E}_{\Delta,n}\mathcal{U}_{n}$, where $\mathcal{U}_{n}=I_{N_{T}+1}\otimes U_{n}$, to a block-diagonal matrix with Toeplitz blocks. For a two-level method with F-relaxation, we obtain

[TABLE]

Using the norm computations from Section 4.3.1, we obtain the following

Lemma 4.4**.**

Let $\Phi$ and $\Phi_{c}$ be $q\times q$ block time-stepping operators, and let $\widetilde{\Phi}_{n}$ and $\widetilde{\Phi}_{c,n}$, $n=1,\ldots,N_{x}$, where $qN_{x}=\dim(\Phi)=\dim(\Phi_{c})$ be the corresponding spatial Fourier symbols. Furthermore, assume that for all $n=1,\ldots,N_{x}$, $\widetilde{\Phi}_{n}$ and $\widetilde{\Phi}_{c,n}$, can be simultaneously diagonalized with eigenvalues $\left\{\lambda_{n,l}\right\}_{l=1,\ldots,q}$ and $\left\{\mu_{n,l}\right\}_{l=1,\ldots,q}$, respectively, and that $|\mu_{n,l}|\neq 1$ for all $n=1,\ldots,N_{x}$, $l=1,\ldots,q$. Then, the $L_{2}$-norm of the coarse-grid iteration matrices, $E_{\Delta}^{F}$ and $E_{\Delta}^{FCF}$, of the two-level methods satisfy

[TABLE]

and

[TABLE]

where $\kappa(U_{n})$ denotes the condition number of $U_{n}$.

Proof 4.5**.**

Let $F$ denote the Fourier matrix that transforms $\Phi$ and $\Phi_{c}$ into block-diagonal form with blocks given by the spatial Fourier symbols, $\left\{\widetilde{\Phi}_{n}\right\}$ and $\left\{\widetilde{\Phi}_{c,n}\right\}$. Furthermore, let $\mathcal{P}$ denote the permutation that reorders the Fourier-transformed block matrix, $\mathcal{F}^{-1}E_{\Delta}\mathcal{F}$, where $\mathcal{F}=I_{N_{T}+1}\otimes F$, into global block-diagonal structure with $N_{x}\times N_{x}$ Toeplitz blocks, $\widetilde{E}_{\Delta,n}$, given in Equation (38) for the two-level method with F-relaxation. Both applied transformations are unitary transformations and, thus,

[TABLE]

The bound for the two-level method with F-relaxation can then be derived as

[TABLE]

with $\mathcal{U}_{n}=I_{N_{T}+1}\otimes U_{n}$, where $U_{n}$ is the transformation defined in Equation (39), for $n=1,\ldots,N_{x}$, and $\mathcal{Q}^{-1}\mathcal{U}_{n}^{-1}\widetilde{E}_{\Delta,n}^{F}\mathcal{U}_{n}\mathcal{Q}=\operatorname{diag}(E_{\Delta,n,l}^{F})_{l=1,\ldots,q}$, with permutation $\mathcal{Q}$ that reorders $\mathcal{U}_{n}^{-1}\widetilde{E}_{\Delta,n}^{F}\mathcal{U}_{n}$ into block-diagonal structure with Toeplitz-blocks, $E_{\Delta,n,l}^{F}$. The bound for the coarse-grid iteration matrix of the two-level method with $FCF$-relaxation is derived analogously.

Remark 4.6**.**

When $|\mu_{n,l}|=1$ for $n\in\{1,\ldots,N_{x}\}$, $l\in\{1,\ldots,q\}$, we obtain $\|E_{\Delta,n,l}^{F}\|_{1}=\|E_{\Delta,n,l}^{F}\|_{\infty}=|\lambda_{n,l}-\mu_{n,l}|N_{T}$ and $\|E_{\Delta,n,l}^{FCF}\|_{1}=\|E_{\Delta,n,l}^{FCF}\|_{\infty}=|\lambda_{n,l}-\mu_{n,l}|(N_{T}-1)|\lambda_{n,l}|^{m}$, and we can easily get analogous bounds to those in Equations (40)-(41), and (35).

4.4 Blending SAMA and reduction analysis

The closeness of SAMA and the two-level reduction analysis motivates applying ideas from SAMA in the two-level reduction analysis and vice versa. In one direction, we can easily make use of the bound commonly used in reduction analysis to improve the computability of the SAMA prediction, in place of the singular value computation that was used in prior work15. This yields

[TABLE]

where $E$ is the iteration matrix using either F- or FCF-relaxation. Similarly, we can compute the SAMA bound using only the coarse representation of the propagators, with $E_{\Delta}^{F}$ or $E_{\Delta}^{FCF}$. Both of these are explored below, in Section 5.1.

In the other direction, we can investigate the differences between the predictions made by reduction analysis on the coarse grid17 or consider the full fine-grid iteration matrix as in SAMA. Here, it is easy to derive the corresponding bounds, of

[TABLE]

with

[TABLE]

For multiple iterations, $k\geq 2$, these bounds become

[TABLE]

with

[TABLE]

4.5 Differences, advantages, and disadvantages of the three tools

The three mode analysis tools presented above are all applied to the MGRIT iteration matrices, $E$, defined in Equations (20) and (21), or to their coarse-level versions, $E_{\Delta}$. Each predicts the worst-case error reduction in the $k$-th iteration, $k\geq 1$, by providing an upper bound or an approximation of $\|E^{k}\|_{2}$ or of $\|E_{\Delta}^{k}\|_{2}$. The analysis of an iteration matrix, $E$, by means of the three tools can be summarized as follows.

Compute the spatial Fourier symbols of $\widetilde{\Phi}_{\boldsymbol{\theta}}$ and $\widetilde{\Phi}_{c,\boldsymbol{\theta}}$ (each of size $q\times q$) of the time integrators on the fine and on the coarse grids for $N_{\boldsymbol{\theta}}$ discrete values of $\boldsymbol{\theta}\in(-\pi,\pi]^{d}$, $d=1,2$. 2. 2.

**LFA

**

1. i)

Compute the space-time Fourier symbol $\widehat{E}(\boldsymbol{\theta},\omega^{(0)})$ (of size $mq\times mq$) for $N_{\omega}$ discrete values of $\omega^{(0)}\in(-\pi/m,\pi/m]$. 2. ii)

For each $k=1,2,\ldots$

   

Compute $\|(\widehat{E}(\boldsymbol{\theta},\omega^{(0)}))^{k}\|_{2}$ for all $(\boldsymbol{\theta},\omega^{(0)})$-tuples and maximize over the $N_{\boldsymbol{\theta}}N_{\omega}$ values.

**SAMA

**

1. i)

Compute the SAMA block $B_{\boldsymbol{\theta}}^{(E)}$ (of size $(N_{t}+1)q\times(N_{t}+1)q)$. 2. ii)

For each $k=1,2,\ldots$

   

Compute $\|(B_{\boldsymbol{\theta}}^{(E)})^{k}\|_{2}$ or $\sqrt{\|(B_{\boldsymbol{\theta}}^{(E)})^{k}\|_{1}\|(B_{\boldsymbol{\theta}}^{(E)})^{k}\|_{\infty}}$ for all $\boldsymbol{\theta}$ and maximize over the $N_{\boldsymbol{\theta}}$ values.

**RA

**

1. i)

Simultaneously diagonalize $\widetilde{\Phi}_{\boldsymbol{\theta}}$ and $\widetilde{\Phi}_{c,\boldsymbol{\theta}}$ (each of size $q\times q$) to obtain eigenvalues $\{\lambda_{\boldsymbol{\theta},l}\}_{l=1,\ldots,q}$ and $\{\mu_{\boldsymbol{\theta},l}\}_{l=1,\ldots,q}$. 2. ii)

For each $k=1,2,\ldots$

   

Compute bound of $\|E_{\Delta}^{k}\|_{2}$ or $\|E^{k}\|_{2}$ by using $\{\lambda_{\boldsymbol{\theta},l}\}$ and $\{\mu_{\boldsymbol{\theta},l}\}$ in formulas for $\|E_{\Delta}^{k}\|_{1}$ and $\|E_{\Delta}^{k}\|_{\infty}$ or for $\|E^{k}\|_{1}$ and $\|E^{k}\|_{\infty}$ for all $\boldsymbol{\theta}$ and $l$ and maximize over the $qN_{\boldsymbol{\theta}}$ values.

The most important advantages and disadvantages of the three mode analysis tools are as follows:

**LFA

**

• $+$

relatively low computational cost

• $-$

good predictions only for large numbers of time steps, since finiteness of time interval is ignored; exactness property not captured

**SAMA

**

• $+$

accurate predictions for scalar and systems cases; takes finiteness of time interval into account; captures exactness property

• $-$

high computational cost, especially for large numbers of time steps

**RA

**

• $+$

ease of computing a bound for small and large numbers of time steps; takes finiteness of time interval into account; captures exactness property when considering powers of iteration matrices as in Remark 4.3

• $-$

bound relies on the assumption of unitary diagonalizability

5 Numerical results

Numerical experiments presented in this section are organized in two parts: first, in Section 5.1, we compare and contrast the three mode analysis tools of Section 4. Secondly, Sections 5.2 and 5.3 are devoted to using appropriate tools for gaining insight into effects of model and algorithmic parameters on the convergence behavior of the methods presented in Section 3.

5.1 Comparing the three mode analysis tools

The three mode analysis tools differ, in particular, in the treatment of the time dimension. While space-time LFA uses a Fourier ansatz in space and time, SAMA couples LFA in space with algebraic computations in time, and the two-level reduction analysis considers error propagation only on the coarse time grid. In this section, we compare and contrast the predictions of the three methods for the two hyperbolic model problems, linear advection in one dimension and incompressible linear elasticity in two dimensions. Here, both problems are discretized on a space-time mesh of size $64\times 64$ or $64^{2}\times 64$, respectively, with spatial mesh size ${\Delta x}=1/2$ and time-step size $\Delta t=1/10$. For the linear advection problem, the flow speed is chosen to be $c=1$ and the material parameters of the elasticity problem are chosen to be $\mu=\rho=1$; the influence of these model parameters on the convergence behavior is considered in Sections 5.2.1 and 5.3, respectively.

In all of the mode analysis tools, we consider a Fourier ansatz in space. Note that for the two-level reduction analysis, this ansatz saves the cost of a computationally expensive solution of an eigenvalue problem, and (for all cases) is equivalent to rigorously analysing the spatial problems with periodic boundary conditions. The results presented here sample the Fourier frequency, $\theta\in(-\pi,\pi]$, or the Fourier frequency pair, $\boldsymbol{\theta}=(\theta_{x},\theta_{y})\in(-\pi,\pi]^{2}$, respectively, on a discrete mesh with spacing $h_{\theta}=\pi/32$. For the space-time LFA predictions, the temporal Fourier base frequency, $\omega_{0}\in(-\pi/m,\pi/m]$, is additionally sampled on a discrete mesh with spacing $h_{\omega}=\pi/32$. The impact of finer Fourier frequency meshes was negligible in the examples considered here. Figure 4 shows the error reduction factors, $\sigma_{\text{LFA}}$, $\sigma_{\text{SAMA}}$, and $\sigma_{\text{RA}}$, defined in Equations (30), (32), and (36)-(37), respectively, for the first 10 two-level iterations with F- and FCF-relaxation with factor-2 temporal coarsening applied to linear advection and linear elasticity. Note that the error reduction factors, $\sigma_{\text{LFA}}$ and $\sigma_{\text{SAMA}}$, of space-time LFA and of SAMA are based on measuring the $L_{2}$-norm of powers of the full iteration matrices, i. e., $\sigma_{\text{LFA}}$ and $\sigma_{\text{SAMA}}$ are upper bounds for the worst-case error reduction at all time points, whereas the error reduction factor, $\sigma_{\text{RA}}$, of the two-level reduction analysis provides bounds on the $L_{2}$-norm of powers of the coarse-grid iteration matrices, i. e., upper bounds for the worst-case error reduction only at C-points. Results show that LFA only predicts the initial convergence behavior, while both SAMA and RA also enable good predictions of short-term and long-term convergence behavior, covering the superlinear convergence of the methods, the effect of non-normality in early iterations, and the exactness property of the algorithms (not shown here). Furthermore, the results demonstrate the difference in predictivity of SAMA and RA for the scalar case and the systems case. While in the scalar case of linear advection, predictions of SAMA and RA are close, in the systems case of linear elasticity, the prediction of RA is pessimistic due to the necessary condition number factor involved in the bound of the $L_{2}$-norm of the operators.

Figure 5 compares error reduction factors predicted by SAMA with experimentally measured error reduction factors for the first 10 two-level iterations with F- and FCF-relaxation with factor-2 temporal coarsening applied to linear elasticity, i. e., for the setting considered in the left plot of Figure 4. As SAMA predicts the worst-case error reduction, MGRIT convergence is measured for three initial conditions, consistenting of either exclusively low or high frequency modes, or for a linear combination of low and high frequency modes. More precisely, the initial condition is chosen as $u_{0}(x)=\varphi(-\theta,x)+\varphi(\theta,x)=2\cos(\theta x)$, with $\theta=\pi/16$ or $\theta=5\pi/8$ and $u_{0}(x)=\varphi(-\pi/8,x)+\varphi(\pi/8,x)+\varphi(-15\pi/16,x)+\varphi(15\pi/16,x)=2\cos(\pi x/8)+2\cos(15\pi x/16)$. Results demonstrate that SAMA predictions provide a good upper bound for the error reduction in each iteration. Note that when considering a random initial space-time guess, convergence is essentially the same for all initial conditions, since the initial space-time error consists of all frequencies. In contrast, convergence is much faster for a low-frequency initial condition when considering a zero initial guess, since only low-frequency errors are present.

5.1.1 Combining SAMA and reduction analysis techniques

The above results show that, in the systems setting, SAMA provides better predictions than RA. However, when computing exact Euclidean operator norms of matrices with size equal to the number of time steps multiplied by the dimension of the spatial LFA symbol, SAMA becomes prohibitively expensive for large numbers of time steps. On the other hand, computing bounds instead of exact Euclidean operator norms makes the reduction analysis computationally tractable. Moreover, considering the iteration matrices only on the coarse grid instead of the full iteration matrices further reduces the cost of RA. Therefore, in this section, we explore combining the SAMA and RA approaches. Again, both model problems are discretized on a space-time mesh of size $64\times 64$ or $64^{2}\times 64$, respectively, with spatial mesh size ${\Delta x}=1/2$ and time-step size $\Delta t=1/10$. Furthermore, spatial Fourier frequencies are sampled on a discrete mesh with spacing $h_{\theta}=\pi/32$ in both dimensions.

In Figures 6 and 7, we explore the effects of considering coarse-grid iteration matrices or full iteration matrices (labeled C-pts or full in the legends) and of computing bounds or exact Euclidean operator norms (labeled 2-norm bound or 2-norm in the legends of the figures). More precisely, in addition to the error reduction factors, $\sigma_{\text{SAMA}}$ and $\sigma_{\text{RA}}$ (solid lines in the figures), also considered in Section 5.1, Figures 6 and 7 show various variants of SAMA- and RA-predicted error reduction factors. In particular, for RA, we consider using full iteration matrices instead of only coarse-grid iteration matrices. For SAMA, three additional error reduction factors are shown: one error reduction factor based on full iteration matrices, but computing bounds of the Euclidean operator norms, and two error reduction factors based on coarse-grid iteration matrices, computing either exact Euclidean operator norms or their bounds. Results using F- and FCF-relaxation show similar relationships between the different variants of error reduction factors. The difference between considering coarse-grid iteration matrices or full iteration matrices is at most a factor of about 1.4 in all cases. Since for a temporal coarsening factor of $m$ this factor is given by $\sqrt{m}$18, for small coarsening factors, considering only coarse-grid iteration matrices gives good estimates of actual error reduction factors that may lie in between the full and the C-point bounds or below both bounds. Furthermore, when using the bound on the 2-norm instead of computing the exact 2-norm, for the advection problem, we observe that error reduction factors increase by a factor of at most 1.7 for linear advection, and by a factor of at most 3.7 in the case of linear elasticity. Thus, for the analysis of MGRIT applied to the linear elasticity probem, considering coarse-grid iteration matrices and using the bound on the 2-norm for the SAMA prediction gives a practical improvement over RA, with predicted error reduction factors being a factor of about five or ten, respectively, smaller than those of RA, and not relying on the (unsatisfied) assumption of unitary diagonalizability.

5.1.2 Predictivity of space-time LFA

Motivated by the observation that long time intervals are needed for LFA to be predictive for parabolic problems, in this section, we investigate whether similar observations apply to hyperbolic problems by comparing space-time LFA predictions with SAMA predictions for the linear advection problem on time intervals of increasing lengths. More precisely, the linear advection problem is discretized on a space-time mesh of size $64\times N_{t}$ using various numbers of time steps, $N_{t}$, with a spatial mesh size of ${\Delta x}=1/2$ and with fixed time-step size $\Delta t=1/10$; again, a flow speed of $c=1$ and factor-2 temporal coarsening are considered. We compute the error reduction factors $\sigma_{\text{LFA}}$ and $\sigma_{\text{SAMA}}$ for the first 20 iterations of the two-level methods. The slopes of the best-fit lines of these error reduction factors as a function of the iterations can be used to compute the predicted average error reductions per iteration, $\overline{\sigma}_{\text{LFA}}$ and $\overline{\sigma}_{\text{SAMA}}$, respectively. Figure 8 shows these predicted average error reduction factors (left) and the difference, $|\overline{\sigma}_{\text{LFA}}-\overline{\sigma}_{\text{SAMA}}|$, of LFA- and SAMA-predicted average error reductions (right) as functions of increasing numbers of time steps, $N_{t}$. Results show that the difference between LFA and SAMA predictions decreases with increasing numbers of time steps. More precisely, while for $N_{t}=128$, LFA-predicted average error reduction factors using F- and FCF-relaxation are about 23% or 90%, respectively, larger than SAMA predictions, for $N_{t}=1024$, LFA predictions are only about 2% or 6%, respectively, larger. Thus, results suggest that space-time LFA is a feasible option for predicting the convergence behavior for large numbers of time steps, usually corresponding to long time intervals as $T=N_{t}\Delta t$. Note that for hyperbolic problems, considering long time intervals may be interesting in practice, in contrast to the case of parabolic problems for which LFA becomes predictive for time intervals that generally are longer than the diffusion time scale.

Interestingly, in results not shown here, we see that space-time LFA is much more predictive of behaviour over the first few MGRIT iterations than it is over later iterations. As a concrete example, in Figure 4, we see that the space-time LFA and SAMA predictions agree nearly perfectly for the first iteration, but diverge after this. When averaging over the first 10 iterations, in place of the first 20 used in Figure 8, we see improvements of about a factor of two in the differences in average error reduction predicted compared to the results shown in Figure 8. These results are offset by correspondingly higher errors in iterations 11-20, where SAMA accurately predicts improvements in convergence factors that are missed by LFA.

5.2 Investigating MGRIT convergence for linear advection

So far, we have only presented results for flow speed $c=1$ and two-level cycling with factor-2 temporal coarsening. In this section, we investigate the effects of other wave speeds, of multilevel cycling, and of other coarsening factors on convergence. In particular, we are interested in answering the question of which frequencies cause slow convergence when MGRIT performance degrades.

5.2.1 Influence of model parameters

The discrete advection problem (2) depends on the factor $\lambda:=(c\Delta t)/{\Delta x}$, which can be seen as an effective CFL number, relating the flow speed $c$ and the discretization parameters ${\Delta x}$ and $\Delta t$. To determine the effect of $\lambda$ on convergence, we consider predicted error reduction for varying flow speeds $c$, while keeping the discretization parameters ${\Delta x}$ and $\Delta t$ fixed. As SAMA provides the most accurate predictions among the three analysis tools and allows insights into error reduction as a function of spatial modes, in this section, we only consider SAMA. Figure 9 shows the error reduction factor $\sigma_{\text{SAMA}}$ for the first 20 two-level iterations with F- and FCF-relaxation and factor-2 temporal coarsening applied to linear advection with flow speeds $c=2$ and $c=.5$ (corresponding to $\lambda=.4$ and $\lambda=.1$), respectively, discretized on a space-time mesh of size $64\times 512$ with spatial mesh size ${\Delta x}=1/2$ and time-step size $\Delta t=1/10$. The right plot in Figure 9 details the error reduction for the first iteration, showing the SAMA-predicted error reduction as a function of the spatial Fourier modes. Convergence with both relaxation strategies is similar, with convergence degrading with increasing effective CFL. We note that the error reduction is slowest for a set of low frequency modes that are close to, but do not include, the constant mode at $\theta=0$, indicating that the convergence rate is limited by the coarse-grid correction process.

5.2.2 Effects of multilevel cycling and coarsening

In Parareal and MGRIT algorithms, computations on the coarsest temporal grid are sequential. Larger temporal coarsening factors reduce the size of the coarsest-grid problem and, thus, the cost of computations. On the other hand, large coarsening factors increase the cost of relaxation. To determine the effect of different coarsening strategies on convergence, we consider two- and three-level MGRIT variants applied to the advection problem discretized on a $64\times 256$ space-time grid with spatial mesh size ${\Delta x}=1/2$ and time-step size $\Delta t=1/10$. The left plot in Figure 10 shows SAMA-predicted error reduction factors for two-level iterations with F- and FCF-relaxation with factor-2 and with factor-4 temporal coarsening. Increasing the coarsening factor leads to slower convergence, especially for F-relaxation. When considering the predicted average error reduction, $\overline{\sigma}_{\text{SAMA}}$, over the first 10 iterations, average error reduction per iteration degrades from 0.13 to 0.31 for F-relaxation and from 0.11 to 0.24 for FCF-relaxation. Note that considering a three-level method with coarsening factors $m$ and $m_{2}$ leads to the same number of time points on the coarsest grid as a two-level method with a coarsening factor of $mm_{2}$. Therefore, we compare convergence of two-level iterations with F- and FCF-relaxation with factor-4 coarsening with convergence of three-level schemes using factor-2 coarsening in both coarsening steps. The right plot in Figure 10 shows error reduction factors for the two-level methods and for three-level V- and F-cycle methods. Neglecting cycle costs, V-cycles lead to slower convergence than the corresponding two-level schemes, while F-cycles converge faster than the two-level schemes. However, taking into account that three-level F-cycles are twice as expensive as three-level V-cycles, the convergence rates of three-level V-cycles are better than those of three-level F-cycles, with an average error reduction over the first 10 iterations of three-level F-cycles with F- or FCF-relaxation of about 0.15 and 0.12, respectively, compared to squared average error reduction rates of three-level V-cycles with F- or FCF-relaxation of about 0.11 and 0.08, respectively. Note that the two-level results in Figure 10 with $m=4$ and FCF-relaxation come close to achieving convergence due to the exactness property, since the coarse-grid problem has only 64 temporal meshpoints, so we expect exact convergence in 32 two-level iterations. This may be the underlying reason for the prominent difference in convergence between F- and FCF-relaxation over the later iterations.

5.3 Investigating MGRIT convergence for linear elasticity

Performing similar experiments for the linear elasticity problem as for the linear advection problem in the previous section is computationally challenging. Additionally to the spatial Fourier symbol being a block matrix as a result of the system structure of the problem and of the mixed finite-element discretization, we consider Fourier frequency pairs on a two-dimensional tensor-product mesh instead of on a one-dimensional mesh. Results in Section 5.1 show that LFA may not be a feasible option, and predictions of RA are pessimistic due to the necessary condition number factor in the computations. SAMA gives a practical improvement over RA, which is also true for SAMA when computing bounds instead of exact Euclidean operator norms. We therefore focus on computing the SAMA bounds given in Equation (42) when analyzing MGRIT convergence.

One remaining drawback of the modified SAMA analysis is that it still becomes expensive for large numbers of time steps. When considering larger coarsening factors, however, we need to increase the number of time steps due to the exactness property of MGRIT. To make the SAMA analysis computationally tractable for larger numbers of time steps, we sample spatial Fourier frequencies on a discrete mesh with spacing $h_{\theta}=\pi/16$. Figure 11 shows that with this mesh spacing, SAMA captures the dominant behavior and shows details in the spectral Fourier domain that are not visible with a coarser mesh spacing of $h_{\theta}=\pi/8$, and are not significantly different when considering a finer mesh spacing of $h_{\theta}=\pi/32$.

In Figure 12, we look at the effects of the coarsening factor on error reduction. The left-hand side of the figure plots error reduction factors for the first 16 iterations of two-level MGRIT with factor-2 and with factor-8 temporal coarsening. While using a coarsening factor of two leads to good convergence behavior of both methods, we see divergence for factor-8 temporal coarsening in the first iterations. Note that the exactness property drives convergence of the iteration with FCF-relaxation and factor-8 coarsening at later iterations. At the right of Figure 12, we plot error reduction factors of the tenth iteration for both schemes (top and bottom rows) and both coarsening strategies (left and right columns) as functions of the Fourier frequency pair $(\theta_{x},\theta_{y})$. Note that the overall “structure” of these plots is quite similar across the different algorithmic parameters (although they are plotted on slightly different axes to best show individual details), with excellent convergence at frequencies both close to and far from the origin in frequency space, and worst convergence achieved for a range of small (but not zero) frequencies. This is quite possibly related to phase errors between the fine- and coarse-scale propagators49; to what extent this mode analysis tool may provide insight into how to cure the poor convergence is a question left for future work.

The discrete elasticity problem, given in Equations (13)-(15), depends on the material parameters $\rho$ and $\mu$ as well as on the discretization parameters ${\Delta x}$ and $\Delta t$. To determine relevant parameter sets for convergence studies, we make the following observations: In (15), the scaling of ${\Delta x}$ in the discrete divergence operator and of $\Delta t$ clearly does not matter at all. In (14), the value of $\Delta t$ matters as an absolute. For standard finite-element discretizations on uniform meshes in two spatial dimensions (as we consider here), the mass matrix can be written as a scaling factor of $(\Delta x)^{2}$ times an operator that is independent of $\Delta x$, while the entries in the stiffness matrix are independent of $\Delta x$. Thus, we can rescale (13) by dividing through by $\rho{\Delta x}^{2}$, and rescaling $p_{i}$ to obtain

[TABLE]

Since the scaling on the two terms involving $K$ differs by a factor of $\Delta t$, there are two natural parameters to consider:$\nu:=\Delta t\mu/\rho{\Delta x}^{2}$ and $\Delta t$.

To perform a thorough set of experiments with these parameters with a reasonable computational complexity, we fix $N_{t}=128$. The left plot in Figure 13 shows that the choice $N_{t}=128$ is reasonable for experimenting with the parameters in the system since this choice ensures that convergence is not affected by the exactness property of MGRIT in the first few iterations. Furthermore, aside from effects of the exactness property for small $N_{t}$, performance of the two-level methods does not depend on the number of time points. Note that varying the number of time steps changes the length of the time interval. This length can also be controlled by varying the time-step size (on a uniform-in-time mesh as we consider here). In the right plot of Figure 13, we consider three different time-step sizes $\Delta t$ for fixed $N_{t}=128$ and $\rho=\mu=\nu=1$. Results show that as $\Delta t$ decreases, there is an increasing initial jump in the worst-case error, but asymptotically, convergence appears to be independent of the time-step size. Furthermore, convergence of the two-level methods with $F$- and with $FCF$-relaxation is similar. After 10 iterations, $FCF$-relaxation gives about a factor of six improvement in error reduction over $F$-relaxation at twice the cost per itertation.

To gain insight into effects of model parameters on convergence, we consider the parameter, $\nu$, factored in the form $\nu=(\Delta t/{\Delta x}^{2})(\mu/\rho)$, i. e., we group discretization and material parameters, and we study its effect on convergence of the two-level methods. More precisely, we look at convergence for fixed $\nu$ by simultaneously varying $\Delta t$ and ${\Delta x}^{2}$ or $\mu$ and $\rho$, respectively, and at error reduction for varying $\nu$ by fixing three of the four parameters and varying the remaining one. Figure 14 shows effects of the material parameters on error reduction when the discretization parameters $\Delta t$ and $\Delta x$ are held constant. The left-hand side of the figure plots error reduction factors for fixed ratio $\mu/\rho=1$, i. e., simultaneously varying $\mu$ and $\rho$ (corresponding to fixed $\nu$), while at the right, error reduction factors are shown for various ratios $\mu/\rho$ for fixed $\mu=1$ (corresponding to varying $\nu$). The plots show that performance of both two-level schemes only depends on the ratio $\mu/\rho$. In the limit of large $\mu/\rho$, convergence is extremely fast, especially for the two-level scheme with $FCF$-relaxation, and unsteady for $F$-relaxation. This behavior is not surprising as this limit corresponds to a stiff material and, thus, the solution is approaching zero, i. e., oscillations are rapidly damped because of the Euler discretization.

The ratio $\mu/\rho$ allows us to determine effects of the parameter $\nu$ on error reduction for large variations in $\nu$. Using the slopes of the best-fit lines as a function of the iterations, we compute the average error reduction per iteration. The right-hand side of Figure 15 shows the average error reduction over iterations two thru 10 as a function of the parameter $\nu$ for $\Delta t=1/4$, $\Delta x=1/2$, and $\mu=1$ fixed, and various values of $\rho$. In all cases, average error reduction is bounded by 0.5, showing good and robust convergence in a reasonable parameter regime. At the left of Figure 15, we plot average error reduction as a function of the time-step size for fixed parameter $\nu=16$. Results demonstrate that convergence does not change until we get to large time steps, again, indicating good and robust convergence.

6 Conclusions

Currently, a key challenge in the development of parallel-in-time algorithms is achieving scalable algorithmic performance for hyperbolic PDEs. While some insight has been gained by both trial-and-error and more systematic computational studies in recent years, predictive analytical tools to aid in this development have been lacking in the literature. Here, we examine and extend mode analysis approaches, long used in the spatial multigrid community, to examine performance of methods from the Parareal/MGRIT class for two hyperbolic model problems. Our extensions lead to tighter bounds on performance that can be computed more efficiently than those existing in the literature. When applied to these model problems, we gain some insight into what poses essential challenges in developing algorithms for hyperbolic PDEs, and what parameter regimes, both physical and computational, are more or less difficult to handle. As identified in other contexts49, 27, convergence of MGRIT appears to be limited by differences between the fine-grid and coarse-grid time propagators for a set of smooth spatial modes near zero frequency, but not at zero frequency. While this is not a unique observation, we point out that the mode analysis realization of this statement has been key to the developments in De Sterck et al. 27, who derive an optimization perspective on the construction of MGRIT coarse operators that is particularly effective for explicit and higher-order discretizations of the linear advection equation. This new approach offers substantially more robust MGRIT performance for the linear advection equation than yet seen in the literature. Another key observation for the problems and discretizations considered in this paper is that, while the FCF-relaxation that is typical in MGRIT is generally more effective than the F-relaxation typical in Parareal, the differences between these approaches in a two-level setting is generally not substantial; when FCF-relaxation works best, F-relaxation has good performance, too, and when F-relaxation is ineffective, FCF-relaxation is not a magic cure (unless the exactness property becomes significant). Note that this conclusion is not universal. For linear advection, FCF-relaxation may be beneficial in the case of higher order Runge-Kutta discretizations, as discussed by Dobrev et al. 17 and De Sterck et al. 27 While this work does not immediately give direction as to how to improve algorithmic performance, it offers predictive tools that can be used to identify and diagnose convergence difficulties and may help in designing or optimizing improved algorithms. Possible avenues for future research clearly include the application of these analysis tools to a broader class of PDEs and discretizations, as well as their extension to consider the case of convergence to the continuum solution, as considered for Parareal by Ruprecht49.

Appendix A Spatial Fourier Symbols for Elasticity Operator

The Fourier symbols of the $Q2$ mass and stiffness matrices using nodal basis functions can be computed using tensor products,

[TABLE]

and

[TABLE]

respectively, with symbols

[TABLE]

of the 1D $Q2$ mass and stiffness matrices, respectively38, 39. The Fourier symbols of the derivative operators, $\widetilde{B}_{x}$ and $\widetilde{B}_{y}$, are given by39

[TABLE]

and

[TABLE]

respectively. The Fourier symbols of the time integrators, $\Phi_{c}$ and $\Phi_{cc}$, on the first and second coarse grid can be derived analogously, simply by adjusting the value of $\Delta t$ in the definition of $\widetilde{\Phi}$.

Acknowledgments

The work of H.D.S and S.M. was partially supported by NSERC Discovery Grants.