Parallel-in-Time Multi-Level Integration of the Shallow-Water Equations on the Rotating Sphere

TL;DR

This paper introduces a multi-level parallel-in-time integration method combining PFASST with spherical harmonics for efficient simulation of shallow-water equations on the sphere, significantly reducing computation time.

Contribution

The paper presents a novel PFASST-based multi-level scheme using spherical harmonics for spatial coarsening, enabling faster parallel-in-time integration of atmospheric PDEs on the sphere.

Findings

PFASST-SH converges upon time refinement.

The coarsening strategy effectively captures high-frequency modes.

The scheme achieves multiple times speedup over serial methods.

Abstract

The modeling of atmospheric processes in the context of weather and climate simulations is an important and computationally expensive challenge. The temporal integration of the underlying PDEs requires a very large number of time steps, even when the terms accounting for the propagation of fast atmospheric waves are treated implicitly. Therefore, the use of parallel-in-time integration schemes to reduce the time-to-solution is of increasing interest, particularly in the numerical weather forecasting field. We present a multi-level parallel-in-time integration method combining the Parallel Full Approximation Scheme in Space and Time (PFASST) with a spatial discretization based on Spherical Harmonics (SH). The iterative algorithm computes multiple time steps concurrently by interweaving parallel high-order fine corrections and serial corrections performed on a coarsened problem. To do that, we design a methodology relying on the spectral basis of the SH to coarsen and interpolate the problem in space. The methods are evaluated on the shallow-water equations on the sphere using a set of tests commonly used in the atmospheric flow community. We assess the convergence of PFASST-SH upon refinement in time. We also investigate the impact of the coarsening strategy on the accuracy of the scheme, and specifically on its ability to capture the high-frequency modes accumulating in the solution. Finally, we study the computational cost of PFASST-SH to demonstrate that our scheme resolves the main features of the solution multiple times faster than the serial schemes.

Tables

Table 1. Table 1: Parameters for the serial SDC and MLSDC-SH schemes. The SDC iteration only involves one sweep on the fine level. The MLSDC-SH iteration is described in Hamon et al., ( 2019 ) , and involves one sweep on the fine level and one sweep on the coarse level.

SDC($M_f+1$, $N_{\text{S}}$)
Parameter	Description
$M_f+1$	SDC nodes on fine level
$N_S$	SDC iterations

Table 2. Table 2: Parameters for the parallel-in-time PFASST-SH scheme.

PFASST($n_{\text{ts}}$, $M_f+1$, $M_c+1$, $N_{\text{PF}}$, $\alpha$)
Parameter	Description
$n_{\text{ts}}$	Time steps solved in parallel
$M_f+1$	SDC nodes on fine level
$M_c+1$	SDC nodes on coarse level
$N_{\text{PF}}$	PFASST iterations
$\alpha$	Spatial coarsening ratio

Table 3. Table 3: Gaussian bump: L ∞ subscript 𝐿 L_{\infty} -error for R norm = 32 subscript 𝑅 norm 32 R_{\textit{norm}}=32 , time step size, and wall-clock time for serial and parallel SDC-based schemes. In the last two columns, S S obs subscript superscript 𝑆 obs 𝑆 S^{\textit{obs}}_{S} denotes the observed speedup achieved with MLSDC-SH and PFASST-SH with respect to SDC, while S ML obs subscript superscript 𝑆 obs ML S^{\textit{obs}}_{\textit{ML}} denotes the speedup achieved with PFASST-SH with respect to MLSDC-SH. The speedups for R norm = 256 subscript 𝑅 norm 256 R_{\textit{norm}}=256 can be read from Fig.  5 .

Integration	$L_{\infty }$-error	Time step	Wall-clock	Observed	Observed
scheme	$R_{\text{norm}}=32$	size (s)	time (s)	speedup $S_S^{\text{obs}}$	speedup $S_{\text{ML}}^{\text{obs}}$
SDC(3,4)	$5.2\times {10}^{-6}$	400	528	-	-
MLSDC(3,2,2,1/2)	$5.2\times {10}^{-6}$	400	354	1.5	-
PFASST(16,3,2,4,1/5)	$4.4\times {10}^{-6}$	128	120	4.1	3.0
SDC(5,8)	$4.6\times {10}^{-11}$	200	3,389	-	-
MLSDC(5,3,4,1/2)	$7.6\times {10}^{-10}$	200	2,269	1.5	-
PFASST(16,5,3,8,4/5)	$3.5\times {10}^{-11}$	100	917	3.7	2.5

Table 4. Table 4: Rossby-Haurwitz wave: L ∞ subscript 𝐿 L_{\infty} -error computed with R norm = 32 subscript 𝑅 norm 32 R_{\textit{norm}}=32 , time step size, and wall-clock time for the serial and parallel SDC-based schemes. S S obs subscript superscript 𝑆 obs 𝑆 S^{\textit{obs}}_{S} denotes the observed speedup achieved with MLSDC-SH and PFASST-SH with respect to SDC, while S ML obs subscript superscript 𝑆 obs ML S^{\textit{obs}}_{\textit{ML}} denotes the speedup achieved with PFASST-SH with respect to MLSDC-SH. The speedups for R norm = 256 subscript 𝑅 norm 256 R_{\textit{norm}}=256 can be read from Fig.  9 .

Integration	$L_{\infty }$-error	Time step	Wall-clock	Observed	Observed
scheme	$R_{\text{norm}}=32$	size (s)	time (s)	speedup $S_S^{\text{obs}}$	speedup $S_{\text{ML}}^{\text{obs}}$
SDC(3,4)	$2.0\times {10}^{-8}$	400	526	-	-
MLSDC(3,2,2,1/2)	$1.7\times {10}^{-8}$	400	366	1.4	-
PFASST(16,3,2,4,1/5)	$1.2\times {10}^{-8}$	100	166	3.2	2.2
SDC(5,8)	$4.6\times {10}^{-12}$	640	1,027	-	-
MLSDC(5,3,4,1/2)	$4.6\times {10}^{-12}$	640	710	1.4	-
PFASST(16,5,3,8,1/2)	$4.6\times {10}^{-12}$	320	198	5.2	3.6

Table 5. Table 5: Unstable barotropic wave: L ∞ subscript 𝐿 L_{\infty} -error for R norm = 64 subscript 𝑅 norm 64 R_{\textit{norm}}=64 , time step size, and wall-clock time for the serial and parallel SDC-based schemes. S S obs subscript superscript 𝑆 obs 𝑆 S^{\textit{obs}}_{S} denotes the observed speedup achieved with MLSDC-SH and PFASST-SH with respect to SDC, while S ML obs subscript superscript 𝑆 obs ML S^{\textit{obs}}_{\textit{ML}} denotes the speedup achieved with PFASST-SH with respect to MLSDC-SH.

Integration	$L_{\infty }$-error	Time step	Wall-clock	Observed	Observed
scheme	$R_{\text{norm}}=64$	size (s)	time (s)	speedup $S_S^{\text{obs}}$	speedup $S_{\text{ML}}^{\text{obs}}$
SDC(5,8)	$3.3\times {10}^{-12}$	240	15,928	-	-
MLSDC(5,3,4,4/5)	$4.2\times {10}^{-12}$	60	46,127	0.3	-
PFASST(16,5,3,16,1)	$6.5\times {10}^{-12}$	240	4,311	3.7	10.7

Equations

$$\frac{\partial \mathbf{U}}{\partial t}=\mathcal{L}(\mathbf{U})+\mathcal{N}(\mathbf{U}),$$

$$\boldsymbol{\mathcal{L}}(\boldsymbol{U})\equiv\left[\begin{array}[]{c}-\bar{\Phi}\delta+\nu\nabla^{2}\Phi^{\prime}\\[5.0pt] \nu\nabla^{2}\zeta\\[5.0pt] -\nabla^{2}\Phi+\nu\nabla^{2}\delta\end{array}\right].$$

$$\boldsymbol{\mathcal{N}}(\boldsymbol{U})=\left[\begin{array}[]{c}-\nabla\cdot(\Phi^{\prime}\boldsymbol{V})\\[5.0pt] -\nabla\cdot(\zeta+f)\boldsymbol{V}\\[0.0pt] \boldsymbol{k}\cdot\nabla\times(\zeta+f)\boldsymbol{V}-\nabla^{2}\bigg{(}\displaystyle\frac{\boldsymbol{V}\cdot\boldsymbol{V}}{2}\bigg{)}\end{array}\right],$$

$$\mathbf{V}=\mathbf{k}\times \nabla \psi +\nabla \chi .$$

$$\zeta$$

$$\delta$$

$$\xi (\lambda ,\mu )=\sum _{r=-R}^R\sum _{s=|r|}^{S(r)}{\xi }_s^r{P}_s^r(\mu )e^{ir\lambda },$$

$$\frac{\partial {\mathbf{\Theta }}_s^r}{\partial t}={\mathbf{L}}_s^r(\mathbf{\Theta })+{\mathbf{N}}_s^r(\mathbf{\Theta }),$$

$$\mathbf{\Theta }=[{\mathbf{\Theta }}_0^0,{\mathbf{\Theta }}_0^1,\dots ,{\mathbf{\Theta }}_{R-1}^R,{\mathbf{\Theta }}_R^R{]}^T.$$

$$\frac{\partial {\mathbf{\Theta }}_s^r}{\partial t}=({\mathbf{F}}_I{)}_s^r(\mathbf{\Theta })+({\mathbf{F}}_E{)}_s^r(\mathbf{\Theta }),$$

$$({\mathbf{F}}_I{)}_s^r$$

$$({\mathbf{F}}_E{)}_s^r$$

$$\frac{\partial \mathbf{\Theta }}{\partial t}(t)$$

$$\mathbf{\Theta }(t^n)$$

$$t^n\equiv t^{n,0}<t^{n,1}<\cdots <t^{n,M}=t^n+\Delta t\equiv t^{n+1}.$$

$${\mathbf{\Theta }}^{m+1,(k+1)}={\mathbf{\Theta }}^n$$

$$\displaystyle+\Delta t\sum_{j=1}^{m+1}\tilde{q}^{I}_{m+1,j}\big{[}\boldsymbol{F}_{I}\big{(}\boldsymbol{\Theta}^{j,(k+1)}\big{)}-\boldsymbol{F}_{I}\big{(}\boldsymbol{\Theta}^{j,(k)}\big{)}\big{]}$$

