# Application of the spectral-element method to the axisymmetric Navier-Stokes equation

A. Fournier, H. -P. Bunge, R. Hollerbach, J. -P. Vilotte

► To cite this version:

A. Fournier, H. -P. Bunge, R. Hollerbach, J. -P. Vilotte. Application of the spectral-element method to the axisymmetric Navier-Stokes equation. Geophysical Journal International, 2004, 156, pp.682-700. 10.1111/j.1365-246X.2004.02149.x. insu-03600288

HAL Id: insu-03600288

https://insu.hal.science/insu-03600288

Submitted on 7 Mar 2022

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
# Application of the spectral-element method to the axisymmetric Navier-Stokes equation

A. Fournier,<fnref n="1" />,<fnref n="3" /> H.-P. Bunge,<fnref n="1" /> R. Hollerbach<fnref n="1" />,<fnref n="2" /> and J.-P. Vilotte<fnref n="3" />

Accepted 2003 August 15. Received 2003 July 24; in original form 2002 August 19

## SUMMARY

We present an application of the spectral-element method to model axisymmetric flows in rapidly rotating domains. The primitive equations are discretized in space with local tensorized bases of high-order polynomials and in time with a second-order accurate scheme that treats viscous and rotational effects implicitly. We handle the pole problem using a weighted quadrature in elements adjacent to the axis of rotation. The resulting algebraic systems are solved efficiently using preconditioned iterative procedures. We validate our implementation through comparisons with analytic and purely spectral solutions to laminar flows in a spherical shell. This axisymmetric tool is the kernel on which complexity will be added subsequently in the long-term prospect of building a parallel spectral-element based geodynamo model.

**Key words:** geodynamo, rotating flows, spectral-element method.

## 1 INTRODUCTION

As the Earth sheds its heat, its interior undergoes large-scale convective motions. Inside its liquid metallic outer core, these motions generate in turn the geomagnetic field, as was originally proposed by Larmor (1919). More than 80 yr after his founding hypothesis it is now widely accepted that thermochemical convection indeed provides enough energy to power the geodynamo (Gubbins & Roberts 1987). Modelling this complex magnetohydrodynamic process is made difficult by the low molecular viscosity of iron under core conditions (Poirier 1988; de Wijs et al. 1998). In fact, the ratio of viscous stresses to the Coriolis force in the force balance of the core, measured by the Ekman number *E*, is very small ($10^{-12}$ at most) resulting in sharp viscous boundary layers (called *Ekman layers*) of a few metres. Thus, we have little hope in the near future of resolving these small length-scales numerically in a computer model of the geodynamo, even if we account for the impressive rise in (parallel) computing power expected in the coming years.

Despite these difficulties great insight into the working of the geodynamo has been gained over the past decade thanks to progress made jointly by laboratory and numerical modellers (Busse 2000). As a matter of fact, Glatzmaier & Roberts (1995) simulated the magnetohydrodynamics of an artificially hyperviscous core and presented the first computer simulation of a geomagnetic field reversal using a 3-D spherical dynamo model. Although far from the appropriate parameter regime, their model produced a magnetic field remarkably similar to the magnetic field of the Earth. This final result led subsequently these and other authors (Glatzmaier & Roberts 1996a; Kuang & Bloxham 1997) to investigate a range of geophysical problems related to the dynamics of the Earth's core,

including the differential rotation of the inner core (Glatzmaier & Roberts 1996b), the angular momentum budget of the Earth (Bloxham 1998), the secular variation of the Earth's magnetic field (Bloxham 2000a) and, in a palaeomagnetic perspective, the validity of the geocentric axial dipole hypothesis (Bloxham 2000b).

From a numerical standpoint, current dynamo models are based on spherical harmonics to describe the horizontal dependency of the variables (Glatzmaier 1984; Kuang & Bloxham 1999; Hollerbach 2000). The method is certainly the most natural one to consider when attacking the problem of modelling the circulation of a convecting (or precessing) Boussinesq liquid metal in spherical geometry (see also Tilgner 1999). For instance, the analytic character of spherical harmonics permits one to perform a poloidal–toroidal decomposition both of the magnetic and the velocity fields, thus satisfying exactly the solenoidal requirements upon these vector fields (Glatzmaier 1984). Moreover, their use leads to a weak numerical dispersion, and they achieve an almost uniform resolution of the spherical surface. They also circumvent the pole problem that arises when using spherical ($r, \theta, \phi$) coordinates. Unfortunately, the main drawback of spherical harmonics originates from their global definition, which requires a rather expensive pseudospectral calculation of the non-linear terms, and consequently gives rise to a difficult processing on parallel computers. As a result, current dynamo simulations are not performed at Ekman numbers smaller than $10^{-4}$ (Christensen et al. 1999) for simulations that span several magnetic diffusion timescales, unless one uses a controversial hyperviscosity (Zhang & Jones 1997; Grote et al. 2000).

Questions remain on the ability of these smooth models to reflect turbulent motions in the Earth's core and to reproduce long-term features of the geomagnetic field, as pointed out by Dormy et al.

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

<footnote n="1">¹Princeton University Department of Geosciences, Princeton, NJ, 08540, USA. E-mail: amf2@alumni.princeton.edu</footnote>

<footnote n="2">²University of Glasgow, Department of Mathematics, Glasgow G12 8QW, UK</footnote>

<footnote n="3">³Département de Modélisation Physique et Numérique, IPGP, 75252 Paris Cedex 05, France</footnote>
(2000). There is hope, however, that if one is able to prescribe a small enough Ekman number, one will reach a parameter regime asymptotically appropriate for the Earth's core. Indeed, from a theoretical standpoint the core has two options as to how to operate its dynamo, commonly referred to as the weak and strong field regimes (Roberts 1978). The dynamo inside the Earth may fluctuate between these states (Zhang & Gubbins 2000), but looking at computer models of the dynamo we have yet to discover how large rotation has to be before a dynamo has the choice between these two distinct regimes. St. Pierre (1993) found that $E = O(10^{-5})$ was sufficiently small to obtain a subcritical strong field dynamo in his plane layer study. However, before applying these results to the real Earth one would have to repeat them in spherical geometry, and vary the Ekman number (and other relevant parameters) enough to be able to determine whether or not there are these two distinct regimes. Indeed, that is precisely the ultimate objective of this work. Nevertheless, St. Pierre's results suggest that the $O(10^{-4})$ Ekman number currently being used may need to be reduced by an order of magnitude before one is even qualitatively in the right regime.

A reduction in Ekman number could be attained by using numerical methods that execute efficiently on modern parallel computers via domain decomposition and explicit message-passing. In fact, domain decomposition methods based on explicit message-passing have already proven to be successful in finite-element models simulating flow inside the Earth's mantle at high convective vigour (Bunge & Baumgardner 1995). Moreover, these methods are well suited to the growing trend of using cost-effective, off-the-shelf PC-clusters in geophysical modelling (Bunge & Dalton 2001). Consequently, our long-term effort aims at developing a numerical dynamo model that retains the accuracy and robustness of spectral methods while performing well on modern parallel computers such as clusters of PCs. Our approach is based upon the use of the spectral-element method (SEM), a variational technique that relies on high-order local shape functions (Patera 1984; Bernardi & Maday 1992). The SEM, in fact, combines the geometrical flexibility of the finite-element method with the exponential convergence and weak numerical dispersion of spectral methods (Maday & Patera 1989). In addition, its local character lends itself naturally to domain decompositions, and allows for non-uniform resolution inside the computational domain, i.e. for grid-refinements in localized regions such as the narrow Ekman boundary layers inside the core. Recent geophysical applications of the SEM include ocean-atmosphere modelling (Taylor et al. 1997; Levin et al. 2000; Giraldo 2001) as well as regional and global seismic wave propagation (Komatitsch & Vilotte 1998; Komatitsch & Tromp 1999; Capdeville et al. 2002; Chaljub et al. 2003). To our knowledge, however, the SEM has not yet been applied to models of deep Earth flows, neither in the mantle nor in the core.

While Chan et al. (2001) already investigated the implementation of a finite-element method to solve the spherical kinematic dynamo problem, we present and validate here the application of the SEM to the Navier-Stokes equation in an axisymmetric, non-magnetic context. This axisymmetric case can readily be generalized to fully 3-D applications by coupling the SEM in the meridional plane with a Fourier expansion in the longitudinal direction. In this so-called Fourier-spectral-element approach (Bernardi et al. 1999), the 3-D problem is broken into a collection of meridional subproblems, which in turn may be parallelized into a number of spatial subdomains. We use cylindrical $(s, \phi, z)$ coordinates and solve for primitive variables. We thus do not rely on the expansion of the velocity in terms of a poloidal and a toroidal field: a poloidal-toroidal decomposition generates high-order differential operators which can in turn lead to a substantial numerical dispersion. We therefore show

explicitly in this paper how the divergence-free requirement on the velocity field is satisfied with our method. We show furthermore, how we handle the singularities at the axis of rotation by using a weighted Gauss-Lobatto quadrature (Bernardi et al. 1999).

The outline of this paper is as follows: Section 2 recalls the system of equations of interest, and its detailed variational treatment is presented in Section 3. We then describe the spatial and temporal discretizations of the variational problem in Sections 4 and 5. The validation of the implementation proceeds by comparing SEM results with analytical solutions for steady and unsteady Stokes problems (Section 6) and with published spectral solutions in a rapidly rotating context (Section 7). The SEM is shown in all cases to exhibit the spectral convergence properties of standard spectral methods and to provide numerical accuracy of better than one per mil relative to the reference solution. A concluding discussion follows in Section 8.

## 2 GOVERNING EQUATIONS

As illustrated in Fig. 1, we are interested in describing the axisymmetric motion of an incompressible Newtonian fluid filling an axisymmetric container of arbitrary meridional shape $\Omega$. The revolution of $\Omega$ around the axis of symmetry $\Gamma$ gives rise to the full 3-D domain $\tilde{\Omega}$. We assume that the rotation rate $\omega$ is constant and that the rotation vector $\omega$ is parallel to $\Gamma$. The unit vector along this axis is denoted by $\hat{\omega}$. Under these conditions, the flow of the fluid is governed by the following non-dimensional equations (e.g. Gubbins & Roberts 1987):

$$ \partial_t \mathbf{u} + 2\hat{\mathbf{z}} \times \mathbf{u} = -\nabla p + E \Delta \mathbf{u} + \mathbf{f} \text{ in } \Omega, \quad (1a) $$

$$ \nabla \cdot \mathbf{u} = 0 \text{ in } \Omega, \quad (1b) $$

<figure><img src="image_3.png" /><figcaption>**Figure 1.** The approach we describe aims primarily at solving the Navier-Stokes equation in spherical/spheroidal shells (right). Its flexibility allows, however, to handle axisymmetric containers of more complicated shape (left) that one could use in a laboratory experiment. In each case, the 3-D domain $\tilde{\Omega}$ follows from the revolution of its meridional section $\Omega$ around its axis of symmetry $\Gamma$. $\partial\Omega$ is the boundary of $\Omega$.</figcaption></figure>
where **u** is the velocity of the fluid, *p* is its pressure augmented of the centrifugal acceleration, and **f** denotes the body forces which include potentially the non-linear interactions. The actual treatment of the non-linearities is beyond the scope of the present paper. However, let us mention that they may be dealt with in an explicit fashion, by absorbing them into **f**. The relative importance of viscous to rotational effects is measured by the non-dimensional Ekman number:

$$E = \frac{\nu}{\omega L^2}, \qquad (2)$$

in which $\nu$ represents the kinematic viscosity of the fluid and $L$ the depth of the container. For problem (1) to be well-posed, we specify boundary conditions $\mathbf{u}_b(t)$ on the domain boundary $\partial\Omega$ (which does not include the intersection of $\Omega$ with $\Gamma$), and conditions on the initial state $\mathbf{u}_0(\mathbf{x})$.

### 3 VARIATIONAL FORMULATION

The spectral-element method, such as the standard finite-element method, relies on the variational formulation of the equations of interest. At any time $t \in [0, T]$, we consider the velocity and pressure field that we denote by $\mathbf{u}_t(\mathbf{x}) = \mathbf{u}(\mathbf{x}, t)$ and $p_t(\mathbf{x}) = p(\mathbf{x}, t)$. Using cylindrical coordinates $(s, \phi, z)$, the three vector components of $\mathbf{u}_t$ will subsequently be indicated by $(u_{t,s}, u_{t,\phi}, u_{t,z})$. The variational formulation of problem (1) is obtained by multiplying eqs (1a) and (1b) with appropriate trial functions and integrating the resulting system over the domain $\Omega$. An elementary volume of integration $d\Omega$ is a torus, obtained by the revolution of a rectangular meridional section of area $dsdz$ around $\Gamma$ (see Fig. 1). It is thus given by

