An optimized Space-Time Multigrid algorithm for parabolic PDEs

TL;DR

This paper enhances a space-time multigrid algorithm for parabolic PDEs by deriving optimal smoothing parameters, proposing a new coarsening strategy, and optimizing the two-cycle process, leading to significant efficiency improvements especially for small time-space ratios.

Contribution

The paper introduces a new coarsening strategy and analytical smoothing parameter optimization for space-time multigrid methods, improving efficiency for parabolic PDEs.

Findings

Optimal smoothing parameter derived for full coarsening

New coarsening strategy simplifies implementation

Efficiency gains observed for small τ/h² ratios

Abstract

We investigate three directions to further improve the highly efficient Space-Time Multigrid algorithm with block-Jacobi smoother introduced in [GanNeu16]. First, we derive an analytical expression for the optimal smoothing parameter in the case of a full space-time coarsening strategy; second, we propose a new and efficient direct coarsening strategy which simplifies the code by preventing changes of coarsening regimes; and third, we also optimize the entire two cycle to investigate if further efficiency gains are possible. Especially, we show that our new coarsening strategy leads to a significant efficiency gain when the ratio $\tau/h^2$ is small, where $\tau$ and $h$ represent the time and space steps. Our analysis is performed for the heat equation in one spatial dimension, using centered finite differences in space and Backward Euler in time, but could be generalized to other situations. We also present numerical experiments that confirm our theoretical findings.

Equations