$$\displaystyle+\Delta t\sum_{j=0}^{M}q_{m+1,j}\boldsymbol{F}\big{(}\boldsymbol{\Theta}^{j,(k)}\big{)}.$$

$$\int_{t^{n}}^{t}\boldsymbol{F}\big{(}\tilde{\boldsymbol{\Theta}}(a)\big{)}da.$$

$$\vec{\mathbf{\Theta }}$$

$$\vec{\mathbf{F}}$$

$$\mathbf{A}(\vec{\mathbf{\Theta }})\equiv \vec{\mathbf{\Theta }}-\Delta t(\mathbf{Q}\otimes {\mathbf{I}}_K)\vec{\mathbf{F}},$$

$$\mathbf{A}(\vec{\mathbf{\Theta }})=1_{M+1}\otimes {\mathbf{\Theta }}^{n,0}.$$

$${\mathbf{A}}_c({\vec{\mathbf{\Theta }}}_c)-{\vec{\mathbf{\tau }}}_c=1_{M_c+1}\otimes {\mathbf{\Theta }}_c^{n,0},$$

$${\vec{\mathbf{\tau }}}_c\equiv {\mathbf{A}}_c({\mathbf{R}}_f^c{\vec{\mathbf{\Theta }}}_f)-{\mathbf{R}}_f^c{\mathbf{A}}_f({\vec{\mathbf{\Theta }}}_f)+{\mathbf{R}}_f^c{\vec{\mathbf{\tau }}}_f,$$

$${\mathbf{\Theta }}_c^{m+1,(k+1)}$$

$$\displaystyle+\Delta t\sum_{j=1}^{m+1}(\tilde{q}^{I}_{m+1,j})_{c}\big{[}\boldsymbol{F}_{I,c}\big{(}\boldsymbol{\Theta}^{j,(k+1)}_{c}\big{)}-\boldsymbol{F}_{I,c}\big{(}\boldsymbol{\Theta}^{j,(k)}_{c}\big{)}\big{]}$$

$$\displaystyle+\Delta t\sum_{j=0}^{M}(q_{m+1,j})_{c}\boldsymbol{F}_{c}\big{(}\boldsymbol{\Theta}^{j,(k)}_{c}\big{)}+\boldsymbol{\tau}^{m+1,(k)}_{c}.$$

$${\vec{\mathbf{\Theta }}}_c={\mathbf{R}}_f^c{\vec{\mathbf{\Theta }}}_f=({\mathbf{R}}^{\text{s}}{)}_f^c({\mathbf{R}}^{\text{t}}{)}_f^c{\vec{\mathbf{\Theta }}}_f,$$

$$({\mathbf{R}}^{\text{t}}{)}_f^c\equiv {\mathbf{\Pi }}_f^c\otimes {\mathbf{I}}_{K_f}\in {\mathbb{R}}^{(M_c+1)K_f\times (M_f+1)K_f},$$

$$({\mathbf{\Pi }}_f^c{)}_{ij}=L_f^{j-1}(t_c^{i-1}),$$

$$({\mathbf{R}}^s{)}_f^c\equiv {\mathbf{I}}_{M_c+1}\otimes {\mathbf{D}}_f^c\in {\mathbb{R}}^{(M_c+1)K_c\times (M_c+1)K_f}.$$

$${\vec{\mathbf{\Theta }}}_f\equiv {\mathbf{P}}_c^f{\vec{\mathbf{\Theta }}}_c=({\mathbf{P}}^t{)}_c^f({\mathbf{P}}^s{)}_c^f{\vec{\mathbf{\Theta }}}_c.$$

$$(\boldsymbol{P}^{s})^{f}_{c}\equiv\big{(}(\boldsymbol{R}^{s})^{c}_{f}\big{)}^{T}\in\mathbb{R}^{(M_{c}+1)K_{f}\times(M_{c}+1)K_{c}}.$$

Full text

We present a multi-level parallel-in-time integration method combining the Parallel Full Approximation Scheme in Space and Time (PFASST) with a spatial discretization based on Spherical Harmonics (SH). The iterative algorithm computes multiple time steps concurrently by interweaving parallel high-order fine corrections and serial corrections performed on a coarsened problem. To do that, we design a methodology relying on the spectral basis of the SH to coarsen and interpolate the problem in space.

The methods are evaluated on the shallow-water equations on the sphere using a set of tests commonly used in the atmospheric flow community. We assess the convergence of PFASST-SH upon refinement in time. We also investigate the impact of the coarsening strategy on the accuracy of the scheme, and specifically on its ability to capture the high-frequency modes accumulating in the solution. Finally, we study the computational cost of PFASST-SH to demonstrate that our scheme resolves the main features of the solution multiple times faster than the serial schemes.

1 Introduction

The accurate simulation of atmospheric flows over long periods of time is one of the critical components in the fields of numerical weather prediction and climate modeling. Despite a rich history of numerical methods in these fields, the development of more accurate and efficient temporal integration methods for weather and climate simulations is still an ongoing challenge. One difficulty is the presence of a wide range of time scales in the equations, including the propagation of fast atmospheric waves, which imposes a severe stability restriction on the time step size of fully explicit integration schemes. Implicit-explicit schemes partly overcome this limitation by treating the stiff terms implicitly, but still involve a very large number of time steps for high-resolution long-range simulations (up to a thousand years for paleoclimate studies). Therefore, parallel-in-time methods are an attractive approach to reduce the time-to-solution by accessing an additional axis of parallelism in the temporal direction.

In this work, we study the performance and accuracy of an implicit-explicit, iterative, multi-level, parallel-in-time integration scheme based on the Parallel Full Approximation Scheme in Space and Time (PFASST, Emmett and Minion,, 2012) in the context of atmospheric flows. The method is applied to the Shallow-Water Equations (SWE) on the rotating sphere, a common two-dimensional proxy that incorporates the horizontal features of the full three-dimensional atmospheric equations. In the classification of Burrage, (1997), the parallelization strategy in PFASST relies on parallelization across the steps, which consists in solving multiple time steps concurrently on different processing units to accelerate the simulations. Other examples of this type of parallelism include Parareal (Lions et al.,, 2001), the Parallel Implicit Time-integration Algorithm (PITA, Farhat and Chandesris,, 2003) and the MultiGrid Reduction in Time (MGRIT, Falgout et al.,, 2014). This approach differs from other classes of methods not considered in this work, namely, parallelization across the method, in which intermediate stage values are computed in parallel (e.g., Butcher,, 1997), and parallelization across the problem, in which the full problem is split into subproblems coupled iteratively (e.g., Gander,, 1999). Parallel-in-time methods based on exponential integration have also been proposed (Gander and Güttel,, 2013; Haut and Wingate,, 2014). They have been applied to the SWE on the plane (Schreiber et al.,, 2017) and on the rotating sphere (Schreiber and Loft,, 2018) using a REXI (Rational Approximation of Exponential Integrators) approach to compute the exponential. Unlike PFASST, these schemes integrate one time step at a time, but perform the integration of a given time step by summing the solutions of fully decoupled systems solved in parallel. This yields a highly accurate and stable integration of the linear(ized) terms that requires an efficient solver for complex-valued linear systems. For nonlinear problems, exponential integration must be combined with other time integrators handling the nonlinear terms using a splitting scheme (Schreiber et al.,, 2019).

The structure of PFASST is similar to that of Parareal as both algorithms combine parallel, relatively expensive updates based on approximate initial conditions and serial inexpensive updates propagating the new initial conditions on the time interval. However, the PFASST solution updates are computed iteratively with Spectral Deferred Corrections (SDC, Dutt et al.,, 2000) instead of direct ODE integrators such as Runge-Kutta schemes in standard Parareal. This makes it possible to amortize the cost of the fine-level SDC updates over multiple PFASST iterations (Minion,, 2011), thereby reducing the computational cost of the fine propagator. The flexibility inherent in the SDC approach also allows for the straightforward construction of implicit-explicit (IMEX) methods of very high order, which has not yet been attained with methods like linear-multistep or Runge-Kutta methods (Minion,, 2003). IMEX schemes are preferred in the current setting over explicit or fully implicit methods.

In PFASST the SDC corrections are interwoven on a hierarchy of levels representing coarse space-time approximations of the fine problem under consideration. This iterative procedure allowing for spatial coarsening results in a larger theoretical parallel efficiency than with the Parareal algorithm. An efficient implementation of the PFASST algorithm has been proposed (Emmett and Minion,, 2014), and the parallel efficiency of the scheme has been demonstrated for some applications, including the N-body problem (Speck et al.,, 2012) and the heat equation (Speck et al.,, 2014). Recently, in Bolten et al., (2017, 2018), the PFASST algorithm has been cast as a multigrid method and its convergence has been studied for synthetic diffusion-dominated and advection-dominated problems. PFASST has also been used in the context of optimal control problems (Götschel and Minion,, 2019).

The space-time coarsening and interpolation strategy is one of the critical determinants of the accuracy and efficiency of the PFASST algorithm. We consider the approach proposed in the serial Multi-Level Spectral Deferred Corrections (MLSDC) scheme in Hamon et al., (2019). It takes advantage of the structure of the spatial discretization based on the Spherical Harmonics (SH) transform, which plays a key role in our parallel-in-time scheme, thereafter referred to as PFASST-SH. Specifically, we have shown in previous work that the SH basis can be truncated, or padded with zeros, to construct consistent spatial restriction and interpolation operators between levels. As in Hamon et al., (2019), we exploit the properties of the SH transform to efficiently solve the implicit systems at each temporal node. This study is relevant for operational numerical weather prediction systems as the SH transform is used in the Integrated Forecast System (IFS) at the European Centre for Medium-Range Weather Forecasts (ECMWF, Wedi et al.,, 2013) and the Global Spectral Model (GSM) at the Japan Meteorological Agency (JMA, Kanamitsu et al.,, 1983).

