# A non-spectral Helmholtz solver for numerical weather prediction models with a mass-based vertical coordinate

Daan Degrauwe, Fabrice Voitus, Piet Termonia

► To cite this version:

Daan Degrauwe, Fabrice Voitus, Piet Termonia. A non-spectral Helmholtz solver for numerical weather prediction models with a mass-based vertical coordinate. Quarterly Journal of the Royal Meteorological Society, 2020, 147 (734), pp.30-44. 10.1002/qj.3902. hal-04950244

HAL Id: hal-04950244

https://cnrs.hal.science/hal-04950244v1

Submitted on 23 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A non-spectral Helmholtz solver for numerical weather prediction models with a mass-based vertical coordinate

Daan Degrauwe<fnref n="1" /> | Fabrice Voitus<fnref n="2" /> | Piet Termonia<fnref n="1" />,<fnref n="3" />

## Correspondence

D. Degrauwe, Royal Meteorological Institute, Ringlaan 3, 1180 Brussels, Belgium.

Email: daan.degrauwe@meteo.be

## Funding information

European Union's Horizon 2020 research and innovation programme, Grant/Award

Number: GA 67162

## Abstract

Although semi-implicit semi-Lagrangian spectral atmospheric models have been very successful for decades, they are believed to face big challenges in the longer term. Foremost, the spectral method relies heavily on data-rich global communications, which may become problematic on future massively parallel machines. This paper investigates how the Helmholtz problem, as it arises in the dynamical core of a semi-implicit non-hydrostatic numerical weather prediction model with a mass-based vertical coordinate and a constant-coefficient reference state, can be solved efficiently without relying on spectral transforms, by using a multigrid-preconditioned iterative solver instead. In the particular case of a limited-area geometry, the convergence rate of this iterative solver can be determined a priori, which allows us to predict the required number of iterations. This knowledge is especially valuable for an atmospheric model that is used for operational weather forecasting, because it guarantees that the model runtime stays constant from one forecast to another. The a priori knowledge of the convergence rate also allows us to choose the parameters of the multi-grid preconditioner optimally. Weak scalability experiments show the superior scalability of this solver with respect to a spectral solver.

## KEYWORDS

Helmholtz problem, Krylov solver, scalability, semi-implicit timestepping, spectral methods

# 1 | INTRODUCTION

The development of numerical weather prediction (NWP) models is constantly being driven by a quest for higher resolutions at a reasonable computing cost. A decisive criterion in the choice of a time integration scheme and the model's timestep is the stability of the model. In this regard, semi-implicit methods (Robert, 1981) have contributed significantly, since such schemes allow for large timesteps without sacrificing stability. Especially the combination of a semi-implicit time integration method with a semi-Lagrangian advection scheme to treat the

nonlinear advective terms, as presented by Robert (1982), has shown to be very successful. The price to pay for this efficiency is that an implicit linear system needs to be solved at every timestep. This linear system takes the shape of a three-dimensional elliptic sign-positive Helmholtz equation. When formulated in spectral space, the solution of the implicit system becomes trivial. Moreover, spectral methods offer an unchallenged accuracy. This has led to the development of semi-implicit semi-Lagrangian spectral models such as the global Integrated Forecast System (IFS; Barros et al., 1995; Wedi et al., 2015) and ARPEGE (Courtier and Geleyn, 1988; Courtier et al., 1991) models,

<footnote n="1">¹Royal Meteorological Institute, Brussels, Belgium</footnote>

<footnote n="2">²CNRM, Météo-France, Toulouse, France</footnote>

<footnote n="3">³Ghent University, Ghent, Belgium</footnote>
and the associated limited-area ALADIN model (ALADIN International Team, 1997; Termonia *et al.*, 2018).<fnref n="1" />

However, some concerns exist on whether the use of spectral methods in atmospheric models can be future-proofed. With the anticipated arrival of exascale computing, the scalability of atmospheric models should be a top research priority (Bauer *et al.*, 2015; Müller *et al.*, 2019). Because spectral transforms rely on the exchange of data throughout the model domain, the scalability of spectral transforms on massively parallel machines is somewhat questionable (Steppeler *et al.*, 2003; Wedi *et al.*, 2013; Kühnlein *et al.*, 2019). It should be mentioned here that, from a theoretical point of view, on a parallel computer system where the bandwidth between different compute nodes can be maintained with increasing system size, the all-to-all communication patterns of spectral methods should scale very well (Zheng and Marguinaud, 2018). Nevertheless, actual tests show that the scalability of spectral transforms is far from perfect in practice (Michalakes *et al.*, 2015).

A second concern about the use of spectral methods in semi-implicit timestepping is that the atmospheric reference state, around which the equations are linearized, needs to be horizontally homogeneous (Thomas *et al.*, 2003; Bénard *et al.*, 2010). While this has not been much of an issue so far because at current resolutions the terms responsible for the fastest waves can still be treated adequately, it is unsure whether this will remain the case in the future. Indeed, going to higher resolutions means that resolved orographic slopes become steeper, and this effect may generate instabilities if not treated implicitly in the model's timestepping (Ikawa, 1988; Skamarock *et al.*, 1997; Kühnlein *et al.*, 2019). The focus of the present paper lies on the scalability aspect of semi-implicit NWP models, in particular those that currently use a spectral solver. The topic of instabilities generated by an explicit treatment of the orography at high resolutions is not dealt with here. However, as will be indicated later, the work presented in this paper can serve as a preparatory step to treat orographic effects in a more stable way.

Different groups of researchers are considering different solutions to these challenges for spectral methods. At the Japan Meteorological Agency (JMA), a limited-area model (Saito *et al.*, 2006) has been developed that is implicit only in the vertical, while horizontal coupling is dealt with in a time-explicit way. At Environment and Climate Change Canada (ECCC), a new dynamical core for the Global Environmental Multiscale model is being developed (Husain *et al.*, 2019) which is characterized by

a novel height-based terrain-following vertical coordinate, and where a Krylov solver is implemented as an alternative to a direct Fourier-based solver for long-term scalability reasons. At the European Centre for Medium-Range Weather Forecasts (ECMWF), alternatives to the current semi-implicit semi-Lagrangian spectral-transform dynamical core of the IFS model are considered. Kühnlein *et al.* (2019) developed a finite-volume-based version of the IFS model with a non-oscillatory forward-in-time advection scheme, and a preconditioned Krylov solver for the Helmholtz equation. Among others, this latter work involves a change of the vertical coordinate from mass-based to height-based in the overall NWP model formulation.

In this paper, a different, less intrusive, approach is followed for the dynamical core of the non-hydrostatic spectral ALADIN model. Instead of redeveloping the dynamical core entirely, it is investigated whether the Helmholtz equation of this model can be solved efficiently without relying on spectral transforms. The advantage of this approach is that the lion's share of the model does not need to be adapted. For instance, the physics parametrizations (Termonia *et al.*, 2018, and references therein) are coupled to the dynamical core with a flux-conservative interface (Catry *et al.*, 2007; Degrauwe *et al.*, 2016) which relies on a mass-based coordinate. Changing the vertical coordinate system would have severe consequences, not only for the physics parametrizations themselves, but also for the way they are coupled to the dynamical core. The idea of changing only the Helmholtz solver of the ALADIN model while keeping the rest of the model intact, builds further on the work by Caluwaerts *et al.* (2015), who have investigated the scientific impact of replacing the spectral solver by a (less accurate) non-spectral solver.

Besides this practical motivation for staying close to the existing dynamics of the ALADIN model, this paper shows that several aspects of this dynamical core turn out to be very favourable for an iterative non-spectral Helmholtz solver as well. Most notably, in this dynamical core, the reference state for the semi-implicit scheme is taken horizontally homogeneous (Bénard *et al.*, 2004), which means that all operators and coefficients are constant in space and time. It will be demonstrated in this paper that, for the specific case of such a constant-coefficient Helmholtz problem, it is possible to determine the convergence rate of an iterative solver *a priori*. The importance of this knowledge can hardly be overestimated. Indeed, for an atmospheric model that is used for operational weather forecasting in 26 European and North African countries (Termonia *et al.*, 2018), it is essential that the forecast finishes within a predefined time window. One simply cannot afford to occasionally have delayed forecasts due to an unpredictable convergence of an iterative Helmholtz

<footnote n="1">¹ARPEGE=Action de Recherche Petite Échelle Grande Échelle;
ALADIN=Aire Limitée Adaptation dynamique Développement
InterNational.</footnote>
solver. It is fair to say that, when the convergence rate is not guaranteed theoretically, one can (and does, as illustrated by Wood et al., 2014) limit the maximum number of iterations, and accept the (small) risk of local loss of accuracy in the case of slow convergence. The main outcome of the present paper is that, thanks to the predictable convergence rate of an iterative solver for the Helmholtz problem as arising in a semi-implicit model with a mass-based vertical coordinate and a constant-coefficient reference state, the required number of iterations is known a priori, as is the runtime of the entire numerical forecast.