$$d\Omega = 2\pi s \, ds \, dz. \qquad (3)$$

Following Bernardi et al. (1999), we define the space of square integrable functions $L_1^2(\Omega)$

$$L_1^2(\Omega) = \left\{ w : \Omega \to \mathbb{R}, \|w\| = \left(\int_{\Omega} w^2 d\Omega\right)^{1/2} < \infty \right\}. \quad (4)$$

We associate the inner product $(\cdot, \cdot)_1$

$$\forall (f, g) \in L_1^2(\Omega) \times L_1^2(\Omega), \quad (f, g)_1 = \int_{\Omega} fg \, d\Omega. \qquad (5)$$

We also introduce the 2-D weighted Sobolev space $H_1^1(\Omega)$ as the subspace of $L_1^2(\Omega)$ containing those functions whose first partial derivatives are also square integrable

$$H_1^1(\Omega) = \{w \in L_1^2(\Omega), \partial_s w \in L_1^2(\Omega), \partial_z w \in L_1^2(\Omega)\}. \qquad (6)$$

To account for boundary conditions, it is necessary to define the subspace of functions in $H_1^1(\Omega)$ which vanish on $\partial\Omega$

$$H_{1\circ}^1(\Omega) = \{w \in H_1^1(\Omega), w = 0 \text{ on } \partial\Omega\}. \qquad (7)$$

In the axisymmetric case considered here, the three components of the velocity $(u_{t,s}, u_{t,\phi}, u_{t,z})$ have to satisfy different conditions on the axis $\Gamma$. Indeed, $u_{t,s}$ and $u_{t,\phi}$ must vanish on $\Gamma$ whereas $u_{t,z}$ must satisfy the symmetry condition $\partial_s u_{t,z} = 0$.

$$u_{t,s} = u_{t,\phi} = 0, \text{ on } \Gamma, \qquad (8a)$$

$$\partial_s u_{t,z} = 0, \text{ on } \Gamma. \qquad (8b)$$

The latter condition is a so-called natural condition, and is automatically satisfied by the solution of the associated variational problem.

However, the nullity condition on $u_{t,s}$ and $u_{t,\phi}$, which is of the essential kind, has to be enforced and requires the introduction of $V_1^1(\Omega)$:

$$V_1^1(\Omega) = \{w \in H_1^1(\Omega), w = 0 \text{ on } \Gamma\}. \qquad (9)$$

Again, the imposition of the boundary conditions on $\partial\Omega$ requires us to define $V_{1\circ}^1(\Omega)$ as

$$V_{1\circ}^1(\Omega) = \{w \in V_1^1(\Omega), w = 0 \text{ on } \partial\Omega\}. \qquad (10)$$

We can now define the space of admissible velocities at any given time $t$

$$\mathbf{H}^1(\Omega) = V_1^1(\Omega) \times V_{1\circ}^1(\Omega) \times H_{1\circ}^1(\Omega), \qquad (11)$$

and the space of velocity trial functions

$$\mathbf{H}_{\circ}^{1}(\Omega) = V_{1\circ}^{1}(\Omega) \times V_{1\circ}^{1}(\Omega) \times H_{1\circ}^{1}(\Omega). \qquad (12)$$

These two spaces, therefore, differ only in that the trial functions have to vanish where the value of the velocity is imposed. As no boundary condition is prescribed on the pressure field, the space of pressure trial functions is the same as the space of pressure basis functions, and consists simply of the space of square integrable functions defined over $\Omega$. This collection of spaces now enables us to recast problem (1) in its equivalent variational form:

For any time $t$ in $[0, T]$ find $(\mathbf{u}_t, p_t)$ in $\mathbf{H}^1(\Omega) \times L_1^2(\Omega)$ with $\mathbf{u}_t - \mathbf{u}_b(t)$ in $\mathbf{H}_{\circ}^1(\Omega)$, such that:

$$\forall \mathbf{v} \in \mathbf{H}_{\circ}^{1}(\Omega), \quad (\partial_t \mathbf{u}_t, \mathbf{v})_1 + (2\hat{\mathbf{z}} \times \mathbf{u}_t, \mathbf{v})_1 + \mathbf{Ea}(\mathbf{u}_t, \mathbf{v}) - d(\mathbf{v}, \mathbf{p}_t) = (\mathbf{f}, \mathbf{v})_1, \qquad (13a)$$

$$\forall q \in L_1^2(\Omega), \quad d(\mathbf{u}_t, q) = 0. \qquad (13b)$$

This problem is a standard saddle-point problem, where eq. (13a) has to be solved for a velocity that satisfies the divergence-free constraint (13b). The pressure field is the Lagrange multiplier associated with this constraint. Here we have introduced the bilinear form $a$, which is the variational equivalent of the Laplacian:

$$a(\mathbf{u}, \mathbf{v}) = a_0(u_s, v_s) + a_0(u_\phi, v_\phi) + a_0(u_z, v_z) + \int_{\Omega} \frac{1}{s^2}(u_s v_s + u_\phi v_\phi) d\Omega, \qquad (14)$$

in which

$$a_0(f, g) = \int_{\Omega} (\partial_s f \partial_s g + \partial_z f \partial_z g) d\Omega. \qquad (15)$$

The divergence/gradient form $d$ is given by

$$d(\mathbf{v}, \mathbf{q}) = \int_{\Omega} \mathbf{q} \left( \partial_s \mathbf{v}_s + \frac{\mathbf{v}_s}{\mathbf{s}} + \partial_z \mathbf{v}_z \right) d\Omega. \qquad (16)$$

Note that both of these forms appear in the variational momentum eq. (13a) after an integration by parts, and that the pressure does not have to be continuous on $\Omega$. Importantly, it can be shown (Bernardi & Maday 1992) that the existence of a unique solution to the saddle-point problem (13) is guaranteed, if a so-called compatibility condition between the velocity and pressure spaces is respected. We will return to this point in more detail in the following section.

### 4 SPECTRAL-ELEMENT METHODOLOGY

In this section, we describe how the weak formulation (13) of the original problem (1) is discretized in space using the spectral-element method. We want to restrict $\mathbf{u}_t$ and $p_t$ in (13) to finite-dimensional spaces $X_h$ and $Y_h$, respectively and denote their discretized version by $\mathbf{u}_{t,h}$ and $p_{t,h}$. As illustrated in the top row of
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Top, the domain $\Omega$ is broken into a collection of $n_e$ non-overlapping spectral elements $\Omega^e$. Each $\Omega^e$ is the image of the reference square $\Lambda^2 = [-1, 1]^2$ under an invertible local mapping $\mathcal{F}^e$. For spatially complex $\Omega$ (a), $\mathcal{F}^e$ is a subparametric transformation, otherwise $\mathcal{F}^e$ is analytical (b). Middle, when $\Omega^e$ is not adjacent to the axis, its local shape functions for velocity and pressure are defined by the tensor product of the Lagrangian interpolants (LI) defined over the family of Gauss-Lobatto Legendre points of order $N$. (c) The 11 velocity LI defined by the GLL points of order 10. (d) The 9 Pressure LI defined by the interior GLL points of order 10. Bottom, when $\Omega^e$ is adjacent to the axis, the discretized velocity must exhibit the proper behaviour when approaching the axis. A weighted quadrature is thus used which incorporates the cylindrical radius in its weight, and has one velocity point strictly on the axis. (e) The 11 velocity LI defined by the weighted GLL points of order 10. Notice the resulting asymmetry in the shape functions in contrast to (c). (f) The 9 pressure LI defined by the interior weighted GLL points of order 10.</figcaption></figure>

Fig. 2, we define these spaces by decomposing the global domain $\Omega$ into a collection of $n_e$ non-overlapping elements $\Omega^e$, such that

$$ \bar{\Omega} = \bigcup_{e=1}^{n_e} \bar{\Omega}^e. \qquad (17) $$

Here, each $\Omega^e$ is the image of a reference square $\Lambda^2 = [-1, +1]^2$ under a local invertible mapping $\mathcal{F}^e: (\xi, \eta) \in \Lambda^2 \Rightarrow (s, z) \in \Omega^e$ with a well-defined inverse. Dealing with a deformed quadrangle enables us to perform a separation of variables $(\xi, \eta)$ and therefore to use a tensorized basis. Figs 2(a) and (b) illustrate the two options we have in practice to implement this mapping. When the shape of the domain $\Omega$ is complex (see Fig. 2a), we use a so-called subpara-

metric mapping (Hughes 1987; Reddy 1993), where the transformation is parametrized by the datum of the images of control points in $\Lambda^2$. When $\Omega$ is simple (e.g. when it is the meridional section of a spheroid, see Fig. 2b), an analytic expression for $\mathcal{F}^e$ is preferred.

In each element, velocity and pressure are approximated locally via a tensorized basis of high-order polynomials (shown in the middle and the bottom panels of Fig. 2), hence the terminology of spectral elements introduced by Patera (1984). To avoid spurious pressure modes, Bernardi & Maday (1988) suggest taking

$$ X_h = \mathbf{H}^1(\Omega) \cap \mathbf{P}_{N,n_e}, \qquad (18) $$

$$ Y_h = L_1^2(\Omega) \cap \mathbb{P}_{N-2,n_e}, \qquad (19) $$
where

$$
\mathbb{P}_{N,n_e} = \{w[\mathcal{F}^e(\xi, \eta)] \mid \Omega^e \in P_N(\xi) \times P_N(\eta), e = 1, n_e\} \quad (20)
$$

and

$$
\mathbf{P}_{N,n_e} = \mathbb{P}_{N,n_e} \times \mathbb{P}_{N,n_e} \times \mathbb{P}_{N,n_e}. \quad (21)
$$

Here, $P_N$ is the space of those polynomials defined over $[-1, 1]$ of degree less or equal to $N$. In other words, each component of the restriction of the velocity in a given element $\Omega^e$ is described in terms of the tensor product of polynomials of order $N$ along the $\xi$ and $\eta$ directions. Definition (18) also requires the velocity to be continuous at the boundary between two elements. For the pressure field, instead, the order of the polynomials is set to $N-2$, and definition (19) does not require pressure to be continuous at the elements boundaries. It can been shown, as in Bernardi & Maday (1992), that the lower degree used to discretize pressure in this so-called $P_N - P_{N-2}$ approach provides a unique discrete solution $(\mathbf{u}_{t,h}, p_{t,h})$ to the problem of interest. Similarly, the discrete space of velocity trial functions is defined as

$$
\mathbf{X}_{\circ,h} = \mathbf{H}_{\circ}^{1}(\Omega) \cap \mathbf{P}_{N,n_e}. \qquad (22)
$$

We now describe in detail the exact nature of the polynomials. The
basis for $P_N$ is related to the Gauss-type quadrature formula used
to evaluate the integrals which appear in the variational formulation
(13). Such integrals can be broken into a sum of elemental integrals,
i.e. one can write

$$
\int_{\Omega} f \, d\Omega = \sum_{e=1}^{n_e} \int_{\Omega^e} f \, d\Omega^e. \qquad (23)
$$

As we use cylindrical coordinates, elements adjacent to the axis of symmetry $\Gamma$ (which we will hereafter refer to as axial elements) have to be distinguished from elements away from the axis. We group the $n_\Gamma$ axial elements in $\Omega_\Gamma$ such that

$$
\overline{\Omega}_{\Gamma} = \bigcup_{e=1}^{n_{\Gamma}} \overline{\Omega}^{e}, \qquad (24)
$$

whereas the non-axial elements are grouped into

$$
\Omega_{\emptyset}, \overline{\Omega}_{\emptyset} = \overline{\Omega} \setminus \overline{\Omega}_{\Gamma} = \bigcup_{e=n_{\Gamma}+1}^{n_e} \overline{\Omega}^e. \quad (25)
$$

When $\Omega^e$ is not axial (such as elements $\Omega^3$, $\Omega^4$, $\Omega^5$, $\Omega^6$ in Fig. 2a, or elements $\Omega^3$, $\Omega^4$ in Fig. 2b), $P_N$ is spanned by the set of Lagrangian interpolants $h_i^N$, $i \in \{0, \dots, N\}$ defined by the $N+1$ Gauss–Lobatto Legendre (GLL) points $\xi_i^N$, $i \in \{0, \dots, N\}$ on $[-1, 1]$. Fig. 2(c) displays this family of polynomials for $N=10$. For the pressure, the basis for $P_{N-2}$ is the set of Lagrangian interpolants $h_i^{N-2}$ defined on the interior GLL nodes $\xi_i^N$, $i \in \{1, \dots, N-1\}$ (see Fig. 2d). We are now in a position to define the quadrature rule over the non-axial elements