We evaluate the properties of the proposed PFASST-SH algorithm using a well-known suite of nonlinear shallow-water test cases that are representative of the horizontal features of atmospheric flows (Williamson et al.,, 1992; Galewsky et al.,, 2004). These wave-propagation examples are challenging for the integration scheme as parallel-in-time methods are known to suffer from convergence difficulties on advection-dominated problems (e.g., Ruprecht,, 2018). In addition, the nonlinear problems considered in this work are characterized by the progressive amplification over time of high-frequency modes not captured by the coarse SDC updates. This can undermine the accuracy of the parallel approach on these modes when spatial coarsening is aggressive. However, the use of damping methods (second-order, fourth-order, or sixth-order diffusion, spectral viscosity) to resolve the spectral blocking problem (Gelb and Gleeson,, 2001) brings into question the relevance of accuracy in the highest frequencies. We investigate this question by using different measures of error on the SH discretization in the experiments. In the numerical examples, we start by evaluating the convergence rate of PFASST-SH as the time step size is decreased. This confirms that spatial coarsening is a key determinant of the accuracy of the parallel scheme, and demonstrates that PFASST-SH can capture the small-scale features of the solution when the spatial coarsening is mild. Then, we measure the computational cost of PFASST-SH to show that it yields a significant reduction in the time-to-solution compared to serial single-level and multi-level SDC schemes for the three numerical examples considered here. A strong scalability test performed on the benchmark proposed in Galewsky et al., (2004) concludes our numerical study of the properties of PFASST-SH.

In Sections 2, 3, and 4, we review some features of the problem – namely, the governing equations, SH transform, and IMEX temporal splitting – already discussed in Hamon et al., (2019) to ensure the completeness of this paper. We present the parallel-in-time numerical scheme based on PFASST-SH in Section 5. We evaluate its performance and accuracy on SWE test cases in Section 6.

2 Mathematical model

In this work, we study the properties of the parallel integration scheme applied to the Shallow-Water Equations (SWE) on the rotating sphere. This simplified mathematical model captures the horizontal features of realistic atmospheric flows and allows us to evaluate the main properties of our numerical scheme – namely, accuracy, performance, and robustness – on a set of well-defined test cases (Williamson et al.,, 1992; Galewsky et al.,, 2004). To overcome the singularity in the velocity field at the poles, we use the vorticity-divergence formulation in physical space (Bourke,, 1972; Hack and Jakob,, 1992). The prognostic variables $\boldsymbol{U}=[\Phi,\,\zeta,\,\delta]^{T}$ are respectively the geopotential, $\Phi$, the vorticity, $\zeta$, and the divergence, $\delta$. The system of governing partial differential equations reads

[TABLE]

where the first group in the right-hand side of (1) contains the stiff terms involved in the linear wave motion induced by gravitational forces and also includes the diffusion term

[TABLE]

The average geopotential, $\bar{\Phi}=g\bar{h}$, is written as the product of the gravitational acceleration by the average height, and $\Phi^{\prime}$ is defined as $\Phi^{\prime}=\Phi-\bar{\Phi}$. The diffusion coefficient is denoted by $\nu$. Although spectral viscosity (Gelb and Gleeson,, 2001) is now used in standard codes like IFS, we employ a second-order diffusion term in the governing equations to stabilize the flow dynamics and reduce the errors caused by under-resolved nonlinear interaction modes. Second-order diffusion remains acceptable here since we are primarily interested in showing that PFASST-SH can resolve the low-frequency and mid-frequency features of the solution significantly faster than serial SDC schemes. This choice also allows for a fair comparison with the results obtained previously on the same test cases with SDC methods (Jia et al.,, 2013; Hamon et al.,, 2019). The second group in the right-hand side of (1) contains all the relatively less stiff, nonlinear terms present in the governing equations

[TABLE]

where $\boldsymbol{k}=[0,0,1]^{T}$. The horizontal velocity vector is $\boldsymbol{V}\equiv\boldsymbol{i}u+\boldsymbol{j}v$, where $\boldsymbol{i}$ and $\boldsymbol{j}$ are the unit vectors in the eastward and northward directions, respectively. The Coriolis force is represented by $f=2\Omega\sin\phi$, where $\Omega$ is the angular rate of rotation, and $\phi$ is the latitude. In Section 4, we use this linear-nonlinear decomposition to define an IMEX temporal splitting that treats the stiff terms implicitly and the less stiff terms explicitly.

To express the velocities as a function of the prognostic variables, $\zeta$ and $\delta$, we first use the Helmholtz theorem which relates $\boldsymbol{V}$ to a scalar stream function, $\psi$, and a scalar velocity potential, $\chi$,

[TABLE]

Using the identities

[TABLE]

the application of the curl and divergence operators to (4) yields $\zeta=\nabla^{2}\psi$ and $\delta=\nabla^{2}\chi$. The Laplace operators can be efficiently inverted using the SH transform to compute the stream function, $\psi$, and the velocity potential, $\chi$, as a function of $\zeta$ and $\delta$. This is discussed in the next section, along with the spatial discretization of $\boldsymbol{\mathcal{L}}$ and $\boldsymbol{\mathcal{N}}$.

3 Spectral transform method

Here, we review the methodology based on the global SH transform used to discretize the governing equations in space. The parallel multi-level scheme proposed in this work heavily relies on the structure of this spatial discretization for the construction of a hierarchy of space-time levels. In the SH scheme, the representation of a function of longitude $\lambda$ and Gaussian latitude $\mu\equiv\sin(\phi)$, $\xi(\lambda,\mu)$, consists of a sum of spherical harmonic basis functions, $P^{r}_{s}(\mu)e^{ir\lambda}$, weighted by the spectral coefficients $\xi^{r}_{s}$,

[TABLE]

where the indices $r$ and $s$ refer to the zonal and total wavenumbers, respectively. In (7), $P^{r}_{s}$ is the normalized associated Legendre polynomial. Without loss of generality, we use a triangular truncation with $S(r)=R$, where $R$ and $S(r)$ are the chosen truncation limits of zonal and total wavenumbers, respectively. We note that truncating the modes corresponding to high-frequency wavenumbers – i.e., reducing $R$ in (7) – yields a consistent coarse representation of the discrete problem. This will be exploited in Section 5 in the integration scheme.

The transformation from physical to spectral space is achieved in two steps. The first step consists in taking the discrete Fourier transform of $\xi(\lambda,\mu)$ in longitude – i.e., over $\lambda$ –, followed by a second step based on the application of the discrete Legendre transformation in latitude. This two-step global transform is applied to (1) to obtain a system of coupled ordinary differential equations involving the prognostic variables in spectral space, $\boldsymbol{\Theta}^{r}_{s}=[\Phi^{r}_{s},\,\zeta^{r}_{s},\,\delta^{r}_{s}]$. Noting that, due to the symmetry of the spectral coefficients, we only have to consider the indices $r\geq 0$, we obtain for $r\in\{0,\dots,R\}$ and $s\in\{r,\dots,R\}$,

[TABLE]

where $\boldsymbol{L}^{r}_{s}$, and $\boldsymbol{N}^{r}_{s}$ are the discrete, spectral representations of the operators defined in (2), and (3). The state variable in spectral space, $\boldsymbol{\Theta}$, is defined as a vector of size $K$ as follows

[TABLE]

An efficient implementation of the global SH transform is described in Temperton, (1991); Rivier et al., (2002); Schaeffer, (2013). The latter reference served as a basis for the developments presented here.

4 Temporal splitting

In atmospheric modeling, the propagation of fast waves – e.g., sound or gravity waves – in the system often imposes a severe stability restriction on the time step size of fully explicit schemes. Fully implicit schemes overcome this stability constraint but require solving expensive global nonlinear systems (Evans et al.,, 2010; Jia et al.,, 2013; Lott et al.,, 2015).

Instead, implicit-explicit (IMEX) schemes only evaluate implicitly the stiff terms involved in the propagation of the fast waves, while less stiff terms are updated explicitly. This strategy reduces the cost of the implicit solves and allows for relatively large stable time steps. In non-hydrostatic atmospheric modeling, dimensional splitting is a commonly used IMEX strategy that only treats implicitly the terms involved in the (fast) vertical dynamics (Ullrich and Jablonowski,, 2012; Weller et al.,, 2013; Giraldo et al.,, 2013; Lock et al.,, 2014; Gardner et al.,, 2018). In Smolarkiewicz et al., (2014), an IMEX scheme for the compressible Euler equations is obtained by treating all terms pertinent to wave motions (namely acoustic, gravity, and rotational) implicitly, but advection explicitly. Alternatively, the approach of Robert et al., (1972); Giraldo, (2005) consists in linearizing the governing PDEs and treating the linearized piece implicitly. The term treated explicitly is obtained by subtracting the linearized part from the nonlinear system.

In this work, we directly discretize the fast linear terms on the right-hand side of (8) implicitly, while the other terms are evaluated explicitly. Specifically, our integration scheme relies on the following splitting, written below in semi-discrete form for $r\in\{0,\dots,R\}$ and $s\in\{r,\dots,R\}$,

[TABLE]

where

[TABLE]

That is, the implicit right-hand side, $(\boldsymbol{F}_{I})^{r}_{s}$, contains the terms representing linear wave motion induced by gravitational forces and the diffusion term. The explicit right-hand side, $(\boldsymbol{F}_{E})^{r}_{s}$, contains a linear harmonic oscillator and the nonlinear terms. As explained in Section 5.4, in the context of the SH spatial discretization, the implicit solve necessary in the IMEX time-stepping approach is inexpensive compared to the cost of the explicit evaluation of nonlinear terms. Hence, the IMEX approach greatly reduces the computational cost per step compared to a fully implicit method by circumventing the need for a global nonlinear solver. The IMEX approach also allows much larger stable time steps than a fully explicit method with little additional computational cost. Next, we describe the parallel integration scheme that is used to advance the semi-discrete system (10) in time.

5 Parallel-in-time integration

We review the fundamentals of the Parallel Full Approximation Scheme in Space and Time (PFASST) algorithm detailed in Emmett and Minion, (2012) after earlier work by Minion, (2011). In the remainder of the paper, an SDC sweep refers to the computation of a correction with respect to a previous iteration. The multi-level scheme consists in coupling iteratively serial coarse SDC sweeps propagating the initial conditions with parallel fine SDC sweeps achieving high-order accuracy at each time step. We start from two building blocks that are key to the accuracy and efficiency of the scheme, namely the SDC sweeps in Section 5.1.1 and the Full Approximation Scheme (FAS) in Section 5.1.2. Then we detail the steps of PFASST and its computational cost in Sections 5.2 to 5.5. We consider a system of coupled ODEs in the generic form