This paper is organized as follows. Section 2 formulates the problem under consideration, that is, the sign-positive Helmholtz equation as arising in a mass-based semi-implicit non-hydrostatic model. A multigrid-based preconditioner for an iterative solver for this Helmholtz equation is presented in section 3. A spectral analysis of the effect of this preconditioner predicts the convergence rate of the iterative Helmholtz solver. It will also be indicated how the a priori knowledge of the convergence rate can serve as a guideline to select the optimal preconditioner parameters. Section 4 shows the results of weak scalability experiments with a spectral solver and two non-spectral solvers. Finally, conclusions are formulated in section 5.

## 2 | THE HELMHOLTZ EQUATION IN A NON-HYDROSTATIC NWP MODEL WITH A MASS-BASED VERTICAL COORDINATE

### 2.1 | Formulation

The non-hydrostatic dynamical core of the spectral semi-implicit ALADIN model was originally presented in Bubnová et al. (1995), and later refined in Bénard et al. (2010). In the current paper, we repeat only the essentials, and we refer to those papers for further details. After a linearization around a constant-coefficient reference state, the time evolution of the prognostic variables in this system is written as

$$ \frac{\partial D}{\partial t} = -R_a G^* \nabla^2 T' + R_a T^* G^* \nabla^2 \hat{q} - R_a T^* \nabla^2 \hat{q} - \frac{R_a T^*}{\pi_s^*} \nabla^2 \pi_s', \quad (1) $$

$$ \frac{\partial d}{\partial t} = - \frac{g^2}{R_a T_e^*} L^* \hat{q}, \quad (2) $$