$$
\forall \Omega^e \in \Omega_\emptyset, \int_{\Omega^e} f \, d\Omega^e \approx \sum_{i,j=0}^{N} \rho_i \rho_j s f[\mathcal{F}^e(\xi_{ij})] |\mathcal{J}^e|(\xi_{ij}), \quad (26)
$$

where the $\rho_i$, $i \in \{0, \dots, N\}$ are the quadrature weights associated
with the Gauss–Lobatto formula of order $N$, $\xi_{ij} = (\xi_i^N, \xi_j^N)$, and
$|\mathcal{J}^e|$ stands for the Jacobian of the mapping $\mathcal{F}^e$.

When $\Omega^e$ is in contact with the axis of symmetry (elements $\Omega^1$ and $\Omega^2$ in Figs 2a and b), a different quadrature must be used to perform the integration in the $\xi$-direction. Indeed, the presence of an $s$ factor in the elementary volume $d\Omega^e$ would lead to an undetermined system of the form ‘0 = 0’, if integrals were to be evaluated on collocation points located on $\Gamma$ (Gerritsma & Phillips 2000). This

and the enforcement of the essential boundary conditions (8a) on
$\Gamma$ favours the use of a weighted Gauss–Lobatto quadrature, which
incorporates the cylindrical radius in its weight. We denote by $\zeta_i^N$,
$i \in \{0, \dots, N\}$ and $\sigma_i$, $i \in \{0, \dots, N\}$, respectively, the nodes and
weights associated with this new quadrature. For any polynomial $\Phi$
in $P_{2N-1}(\Lambda)$, we then have

$$
\int_{\Lambda} \Phi(\xi)(1+\xi) d\xi = \sum_{i=0}^{N} \sigma_i \Phi(\zeta_i^N). \quad (27)
$$

In the $\xi$-direction, a basis for $P_N$ is thus the set of Lagrangian
interpolants $l_i^N$, $i \in \{0, \dots, N\}$ defined by the $\zeta_i^N$, $i \in \{0, \dots, N\}$,
and a basis for $P_{N-2}$ is the set of Lagrangian interpolants $l_i^{N-2}$
defined by the $\zeta_i^N$, $i \in \{1, \dots, N-1\}$. We show these two bases in
Figs 2(e) and (f), respectively, for a polynomial order $N=10$. In the
$\eta$ direction, for which no underdetermination is expected, we retain
the quadrature rule and the related basis that we defined previously
for non-axial elements. In summary, the following integration rule
applies for elements adjacent to $\Gamma$

$$
\forall \Omega^e \in \Omega_\Gamma, \int_{\Omega^e} f \, d\Omega^e \approx \sum_{i,j=0}^{N} \sigma_i \rho_j \frac{s(\zeta_{ij}^N)}{1 + \zeta_i^N} f[\mathcal{F}^e(\zeta_{ij}^N)] |\mathcal{J}^e|(\zeta_{ij}^N) \quad (28)
$$

where $\zeta_{ij} = (\zeta_i^N, \zeta_j^N)$. The apparent singularity involving the term $s(\zeta_i^N, \zeta_j^N)/(1 + \zeta_i^N)$ when $\zeta_i^N = -1$, or, equivalently, when $s = 0$ is easily removed by the application of L'Hospital's rule. Further details concerning the quadrature formulae can be found in Appendix A, or to a greater extent in (Bernardi et al. 1999 chapters IV and VI). Note that in any situation, the basis for the velocity is continuous across subdomain boundaries, while the basis for the pressure is not. The discrete velocity field is given by

$$
\mathbf{u}_{t,h}[\mathbf{x}(\xi, \eta)] = \sum_{e=1}^{n_{\Gamma}} \sum_{i,j=0}^{N} \left( u_{t,s}^{\text{eij}} u_{t,\phi}^{\text{eij}} u_{t,z}^{\text{eij}} \right) l_i^N(\xi) h_j^N(\eta) \\
+ \sum_{e=n_{\Gamma}+1}^{n_e} \sum_{i,j=0}^{N} \left( u_{t,s}^{\text{eij}} u_{t,\phi}^{\text{eij}} u_{t,z}^{\text{eij}} \right) h_i^N(\xi) h_j^N(\eta). \quad (29)
$$

Here, the $(u_{t,s}^{\text{eij}}, u_{t,\phi}^{\text{eij}}, u_{t,z}^{\text{eij}})$ are the nodal velocities at the collocation points in the $\text{eth}$ element and $\mathbf{x} = (s, z)$ is the meridional position vector. Likewise, the discrete pressure reads:

$$
p_{t,h}[\mathbf{x}(\xi, \eta)] = \sum_{e=1}^{n_{\Gamma}} \sum_{i,j=1}^{N-1} p_i^{\text{eij}} l_i^{N-2}(\xi) h_j^{N-2}(\eta) \\
+ \sum_{e=n_{\Gamma}+1}^{n_e} \sum_{i,j=1}^{N-1} p_i^{\text{eij}} h_i^{N-2}(\xi) h_j^{N-2}(\eta). \quad (30)
$$

In the remainder of this paper, $\mathbf{u}_t = (\mathbf{u}_{t,s}, \mathbf{u}_{t,\phi}, \mathbf{u}_{t,z})$ will be the vector of velocity unknowns, and $\mathbf{p}_t$ the vector of pressure unknowns. Fig. 3 displays a simple spherical mesh showing the collocation points associated with velocity and pressure.

The spatial discretization of problem (13) proceeds by specifying
the trial functions. We follow a classical Galerkin approach, and
build $X_{\circ,h}$ and $Y_h$ with the nodal shape functions associated with
the velocity and pressure degrees of freedom, respectively. This
leads to the semi-discrete version of problem (13):

Find at any time $t \in [0, T]$ the solution $(\mathbf{u}_t, \mathbf{p}_t)$ of

$$
M \partial_t \mathbf{u}_t + C \mathbf{u}_t + E K \mathbf{u}_t - \mathbf{D}^T \mathbf{p}_t = M \mathbf{f}_t, \quad (31a)
$$

$$
-D \mathbf{u}_t = 0. \qquad (31b)
$$

In this system, **M** is the diagonal mass matrix, **C** is the Coriolis antisymmetric matrix, **K** is the stiffness matrix, and **D**/ **D**^T denotes
<figure><img src="image_5.png" /><figcaption>**Figure 3.** SEM mesh for spherical shell geometry with $n_e = 4$ spectral elements of order $N = 4$ used for analytic and numerical benchmarks (see the text). Velocity and pressure nodes are represented by black triangles and white circles, respectively. Axial elements $\Omega^1$ and $\Omega^2$ resort to a weighted Gauss-Lobatto quadrature (see the text) resulting in a different spacing of nodes in latitude relative to non-axial elements $\Omega^3$ and $\Omega^4$, the nodes of which are the images of the standard GLL points.</figcaption></figure>

the divergence/gradient matrix. On the right-hand side, $\mathbf{f}_t$ denotes the forcing vector. An extensive derivation of system (31), together with a detailed description of the various matrices is given in Appendix B. It is worthwhile to mention that these matrices are not stored, except for the diagonal mass matrix. Instead, because of the tensorized formulation, the result of their actions on nodal vectors is directly computed and assembled. Storing the stiffness matrix and applying it to a nodal field would require $O(n_e N^4)$ operations. Instead, the resulting field can be computed in $O(n_e N^3)$ operations, along with a significant reduction in memory requirements.

## 5 TIME DISCRETIZATION

Having presented the spatial discretization of (13), we are now ready to specify how time marches on. We break the interval $[0, T]$ into segments of equal length $\Delta t$, and denote by $\mathbf{u}_n$ and $\mathbf{p}_n$ the value of $\mathbf{u}_t$ and $\mathbf{p}_t$ at $t = t_n = n \Delta t$. The time derivative in (31a) is approximated via a second-order backward differentiation formula (BDF2):

$$ \partial_t \mathbf{u}_n = \frac{3\mathbf{u}_n - 4\mathbf{u}_{n-1} + \mathbf{u}_{n-2}}{2\Delta t}. \quad (32) $$

At each time step $t_n$, we have to solve a modified Stokes problem of the form:

$$ \mathbf{H}\mathbf{u}_n - \mathbf{D}^T \mathbf{p}_n = \mathbf{M}\mathbf{t}_n, \quad (33a) $$

$$ -\mathbf{D}\mathbf{u}_n = 0, \quad (33b) $$

in which:

$$ \mathbf{H} = \frac{3}{2\Delta t}\mathbf{M} + \mathbf{C} + \mathbf{E}\mathbf{K} \quad (34) $$

is a Helmholtz operator modified by the addition of the effects due to rotation, and $\mathbf{t}_n = \mathbf{f}_n + (4\mathbf{u}_{n-1} - \mathbf{u}_{n-2})/(2\Delta t)$.

Our strategy to invert the coupled system (33) follows a so-called Operator Integrated Factor (OIF) splitting scheme, originally introduced by Maday *et al.* (1990). This is a modified version of the more standard Uzawa algorithm (Arrow *et al.* 1958), which we describe briefly here. Problem (33) is the discrete version of the original saddle-point problem (13). In order to apply an Uzawa method, one would split (33) and solve first for the pressure field $\mathbf{p}_n$. Indeed, multiplying (33a) by $\mathbf{D}\mathbf{H}^{-1}$ and using the discrete incompressibility condition (33b) leads to the following elliptic system:

$$ \mathbf{D}\mathbf{H}^{-1}\mathbf{D}^T \mathbf{p}_n = -\mathbf{D}\mathbf{H}^{-1}\mathbf{M}\mathbf{t}_n. \quad (35) $$

Once $\mathbf{p}_n$ is known, it can be used in (33a) to compute the velocity field $\mathbf{u}_n$. Note, however, that the size of problem (35) precludes a direct solution, and that each iteration would require that one inverts **H** (iteratively as well), resulting in an expensive procedure. The scheme proposed by Maday *et al.* (1990) overcomes this problem by relying on the fact that the SEM mass matrix is diagonal, and therefore straightforward to invert. Following Couzy (1995), we write (33) in the equivalent matrix form

$$ \begin{bmatrix} \mathbf{H} & -\mathbf{D}^T \\ -\mathbf{D} & 0 \end{bmatrix} \begin{bmatrix} \mathbf{u}_n \\ \mathbf{p}_n \end{bmatrix} = \begin{bmatrix} \mathbf{M}\mathbf{t}_n \\ 0 \end{bmatrix} \quad (36) $$

and introduce the auxiliary matrix **Q** to rewrite the Stokes system in the following way:

$$ \begin{bmatrix} \mathbf{H} & -\mathbf{HQ}\mathbf{D}^T \\ -\mathbf{D} & 0 \end{bmatrix} \begin{bmatrix} \mathbf{u}_n \\ \delta\mathbf{p} \end{bmatrix} = \begin{bmatrix} \mathbf{M}\mathbf{t}_n + \mathbf{D}^T \mathbf{p}_{n-1} \\ 0 \end{bmatrix} + \begin{bmatrix} \mathbf{r} \\ 0 \end{bmatrix}, \quad (37) $$

where $\delta \mathbf{p} = \mathbf{p}_n - \mathbf{p}_{n-1}$ is the pressure increment, and the residual term is

$$ \mathbf{r} = -(\mathbf{H}\mathbf{Q} - \mathbf{I})\mathbf{D}^T \delta\mathbf{p}, \quad (38) $$

in which **I** is the identity matrix. If $\mathbf{Q} = \mathbf{H}^{-1}$, we retrieve the standard (expensive) Uzawa system. On the other hand, taking $\mathbf{Q} = \frac{\Delta t}{3/2}\mathbf{M}^{-1}$ is a computationally convenient choice, as **M** is diagonal. It can be shown that, in this case, neglecting **r** in (37) leads to a method which is formally second-order accurate in time (Fischer 1997 and references therein), and therefore does not affect the overall accuracy of the time scheme. This is the option we retain. Dropping the residual and carrying out a round of block Gaussian eliminations leads to the reformulated Stokes problem

$$ \begin{bmatrix} \mathbf{H} & -\frac{\Delta t}{3/2}\mathbf{H}\mathbf{M}^{-1}\mathbf{D}^T \\ 0 & \mathbf{E} \end{bmatrix} \begin{bmatrix} \mathbf{u}_n \\ \delta\mathbf{p} \end{bmatrix} = \begin{bmatrix} \mathbf{M}\mathbf{t}_n + \mathbf{D}^T \mathbf{p}_{n-1} \\ \mathbf{g} \end{bmatrix}, \quad (39) $$

where