[TABLE]

where $\boldsymbol{F}_{I}$ and $\boldsymbol{F}_{E}$ are the implicit and explicit right-hand sides, respectively, with $\boldsymbol{F}=\boldsymbol{F}_{I}+\boldsymbol{F}_{E}$ and $\boldsymbol{\Theta}(t)$ is the state variable at time $t$.

5.1 Ingredients of the time integration scheme

5.1.1 IMplicit-EXplicit Spectral Deferred Correction (IMEX SDC)

SDC methods have been presented in Dutt et al., (2000) and then generalized to methods with different temporal splittings in Minion, (2003); Bourlioux et al., (2003); Layton and Minion, (2004). In the context of fast-wave slow-wave problems, the properties of IMEX SDC schemes have been studied in Ruprecht and Speck, (2016). In SDC methods, the interval $[t^{n},t^{n+1}]$ is decomposed into $M$ subintervals using $M+1$ Gauss-Lobatto temporal nodes, such that

[TABLE]

In the remainder of this paper, we use the shorthand notation $t^{m}=t^{n,m}$. We denote by $\boldsymbol{\Theta}^{m+1,(k+1)}$ the approximate solution at node $m+1$ and at sweep $(k+1)$. The SDC scheme applied to the implicit-explicit temporal splitting (13) iteratively improves the accuracy of the approximation based on the discrete correction equation

[TABLE]

The coefficients $q_{m+1,j}$ are chosen to be the Gauss-Lobatto quadrature points, such that the term on the third line of (16) is a high-order approximation of the integral

[TABLE]

Finally, the coefficients $\tilde{q}^{E}_{m+1,j}$ correspond to forward-Euler time stepping, while the weights $\tilde{q}^{I}_{m+1,j}$ are chosen to be the coefficients of the upper triangular matrix in the LU decomposition of $\boldsymbol{Q}^{T}$, where $\boldsymbol{Q}=\{q_{ij}\}\in\mathbb{R}^{(M+1)\times(M+1)}$. We refer to Weiser, (2015) for a proof that this choice leads to fast convergence of the iterative process to the fixed-point solution for stiff problems, and to Hamon et al., (2018) for numerical examples illustrating the improved convergence. Each pass of the discrete version of the update equation (16), referred to as sweep, increases the formal order of accuracy by one until the order of accuracy of the quadrature applied to the third integral is reached (Christlieb et al.,, 2009). We mention here that our implementation splits the implicit step of (16) into two substeps, with a substep for the physical linear terms followed by a substep for the artificial diffusion terms. This additional splitting is described in Bourlioux et al., (2003) and will not be discussed in the present paper for brevity.

Using the compact notation of Bolten et al., (2017), we introduce the space-time vectors $\vec{\boldsymbol{\Theta}}\in\mathbb{C}^{(M+1)K}$ and $\vec{\boldsymbol{F}}\in\mathbb{C}^{(M+1)K}$ are such that

[TABLE]

We also define the operator

[TABLE]

where $\otimes$ denotes the Kronecker product and $\boldsymbol{I}_{K}\in\mathbb{R}^{K\times K}$ is the identity matrix. $\boldsymbol{1}_{M+1}\in\mathbb{R}^{M+1}$ is a vector of ones. Using these definitions, the integration scheme (16) can be written as an iterative solution method for the collocation problem defined by

[TABLE]

In the following sections, we describe how the SDC sweeps are applied iteratively on coarse and fine representations of the problem to obtain a parallel-in-time integration scheme.

5.1.2 Full Approximation Scheme (FAS)

The PFASST algorithm requires the definition of a coarse approximation in space and in time of the discrete problem of interest. We refer to the former as the coarse level $\ell=c$, and to the latter as the fine level $\ell=f$. The space-time coarsening methodology used to construct the coarse space is reviewed in Section 5.3. On these two levels, the PFASST algorithm computes multiple time steps concurrently by interweaving parallel high-order SDC sweeps on the fine level and serial SDC sweeps performed on the coarse level to propagate the updated initial conditions between time steps (see Section 5.2). We do not consider the case in which PFASST involves more than one coarse level. To ensure that the fine and coarse levels are properly coupled, we employ the Full Approximation Scheme (FAS, Brandt,, 1977), described next. In FAS, the discrete equations solved on the coarse level are modified with the introduction of a correction term, $\vec{\boldsymbol{\tau}}_{c}$. The coarse collocation problem with the correction term becomes

[TABLE]

where $\boldsymbol{A}_{c}$ is the approximation of $\boldsymbol{A}$ at the coarse level. The FAS correction term at the coarse level is defined as

[TABLE]

with, for the two-level case, $\vec{\boldsymbol{\tau}}_{f}=\boldsymbol{0}$ on the fine level. In (22) and (23), $\vec{\boldsymbol{\Theta}}_{\ell}\in\mathbb{C}^{(M_{\ell}+1)K_{\ell}}$ and $\vec{\boldsymbol{F}}_{\ell}\in\mathbb{C}^{(M_{\ell}+1)K_{\ell}}$ are respectively the space-time vector and right-hand side at level $\ell$. $K_{\ell}$ represents the total number of spectral coefficients in (9) on level $\ell$. The matrix $\boldsymbol{R}^{c}_{f}\in\mathbb{R}^{(M_{c}+1)K_{c}\times(M_{f}+1)K_{f}}$ is the linear restriction operator from the fine level to the coarse level that consists in truncating the spectral coefficients corresponding to the high-frequency modes in the SH basis. With this modification of the coarse equations, we note that the restriction of the fine solution, $\boldsymbol{R}^{c}_{f}\vec{\boldsymbol{\Theta}}_{f}$, is a solution of the coarse problem. On the coarse problem (22), the modified SDC update for temporal node $m+1$ at sweep $(k+1)$ is

[TABLE]

5.2 Parallel Full Approximation Scheme in Space and Time (PFASST) algorithm

We present the steps of the PFASST algorithm presented in Emmett and Minion, (2012), with the improved communication pattern of Emmett and Minion, (2014). The presentation of the algorithm is done for the two-level case that we exclusively consider in this work. We consider that $n_{\textit{ts}}$ processors are available to solve a block of $n_{\textit{ts}}$ time steps in parallel and we assume that the $n^{\text{th}}$ time step is assigned to processor, or set of computing resources, $\mathcal{P}_{n}$. We denote by $\boldsymbol{\Theta}^{n,m,(k)}$ the approximate solution at time step $n$, at SDC node index $m$, and at PFASST iteration $(k)$. In compact form, $\vec{\boldsymbol{\Theta}}^{n,(k)}$ is the space-time vector containing the approximate solution at time step $n$, at all temporal nodes, and PFASST iteration $(k)$. This algorithm, illustrated in Fig. 1 on a block of three time steps solved in parallel, starts with a prediction step, followed by a sequence of iterations to correct the prediction.

5.2.1 Prediction step

When the algorithm starts, the only known initial condition – that is, the solution at SDC node $m=0$ on a given time step – is at the first time step of the block, denoted by $\Delta t^{0}$, and owned by processor $\mathcal{P}_{0}$. The prediction procedure aims at using serial coarse SDC sweeps to generate an approximate fine solution for all the SDC nodes of the time steps of the block.

At the beginning of the prediction step, processor $\mathcal{P}_{0}$ sends its initial condition to the other processors. Then, each processor restricts in space the fine initial condition received from $\mathcal{P}_{0}$ using the methodology described in Section 5.3. After that, processor $\mathcal{P}_{n}$ ($0\leq n\leq n_{\textit{ts}}-1$) performs $n+1$ relatively inexpensive coarse sweeps with an updated initial condition (SDC node $m=0$) received from processor $\mathcal{P}_{n-1}$ before each sweep. A coarse sweep consists in applying (24) at each coarse node of its time interval. Finally, each processor interpolates the solution resulting from the coarse sweeps to the fine level to generate the fine prediction for the PFASST iterations. The fine implicit and explicit right-hand sides are then reevaluated using this fine prediction. We refer the reader to Emmett and Minion, (2012) for a detailed description of the prediction step.

This prediction procedure does not require fine sweeps and is therefore inexpensive compared to a PFASST iteration. We note that the prediction step introduces a load imbalance in the algorithm since the processors do not perform the same number of coarse sweeps. This imbalance, illustrated in Fig. 1, remains small compared to the total cost of the algorithm. Without any predictor, the load imbalance would still appear in the first iteration after the coarse sweep due to its serial nature (see below). The idea of the predictor is to use this fact to improve the initial guess for each time step.

5.2.2 PFASST iteration

After the prediction step, a sequence of iterations is used to improve the quality of the predicted solution and achieve high-order accuracy on the time intervals with expensive, but parallel, fine SDC sweeps. The PFASST iteration is illustrated in Fig. 1 and written in pseudocode in Algorithm 1. We use the PFASST algorithm with a fixed number of iterations. The behavior of the algorithm with a variable number of iterations and a convergence check to stop the iterations will be studied in future work.

We describe below the steps taken by processor $\mathcal{P}_{n}$ ($0\leq n\leq n_{\textit{ts}}-1$) during one PFASST iteration. The implementation of the algorithm on a high-performance computing platform as well as the optimization of the communication pattern is discussed in Emmett and Minion, (2014).

1. A)

The iteration starts with a sweep on the fine level. This procedure consists in applying the discrete correction (16) at each fine SDC node and is the most computationally expensive step in the iteration. For the last iteration of a time step, we skip Steps B, C, and D and we stop Algorithm 1 here. This approach is more robust than stopping the iterations after the coarse calculations and interpolations (Step D) when the coarsening strategy is very aggressive. If this is not the last iteration, the updated fine value resulting from this sweep at the last SDC node ($m=M_{f}$) is sent to the next processor, $\mathcal{P}_{n+1}$. 2. B)

The fine approximate solution is restricted in time and in space to the coarse level. Then, the coarse implicit and explicit right-hand sides are reevaluated at all the coarse SDC nodes using this restricted approximated solution. We also compute the FAS correction. These quantities will be used to compute the coarse correction during the sweep of Step C. Importantly, we also save the solution and the right-hand sides at all the SDC nodes to interpolate in space and in time the coarse change in these quantities to the fine level at Step D. 3. C)