$$ \frac{\partial T'}{\partial t} = - \frac{R_a T^*}{C_{va}} (D + d), \quad (3) $$

$$ \frac{\partial \hat{q}}{\partial t} = S^* D - \frac{C_{pa}}{C_{va}} (D + d), \quad (4) $$

$$ \frac{\partial \pi_s'}{\partial t} = -\pi_s^* N^* D, \quad (5) $$

where D is the horizontal divergence, d a vertical divergence variable, $T' = T - T^*$ the temperature deviation w.r.t. a constant reference temperature $T^*$, $\hat{q} = \ln(p/\pi)$ a variable representing the deviation of the true pressure p w.r.t. the hydrostatic pressure $\pi$, and $\pi_s' = \pi_s - \pi_s^*$ the deviation of the surface pressure w.r.t. a constant reference surface pressure $\pi_s^*$. For the meaning of the other variables and the operators $G^*$, $S^*$, $N^*$ and $L^*$, we refer to Bénard et al. (2010). After time-discretization with a semi-implicit scheme with timestep $\Delta t$ and a vertical discretization as described in Bubnová et al. (1995), these equations become

$$ D - \delta t R_a T^* \nabla^2 [(G^* - I)\hat{q} - G^*(T'/T^*) - (\pi_s'/\pi_s^*)] = \tilde{D}, \quad (6) $$

$$ d - \delta t \left( - \frac{g^2}{R_a T_e^*} L^* \hat{q} \right) = \tilde{d}, \quad (7) $$

$$ T' + \delta t \frac{RT^*}{C_{va}} (D + d) = \tilde{T}, \quad (8) $$

$$ \hat{q} - \delta t \left[ S^* D - \frac{C_{pa}}{C_{va}} (D + d) \right] = \tilde{\hat{q}}, \quad (9) $$

$$ \pi_s' + \delta t \pi_s^* N^* D = \tilde{\pi}_s, \quad (10) $$

where $\delta t = \Delta t/2$. The prognostic variables D, d, T', $\hat{q}$ and $\pi_s'$ on the left-hand side are all considered at time $t + \Delta t$, and the right-hand sides are available explicitly at time t. It should be noted that the vertical operators $G^*$, $L^*$, $S^*$ and $N^*$ have been replaced by their respective vertically discrete analogue matrices $G^*$, $L^*$, $S^*$ and $N^*$.

A characteristic feature of the system (6)–(10) is that all coefficients are constant horizontally and in time. Although such a constant-coefficient system leads to residual nonlinear terms that are evaluated explicitly in time (Bénard et al., 2004), the use of a spectral Helmholtz solver necessitates this choice.

By further manipulation (detailed in the Appendix), the preceding system (6)–(10) can be reduced to a single Helmholtz equation in D:

$$ (I - \delta t^2 \nabla^2 \mathbf{B}_D^*) D = D^{\bullet\bullet} \quad (11) $$

where $\mathbf{B}_D^*$ is a square matrix with dimension the number of vertical levels, and the expression for the right-hand side $D^{\bullet\bullet}$ is given in the Appendix. By construction, $\mathbf{B}_D^*$ is positive definite.
<figure><img src="image_2.png" /></figure>

It is worth mentioning that this elimination procedure, requiring only the solution of a single Helmholtz equation, differs from the original formulation in Bubnová et al. (1995), which required the solution of two Helmholtz equations.

## 2.2 | Vertical decoupling

Since the coefficients in the 3D Helmholtz system (11) are horizontally constant, it is possible to decouple the system in the vertical dimension by considering the eigendecomposition of the matrix $B_D^*$:

$$B_D^* V_B = V_B \Lambda_B, \quad (12)$$

where $V_B$ is a matrix whose columns are the eigenvectors and $\Lambda_{bfB}$ is a diagonal matrix with the eigenvalues $\lambda_\ell$. Substitution in the 3D Helmholtz equation (11) then decouples the system into a set of N (the number of vertical levels) 2D Helmholtz equations:

$$(1 - c_{\ell}^2 \delta t^2 \nabla^2) \psi_{\ell} = \tilde{\psi}_{\ell} \quad \text{for } \ell = 1, \dots, N, \quad (13)$$

where $c_\ell = \sqrt{\lambda_\ell}$ and $\tilde{\psi}_\ell$ is the $\ell$th row of the matrix $V_B^{-1} D^{**}$.

## 2.3 | Spectral solver

A spectral solver looks for the solution of (13) in the space spanned by a set of orthogonal eigenfunctions of the Laplace operator. For a spherical geometry, these eigenfunctions are the spherical harmonics, while for a biperiodic rectangular domain, Fourier modes are used. In this paper, we will restrict the analysis and experiments to the limited-area case, that is, we will always consider a rectangular Cartesian grid. At the boundaries, biperiodic boundary conditions are applied. This corresponds to what is done in the ALADIN model: during the solution of the Helmholtz problem, which is done in spectral space, biperiodic boundary conditions are assumed. At a later moment in the model timestep, these (unphysical) boundary conditions are replaced by physical boundary conditions through Davies relaxation to a large-scale model (Machenhauer and Haugen, 1987; Degrauwe et al., 2012; Termonia et al., 2012). In other words, the assumption of periodic boundary conditions throughout this paper does not make the conclusions invalid for an actual atmospheric model.

In a bi-Fourier spectral space, the solution of a 2D Helmholtz equation such as (13) is trivial. Indeed, considering the bi-Fourier decomposition

of the RHS as

$$\tilde{\psi}_{\ell} = \sum_{k_x, k_y} \hat{\psi}_{k_x, k_y} \exp[i(k_x x + k_y y)], \quad (14)$$

then the solution of (13) is given by

$$\psi_{\ell} = \sum_{k_x, k_y} \frac{\hat{\psi}_{k_x, k_y}}{1 + c_{\ell}^2 \delta t^2 (k_x^2 + k_y^2)} \exp[i(k_x x + k_y y)]. \quad (15)$$

## 2.4 | Properties of the constant-coefficient 2D Helmholtz problem

Discretizing the constant-coefficient 2D Helmholtz Equation (13) with second-order finite differences on a Cartesian grid with uniform resolution $\Delta x = \Delta y$ and dimensions $n_x \times n_y$ (where $n_x$ and $n_y$ are assumed to be even for the sake of simplicity) turns this partial differential equation into a linear system

$$A \psi_{\ell} = \tilde{\psi}_{\ell}, \quad (16)$$

where $A$ is a $n_x n_y \times n_x n_y$ matrix. When solving this linear system with an iterative solver, the convergence rate of the solver is determined by the spectrum $\sigma_A$ of the matrix $A$. The eigenvalues of this matrix are

$$\lambda = 1 + \frac{1}{2} \mu_{\ell}^2 [2 - \cos(k_x \Delta x) + \cos(k_y \Delta y)], \quad (17)$$

where $k_x = 2\pi/(n_x \Delta x)$, $\{0, \dots, n_x - 1\}$ and $k_y = 2\pi/(n_y \Delta y)$, $\{0, \dots, n_y - 1\}$ are wavenumbers in both directions, and $\mu_\ell = c_\ell \Delta t / \Delta x$ is the wave Courant number corresponding to the vertical mode with index $\ell$.

The spectrum of $A$ is then $\sigma_A = [1, 1 + 2\mu_\ell^2]$. As a consequence, the convergence rate of an iterative solver is determined by the Courant number $\mu_\ell$.

## 3 | A MULTIGRID-PRECONDITIONED ITERATIVE SOLVER

As mentioned in the previous section, the solution of a 2D Helmholtz Equation (13) is obtained straightforwardly in spectral space. However, if one wants to avoid the spectral transforms, a different solver for the 2D Helmholtz equation is required. In this section, the application of an iterative solver to this problem is considered. It should be noted that the numerical problems arising from implicit timestepping in atmospheric models are typically very ill-conditioned because of the anisotropy: the horizontal
<figure><img src="image_3.png" /></figure>

resolution is usually orders of magnitude coarser than the vertical resolution, especially near the surface. For this reason, it is essential to apply a preconditioner to the system, in order to enable the (rapid) convergence of the iterative solver. Most preconditioners that have been developed so far for atmospheric models focus on this anisotropy aspect of the problem, that is, they take into account the strong vertical coupling while neglecting the off-diagonal terms that are responsible for the weaker horizontal coupling between gridpoints, thus leading to the so-called line relaxation preconditioners (Skamarock et al., 1997; Thomas et al., 1997; 2000; Qaddouri and Côté, 2003; Müller and Scheichl, 2014; Kühnlein et al., 2019, and references therein).

At this point, we come to a crucial difference – from the iterative solver point of view – between the formulation in a mass-based constant-coefficient model such as the ALADIN model and other atmospheric models that use a different vertical coordinate and/or a non-constant reference state. Whereas other models need solvers or preconditioners that deal appropriately with the strong vertical coupling, this is not an issue in a mass-based constant-coefficient model, because of the vertical eigen-mode decomposition. Some researchers (Börm and Hiptmair, 2001; Thomas et al., 2003; Dedner et al., 2016) have already explored the route of using an eigenmode decomposition in the stiffest dimension to deal with the challenge of directional anisotropy. In that research, the eigenmode decomposition is used to perform a theoretical analysis of the stiffness problem. In the present paper, this approach is taken a step further, and it will be shown how this vertical decomposition not only solves the problem of anisotropy, but also offers the opportunity to develop a highly efficient and robust multigrid preconditioner which deals appropriately with the horizontal coupling in the Helmholtz problem.

## 3.1 | The multigrid preconditioner

A good preconditioner is characterized by two properties. First, it must resemble the original system, and second it must be easy to invert. The basic idea behind multigrid preconditioning is that solving a problem numerically is easier at a lower resolution for two reasons (Sandbach et al., 2015). Besides the obvious reduction of the size of the problem, a lower resolution also leads to a better condition number for the sign-positive Helmholtz equation. In Section 2.4 it was shown how the conditioning (and thus the convergence speed of iterative methods) of a 2D Helmholtz problem is determined by the Courant number $\mu_c$, which in turn is directly related to the grid resolution.

The use of multigrid techniques in atmospheric modelling is far from new, with Fulton et al. (1986) already giving a review more than thirty years ago. Also more recently, several researchers (Chen and Sun, 2001; Buckeridge and Scheichl, 2010; Müller and Scheichl, 2014; Dedner et al., 2016; Mitchell and Müller, 2016) have investigated the virtues of multigrid-based solvers and preconditioners. The reason for this interest in multigrid techniques is that it can be shown that they are algorithmically optimal, that is, the number of operations is proportional to the number of gridpoints (Briggs et al., 2000). This contrasts with, for example, spectral methods, for which the number of operations increases more quickly than the number of gridpoints (Wedi et al., 2013).

The most basic form of a multigrid preconditioner (or solver), is based on a so-called two-grid cycle, which consists of the following steps in order to solve the system $Ax = b$:

1. A number $v_1$ of pre-smoothing steps with a stationary method, such as Jacobi or Gauss-Seidel relaxation. Denoting with $P_s$ the preconditioner describing this pre-smoother, one writes

$$x \rightarrow x + P_s^{-1}(b - Ax) \quad (v_1 \text{ times}). \quad (18)$$

2. Coarsening (restriction) of the residual $r = b - Ax$ from the fine grid to the coarse grid:

$$r^c = I^c r, \quad (19)$$

where $I^c$ is an operator describing the coarsening.

3. Solution of the residual problem at the coarser resolution:

$$\delta x^c = (A^c)^{-1} r^c, \quad (20)$$

where $A^c$ represents the discretized problem at the coarser resolution.

4. Refining (prolongation) of the coarse increment to the fine grid:

$$\delta x = I^f \delta x^c, \quad (21)$$

where $I^f$ is an operator describing the refining.

5. Adding the increment to the solution:

$$x \rightarrow x + \delta x. \quad (22)$$

6. A number $v_2$ of post-smoothing steps. As for the pre-smoothing, one writes

$$x \rightarrow x + P_s^{-1}(b - Ax) \quad (v_2 \text{ times}). \quad (23)$$
<figure><img src="image_4.png" /></figure>

In this paper, a red–black Gauss–Seidel smoother (Saad, 2003) will be used. Other smoothers have been tested, but turned out to be more costly in terms of communications versus performance. The coarse grid is assumed to have half the resolution of the fine grid. The restriction is done with full-weighting coarsening (Saad, 2003), which writes for the 1D case

$$r_j^c = \frac{1}{4}(r_{2j-1} + 2r_{2j} + r_{2j+1}). \quad (24)$$

The prolongation is done with linear interpolation. Note that these choices make the operators $I^c$ and $I^f$ adjoint (i.e., $I^c = (I^f)^T$), a property that will benefit the behaviour of the Krylov solver (see below).

If the direct solution of the coarse problem (20) is still excessively costly, a recursive strategy can be applied, where the coarse solution itself is obtained with a further multigridding approach. As such, the problem is reduced (each time improving its conditioning) until it can be solved with a few iterations of a stationary solver. In this work, this is also done with red–black Gauss–Seidel relaxations. The number of iterations at the coarsest resolution will be denoted $v_0$.

In summary, four parameters need to be chosen for the multigrid preconditioner: the number of pre-smoothing steps $v_1$, the multigrid recursion depth $d$, the number of bottom iterations $v_0$, and the number of post-smoothing steps $v_2$. The choice of a preconditioner (or its parameters) is often considered to be a difficult task, which is carried out based on intuition and/or experimentation (Mitchell and Müller, 2016). However, for our specific case of a vertically decoupled Helmholtz system of a constant-coefficient mass-based atmospheric model on a biperiodic rectangular domain, as expressed in Equation (13), it is shown in the next sections that it is possible to determine the optimal preconditioner parameters, that is, those parameter values that minimize the spectrum of the preconditioned system $AP^{-1}$, where $P$ is the matrix representation of the entire multigrid preconditioner.

## 3.2 | Spectral analysis of the multigrid preconditioner

In this section, a spectral analysis of the effect of a multigrid preconditioner for a constant-coefficient Helmholtz equation is performed. For simplicity, the one-dimensional case is studied instead of the two-dimensional case. The generalization from one to two spatial dimensions is straightforward. As indicated before, periodic boundary conditions are assumed for the Helmholtz problem. Not only is this necessary to be able to perform the spectral

analysis of this section, it also corresponds to what actually happens inside the dynamical core of the limited-area model ALADIN.

Consider the effect of the different steps (18)–(23) in a two-level multigrid preconditioner on a harmonic wave $\psi_k = e^{ikx}$ with wavenumber $k$:

* Applying the 1D Helmholtz operator $A = 1 - (c^2\Delta t^2/4)\nabla^2$, discretized with second-order centred finite differences $(\nabla^2 f_j = (f_{j-1} - 2f_j + f_{j+1})/\Delta x^2)$, only results in an amplitude change, because harmonic waves are eigenfunctions of the discretized Laplacian operator:

$$A\psi_k = \left[1 - \frac{1}{2}\mu^2(\cos k\Delta x - 1)\right]\psi_k; \quad (25)$$

* Applying a smoothing step with a red–black Gauss–Seidel iteration not only yields a wave with the same wavenumber $k$, but also one with a wavenumber $k + \pi/\Delta x$ (remapped to the wavenumber interval $2\pi/(n_x\Delta x)\{0, \dots, n_x - 1\}$ if necessary):

$$P_s^{-1}\psi_k = \frac{1}{1 + \mu^2/2} \left(1 - \frac{\mu^2 \cos k\Delta x}{4(1 + \mu^2/2)}\right)\psi_k - \frac{\mu^2 \cos k\Delta x}{4(1 + \mu^2/2)^2}\psi_{k+\pi/\Delta x}; \quad (26)$$

* Applying a full-weighting coarsening again projects a harmonic function onto itself:

$$I^c\psi_k = \frac{1 + \cos k\Delta x}{2}\psi_k; \quad (27)$$

* A linear interpolation projects the harmonic wave with wavenumber $k$ on two waves with wavenumbers $k$ and $k + \pi/\Delta x$:

$$I^f\psi_k = \frac{1 + \cos k\Delta x}{2}\psi_k + \frac{1 - \cos k\Delta x}{2}\psi_{k+\pi/\Delta x}; \quad (28)$$

* The coarser grid problem is itself a Helmholtz problem with Courant number $\mu/2$, for which harmonic waves are the eigenfunctions:

$$(A^c)^{-1}\psi_k = \frac{1}{1 - (\mu^2/8)(\cos 2k\Delta x - 1)}\psi_k. \quad (29)$$

It is observed that each individual step of the multigrid preconditioner projects a harmonic wave with wavenumber $k$ onto two waves with respective wavenumbers $k$ and $k + \pi/\Delta x$. Keeping in mind that $\exp[i(k + 2\pi/\Delta x)x] \equiv \exp(ikx)$, it follows that also the wave with wavenumber $k + \pi/\Delta x$ is projected onto two waves with respective wavenumbers $k$ and $k + \pi/\Delta x$. Since this property holds
<figure><img src="image_5.png" /></figure>

for every single step (smoothing, restriction, coarse solver, prolongation) in the multigrid preconditioning procedure, it also holds for the aggregated operator $\mathbf{AP}^{-1}$. One can thus formally write

$$ \mathbf{AP}^{-1}\psi_k = \alpha_{11}\psi_k + \alpha_{21}\psi_{k+\pi/\Delta x}, \quad (30) $$

$$ \mathbf{AP}^{-1}\psi_{k+\pi/\Delta x} = \alpha_{12}\psi_k + \alpha_{22}\psi_{k+\pi/\Delta x}, \quad (31) $$

where the coefficients $\alpha_{11}$, $\alpha_{12}$, $\alpha_{21}$ and $\alpha_{22}$ can be, in principle, analytically determined from Equations (25)–(29).

In other words, any linear combination of two waves $\exp(ikx)$ and $\exp[i(k + \pi/\Delta x)x]$ will be projected by the operator $\mathbf{AP}^{-1}$ onto a (different) linear combination of those same two harmonic waves. As a consequence, the eigenvectors of the preconditioned system $\mathbf{AP}^{-1}$ can be expressed as linear combinations of two waves with wavenumbers $k$ and $k + \pi/\Delta x$. Furthermore, applying a Rayleigh–Ritz method (Saad, 2011) with a vector basis consisting of these two waves will yield the exact corresponding eigenvalues. Indeed, considering the orthonormal vector basis $\mathbf{V}_k = [\psi_k \ \psi_{k+\pi/\Delta x}]$, one can write

$$ \mathbf{V}_k^\mathrm{T}\mathbf{AP}^{-1}\mathbf{V}_k = \alpha_k, \quad (32) $$

where the matrix $\alpha_k$ is given by

$$ \alpha_k = \begin{pmatrix} \alpha_{11} & \alpha_{12} \\ \alpha_{21} & \alpha_{22} \end{pmatrix}. \quad (33) $$

The eigenvalues of $\mathbf{AP}^{-1}$ corresponding to the eigenvectors that lie in the space spanned by $\mathbf{V}_k$ are exactly equal to the eigenvalues of the matrix $\alpha_k$.

For every set of parameters describing the system and the preconditioner, it is then possible to find the wavenumber $k_1$ that maximizes the largest eigenvalue of $\alpha_k$, and to find another wavenumber $k_2$ that minimizes the smallest eigenvalue of $\alpha_k$. These two extreme eigenvalues will be the extreme eigenvalues of the full system $\mathbf{AP}^{-1}$, and they are independent of the grid dimensions. This implies grid-independent convergence of the iterative solver.

For a multigrid method with three grid levels, the same reasoning can be followed, except that the Rayleigh–Ritz vector basis should now consist of four harmonic waves with wavenumbers $k + (\pi/2\Delta x)\{0, 1, 2, 3\}$, because of the smoothing, reduction and prolongation at the second grid level, and the matrix $\alpha_k$ will become a $4 \times 4$ matrix. The same methodology can be extended to more multigrid levels.

Generalizing from the 1D case to the 2D case is straightforward, although performing the calculations analytically quickly becomes intractable, especially the calculation of the eigenvalues of the matrix $\alpha_k$. However, what is

<figure><img src="image_6.png" /><figcaption>FIGURE 1 Minimum (solid blue) and maximum (dashed green) eigenvalue of the matrix $\mathbf{AP}^{-1}$ as a function of the Courant number. A represents the 2D Helmholtz system as defined in Equation (13); P represents a multigrid preconditioner with parameters $d=2$, $v_0 = 2$, $v_1 = 1$, $v_2 = 1$. The convergence rates of a preconditioned Richardson solver (short dashed red) and a preconditioned Krylov solver GCR(k=1) (dotted cyan) are also shown. [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

important is that the eigenvectors of $\mathbf{AP}^{-1}$ can still be expressed as linear combinations of a limited set of harmonic waves, and that the corresponding eigenvalues can still be obtained (numerically) from a low-dimensional Rayleigh–Ritz calculation, independent of the grid dimensions.

To summarize, it is possible to determine the minimum eigenvalue $\lambda_{\min}$ and maximum eigenvalue $\lambda_{\max}$ of the system $\mathbf{AP}^{-1}$ for a given Courant number $\mu$ and given preconditioner parameter values $d$, $v_0$, $v_1$ and $v_2$. These extreme eigenvalues are not affected by the grid dimensions, so they can be calculated once and for all. Moreover, they can be efficiently determined with a low-dimensional Rayleigh–Ritz calculation. This is illustrated in Figure 1, which shows the minimum and maximum eigenvalue (for the 2D case) as a function of the Courant number, for given multigrid preconditioner parameter values $d=2$, $v_0 = 2$, $v_1 = 1$, $v_2 = 1$. The points where the slope of the curves changes discontinuously correspond to a swapping of the eigenmode yielding the minimum or maximum eigenvalue.

## 3.3 | Convergence rate of multigrid preconditioned solvers

Once the extreme eigenvalues $\lambda_{\min}$ and $\lambda_{\max}$ of the preconditioned system have been calculated for given preconditioner parameters, it becomes possible to estimate the convergence rate of an iterative solver for the 2D
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>FIGURE 2 Comparison between predicted (solid lines) and measured (symbols) convergence of a Richardson solver (blue) and a Krylov solver (green) applied to a 2D Helmholtz problem with Courant number $\mu = 1.44$. The Krylov solver is a truncated generalized conjugate residual solver (GCR(k=1)). Both solvers are preconditioned with a multigrid preconditioner with parameters $d=1$, $v_0 = 2$, $v_1 = 1$ and $v_2 = 1$ [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

Helmholtz problem. Indeed, when considering solving the system with Richardson iterations (Saad, 2003), the convergence rate is given by

$$ \rho_R(\mu, d, v_0, v_1, v_2) = \max\{|1 - \lambda_{\min}|, |1 - \lambda_{\max}|\}. \quad (34) $$

A strict upper bound for the convergence rate of a Krylov solver is given by (Saad, 2003)

$$ \rho_K(\mu, d, v_0, v_1, v_2) = \frac{\sqrt{\lambda_{\max}/\lambda_{\min}} - 1}{\sqrt{\lambda_{\max}/\lambda_{\min}} + 1}. \quad (35) $$

These convergence rates are shown in Figure 1.

To verify if the convergence rates, as estimated from the extreme eigenvalues, correspond to the actual convergence rates, Figure 2 shows a comparison between the predicted and measured convergence for a 2D Helmholtz system solved with Richardson iterations and with a Krylov solver. This confirms that the estimated convergence rates are very accurate.

When atmospheric models are used for operational weather forecasting, it is necessary that the model run finishes in time. The importance of being able to determine the convergence rate of the iterative solver *a priori* (i.e., before actually starting the forecast run) therefore cannot be overestimated. Indeed, this knowledge allows us to determine the required number of iterations to reach a given precision, which in turn leads to very predictable and stable runtimes.

An additional advantage of knowing the required number of iterations *a priori* is that it is not even necessary

to test for convergence during the iterative solving. This is relevant because convergence tests are usually based on residual or increment norms, which require a (admittedly very small-volume) global communication (reduction operation). On massively parallel machines, the overhead associated with such a global reduction operation may become non-negligible (Zheng and Marguinaud, 2018), as it requires a synchronization between all processes. Of course, Krylov methods involve taking scalar products, which require the same type of global reduction operations and the resulting synchronization. But when using a preconditioned Richardson solver, it thus becomes possible to develop an iterative Helmholtz solver that does not require any global synchronization.

## 3.4 | Selection of optimal multigrid preconditioner parameters

Besides a guaranteed robustness of the solver, the fact that the convergence speed of the iterative Helmholtz solver can be determined *a priori*, without actually running the iterative solver, has another advantage: it becomes possible to compare different values for the preconditioner parameters, and to pick those parameters that lead to the fastest convergence.

Figure 3 shows the convergence speed (expressed as the number of digits gained per iteration) as a function of the Courant number when using a Krylov solver with a multigrid preconditioner with various depth values. It is observed that multigridding is not necessary ($d=0$) for small values of the Courant number. For higher values of the Courant number, however, the convergence speed without multigridding becomes quickly worse, and multigridding accelerates the solver significantly. A depth of $d=1$ gives fast convergence up to a Courant number of about 2.5, for even larger values, a higher value for the depth is required.

These observations make sense if one considers the physical meaning of the Courant number (Sandbach et al., 2015). A low value means that disturbances travel slowly in the system. Therefore a stationary preconditioner such as red-black Gauss-Seidel relaxation can be expected to perform well. However, for a system with a higher Courant number, the region of influence of a disturbance during a single timestep becomes larger, and the solver will require that information is exchanged over a larger distance in order to obtain an accurate solution. Exchanging information over such a larger distance is what multigrid techniques achieve through considering the problem at a coarser resolution.

Of course, when comparing different preconditioner settings, their relative cost should be taken into account.
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>**FIGURE 3** Convergence speed of a multigrid-preconditioned Krylov solver for the 2D Helmholtz system, as a function of the Courant number. The multigrid preconditioner has a depth $d=0$ (solid blue), $d=1$ (dashed green), $d=2$ (short-dashed red), $d=3$ (dotted cyan) or $d=4$ (dash-dotted magenta). The other parameters are $v_0 = 2$, $v_1 = 1$, $v_2 = 1$ [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

It is to be expected that increasing the number of pre- or post-smoothing steps will lead to faster convergence, but also to a higher cost in terms of computations and communications. If one wants to determine the preconditioner settings that minimize the number of communications, one could look at the number of digits gained *per halo exchange*, where the required communications for the different steps in the preconditioner (pre-smoothing, restriction, coarse solve, prolongation and post-smoothing) are taken into account. Figure 4 shows this measure, again as a function of the Courant number, for various parameter combinations. As indicated on the figure, the optimal parameter values depend on the Courant number, where more multigrid levels are necessary for higher Courant values. The optimal number of pre- and post-smoothing turns out to be $v_1 = v_2 = 1$ throughout the range of Courant numbers. The optimal number of relaxation iterations at the coarsest resolution $v_0$ varies between 1 and 3. It is also remarkable that the increased convergence speed of multigridding with a single level ($d=1$) does not outweigh the cost of the additional communications: for Courant numbers around $\mu = 2$, it turns out to be better to use a cheaper but slower-converging Gauss-Seidel preconditioner than a faster but more expensive multigrid preconditioner.

**3.5 | Application in a decoupled 3D system**

As explained in Section 2.2, the vertical decoupling of the 3D Helmholtz system leads to a set of 2D Helmholtz

<figure><img src="image_11.png" /><figcaption>**FIGURE 4** Convergence of a multigrid-preconditioned Krylov solver applied to a 2D Helmholtz problem. The grey lines represent different combinations of preconditioner parameters $\{d, v_0, v_1, v_2\} \in [0:4] \times [1:4] \times [0:2] \times [0:2]$. The optimal parameters (depending on the Courant number) are indicated by the coloured lines and the associated labels</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**FIGURE 5** Courant numbers of vertical eigenmodes for a system with horizontal resolution 1 km, timestep 30 s, and a vertical discretization in 90 levels [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

systems, each with its own wave Courant number. As an example, Figure 5 shows the Courant numbers of the vertical eigenmodes for a system with horizontal resolution $\Delta x = 1$ km, a timestep of $\Delta t = 30$ s, and a vertical discretization with 90 levels. It is noticed that only a few vertical eigenmodes have a large Courant number.

When applying the multigrid preconditioner described above, it is possible to choose the preconditioner parameters for each vertical eigenmode separately. As such, only the few eigenmodes with large Courant numbers will require many iterations or a large value for the multigridding depth. For the other eigenmodes, fewer and cheaper iterations will suffice. This observation was already made
by Qaddouri and Côté (2003). As such, the number of computations and communications can be significantly reduced with respect to an approach where the same preconditioner parameters are used for all vertical eigenmodes. This strategy can only be applied if a vertical decomposition is possible, which is a consequence of the constant-coefficient reference state.

To assess the significance of this additional advantage of the vertical decoupling of the 3D Helmholtz system, consider the 90-level vertical discretization of Figure 5. The highest Courant number in this system is $\mu_{max} = 11.24$. The optimal preconditioner settings (in terms of communications) corresponding to this Courant number can be read from Figure 4, and turns out to be $d = 3, v_0 = 2, v_1 = 1$ and $v_2 = 1$. The number of halo exchanges for a single iteration of the Krylov solver with these preconditioner settings is 5.5. The required number of iterations to reach a precision of $10^{-8}$ can be deduced from Figure 3 and is equal to 6. Therefore, solving the Helmholtz system using these preconditioner settings for all 90 vertical eigenmodes would require a total of $5.5 \times 6 \times 90 = 2970$ halo exchanges. On the other hand, by using optimized preconditioner settings for each vertical eigenmode separately, the total number of halo exchanges is reduced to 548. This represents a saving in the communication volume by a factor of more than 5.

# 4 | SCALABILITY OF THE MULTIGRID-PRECONDITIONED ITERATIVE SOLVERS

## 4.1 | Experiment set-up

The scalability of three solvers for the Helmholtz equation in a constant-coefficient mass-based non-hydrostatic atmospheric model is investigated. The first solver is the spectral solver that is currently used in the ALADIN NWP model, as described in Section 2.3. The second solver is the truncated generalized conjugate residual method (GCR(k)), as introduced by Eisenstat *et al.* (1983). This Krylov solver is commonly used in atmospheric models (Skamarock *et al.*, 1997; Thomas *et al.*, 2003; Smolarkiewicz and Szmelter, 2011) because, as opposed to the well-known conjugate gradient algorithm, GCR(k) does not require the linear system to be symmetric. The reason this solver is used in this work is that, even though the Helmholtz system (11) itself is symmetric, the multigrid preconditioner is not. The untruncated version of the GCR algorithm requires the storage of Krylov base vectors for successive iterations. To avoid excessive memory consumption, a truncated version of the algorithm was

developed where the number of base vectors is limited to a fixed number *k*. Because our preconditioned system will be close to symmetric, it is to be expected that the convergence speed is quite close to that of the conjugate gradient algorithm (Eisenstat *et al.*, 1983), even when the truncation parameter is set to *k* = 1. The third solver under investigation is a Richardson solver. As explained before, a specific property of this solver is that it does not require any kind of global communications or system-wide synchronization. Both non-spectral solvers are fitted with the preconditioner described in the previous section. The parameters of the preconditioner are chosen for each vertical eigenmode separately, so as to minimize the number of communications. As for the required target precision, several values are found in literature (Buckeridge and Scheichl, 2010; Müller and Scheichl, 2014; Mitchell and Müller, 2016). For this study, a reduction of the residual error by $10^8$ is chosen, which is more than sufficient for numerical weather prediction purposes.

Rather than testing the different solvers inside a fully-fledged atmospheric model, the scalability experiments in this section are performed with a dedicated computer program that focuses on the solution of the implicit system as given in Equations (6)–(10). More concretely, this program carries out the following steps:

1. [Only for the spectral solver] forward spectral transform of prognostic variables $D, d, T', \hat{q}$ and $\pi'$ at time instant $t$;

2. Calculation of the RHS terms in Equations (6)–(10): $\tilde{D}, \tilde{d}, \tilde{T}, \tilde{q}$ and $\tilde{\pi}_s$;

3. Calculation of the RHS term of the 3D Helmholtz Equation (11): $D^{**}$;

4. Projection of $D^{**}$ on vertical eigenmodes: $\tilde{\psi} = V_B^{-1}D^{**}$;

5. Solution of a 2D Helmholtz Equation 13 for each vertical eigenmode, with the spectral solver, the Krylov solver, or the Richardson solver;

6. Inverse projection from vertical eigenmodes: $D = V_B\psi$;

7. Substitution of $D$ in Equations (7)–(10) to obtain the other prognostic variables $d, T', \hat{q}$ and $\pi'$ at time instant $t + \Delta t$;

8. [Only for the spectral solver] inverse spectral transform of prognostic variables $D, d, T', \hat{q}$ and $\pi'$.

It is worth mentioning that the projection to and from vertical eigenmodes can be done columnwise, so steps 4 and 6 do not require any horizontal communications.

The experiments will run with a vertical discretization in 96 levels, and a ratio of horizontal resolution over timestep of $\Delta x/\Delta t = 1 \text{ km} \cdot (30 \text{ s})^{-1}$. The values of the variables at time instant $t$ are random numbers, and the initial guess for the iterative solvers is taken equal to zero.
## 4.2 | Scalability results

The scalability is tested by means of a weak scalability experiment where the problem size increases proportionally to the processor count. This means that, as the number of processors increases, the spatial resolution is increased, such that the total domain extent remains constant, and such that each compute core is always attributed a 32 × 32 horizontal subdomain. The time step $\Delta t$ is scaled proportionally to the horizontal resolution in order to maintain the ratio $\Delta x/\Delta t = 1\,\text{km} \cdot (30\,\text{s})^{-1}$. In the case of perfect scalability, the run time of these experiments should not increase with increasing processor count.

Runs are performed on two machines. The first machine (beaufix) is a Bull cluster at Météo-France, where the compute nodes consist of two Intel Xeon E5-2698 v4 processors with 20 physical cores each. The interconnect network is an Infiniband-based non-blocking (1:1) fat tree topology. The second machine (cca) is a Cray cluster at ECMWF, consisting of Intel Xeon E5-2695 v4 processors with 18 physical cores each. This cluster uses a Cray Aries interconnect in a dragonfly topology.

Figure 6 shows the run times on beaufix. Both pure MPI and hybrid MPI-OpenMP runs were performed. It is observed that, at low processor counts, the three solvers perform quite comparably. However, at higher processor counts, the scalability of the GCR(k) solver and the Richardson solver is much better than that of the spectral solver. The performance of the two non-spectral solvers is almost identical. Although the GCR(k) solver requires fewer iterations to reach the target precision, the cost per iteration is lower for the Richardson solver. The same observations are made from the runs on cca, as shown in Figure 7.

Some important remarks are in place when considering these results. First, these experiments do not measure the scalability of an entire atmospheric model, but only that of the Helmholtz solver. Since the other parts of an atmospheric model usually scale much better, for example the physics parametrizations that are performed in grid-columns, the overall scalability will be much better than is shown here. Concretely, in the operational runs performed with ALADIN at Météo-France, the spectral Helmholtz solver only accounts for about 10% of the total runtime.

A second remark is that the two machines used for the experiments were purchased in order to run spectral models effectively. This means that they were designed to provide good all-to-all communications bandwidth (which, of course, comes at a financial cost). Since non-spectral solvers do not rely on all-to-all communications, the design (and cost) of future machines could be adapted accordingly.

<figure><img src="image_13.png" /><figcaption>**FIGURE 6** Wall-clock times of the weak-scalability experiment on beaufix: spectral solver (solid blue), GCR(k) solver (dashed green) and Richardson solver (short-dashed red). Circle symbols denote pure MPI experiments, and cross symbols denote hybrid MPI-OpenMP experiments with 20 threads [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**FIGURE 7** Wall-clock times of the weak-scalability experiment on cca: spectral solver (solid blue), GCR(k) solver (dashed green) and Richardson solver (short-dashed red). Circle markers denote pure MPI experiments; cross markers denote hybrid MPI-OpenMP experiments with 18 threads [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

## 5 | CONCLUSIONS

This paper examines the scalability of atmospheric models with a semi-implicit timestepping scheme using a constant-coefficient reference state and a mass-based vertical coordinate. The Helmholtz problem that appears in such models is commonly solved with a spectral transform method. However, the scalability of spectral transforms on future massively-parallel machines may become problematic, because they involve data-rich domain-wide communications. This paper investigates the viability of an iterative non-spectral solver as an alternative to the existing spectral Helmholtz solver of the semi-implicit ALADIN
model, without touching the other features of this model,
most notably its mass-based vertical coordinate system.

The application of iterative solvers in atmospheric
modelling is far from new. The specificity of the work
presented in this paper is that the system under consid-
eration employs a constant-coefficient reference state for
the semi-implicit scheme. This allows for a vertical decou-
pling of the three-dimensional Helmholtz problem into a
set of two-dimensional Helmholtz problems. A direct con-
sequence is that the anisotropy due to the difference in
horizontal and vertical resolution, which is an aspect that
needs proper attention in many semi-implicit atmospheric
models, is not an issue for our constant-coefficient model.
As a result, instead of focusing on the strong vertical cou-
pling, the multigrid-based preconditioner considered in
this paper appropriately treats the horizontal coupling in
the Helmholtz system. In theory, the techniques presented
in this paper could also be applied in an atmospheric
model with a different vertical coordinate, as long as the
model uses a constant-coefficient reference state.

A possible concern with iterative solvers is that con-
vergence occasionally may be slow, for example under
specific weather conditions. For an atmospheric model
such as ALADIN, which is used for operational weather
forecasting, this would pose a significant problem. How-
ever, it is shown in this paper that, for the specific case of
a constant-coefficient Helmholtz problem on a Cartesian
grid (exactly as it arises in the ALADIN dynamical core),
the convergence rate of a multigrid-preconditioned iter-
ative solver is guaranteed, and can even be determined a
priori with a low-dimensional Rayleigh-Ritz calculation.
In other words, the number of iterations required to reach
a predefined precision is known before even starting the
model run, which makes the runtime of the model per-
fectly predictable. As such, the multigrid-preconditioned
Helmholtz solver presented in this paper is, in com-
bination with the constant-coefficient, mass-based
semi-implicit dynamics of ALADIN, extremely well suited
for use in an operational context.

This paper also shows how the *a priori* knowledge
of the convergence rate of the iterative Helmholtz solver
can be used to select the parameters of the multigrid
preconditioner in an optimal way. By choosing optimal
parameters for each vertical eigenmode separately, the
efficiency of the iterative solver is greatly improved. Two
multigrid-preconditioned iterative solvers were tested by
means of weak scalability experiments: a Krylov solver
(GCR(k)) and a Richardson solver. These experiments con-
firm the superior scalability of both iterative solvers with
respect to the spectral solver.

A major strength of non-spectral solvers lies in their
ability to deal with problems that are not necessarily
horizontally homogeneous (Smolarkiewicz and Margolin,

1994; Skamarock *et al.*, 1997). For instance, treating in
an implicit way the orographic terms that typically result
from the use of terrain-following vertical coordinates is
simply not possible with a spectral solver, while it is fea-
sible with an iterative gridpoint solver. Thus, by restrict-
ing the work in this paper to the application of iterative
solvers on a constant-coefficient problem, we may not yet
fully exploit the power of these solvers. However, as indi-
cated by many researchers (Barros, 1991; Thomas *et al.*,
2003; Dedner *et al.*, 2016; Mitchell and Müller, 2016), a
constant-coefficient solver may act as an excellent precon-
ditioner for a more complicated system where the coeffi-
cients are not horizontally constant. For example, Thomas
*et al.* (2003) show that a preconditioner without orog-
raphy performs well in the case of complex orography,
despite the substantial departure of the preconditioner
from the governing Helmholtz operator. It is along this line
that our research will continue for the Helmholtz solver
of the ALADIN model. It can be expected that the con-
vergence with a constant-coefficient preconditioner for a
non-constant-coefficient problem may be slower.

Although the results in this paper are quite promising,
significant challenges still need to be overcome before the
ALADIN model can be run entirely without spectral trans-
forms. Indeed, not only is the Helmholtz system solved in
spectral space, but also the calculation of spatial deriva-
tives and high-order diffusion are done spectrally. While
these operations could be done in gridpoint space, for
example by means of finite differences, the effect on the
accuracy should be closely investigated. The second-order
finite differences used in this paper may be insufficient in
this respect. Nevertheless, the calculations of this paper
can be repeated for higher-order finite-difference opera-
tors. While the numerical results (e.g., the convergence
rates and the required number of iterations) may change
for higher-order discretizations, the overall conclusions of
this paper will remain valid.

ACKNOWLEDGMENTS

This work was supported by the European UnionáĂŹs
Horizon 2020 research and innovation programme (GA
67162). The authors wish to thank the two anonymous
reviewers who helped to improve the manuscript.

ORCID

Daan Degrauwe https://orcid.org/0000-0002-7340-0510

Fabrice Voitus https://orcid.org/0000-0003-4920-3063

REFERENCES

ALADIN International Team (1997) The ALADIN project: mesoscale
modelling seen as a basic tool for weather forecasting and atmo-
spheric research. *WMO Bulletin*, 46(4), 317–324
Barros, S.R.M. (1991) Multigrid methods for two- and three-dimensional Poisson-type equations on the sphere. *Journal of Computational Physics*, 92(2), 313-348. https://doi.org/10.1016/0021-9991(91)90213-5

Barros, S.R.M., Dent, D., Isaksen, L., Robinson, G., Mozdzynski, G. and Wollenweber, F. (1995) The IFS model: a parallel production weather code. *Parallel Computing*, 21(10), 1621-1638. https://doi.org/10.1016/0167-8191(96)80002-0

Bauer, P., Thorpe, A.J. and Brunet, G. (2015) The quiet revolution of numerical weather prediction. *Nature*, 525, 47-55. https://doi.org/10.1038/nature14956

Bénard, P., Laprise, R., Vivoda, J. and Smolíková, P. (2004) Stability of leapfrog constant-coefficients semi-implicit schemes for the fully elastic system of Euler equations: flat-terrain case. *Monthly Weather Review*, 132(5), 1306-1318. https://doi.org/10.1175/1520-0493(2004)132<1306:SOLCSS>2.0.CO;2

Bénard, P., Vivoda, J., Mašek, J., Smolíková, P., Yessad, K., Smith, C., Brožková, R. and Geleyn, J.-F. (2010) Dynamical kernel of the Aladin-NH spectral limited-area model: revised formulation and sensitivity experiments. *Quarterly Journal of the Royal Meteorological Society*, 136, 155-169. https://doi.org/10.1002/qj.522

Börm, S. and Hiptmair, R. (2001) Analysis of tensor product multigrid. *Numerical Algorithms*, 26(3), 219-234. https://doi.org/10.1023/A:1016686408271

Briggs, W.L., Van Emden, H. and McCormick, S.F. (2000) *A Multigrid Tutorial* (2nd Edition). SIAM, Philadelphia, PA.

Bubnová, R., Hello, G., Bénard, P. and Geleyn, J.-F. (1995) Integration of the fully elastic equations cast in the hydrostatic pressure terrain-following coordinate in the framework of the ARPEGE/ALADIN NWP system. *Monthly Weather Review*, 123(2), 515-535. https://doi.org/10.1175/1520-0493(1995)123<0515:TOTFEE>2.0.CO;2

Buckeridge, S. and Scheichl, R. (2010) Parallel geometric multigrid for global weather prediction. *Numerical Linear Algebra with Applications*, 17(2-3), 325-342. https://doi.org/10.1002/nla.699

Caluwaerts, S., Degrauwe, D., Termonia, P., Voitus, F., Bénard, P. and Geleyn, J.-F. (2015) Importance of temporal symmetry in spatial discretization for geostrophic adjustment in semi-implicit Z-grid schemes. *Quarterly Journal of the Royal Meteorological Society*, 141, 128-138. https://doi.org/10.1002/qj.2344

Catry, B., Geleyn, J.-F., Tudor, M., Bénard, P. and Trojáková, A. (2007) Flux-conservative thermodynamic equations in a mass-weighted framework. *Tellus A: Dynamic Meteorology and Oceanography*, 59(1), 71-79. https://doi.org/10.1111/j.1600-0870.2006.00212.x

Chen, S.-H. and Sun, W.-Y. (2001) Application of the multigrid method and a flexible hybrid coordinate in a nonhydrostatic model. *Monthly Weather Review*, 129(11), 2660-2676. https://doi.org/10.1175/1520-0493(2001)129<2660:AOTMMA>2.0.CO;2

Courtier, P. and Geleyn, J.-F. (1988) A global numerical weather prediction model with variable resolution: application to the shallow-water equations. *Quarterly Journal of the Royal Meteorological Society*, 114, 1321-1346. https://doi.org/10.1002/qj.49711448309

Courtier, P., Freydier, C., Geleyn, J.-F., Rabier, F. and Rochas, M. (1991). The ARPEGE project at Météo-France, pp. 193-231 in Proceedings of 1991 ECMWF Seminar on Numerical Methods in Atmospheric Models. ECMWF, Reading, UK.

Dedner, A., Müller, E. and Scheichl, R. (2016) Efficient multigrid preconditioners for atmospheric flow simulations at high aspect

ratio. *International Journal for Numerical Methods in Fluids*, 80(1), 76-102. https://doi.org/10.1002/fld.4072

Degrauwe, D., Caluwaerts, S., Voitus, F., Hamdi, R. and Termonia, P. (2012) Application of Boyd's periodization and relaxation method in a spectral atmospheric limited-area model. Part II: accuracy analysis and detailed study of the operational impact. *Monthly Weather Review*, 140(10), 3149-3162. https://doi.org/10.1175/MWR-D-12-00032.1

Degrauwe, D., Seity, Y., Bouyssel, F. and Termonia, P. (2016) Generalization and application of the flux-conservative thermodynamic equations in the AROME model of the ALADIN system. *Geoscientific Model Development*, 9(6), 2129-2142. https://doi.org/10.5194/gmd-9-2129-2016

Eisenstat, S., Elman, H. and Schultz, M. (1983) Variational iterative methods for nonsymmetric systems of linear equations. *SIAM Journal on Numerical Analysis*, 20(2), 345-357. https://doi.org/10.1137/0720023

Fulton, S.R., Ciesielski, P.E. and Schubert, W.H. (1986) Multigrid methods for elliptic problems: a review. *Monthly Weather Review*, 114(5), 943-959. https://doi.org/10.1175/1520-0493(1986)114<0943:MMFEPA>2.0.CO;2

Husain, S.Z., Girard, C., Qaddouri, A. and Plante, A. (2019) A new dynamical, Core of the Global Environmental Multiscale (GEM) model with a height-based terrain-following vertical coordinate. *Monthly Weather Review*, 147(7), 2555-2578. https://doi.org/10.1175/MWR-D-18-0438.1

Ikawa, M. (1988) Comparison of some schemes for nonhydrostatic models with orography. *Journal of the Meteorological Society of Japan. Series II*, 66(5), 753-776. https://doi.org/10.2151/jmsj1965.66.5_753

Kühnlein, C., Deconinck, W., Klein, R., Malardel, S., Piotrowski, Z.P., Smolarkiewicz, P.K., Szmelter, J. and Wedi, N.P. (2019) Fvm 1.0: a nonhydrostatic finite-volume dynamical core for the IFS. *Geoscientific Model Development*, 12(2), 651-676. https://doi.org/10.5194/gmd-12-651-2019

Machenhauer, B. and Haugen, J.E. (1987). Test of a spectral limited-area shallow-water model with time-dependent lateral boundary conditions and combined normal mode/semi-Lagrangian time integration schemes, pp. 361-377 in ECMWF Workshop Proceedings: Techniques for Horizontal Discretization in Numerical Weather Prediction Models. ECMWF, Reading, UK.

Michalakes, J., Govett, M., Benson, R., Black, T., Juang, H.-M.H., Reinecke, A., Skamarock, W., Duda, M., Henderson, T., Madden, P., Mozdzynski, G. and Vasic, R. (2015). AVEC report: NGGPS level-1 benchmarks and software evaluation. https://repository.library.noaa.gov/view/noaa/18654; accessed 12 September 2020.

Mitchell, L. and Müller, E.H. (2016) High level implementation of geometric multigrid solvers for finite-element problems: applications in atmospheric modelling. *Journal of Computational Physics*, 327, 1-18. https://doi.org/10.1016/j.jcp.2016.09.037

Müller, A., Deconinck, W., Kühnlein, C., Mengaldo, G., Lange, M., Wedi, N., Bauer, P., Smolarkiewicz, P.K., Diamantakis, M., Lock, S.-J., Hamrud, M., Saarinen, S., Mozdzynski, G., Thiemert, D., Glinton, M., Bénard, P., Voitus, F., Colavolpe, C., Marguinaud, P., Zheng, Y., Van Bever, J., Degrauwe, D., Smet, G., Termonia, P., Nielsen, K.P., Sass, B.H., Poulsen, J.W., Berg, P., Osuna, C., Fuhrer, O., Clement, V., Baldauf, M., Gillard, M., Szmelter, J., O'Brien, E., McKinstry, A., Robinson, O., Shukla, P., Lysaght, M.,
Kulczewski, M., Ciznicki, M., Piatek, W., Ciesielski, S., Blazewicz, M., Kurowski, K., Procyk, M., Spychala, P., Bosak, B., Piotrowski, Z.P., Wyszogrodzki, A., Raffin, E., Mazauric, C., Guibert, D., Douriez, L., Vigouroux, X., Gray, A., Messmer, P., Macfaden, A.J. and New, N. (2019) The ESCAPE project: energy-efficient scalable algorithms for weather prediction at exascale. Geoscientific Model Development, 12(10), 4425–4441. https://doi.org/10.5194/gmd-12-4425-2019

Müller, E.H. and Scheichl, R. (2014) Massively parallel solvers for elliptic partial differential equations in numerical weather and climate prediction. Quarterly Journal of the Royal Meteorological Society, 140, 2608–2624. https://doi.org/10.1002/qj.2327

Qaddouri, A. and Côté, J. (2003). Preconditioning for an iterative elliptic solver on a vector processor, pp. 343–353 in High Performance Computing for Computational Science — VECPAR 2002. Palma, J.M.L.M., Sousa, A.A., Dongarra, J., Hernández, V. (eds), Springer, Berlin.

Robert, A. (1981) A stable numerical integration scheme for the primitive meteorological equations. Atmosphere-Ocean, 19(1), 35–46. https://doi.org/10.1080/07055900.1981.9649098

Robert, A. (1982) A semi-Lagrangian and semi-implicit numerical integration scheme for the primitive meteorological equations. Journal of the Meteorological Society of Japan, 60(1), 319–325. https://doi.org/10.2151/jmsj1965.60.1_319

Saad, Y. (2003) Iterative Methods for Sparse Linear Systems. SIAM, Philadelphia, PA.

Saad, Y. (2011) Numerical Methods for Large Eigenvalue Problems. SIAM, Philadelphia, PA.

Saito, K., Fujita, T., Yamada, Y., Ishida, J.-I., Kumagai, Y., Aranami, K., Ohmori, S., Nagasawa, R., Kumagai, S., Muroi, C., Kato, T., Eito, H. and Yamazaki, Y. (2006) The operational JMA non-hydrostatic mesoscale model. Monthly Weather Review, 134(4), 1266–1298. https://doi.org/10.1175/MWR3120.1

Sandbach, S., Thuburn, J., Vassilev, D. and Duda, M.G. (2015) A semi-implicit version of the MPAS-atmosphere dynamical core. Monthly Weather Review, 143(9), 3838–3855. https://doi.org/10.1175/MWR-D-15-0059.1

Skamarock, W.C., Smolarkiewicz, P.K. and Klemp, J.B. (1997) Preconditioned conjugate-residual solvers for Helmholtz equations in nonhydrostatic models. Monthly Weather Review, 125(4), 587–599. https://doi.org/10.1175/1520-0493(1997)125<0587:PCRFH>2.0.CO;2

Smolarkiewicz, P.K. and Margolin, L.G. (1994) Variational Elliptic Solver for Atmospheric Applications. Technical Report LA-12712-MS, Los Alamos National Laboratory, Santa Fe, NM.

Smolarkiewicz, P.K. and Szmelter, J. (2011) A nonhydrostatic unstructured-mesh soundproof model for simulation of internal gravity waves. Acta Geophysica, 59(6). https://doi.org/10.2478/s11600-011-0043-z

Steppeler, J., Hess, R., Schättler, U. and Bonaventura, L. (2003) Review of numerical methods for nonhydrostatic weather prediction models. Meteorology and Atmospheric Physics, 82(1), 287–301. https://doi.org/10.1007/s00703-001-0593-8

Termonia, P., Voitus, F., Degrauwe, D., Caluwaerts, S. and Hamdi, R. (2012) Application of Boyd’s periodization and relaxation method in a spectral atmospheric limited-area model. Part I: implementation and reproducibility tests. Monthly Weather Review, 140(10), 3137–3148. https://doi.org/10.1175/MWR-D-12-00033.1

Termonia, P., Fischer, C., Bazile, E., Bouyssel, F., Brožková, R., Bénard, P., Bochenek, B., Degrauwe, D., Derková, M., El Khatib,

R., Hamdi, R., Mašek, J., Pottier, P., Pristov, N., Seity, Y., Smolíková, P., Španiel, O., Tudor, M., Wang, Y., Wittmann, C. and Joly, A. (2018) The ALADIN system and its canonical model configurations AROME cy41t1 and ALARO cy40t1. Geoscientific Model Development, 11(1), 257–281. https://doi.org/10.5194/gmd-11-257-2018

Thomas, S., Girard, C., Doms, G. and Schättler, U. (2000) Semi-implicit scheme for the DWD Lokal-Modell. Meteorology and Atmospheric Physics, 73(1), 105–125. https://doi.org/10.1007/s007030050068

Thomas, S.J., Malevsky, A.V., Desgagné, M., Benoit, R., Pellerin, P. and Valin, M. (1997) Massively parallel implementation of the mesoscale compressible community model. Parallel Computing, 23(14), 2143–2160. https://doi.org/10.1016/S0167-8191(97)00105-1

Thomas, S.J., Hacker, J.P., Smolarkiewicz, P.K. and Stull, R.B. (2003) Spectral preconditioners for nonhydrostatic atmospheric models. Monthly Weather Review, 131(10), 2464–2478. https://doi.org/10.1175/1520-0493(2003)131<2464:SPFNAM>2.0.CO;2

Wedi, N., Bauer, P., Deconinck, W., Diamantakis, M., Hamrud, M., Kühnlein, C., Malardel, S., Mogensen, K., Mozdzynski, G. and Smolarkiewicz, P. (2015). The modelling infrastructure of the Integrated Forecasting System: recent advances and future challenges. Technical Report 760, ECMWF, Reading, UK.

Wedi, N.P., Hamrud, M. and Mozdzynski, G. (2013) A fast spherical harmonics transform for global NWP and climate models. Monthly Weather Review, 141(10), 3450–3461. https://doi.org/10.1175/MWR-D-13-00016.1

Wood, N., Staniforth, A., White, A.A., Allen, T., Diamantakis, M., Gross, M., Melvin, T., Smith, C., Vosper, S.B., Zerroukat, M. and Thuburn, J. (2014) An inherently mass-conserving semi-implicit semi-Lagrangian discretization of the deep-atmosphere global non-hydrostatic equations. Quarterly Journal of the Royal Meteorological Society, 140, 1505–1520. https://doi.org/10.1002/qj.2235

Zheng, Y. and Marguinaud, P. (2018) Simulation of the performance and scalability of message passing interface (MPI) communications of atmospheric models running on exascale supercomputers. Geoscientific Model Development, 11(8), 3409–3426. https://doi.org/10.5194/gmd-11-3409-2018

**How to cite this article:** Degrauwe D, Voitus F, Termonia P. A non-spectral Helmholtz solver for numerical weather prediction models with a mass-based vertical coordinate. *Q J R Meteorol. Soc.* 2021;147:30–44. https://doi.org/10.1002/qj.3902

## APPENDIX

Reduction of the linear system to a single Helmholtz equation Substitution of Equations (8), (9) and (10) in Equations (6) and (7) yields

$$[\mathbf{I} - \delta t^2 c_*^2 \nabla^2 (\mathbf{I} + C_{va}/C_{pa}(\mathbf{G}^* \mathbf{S}^* - \mathbf{G}^* - \mathbf{S}^* + \mathbf{N}^*))]D - \delta t^2 c_*^2 \nabla^2 (\mathbf{I} - C_{va}/C_{pa}\mathbf{G}^*)d = D^*, \quad (A1)$$
$$ \left[ \mathbf{I} - \delta t^2 \frac{c_*^2}{rH_*^2} \mathbf{L}^* \right] d - \delta t^2 \frac{c_*^2}{rH_*^2} \mathbf{L}^* (\mathbf{I} - C_{va}/C_{pa} \mathbf{S}^*) D = d^*, \quad (A2) $$

where $\mathbf{I}$ denotes the identity vertical matrix, $c_*^2 = (C_{pa}/C_{va})R_a T^*$ and $H_* = R_a T^*/g$. The right-hand sides of (A1) and (A2) are given respectively by

$$ D^* = \tilde{D} - \delta t R T^* \nabla^2 [(\mathbf{I} - \mathbf{G}^*) \tilde{\hat{q}} + \mathbf{G}^* (\tilde{T}/T^*) + (\tilde{\pi}_s/\pi_s^*)], \quad (A3) $$

$$ d^* = \tilde{d} - \delta t \frac{R T^*}{r H_*^2} \mathbf{L}^* \tilde{\hat{q}}. \quad (A4) $$

For convenience, the vertical Helmholtz operator is introduced as

$$ \mathbf{H}^* = \mathbf{I} - \delta t^2 \frac{c_*^2}{rH_*^2} \mathbf{L}^*. $$

Since this latter operator is a constant-coefficient positive-definite matrix (by construction), its inverse

matrix (denoted $\mathbf{H}^{*-1}$) exists, and only needs to be computed once.

The vertical divergence $d$ is then related to the horizontal divergence $D$ through

$$ d = \mathbf{H}^{*-1} \left( d^* + \delta t^2 \frac{c_*^2}{rH_*^2} \mathbf{L}^* (\mathbf{I} - C_{va}/C_{pa} \mathbf{S}^*) D \right). \quad (A5) $$

Substituting in (A1) finally reduces the system to a single Helmholtz equation

$$ (\mathbf{I} - \delta t^2 \nabla^2 \mathbf{B}_D^*) D = D^{**}, \quad (A6) $$

where the matrix $\mathbf{B}_D^*$ is given by

$$ \mathbf{B}_D^* = \frac{R_a^2 T^*}{C_{pa}} \left( \mathbf{G}^* \mathbf{S}^* + \frac{C_{pa}}{R_a} \mathbf{N}^* \right) + c_*^2 (\mathbf{I} - C_{va}/C_{pa} \mathbf{G}^*) \mathbf{H}^{*-1} (\mathbf{I} - C_{va}/C_{pa} \mathbf{S}^*) \quad (A7) $$

and the right-hand side is given by

$$ D^{**} = D^{*} + \delta t^2 c_*^2 \nabla^2 (\mathbf{I} - C_{va}/C_{pa} \mathbf{G}^*) \mathbf{H}^{*-1} d^*. \quad (A8) $$