$$ \mathbf{E} = \frac{\Delta t}{3/2}\mathbf{D}\mathbf{M}^{-1}\mathbf{D}^T, \quad (40) $$

and

$$ \mathbf{g} = -\mathbf{D}\mathbf{H}^{-1}(\mathbf{M}\mathbf{t}_n + \mathbf{D}^T \mathbf{p}_{n-1}). \quad (41) $$

$\mathbf{E}$ is directly proportional to $\mathbf{D}\mathbf{M}^{-1}\mathbf{D}^T$, also known as the pseudo-Laplacian operator (Maday *et al.* 1993), and $\mathbf{g}$ is the so-called inhomogeneity. To summarize, the procedure we follow at each time step consists of first computing $\mathbf{g}$, by inverting the modified Helmholtz operator **H**. In other words, we treat viscous and rotational effects implicitly

$$ \mathbf{g} = -\mathbf{D}\mathbf{H}^{-1}(\mathbf{M}\mathbf{t}_n + \mathbf{D}^T \mathbf{p}_{n-1}) = -\mathbf{D}\mathbf{u}^*, \quad (42) $$
where $\mathbf{u}^*$ can be interpreted as a first guess for the velocity. $\mathbf{H}$, which is not symmetric, is inverted iteratively, using a preconditioned stabilized biconjugate gradient algorithm (Van Der Vorst 1992). The preconditioner is of the element-by-element kind (Wathen 1989), and proves to be efficient enough, as $\mathbf{H}$ is diagonally dominant. The pressure increment $\delta \mathbf{p} = \mathbf{p}_n - \mathbf{p}_{n-1}$ that enforces the incompressibility constraint is then calculated by inverting $\mathbf{E}$

$$ \delta \mathbf{p} = \mathbf{E}^{-1} \mathbf{g}. \qquad (43) $$

As $\mathbf{E}$ is symmetric (see eq. 40), eq. (43) is solved iteratively by means of a preconditioned conjugate gradient algorithm. The final velocity $\mathbf{u}_n$ follows from

$$ \mathbf{u}_n = \frac{\Delta t}{3/2} \mathbf{M}^{-1} \mathbf{D}^T \delta \mathbf{p} + \mathbf{u}^*. \qquad (44) $$

This splitting approach is similar to classical splitting techniques, such as the fractional time step method originally devised by Chorin (1968). It differs nevertheless, in that the splitting is effected in the discrete form of the equations. Unlike a fractional step method, no additional pressure boundary conditions are thus prescribed, and no temporal error is introduced. We should note that inconsistent pressure boundary conditions tend to create so-called *divergence boundary layers*, located near the domain boundary $\partial \Omega$ (Blair Perot 1993; Tomboulides 1993, and references therein). As rotating fluids embedded in a container with rigid boundaries tend to generate sharp boundary layers that can in turn influence the bulk flow (Greenspan 1990), we would rather avoid to generate inconsistent boundary layers. Our strategy permits this, albeit at a somewhat larger cost than standard splitting schemes. Indeed, as pointed out by Maday et al. (1993), the pseudo-Laplacian involved in (43) has a much worse condition than the standard Laplacian that follows from a fractional step approach. It is therefore more difficult to invert iteratively. This problem can, however, be alleviated using an additive overlapping Schwarz preconditioner, which we describe in Appendix C.

Also, note that the examples that follow correspond to linear problems. The implicit technique described above being unconditionally stable, there is no stability constraint on the time step size—this is precisely why one tries to treat as many terms as possible in an implicit fashion. In a non-linear situation, however, the explicit treatment of the non-linear terms (following for instance an Adams-Bashforth formula) controls the maximum value of the time step that one can choose. The reader is referred to the book by Deville et al. (2002), chapters 3 and 6, for an extensive treatment of this issue in the spectral-element framework.

## 6 SEM VERSUS ANALYTIC SOLUTIONS: STEADY AND UNSTEADY STOKES PROBLEMS

We verify the accuracy of our implementation of the SEM by comparing it to a set of analytical solutions in a spherical shell configuration: the outer ($r_o$) and inner ($r_i$) radii are chosen such that $r_i/r_o$ is equal to 1/3. The basic idea behind our analytical tests is to define a simple reference divergence-free velocity, and to compute the body force that ensures conservation of momentum. In other words, we solve the forward problem, where a known velocity field is used to analytically compute the forcing of the right-hand side, and we then use this forcing as an input to our SEM code, in order to retrieve the velocity field numerically.

### 6.1 Steady stokes problem

In a first series of tests, we disregard inertia and the effects of rotation to consider a steady Stokes problem. The goal of this test is twofold: first, we wish to verify that the proper spaces are used to discretize velocity and pressure in our implementation of the $P_N - P_{N-2}$ approach, that is, we wish to confirm that the SEM velocity is indeed divergence-free. Second, we also wish to retrieve the classical spectral convergence properties of spectral methods.

The steady Stokes problem reads:

$$ \Delta \mathbf{u} - \nabla p + b \hat{\mathbf{r}} = \mathbf{0} \quad \text{in } \Omega, \qquad (45a) $$

$$ \nabla \cdot \mathbf{u} = \mathbf{0} \quad \text{in } \Omega, \qquad (45b) $$

$$ \mathbf{u} = \mathbf{0} \quad \text{on } \partial \Omega. \qquad (45c) $$

Note that $\hat{\mathbf{r}}$ is the unit vector in the radial direction and that the prescribed forcing $b\hat{\mathbf{r}}$ we seek in eq. (45a) is purely radial (it could be interpreted as an imposed buoyancy force). To define the analytical reference solution $\mathbf{u}_a$, we start by making the standard poloidal-toroidal decomposition of the velocity (e.g. appendix B of Dahlen & Tromp 1998):

$$ \mathbf{u}_a = \nabla \times (E^s \hat{\mathbf{r}}) + \nabla \times \nabla \times (F^s \hat{\mathbf{r}}), \qquad (46) $$

where $E^s$ and $F^s$ are the toroidal and poloidal fields, respectively, and where the superscript 's' stands for 'steady'. Using this expansion, we automatically satisfy eq. (45b) with our reference velocity solution. Each field is then sought in terms of zonal spherical harmonics

$$ \{E^s, F^s, b, p\} = \sum_l \{E_l^s, F_l^s, b_l, p_l\} (r) L_l(\cos \theta), \qquad (47) $$

in which $L_l$ is the Legendre polynomial of degree $l$. As the problem of interest is linear, we can consider one harmonic at a time. The radial components of the first and second curls of (45a) reduce to

$$ E_l^s = 0, \qquad (48) $$

$$ \left[ \frac{d^2}{dr^2} - \frac{l(l+1)}{r^2} \right]^2 F_l^s = b_l. \qquad (49) $$

By seeking a purely radial forcing, the toroidal field is identically zero. As far as the poloidal field is concerned, since (49) is a fourth-order equation, we need four boundary conditions, two each at $r_i$ and $r_o$. The no-slip boundary conditions imply that

$$ F_l^s = \frac{d}{dr} F_l^s = 0 \quad \text{at } r = r_i, r_o. \qquad (50) $$

The procedure for our test is then as follows:

(1) choose an expression for $F_l^s$ that matches the boundary conditions (50);

(2) solve eq. (49) analytically for the appropriate forcing $b_l$;

(3) use this forcing as an input for the SEM code;

(4) solve the Stokes problem using the SEM, starting from a zero initial guess for velocity and pressure;

(5) quantify the accuracy of the numerical solution $\mathbf{u}_h$ with respect to the analytical solution $\mathbf{u}_a$.

The Stokes problem is solved here with a standard Uzawa algorithm (Arrow et al. 1958), and the mesh we use is represented in Fig. 3. It consists of $n_e = 4$ spectral elements. Note that we also vary the polynomial order $N$ from 4 to 12 in our test, and that Fig. 3 corresponds to the coarsest mesh with $N = 4$. Depending on the
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Left, log-log plot of relative error (in an $L_1^2$ sense) for the steady Stokes problem (see the text) as function of polynomial order $N$ for harmonic degrees $l = 1, 3$ and $5$ (top to bottom). Note spectral convergence as $N$ increases. Right, SEM solution $\mathbf{u}_h$ for the same harmonic degrees obtained using $N = 11$. The analytic reference solution is not shown here, at it is indistinguishable from the SEM solution.</figcaption></figure>

spherical harmonic degree $l$ of the input velocity field, we either enforce a zero vertical velocity (when $l$ is even), or a zero radial velocity (when $l$ is odd) at the equator. Results for the $l = 1, 3, 5$ harmonics are displayed in Fig. 4. In each case the relative error

$$ \|\mathbf{e}\| = \left[ \frac{\int_{\Omega}(\mathbf{u}_h - \mathbf{u}_a)^2 d\Omega}{\int_{\Omega}\mathbf{u}_a^2 d\Omega} \right]^{1/2} \quad (51) $$

is very small (below one per cent for all cases with $N > 5$). Moreover, it decreases exponentially with the polynomial order $N$. It therefore exhibits the expected spectral convergence properties of classical spectral methods. Indeed, when we increase the polynomial order $N$, we find that the accuracy of the numerical solution is only limited

by the regularity of the solution sought. This behaviour validates our implementation of the $P_N - P_{N-2}$ method, and furthermore confirms that no spurious pressure modes are present that would prevent the velocity from being divergence-free. Meeting this *sine qua non* requirement enables us to turn our attention to time-dependent problems.

## 6.2 Unsteady Stokes problem

We now assess the temporal error of the time-marching scheme. The procedure is identical to the one we followed in the previous subsection, save that we introduce temporal variations. In other words, over the time interval $[0, T]$ we now consider an unsteady Stokes
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Top, SEM mesh (left) with $n_e = 4$ spectral elements of order $N = 11$ used for the unsteady Stokes problem (see the text) together with snapshot (right) of harmonic degree 1 velocity field. Bottom, log-log variation of relative error $||e||(T)$ versus time step size for five different values of $\Delta t$ (white circles). The black fit line corresponds to an exponent of 1.94, which is close to the expected theoretical value (2) for a second-order accurate time-scheme.</figcaption></figure>

problem of the form:

$$ \Delta \mathbf{u} - \nabla p + b\hat{\mathbf{r}} = \partial_t \mathbf{u} \quad \text{in } \Omega \quad \forall t \in [0, T] \qquad (52a) $$

$$ \nabla \cdot \mathbf{u} = 0 \quad \text{in } \Omega \quad \forall t \in [0, T] \qquad (52b) $$

$$ \mathbf{u} = \mathbf{0} \quad \text{on } \partial\Omega \quad \forall t \in [0, T], \qquad (52c) $$

supplemented by the initial condition $\mathbf{u} = \mathbf{0}$ at $t = 0$. In an attempt to focus our attention exclusively on temporal errors, we seek to ensure that spatial errors are negligible in this benchmark. To this end, we consider zonal harmonic $l = 1$ and choose a mesh of $n_e = 4$ elements having a rather high polynomial order of $N = 11$ (see Fig. 5, top left). We recall that this fine mesh resulted in a spatial error of $3.2 \times 10^{-7}$ in our earlier steady benchmark case (Fig. 4, top left). The negligible spatial error guarantees that our solution will be dominated by temporal error due to the time-marching scheme. The reference velocity field $\mathbf{u}_a$, and the forcing to conserve momentum in eq. (52a) are determined as before. While the toroidal component of $\mathbf{u}_a$ is still zero, we define its unsteady poloidal component as:

$$ F_l^u(r, t) = \sin(2\pi t/\tau) F_l^s(r), \qquad (53) $$

meaning that we let the steady-state solution from the previous sub-section oscillate with some arbitrary period $\tau$. The time-dependent force field consistent with this velocity can be used again as an input to the SEM code. After a time $T$ larger than $\tau$, we evaluate the

normalized error

$$ \|\mathbf{e}\|(T) = \left\{ \frac{\int_{\Omega} [\mathbf{u}_h(\mathbf{x}, T) - \mathbf{u}_a(\mathbf{x}, T)]^2 d\Omega}{\int_{\Omega} \mathbf{u}_a^2(\mathbf{x}, T) d\Omega} \right\}^{1/2} \qquad (54) $$

We repeat this procedure for various values of the time step $\Delta t$, and display the results in the bottom curve of Fig. 5. The error level (always above $10^{-5}$) is dominated by temporal error, as expected. The largest $\Delta t$ has a value equal to $\tau/(10\pi)$, while smaller $\Delta t$s follow a geometrical sequence of common ratio 1/2. The error level is proportional to the time step size with a power close enough to the expected value (1.94 versus 2) to confirm that neglecting the residual term in eq. (37) does not affect the overall order 2 accuracy of the time-scheme.