Processor $\mathcal{P}_{n}$ receives the updated coarse initial condition from $\mathcal{P}_{n-1}$. To reflect this change in initial condition, the coarse right-hand sides at the first SDC node of this time interval are reevaluated. At this point, we are ready to perform a coarse sweep. This sweep yields a new coarse value at the last SDC node of time step $n$ that is sent to $\mathcal{P}_{n+1}$ to conclude Step C. 4. D)

We interpolate the coarse solution update in space and in time to the fine level. Then, instead of reevaluating the right-hand sides on the fine level, we interpolate the coarse change in right-hand sides to the fine level using the data saved at Step B. The same choice was made in Hamon et al., (2019). Our numerical tests show that avoiding these fine function evaluations reduces the computational cost without significantly undermining the accuracy of the scheme. After the interpolation, the new fine initial condition is received from $\mathcal{P}_{n-1}$. For consistency, we conclude the PFASST iteration by the application of the interpolated coarse increment computed for SDC node $m=0$ to this new initial condition. This initial condition will be used at the next iteration performed by $\mathcal{P}_{n}$.

For the special case of $n=0$, the processor in charge of the first time step, $\mathcal{P}_{0}$, follows Algorithm 1 but skips the ReceiveInitialCondition steps. Next, we describe the restriction and interpolation operators employed in Steps B and D, respectively.

5.3 Construction of the coarse level

The coarsening strategy used to construct the space-time level $\ell=c$ has a strong impact on the performance and accuracy of the PFASST algorithm. The goal is to obtain a coarse problem on which the SDC sweeps are inexpensive but can still capture most of the features of the fine solution. This is a challenge for nonlinear shallow-water problems whose temporal evolution is characterized by the progressive amplification of the modes corresponding to high-frequency wavenumbers that cannot be easily resolved by the coarse sweeps. Therefore, the presence of these high-frequency modes in the fine solution will inherently limit the spatial coarsening that can be applied to this class of problems.

Here, we adopt the strategy proposed in Hamon et al., (2019) to construct the coarse level and transfer data between levels. The spatial coarsening strategy is performed entirely in spectral space and acts directly on the spectral basis. This approach is based on mode truncation and zero-padding and does not introduce spurious modes in the approximated solution. The restriction of the approximate solution and right-hand sides from the fine level to the coarse level is performed in two steps. This operation, $\boldsymbol{R}^{c}_{f}$, can be decomposed into a restriction in time followed by a restriction in space as follows

[TABLE]

where $(\boldsymbol{R}^{\textit{t}})^{c}_{f}$ and $(\boldsymbol{R}^{\textit{s}})^{c}_{f}$ represents the temporal and spatial operator, respectively. As explained above, $\vec{\boldsymbol{\Theta}}_{\ell}\in\mathbb{C}^{(M_{\ell}+1)K_{\ell}}$ is the space-time vector storing the state of the system at level $\ell$, where $K_{\ell}$ denotes the number of spectral coefficients used in (9) at level $\ell$ and $M_{\ell}+1$ is the number of SDC nodes at level $\ell$.

In (25), the restriction operator in time is defined using the Kronecker product as

[TABLE]

where $\boldsymbol{I}_{K_{f}}\in\mathbb{R}^{K_{f}\times K_{f}}$ is the identity matrix, and $\boldsymbol{\Pi}^{c}_{f}\in\mathbb{R}^{(M_{c}+1)\times(M_{f}+1)}$ is the rectangle matrix employed to interpolate a scalar function from the fine temporal discretization to the coarse temporal discretization. This matrix is defined using the Lagrange polynomials $L^{j}_{f}$ on the fine temporal discretization as follows

[TABLE]

using the SDC node $i-1$ at the coarse level, denoted by $t^{i-1}_{c}$. In the numerical examples, we will limit the analysis to the special cases of two, three, and five Gauss-Lobatto nodes. For this choice, applying this restriction operator in time is equivalent to performing pointwise injection. The restriction in space requires choosing the number of spectral coefficients that can be represented on the coarse level, denoted by $K_{c}$. Then, we truncate the spectral representation of the primary variables (9) in the SH transform to remove $K_{f}-K_{c}$ spectral coefficients corresponding to the high-frequency spatial features of the approximate solution. This is achieved by applying the matrix

[TABLE]

In (28), $\boldsymbol{D}^{c}_{f}\in\mathbb{R}^{K_{c}\times K_{f}}$ is a rectangle truncation matrix defined as such that $(\boldsymbol{D}^{c}_{f})_{ii}=1$, and $(\boldsymbol{D}^{c}_{f})_{ij}=0$ whenever $i\neq j$.

The interpolation procedure used to transfer the approximate solution from the coarse level to the fine level is also defined in two steps. It begins with the application of the interpolation operator in space, $(\boldsymbol{P}^{s})^{f}_{c}$ and continues with the application of the interpolation operator in time, $(\boldsymbol{P}^{t})^{f}_{c}$,

[TABLE]

The interpolation operator in space is defined as the transpose of the restriction operator in space:

[TABLE]

These operations are realized in an efficient way by padding the spectral representation of the primary variables at the coarse level with $K_{f}-K_{c}$ zeros. Finally, the interpolation operator in time is analogous to (26) and reads

[TABLE]

where the rectangle interpolation matrix $\boldsymbol{\Pi}^{f}_{c}$ is constructed with the Lagrange polynomials $L^{j}_{c}$ on the coarse temporal discretization. For two, three, and five Gauss-Lobatto nodes, (31) amounts to performing pointwise injection at the fine nodes that correspond to the coarse nodes, and then polynomial interpolation to compute the solution at the remaining fine nodes. In Section 6, we will show that the accuracy of PFASST-SH is heavily dependent on the choice of the spatial coarsening ratio in the coarsening step, since this parameter determines the range of spectral modes that can be captured on the coarse level.

5.4 Solver for the implicit systems

The parallel-in-time integration scheme entails solving implicit systems at each SDC node such that $m>0$. They are in the form

[TABLE]

In (32), $\boldsymbol{b}$ is obtained from (16) or (24). For simplicity, we have dropped the subscripts denoting the space-time levels. The structure of the linear right-hand side in (11) is determined by the spatial discretization and temporal splitting described above. As in Hamon et al., (2019), this structure is exploited to circumvent the need for a linear solver and efficiently get the updated solution in spectral space via local updates only. We refer the reader to Hamon et al., (2019) for a presentation of this solution strategy.

5.5 Computational cost of PFASST-SH

Next we compare the theoretical computational cost of PFASST-SH to those of the serial SDC and MLSDC-SH schemes (Hamon et al.,, 2019). The PFASST-SH integration scheme with a block of $n_{\textit{ts}}$ time steps solved in parallel – which is assumed here to correspond to the number of processors –, $M_{f}+1$ temporal nodes on the fine level, $M_{c}+1$ temporal nodes on the coarse level, $N_{\textit{PF}}$ iterations, and a spatial coarsening ratio between the number of coarse zonal wavenumbers and the number of fine zonal wavenumbers, $\alpha=R_{c}/R_{f}$, is denoted by PFASST($n_{\textit{ts}}$, $M_{f}+1$, $M_{c}+1$, $N_{\textit{PF}}$, $\alpha$). We refer to the single-level SDC scheme with $M_{f}+1$ temporal nodes and $N_{S}$ fine sweeps as SDC($M_{f}+1$, $N_{S}$). We denote by MLSDC($M_{f}+1$, $M_{c}+1$, $N_{ML}$, $\alpha$) the MLSDC-SH scheme with $M_{f}+1$ nodes on the fine level, $M_{c}+1$ nodes of the coarse level, $N_{\textit{ML}}$ iterations, and a spatial coarsening ratio of $\alpha$. This notation for the PFASST-SH, SDC, and MLSDC-SH schemes is summarized in Tables 1 and 2.

To approximate the computational cost of the integration schemes, we consider a block of $n_{\textit{ts}}$ time steps of individual size $\Delta t_{\textit{PF}}$ solved in parallel with PFASST-SH. $C^{s}_{\ell}$ denotes the cost of solving the implicit system (32) on level $\ell$. $C^{\textit{fi}}_{\ell}$ (respectively, $C^{\textit{fe}}_{\ell}$) is the cost of evaluating the implicit (respectively, explicit) right-hand side on level $\ell$. We decompose the cost of PFASST($n_{\textit{ts}}$, $M_{f}+1$, $M_{c}+1$, $N_{\textit{PF}}$, $\alpha$) into the cost of the prediction, and the cost of the PFASST-SH iterations. The cost of the prediction is

[TABLE]

where the term in (33a) represents the cost of re-evaluating the coarse right-hand sides at all the temporal nodes after the restriction in Step B. The term in (33b) corresponds to the cost of re-evaluating the right-hand side at the first node once the new initial condition has been received. Finally, the cost of the $n_{\textit{ts}}$ serial coarse sweeps performed during the prediction is in (33c).

The cost of the PFASST iteration detailed in Algorithm 1 reads

[TABLE]

The term in (34a) represents the cost of the fine sweep (Step A). The term in (34b) is the cost of re-evaluating the coarse right-hand sides at all nodes except the first one after the restriction (Step B). The term in (34c) represents the cost of re-evaluating the coarse right-hand sides at the first node after receiving the new initial condition in Step C. The cost of the coarse sweep is in (34d), and the cost of re-evaluating the fine right-hand side after receiving the new initial condition is in (34e).

Assuming that $N_{\textit{PF}}$ iterations are performed and denoting the communication costs of the full algorithm – i.e., prediction and iterations – by $C^{\textit{comm}}$, we obtain

[TABLE]

This derivation assumes that the cost of computing the FAS correction at the end of Step B (a linear combination of function values already accounted for) is negligible. The computational cost of the serial SDC and MLSDC-SH are considered in Hamon et al., (2019). Briefly, using the same block of $n_{\textit{ts}}$ time steps of individual size $\Delta t_{\textit{PF}}=\Delta t_{\textit{SDC}}$ the serial SDC($M_{f}+1$, $N_{S}$) has a computational cost of

[TABLE]