$$\{\begin{array}{ll}{\partial }_{t}u(x,t)={\partial }_{xx}u(x,t)+f(x,t),& (x,t)\in \Omega \times (0,T),\\ u(x,t)=0,& x\in \partial \Omega ,\\ u(x,0)=u_0(x),& x\in \Omega .\end{array}$$

$$\{\begin{array}{ll}{\partial }_t\mathbf{u}(t)=A_h\mathbf{u}(t)+\mathbf{f}(t),& t\in (0,T),\\ \mathbf{u}(0)={\mathbf{u}}_0,& \end{array}$$

$$\{\begin{array}{ll}{Q}_{\tau ,h}{\mathbf{u}}_{n+1}=B_{\tau ,h}{\mathbf{u}}_n+{\mathbf{f}}_{n+1}& n=0,\dots ,N_t-1,\\ {\mathbf{u}}_0={\mathbf{u}}_0,& \end{array}$$

$$\underset{=:L_{\tau ,h}}{\underbrace{\left(\begin{array}{cccc}{Q}_{\tau ,h}& & & \\ -B_{\tau ,h}& Q_{\tau ,h}& & \\ & \ddots & \ddots & \\ & & -B_{\tau ,h}& Q_{\tau ,h}\end{array}\right)}}\underset{=:\mathbf{u}}{\underbrace{\left(\begin{array}{c}{\mathbf{u}}_1\\ {\mathbf{u}}_2\\ ⋮\\ {\mathbf{u}}_{N_t}\end{array}\right)}}=\underset{=:\mathbf{f}}{\underbrace{\left(\begin{array}{c}{\mathbf{f}}_1+B_{\tau ,h}{\mathbf{u}}_0\\ {\mathbf{f}}_2\\ ⋮\\ {\mathbf{f}}_{N_t}\end{array}\right)}}.$$

$$S_{\tau ,h}^{\nu }=[I_{\tau ,h}-\omega D_{\tau ,h}^{-1}{L}_{\tau ,h}{]}^{\nu },$$

$$S_{\tau ,h}=[I_{\tau }\otimes (1-\omega )I_h-\omega {\Gamma }_{-1,\tau }\otimes Q_{\tau ,h}^{-1}{B}_{\tau ,h}].$$

$${\Theta }_t$$

$${\Theta }_x$$

$$\hat{u}({\theta }_t,{\theta }_x):=\frac{1}{N_t}\frac{1}{N_x}\sum _{n=1}^{N_t}\sum _{j=1}^{N_x}{u}_{n,j}{e}^{-in{\theta }_t}{e}^{-ij{\theta }_x}.$$

$$\mathbf{u}=\sum _{{\theta }_t\in {\Theta }_t}\sum _{{\theta }_x\in {\Theta }_x}\hat{u}({\theta }_t,{\theta }_x)\mathbf{\Phi }({\theta }_t,{\theta }_x).$$

$${\mu }_S(\omega )=\max \{|{\hat{S}}_{\tau ,h}(\omega ;{\theta }_t,{\theta }_x)|;({\theta }_t,{\theta }_x)\in {\Theta }^{high}\},$$

$${\omega }^{\ast }=\underset{\omega \in (0,1]}{\mathrm{argmin}}{\mu }_S(\omega ).$$

$${\hat{S}}_{\tau ,h}(\omega ;{\theta }_t,{\theta }_x)=1-\omega +\omega e^{-i{\theta }_t}R(z({\theta }_x)),$$

$$\begin{gathered} |{\hat{S}}_{\tau ,h}(\omega ;{\theta }_t,{\theta }_x){|}^2=(1-\omega {)}^2+2\omega (1-\omega )\cos {\theta }_{t}R(z({\theta }_x)) \\ +{\omega }^2{R}^2(z({\theta }_x)). \end{gathered}$$

$$R(z)=\frac{1}{1-z},$$

$$|{\hat{S}}_{\tau ,h}(\omega ;{\theta }_t,{\theta }_x){|}^2=(1-\omega {)}^2+\frac{2\omega (1-\omega )\cos {\theta }_t}{c_x}+\frac{{\omega }^2}{c_x^2}.$$

$$c:=1-2\sigma [\cos \frac{\pi }{2}-1]=1+2\sigma .$$

$${\theta }_t^{\ast }=\{\begin{array}{ll}0& \text{if no time coarsening,}\\ \pi /2& \text{if factor 2 time coarsening,}\\ \pi /4& \text{if factor 4 time coarsening.}\end{array}$$

$$2\omega (1-\omega )\cos {\theta }_{t}R(z({\theta }_x)).$$

$$({\theta }_t^{\ast },{\theta }_x^{\ast })=(\pi /2,0).$$

$$({\theta }_t^{\ast },{\theta }_x^{\ast })=(0,\pi /2).$$

$$({\theta }_t^{\ast },{\theta }_x^{\ast })=\{\begin{array}{ll}(0,\pi /2)& \text{if}(c,\omega )\in {\mathcal{D}}_1,\\ (\pi /2,0)& \text{otherwise,}\end{array}$$

$${\mathcal{D}}_1:=\left\{c\in [1,+\infty ),0<\omega \le \frac{2c}{c^2+2c-1}\right\}.$$

$$({\theta }_t^{\ast },{\theta }_x^{\ast })=\{\begin{array}{ll}(0,\pi /2)& \text{if}(c,\omega )\in {\mathcal{D}}_2,\\ (\pi /4,0)& \text{otherwise,}\end{array}$$

$${\mathcal{D}}_2:=\left\{c\in [1,\sqrt{2}],0<\omega \le \frac{\sqrt{2}{c}^2-2c}{(\sqrt{2}-1)c^2-2c+1}\right\}.$$

$$|{\hat{S}}_{\tau ,h}(\omega ;\pi /2,0){|}^2\le |{\hat{S}}_{\tau ,h}(\omega ;0,\pi /2){|}^2.$$

$$|{\hat{S}}_{\tau ,h}(\omega ;\pi /4,0){|}^2\le |{\hat{S}}_{\tau ,h}(\omega ;0,\pi /2){|}^2.$$

$${\omega }^{\ast }=\frac{1}{2}.$$

$${\omega }^{\ast }=1.$$

$${\omega }^{\ast }=\{\begin{array}{ll}\frac{1}{2}& \text{if}\sigma >\frac{1}{\sqrt{2}}\\ \frac{2c}{c^2+2c-1}& \text{otherwise.}\end{array}$$

$${\omega }^{\ast }=\{\begin{array}{ll}\frac{1}{2}& \text{if}\sigma >\frac{\sqrt{2}-2+\sqrt{2-\sqrt{2}}}{2}\\ \frac{\sqrt{2}{c}^2-2c}{(\sqrt{2}-1)c^2-2c+1}& \text{otherwise.}\end{array}$$

Full text

††journal: International Journal of Computing and Visualization in Science and Engineering

1 Introduction

During the last twenty years, faster computing times have been achieved through the increase of the number of processing cores available rather than the development of faster clock speeds. However, as spatial parallelism techniques saturate, developing parallel techniques for the time dimension becomes important. Popular examples of such methods include Parareal lions2001resolution and its multilevel version: Multigrid Reduction-in-Time (MGRIT) falgout2014parallel . For a review on Parallel-in-time methods, the reader is referred to gander201550 ; ong2020applications ; falgout2017multigrid . Yet, many of the methods studied in these references coarsen the problem only in the time dimension. Therefore they are not truly scalable when applied to space-time problems. In order to overcome this issue, new methods that also include spatial coarsening have been developed, see e.g. horton1995space ; fischer2005parareal ; neumuller2013space ; gander2016analysis ; ruprecht2014convergence ; schroder2021 ; gander2022 .

Among these space-time approaches, we are interested in a specific, efficient Space-Time Multigrid algorithm (STMG), originally introduced in gander2016analysis based on a block-Jacobi smoother. This method, which benefits from excellent scalability properties, has been successfully applied to several problems neumuller2013space ; neumuller2018fully ; neumuller2019parallel .

In gander2016analysis , one of the goals of the authors was to develop a space-time multigrid method based only on standard components that always permits coarsening in the time direction, in order to overcome a limitation in the early parabolic multigrid method Hackbusch:1984:PMG , leading to a faster and more scalable multigrid algorithm. In doing so, in gander2016analysis , emphasis was put on semi-coarsening in time during the optimization of the smoother with respect to the damping parameter. The analysis then reveals that the ratio $\sigma:=\frac{\tau}{h^{2}}$, where $\tau$ and $h$ denote the time and space steps, plays a key role in the convergence of the method. Especially, in order to guarantee fast convergence, specific coarsening strategies are needed for this ratio to stay above a given threshold. The authors in gander2016analysis then propose a strategy based on alternating semi-coarsening in time, and full space-time coarsening.

Our goal here is to further optimize the STMG algorithm from gander2016analysis in three ways. First, we optimize the damping parameter of the block-Jacobi smoother for the case of full space-time coarsening. Second, we propose a new coarsening strategy to control the ratio $\sigma$, which results in a simpler algorithm while maintaining efficiency. And third, we investigate if optimizing the entire two-level process like in lucero2020optimization can lead to further gains. The whole analysis is conducted in the specific case of the heat equation in one spatial dimension discretized using centered finite differences in space and the Backward Euler method in time, but our approach could be generalized to higher dimensions and other equations.

We begin with introducing the continuous and discrete versions of the model problem in Section 2. Then, in Section 3, we present the STMG algorithm in our specific context. In Section 4, we review some useful notations and results from Local Fourier Analysis. Next, we analyze the performance of the smoother and derive the optimal damping parameter for full space-time coarsening in Section 5. Section 6 is dedicated to the comparison of different coarsening strategies preventing $\sigma$ from becoming too small, as well as to the optimization of the whole two-level process. Finally, numerical results are presented in Section 7.

2 Model problem

We consider as our parabolic model problem the one-dimensional heat equation. Let $T>0$, the space-time domain is given by $\Omega\times(0,T)$ where $\Omega:=(0,1)$. Given a source term $f$ and an initial condition $u_{0}$, the strong formulation of the problem reads

[TABLE]

By (lions2012non, , Chapter 4), we know that if $f\in L^{2}(0,T;H^{-1}(\Omega))$ and $u_{0}\in L^{2}(\Omega)$, then there exists a unique solution $u\in L^{2}(0,T;H^{1}_{0}(\Omega))\cap C^{0}(0,T;L^{2}(\Omega))$ to problem (1). Under such regularity assumptions on the data, we may discretize the problem in space using a step $h$, which leads to the matrix differential equation

[TABLE]

where $A_{h}\in\mathbb{R}^{N_{x}\times N_{x}}$ is the discrete Laplace operator and $N_{x}$ is the number of space steps. Problem (2) can then be discretized in time with a time step $\tau$ using a Runge-Kutta method. Let $R(z)=:Q^{-1}(z)B(z)$ be the associated stability function wanner1996solving , where $Q$ and $B$ are given polynomials. For simplicity, we introduce $Q_{\tau,h}:=Q(\tau A_{h})$ and $B_{\tau,h}:=B(\tau A_{h})$. Using these notations, we finally obtain

[TABLE]

where $N_{t}$ denotes the number of time steps. The linear system (3) can then be rewritten as a matrix system,

[TABLE]

3 Space-Time Multigrid

System (4) can be solved using the STMG algorithm. It follows the same steps as traditional multigrid algorithms trottenberg2000multigrid , which yields in the case of one simple V-cycle:

pre-smoothing ($\nu_{1}$ steps); 2. 2.

computation of the residual and restriction to the coarse grid; 3. 3.

coarse grid solve; 4. 4.

prolongation to the fine grid and correction; 5. 5.

post-smoothing ($\nu_{2}$ steps).

In the following subsections, we will describe each of the different components of the algorithm in detail.

3.1 Smoothing

Each iteration of STMG starts with $\nu_{1}$ pre-smoothing steps and ends with $\nu_{2}$ post-smoothing steps. Smoothing is used to damp highly oscillating modes in the error so that it can be properly represented on the coarse grid. Following gander2016analysis , we use a damped block-Jacobi smoother,

[TABLE]

where $\omega$ is the damping parameter and $D_{\tau,h}$ is a block-diagonal square matrix of size $N_{t}N_{x}$ composed of $N_{t}$ times the block $Q_{\tau,h}$ on the diagonal.

Remark 1

Let $\Gamma_{-1,\tau}$ be the $N_{t}\times N_{t}$ matrix whose entries are ones on the lower-diagonal and zeros elsewhere. Using the definition of $L_{\tau,h}$ in (4), the block-Jacobi smoother is given by

[TABLE]

Remark 2

For the smoothing iterative procedure to be convergent, we need the spectral radius of the smoother $S_{\tau,h}$ to verify $\rho(S_{\tau,h})=|1-\omega|<1$, that is $\omega\in(0,2)$. In the following, we will restrict ourselves to case $\omega\in(0,1]$ in order to avoid non-uniform convergence, see gander2016analysis .

3.2 Transfer operators

To go from the coarse to the fine grid and back, we need restriction and prolongation operators for both space and time grids, which we generically denote by $R_{h}$, $P_{h}$ and $R_{\tau}$, $P_{\tau}$. Their definitions depend on the type of coarsening we do, and are sometimes identities:

semi-coarsening in space $x$: $P_{\tau}=R_{\tau}=I_{\tau}$; 2. 2.

semi-coarsening in time $t$: $P_{h}=R_{h}=I_{h}$; 3. 3.

full-coarsening in space and time; 4. 4.

new coarsening strategy (see details below).

Unless specified otherwise, when the transfer operators are not equal to the identity, the prolongation operator is taken as the linear interpolation and the restriction operator is computed using the full-weighting technique. Note that using different transfer operators would result in a similar analysis.

The approach used in gander2016analysis is to coarsen by the same factor in space and time at each level. At some point, the smoothing operator does not smooth very well (this will be further discussed in Section 5) and we switch to doing semi-coarsening in time. To simplify the analysis, this technique will consist of a factor $2$ space-time coarsening followed by a semi-coarsening in time (see the ”Original strategy” in Figure 1). This is a slower technique than the one initially proposed but it gives us a basis for comparison.

The coarsening technique we propose consists in coarsening by a factor $4$ in time and $2$ in space, therefore maintaining the ratio $\sigma=\tau/h^{2}$ constant (see the ”New strategy” in Figure 1). This will yield a computationally cheaper strategy, while maintaining the convergence properties of the original method, therefore making the STMG algorithm even more efficient.

3.3 Coarse-grid solve

The restricted residual is solved on the coarser grid, meaning that we either perform a direct solve (two-level method) or we repeat steps 1 to 5 on this coarser grid until we reach a grid for which a direct solve is considered to be cheap (multi-level method).

For the coarse problem, we have two options: either use a Galerkin approximation of the operator $L_{\tau,h}$, that is e.g. $L_{2\tau,2h}=(R_{\tau}\otimes R_{h})L_{\tau,h}(P_{\tau}\otimes P_{h})$, or a rediscretization with coarser space and/or time steps, in which case $L_{2\tau,2h}=I_{2\tau}\otimes Q_{2\tau,2h}+\Gamma_{-1,2\tau}\otimes B_{2\tau,2h}$. In our analysis, we will consider rediscretization.

4 Analysis techniques

The analysis performed in this article relies on Local Fourier Analysis (LFA), first introduced in brandt1977multi and made rigorous in brandt1994rigorous . For an introduction, we refer the reader to the well written books (trottenberg2000multigrid, , Chapter 3) and (hackbusch1985multi, , Section 8.2). To this end, we will assume that our problem is periodic in time and space.

First, let us recall the Fourier mode decomposition for the vector $\boldsymbol{u}\in\mathbb{R}^{N_{t}N_{x}}$. All the following definitions can be found in gander2016analysis .

Definition 1

Let $\boldsymbol{\Phi}:(-\pi,\pi]^{2}\to\mathbb{R}^{N_{t}N_{x}}$ be defined by $\Phi_{n,j}(\alpha,\beta):=e^{in\alpha}e^{ij\beta}$, for $1\leq n\leq N_{t}$ and $1\leq j\leq N_{x}$. The vector $\boldsymbol{\Phi}(\theta_{t},\theta_{x})$ is called the Fourier mode with frequency $(\theta_{t},\theta_{x})\in\Theta_{t}\times\Theta_{x}=:\Theta$ where

[TABLE]

Definition 2

Let $\boldsymbol{u}=(\boldsymbol{u}_{1},\boldsymbol{u}_{2},\dots,\boldsymbol{u}_{N_{t}})\in\mathbb{R}^{N_{t}N_{x}}$. For each frequency $(\theta_{t},\theta_{x})\in\Theta$, the Fourier coefficient $\hat{u}$ associated to $(\theta_{t},\theta_{x})$ is defined by

[TABLE]

Using these notations, the Fourier mode decomposition of $\boldsymbol{u}$ reads

[TABLE]

The space-time frequency domain $\Theta$ can be decomposed as the disjoint union of low frequency and high frequency domains $\Theta=\Theta^{low}\cup\Theta^{high}$. The definitions of these low/high frequency domains depend on the type of coarsening considered,

• •

semi-coarsening in time:

$\Theta^{low}=\Theta\cap(-\pi/2,\pi/2]\times(-\pi,\pi]$;

• •

semi-coarsening in space:

$\Theta^{low}=\Theta\cap(-\pi,\pi]\times(-\pi/2,\pi/2]$;

• •

full space-time coarsening:

$\Theta^{low}=\Theta\cap(-\pi/2,\pi/2]\times(-\pi/2,\pi/2]$;

• •

new coarsening strategy:

$\Theta^{low}=\Theta\cap(-\pi/4,\pi/4]\times(-\pi/2,\pi/2]$.

In each case, we directly deduce $\Theta^{high}=\Theta\setminus\Theta^{low}$, see Figure 2.

5 Analysis of the smoother

In this section, we begin with the analysis for semi-coarsening in space and time, and we deduce the optimal damping parameter $\omega^{*}$ for full space-time coarsening as well as for our new coarsening strategy (see Figure 1). Then, for the full space-time and new coarsening strategies, we compare the performance of the smoother with our optimized damping parameter to the one with the original choice $\omega^{*}=0.5$ made in gander2016analysis .

To compute the optimal damping parameter, we first need to consider the smoothing factor defined as

[TABLE]

where $\Theta^{high}$ is the set of high frequencies defined in Section 4 and $\hat{S}_{\tau,h}$ is the Fourier symbol associated to (5). Thus, computing $\mu_{S}$ is equivalent to finding the modes $(\theta_{t}^{*},\theta_{x}^{*})\in\Theta^{high}$ that degrade the smoothing factor the most, i.e., that maximize $|\hat{S}_{\tau,h}(\omega;\theta_{x},\theta_{t})|$. Now, given an expression of $\mu_{S}$, the optimal damping parameter is obtained by solving the optimization problem

[TABLE]

In the case of a discretization using centered finite differences in space and a generic Runge Kutta method in time, the Fourier coefficient $\hat{S}_{\tau,h}$ of (5) is given by

[TABLE]

where $z(\theta_{x})$ is the Fourier coefficient of the matrix $L_{h}$ in (2) scaled by $\tau$. To ease computations, we consider

[TABLE]

In particular, we will choose $R$ to be the stability function associated to the Backward Euler method,

[TABLE]

and $z(\theta_{x})=2\tau[\cos\theta_{x}-1]/h^{2}=2\sigma[\cos\theta_{x}-1]$ (from a centered spatial discretization).

Denote by $c_{x}:=1-2\sigma[\cos\theta_{x}-1]$. Then, the expression (11) can then explicitly be written as

[TABLE]

When $\theta_{x}=\frac{\pi}{2}$, we denote $c_{x}$ by

[TABLE]

Remark 3

As the term (13) is symmetric with respect to $\theta_{x}$ (due to the symmetry of the spatial discretization) and $\theta_{t}$, we can restrict our arguments to the high frequencies that lie in the upper right quadrant.

Lemma 1

Let $R$ be the stability function of a Runge-Kutta method such that $R(z)\geq 0$ for all $z\in\mathbb{R}$. For $\omega\in(0,1]$, the time mode that degrades the smoother the most is

[TABLE]

Proof

The term (11) can be separated into three terms, the first and third do not depend on $\theta_{t}$ and are positive. We only need to maximize the second term

[TABLE]

As $R(z)\geq 0$, we only need to focus on maximizing $\cos\theta_{t}$. For $\theta_{t}\in[0,\pi]$ (no coarsening), it is obtained for $\theta_{t}^{*}=0$. For $\theta_{t}\in[\pi/2,\pi]$ (factor $2$ coarsening in time), we have that $\theta_{t}^{*}=\pi/2$. For $\theta_{t}\in[\pi/4,\pi]$ (factor $4$ coarsening in time), we have that $\theta_{t}^{*}=\pi/4$.

Lemma 2

Let $R$ be the stability function of a Runge-Kutta method such that $R$ is a monotonically increasing function. For $\omega\in(0,1]$, the spatial mode that degrades the smoother the most is the one that maximizes $z(\theta_{x})$.

Proof

Proceeding as previously, we look at the terms in (11), the first one is a constant, the second one depends on $\theta_{t}$ and $\theta_{x}$ and the third one depends on $\theta_{x}$. By Lemma 1, $\cos\theta_{t}^{*}\geq 0$, thus in order to maximize the term (11), the term $R(z(\theta_{x}))$ needs to be maximized. Since $R$ is a monotonically increasing function, we need to maximize $z(\theta_{x})$.

Lemma 3

Let us consider the Backward Euler method whose stability function is given by (12). Then we are able to compute the modes $(\theta_{t}^{*},\theta_{x}^{*})$ for which the term (13) is the largest for each of the following coarsening strategies.

1. (i)

When considering semi-coarsening in time,

[TABLE] 2. (ii)

When considering semi-coarsening in space,

[TABLE] 3. (iii)

When considering full space-time coarsening,

[TABLE]

where the region $\mathcal{D}_{1}\subset(1,+\infty)\times(0,1]$ is defined by

[TABLE] 4. (iv)

When considering the new coarsening strategy,

[TABLE]

where the region $\mathcal{D}_{2}\subset(1,+\infty)\times(0,1]$ is defined by

[TABLE]

Proof

(i) If no coarsening in space is used, the high frequency spatial modes are in the set $[0,\pi]$. Thus by Lemma 2, we need to maximize $z(\theta_{x})=2\sigma[\cos\theta_{x}-1]$ which is achieved for $\theta_{x}^{*}=0$, and we get from Lemma 1 that $\theta_{x}^{*}=\pi/2$.

(ii) In the case of semi-coarsening in space, this maximum is achieved for $\theta_{x}^{*}=\pi/2\in[\pi/2,\pi]$. Besides, Lemma 1 yields $\theta_{t}^{*}=0$.

(iii) In the case of full space-time coarsening, we can see the domain of high frequencies as the union of the high frequencies associated to semi-coarsening in time and to semi-coarsening in space (see Figure 1). Thus, depending on the values of $\omega$ and $\sigma$, the modes $(\theta_{t}^{*},\theta_{x}^{*})$ that degrade the smoother the most are either those for semi-coarsening in time $(\pi/2,0)$ (we say that time dominates) or those for semi-coarsening in space $(0,\pi/2)$ (we say that space dominates). More specifically, we have that space dominates for values of $\omega\in(0,1]$ such that

[TABLE]

Solving this inequality for $\omega$ and $\sigma$ leads to $(c,\omega)\in\mathcal{D}_{1}$. In Figure 3 (top), this region $\mathcal{D}_{1}$ where space dominates is represented in orange, and the region where time dominates is represented in blue. The previous inequality becomes an equality when $\omega=\frac{2c}{c^{2}+2c-1}$, see the grey curve in Figure 3.

(iv) For the new coarsening strategy, we need to compute the optimal modes for semi-coarsening in time with a factor of $4$. It is given by $(\theta_{t}^{*},\theta_{x}^{*})=(\pi/4,0)$ using Lemmas 1 and 2. So this time, we have that $(\theta_{t}^{*},\theta_{x}^{*})=(\pi/4,0)$ when time dominates and $(0,\pi/2)$ when space dominates. As in the previous case, we know that space dominates for $\omega\in(0,1]$ such that

[TABLE]

Here, the solution is given by $(c,\omega)\in\mathcal{D}_{2}$, see Figure 3 (bottom), and some technical computations show that equality is reached for $\omega=\frac{\sqrt{2}c^{2}-2c}{(\sqrt{2}-1)c^{2}-2c+1}$.

Now that we have found the pairs $(\theta_{t}^{*},\theta_{x}^{*})$ for each coarsening strategy, we are able to compute $\mu_{S}$ in each case. The next step is to figure out the optimal values for the damping parameter $\omega$ by solving (9).

Theorem 5.1

Consider the heat equation (1) discretized using a second order centered finite difference scheme in space and the Backward Euler method in time, which we solve with the Space-Time Multigrid algorithm. Then, we are able to compute the optimal damping parameter $\omega^{*}$ for the following coarsening strategies.

1. (i)

When considering semi-coarsening in time,

[TABLE] 2. (ii)

When considering semi-coarsening in space,

[TABLE] 3. (iii)

When considering full space-time coarsening,

[TABLE] 4. (iv)

When considering the new coarsening strategy,

[TABLE]

Proof

(i) For semi-coarsening in time, to obtain the optimal values $\omega^{*}$, we consider $\mu_{S}$ as defined in (8), then compute $d\mu_{S}/d\omega$ and solve $d\mu_{S}/d\omega=0$. If we coarsen by a factor $2$,

[TABLE]

where solving for $\omega$ yields $\omega^{*}=1/2$. Now, if we coarsen by a factor $4$,

[TABLE]

where solving for $\omega$ yields again $\omega^{*}=1/2$.

(ii) In the case of semi-coarsening in space, we have

[TABLE]

We notice that $\mu_{S}$ is a parabola with a minimum at $\omega_{0}=\frac{c}{c-1}$, which is outside the domain $(0,1]$. Therefore it is a decreasing function of $\omega$ on that domain, and it follows that $\omega^{*}=1$.

(iii) In the case of full space-time coarsening, we can observe two regimes (see Lemma 3 and Figure 3): one where time dominates, and the other one where space dominates. Also, we know that $\omega\mapsto|\hat{S}_{\tau,h}(\omega;0,\pi/2)|^{2}$ and $\omega\mapsto|\hat{S}_{\tau,h}(\omega;\pi/2,0)|^{2}$ are parabolas attaining their minima on $(0,1]$ at $\omega_{0}=1$ and $\omega_{0}=1/2$. Therefore, for a given $\sigma$, there are two possible behaviours for the continuous function $\omega\mapsto\mu_{S}(\omega)$, see Figure 3 (top). First, if $\sigma>1/\sqrt{2}$, then $\mu_{S}$ is decreasing on $[0,1/2]$ and increasing on $[1/2,1]$, thus $\omega^{*}=1/2$. Second, if $\sigma<1/\sqrt{2}$, then $\mu_{S}$ is decreasing on $[0,\frac{2c}{c^{2}+2c-1}]$ and increasing on $[\frac{2c}{c^{2}+2c-1},1]$, thus $\omega^{*}=\frac{2c}{c^{2}+2c-1}$.

(iv) In the case of the new coarsening strategy, we proceed exactly in the same way since we know that $\omega\mapsto|\hat{S}_{\tau,h}(\omega;0,\pi/4)|^{2}$ is also a parabola attaining its minimum at $\omega_{0}=1/2$. Here again, for a given $\sigma$, there are two distinct cases, see Figure 3 (bottom). First, if $\sigma>\frac{\sqrt{2}-2+\sqrt{2-\sqrt{2}}}{2}$, then $\mu_{S}$ is decreasing on the interval $[0,1/2]$ and increasing on the interval $[1/2,1]$, thus $\omega^{*}=1/2$. Second, if $\sigma<\frac{\sqrt{2}-2+\sqrt{2-\sqrt{2}}}{2}$, then $\mu_{S}$ is decreasing on $[0,\frac{\sqrt{2}c^{2}-2c}{(\sqrt{2}-1)c^{2}-2c+1}]$ and increasing on $[\frac{\sqrt{2}c^{2}-2c}{(\sqrt{2}-1)c^{2}-2c+1},1]$, thus $\omega^{*}=\frac{\sqrt{2}c^{2}-2c}{(\sqrt{2}-1)c^{2}-2c+1}$.

Remark 4

It can be checked that the optimal damping parameters for full space-time coarsening as well as for the new strategy belong to $[1/2,1]$ for $\sigma\in[0,1/\sqrt{2}]$ and $\sigma\in[0,\frac{\sqrt{2}-2+\sqrt{2-\sqrt{2}}}{2}]$ respectively.

Figures 4(a) and 4(b) show a comparison between the smoothing factors corresponding to the choice of the damping parameter in gander2016analysis ($\omega=1/2$) and the ones proposed in Theorem 5.1 ($\omega=\omega^{*}$) for the original and the new coarsening strategies. In each case, as expected, choosing $\omega=\omega^{*}$ leads to a smaller smoothing factor, although the gain does not seem very important for the new strategy. In order to better illustrate this gain, we plotted in Figure 4(c) the ratio $\frac{\ln[\mu_{S}(\omega^{*})]}{\ln[\mu_{S}(0.5)]}$ as a function of $\sigma$, which we call the efficiency. This quantity is an approximation of the ratio $\frac{N^{default}}{N^{optimal}}$ where $N^{default}$ and $N^{optimal}$ are the numbers of smoothing iterations it takes to reach a given error if we take $\omega=1/2$ and $\omega=\omega^{*}$. In Figure 4(c), we see that considering an optimal damping parameter can give up to a factor $2$ efficiency with respect to the value proposed in gander2016analysis . This means that an iteration using the optimal damping parameter would give a similar smoothing of the solution as two iterations using the original one. Additionally, using an optimal damping parameter does not give any overhead: in any case, the value of $\sigma$ needs to be checked at each level to decide if semi-coarsening in time or full space-time coarsening is more appropriate.

Remark 5

Taking a look at the smoothing factors displayed in Figures 4(a) and 4(b), one may think that the full space-time coarsening strategy is much more efficient than the new coarsening strategy. Actually, comparing these values is not really fair as the two strategies do not correspond to the same grids. Indeed, the first one is for coarsening from $(\tau,h)$ to $(2\tau,2h)$ while the second one is for coarsening from $(\tau,h)$ to $(4\tau,2h)$. Therefore $\mu_{S}$ in Figure 4(b) is computed by taking the maximum over a set of high frequencies $\Theta^{high}$ larger than the one in the case of Figure 4(a), which explains why the values are significantly higher.

6 Comparison of different multigrid strategies

In this section, we are interested in the study of Space-Time Multigrid strategies for solving (1). The most natural strategy would be to build a multigrid cycle based on full space-time coarsening. This choice does not alter the geometric properties of the mesh and the gain obtained by each coarsening step is significant in terms of computational cost since both space and time dimensions are divided by two. However, if the ratio $\sigma$ is given by $\sigma^{0}:=\tau/h^{2}$ at the finest level, then at the second level it becomes $\sigma=\sigma^{0}/2$, and at the $n$-th level we end up with $\sigma=\sigma^{0}/2^{n-1}$. As illustrated in Figures 3 and 4, the smoothing factor $\mu_{S}$ tends to deteriorate when $\sigma$ gets too small. Thus, if we use too many levels with this strategy, we will end up with very inefficient smoothing steps, which might slow down the whole multigrid algorithm.

For the reasons mentioned above, we thus want to study strategies such that the ratio $\sigma$ remains close to its value $\sigma^{0}$ at the finest level along all different levels. Especially, we focus on the two strategies presented in Figure 1. In the first one (called original strategy), inspired by the alternating strategy introduced by Neumüller and Gander in (gander2016analysis, , Section 5), the simplest V-cycle consists of one full space-time coarsening, then one semi-coarsening in time. This leads to a three-level method, starting from the level $(\tau,h)$ and going down to the coarsest level $(4\tau,2h)$. In the second strategy (called new strategy), which we present for the first time in this manuscript, the simplest V-cycle consists of a one step coarsening from the level $(\tau,h)$ directly to the level $(4\tau,2h)$. Therefore it leads to a standard two-level method.

In what follows, we will derive the expressions of the multigrid iteration operators associated with these different strategies, and we will estimate their spectral radii in the case where the damping parameters of the smoothers are chosen according to Theorem 5.1. Then, we will be interested in the action of these multigrid operators on low frequency modes. Finally, we will try to find damping parameters that optimize the whole multi-level method instead of optimizing only the smoother.

6.1 Original strategy

For the simplest three-level V-cycle illustrated in Figure 1 (left), we consider $\nu_{1}$ and $\nu_{2}$ pre/post-smoothing steps at the finest level and $\eta_{1}$ and $\eta_{2}$ pre/post-smoothing steps at the level $(2\tau,2h)$.

Notation The operator associated with the space-time prolongation from a given coarse grid $(m\tau^{\prime},nh^{\prime})$ to the finer grid $(\tau^{\prime},h^{\prime})$ will be denoted by $P_{m,n}$, for $m,n\in\{1,2\}$. This notation does not take into account the grid $(\tau^{\prime},h^{\prime})$ considered therefore the definition of $P_{m,n}$ depends on the context, e.g. for $m=n=2$:

[TABLE]

Similarly, we denote the corresponding restriction operator from the grid $(\tau^{\prime},h^{\prime})$ to the coarser grid $(m\tau^{\prime},nh^{\prime})$ by $R_{m,n}=\frac{1}{2}P_{m,n}^{T}$.

Using the notations introduced above, we can write the expression of the multigrid iteration operator as

[TABLE]

where $\widetilde{L_{2\tau,2h}^{-1}}$ is an approximation of the inverse of $L_{2\tau,2h}$ given by

[TABLE]

Remark 6

The definition of $\widetilde{L_{2\tau,2h}^{-1}}$ reflects the fact that, instead of solving exactly on the grid $(2\tau,2h)$ with $L_{2\tau,2h}^{-1}$, we solve approximately using one V-cycle on a coarser grid starting with zero initial approximation.

In order to study the properties of the operator $M^{o}$, we compute its action on the Fourier modes $\boldsymbol{\Phi}(\theta_{t},\theta_{x})$, for all $(\theta_{t},\theta_{x})\in\Theta$.

Preliminary notions.

Unlike the smoothing operator, the multigrid iteration operator cannot be diagonalized by the Fourier modes. Indeed, for a given fine grid low frequency mode, when working on the coarse grid, it is not possible to distinguish this mode from a corresponding fine grid high frequency mode (often referred to as companion mode). This phenomenon is called aliasing. In terms of operator, when studying the action of $M^{o}$ on the fine grid modes, the aliasing phenomenon results in a coupling between low frequency modes and their high frequency companion modes. However, as we will show next, it is still possible to block diagonalize $M^{o}$ by gathering the fine grid modes in groups of companion modes, called spaces of harmonics. In order to explain this in detail, let us begin with introducing two mapping operators

Definition 3

Let $\gamma_{2}:(-\frac{\pi}{2},\frac{\pi}{2}]\to(-\pi,-\frac{\pi}{2}]\cup(\frac{\pi}{2},\pi]$ and $\gamma_{4}:(-\frac{\pi}{4},\frac{\pi}{4}]\to(-\frac{\pi}{2},-\frac{\pi}{4}]\cup(\frac{\pi}{4},\frac{\pi}{2}]$ be defined as

[TABLE]

where we have chosen the convention $\text{sign}(0)=-1$.

For the reader’s convenience, each low time frequency in $(-\frac{\pi}{4},\frac{\pi}{4}]$ will be denoted by $\check{\check{\theta}}_{t}$, and each low space frequency in $(-\frac{\pi}{2},\frac{\pi}{2}]$ will be denoted by $\check{\theta}_{x}$. Then, for each low time frequency, we will denote by $\hat{\check{\theta}}_{t}:=\gamma_{2}(\check{\check{\theta}}_{t})$, $\check{\hat{\theta}}_{t}:=\gamma_{4}(\check{\check{\theta}}_{t})$ and $\hat{\hat{\theta}}_{t}:=\gamma_{2}\circ\gamma_{4}(\check{\check{\theta}}_{t})$. Using the same notation for space frequencies, we get $\hat{\theta}_{x}:=\gamma_{2}(\check{\theta}_{x})$.

In the present section, the coarsest grid in the cycles studied is $(4\tau,2h)$, which means that the low frequency domain, corresponding to the frequencies that can be represented on this coarse grid, is given by $\Theta^{low}=(-\pi/4,\pi/4]\times(-\pi/2,\pi/2]$. Using the frequency mapping operators defined above, it is possible to rewrite the Fourier decomposition (7) as a sum on the low frequencies only.

Lemma 4

Let $\boldsymbol{u}=(\boldsymbol{u}_{1},\boldsymbol{u}_{2},\dots,\boldsymbol{u}_{N_{t}})\in\mathbb{R}^{N_{t}N_{x}}$. Then the vector $\boldsymbol{u}$ can be written as

[TABLE]

where we have used the notation

[TABLE]

Proof

The proof simply relies on the fact that the operators $\gamma_{2}:(-\frac{\pi}{2},\frac{\pi}{2}]\to(-\pi,-\frac{\pi}{2}]\cup(\frac{\pi}{2},\pi]$ and $\gamma_{4}:(-\frac{\pi}{4},\frac{\pi}{4}]\to(-\frac{\pi}{2},-\frac{\pi}{4}]\cup(\frac{\pi}{4},\frac{\pi}{2}]$ are one-to-one.

This result means that each vector $u$ can be expressed as a linear combination, over all $(\check{\check{\theta}}_{t},\check{\theta}_{x})\in\Theta^{low}$, of elements in the space of harmonics associated to the fine grid defined by

[TABLE]

For simplicity, we have chosen the notation $(\theta_{t},\theta_{x})$ to refer to the low frequency $(\check{\check{\theta}}_{t},\check{\theta}_{x})$. In the same spirit, we define the Fourier spaces associated with the coarser grids $(2\tau,2h)$ and $(4\tau,2h)$ as

[TABLE]

Especially, for each low frequency $(\theta_{t},\theta_{x})$, the eight fine grid modes spanning the space of harmonics get aliased to the same coarse grid mode $\boldsymbol{\Phi}(4\theta_{t},2\theta_{x})$ spanning the Fourier space associated to the coarsest grid. Therefore, in order to study the convergence properties of the operator $M^{o}$, we will focus on its action on these eight fine grid modes for each low frequency.

Action of $S_{\tau,h}$.

For each mode $\boldsymbol{\Phi}(\theta_{t},\theta_{x})$, we have already seen in (10) that

[TABLE]

For the reader’s convenience, since $\omega$ is fixed here, we have dropped the dependence on $\omega$ in $\hat{S}_{\tau,h}$. Now, we can deduce from (18) that the smoother maps the space of harmonics to itself. Furthermore, the action of the smoother on this space of harmonics takes the form of a diagonal matrix in $\mathbb{C}^{8\times 8}$:

[TABLE]

[TABLE]

Action of $L_{\tau,h}$.

For each mode $\boldsymbol{\Phi}(\theta_{t},\theta_{x})$, we have

[TABLE]

This yields the following operator acting on the space of harmonics:

[TABLE]

Action of $R_{2,2}$.

In the context of the restriction from the fine grid $(\tau,h)$, for each mode $\boldsymbol{\Phi}(\theta_{t},\theta_{x})$, we have

[TABLE]

Since $2\gamma_{2}(\theta)\equiv 2\theta$ for all $\theta\in(-\pi/2,\pi/2]$, it is clear from the previous equation that the first four modes in (17) will be turned into the same “coarser” mode $\boldsymbol{\Phi}(2\check{\theta}_{t},2\theta_{x})$. Similarly, the last four modes in (17) will be turned into $\boldsymbol{\Phi}(2\hat{\theta}_{t},2\theta_{x})$. In matrix form, this takes the form of the operator

[TABLE]

Action of $P_{2,2}$.

Since the prolongation operator can be written as $P_{2,2}=2R_{2,2}^{T}$, we deduce that its action on the Fourier space associated with the grid $(2\tau,2h)$ can be modeled by the matrix operator

[TABLE]

The next step is to express in a similar way the action of the approximate inverse $(\widetilde{L_{2\tau,2h}^{-1}})$. In order to do this, we first need to study the action of the operators involved in its expression.

Action of $L_{2\tau,2h}$.

In the case where we consider the rediscretization approach rather than the Galerkin approach, we have $L_{2\tau,2h}=I_{2\tau}\otimes Q_{2\tau,2h}+\Gamma_{-1,2\tau}\otimes B_{2\tau,2h}$. Therefore we obtain the following operator acting on the Fourier space:

[TABLE]

where we have used the notation

[TABLE]

Action of $S_{2\tau,2h}$.

Adapting the definition of the damped block Jacobi smoother (5) to the grid $(2\tau,2h)$, we get

[TABLE]

where $D_{2\tau,2h}:=\text{diag}\{Q_{2\tau,2h}\}$. We may now use this expression to study the action of this smoothing operator on the Fourier space, which leads to

[TABLE]

where we have used the notation

[TABLE]

Action of $R_{2,1}$.

In the context of the grid $(2\tau,2h)$, we have $R_{2,1}=R_{2\tau\to 4\tau}\otimes I_{2h}$. Therefore we obtain the following matrix operator acting on the Fourier space:

[TABLE]

Action of $P_{2,1}$.

As in the case of $P_{2,2}$, we can deduce the representation of $P_{2,1}$ in Fourier space from the one of $R_{2,1}$

[TABLE]

Action of $L_{4\tau,2h}$.

Considering again the rediscretization approach, we proceed as in the case of $L_{2\tau,2h}$, which yields here the operator

[TABLE]

where we have used the notation

[TABLE]

Now we have studied the action of each operator involved in the expression of $M^{o}$, it is clear that the operator $\widehat{M}^{o}$ representing the action of $M^{o}$ on the space of harmonics verifies the following mapping property

[TABLE]

Moreover, for each low frequency $(\theta_{t},\theta_{x})$, we are now in a position to compute the spectral radius $\rho^{o}(\theta_{t},\theta_{x})$ of this $8\times 8$ complex-valued matrix. Especially, we are interested in the convergence factor of the multigrid method obtained with this coarsening strategy, i.e. we would like to compute

[TABLE]

6.2 New strategy

For the simplest two-level V-cycle illustrated in Figure 1 (right), with $\nu_{1}$ and $\nu_{2}$ pre/post-smoothing steps, we get the multigrid iteration operator

[TABLE]

where $P_{4,2}:=P_{2,2}P_{2,1}$ and $R_{4,2}:=R_{2,1}R_{2,2}=\frac{1}{4}P_{4,2}^{T}$.

Similarly to the analysis in Subsection 6.1, we need to compute the Fourier symbols for each of the operators in (22). For the operators $S_{\tau,h}$, $L_{\tau,h}$ and $L_{2\tau,2h}$ please refer to the discussion in Subsection 6.1.

Action of $R_{4,2}$.

Given the definition of $R_{4,2}$ above, we have in the context of a fine grid $(\tau,h)$ that $R_{4,2}=(R_{2\tau\to 4\tau}R_{\tau\to 2\tau})\otimes R_{\tau\to 2\tau}$. By writing the explicit form of the operator $R_{2\tau\to 4\tau}R_{\tau\to 2\tau}$, we obtain

[TABLE]

where we have used the trigonometric relation (obtained through Werner’s formula)

[TABLE]

We can therefore represent the action of this restriction operator on the space of harmonics by

[TABLE]

Action of $P_{4,2}$.

By the definition of $P_{4,2}$, we may represent it by

[TABLE]

Remark 7

If $\eta_{1}=\eta_{2}=0$ in the original coarsening strategy, then both coarsening strategies are equivalent. Indeed, in this case the approximation (15) of the inverse of $L_{2\tau,2h}$ becomes

[TABLE]

Then, since $P_{4,2}=P_{2,2}P_{2,1}$ and $R_{4,2}=R_{2,1}R_{2,2}$, it follows that $M^{o}=M^{n}$.

As in the case of the previous strategy, it follows from these results that the operator $\widehat{M}^{n}$ representing the action of $M^{n}$ on the space of harmonics verifies the same mapping property

[TABLE]

We may now study the following quantity, which depends on $\sigma$,

[TABLE]

where for each frequency $(\theta_{t},\theta_{x})$, $\rho^{n}(\theta_{t},\theta_{x})$ denotes the spectral radius of $\widehat{M}^{n}$.

For $\nu_{1}=\eta_{1}=3$ pre-smoothing and $\nu_{2}=\eta_{2}=3$ post-smoothing steps, we plotted in Figures 5 and 6 the dependence of $\bar{\rho}^{o}$ and $\bar{\rho}^{n}$ on $\sigma$.

In Figure 5, we see that for $\omega=1/2$, $\bar{\rho}^{o}$ is always smaller than $\bar{\rho}^{n}$, which makes sense since there is an additional intermediate smoothing step in the original strategy (at the level $(2\tau,2h)$) compared to the new strategy. Moreover, we also see that the convergence factors are very similar except for values of $\sigma$ approximately in $(0.02,400)$. This corresponds to the range of values of $\sigma$ for which the intermediate smoothing step has the highest impact. Now in Figure 6, we plotted the two spectral radii again but this time choosing the damping parameter $\omega=\omega^{*}$ as prescribed by Theorem 5.1. For intermediate and high values of $\sigma$, the remarks made in the case $\omega=\frac{1}{2}$ remain valid. But for low values of $\sigma$, more specifically for $\sigma\leq 0.04$, we see that $\bar{\rho}^{n}$ is significantly smaller than $\bar{\rho}^{o}$. This means that for such values of $\sigma$, the new strategy is not only cheaper than the original strategy but also faster. One possible explanation for this surprisingly nice property is that, when $\sigma$ is small, $\omega^{*}$ is close to the optimal damping parameter for the multi-level method. The reader is referred to Subsection 6.4 for a more detailed discussion of this issue.

6.3 Action of the multigrid operators on $\Theta^{low}$

In order to illustrate the action of the two-grid operators $M^{o}$ and $M^{n}$ on low frequency modes, we begin with choosing a simple linear combination of all low frequency modes as input

[TABLE]

For any combination of the space-time frequency modes

[TABLE]

with coefficients $v(\theta_{t},\theta_{x})\in\mathbb{C}$, we choose to represent it graphically by the moduli of its coefficients $|v|$ defined on the space-time frequency space $\Theta$. Following this choice, we plotted the graph of $\boldsymbol{V}_{0}$ in Figure 7, illustrating that $|v_{0}|$ takes the value 1 if $(\theta_{t},\theta_{x})\in\Theta^{low}$ and 0 otherwise. For $\nu_{1}=\eta_{1}=3$ pre-smoothing and $\nu_{2}=\eta_{2}=3$ post-smoothing steps, we have computed $\widehat{M}^{o}\boldsymbol{V}_{0}$ and $\widehat{M}^{n}\boldsymbol{V}_{0}$ for $\omega=\omega^{*}$ and different values of $\sigma$, and plotted the results in Figures 8 and 9. According to these results, it appears that the intermediate smoothing step makes the original strategy more efficient than the new strategy when applied to the low frequency input $\boldsymbol{V}_{0}$ for both values $\sigma=1$ and $\sigma=10^{-2}$. Especially, the original strategy is much more efficient to damp space-time frequency modes $\boldsymbol{\Phi}(\theta_{t},\theta_{x})$ for which $\theta_{t}\in[-\frac{3\pi}{4},-\frac{\pi}{4}]\cup[\frac{\pi}{4},\frac{3\pi}{4}]$. An explanation for this is that the intermediate smoothing step is related to the coarsening from the finer grid $(2\tau,2h)$ to the coarser grid $(4\tau,2h)$. Therefore it is optimized to damp efficiently associated high frequency modes, i.e. modes for which $\theta_{t}\in[-\frac{\pi}{2},-\frac{\pi}{4}]\cup[\frac{\pi}{4},\frac{\pi}{2}]$. This damping effect also impacts the corresponding companion modes with respect to the coarsening $(\tau,h)\to(2\tau,2h)$, i.e. modes for which $\theta_{t}\in[-\frac{3\pi}{4},-\frac{\pi}{2}]\cup[\frac{\pi}{2},\frac{3\pi}{4}]$.

For $\sigma=10^{-2}$, we see in Figure 8 that the coefficients of $\widehat{M}^{n}\boldsymbol{V}_{0}$ are maximal near the boundaries of $\Theta^{low}$. Thus for the original strategy, the coefficients of $\widehat{M}^{o}\boldsymbol{V}_{0}$ near the boundaries $\left\{\theta_{t}=\pm\frac{\pi}{4}\right\}$ are well damped thanks to the effect of the intermediate smoothing step but the ones near the boundaries $\left\{\theta_{x}=\pm\frac{\pi}{2}\right\}$ are similar to the ones for the new strategy. Hence the global performances of both strategies are rather similar when applied to the input $\phi_{0}$.

Besides, in the case where $\sigma=1$, Figure 9 shows that the coefficients of $\widehat{M}^{n}\boldsymbol{V}_{0}$ are maximal for the modes ${(\theta_{t},\theta_{x})=\left(\pm\frac{\pi}{4},0\right)}$. So the damping effect of the intermediate smoothing step is more important than in the case $\sigma=10^{-2}$, and we end up with maximal coefficients of $\widehat{M}^{o}\boldsymbol{V}_{0}$ about half of the ones of $\widehat{M}^{n}\boldsymbol{V}_{0}$. Thus the original strategy is much more efficient than the new one in this case.

6.4 Optimizing the damping parameter for the multi-level method

In the previous subsections (except from Figure 5), we had chosen for each strategy a damping parameter equal to $\omega^{*}$ according to Theorem 5.1 in order to maximize the efficiency of the smoothing step. However, this choice might not be optimal when considering the multi-level strategies presented above. Indeed, there may exist values of $\omega$ such that the spectral radii $\bar{\rho}^{o}$ and $\bar{\rho}^{n}$ are minimal, for given values $\nu_{1}$, $\nu_{2}$, $\eta_{1}$ and $\eta_{2}$. Let us denote these values of damping parameter by $\omega_{opt}^{o}$ and $\omega_{opt}^{n}$. By definition, we have

[TABLE]

Since the expressions of $\bar{\rho}^{o}$ and $\bar{\rho}^{o}$ are already quite intricate, we do not seek to find the analytical expressions of $\omega_{opt}^{o}$ and $\omega_{opt}^{n}$. Instead, we compute them numerically as a function of $\sigma$, for fixed numbers of smoothing steps $\nu_{1}$, $\nu_{2}$, $\eta_{1}$ and $\eta_{2}$.

In Figures 10 and 11, we plotted the values of $\omega_{opt}^{o}$ and $\omega_{opt}^{n}$ with respect to $\sigma$ for $\nu_{1}=\nu_{2}=\eta_{1}=\eta_{2}=3$, and we compared them to $\omega^{*}$ for each strategy. Then, in Figure 12, we have compared $\bar{\rho}^{o}$ and $\bar{\rho}^{n}$ when taking the optimal damping parameters $\omega_{opt}^{o}$ and $\omega_{opt}^{n}$ computed numerically. We see that for these optimal versions, the original strategy always performs better than the new one, which was expected since it contains additional intermediate pre/post-smoothing steps. Thus when taking the best out of each strategy, we obtain a more efficient (though more expensive) algorithm with the original strategy.

7 Numerical experiments

To illustrate our theory, we compute the solution to the heat equation (1) on the time interval $[0,0.1]$ with the source term

[TABLE]

and initial value ${u_{0}(x)=0}$. We start the iterative process with a random initial guess.

We solve the problem in two different regimes: first for $\sigma$ small (Figure 13, where both coarsening strategies with optimal parameters have an efficiency close to $2$), then for $\sigma$ large (Figure 14, where both the default and optimal strategies share the same smoothing factor $\omega^{*}=0.5$). In all cases, $\nu_{1}=\eta_{1}=3$ pre-smoothing and $\nu_{2}=\eta_{2}=3$ post-smoothing steps were performed. The error is measured in the norm $L^{\infty}(0,T;L^{2}(\Omega))$.

For $\sigma\approx 1.56\times 10^{-1}$, we observe in Figure 13(a) that for the default parameter $\omega=0.5$, the convergence is quite slow for both strategies. We also see that the errors for the two strategies are very close (the red line is close to the blue line). Besides, for this value of $\sigma$ and the original coarsening strategy, choosing the damping parameter in an optimal fashion (see Theorem 5.1) yields some gains in terms of efficiency. With this "optimal" version of the algorithm, we also observe that after 5 iterations of the new strategy, we get an error similar to the one when $\omega=0.5$ is used for 10 iterations. In terms of computing time, this is a gain of approximately $30\%$. However, it is worth noting that even further gains can be achieved by computing an optimal two-level damping parameter. In the considered regime, the new coarsening strategy is supposed to be much worse than the original one (see Figures 5 and 6), which is what we observe for all damping parameter choices. So using the new coarsening strategy is more beneficial when the value of $\sigma$ is larger, which is what we encounter in practical situations.

From Figure 14(a), we notice that in the regime where $\sigma$ is large, the convergence rate of the new coarsening strategy is slightly lower than for the original strategy, which is to be expected given the values of the convergence factors in Figures 5 and 6. However, when plotting the error as a function of the computational time (see Figure 14(b)), we observe a quicker convergence for the new coarsening strategy. This is a very interesting optimization as the original STMG method presented in gander2016analysis to which we compare our new strategy is considered to be extremely efficient in this regime. The results have been obtained using a sequential code and, in practice, the smoothing procedure would be computed in parallel, making it cheaper. Numerically, we observe a benefit to using the new strategy as soon as it is approximately $15\%$ cheaper than the original method in computational cost per iteration.

8 Conclusion

We derived optimal smoothing factors for the Space-Time Multigrid algorithm proposed in gander2016analysis , demonstrated their efficiency compared to the original method and verified it numerically. Additionally, the original approach suggests to do a coarsening by a factor $2$ in space and in time as long as $\tau/h^{2}>1/\sqrt{2}$ holds, and then to proceed to do coarsening in time only. In order to prevent this change of regimes, we suggested to coarsen by a factor $2$ in space and $4$ in time. Our theoretical analysis revealed that this new strategy leads to a multi-level algorithm which has a comparable convergence factor to the original strategy, except for intermediate values of the ratio $\sigma$, where it is slightly bigger. However, this new strategy is always competitive in practice as it is cheaper than the original one. Indeed, since it does not require any intermediate smoothing steps, it is slightly slower for intermediate values of $\sigma$, and it is faster for larger values of $\sigma$, where the original algorithm is known to perform extremely well. Finally, we also designed a numerical technique to optimize the multi-level method with respect to the damping parameter, therefore leading to an optimal version of the STMG algorithm in this context.

As far as future work is concerned, it would be interesting to extend the multigrid analysis to the case of higher order Runge Kutta methods for the time discretization. Another natural extension would be to the heat equation in higher dimension. Since the size of the problem is then bigger, we expect that the acceleration obtained with our new strategy would be even more beneficial.

Acknowledgements

This project has received funding from the Federal Ministry of Education and Research and the European High-Performance Computing Joint Undertaking (JU) under grant agreement No 955701, Time-X. The JU receives support from the European Union’s Horizon 2020 research and innovation program and Belgium, France, Germany, Switzerland.