# 7 SEM VERSUS EXISTING NUMERICAL SOLUTIONS: THE PROUDMAN-STEWARTSON PROBLEM

## 7.1 Description

We conclude our presentation of the SEM by applying it to a simple flow problem more relevant to geophysical situations. As shown in Fig. 6(a), we consider a reference frame rotating at a constant rate $\omega$, where flow is induced inside a spherical shell by the super rotation of the inner sphere. We assume that the effects of rotation dominate the viscous effects, which corresponds to $E \ll 1$ in eq. (1a). When the super rotation is small enough, the solution is steady and axisymmetric (Proudman 1956). Moreover, away from viscous boundary layers, the velocity must obey the Taylor-Proudman theorem, i.e. it must be invariant along the axis of rotation

$$ \partial_t \mathbf{u} = \mathbf{0}. \qquad (55) $$

The Taylor-Proudman theorem leads to different flow regimes inside and outside of an imaginary cylinder $C$, that circumscribes the inner sphere and is aligned parallel to the axis of rotation. This cylinder, commonly referred to as the *tangent cylinder*, is represented by a dotted line in Fig. 6(a). Outside of $C$ a fluid particle is insensitive to the super rotation of the inner sphere. It therefore stays at rest with respect to the background rotation $\omega$. Inside of $C$, however, a fluid particle senses the super rotation of the inner sphere, and is entrained in its direction. The background rotation induces via the Coriolis force a meridional circulation that is controlled by pumping and suction inside the viscous Ekman boundary layers located at the inner and the outer shell boundaries (see Fig. 6a). The circulation is completed alongside of $C$, where a viscous shear layer (referred to as the *Stewartson layer*) accommodates the angular velocity jump between regions inside and outside $C$.

This classic kinematic flow problem in its asymptotic form was originally proposed by Proudman (1956). It was developed later by Stewardson (1957, 1966), who derived the exact structure of the shear layers alongside of $C$. Since then it has been treated numerically by several groups (Hollerbach 1994; Dormy et al. 1998). In fact, though simple in essence, this so-called 'Proudman-Stewartson problem' exhibits the essential features of flows dominated by rotation. Its linear character, moreover, makes it an ideal test case to verify the numerical accuracy of our method, leaving aside complications that inevitably appear when one introduces non-linear effects. The Proudman-Stewartson problem is also quite challenging numerically, as one must resolve the narrow Ekman and Stewartson layers that characterize the circulation when $E \ll 1$.

To summarize, we want to study numerically the following set of equations:

$$ \partial_t \mathbf{u} + 2\hat{\mathbf{z}} \times \mathbf{u} = -\nabla p + E \Delta \mathbf{u} \quad \text{in } \Omega \qquad (56a) $$
<table><caption>**Table 1.** Summary of the Proudman–Stewartson problem results. Relative difference $||d||$ of SEM and spectral solution of Hollerbach (1994) for Ekman numbers $10^{-2}$, $10^{-3}$ and $10^{-4}$. $n_e$ is the number of elements, $N$ the polynomial order. Memory requirement is proportional to number of Gauss–Lobatto Legendre points ($N_{nodes}$), while computational cost scales approximately with $n_eN^3$. Note that $||d||$ is in all cases well below the one per mil level.</caption><thead><tr><th>E</th><th>n<sub>e</sub></th><th>N</th><th>N<sub>nodes</sub></th><th>n<sub>e</sub>N<sup>3</sup></th><th>||d||</th></tr></thead><tbody><tr><td>10<sup>-2</sup></td><td>4</td><td>11</td><td>529</td><td>5324</td><td>0.000 002</td></tr><tr><td>10<sup>-3</sup></td><td>9</td><td>11</td><td>1156</td><td>10 648</td><td>0.000 167</td></tr><tr><td>10<sup>-4</sup></td><td>140</td><td>8</td><td>9153</td><td>71 680</td><td>0.000 119</td></tr></tbody></table>

$$ \nabla \cdot \mathbf{u} = 0 \text{ in } \Omega \qquad (56b) $$

$$ \mathbf{u}(\mathbf{x}, t) = s\hat{\phi} \text{ at } r = r_i \qquad (56c) $$

$$ \mathbf{u}(\mathbf{x}, t) = \mathbf{0} \text{ at } r = r_o \qquad (56d) $$

$$ \mathbf{u}(\mathbf{x}, 0) = \mathbf{0} \text{ in } \Omega, \qquad (56e) $$

for values of the Ekman number ranging from $E = 10^{-2}$ to $10^{-4}$.

## 7.2 Reference numerical solution

First, we describe the reference numerical solution together with the changes in the physics of the solution, as we go from a slowly rotating system to a system rotating more rapidly. Our numerical reference solution is the spectral solution published by Hollerbach (1994). It is displayed in Fig. 6 for three values of $E$, $10^{-2}$ (Fig. 6b), $10^{-3}$ (Fig. 6c) and $10^{-4}$ (Fig. 6d). Hollerbach (1994) computed solutions

to the steady problem, by solving system (56) withdrawing the time-derivative in eq. (56a). His spectral method satisfies the solenoidal constraint exactly, as the meridional circulation is described in terms of a stream function, and in his approach the field variables are expanded by means of Chebyshev and Legendre polynomials in the radial and angular direction, respectively. His published results include solutions for an Ekman number as low as $10^{-5}$. These results were subsequently confirmed by Dormy *et al.* (1998), who solved the Stewartson problem for $E$ as small as $2.37 \times 10^{-8}$, using Legendre polynomials in the angular direction together with a finite difference scheme in radius.

Looking at Fig. 6(b), i.e. at the stream function and angular velocity of the flow at $E = 10^{-2}$, we note that the solution is still essentially viscous and smooth. When $E$ is decreased by an order of magnitude (Fig. 6c), Ekman layers form in the vicinity of the shell boundaries. Away from these Ekman layers angular velocity contours and meridional circulation tend to align themselves with the axis of rotation, although the Taylor–Proudman theorem is not fully satisfied as viscous effects are still noticeable. When we decrease $E$

<figure><img src="image_8.png" /><figcaption>**Figure 6.** (a) Sketch (left) of the flow in the Proudman–Stewartson problem (see the text), in a rapidly rotating spherical shell. The flow is induced by a slight super rotation $\Delta\omega$ of the inner sphere. Motion is generated inside the tangent cylinder $C$ only, as the fluid located outside $C$ is insensitive to the super rotation by virtue of the Taylor–Proudman theorem. Meridional circulation inside $C$ is controlled by Ekman pumping/suction occurring at both boundaries, and is completed alongside $C$. A shear layer is created along $C$ to accommodate the shear between the two regions. Numerical solutions (right) to the Proudman–Stewartson flow problem, computed by Hollerbach (1994), for $E = 10^{-2}$, $10^{-3}$ and $10^{-4}$ (b, c, d, top to bottom). On the left contour plots of the meridional stream function are shown, while the right shows contour plots of the angular velocity. Note the transition from a smooth, essentially viscous solution (b) to a solution dominated by rotation (d).</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** SEM solutions to the Proudman–Stewartson problem (see the text) for $E = 10^{-2}$, $10^{-3}$, and $10^{-4}$ (top to bottom). From left to right, computational mesh with number of elements ($n_e$) and polynomial order ($N$); velocity map of meridional circulation ($u_s, u_z$); and contour plot of angular velocity $u_\phi/s$, with contour intervals of 1/15. From top to bottom, $100|u_s|_{\max} = 4.67, 5.02, 6.13$ and $100|u_z|_{\max} = 6.05, 7.60, 7.92$. Note that the relative difference $\|d\|$ of SEM versus Hollerbach (indicated at centre of the shell) is well below one per mil in all cases.</figcaption></figure>

further to a value of $E = 10^{-4}$, the Ekman layers sharpen significantly and the Stewartson shear layer starts to develop (Fig. 6d, right). For this value of the Ekman number $E$, the circulation inside the tangent cylinder follows the Taylor–Proudman theorem quite well (Fig. 6d, left). This last case is the most challenging numerically.

## 7.3 SEM solution to the Stewartson problem

We now turn our attention to the SEM solution of the Stewartson problem for the same three Ekman numbers considered before. Our

time-dependent calculations start with a fluid initially at rest. We then time step our code until the flow settles to steady-state, and compare this steady-state solution $u_h$ to the solutions obtained by Hollerbach. To compute the relative difference between the two numerical solutions:

$$ \|d\| = \left[ \frac{\int_{\Omega} (u_h - u_s)^2 d\Omega}{\int_{\Omega} u_s^2 d\Omega} \right]^{1/2}, \qquad (57) $$

we project the purely spectral solution on the spectral-element grid, and consider agreement between the two numerical approaches as acceptable at levels better than $\|d\| < 10^{-3}$.
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Relative difference $||d||$ of SEM versus Hollerbach (1994) spectral solution to Proudman-Stewartson problem (see the text) for $E = 10^{-3}$ as function of polynomial order $N$ for computational mesh with and without grid-refinement, as indicated by inlet-figures. The mesh with uniform resolution consists of $n_e = 4$ spectral elements, while the non-uniform mesh includes $n_e = 16$ elements, concentrated in Ekman layers at both shell boundaries. The $n_e = 16$ mesh (bottom curve) shows lower overall difference, as expected, due to its larger number of elements. Note, however, that this mesh is characterized by spectral convergence up to the highest polynomial order ($N = 12$) employed, while the uniform resolution mesh (top curve) shows a flattening of the convergence, indicating that mesh-refinement is well suited to describe the Ekman layers and results in an overall faster convergence rate.</figcaption></figure>

The SEM results (angular velocity and meridional circulation) for $E = 10^{-2}$, $10^{-3}$ and $10^{-4}$ are displayed on the right-hand side of Fig. 7 and summarized in Table 1. Recall that unlike Hollerbach (1994), we do not describe the meridional circulation by means of a stream-function in our method, and that we solve instead directly for the primitive variables ($u_s, u_z$). We also show on the left-hand side of Fig. 7 the computational mesh used to obtain the solution.

For the large Ekman number case ($E = 10^{-2}$, Figs 7a-c) four spectral elements are needed to retrieve the solution. Setting the polynomial order to 11 leads by virtue of the spectral convergence properties of the SEM (see Section 6.1) to excellent agreement with Hollerbach's solution, as the relative difference between the two methods for this smooth solution is $||d|| = 0.000002$. For the medium and low Ekman number cases ($E = 10^{-3}$ and $10^{-4}$, Figs 7d-i), we reproduce Hollerbach's solution with our required level of accuracy ($||d|| < 10^{-3}$) at the expense of an increased numerical resolution (see Table 1). In the latter case, we use 140 elements of polynomial order 8 (see Fig. 7g). These elements are gently squeezed toward the boundaries of the shell. If one were to squeeze them further, the contrast in aspect ratio between neighbouring elements and the resulting anisotropy in the computational mesh would lead to a deterioration of condition in the algebraic system (39) and thus to a larger computational cost. We do not consider smaller Ekman number solutions here, because their computational requirements (integration time and computer memory) would be prohibitive for single processor runs using standard personal computers. However, we are currently parallelizing our code to enable computations at lower Ekman numbers.

## 7.4 Adaptivity and enhanced convergence

We complete our study of the SEM by examining its potential for grid-refinement. Indeed, the local character of the SEM allows for a

non-uniform paving of the domain. We made use of this property in our study of the Stewartson problem by radially squeezing spectral elements in the vicinity of the Ekman boundary layers. More generally, the adaptivity of grid-based methods, such as the SEM, is well suited to capture strong gradients by providing increased numerical solution in regions containing fine-scale structures.

To illustrate this property in a quantitative way, we consider the solution of the Stewartson problem at $E = 10^{-3}$ obtained with two different meshes (portrayed in Fig. 8) for a number of polynomial orders $N$. As we have seen before, the value of $E$ in this intermediate case is small enough to permit the development of strong Ekman layers. We want to quantify how well we succeed in describing these localized features of the Proudman-Stewartson solution by using either an almost uniform mesh of $n_e = 4$ elements, or a non-uniform mesh having $n_e = 16$ spectral elements where the exterior elements are squeezed by a factor of three relative to interior elements. We expect, for the same polynomial order $N$, that our mesh with local grid-refinement should yield better results than our mesh of almost uniform grid resolution due to its larger number of elements. However, more importantly, we also expect that our ability to adapt the mesh geometry to conform to the inherent flow structure, with narrow Ekman layers concentrated near the shell boundaries, should furthermore lead to a faster convergence rate when we increase the polynomial order. This behaviour is indeed confirmed by the two curves in Fig. 8 showing the accuracy of our solution as a function of the polynomial order $N$. The top curve (for the uniform mesh) tends to flatten as we increase the polynomial order above 10, while the bottom curve (for the non-uniform mesh) does not show any sign of flattening.