where the term in (36a) represents the cost of recomputing the right-hand sides at the beginning of each time step – not accounted for in Hamon et al., (2019) because it is the same for SDC and MLSDC-SH – and the term in (36b) is the cost of the sweeps. To compare (35) and (36b), we express the cost of the coarse operators as functions of the fine grid quantities and the spatial coarsening ratio. To this end, we assume that the cost of the operators is proportional to the number of spectral coefficients in (9), denoted by $K_{\ell}$. We know that $K_{\ell}=3R_{\ell}(R_{\ell}+1)/2$, where $R_{\ell}$ is the number of zonal wavenumbers on level $\ell$, since there are three primary variables individually represented with $R_{\ell}(R_{\ell}+1)/2$ spectral coefficients. This yields

[TABLE]

where we have used the definition of the spatial coarsening ratio $\alpha=R_{c}/R_{f}$ and the fact that $(R_{f}+1/\alpha)/(R_{f}+1)\approx 1$ for the parameters considered in the numerical examples of Section 6 – i.e., $R_{f}\geq 256$ and $\alpha\geq 1/5$. We obtain from an analogous derivation that $C^{fi}_{c}\approx\alpha^{2}C^{fi}_{f}$ and $C^{fe}_{c}\approx\alpha^{2}C^{fe}_{f}$. For simplicity, we further assume that the cost of evaluating each right-hand side is equal to that of solving the implicit system – even though, in practice, evaluating the nonlinear term requires applying the SH transform and is therefore more expensive than the other operations. This assumption reads

[TABLE]

Then, we use (37) and (38) to eliminate $C^{s}_{f}$ and $C^{s}_{c}$ from $C^{\textit{PFASST}(n_{\textit{ts}},M_{f}+1,M_{c}+1,N_{\textit{PF}},\alpha)}$ and $C^{\textit{SDC}(M_{f}+1,N_{\textit{S}})}$. After dividing both $C^{\textit{PFASST}(n_{\textit{ts}},M_{f}+1,M_{c}+1,N_{\textit{PF}},\alpha)}$ and $C^{\textit{SDC}(M_{f}+1,N_{\textit{S}})}$ by $3n_{\textit{ts}}M_{f}$, we obtain a compact form of the theoretical speedup of PFASST-SH with respect to SDC

[TABLE]

where the coefficients are defined as $a=2/3$, $b=(3M_{f}+2)/(3n_{\textit{ts}}M_{f})$, $c=(5M_{c}+2)/(3n_{\textit{ts}}M_{f})$, and d=\big{(}2M_{c}+(3M_{c}+2)n_{\textit{ts}}+2\big{)}/(3n_{\textit{ts}}M_{f}).

The above derivation assumes that PFASST-SH and SDC can achieve the same accuracy for the same time step size, which is only true if and when the residuals in the PFASST-SH iterations converge to the same level as as the SDC scheme. Whenever this is not the case, we rescale the expression of (39) by the ratio $\Delta t_{\textit{PF}}/\Delta t_{\textit{SDC}}$. The theoretical speedup of PFASST-SH with respect to MLSDC-SH, whose derivation is omitted, reads

[TABLE]

In (40), we have introduced $e=(3M_{f}+2)/(5n_{\textit{ts}}M_{c})$, $f=(5M_{c}+2)/(5n_{\textit{ts}}M_{c})$, and g=\big{(}2M_{c}+(3M_{c}+2)n_{\textit{ts}}+2\big{)}/(5n_{\textit{ts}}M_{c}). We will compare the theoretical speedups obtained in this section with the observed speedups based on wall-clock time in the next section.

6 Numerical examples

We are ready to analyze the accuracy, stability, and computational cost of the proposed parallel-in-time integration method to determine the conditions in which PFASST-SH reduces the time-to-solution. We use numerical examples of increasing complexity to successively evaluate the sensitivity of these PFASST-SH properties to the time step size, the number of iterations, and the number of processors. The examples are performed using the finest spatial resolution given by $R_{f}=256$ in (7), which is sufficient to generate complex flow dynamics over the time integration window. The error is computed with the following norm on the spectral coefficients

[TABLE]

where $\phi$ is the approximate solution and $\phi_{\textit{ref}}$ is the reference solution. We note that in (41), the norm of the difference between the computed solution and the reference solution is normalized by the norm of the reference solution. If we set $R_{\textit{norm}}=R_{f}$ in (41), $||\cdot||_{\infty,\,R_{f}}$ is the standard $L_{\infty}$-norm in spectral space that measures the ability of PFASST-SH to capture the fully detailed solution. If $R_{\textit{norm}}<R_{f}$, $||\cdot||_{\infty,\,R_{\textit{norm}}}$ is a semi-norm that assesses the accuracy of PFASST-SH for the large-scale features represented by the slow modes of the spectral basis. As mentioned in Section 1, the use of damping methods to resolve the spectral blocking problem – here, second-order artificial diffusion – implies that the accuracy on the modes corresponding to high-frequency wavenumbers is not of primary interest. Hence, in the numerical examples, we consider $R_{\textit{norm}}=32$, $R_{\textit{norm}}=64$, and $R_{\textit{norm}}=R_{f}=256$.

The studies presented next have been conducted on the Cori supercomputer at Lawrence Berkeley National Laboratory. Even though the following test cases are run without any parallelization in space for simplicity, our parallel-in-time method is meant to allow access to an additional axis of parallelism whenever spatial parallelization is saturated.

6.1 Propagation of a Gaussian dome

We first consider the propagation of a Gaussian bump on the sphere. This example is adapted from Swarztrauber, (2004) and was used in Hamon et al., (2019) to assess the performance of MLSDC-SH. The initial velocities are set to zero ($u=v=0$), and a Gaussian dome centered at $\lambda_{c}=\pi$ and $\phi_{c}=\pi/4$ is introduced in the geopotential field, such that the height is given by

[TABLE]

where $a$ denotes the earth radius. The distance $d$ is defined as

[TABLE]

with

[TABLE]

We use the same values as in Swarztrauber, (2004) for the Earth radius, the gravitational acceleration, and the angular rate of rotation $\Omega$ involved in the Coriolis force. In (42), we set $\bar{h}=29,400\,\text{m}$ and $A=6,000\,\text{m}$. We simulate $102,400\,\text{s}$ ($28\,\text{hours}$, $26\,\text{min}$, $40\,\text{s}$) of propagation with a spatial resolution $R_{f}=256$ and a diffusion coefficient $\nu=10^{5}\,\text{m}^{2}.\text{s}^{-1}$. PFASST-SH is run with 16 processors, using one time step per compute node. The geopotential field obtained using PFASST($n_{\textit{ts}}=16$, $M_{f}+1=5$, $M_{c}+1=3$, $N_{\textit{PF}}=8$, $\alpha=4/5$) – i.e., with five SDC nodes on the fine level, three SDC nodes on the coarse level, eight iterations, and a spatial coarsening ratio $\alpha=4/5$ – with a time step size $\Delta t=100\,\text{s}$ is shown in Fig. 2. This figure also shows the difference between this PFASST-SH solution and the reference solution obtained using SDC($M_{f}+1=5$, $N_{S}=8$) based on the same spatial resolution and $\Delta t_{\textit{ref}}=60\,\text{s}$, illustrating the fact that PFASST-SH can accurately reproduce the reference geopotential map without numerical artifacts. The spectrum of the reference solution is in Fig. 3.

To gain an understanding of the accuracy and stability of PFASST-SH, we perform a study of the error norm as a function of time step size for a fixed number of iterations. The geopotential error with respect to the reference solution is shown as a function of time step size in Fig. 4. In this figure and in the numerical results shown in the remainder of this paper, we consider a wide range of normalized errors – from $10^{-4}$ to $10^{-11}$ – to compare PFASST-SH and SDC in a regime where the SDC schemes achieve their formal order of accuracy – up to eighth order for SDC(5,8). We focus on PFASST($n_{\textit{ts}}=16$, $M_{f}+1=3$, $M_{c}+1=2$, $N_{\textit{PF}}=4$, $\alpha$) and PFASST($n_{\textit{ts}}=16$, $M_{f}+1=5$, $M_{c}+1=3$, $N_{\textit{PF}}=8$, $\alpha$) with four and eight iterations, respectively. We explain at the end of this section how this fixed number of iterations, $N_{\textit{PF}}$, is determined. We consider the spatial coarsening ratios $\alpha=1/5,\,1/2,\,\text{and}\,4/5$. The range of error magnitudes observed here is consistent with the results in the numerical examples of Jia et al., (2013).

Figure 4 illustrates two key properties of PFASST-SH. First, considering stability, we observe that the PFASST-SH schemes can take relatively large time steps, on the order of $10^{2}\,\text{s}$, but we also note that their largest time step size decreases when the spatial coarsening ratio, $\alpha$, is increased. The largest time steps taken by PFASST-SH are about ten times smaller than with the corresponding serial SDC schemes. Second, we see that the PFASST-SH accuracy heavily depends on the choice of a spatial coarsening ratio, as this parameter determines the portion of the spectrum of the solution that can be resolved on the coarse level (see Fig. 3).

This is because for an aggressive spatial coarsening – i.e., $\alpha$ close to zero –, the main features of the solution cannot be captured on the coarse level and as a result, the PFASST-SH error stagnates at a relatively large magnitude. Conversely, for a milder spatial coarsening – with $\alpha$ close to one –, higher-frequency modes are included on the coarse level. This significantly improves the accuracy of PFASST-SH, but also imposes a more severe stability limit on the time step size.

Finally, for $\alpha=1/2$ and $\alpha=4/5$, we observe that the PFASST-SH accuracy better matches that of the SDC schemes when we only consider the low-frequency modes in the computation of the error ($R_{\textit{norm}}=32$). This means that PFASST-SH can efficiently resolve the spectral coefficients corresponding to the large-scale features of the solution, at the expense of the accuracy for the small-scale features.

We now evaluate the computational cost of the parallel-in-time scheme. In Fig. 5, the geopotential error norm is shown as a function of wall-clock time for the same PFASST-SH schemes as in Fig. 4. We observe in Table 3 that for their largest time step size, these PFASST-SH schemes are consistently faster than the corresponding serial SDC schemes for the range of spatial coarsening ratios considered here.