## 8 DISCUSSION

In this paper we have presented the application of the spectral-element method to model axisymmetric flows in a rapidly rotating reference frame. Spatial discretization relies on breaking the physical domain into a collection of $n_e$ non-overlapping elements, and using a local tensorized basis of polynomials of high order $N$. The temporal discretization is achieved by means of a second-order accurate backward differentiation scheme for the time-derivative, and we follow an operator integrated factor (OIF) approach to treat the resulting modified Stokes problem at each time step. Our use of the OIF strategy is motivated by the fact that the mass matrix arising from the spatial discretization is diagonal. Our implementation of the OIF, relative to a standard implementation of the OIF (Maday et al. 1990) where only viscous effects are treated, is modified such as to account for the effects of the Coriolis force. Hence both viscous and Coriolis effects are treated fully implicitly in our approach. The computation of the pressure increment is accomplished using an efficient additive overlapping Schwarz preconditioner.

We have validated our implementation of the method by comparing it to analytical and published spectral reference solutions of axisymmetric laminar flows in a spherical shell. Our steady-state analytical benchmark solutions demonstrate that no spurious pressure mode exists due to the discretization of the pressure using a lower polynomial order ($N-2$) relative to the polynomial order ($N$) that is used for the velocity. The $P_N - P_{N-2}$ approach for velocity and pressure therefore guarantees the solenoidal character of the velocity field. These results furthermore served to demonstrate the spectral convergence properties of the SEM. The unsteady analytic benchmark solution allowed us finally to verify that the relatively unusual time discretization strategy adopted here is second-order accurate.
<figure><img src="image_11.png" /><figcaption>**Figure 9.** 3-D spherical mesh for mixed Fourier–spectral-element approach (see the text). The mesh consists of 32 spectral elements in the meridional plane together with 32 equally spaced points in longitude. The spatial domain decomposition (indicated by grey-scale) for parallel processing is illustrated for 16 processors handling two spectral elements each.</figcaption></figure>

We have used published spectral solutions to the kinematic Proudman–Stewartson problem (Hollerbach 1994) at Ekman numbers of $10^{-2}$, $10^{-3}$ and $10^{-4}$ to assess the efficacy of the SEM in a rapidly rotating context. These solutions, which exhibit some of the essential features of flows dominated by rotation, confirm consistency between the two different numerical techniques at levels of better than one per mil. In fact, our calculations for the Stewartson problem demonstrate for all practical purpose that the SEM solutions are identical to the solution published by Hollerbach. This close correspondence between SEM and spectral solution shows that our inclusion into the SEM of the effects of a dominant Coriolis force, a relatively novel extension to standard applications of the SEM, does not affect the overall accuracy of the method.

We have shown that the local character of the SEM is well suited for a non-uniform paving of the computational domain, i.e. for grid-refinement in regions of the flow characterized by strong local gradients such as Ekman boundary layers at the core–mantle boundary (CMB) and the inner-core boundary (ICB). In fact, exploiting mesh-refinement in the resolution of the Proudman–Stewartson problem leads to a substantial improvement in the convergence properties of the method. Given the high numerical accuracy of the SEM, its local and hence inherently parallel character, as well as its adaptivity, these early results are encouraging in the long-term prospect of building a parallel spectral-element based model of the geodynamo.

The extension of the present kernel to 3-D problems can be achieved through a mixed Fourier–spectral-element approach (Bernardi et al. 1999), where the longitudinal dependence of the variables is expanded into Fourier series. The approach, effectively, breaks the 3-D problem into a set of meridional problems quite similar to those presented in this paper. The special treatment of spectral elements adjacent to the axis of rotation we introduced is well suited for this technique. Moreover, the Schwarz preconditioner we use here to compute the pressure increment at each time

step lends itself naturally to computationally demanding problems in 3-D applications.

The local character of the SEM allows for a straightforward spatial domain decomposition in the meridional plane using active message-passing to communicate among subdomains, i.e. for interprocessor communication. The strategy is illustrated in Fig. 9, where we show an example of a 3-D mesh, grey-scale coded for parallelization among processors. This approach should allow us, in the future, to take advantage of the gain in computational speed on modern parallel computers.

Our method should carry over to the implementation of the magnetic component of the geodynamo problem. For example, imposing the divergence-free constraint on the magnetic field **B** can be accomplished using the same approach we used to enforce the divergence-free constraint on the velocity field **u**. Likewise, accounting for the effects of the inner core (Hollerbach & Jones 1993) can be accomplished by extending our computational domain beyond the ICB radially inwards to include the solid inner core. The main complexity related to the magnetic induction equation lies in imposing the magnetic boundary conditions at the CMB, where the magnetic field has to be connected to an exterior potential field. This matching condition translates to an elegant analytic boundary condition in the context of spherical harmonics (Glatzmaier 1984). It is, however, less amenable to a local method, and it is yet to find out how this matching is to be accomplished. An alternative to this matching has been presented by Chan et al. (2001), who meshed the external space up to a few core radii and used approximate magnetic boundary conditions on the outer surface of the computational domain. We are currently investigating different strategies regarding the implementation of the induction equation and we hope to benefit, in a non-linear context, from the weak numerical dispersion of the SEM.

## ACKNOWLEDGMENTS

AF thanks E. Chaljub for his insightful comments over the course of this work (along with a careful reading of an earlier version of the manuscript), G. Glatzmaier for his comments on an earlier version of the manuscript, E. Dormy and Y. Maday for fruitful discussions, and M. Azaïez for inviting him to the EUROMECH-ERCOFTAC colloquium 'Spectral Methods and Time Stepping Schemes For Incompressible Flows in Complex Geometries' held in Toulouse in 2001 October.

Plots and analyses were made using the following freely available software: the generic mapping tools (Wessel & Smith 1991), gnuplot (www.gnuplot.info/), xmgrace (http://plasma-gate.weizmann.ac.il/Grace/) and the pstricks package (www.tug.org/applications/PSTRicks/).

This work was supported by NSF grant EAR-01-13050.

## REFERENCES

* Arrow, K., Hurwicz, L. & Uzawa, H., 1958. *Studies in Nonlinear Programming*, Stanford University Press, Stanford.

* Bernardi, C. & Maday, Y., 1988. A collocation method over staggered grids for the Stokes problem, *Int. J. Num. Meth. Fluids*, **8**, 537–557.

* Bernardi, C. & Maday, Y., 1992. *Approximations Spectrales de Problèmes aux Limites Elliptiques*, *Mathématiques & Applications*, **Vol. 10**, Springer-Verlag, Paris.

* Bernardi, C., Dauge, M. & Maday, Y., 1999. *Spectral Methods for Axisymmetric Domains*, *Series in Applied Mathematics*, **Vol. 3**, Gauthier-Villars, Paris, numerical algorithms and tests due to Mejdi Azaïez.

* Blair Perot, J., 1993. An analysis of the fractional step method, *J. Comput. Phys.*, **108**, 51–58.
Bloxham, J., 1998. Dynamics of angular momentum of the Earth's core, *Annu. Rev. Earth Planet. Sci.*, **26**, 501–517.

Bloxham, J., 2000a. The effect of thermal core–mantle interactions on the palaeomagnetic secular variation, *Phil. Trans. R. Soc. Lond. A*, **358**, 1171–1179.

Bloxham, J., 2000b. Sensitivity of the geomagnetic axial dipole to thermal core–mantle interactions, *Nature*, **405**, 63–65.

Bunge, H.P. & Baumgardner, J.R., 1995. Mantle convection modeling on parallel virtual machines, *Comp. Phys.*, **9**, 207–215.

Bunge, H.-P. & Dalton, M., 2001. Building a high-performance linux cluster for large-scale geophysical modeling, in *Linux Clusters: the HPC Revolution*, NCSA Conf. Proc.

Busse, F., 2000. Homogeneous dynamos in planetary cores and in the laboratory, *Annu. Rev. Fluid. Mech.*, **32**, 383–408.

Capdeville, Y., Chaljub, E., Vilotte, J.-P. & Montagner, J.-P., 2002. Coupling the spectral element method with a modal solution for elastic wave propagation in global Earth models, *Geophys. J. Int.* **152**, 34–67.

Chaljub, E., Capdeville, Y. & Vilotte, J.-P., 2003. Solving elastodynamics in a fluid–solid heterogeneous sphere: a parallel spectral element approximation on non-conforming grids, *J. Comput. Phys.*, **152**, 457–491.

Chan, K.H., Zhang, K., Zou, J. & Schubert, G., 2001. A non-linear, 3-D spherical $\alpha^2$ dynamo using a finite element method, *Phys. Earth planet. Inter.*, **128**, 35–50.

Chorin, A.J., 1968. Numerical solution of the Navier–Stokes equations, *Math. Comput.*, **22**, 745–762.

Christensen, U., Olson, P. & Glatzmaier, G.A., 1999. Numerical modeling of the geodynamo: a systematic parameter study, *Geophys. J. Int.*, **138**, 393–409.

Couzy, W., 1995. Spectral element discretization of the unsteady Navier–Stokes equations and its iterative solution on parallel computers, *PhD thesis*, École Polytechnique Fédérale de Lausanne.

Dahlen, F.A. & Tromp, J., 1998. *Theoretical Global Seismology*, Princeton University Press, Princeton.

Deville, M.O., Fischer, P.F. & Mund, E.H., 2002. *High-Order Methods for Incompressible Fluid Flow*, Cambridge Monographs on Applied and Computational Mathematics, Vol. 9, Cambridge University Press, Cambridge.

de Wijs, G.A., Kress, G., Vocado, L., Dobson, D., Alfè, D., Gillan, M.J. & Price, G.D., 1998. The viscosity of liquid iron at the physical conditions of the Earth's core, *Nature*, **392**, 805–807.

Dormy, E., Cardin, P. & Jault, D., 1998. MHD flow in a slightly differentially rotating spherical shell, with conducting inner core, in a dipolar magnetic field, *Earth planet. Sci. Lett.*, **160**, 15–30.

Dormy, E., Valet, J.-P. & Courtillot, V., 2000. Numerical models of the geodynamo and observational constraints, *Geochem. Geophys. Geosyst.*, **1**(62).

Fischer, P.F., 1997. An overlapping Schwarz method for spectral element solution of the incompressible Navier–Stokes equations, *J. Comput. Phys.*, **133**, 84–101.

Fischer, P.F., Miller, N. & Tufo, H., 2000. An overlapping Schwarz method for spectral element simulation of three-dimensional incompressible flows, in *Parallel Solution of Partial Differential Equations*, pp. 159–181, eds Bjørstad, P. & Luskin, M., Springer-Verlag, Berlin.

Gerritsma, M. & Phillips, T., 2000. Spectral element methods for axisymmetric Stokes problems, *J. Comput. Phys.*, **164**, 81–103.

Giraldo, F.X., 2001. A spectral element shallow water model on spherical geodesic grids, *Int. J. Num. Meth. Fluids*, **35**, 869–901.

Glatzmaier, G.A., 1984. Numerical simulations of stellar convective dynamos I. The model and method, *J. Comput. Phys.*, **55**, 461–484.

Glatzmaier, G.A. & Roberts, P.H., 1995. A three-dimensional self-consistent computer simulation of a geomagnetic reversal, *Nature*, **377**, 203–209.

Glatzmaier, G.A. & Roberts, P.H., 1996a. An anelastic evolutionary geodynamo simulation driven by compositional and thermal convection, *Physica D*, **97**, 81–94.

Glatzmaier, G.A. & Roberts, P.H., 1996b. Rotation and magnetism of Earth's inner core, *Science*, **274**, 1887–1891.

Greenspan, H.P., 1990. *The Theory of Rotating Fluids*, 2nd edn, Breukelen Press, Brookline.

Grote, E., Busse, F.H. & Tilgner, A., 2000. Effects of hyperdiffusivities on dynamo simulations, *Geophys. Res. Lett.*, **27**, 2001–2004.

Gubbins, D. & Roberts, P.H., 1987. Magnetohydrodynamics of the Earth's core, in *Geomagnetism*, Vol. 2, ed Jacobs, J.A., Academic Press, London.

Hollerbach, R., 1994. Magnetohydrodynamic Ekman and Stewartson layers in a rotating spherical shell, *Proc. R. Soc. Lond., A*, **444**, 333–346.

Hollerbach, R., 2000. A spectral solution of the magneto-convection equations in spherical geometry, *Int. J. Num. Meth. Fluids*, **32**, 773–797.

Hollerbach, R. & Jones, C., 1993. Influence of the Earth's inner core on geomagnetic fluctuations and reversals, *Nature*, **365**, 541–543.

Hughes, T.J.R., 1987. *The Finite Element Method: Linear Static and Dynamic Finite Element Analysis*, Prentice-Hall, Englewood Cliffs.

Komatitsch, D. & Tromp, J., 1999. Introduction to the spectral element method for three-dimensional seismic wave propagation, *Geophys. J. Int.*, **139**, 806–822.

Komatitsch, D. & Vilotte, J.-P., 1998. The spectral element method: an effective tool to simulate the seismic response of 2D and 3D geological structures, *Bull. seism. Soc. Am.*, **88**, 368–392.

Kuang, W. & Bloxham, J., 1997. An Earth-like numerical dynamo model, *Nature*, **389**, 371–374.

Kuang, W. & Bloxham, J., 1999. Numerical modeling of magnetohydrodynamic convection in a rapidly rotating spherical shell: weak and strong field dynamo action, *J. Comput. Phys.*, **51**, 51–81.

Larmor, J., 1919. Possible rotational origin of magnetic fields of Sun and Earth, *Elec. Rev.*, **85**, 412.

Levin, J.G., Iskandarani, M. & Haidvogel, D.B., 2000. A nonconforming spectral element ocean model, *Int. J. Num. Meth. Fluids*, **34**, 495–525.

Maday, Y. & Patera, A.T., 1989. Spectral element methods for the incompressible Navier–Stokes equations, *State-of-the-Art Surveys on Computational Mechanics*, pp. 71–143, ASME.

Maday, Y., Patera, A.T. & Rønquist, E.M., 1990. An operator-integration-factor splitting method for time-dependent problems: applications to incompressible fluid flow, *J. Sci. Comp.*, **5**, 263–292.

Maday, Y., Meiron, D., Patera, A.T. & Rønquist, E.M., 1993. Analysis of iterative methods for the steady and unsteady Stokes problem: application to spectral element discretizations, *SIAM J. Sci. Comput.*, **14**, 310–337.

Patera, A.T., 1984. A spectral element method for fluid dynamics: laminar flow in a channel expansion, *J. Comput. Phys.*, **54**, 468–488.

Poirier, J.-P., 1988. Transport properties of liquid metals and viscosity of the Earth's core, *Geophys. J. R. astr. Soc.*, **92**, 99–105.

Proudman, I., 1956. The almost-rigid rotation of viscous fluid between concentric spheres, *J. Fluid. Mech.*, **1**, 505–516.

Reddy, J.N., 1993. *An Introduction to the Finite Element Method, Engineering Mechanics* McGraw-Hill, New York.

Roberts, P.H., 1978. Magnetoconvection in a rapidly rotating fluid, in *Rotating Fluid in Geophysics*, pp. 421–435, eds Roberts, P.H. & Soward, A.M., Academic, New York.

Smith, B., Bjørstad, P. & Gropp, W., 1996. *Domain Decomposition*, Cambridge University Press, Cambridge.

St. Pierre, M.G., 1993. The strong field branch of the Childress–Soward dynamo, in *Solar and Planetary Dynamos*, pp. 295–302, eds Proctor, M.R.E., Matthews, P.C. & Rucklidge, A.M., Cambridge University Press, Cambridge.

Stewartson, K., 1957. On almost rigid rotations, *J. Fluid. Mech.*, **3**, 299–303.

Stewartson, K., 1966. On almost rigid rotations part II, *J. Fluid. Mech.*, **26**, 131–144.

Taylor, M., Tribbia, J. & Iskandarani, M., 1997. The spectral element method for the shallow water equations on the sphere, *J. Comput. Phys.*, **130**, 92–108.

Tilgner, A., 1999. Spectral methods for the simulation of incompressible flows in spherical shells, *Int. J. Num. Meth. Fluids*, **30**, 713–724.

Tomboulides, A.G., 1993. Direct and large-eddy simulation of wake flows: flow past a sphere, *PhD thesis*, Princeton University.

Van Der Vorst, H., 1992. BI-CGSTAB: a fast and smoothly converging variant of BI-CG for the solution of nonsymmetric linear systems, *SIAM J. Sci. Statis. Comput.*, **13**, 631–644.
Wathen, A.J., 1989. An analysis of some element-by-element techniques, *Comp. Meth. Appl. Mech. Eng.*, **74**, 271–287.

Wessel, P. & Smith, W., 1991. Free software helps map and display data, *Trans. Am. Geophys. Union*, **72**, 441–445.

Zhang, K. & Gubbins, D., 2000. Scale disparities and magnetohydrodynamics in the Earth's core, *Phil. Trans. R. Soc. Lond. A*, **358**, 899–920.

Zhang, K. & Jones, C.A., 1997. The effect of hyperviscosity on geodynamo models, *Geophys. Res. Lett.*, **24**, 2869–2872.

# APPENDIX A: QUADRATURE FORMULAE AND POLYNOMIAL INTERPOLATION

## A1 Orthogonal polynomials in $L^2(\Lambda)$

Let us denote the Legendre polynomials of order $N$ with $L_N$. The Legendre polynomials are orthogonal in $L^2(\Lambda)$, that is

$$ \forall (N_1, N_2), \quad \int_{\Lambda} L_{N_1}(\xi)L_{N_2}(\xi)d\xi = \begin{cases} 0 & \text{if } N_1 \neq N_2 \\ 1/(N_1 + 1/2) & \text{if } N_1 = N_2 \end{cases} \qquad (A1) $$

They satisfy the conditions $L_N(1) = 1$ and $L_N(-1) = (-1)^N$. Each $L_N$ satisfies the following differential equation:

$$ [(1 - \xi^2)L'_{N}]' + N(N+1)L_N = 0. \qquad (A2) $$

Legendre polynomials are computed by means of the induction formula

$$ \begin{cases} L_0(\xi) = 1 \text{ and } L_1(\xi) = \xi, \\ (N+1)L_{N+1}(\xi) = (2N+1)\xi L_N(\xi) - N L_{N-1}(\xi), \quad N > 1. \end{cases} \qquad (A3) $$

## A2 Standard Gauss–Lobatto Legendre formula

We recall here the main properties of the standard Gauss–Lobatto Legendre formula. Let us set $\xi_i^N = -1$ and $\xi_i^N = 1$. Then there exists a unique set of $N-1$ nodes $\xi_i^N$, $1 \le i \le N-1$ in $\Lambda$ and of $N+1$ weights $\rho_i^N$, $0 \le i \le N$, such that the following exactness property holds

$$ \forall \Phi \in P_{2N-1}(\Lambda), \quad \int_{-1}^{1} \Phi(\xi) d\xi = \sum_{i=0}^{N} \rho_i^N \Phi(\xi_i^N). \qquad (A4) $$

The $\xi_i$, $1 \le i \le N-1$ are the zeros of $L'_N$ and the $\rho_i^N$ can be expressed as follows

$$ \rho_i^N = \frac{2}{N(N+1)L_N^2(\xi_i^N)}, \quad 0 \le i \le N. \qquad (A5) $$

A basis for $P_N$ is made of the Lagrangian interpolants $h_i^N$, $0 \le i \le N$ given by

$$ h_0^N(\xi) = (-1)^{N-1} \frac{(1-\xi)L'_N(\xi)}{N(N+1)}, \qquad (A6) $$

$$ h_i^N(\xi) = \frac{1}{N(N+1)L_N(\xi_i^N)} \frac{(1-\xi^2)L'_N(\xi)}{\xi_i^N - \xi}, \quad 1 \le i \le N-1, \qquad (A7) $$

$$ h_N^N(\xi) = \frac{(1 + \xi)L'_N(\xi)}{N(N+1)}. \qquad (A8) $$

The derivatives of these interpolants, which appear for instance in the divergence and the Laplacian bilinear forms can be estimated using eq. (A2). We can also derive the pressure basis functions as the set of Lagrangian interpolants $h_i^{N-2}$ defined by the interior Gauss–Lobatto points $\xi_i^N$, $1 \le i \le N-1$

$$ h_i^{N-2}(\xi) = \frac{1 - \xi_i^{N-2}}{N(N+1)L_N(\xi_i^N)} \frac{L'_N(\xi)}{\xi_i^N - \xi}, \quad 1 \le i \le N-1. \qquad (A9) $$

## A3 Orthogonal polynomials in $L_1^2(\Lambda)$

The weighted quadrature we use is based upon a class of polynomials $M_N$ defined by

$$ M_N(\xi) = \frac{L_N(\xi) + L_{N+1}(\xi)}{1 + \xi}, \quad N \ge 0. \qquad (A10) $$

They are orthogonal in $L_1^2(\Lambda)$, that is with the weighted measure $(1 + \xi) d\xi$

$$ \forall (N_1, N_2), \quad \int_{\Lambda} M_{N_1}(\xi) M_{N_2}(\xi)(1 + \xi) d\xi = \begin{cases} 0 & \text{if } N_1 \neq N_2 \\ 2/(N_1 + 1) & \text{if } N_1 = N_2 \end{cases} \qquad (A11) $$
Each $M_N$ satisfies $M_N(1) = 1$ and the differential equation

$$[(1 + \xi)^2(1 - \xi)M'_{N}] + N(N + 2)(1 + \xi)M_N = 0. \quad (A12)$$

The induction formula on the $M_N$ is

$$\left\{ \begin{array}{l} M_0(\xi) = 1 \text{ and } M_1(\xi) = \frac{1}{2}(3\xi - 1), \\ \frac{N+2}{2N+3} M_{N+1}(\xi) = \left(\xi - \frac{1}{(2N+1)(2N+3)}\right) M_N(\xi) - \frac{N}{2N+1} M_{N-1}(\xi), \quad N > 1. \end{array} \right. \quad (A13)$$

## A4 Weighted Gauss-Lobatto Legendre (WGLL) formula

We can now define a weighted quadrature formula of the Gauss-Lobatto kind. Let us set $\xi_0^N = -1$ and $\xi_N^N = 1$. Then there exists a unique set of $N-1$ nodes $\xi_i^N$, $1 \le i \le N-1$ in $\Lambda$ and of $N+1$ weights $\sigma_i^N$, $0 \le i \le N$, such that the following exactness property holds

$$\forall \Phi \in P_{2N-1}(\Lambda), \quad \int_{\Lambda} \Phi(\xi)(1+\xi) d\xi = \sum_{i=0}^{N} \sigma_i^N \Phi(\xi_i^N). \quad (A14)$$

The $\xi_i^N$, $1 \le i \le N-1$ are the zeros of $M'_N$ and the $\sigma_i$ can be expressed as follows

$$\sigma_0^N = \frac{8}{N(N+2)M_N^2(-1)}, \quad \text{and} \quad (A15)$$

$$\sigma_i^N = \frac{4}{N(N+2)M_N^2(\xi_i^N)}, \quad 1 \le i \le N. \quad (A16)$$

Again, under these circumstances, we can set as a basis for $P_N(\Lambda)$ the Lagrangian interpolants $I_i^N$, $0 \le i \le N$

$$I_0^N(\xi) = \frac{(\xi - 1)M'_N(\xi)}{2M'_N(-1) + N(N+2)M_N(-1)}, \quad (A17)$$

$$I_i^N(\xi) = \frac{1}{N(N+2)M_N(\xi_i^N)} \frac{(1-\xi^2)M'_N(\xi)}{\xi_i^N - \xi}, \quad 1 \le i \le N-1, \quad (A18)$$

$$I_N^N(\xi) = \frac{(1 + \xi)M'_N(\xi)}{N(N + 2)}. \quad (A19)$$

The pressure basis functions $I_i^{N-2}$, $1 \le i \le N-1$ are the Lagrangian interpolants on the interior nodes

$$I_i^{N-2}(\xi) = \frac{1 - \xi_i^{N2}}{N(N+2)M_N(\xi_i^N)} \frac{M'_N(\xi)}{\xi_i^N - \xi}, \quad 1 \le i \le N-1. \quad (A20)$$

## APPENDIX B: DERIVATION OF THE ALGEBRAIC SYSTEM

Here we explain in detail how the algebraic system (31) resulting from the spatial discretization by spectral elements is derived. Let us recall that the expansion of the velocity using the elemental basis functions anchored at the Gauss-Lobatto points write

$$\mathbf{u}_{i,h}[s(\xi, \eta), z(\xi, \eta)] = \sum_{e=1}^{n_{\Gamma}} \sum_{i=0}^{N} \sum_{j=0}^{N} (u_{i,s}^{\text{eij}}, u_{t,\phi}^{\text{eij}}, u_{t,z}^{\text{eij}}) I_i^N(\xi) h_j^N(\eta) + \sum_{e=n_{\Gamma}+1}^{n_e} \sum_{i=0}^{N} \sum_{j=0}^{N} (u_{t,s}^{\text{eij}}, u_{t,\phi}^{\text{eij}}, u_{t,z}^{\text{eij}}) h_i^N(\xi) h_j^N(\eta). \quad (B1)$$

The $(u_{i,s}^{\text{eij}}, u_{t,\phi}^{\text{eij}}, u_{t,z}^{\text{eij}})$ are the nodal velocities at the collocation points in the $e$th element, and $\Omega_{\Gamma}(\Omega_{\partial\Gamma})$ refers to the collection of elements which are (not) adjacent to the axis $\Gamma$.

The triplet $(e, i, j)$ corresponds to a local elemental ordering of the nodes. It is useful to define a global numbering of the velocity nodes $\text{num}_v$:

$$\text{num}_v: \{1, 2, 3\} \times \{1, \dots, n_e\} \times \{0, \dots, N\}^2 \to \{1, \dots, N_v\} \\ (\alpha, e, i, j) \mapsto \text{num}_v(\alpha, e, i, j), \quad (B2)$$

in which $\alpha$ refers to a component of the velocity field and $N_v$ is the total number of degrees of freedom for the velocity field ($N_v = N_{v,s} + N_{v,\phi} + N_{v,z}$). Note that $\text{num}_v$ is a surjection as some nodes can belong to more than one element. $N_v$ can be estimated only when the exact topology of the mesh is known.
The discrete pressure is expanded according to

$$p_{t,h}[s(\xi, \eta), z(\xi, \eta)] = \sum_{e=1}^{n_{\Gamma}} \sum_{i=1}^{N-1} \sum_{j=1}^{N-1} p_t^{\text{eij}} l_i^{N-2}(\xi) h_j^{N-2}(\eta) + \sum_{e=n_{\Gamma}+1}^{n_e} \sum_{i=1}^{N-1} \sum_{j=1}^{N-1} p_t^{\text{eij}} h_i^{N-2}(\xi) h_j^{N-2}(\eta), \quad (B3)$$

and we can define equivalently a global numbering of the pressure nodes num<sub>p</sub>x:

$$\text{num}_p : \{1, \dots, n_e\} \times \{1, \dots, N-1\}^2 \to \{1, \dots, N_p\} \\ (e, i, j) \mapsto \text{num}_p(e, i, j), \quad (B4)$$

As pressure is defined element-wise, num<sub>p</sub> is a bijection. The total number of pressure degrees of freedom is N<sub>p</sub> = n<sub>e</sub>(N - 1)<sup>2</sup>.

Following a Galerkin procedure (Reddy 1993), the trial spaces X<sub>o,h</sub> and Y<sub>h</sub> are built with the nodal shape functions associated with the velocity and pressure degrees of freedom, respectively.

Velocity trial functions Ψ<sub>I</sub>, 1 ≤ I ≤ N<sub>v</sub> can be conveniently chosen so that, when one computes the integrals involved in problem (13), one obtains a set of three scalar equations, one for each component of the original momentum eq. (1a). In other words, we set

$$\Psi_I = \begin{cases} (\Psi_{I_s}, 0, 0), & I_s = I - \Delta N_s \quad \text{if } 1 \le I \le N_{v,s} \\ (0, \Psi_{I_\phi}, 0), & I_\phi = I - \Delta N_\phi, \quad \text{if } N_{v,s} + 1 \le I \le N_{v,s} + N_{v,\phi}, \\ (0, 0, \Psi_{I_z}), & I_z = I - \Delta N_z, \quad \text{if } N_{v,s} + N_{v,\phi} + 1 \le I \le N_v, \end{cases} \quad (B5)$$

where ΔN<sub>s</sub> = 0, ΔN<sub>φ</sub> = N<sub>v,s</sub> and ΔN<sub>z</sub> = N<sub>v,s</sub> + N<sub>v,φ</sub>.

For any α ∈ {1, 2, 3}, the scalar function Ψ<sub>Iα</sub> is defined by its restrictions on the collection of the subdomains F<sub>v</sub>(I<sub>α</sub>) which contain the I<sub>α</sub>th nodal value for component α. Let us define

$$\forall I_{\alpha} \in \{1, \dots, N_{v,\alpha}\}, \\ F_v(I_{\alpha}) = \{e \in \{1, \dots, n_e\}, \exists!(i_{ae}, j_{ae}) \in \{0, N\}^2, \text{num}_v(\alpha, e, i_{ae}, j_{ae}) = I_{\alpha} + \Delta N_{\alpha}\}. \quad (B6)$$

The definition of the velocity trial function can now be properly written in terms of its elemental restrictions

$$\forall e \in F_v(I_\alpha), \Psi_{I_\alpha} |_{\Omega_e} = \begin{cases} l_{i_{ae}}^N(\xi) h_{j_{ae}}^N(\eta) & \text{if } \Omega_e \in \Omega_\Gamma, \\ h_{i_{ae}}^N(\xi) h_{j_{ae}}^N(\eta) & \text{if } \Omega_e \in \Omega_\phi. \end{cases} \quad (B7)$$

The pressure test functions are on the contrary defined element-wise. We can define similarly F<sub>p</sub>(I) as

$$\forall I \in \{1, \dots, N_p\}, \\ F_p(I) = \{e \in \{1, \dots, n_e\}, \exists!(i_e, j_e) \in \{1, N-1\}^2, \text{num}_p(e, i_e, j_e) = I\}. \quad (B8)$$

As num<sub>p</sub> is a bijection, F<sub>p</sub>(I) is a singlet. Again, a pressure trial function Φ<sub>I</sub> is defined by its restriction on F<sub>p</sub>(I)

$$\forall I \in \{1, \dots, N_p\}, \Phi_I(\xi, \eta) = \begin{cases} l_{F_p(I)}^{N-2}(\xi) h_{j_{F_p(I)}}^{N-2}(\eta) & \text{if } \Omega_{F_p(I)} \in \Omega_\Gamma, \\ h_{F_p(I)}^{N-2}(\xi) h_{j_{F_p(I)}}^{N-2}(\eta) & \text{if } \Omega_{F_p(I)} \in \Omega_\phi. \end{cases} \quad (B9)$$

Having defined both the shape and trial spaces for velocity and pressure, we can use the quadrature rules (26) and (28) to compute the integrals involved in the variational formulation (13). The semi-discrete problem then writes: at any time t ∈ [0, T], find (u<sub>t</sub>, p<sub>t</sub>) solution of

$$M \partial_t u_t + C u_t + E K u_t - D^T p_t = M f_t, \quad (B10a)$$

$$-Du_t = 0. \quad (B10b)$$

M is the N<sub>v</sub> × N<sub>v</sub> mass matrix which can be expressed as

$$M = \begin{bmatrix} M_s & 0 & 0 \\ 0 & M_\phi & 0 \\ 0 & 0 & M_z \end{bmatrix}, \quad (B11)$$

where, for each component, the N<sub>v,α</sub> × N<sub>v,α</sub> mass matrix M<sub>α</sub> is defined as

$$[M_{\alpha}]_{I_{\alpha} J_{\alpha}} = (\Psi_{I_{\alpha}}, \Psi_{J_{\alpha}})_1, \quad (I_{\alpha}, J_{\alpha}) \in \{1, \dots, N_{v,\alpha}\}^2. \quad (B12)$$

As the basis functions are defined by the Lagrange interpolants upon the quadrature points, the mass matrix has the remarkable property of being diagonal.

$$[M_{\alpha}]_{I_{\alpha} J_{\alpha}} = \delta_{I_{\alpha} J_{\alpha}} \sum_{e \in F_v(I_{\alpha})} \begin{cases} \sigma_{i_{ae}} \rho_{j_{ae}} \frac{s(\zeta_{i_{ae}}^N, \xi_{j_{ae}}^N)}{1 + \zeta_{i_{ae}}^N} |\mathcal{J}^e| (\zeta_{i_{ae}}^N, \xi_{j_{ae}}^N) & \text{if } \Omega_e \in \Omega_\Gamma, \\ s(\xi_{i_{ae}}^N, \xi_{j_{ae}}^N) \rho_{i_{ae}} \rho_{j_{ae}} |\mathcal{J}^e| (\xi_{i_{ae}}^N, \xi_{j_{ae}}^N) & \text{if } \Omega_e \in \Omega_\phi, \end{cases} \quad (B13)$$
where $\delta$ is the Kronecker symbol.

C is the $N_v \times N_v$ Coriolis antisymmetric matrix,

$$ C = \begin{bmatrix} 0 & -2M_s & 0 \\ 2M_\phi & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix}. \qquad (B14) $$

K is the $N_v \times N_v$ stiffness matrix

$$ K = \begin{bmatrix} K_s + W_s & 0 & 0 \\ 0 & K_\phi + W_\phi & 0 \\ 0 & 0 & K_z \end{bmatrix}, \qquad (B15) $$

in which we have used the following notation:

$$ [K_\alpha]_{I_\alpha J_\alpha} = a_0(\Psi_{I_\alpha}, \Psi_{J_\alpha}), \quad (I_\alpha, J_\alpha) \in \{1, \dots, N_{v,\alpha}\}^2, \qquad (B16) $$

and

$$ [W_\alpha]_{I_\alpha J_\alpha} = \left( \frac{\Psi_{I_\alpha}}{s}, \frac{\Psi_{J_\alpha}}{s} \right)_1, \quad (I_\alpha, J_\alpha) \in \{1, \dots, N_{v,\alpha}\}^2. \qquad (B17) $$

We are using high-order polynomials (typically $N \ge 6$), and the long-range interactions between nodes makes $K_\alpha$ dense. We have also introduced the $N_v \times N_p$ rectangular gradient matrix $D^T$ given by

$$ D^T = \begin{bmatrix} D_s^T \\ 0 \\ D_z^T \end{bmatrix}. \qquad (B18) $$

The two non-zero components of this matrix are

$$ [D_s^T]_{I_s J} = \left( \partial_s \Psi_{I_s} + \frac{\Psi_{I_s}}{s}, \Phi_J \right)_1, \quad (I_s, J) \in \{1, \dots, N_{v,s}\} \times \{1, \dots, N_p\}, \qquad (B19) $$

and

$$ [D_z^T]_{I_z J} = (\partial_z \Psi_{I_z}, \Phi_J)_1, \quad (I_z, J) \in \{1, \dots, N_{v,z}\} \times \{1, \dots, N_p\}. \qquad (B20) $$

## APPENDIX C: A MULTILEVEL ELLIPTIC SOLVER BASED UPON AN OVERLAPPING SCHWARZ METHOD

The pressure increment at each time step is computed by inverting

$$ E = \frac{\Delta t}{3/2} DM^{-1}D^T. \qquad (C1) $$

<figure><img src="image_12.png" /><figcaption>Figure C1. Typical iteration count for the calculation of the pressure increment $\delta p$ in the case of the Proudman-Stewartson problem. The poor conditioning of the pseudo-Laplacian to invert leads to a very slow convergence when no preconditioner is used in the conjugate gradient algorithm (circles). The local component of the preconditioner removes efficiently the high-frequency content of the residual field (squares). Adding its coarse component removes the large-scale components of the error and leads to an extra factor of two reduction of the iteration counts (triangles).</figcaption></figure>
This matrix is symmetric positive definite, and is inverted by means of a preconditioned conjugate gradient algorithm. The preconditioner $P^{-1}$ is an additive Schwarz preconditioner of the form (Fischer 1997; Fischer et al. 2000)

$$P^{-1} = R_0^T E_0^{-1} R_0 + \sum_{e=1}^{n_e} R_e^T E_e^{-1} R_e. \qquad (C2)$$

It is the sum of a global coarse grid operator (subscript 0) and local subdomains operators (subscript e). $R_{0,e}$ and $R_{0,e}^T$ denote restriction and extension operators. The method has a natural parallel aspect in that the subdomains problems can be solved independently. It is based upon the same ideas as the more classical multigrid approach. The local Laplacian operators $E_e$ are defined over overlapping regions centred on each spectral element and aim at removing the high-frequency components of the residual field. The coarse grid operator, $E_0$, is constructed as the linear finite-element Laplacian derived from triangles which vertices are coincident with the vertices of the spectral element. It aims at removing the large scale component of the residual field, hence its fundamental global character, which theoretically makes the iteration count independent of $n_e$ (e.g. Smith et al. 1996). Both local and coarse problems are small enough so that the $E_e$ and $E_0$ are factorized and inverted using standard linear algebra libraries.

The efficiency of the preconditioner is illustrated in Fig. C1, for a typical calculation of the pressure increment. With respect to a non-preconditioned case, the iteration count is decreased roughly by a factor of ten. The coarse grid solver contributes in itself to a decrease by a factor of two, and has a very modest computational cost.