To achieve a normalized error norm of about $10^{-6}$ in the geopotential variable, PFASST($n_{\textit{ts}}=16$, $M_{f}+1=3$, $M_{c}+1=2$, $N_{\textit{PF}}=4$, $\alpha=1/5$) has to take a time step size smaller or equal to $\Delta t=128\,\text{s}$, whereas SDC($n_{\textit{ts}}=3$, $N_{S}=4$) and MLSDC($M_{f}+1=3$, $M_{c}+1=2$, $N_{\textit{ML}}=2$, $\alpha=1/2$) can achieve the same accuracy with a time step size $\Delta t=400\,\text{s}$. Still, for this relatively low accuracy, PFASST-SH yields a speedup of $S^{\textit{obs}}_{S}=4.1$ (respectively, $S^{\textit{obs}}_{\textit{ML}}=3.0$) compared to SDC($M_{f}+1=3$, $N_{S}=4$) (respectively, MLSDC($M_{f}+1=3$, $M_{c}+1=2$, $N_{\textit{ML}}=2$, $\alpha=1/2$)). These observed speedups are consistent with the values derived in Section 5.5. With a milder coarsening in space, the accuracy of the parallel-in-time algorithm improves, but the speedup is slightly reduced. Specifically, for a normalized error norm of about $10^{-9}$, PFASST($n_{\textit{ts}}=16$, $M_{f}+1=5$, $M_{c}+1=3$, $N_{\textit{PF}}=8$, $\alpha=4/5$), achieves a speedup of $S^{\textit{obs}}_{S}=3.7$ (respectively, $S^{\textit{obs}}_{\textit{ML}}=2.5$) compared to SDC($M_{f}+1=5$, $N_{S}=8$) (respectively, MLSDC($M_{f}+1=5$, $M_{c}+1=3$, $N_{\textit{ML}}=4$, $\alpha=1/2$)). We note that the MLSDC(5,3,4,1/2) error is about 20 times larger than the SDC(5,8) error for the same time step size since the four coarse sweeps do not resolve well the high-frequency modes of the solution. This is documented in the discussion of MLSDC-SH found in Hamon et al., (2019).

To conclude this section, we motivate the choice of the number of PFASST-SH iterations by studying the sensitivity of the algorithm to $N_{\textit{PF}}$ for a fixed time step size ($\Delta t=100\,s$). These results are shown in Fig. 5. We see that for a small $N_{\textit{PF}}$ (for instance, $N_{\textit{PF}}<6$ for PFASST($n_{\textit{ts}}=16$, $M_{f}+1=5$, $M_{c}+1=3$, $\cdot$, $\alpha=4/5$)), performing an additional iteration slightly increases the computational cost, but significantly reduces the error. Instead, for a larger $N_{\textit{PF}}$ (larger than the threshold $N_{\textit{PF}}=6$ for PFASST($n_{\textit{ts}}=16$, $M_{f}+1=5$, $M_{c}+1=3$, $\cdot$, $\alpha=4/5$)), increasing the number of iterations still increases the computational cost, but does not reduce the error anymore. This remaining error persists because the high-frequency modes cannot be captured on the coarse level.

We used this sensitivity analysis to choose the fixed $N_{\textit{PF}}$ in the refinement studies on the time step size (Figs. 4 and 5). Specifically, we set $N_{PF}$ slightly above the threshold mentioned in the previous paragraph to minimize the number of iterations that increase the computational cost without reducing the error. This methodology yielded $N_{\textit{PF}}=4$ for PFASST($n_{\textit{ts}}=16$, $M_{f}+1=3$, $M_{c}+1=2$, $\cdot$, $\alpha=4/5$) and $N_{\textit{PF}}=8$ for PFASST($n_{\textit{ts}}=16$, $M_{f}+1=5$, $M_{c}+1=3$, $\cdot$, $\alpha=4/5$). A more flexible way to choose $N_{\textit{PF}}$ would rely on a convergence check to stop the iterations, which we will implement in future work.

6.2 Rossby-Haurwitz wave

We continue the analysis of the PFASST-SH algorithm with the Rossby-Haurwitz wave test case. The initial velocity field is non-divergent, and the initial geopotential distribution, given in Williamson et al., (1992), is such that the temporal derivative of the divergence is zero. The diffusion coefficient is set to $10^{5}\,\text{m}^{2}.\text{s}^{-1}$. We simulate $102,400\,\text{s}$ ($28\,\text{hours}$, $26\,\text{min}$, $40\,\text{s}$) of propagation with a spatial resolution given by $R_{f}=256$. We note that when PFASST-SH is used with an aggressive coarsening ratio and/or large time steps to simulate the problem over a longer period of time, the wave becomes unstable, as documented for standard time integration schemes in Thuburn and Li, (2000); Bénard, (2019). The spectrum of the reference solution obtained using SDC($M_{f}+1=5$, $N_{S}=8$) with the same spatial resolution and a time step size $\Delta t_{\textit{ref}}=100\,\text{s}$ is shown in Fig. 7. As in Section 6.1, we see that even though relatively large spectral coefficients are truncated on the coarse level for $\alpha=1/5$, they are much smaller for $\alpha=4/5$.

We first study the error norm as a function of time step size for a fixed number of iterations. We consider the same PFASST-SH schemes with the same number of iterations as in the previous test case, namely, PFASST($n_{\textit{ts}}=16$, $M_{f}+1=3$, $M_{c}+1=2$, $N_{\textit{PF}}=4$, $\alpha$) and PFASST($n_{\textit{ts}}=16$, $M_{f}+1=5$, $M_{c}+1=3$, $N_{\textit{PF}}=8$, $\alpha$). Here, the fixed number of iterations, $N_{\textit{PF}}$, has been chosen using the approach described at the end of the previous section.

The results shown in Fig. 8 are consistent with those of the previous test case. Specifically, the largest stable time steps achieved with PFASST-SH decrease when the spatial coarsening ratio is close to one, and are still smaller than with the serial SDC schemes. In addition, the error obtained with the PFASST-SH schemes first decreases as the time step size is reduced, and then stagnates for small time steps. As discussed earlier, the error norm achieved for small time steps depends on the spatial coarsening ratio as this parameter determines the fraction of the spectrum that can be represented on the coarse level. Figure 8 shows that the error norm remains large with aggressive spatial coarsening, but essentially vanishes for mild coarsening in space. We observe that PFASST($n_{\textit{ts}}=16$, $M_{f}+1=5$, $M_{c}+1=3$, $N_{\textit{PF}}=8$, $\alpha=4/5$) is as accurate as the serial SDC($M_{f}+1=5$, $N_{S}=8$) for all the time step sizes considered here. PFASST($n_{\textit{ts}}=16$, $M_{f}+1=5$, $M_{c}+1=3$, $N_{\textit{PF}}=8$, $\alpha=1/2$) is less accurate than PFASST(16,5,3,8,4/5) for $R_{\textit{norm}}=256$, but the two schemes have the same accuracy for $R_{\textit{norm}}=32$.

Next, we assess the computational cost of these PFASST-SH schemes by considering the error norm as a function of wall-clock time. As in the previous test case, we note in Fig. 9 that for their largest stable time steps, the PFASST-SH schemes are more efficient than the corresponding serial SDC schemes. We also see that for a given accuracy, the difference in wall-clock time is larger when the error is computed with $R_{\textit{norm}}=32$, which suggests that the PFASST-SH algorithm resolves the main features of the solution significantly faster than the serial schemes.

These results are summarized in Table 4. We observe that PFASST($n_{\textit{ts}}=16$, $M_{f}+1=3$, $M_{c}+1=2$, $N_{\textit{ML}}=4$, $\alpha=1/5$) can achieve a geopotential error norm of $2\times 10^{-8}$ for a time step size of $\Delta t=100\,\text{s}$, which is a quarter of the time step size required by SDC($M_{f}+1=3$, $N_{S}=4$) and MLSDC($M_{f}+1=3$, $M_{c}+1=2$, $N_{\textit{ML}}=2$, $\alpha=1/2$) to reach this tolerance. The reduced stability of the parallel-in-time scheme explains why the speedups obtained with PFASST($n_{\textit{ts}}=16$, $M_{f}+1=3$, $M_{c}+1=2$, $N_{\textit{PF}}=4$, $\alpha=1/5$) are slightly smaller than for the previous test case. We observe a speedup of $S^{\textit{obs}}_{S}=3.2$ compared to SDC($M_{f}+1=3$, $N_{S}=4$) and of $S^{\textit{obs}}_{ML}=2.2$ compared to MLSDC($M_{f}+1=3$, $M_{c}+1=2$, $N_{\textit{ML}}=2$, $\alpha=1/2$), which is close to the theoretical speedup of Section 5.5. This remark also applies to the speedup obtained with PFASST($16$, $M_{f}+1=5$, $M_{c}+1=3$, $N_{\textit{PF}}=8$, $\alpha=1/2$).

6.3 Nonlinear evolution of an unstable barotropic wave

We conclude the analysis of PFASST-SH by considering the initial condition proposed by Galewsky et al., (2004). This test case consists of a stationary zonal jet perturbed by the introduction of a Gaussian bump in the geopotential field. This leads to the propagation of gravity waves followed by the development of complex vortical dynamics. These processes operate on multiple time scales and are representative of the horizontal aspects of atmospheric flows. We run the simulations for 144 hours using a diffusion coefficient of $\nu=10^{5}\,\text{m}^{2}.\text{s}^{-1}$ as in Galewsky et al., (2004).

We have shown in Hamon et al., (2019) that this is a particularly challenging test for the MLSDC-SH scheme. This is due to the fast amplification of small-scale features that significantly undermines the accuracy of the multi-level integration scheme when the spatial coarsening strategy is too aggressive. Based on a result from our previous work on MLSDC-SH, we only consider coarsening ratios closer to one – i.e., $\alpha=R_{c}/R_{f}\geq 3/5$ – in this example. This test case is therefore well suited to assess the robustness and efficiency of PFASST-SH when the coarse corrections are relatively expensive due to the limited spatial coarsening.

The vortices generated with PFASST-SH are in Fig. 10 for a resolution given by $R_{f}=256$, along with the error in the vorticity field with respect to a reference solution computed with SDC($M_{f}+1=5$, $N_{S}=8$). This reference solution relies on a time step size $\Delta t_{\textit{ref}}=60\,\text{s}$ and the same spatial resolution. Figure 10 shows that PFASST-SH can accurately match the vorticity field generated with SDC($M_{f}+1=5$, $N_{S}=8$), with small errors concentrated in the vortical structures of the flow. The spectrum of the reference is shown in Fig. 11. Comparing the vorticity spectrum with those of Figs. 3 and 7 highlights the amplification of large spectral coefficients corresponding to high-frequency features in this test case.

We now proceed to the analysis of the evolution of the error norm as a function of time step size in Fig. 12. We focus on PFASST($n_{\textit{ts}}=16$, $M_{f}+1=3$, $M_{c}+1=2$, $N_{\textit{PF}}=4$, $\alpha$) and PFASST($n_{\textit{ts}}=16$, $M_{F}+1=5$, $M_{c}+1=3$, $N_{\textit{PF}}=8$, $\alpha$) with four and eight iterations, respectively. We have chosen these values of $N_{\textit{NF}}$ to be consistent with the parameters used in the other test cases. For this example characterized by the presence of large spectral coefficients associated with high-frequency modes, PFASST-SH is accurate only if the spatial coarsening ratio is close to one. But, for a small value of $\alpha$, most of the spectrum is represented on the coarse level, which imposes a more severe stability restriction on the time step size. This is particularly the case for PFASST($n_{\textit{ts}}=16$, $M_{f}+1=3$, $M_{c}+1=2$, $N_{\textit{PF}}=4$, $\alpha$), as this scheme is only stable for $\Delta t\leq 15\,\text{s}$. This is almost two orders of magnitude smaller than the largest time step sizes achieved by the serial schemes.

With more SDC nodes on the fine and coarse levels, PFASST($n_{\textit{ts}}=16$, $M_{F}+1=5$, $M_{c}+1=3$, $N_{\textit{PF}}=8$, $\alpha$) can still achieve stable time step sizes as large as $240\,\text{s}$ for $\alpha=3/5$, which is only one order of magnitude smaller than the largest time step size taken by SDC($M_{f}+1=5$, $N_{S}=8$). However, we see that reducing the time step size for a fixed number of eight iterations in PFASST($n_{\textit{ts}}=16$, $M_{F}+1=5$, $M_{c}+1=3$, $N_{\textit{PF}}=8$, $\alpha$) does not reduce the error. We show later in this section that for the Galewsky benchmark test, PFASST-SH performs better with a larger number of iterations.

Figure 13 illustrates the computational cost of PFASST-SH by showing the error norm as a function of wall-clock time. Due to its limited stability, PFASST($n_{\textit{ts}}=16$, $M_{f}+1=3$, $M_{c}+1=2$, $N_{\textit{PF}}=4$, $\alpha$) is not competitive for this test case, as it is consistently slower than the serial SDC schemes for all time step sizes. Conversely, for their largest stable time step, the PFASST($n_{\textit{ts}}=16$, $M_{f}+1=5$, $M_{c}+1=3$, $N_{\textit{PF}}=8$, $\alpha$) schemes are more efficient than both SDC($M_{f}+1=5$, $N_{S}=8$) and SDC($M_{f}+1=3$, $N_{S}=4$). The speedup is slightly larger when the error norm only accounts for the low-frequency spectral coefficients ($R_{\textit{norm}}=64$). This is because for $\alpha<1$, the high-frequency features cannot be resolved by PFASST-SH when the number of iterations is small compared to the number of time steps.

In Figs. 12 and 13, the error norm obtained with the PFASST-SH schemes stagnates at a relatively large magnitude. Here, we study the sensitivity of PFASST-SH to the number of iterations, $N_{\textit{PF}}$ using the same methodology as in Fig. 5 to show that this limitation can be overcome by increasing $N_{\textit{PF}}$. In Fig. 14, the PFASST-SH results are obtained with a fixed time step size $\Delta t=240\,\text{s}$, a spatial coarsening ratio $\alpha\in\{3/5,4/5,19/20,1\}$ and different values of the number of iterations, $N_{\textit{PF}}$. We now focus on PFASST($n_{\textit{ts}}=16$, $M_{f}+1=5$, $M_{c}+1=3$, $N_{\textit{PF}}$, $\alpha$) with $N_{PF}\geq 8$ as this family of schemes yields consistent speedups compared to the serial SDC schemes.

Figure 14 shows that increasing the number of PFASST-SH iterations reduces the error in the spectral coefficients for both $R_{\textit{norm}}=256$ and $R_{\textit{norm}}=64$ until a threshold is reached at approximately $N_{\textit{PF}}=16$. We note that for this challenging test case, the threshold is much higher than for the Gaussian test case of Section 6.1. This error reduction for $N_{\textit{PF}}<16$ takes place because increasing $N_{\textit{PF}}$ allows to progressively propagate more information on high-frequency modes to the last processor.

In terms of computational cost, performing more PFASST-SH iterations only slightly increases the wall-clock time. This is because in PFASST-SH, the cost of the parallel fine sweeps is amortized over multiple iterations. Therefore, unlike Parareal, PFASST-SH can still achieve a speedup when the number of iterations, $N_{\textit{PF}}$, is close or equal to the number of processors, $n_{\textit{ts}}$. In fact, we see in Fig. 14 that increasing $N_{\textit{PF}}$ leads to an increase in the speedup compared to the serial schemes.

The speedups observed for this experiment are in Table 5. We highlight that by running the PFASST-SH simulations for different numbers of iterations and spatial coarsening ratios, we found a near-optimal set of parameters for this challenging test case. In the best configuration, PFASST-SH uses no coarsening in space ($\alpha=1$) and performs 16 iterations ($N_{\textit{PF}}=n_{\textit{ts}}$). This yields an accurate parallel-in-time scheme with a good stability. Specifically, for an error norm of $3\times 10^{-12}$, PFASST($n_{\textit{ts}}=16$, $M_{f}+1=5$, $M_{c}+1=3$, $N_{\textit{PF}}=16$, $\alpha=1$) can take the same time step size of $240\,\text{s}$ as SDC($M_{f}+1=5$, $N_{S}=8$), leading to a speedup of 3.7. For this error norm, MLSDC($M_{f}+1=5$, $M_{c}+1=3$, $N_{\textit{ML}}=4$, $\alpha=4/5$) is more expensive than the single-level SDC($M_{f}+1=5$, $N_{S}=8$), as shown in Hamon et al., (2019).

Finally, we run a strong scalability study to assess the evolution of the accuracy and efficiency of PFASST-SH when the number of processors is increased from four to 64. The results of this study run with the same resolution ($R_{f}=256$) and physical parameters as in the previous paragraphs are in Fig. 15. We simulate 144 hours with a time step of $\Delta t=225\,s$. Figure 15 shows that we can significantly reduce the wall-clock time by increasing the number of processors, and therefore, the number of time steps solved in parallel. Specifically, the wall-clock time of PFASST($n_{\textit{ts}}=16$, $M_{f}+1=5$, $M_{c}+1=3$, $N_{\textit{PF}}=30$, $\alpha=3/5$) is reduced by a factor of 7.8 (respectively, 13.7) when the number of processors is increased from four to 32 (respectively, 64). This is an encouraging result although it does not prove that PFASST-SH can improve the scalability of global spectral models – this would require performing the same study on a large-scale problem on which parallelism in space is saturated. But, Fig. 15 shows that increasing the number of processors causes a deterioration of the error as PFASST-SH becomes unstable. This is a result that will be improved in future work with the research directions discussed in the conclusion.

7 Conclusion and future work

We have proposed an iterative, multilevel, parallel-in-time integration method for the Shallow-Water Equations (SWE) on the rotating sphere. Our approach combines the global Spherical Harmonics (SH) transform with the Parallel Full Approximation Scheme in Space and Time (PFASST) to construct a robust and efficient numerical scheme for the nonlinear wave-propagation problems arising from the SWE. The method computes multiple time steps in parallel by applying a sequence of concurrent fine corrections to iteratively improve an initial condition propagated on the coarse level. A key feature of this algorithm is the coarsening and interpolation procedure designed in Hamon et al., (2019) to accurately transfer the approximate solution between levels without introducing spurious modes.

We have studied the properties of PFASST-SH using a suite of standard nonlinear test cases designed for the development of dynamical cores for weather and climate simulations. These numerical tests illustrate the importance of the spatial coarsening ratio, $\alpha$, for the accuracy and computational cost of the scheme. Aggressive spatial coarsening yields a PFASST-SH scheme that only partially resolves the small-scale features of the solution, and that is therefore inexpensive but also relatively inaccurate compared to the serial single-level and multi-level SDC schemes. Conversely, milder spatial coarsening leads to a more costly, but significantly more accurate scheme that captures the full spectrum of the solution. Our results show that PFASST-SH remains robust and efficient for the challenging test case proposed in Galewsky et al., (2004) and characterized by the progressive amplification of large high-frequency modes in the spectrum of the solution. For this numerical experiment, we have demonstrated that our parallel-in-time algorithm can both fully represent the fine solution of the problem and achieve a speedup compared to the serial single-level and multi-level schemes.

Future work includes combining PFASST with Semi-Lagrangian (SL) methods, which have exhibited excellent stability properties in the context of a Parareal-based predictor-corrector scheme (Schmitt et al.,, 2018). This would aim at increasing the largest stable time step size taken by PFASST to further reduce the cost of the parallel-in-time integration method. Another research direction consists in exploiting the potential of exponential integrators in PFASST, as done before using a Parareal approach (Gander and Güttel,, 2013) and a serial SDC approach (Buvoli,, 2015). Exponential integration could advantageously replace the backward-Euler step in the PFASST iteration to reduce the dispersion errors generated by the current scheme. Finally, future work also includes extending these results to the full atmospheric dynamics – i.e., accounting for the horizontal and vertical components – to demonstrate the applicability of PFASST in dynamical cores and target realistic, larger-scale weather and climate simulations.

8 Acknowledgements

The work of François Hamon and Michael Minion was supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Applied Mathematics program under contract number DE-AC02005CH11231. The simulations were performed at the National Energy Research Scientific Computing Center (NERSC), a DOE Office of Science User Facility supported by the U.S. DOE under Contract No. DE-AC02-05CH11231. Martin Schreiber received funding from NCAR for a research stay in summer 2017 at the Mesa Labs which also led to contributions to the present paper.