00footnotetext: © Crown copyright 2020. Reproduced with the permission of the
Controller of HMSO.

# Hybridised multigrid preconditioners for a compatible finite element
dynamical core

Jack D. Betteridge  Imperial College, London, United Kingdom  Colin J. Cotter
Imperial College, London, United Kingdom  Thomas H. Gibson  National Center
for Supercomputing Applications, Urbana, IL, USA  Advanced Micro Devices,
Inc., Austin, TX, United States  Matthew J. Griffith  University of Bath, UK
European Centre for Medium Range Weather Forecasts, Reading, UK  Thomas Melvin
Met Office, Exeter, United Kingdom  Eike H. Müller  University of Bath, UK
Email<EMAIL_ADDRESS>
###### Abstract

Compatible finite element discretisations for the atmospheric equations of
motion have recently attracted considerable interest. Semi-implicit
timestepping methods require the repeated solution of a large saddle-point
system of linear equations. Preconditioning this system is challenging since
the velocity mass matrix is non-diagonal, leading to a dense Schur complement.
Hybridisable discretisations overcome this issue: weakly enforcing continuity
of the velocity field with Lagrange multipliers leads to a sparse system of
equations, which has a similar structure to the pressure Schur complement in
traditional approaches. We describe how the hybridised sparse system can be
preconditioned with a non-nested two-level preconditioner. To solve the coarse
system, we use the multigrid pressure solver that is employed in the
approximate Schur complement method previously proposed by the some of the
authors. Our approach significantly reduces the number of solver iterations.
The method shows excellent performance and scales to large numbers of cores in
the Met Office next-generation climate- and weather prediction model LFRic.
keywords: Numerical methods and NWP, Dynamics, Atmosphere, Multigrid, Parallel
Scalability, Linear Solvers, Finite Element Discretisation, Hybridisation

##  1 Introduction

Over the past decade, finite element discretisations of the governing
equations for atmospheric motion have become increasingly popular. A majority
of the focus has been centred on high-order continuous Galerkin (CG), spectral
element, and discontinuous Galerkin (DG) methods (for an overview, see [1, 2,
3, 4, 5, 6, 7]). More recently, a theory of mixed-finite element
discretisations, known as mimetic or compatible finite element methods [8, 9,
10, 11, 12, 13], has been developed within the context of simulating
geophysical fluid dynamics. Compatible finite element elements are constructed
from discrete spaces that satisfy an $L^{2}$ de-Rham co-homology in which
differential operators surjectively map from one space to another [14, 15,
16]. Discretisations arising from this framework have a long history in both
numerical analysis and applications ranging from structural mechanics to
porous media flows, see [17, 18] for a summary of contributions. The embedded
property of these discrete spaces defines a finite element extension of the
Arakawa C-grid staggered finite difference approaches [19], while avoiding the
need for discretising on structured, orthogonal meshes (such as traditional
latitude-longitude meshes). The realization made by Cotter and Shipton in [8]
that compatible finite element methods respect essential geostrophic balance
relations for steady-state solutions of the rotating shallow water equations,
combined with other properties making them analogous to the Arakawa C-grid,
has galvanized an effort to develop similar techniques for numerical weather
prediction. This framework of compatible finite elements is a core area of
research in the UK Met Office’s effort to redesign an operational dynamical
core, known as the Gung-Ho project. A key result of this project has been the
development of the LFRic model [20, 21], which solves the compressible Euler
equations for a perfect gas in a rotating frame. The model uses a compatible
finite element discretisation on a hexahedral grid based on tensor product
spaces derived from the classical Raviart-Thomas (RT) method [22].

###  1.1 Semi-implicit time discretisation of the compressible Euler equations

Atmospheric flows are characterised by a separation of scales: although fast
acoustic and inertia-gravity waves carry very little energy, they need to be
represented in the model since they drive slower large-scale balanced flow
through non-linear interactions. As a result, the timestep size in a model
based on fully explicit timestepping would be prohibitively small due to the
CFL condition imposed by the fast acoustic waves. To overcome this issue, the
LFRic model employs a semi-implicit time discretisation which treats the fast
waves implicitly in all spatial directions. This allows simulations with a
significantly larger timestep size, which can potentially reduce the overall
model runtime.

The drawback of semi-implicit methods is that they require the solution of a
non-linear system in every model timestep. Since the implicit system is only
mildly non-linear, it can be solved efficiently with a quasi-Newton-type
iteration, which in turn requires the repeated solution of a linear system.
This linear solve is one of the computational bottlenecks of the LFRic model.
While compatible finite element discretisations have several useful properties
for geophysical flow applications, the efficient implementation of solver
algorithms for the innermost linear system poses a significant challenge. This
is primarily due to the fact that, unlike for standard finite difference
discretisations, the linear system takes the form of an indefinite mixed
problem with saddle-point structure. To deliver forecasts under strict
operational time constraints, solvers must efficiently scale to many cores on
distributed-memory clusters while being algorithmically robust.

###  1.2 Hybridisation and related work on preconditioners

There are a number of possible approaches for handling the linear saddle-point
system arising from compatible finite element discretisations. This includes
H(div)-multigrid [23], requiring complex overlapping Schwarz smoothers;
auxiliary space multigrid [24]; and global Schur-complement preconditioners,
requiring inverting both a velocity mass matrix and an elliptic Schur-
complement operator. So far, a pragmatic solution has been employed in LFRic:
the saddle-point system is preconditioned with an approximate Schur-complement
method which is obtained by lumping the velocity mass matrix [25, 26]. The
resulting elliptic Schur-complement system in the piecewise constant pressure
space is solved with a bespoke multigrid algorithm that exploits the tensor-
product structure of the function spaces, following [27]. The drawback of this
approach is that it requires the repeated, expensive application of the
saddle-point operator for all four prognostic variables in a preconditioned
Krylov-iteration.

To address this issue, we focus on a technique known as hybridisation [28, 29]
in this paper. Hybridisation utilises element-wise static condensation to
obtain a reduced problem on the mesh interfaces. Once solved, recovery of the
prognostic variables is achieved by solving purely local linear systems.
Although the reduced, sparse system is smaller than the original saddle-point
system for the prognostic variables, it has to be solved efficiently. Non-
nested algorithms [30] have been proposed for the solution of linear systems
that arise from the hybridisation of mimetic finite element discretisations.
The authors of [31] describe a non-nested multigrid method for the solution of
the hybridisable mixed formulation of a second order elliptic boundary value
problem, while [32] use a similar approach for the DG discetisation of the
same equation. Both papers employ a conforming coarse level space based on
piecewise linear elements; this space allows the construction of a standard
multigrid method with $h$-coarsening. The hybridised method has also been
employed to solve the linear system that arises in an IMEX-DG discretisation
of the shallow water equations in [33]. While the authors of [33] use a direct
solver for the resulting elliptic system on the mesh interfaces, a non-nested
multigrid algorithm similar to [32] is used for the IMEX-DG discretisation of
the shallow water equations in [34].

###  1.3 Main contributions

In this paper, we describe how a hybridisable discretisation of the
compressible Euler equations in LFRic allows the efficient solution of the
saddle-point system for the prognostic variables if it is combined with a non-
nested multigrid preconditioner. Our work is motivated by two scientific
questions:

  1. 1.

Can hybridisation be used to solve the saddle-point system in LFRic more
efficiently?

  2. 2.

What performance gains can be achieved by using bespoke non-nested multigrid
preconditioners?

The new contributions of our work are:

  1. i.

We decribe how the exact elimination of density and potential temperature
allows the construction of a reduced saddle-point system for Exner pressure
and momentum; solving this system with the approximate Schur-complement
approach with a pressure multigrid preconditioner improves on the results in
[25] and forms the basis for our hybridisable discretisation.

  2. ii.

We derive a hybridisable discretisation of the compressible Euler equations
and use static condensation to construct an elliptic system for the unknowns
on the mesh skeleton.

  3. iii.

We introduce a non-nested multigrid algorithm for the solution of this
elliptic system; since the unknowns on the mesh skeleton approximate pressure,
we use the pressure multigrid algorithm from [25] to solve the piecewise
constant pressure equation on the coarse level of the multigrid hierarchy.

  4. iv.

We compare the performance of different solver setups to demonstrate the
impact of hybridisation and the use of a multigrid preconditioner, thus
providing empirical evidence to answer our two central scientific questions
above.

  5. v.

We demonstrate the excellent parallel scalability of our implementation on up
to 13824 cores.

The choice of a piecewise constant coarse level space in the non-nested
multigrid preconditioner is mainly motivated by practical considerations,
since it allows the re-use of the efficient pressure multigrid algorithm in
[25]. Although a rigorous theoretical justification of this approach is beyond
the scope of this paper, our numerical results show that it works very well in
practice. While for the present, lowest order discretisation employed in
LFRic, our hybridised solver is not able to beat the optimised approximate
Schur-complement approach in [25], we will argue in Section 4 that this
picture is likely to change for higher-order discretisations. This is due to
the dramatic increase in condition number and problem size for the resulting
mixed operators. In situations like this, it has been show that static
condensation approachs can provide significant benefit [35]. We therefore
believe that the hybridised solver introduced in this paper has significant
potential to accelerate semi-implicit compatible finite element models for
atmospheric fluid dynamics.

#### Structure.

The rest of this paper is organised as follows: in Section 2 we describe the
mixed finite element discretisation of the compressible Euler equations that
is used in LFRic and show how potential temperature and density can be
eliminated to obtain a smaller linear system for Exner pressure and velocity
only. We also explain how hybridisation allows the construction of a sparse
Schur complement system. After reviewing the approximate Schur-complement
solver and pressure multigrid algorithm from [25], a new non-nested multigrid
algorithm for solving this system is introduced in this section. Following a
description of our implementation and the model setup, we compare the
performance of different solver algorithms for a baroclinic test case in
Section 3, which also contains results from massively parallel scaling runs.
We conclude and outline avenues for further work in Section 4. Further
technical details on the LFRic implementation and parameter tuning can be
found in the appendices.

##  2 Methods

###  2.1 Model equations

To obtain an update for the non-linear outer iteration in a semi-implicit
timestepping scheme we need to compute increments in Exner pressure
$\Pi^{\prime}$, potential temperature $\theta^{\prime}$, density
$\rho^{\prime}$ and velocity $\bm{u}^{\prime}$. This can be achieved by
linearising the compressible Euler equations around a suitable reference state
and solving the following four equations:

$\displaystyle\left[\bm{u}^{\prime}+\tau_{u}\Delta
t\left(\mu\widehat{\bm{z}}(\widehat{\bm{z}}\cdot\bm{u}^{\prime})+2\bm{\Omega}\times\bm{u}^{\prime}\right)\right]+\tau_{u}\Delta
tc_{p}\left(\theta^{\prime}\nabla\Pi^{*}+\theta^{*}\nabla\Pi^{\prime}\right)$
$\displaystyle=$ $\displaystyle\bm{r}_{u},$ (1a)
$\displaystyle\rho^{\prime}+\tau_{\rho}\Delta
t\nabla\cdot\left(\rho^{*}\bm{u}^{\prime}\right)$ $\displaystyle=$
$\displaystyle r_{\rho},$ (1b)
$\displaystyle\theta^{\prime}+\tau_{\theta}\Delta
t\bm{u}^{\prime}\cdot\nabla\theta^{*}$ $\displaystyle=$ $\displaystyle
r_{\theta},$ (1c)
$\displaystyle\frac{1-\kappa}{\kappa}\frac{\Pi^{\prime}}{\Pi^{*}}-\frac{\rho^{\prime}}{\rho^{*}}-\frac{\theta^{\prime}}{\theta^{*}}$
$\displaystyle=$ $\displaystyle r_{\Pi}.$ (1d)

Here $\Delta t$ is the timestep size, $\tau_{u},\tau_{\rho},\tau_{\theta}$ are
relaxation parameters, the unit normal in the vertical direction is denoted by
$\widehat{\bm{z}}$, and we have $\mu=1+\tau_{u}\tau_{\theta}\Delta t^{2}N^{2}$
with the Brunt-Väisälä frequency $N^{2}$. The right-hand sides $\bm{r}_{u}$,
$r_{\rho}$, $r_{\theta}$, and $r_{\Pi}$ are residuals obtained from the
linearisation process. In LFRic the reference states $\theta^{*}$, $\rho^{*}$
and $\Pi^{*}$ are the fields at the previous timestep. The expressions on the
right-hand side depend on the current non-linear iterate, but their exact form
is not relevant for the following derivation. The equations have to be solved
in a thin spherical shell $\Omega$, with suitable boundary conditions at the
top and bottom of the atmosphere. For a further discussion of the equations we
refer the reader to [36] and [25, Section 3.3]. It should be noted that in
[25], Eqs. (1a) - (1d) are discretised to obtain a $4\times 4$ block matrix
system, which is solved with a suitable preconditioned Krylov subspace method.
Since each application of the $4\times 4$ block matrix is expensive, we pursue
a different approach here: by eliminating $\theta^{\prime}$ and
$\rho^{\prime}$ first, we obtain a smaller $2\times 2$ block system for
velocity $\bm{u}^{\prime}$ and Exner pressure $\Pi^{\prime}$ only. This system
is then solved with a preconditioned Krylov subspace method or the hybridised
approach described in Section 2.4, followed by the reconstruction of
$\theta^{\prime}$ and $\rho^{\prime}$ from $\bm{u}^{\prime}$ and
$\Pi^{\prime}$. Since the elimination and recovery of $\theta^{\prime}$ and
$\rho^{\prime}$ is only necessary once per linear solve, this increases
efficiency.

#### Analytic elimination of potential temperature.

Eq. (1c) can be used to eliminate $\theta^{\prime}$ from Eqs. (1a), (1b) and
(1d), leading to the following system of three equations:

$\displaystyle\left[\bm{u}^{\prime}+\tau_{u}\Delta
t\left(\mu\widehat{\bm{z}}(\widehat{\bm{z}}\cdot\bm{u}^{\prime})+2\bm{\Omega}\times\bm{u}^{\prime}\right)-\tau_{u}\tau_{\theta}\Delta
t^{2}c_{p}\left(\bm{u}^{\prime}\cdot\nabla\theta^{*}\right)\nabla\Pi^{*}\right]+\tau_{u}\Delta
tc_{p}\theta^{*}\nabla\Pi^{\prime}$ $\displaystyle=$
$\displaystyle\overline{\bm{r}}_{u}\equiv\bm{r}_{u}-\tau_{u}\Delta
tc_{p}\nabla\Pi^{*}r_{\theta},$ (2a)
$\displaystyle\rho^{\prime}+\tau_{\rho}\Delta
t\nabla\cdot\left(\rho^{*}\bm{u}^{\prime}\right)$ $\displaystyle=$
$\displaystyle r_{\rho},$ (2b)
$\displaystyle\frac{1-\kappa}{\kappa}\frac{\Pi^{\prime}}{\Pi^{*}}-\frac{\rho^{\prime}}{\rho^{*}}+\tau_{\theta}\Delta
t\frac{\bm{u}^{\prime}\cdot\nabla\theta^{*}}{\theta^{*}}$ $\displaystyle=$
$\displaystyle\overline{r}_{\Pi}\equiv r_{\Pi}+\frac{r_{\theta}}{\theta^{*}}.$
(2c)

###  2.2 Mixed finite element discretisation

To discretise Eqs. (2a) - (2c), we proceed as in [25, Section 3.2] and
introduce finite element spaces on a hexahedral grid $\Omega_{h}$ covering the
domain $\Omega$: $\mathbb{W}_{2}$ is the Raviart-Thomas [37] space $RT_{p}$
and $\mathbb{W}_{3}=Q_{p}^{\text{DG}}$ is the scalar discontinuous Galerkin
space, where in both cases $p$ is the order of the discretisation. The
velocity space
$\mathbb{W}_{2}=\mathbb{W}_{2}^{\text{h}}\oplus\mathbb{W}_{2}^{\text{z}}$ can
be written as the direct sum of a space $\mathbb{W}_{2}^{h}$, which only
contains vectors pointing in the horizontal, and a space
$\mathbb{W}_{2}^{\text{z}}$ with purely vertical fields (see [25, Figure 1]).
Recall further that the space $\mathbb{W}_{2}^{\text{z}}$ is continuous in the
vertical direction, but discontinuous in the horizontal direction, while
$\mathbb{W}_{2}^{\text{h}}$ is horizontally continuous and vertically
discontinuous. To represent the reference profile for potential temperature, a
third space $\mathbb{W}_{\theta}$ is introduced. $\mathbb{W}_{\theta}$ is the
scalar-valued equivalent of $\mathbb{W}_{2}^{\text{z}}$ and has the same
continuity, which is analogous to using a Charney-Phillips staggered vertical
discretisation. A further discussion of the function spaces can be found in
[20].

We now assume that $\bm{u}^{\prime}\in\mathbb{W}_{2}$,
$\Pi^{\prime},\rho^{\prime},\Pi^{*},\rho^{*}\in\mathbb{W}_{3}$ and
$\theta^{*}\in\mathbb{W}_{\theta}$. Multiplying Eqs. (2a) - (2c) by suitable
test functions $\bm{w}\in\mathbb{W}_{2}$ and $\phi,\psi\in\mathbb{W}_{3}$ and
integrating over $\Omega_{h}$ leads to the following weak form of the problem:

$\displaystyle\breve{m}_{2}(\bm{w},\bm{u}^{\prime})-q_{22}(\bm{w},\bm{u}^{\prime})+g(\bm{w},\Pi^{\prime})$
$\displaystyle=r_{u}(\bm{w}),$ (3) $\displaystyle
d(\phi,\bm{u}^{\prime})+m_{3}(\phi,\rho^{\prime})$
$\displaystyle=r_{\rho}(\phi),\qquad\text{for all $\bm{w}\in\mathbb{W}_{2}$
and all $\phi,\psi\in\mathbb{W}_{3}$.}$ $\displaystyle
q_{32}(\psi,\bm{u}^{\prime})-m_{3}^{\rho}(\psi,\rho^{\prime})+m_{3}^{\Pi}(\psi,\Pi^{\prime})$
$\displaystyle=r_{\Pi}(\psi),$

Writing $\langle\cdot\rangle_{\Omega_{h}}$ for integrals over the mesh
$\Omega_{h}$ and $\langle\cdot\rangle_{\mathcal{E}_{h}}$ for integrals over
the mesh skeleton $\mathcal{E}_{h}$, i.e.

$\displaystyle\langle f\rangle_{\Omega_{h}}$
$\displaystyle=\sum_{K\in\Omega_{h}}\int_{K}f\;dx,$ $\displaystyle\langle
g\rangle_{\mathcal{E}_{h}}$ $\displaystyle=\sum_{\partial
K\in\mathcal{E}_{h}}\int_{\partial K}g\;dS,$

the eight bilinear forms on the left-hand side of Eq. (3) are

$\displaystyle\breve{m}_{2}(\bm{w},\bm{u}^{\prime})$
$\displaystyle=\left\langle\bm{w}\cdot\left(\bm{u}^{\prime}+\tau_{u}\Delta
t\left(\mu\widehat{\bm{z}}(\widehat{\bm{z}}\cdot\bm{u}^{\prime})+2\bm{\Omega}\times\bm{u}^{\prime}\right)\right)\right\rangle_{\Omega_{h}},$
(4) $\displaystyle q_{22}(\bm{w},\bm{u}^{\prime})$
$\displaystyle=\tau_{u}\tau_{\theta}\Delta t^{2}\left\langle
c_{p}(\widehat{\bm{z}}\cdot\bm{w})\partial_{z}\widehat{\Pi}^{*}(\widehat{\bm{z}}\cdot\bm{u}^{\prime})\partial_{z}\theta^{*}\right\rangle_{\Omega_{h}},$
$\displaystyle g(\bm{w},\Pi^{\prime})$ $\displaystyle=\tau_{u}\Delta
tc_{p}\left(\left\langle[\\![\theta^{*}\bm{w}]\\!]\\{\\!\\!\\{\Pi^{\prime}\\}\\!\\!\\}\right\rangle_{\mathcal{E}_{h}}-\left\langle\nabla\cdot(\theta^{*}\bm{w})\Pi^{\prime}\right\rangle_{\Omega_{h}}\right),$
$\displaystyle d(\phi,\bm{u}^{\prime})$ $\displaystyle=\tau_{\rho}\Delta
t\left\langle\phi\nabla\cdot(\rho^{*}\bm{u}^{\prime})\right\rangle_{\Omega_{h}},$
$\displaystyle m_{3}(\phi,\rho^{\prime})$
$\displaystyle=\left\langle\phi\rho^{\prime}\right\rangle_{\Omega_{h}},$
$\displaystyle q_{32}(\psi,\bm{u}^{\prime})$
$\displaystyle=\tau_{\theta}\Delta
t\left\langle\psi\frac{\partial_{z}\theta^{*}}{\theta^{*}}(\widehat{\bm{z}}\cdot\bm{u}^{\prime})\right\rangle_{\Omega_{h}},$
$\displaystyle m_{3}^{\rho}(\psi,\rho^{\prime})$
$\displaystyle=\left\langle\frac{\psi\rho^{\prime}}{\rho^{*}}\right\rangle_{\Omega_{h}},$
$\displaystyle m_{3}^{\Pi}(\psi,\Pi^{\prime})$
$\displaystyle=\frac{1-\kappa}{\kappa}\left\langle\frac{\psi\Pi^{\prime}}{\Pi^{*}}\right\rangle_{\Omega_{h}}.$

Since it is important for the discussion of hybridisation in Section 2.4, we
explain the derivation of the weak form $g(\cdot,\cdot)$ in more detail. For
this, consider the term $\theta^{*}\nabla\Pi^{\prime}$ in Eq. (1a). The
function space $\mathbb{W}_{3}$ is discontinuous and the derivative of
$\Pi^{\prime}$ can not be evaluated directly. Instead, we integrate the term
by parts after it has been multiplied by a test function
$\bm{w}\in\mathbb{W}_{2}$. Considering a single cell $K$ of the mesh, this
leads to two integrals:

$\int_{K}\theta^{*}\bm{w}\cdot\nabla\Pi^{\prime}\;dx=-\int_{K}\nabla\cdot\left(\theta^{*}\bm{w}\right)\Pi^{\prime}\;dx+\int_{\partial
K}\left(\theta^{*}\bm{w}\cdot\bm{n}\right)\Pi^{\prime}|_{\partial K}\;dS,$ (5)

where $\bm{n}$ is the unit outward normal vector at each point of the cell
surface $\partial K$. In the surface integral, we need to decide how to
evaluate the field $\Pi^{\prime}|_{\partial K}$, which is not well-defined
since $\mathbb{W}_{3}$ is discontinuous. The natural choice is to take the
average
$\\{\\!\\!\\{\Pi^{\prime}\\}\\!\\!\\}:=\frac{1}{2}\left(\Pi^{\prime}_{+}+\Pi^{\prime}_{-}\right)$
of the fields in the two cells $K^{+}$, $K^{-}$ that are adjacent to a
particular facet. Summing the right-hand side of Eq. (5) over all cells of the
mesh then leads to the expression in $g(\cdot,\cdot)$, namely

$-\left\langle\nabla\cdot\left(\theta^{*}\bm{w}\right)\Pi^{\prime}\right\rangle_{\Omega_{h}}+\left\langle[\\![\theta^{*}\bm{w}]\\!]\\{\\!\\!\\{\Pi^{\prime}\\}\\!\\!\\}\right\rangle_{\mathcal{E}_{h}}.$
(6)

In this expression
$[\\![\theta^{*}\bm{w}]\\!]\equiv\theta^{*}_{+}(\bm{w}_{+}\cdot\bm{n}_{+})+\theta^{*}_{-}(\bm{w}_{-}\cdot\bm{n}_{-})$
and the subscripts “$+$” and “$-$” are again used to identify quantities
defined on the two cells that are adjacent to a particular facet. Note that
although here $\bm{w}_{+}\cdot\bm{n}_{+}=-\bm{w}_{-}\cdot\bm{n}_{-}$ as
$\bm{n}_{+}=-\bm{n}_{-}$ and the normal components of the velocity field are
continuous across facets, this is no longer true in our hybridisable
formulation (see Section 2.4).

To obtain the bilinear forms $q_{22}(\cdot,\cdot)$ and $q_{32}(\cdot,\cdot)$,
the following approximations have been made: when evaluating $\nabla\Pi^{*}$
in Eq. (2a), the reference profile is projected into the function space
$\mathbb{W}_{\theta}$ to obtain $\widehat{\Pi}^{*}\in\mathbb{W}_{\theta}$.
Since $\mathbb{W}_{\theta}$ is only continuous in the vertical direction, the
horizontal components of the derivatives $\nabla\theta^{*}$ and
$\nabla\Pi^{*}$ in Eqs. (2a) and (2c) are dropped to obtain
$\widehat{\bm{z}}\partial_{z}\theta^{*}$ and
$\widehat{\bm{z}}\partial_{z}\widehat{\Pi}^{*}$. This is a sensible
approximation because the reference profiles vary predominantly in the
vertical direction.

The linear forms on the right-hand side of Eq. (3) are given by
$r_{u}(\bm{w})=\left\langle\bm{w}\cdot\overline{\bm{r}}_{u}\right\rangle_{\Omega_{h}}$,
$r_{\rho}(\phi)=\left\langle\phi r_{\rho}\right\rangle_{\Omega_{h}}$ and
$r_{\Pi}(\psi)=\left\langle\psi\overline{r}_{\Pi}\right\rangle_{\Omega_{h}}$.

In general, the integrals in Eq. (4) are approximated by $n$-point Gaussian
quadrature rules (with $n=p+3$ for an element of order $p$), in each
coordinate direction for the volume integrals
$\langle\cdot\rangle_{\Omega_{h}}$ and in two dimensions restricted to each
facet for the facet integrals $\langle\cdot\rangle_{\mathcal{E}_{h}}$. This
value of $n$ is chosen as it is generally sufficient to integrate the terms in
Eq. (4) exactly on an affine mesh. However, since it desirable for the
equation of state $r_{\Pi}$ in Eq. (2c) to hold exactly, the bilinear forms
$q_{32}$, $m_{3}^{\rho}$ and $m_{3}^{\Pi}$ in Eq. (4) are instead evaluated at
the nodal points of the $\mathbb{W}_{3}$ space, such that the diagnostic
relationship $r_{\Pi}$ between $\rho,\,\theta,\,\Pi$ holds exactly at these
nodal points. For $p=0$ elements this can be simply enforced by using a one-
point quadrature rule with the quadrature point located in the cell centre.
Since for $p=0$ the test functions $\psi\in\mathbb{W}_{3}$ are constant, this
is equivalent to enforcing Eq. (2c) in strong form at the quadrature points.

#### Matrix formulation.

Choosing suitable basis functions $\bm{v}_{j}(\bm{\chi})\in\mathbb{W}_{2}$,
$\sigma_{j}(\bm{\chi})\in\mathbb{W}_{3}$ that depend on the spatial coordinate
$\bm{\chi}\in\Omega_{h}$, the fields $\bm{u}^{\prime}$, $\rho^{\prime}$ and
$\Pi^{\prime}$ can be expanded as

$\displaystyle\bm{u}^{\prime}(\bm{\chi})$
$\displaystyle=\sum_{j}\widetilde{u}^{\prime}_{j}\bm{v}_{j}(\bm{\chi})\in\mathbb{W}_{2},$
$\displaystyle\rho^{\prime}(\bm{\chi})$
$\displaystyle=\sum_{j}\widetilde{\rho}^{\prime}_{j}\sigma_{j}(\bm{\chi})\in\mathbb{W}_{3},$
$\displaystyle\Pi^{\prime}(\bm{\chi})$
$\displaystyle=\sum_{j}\widetilde{\Pi}^{\prime}_{j}\sigma_{j}(\bm{\chi})\in\mathbb{W}_{3},$
(7)

where for each quantity $a^{\prime}$ the corresponding dof (degrees of
freedom)-vector is given by $\widetilde{a}^{\prime}=[a_{1},a_{2},\dots]$. This
allows expressing the weak form in Eq. (3) as a $3\times 3$ matrix equation
for the dof-vectors $\widetilde{u}^{\prime}$, $\widetilde{\rho}^{\prime}$ and
$\widetilde{\Pi}^{\prime}$:

$\begin{pmatrix}\breve{M}_{2}-Q_{22}&0&G\\\ D&M_{3}&0\\\
Q_{32}&-M_{3}^{\rho}&M_{3}^{\Pi}\end{pmatrix}\begin{pmatrix}\widetilde{u}^{\prime}\\\
\widetilde{\rho}^{\prime}\\\
\widetilde{\Pi}^{\prime}\end{pmatrix}=\begin{pmatrix}\mathcal{R}_{u}\\\
\mathcal{R}_{\rho}\\\ \mathcal{R}_{\Pi}\end{pmatrix}.$ (8)

The entries of the matrices $\breve{M}_{2}$, $Q_{22}$, $G$, $D$, $M_{3}$,
$Q_{32}$, $M_{3}^{\rho}$ and $M_{3}^{\Pi}$ are obtained by evaluating the
bilinear forms in Eq. (4) for the basis functions, for example
$(\breve{M}_{2})_{ij}=\breve{m}_{2}(\bm{v}_{i},\bm{v}_{j})$ and
$G_{ij}=g(\bm{v}_{i},\sigma_{j})$. Similarly, the components of the residual
(dual) vectors $\mathcal{R}_{u}$, $\mathcal{R}_{\rho}$ and $\mathcal{R}_{\Pi}$
on the right-hand side of Eq. (8) are obtained by evaluating the linear forms
$r_{u}$, $r_{\rho}$ and $r_{\Pi}$ for the basis functions, for example
$(\mathcal{R}_{u})_{i}=r_{u}(\bm{v}_{i})$.

#### Elimination of density.

Finally, the density $\widetilde{\rho}^{\prime}$ is eliminated algebraically
from Eq. (8) by using
$\widetilde{\rho^{\prime}}=M_{3}^{-1}(\mathcal{R}_{\rho}-D\widetilde{u^{\prime}})$
to obtain a $2\times 2$ matrix system for velocity $\widetilde{u}^{\prime}$
and Exner pressure $\widetilde{\Pi}^{\prime}$ only:

$\begin{pmatrix}\breve{M}_{2}-Q_{22}&G\\\
Q_{32}+M_{3}^{\rho}M_{3}^{-1}D&M_{3}^{\Pi}\end{pmatrix}\begin{pmatrix}\widetilde{u}^{\prime}\\\
\widetilde{\Pi}^{\prime}\end{pmatrix}=\begin{pmatrix}\mathcal{R}_{u}\\\
\overline{\mathcal{R}}_{\Pi}\end{pmatrix}\equiv\begin{pmatrix}\mathcal{R}_{u}\\\
\mathcal{R}_{\Pi}+M_{3}^{\rho}M_{3}^{-1}\mathcal{R}_{\rho}\end{pmatrix}$ (9)

In the next section, we review a suitable preconditioned Krylov subspace
method for solving Eq. (9) (namely the approximate Schur complement pressure
multigrid method in [25]) and then introduce an alternative solution method
based on hybridisation in Section 2.4.

###  2.3 Pressure multigrid for approximate Schur complement solver

The traditional approach for solving Eq. (9), which is pursued in the ENDGame
dynamical core [36], is to eliminate the velocity $\widetilde{u}^{\prime}$ and
then solve the resulting elliptic Schur complement equation for the Exner
pressure increment $\widetilde{\Pi}^{\prime}$. For the finite element
discretisation considered here, however, the matrix $\breve{M}_{2}-Q_{22}$ is
not diagonal. This makes this approach unfeasible since the resulting Schur
complement would be dense because it contains the inverse of
$\breve{M}_{2}-Q_{22}$. As in [25], we address this issue by constructing an
approximate Schur complement which is sparse, and then using this to
precondition a Krylov-subspace iteration over the mixed system in Eq. (9). To
construct the approximate Schur complement, we replace $\breve{M}_{2}-Q_{22}$
by a suitable lumped, diagonal approximation $\mathring{M}_{2}$, which results
in the following $2\times 2$ system:

$\begin{pmatrix}\mathring{M}_{2}&G\\\
Q_{32}+M_{3}^{\rho}M_{3}^{-1}D&M_{3}^{\Pi}\end{pmatrix}\begin{pmatrix}\widetilde{u}^{\prime}\\\
\widetilde{\Pi}^{\prime}\end{pmatrix}=\begin{pmatrix}\mathcal{R}_{u}\\\
\overline{\mathcal{R}}_{\Pi}\end{pmatrix}.$ (10)

Eq. (10) has a sparse Schur complement and can be solved as follows:

Step 1:

Compute the modified right-hand side

$\mathcal{B}=\overline{\mathcal{R}}_{\Pi}-\left(Q_{32}+M_{3}^{\rho}M_{3}^{-1}D\right)\mathring{M}_{2}^{-1}\mathcal{R}_{u}$
(11)

Step 2:

Solve the sparse elliptic system

$H\widetilde{\Pi}^{\prime}=\mathcal{B}\qquad\text{with}\quad
H=-\left(Q_{32}+M_{3}^{\rho}M_{3}^{-1}D\right)\mathring{M}_{2}^{-1}G+M_{3}^{\Pi}$
(12)

Step 3:

Recover the velocity using

$\widetilde{u}^{\prime}=\mathring{M}_{2}^{-1}(\mathcal{R}_{u}-G\widetilde{\Pi}^{\prime})$
(13)

To (approximately) solve the elliptic system in Eq. (12) we use a tensor-
product multigrid [27] approach. Due to the high aspect ratio of the grid
($\Delta z\ll\Delta x$), the Helmholtz operator is dominated by couplings
within each vertical column. A suitable operator $\widehat{H}_{z}\approx H$,
which captures these columnwise couplings, is constructed as in [25, Section
3.6] and used as a block-Jacobi smoother in the multigrid iteration:

$\widetilde{\Pi}^{\prime}\mathrel{\reflectbox{$\mapsto$}}\widetilde{\Pi}^{\prime}+\omega\widehat{H}_{z}^{-1}\left(\mathcal{B}-H\widetilde{\Pi}^{\prime}\right),$
(14)

where $\omega$ is a relaxation factor. Since the matrix $\widehat{H}_{z}$ is
block-diagonal, it can be inverted independently in each vertical column with
the Thomas algorithm [38]. We refer the reader to [25] for further details,
such as the construction of suitable prolongation and restriction operators.

###  2.4 Hybridisable formulation

The main drawback of the approximate Schur complement method in Section 2.3 is
that it still requires an iterative scheme to be applied to the mixed operator
in Eq. (9) which is expensive: at lowest order, on average four unknowns are
stored per grid cell to represent the pressure and velocity unknowns (one
pressure unknown per cell plus one shared velocity unknown for each of the six
facets). The mixed operator is stored and applied in the form of local
stiffness matrices and since there are seven unknowns per grid cell (one
pressure unknown and six velocity unknowns), each of these matrices has
$7^{2}=49$ unknowns, which makes the operator application expensive. The
hybridisable approach presented here seeks to address these issues. Recall
that the continuity of the velocity space prevents the construction of a
sparse Schur complement. We therefore introduce a broken function space
$\mathbb{W}_{2}^{b}$ for the velocity field such that in each grid cell
$\mathbb{W}_{2}^{b}$ is constructed from the same local basis elements as
$\mathbb{W}_{2}$, but — in contrast to $\mathbb{W}_{2}$ — does not necessarily
have a continuous normal component across neighbouring cells, see Figure 1.

Figure 1: Function space degrees of freedom (DOFs) used for hybridisable
formulation. The original velocity function space $\mathbb{W}_{2}$ (left) is
replaced by the broken function space $\mathbb{W}_{2}^{b}$ (centre). The trace
function space $\mathbb{W}_{2}^{\text{tr}}$ (right) can be decomposed into
$\mathbb{W}_{2}^{\text{tr,h}}$ and $\mathbb{W}_{2}^{\text{tr,z}}$; compare to
the decomposition of $\mathbb{W}_{2}$ in [25, Figure 1]. Arrows indicate DOFs
evaluating normal components of vectors, and points indicate evaluation of
scalars.

We expand the broken velocity field $\bm{u}^{\prime b}\in\mathbb{W}_{2}^{b}$
as in Eq. (7)

$\bm{u}^{\prime b}(\bm{\chi})=\sum_{j}\widetilde{u}^{\prime
b}_{j}\bm{v}^{b}_{j}(\bm{\chi})\in\mathbb{W}_{2}^{b},$ (15)

where now each basis function $\bm{v}^{b}_{j}(\bm{\chi})$ is non-zero only in
a single grid cell $K$. More specifically, each $\bm{v}^{b}_{j}(\bm{\chi})$ is
obtained by multiplying a Raviart Thomas basis function with the indicator
function $\chi_{K}(\bm{\chi})$ which satifies $\chi_{K}(\bm{\chi})=1$ for
$\bm{\chi}\in K$ and $\chi_{K}(\bm{\chi})=0$ otherwise. To weakly enforce the
continuity of the velocity fields, Lagrange multipliers are introduced. This
leads to an additional equation, which, if solved exactly, ensures that the
velocity field is projected into the subspace
$\mathbb{W}_{2}\subset\mathbb{W}_{2}^{b}$ such that $\mathbb{W}_{2}$ has
continuous normal components across facets.

Replacing $\mathbb{W}_{2}$ by $\mathbb{W}_{2}^{b}$ affects four of the
bilinear forms in Eq. (4). The forms $\breve{m}_{2}(\cdot,\cdot)$,
$d(\cdot,\cdot)$ and $q_{32}(\cdot,\cdot)$ are volume integrals, and
evaluating them for the basis functions $\bm{v}_{j}^{b}$ leads to three
matrices $\breve{M}_{2}^{b}$, $D^{b}$ and $Q_{32}^{b}$ with, for example
$D^{b}_{ij}=d(\sigma_{i},\bm{v}_{j}^{b})$. It is important to note that with
this choice of basis all three matrices are block-diagonal, coupling only
unknowns in a single grid cell. This also implies that $\breve{M}_{2}^{b}$,
$D^{b}$ and $Q_{32}^{b}$ can be inverted cell-wise. The bilinear form
$g(\cdot,\cdot)$ requires more thought. Instead of constructing it by using
the average $\\{\\!\\!\\{\Pi^{\prime}\\}\\!\\!\\}$ in the facet integral in
Eq. (6), we introduce a new field $\lambda\in\mathbb{W}_{2}^{\text{tr}}$ which
is an approximation of the Exner pressure increment in the continuous
equations. The trace space $\mathbb{W}_{2}^{\text{tr}}$ is a scalar valued
function space with support on the skeleton $\mathcal{E}_{h}$. It is piecewise
continuous, with each function being represented by a polynomial of degree $p$
on an individual facet. Consequently, the facet integral in Eq. (6) becomes
$\left\langle\right[\\![\theta^{*}\bm{w}]\\!]\lambda\rangle_{\mathcal{E}_{h}}$
and the bilinear form $g(\bm{w},\Pi^{\prime})$ in Eq. (4) is replaced by

$g^{b}(\bm{w}^{b},\Pi^{\prime})+k(\bm{w}^{b},\lambda)\qquad\text{for
$\bm{w}^{b}\in\mathbb{W}_{2}^{b}$, $\Pi\in\mathbb{W}_{3}$ and
$\lambda\in\mathbb{W}_{2}^{\text{tr}}$}$ (16)

with

$\displaystyle g^{b}(\bm{w}^{b},\Pi^{\prime})$ $\displaystyle=-\tau_{u}\Delta
tc_{p}\left\langle\nabla\cdot(\theta^{*}\bm{w}^{b})\Pi^{\prime}\right\rangle_{\Omega_{h}},$
(17) $\displaystyle k(\bm{w}^{b},\lambda)$ $\displaystyle=\tau_{u}\Delta
tc_{p}\left\langle[\\![\theta^{*}\bm{w}^{b}]\\!]\lambda\right\rangle_{\mathcal{E}_{h}}.$

Finally, continuity of the normal components of the velocity field is enforced
weakly through the equation

$c(\bm{u}^{\prime b},\nu)=\left\langle[\\![\bm{u}^{\prime
b}]\\!]\nu\right\rangle_{\mathcal{E}_{h}}=0\qquad\text{for all
$\nu\in\mathbb{W}_{2}^{\text{tr}}$}.$ (18)

This is the crucial step that makes the hybridisable method consistent with
the mixed method that it is based on, since it ensures that the resulting
velocity field is H(div) conforming.

As above, matrix representations $G^{b}$, $K$ and $C$ of the bilinear forms
$g^{b}(\cdot,\cdot)$, $k(\cdot,\cdot)$ and $c(\cdot,\cdot)$ are obtained by
evaluating them for the basis functions. This leads to the following matrix
equation

$\begin{pmatrix}\breve{M}^{b}_{u}-Q_{22}^{b}&G^{b}&K\\\
Q_{32}^{b}+M_{3}^{\rho}M_{3}^{-1}D^{b}&M_{3}^{\Pi}&0\\\
C^{T}&0&0\end{pmatrix}\begin{pmatrix}\widetilde{u}^{\prime b}\\\
\widetilde{\Pi}^{\prime}\\\
\widetilde{\lambda}\end{pmatrix}=\begin{pmatrix}\mathcal{R}^{b}_{u}\\\
\overline{\mathcal{R}}_{\Pi}\\\ 0\end{pmatrix}.$ (19)

Since $\mathbb{W}_{2}\subset\mathbb{W}_{2}^{b}$, the vector
$\mathcal{R}_{u}^{b}$ on the right-hand side is readily obtained from
$\mathcal{R}_{u}$ with a simple basis transformation. The $3\times 3$ system
in Eq. (19) is larger than the original $2\times 2$ system in Eq. (9), both
due to the fact that $\mathbb{W}_{2}^{b}$ has more unknowns than
$\mathbb{W}_{2}$ and since the Lagrange multipliers $\lambda$ need to be
solved for as well. However, the $3\times 3$ matrix in Eq. (19) has the
crucial advantage that its upper left $2\times 2$ submatrix has a block-
diagonal structure, with each block only coupling the unknowns in a single
grid cell. This implies that it can be inverted very efficiently by inverting
small, cell-local matrices. Furthermore, eliminating velocity and pressure
from Eq. (19) results in a (block-)sparse Schur complement for
$\widetilde{\lambda}$, and the velocity $\widetilde{u}^{\prime}$ and Exner
pressure $\widetilde{\Pi}^{\prime}$ can be recovered with purely local
operations if the dof-vector $\widetilde{\lambda}$ is known. To make this more
explicit, we rewrite Eq. (19) as

$\begin{pmatrix}\mathbf{A}&\mathbf{K}\\\
\mathbf{C}^{T}&0\end{pmatrix}\begin{pmatrix}\widetilde{\mathbf{X}}^{\prime}\\\
\widetilde{\lambda}\end{pmatrix}=\begin{pmatrix}\bm{\mathcal{R}}_{X}\\\
0\end{pmatrix},$ (20)

with

$\displaystyle\mathbf{A}$ $\displaystyle\equiv$
$\displaystyle\begin{pmatrix}\breve{M}^{b}_{u}-Q_{22}^{b}&G^{b}\\\
Q_{32}^{b}+M_{3}^{\rho}M_{3}^{-1}D^{b}&M_{3}^{\Pi}\end{pmatrix},$ (21a)
$\displaystyle\mathbf{K}$ $\displaystyle\equiv$
$\displaystyle\left[K,0\right]^{T},$ (21b) $\displaystyle\mathbf{C}$
$\displaystyle\equiv$ $\displaystyle\left[C,0\right]^{T}`,$ (21c)
$\displaystyle\widetilde{\mathbf{X}}^{\prime}$ $\displaystyle\equiv$
$\displaystyle\left[\widetilde{u}^{\prime
b},\widetilde{\Pi^{\prime}}\right]^{T},$ (21d)
$\displaystyle\bm{\mathcal{R}}_{X}$ $\displaystyle\equiv$
$\displaystyle\left[\mathcal{R}_{u}^{b},\overline{\mathcal{R}}_{\Pi}\right]^{T},$
(21e)

where bold operators and vectors indicate objects that act on a combination of
function spaces. Eliminating $\widetilde{\mathbf{X}}^{\prime}$ from (20)
results in an equation for the Lagrange multipliers $\widetilde{\lambda}$

$S\widetilde{\lambda}=\mathcal{B}_{\lambda}$ (22)

with

$\displaystyle S$ $\displaystyle\equiv$
$\displaystyle\mathbf{C}^{T}\mathbf{A}^{-1}\mathbf{K}\qquad\text{and}$ (23)
$\displaystyle\mathcal{B}_{\lambda}$ $\displaystyle\equiv$
$\displaystyle\mathbf{C}^{T}\mathbf{A}^{-1}\bm{\mathcal{R}}_{X}.$ (24)

Since $\mathbf{A}^{-1}$ is block-diagonal and only couples unknown within a
single grid cell, the Schur complement matrix $S$ in Eq. (23) is (block-)
sparse and can be applied cheaply while iteratively solving Eq. (22) with a
Krylov subspace method. Observe further that instead of iterating in the
pressure-velocity space as in Eq. (9), we now only need to iterate in the
smaller trace space: at lowest order, we need to store on average three trace
unknowns per cell (one shared trace unknown for each of the six facets). The
right-hand side $\mathcal{B}_{\lambda}$ in Eq. (24), which can be interpreted
as the matrix-representation of a one-form in the dual space
$\mathbb{W}_{2}^{\text{tr}*}$, can be constructed using purely local
operations. Once the solution $\widetilde{\lambda}$ of Eq. (22) is known, the
Exner pressure $\widetilde{\Pi}^{\prime}$ and (broken) velocity increment
$\widetilde{u}^{\prime b}$ can be obtained via the local operation

$\widetilde{\mathbf{X}}^{\prime}=\mathbf{A}^{-1}\left(\bm{\mathcal{R}}_{X}-\mathbf{K}\widetilde{\lambda}\right).$
(25)

###  2.5 Non-nested multigrid

To solve Eq. (22) we use a non-nested multigrid algorithm similar to that
described in [31]. The central idea is the following: as for any multigrid
algorithm, first a small number of smoother iterations are applied to reduce
high-frequency error components. Since the Lagrange multiplier field
$\lambda\in\mathbb{W}_{2}^{\text{tr}}$ approximates the Exner pressure
increment, we then restrict the residual to the pressure space
$\mathbb{W}_{3}$ and solve the resulting coarse level equation with the
pressure multigrid algorithm described in Section 2.3, before prolongating the
coarse level correction back to $\mathbb{W}_{2}^{\text{tr}}$. The same
algorithm has been proposed for the mixed Poisson problem in [31], with the
crucial difference that the authors of [31] employ a conforming coarse level
space instead of the discontinuous space used in this work. We begin by
discussing an efficient smoother algorithm for Eq. (22).

#### Trace space smoother

As for the pressure multigrid algorithm in Section 2.3, the smoother in the
trace space $\mathbb{W}_{2}^{\text{tr}}$ needs to take account of the strong
anisotropy caused by the high aspect ratio of the grid ($\Delta z\ll\Delta
x$). First note that the trace space $\mathbb{W}_{2}^{\text{tr}}$ can be
written as the direct sum
$\mathbb{W}_{2}^{\text{tr}}=\mathbb{W}_{2}^{\text{tr,h}}\oplus\mathbb{W}_{2}^{\text{tr,z}}$.
As shown in Figure 1, the degrees of freedom of $\mathbb{W}_{2}^{\text{tr,h}}$
are associated with vertical facets (aligned parallel to the vertical
direction $\widehat{\bm{z}}$) while the degrees of freedom of
$\mathbb{W}_{2}^{\text{tr,z}}$ are associated with horizontal facets (aligned
perpendicular to the vertical direction $\widehat{\bm{z}}$). Similar to
$\mathbb{W}_{2}^{\text{z}}$, the function space $\mathbb{W}_{2}^{\text{tr,z}}$
is continuous in the vertical but discontinuous in the horizontal direction,
while $\mathbb{W}_{2}^{\text{tr,h}}$ is horizontally continuous and vertically
discontinuous.

The decomposition of $\mathbb{W}_{2}^{\text{tr}}$ implies that the dof-vector
$\widetilde{\lambda}=(\widetilde{\lambda}^{\text{h}},\widetilde{\lambda}^{\text{z}})$
can be partitioned such that the unknowns associated with
$\mathbb{W}_{2}^{\text{tr,h}}$ are collected in the vector
$\widetilde{\lambda}^{\text{h}}$, while the unknowns associated with
$\mathbb{W}_{2}^{\text{tr,z}}$ are stored in $\widetilde{\lambda}^{\text{z}}$.
Further, since $\mathbb{W}_{2}^{\text{tr,z}}$ is discontinuous in the
horizontal direction, we can write
$\widetilde{\lambda}^{\text{z}}=(\widetilde{\lambda}^{\text{z}}_{1},\widetilde{\lambda}^{\text{z}}_{2},\dots,\widetilde{\lambda}^{\text{z}}_{n_{\text{col}}})$
where $\widetilde{\lambda}^{\text{z}}_{j}$ is the dof-vector of all unknowns
associated with a particular column $j$. With this ordering, the matrix $S$ in
Eq. (22) has the following block-structure:

$S=\begin{pmatrix}S^{\text{hh}}&S^{\text{hz}}\\\
S^{\text{zh}}&S^{\text{zz}}\end{pmatrix}\qquad\text{with}\qquad
S^{\text{zz}}=\text{blockdiag}\left(S^{\text{zz}}_{1},S^{\text{zz}}_{2},\dots,S^{\text{zz}}_{n_{\text{col}}}\right)$
(26)

Due to the strong vertical anisotropy, the entries of the matrix
$S^{\text{zz}}$ are much larger than all other entries in $S$. Furthermore,
the block-diagonal matrix $S^{\text{zz}}$ can be inverted by independently
solving a tridiagonal system in each vertical column of the grid.

The smoother is written down explicitly in Algorithm 1, and has a very similar
structure to the block-Jacobi method in Eq. (14) that is used as a smoother
for the pressure multigrid algorithm in Section 2.3. While we use a block-
Jacobi iteration for the degrees of freedom associated with
$\mathbb{W}_{2}^{\text{tr,z}}$, a simple point Jacobi procedure is used for
the unknowns associated with $\mathbb{W}_{2}^{\text{tr,h}}$. Recall that
solving a system of the form

$S_{j}^{\text{zz}}\widetilde{\delta\lambda}_{j}=(\mathcal{R}_{\lambda}^{\text{z}})_{j}$
(27)

in line 5 of Algorithm 1 with the Thomas algorithm [38] will incur a cost of
$\mathcal{O}(L)$ where $L$ is the number of vertical layers of the grid (at
higher order $S^{\text{zz}}_{j}$ is a banded matrix with bandwidth $b$, which
can be solved in $\mathcal{O}(bL)$ time if we ignore the initial overhead of
the banded factorisation).

Algorithm 1
$\textsf{LineSmooth}(\widetilde{\lambda};S,\mathcal{B}_{\lambda};\omega,n_{\text{smooth}})$:
update $\widetilde{\lambda}$ with $n_{\text{smooth}}$ (block-) Jacobi
iterations with relaxation factor $\omega$

1:for $k=1,2,\dots,n_{\text{smooth}}$ do

2: Compute
$\mathcal{R}_{\lambda}\mathrel{\reflectbox{$\mapsto$}}\mathcal{B}_{\lambda}-S\widetilde{\lambda}$

3: Partition
$\mathcal{R}_{\lambda}=(\mathcal{R}_{\lambda}^{\text{h}},\mathcal{R}_{\lambda}^{\text{z}})$
and
$\mathcal{R}_{\lambda}^{\text{z}}=((\mathcal{R}_{\lambda}^{\text{z}})_{1},(\mathcal{R}_{\lambda}^{\text{z}})_{2},\dots,(\mathcal{R}_{\lambda}^{\text{z}})_{n_{\text{col}}})$
in the same way as $\widetilde{\lambda}$

4: for $j=1,2,\dots,n_{\text{col}}$ do

5: Solve the tridiagonal system
$S_{j}^{\text{zz}}\widetilde{\delta\lambda}_{j}=(\mathcal{R}_{\lambda}^{\text{z}})_{j}$
for $\widetilde{\delta\lambda}_{j}$

6: Update
$\widetilde{\lambda}^{\text{z}}_{j}\mathrel{\reflectbox{$\mapsto$}}\widetilde{\lambda}^{\text{z}}_{j}-\omega\widetilde{\delta\lambda}_{j}$

7: end for

8: Update
$\widetilde{\lambda}^{\text{h}}\mathrel{\reflectbox{$\mapsto$}}\widetilde{\lambda}^{\text{h}}-\omega\left(\text{diag}(S^{\text{hh}})\right)^{-1}\mathcal{R}_{\lambda}^{\text{h}}$

9:end for

Of course, the smoother in Algorithm 1 could be used directly to precondition
a Krylov-subspace method for the solution of
$S\widetilde{\lambda}=\mathcal{B}_{\lambda}$, and in our numerical results in
Section 3 we also investigate this setup as a simpler reference method. As
expected for a single-level preconditioner, the number of Krylov iterations
turns out to be very large, which leads to a larger overall solution time. The
hierarchical algorithm that we describe next overcomes this issue.

#### Non-nested two level algorithm

A single iteration of the non-nested two-level method is described in
Algorithm 2, where
$\textsc{LineSmooth}(\widetilde{\lambda};S,\mathcal{B}_{\lambda},\omega,n_{\text{smooth}})$
corresponds to $n_{\text{smooth}}$ applications of the linesmoother in
Algorithm 1. Crucially, the operator $H$ in the coarse level system

$H\widetilde{\delta\Pi}=\overline{\mathcal{B}}$ (28)

in line 6 of Algorithm 2 is the operator defined in Eq. (12), and we can
therefore solve Eq. (28) with the pressure multigrid algorithm from Section
2.3.

Algorithm 2 Non-nested two-level cycle for approximately solving
$S\widetilde{\lambda}=\mathcal{B}_{\lambda}$

1:Set $\widetilde{\lambda}\mapsto 0$

2:Presmooth:
$\textsc{LineSmooth}(\widetilde{\lambda};S,\mathcal{B}_{\lambda};\omega,n_{\text{presmooth}})$
(see Algorithm 1)

3:Compute residual:
$\mathcal{R}_{\lambda}\mathrel{\reflectbox{$\mapsto$}}\mathcal{B}_{\lambda}-S\widetilde{\lambda}$

4:Restrict residual:
$\mathcal{B}\mathrel{\reflectbox{$\mapsto$}}\textsc{Restrict}\left(\mathcal{R}_{\lambda}\right)$

5:Scale right-hand side:
$\overline{\mathcal{B}}\mathrel{\reflectbox{$\mapsto$}}\Gamma\mathcal{B}$ (see
Eq. (31))

6:Solve pressure system $H\widetilde{\delta\Pi}=\overline{\mathcal{B}}$ for
$\widetilde{\delta\Pi}$.

7:Prolongate correction and update solution:
$\widetilde{\lambda}\mathrel{\reflectbox{$\mapsto$}}\widetilde{\lambda}+\textsc{Prolong}\left(\widetilde{\delta\Pi}\right)$

8:Postsmooth:
$\textsc{LineSmooth}(\widetilde{\lambda};S,\mathcal{B}_{\lambda};\omega,n_{\text{postsmooth}})$
(see Algorithm 1)

#### Intergrid operators

To prolongate the solution from the pressure space $\mathbb{W}_{3}$ to the
trace space $\mathbb{W}_{2}^{\text{tr}}$ we use the fact that
$\lambda\in\mathbb{W}_{2}^{\text{tr}}$ is an approximation of the Exner
pressure increment $\Pi^{\prime}\in\mathbb{W}_{3}$. Hence, the natural
prolongation
$\mathfrak{P}:\mathbb{W}_{3}\rightarrow\mathbb{W}_{2}^{\text{tr}}$ is the
averaging procedure

$\mathfrak{P}:\Pi^{\prime}\mapsto\lambda=\\{\\!\\!\\{\Pi^{\prime}\\}\\!\\!\\}.$
(29)

The corresponding linear operator acting on the dof-vector
$\widetilde{\Pi}^{\prime}$ can be written as
$\textsf{Prolong}\left(\widetilde{\Pi}^{\prime}\right)\equiv
P\widetilde{\Pi}^{\prime}$. At lowest order, the application of the matrix $P$
simply corresponds to averaging the values of the unknowns on neighbouring
cells to obtain the trace-unknown of the facet between those two cells.
Restricting the dof-vector $\widetilde{\lambda}$ of a field in the trace-space
$\mathbb{W}_{2}^{\text{tr}}$ to the dof-vector of a field in the pressure
space $\mathbb{W}_{3}$ then corresponds to the multiplication with the
transpose matrix, i.e.
$\textsf{Restrict}\left(\widetilde{\lambda}\right)\equiv
P^{T}\widetilde{\lambda}$.

#### Coarse level operator

Given a fine level matrix $S$, a prolongation matrix $P$ and a restriction
matrix $R=P^{T}$, the canonical way of constructing a suitable coarse level
operator is to form the Galerkin product $H_{\text{G}}=P^{T}SP$. This leads to
the coarse level equation
$H_{G}\widetilde{\delta\Pi}=\mathcal{B}=R^{T}\mathcal{R}_{\lambda}$.
Unfortunately, the triple matrix-product required in the construction of
$H_{G}$ is expensive and $H_{G}$ has a large stencil. To address these issues
we use a different approach here. We know that $H$ as defined in Eq. (12)
(which only contains couplings between adjacent cells) is a good approximation
to $H_{G}$, in the sense that the action of $H_{G}$ on smooth fields is
similar to the action of $H$, up to some scaling factor $\Gamma$, which
appears in line 5 of Algorithm 2. More precisely,

$\Gamma H_{G}\widetilde{\Pi}_{\text{smooth}}\approx
H\widetilde{\Pi}_{\text{smooth}}$ (30)

where $\Gamma$ is a diagonal matrix and $\widetilde{\Pi}_{\text{smooth}}$ is
the dof-vector of a $\mathbb{W}_{3}$ field that varies slowly in space. To
find the entries of $\Gamma$, we consider the constant field
$\widetilde{\Pi}_{\text{const}}$ with $(\widetilde{\Pi}_{\text{const}})_{j}=1$
for all $j$. Then $\Gamma_{jj}$ is given by the ratio of the $j$-th row-sums
of $H$ and $H_{G}$:

$\Gamma_{jj}=\frac{\sum_{k}\left(H\right)_{jk}}{\sum_{k}\left(H_{G}\right)_{jk}}=\frac{(H\widetilde{\Pi}_{\text{const}})_{j}}{(H_{G}\widetilde{\Pi}_{\text{const}})_{j}}=\frac{(H\widetilde{\Pi}_{\text{const}})_{j}}{(P^{T}SP\widetilde{\Pi}_{\text{const}})_{j}}.$
(31)

Evaluating the expressions $H\widetilde{\Pi}_{\text{const}}$ and
$P^{T}SP\widetilde{\Pi}_{\text{const}}$ that are needed in the numerator and
denominator of the rightmost fraction in Eq. (31) requires only
multiplications by the sparse matrices $S$, $H$, $P$ and $P^{T}$.

The coarse level problem $H\widetilde{\delta H}=\overline{\mathcal{B}}$ in
line 6 of Algorithm 2 is solved with the pressure multigrid algorithm
described in Section 2.3.

##  3 Results

###  3.1 Implementation

The non-nested multigrid solver described in Section 2.5 has been implemented
in the full LFRic model, which uses a separation of concerns philosophy to
orchestrate the efficient parallel execution of computational kernels [21]. In
our case, these kernels correspond, for example, to the application of the
linear operator $S$ defined in Eq. (23) in a vertical column, or a single
trididiagonal solve in Algorithm 1. From the discussion in Section 2 it should
be evident that the solver architecture is highly non-trivial, and it is
crucial to have a simple yet robust mechanism for independently testing and
swapping different solver components. This allows a comparison to be made to
alternative approaches, such as the approximate Schur complement method in
[25], to ultimately identify the most promising solver configuration. To
realise this, a hierarchy of Fortran 2003 derived types has been developed
that conform to common interfaces: following the same modular design as other
libraries [39, 40, 41, 42], the “Solver API” in LFRic [25, 21] allows the
implementation of iterative solvers, linear operators and preconditioners
which can be plugged together by the algorithm developer. Further details on
the modular implementation of the solver algorithms discussed in this paper
can be found in Appendix A.

###  3.2 Solver configurations

This paper introduces two innovations: the solution of the linear system in
Eqs. (1a) - (1d) with the help of a hybridised discretisation and efficient
preconditioning of the hybridised system in Eq. (22) with a non-nested
multigrid algorithm. To assess the impact of each of these changes on model
performance, we compare three different solver setups:

Configuration 1 (Pressure multigrid):

The current non-hybridised solver introduced in [25], preconditioning the
mixed system in Eq. (9) with the approximate Schur complement algorithm that
solves the pressure system with multigrid as described in Section 2.3. This
configuration can be considered as a baseline for our comparisons.

Configuration 2 (Hybridised single-level):

The hybridised solver as described in Section 2.4, preconditioning the trace
system in Eq. (22) with a few iterations of the block-Jacobi smoother in
Algorithm 1.

Configuration 3 (Hybridised multigrid):

The hybridised solver as described in Section 2.4, preconditioning the trace
system in Eq. (22) with the non-nested multigrid algorithm in Algorithm 2.

The gains in solver performance achieved by hybridising the equations can be
assesed by comparing configurations 1 and 3, using robust multigrid solvers
for both discretisations. However, given the complexity of the non-nested
multigrid algorithm in Algorithm 2, it is worth assessing whether the
hybridised solver, using a much simpler single-level preconditioner, is
competitive compared to the non-hybridised solver. As in [25], we therefore
compare our hybridised multigrid preconditioner (Configuration 3) to the
block-Jacobi preconditioner (Configuration 2).

For the pressure multigrid preconditioner (Configuration 1), the mixed system
in Eq. (9) is solved with BiCGStab [43, 44], reducing the two-norm of the
residual by six orders of magnitude. The same BiCGStab iteration and relative
residual reduction is used to solve the system in Eq. (22) for the two
hybridised solvers (Configurations 2 and 3).

For the hybridised single-level preconditioner (Configuration 2), we can apply
different numbers of smoothing steps, with a trade-off between the total
number of BiCGStab iterations and the time per iteration. Thus, we consider
three different setups, with $n_{\text{smooth}}=$ 1, 2 and 3 in Algorithm 1.
Note that using $n_{\text{smooth}}=3$ then uses the equivalent number of
smoothing steps as the hybridised multigrid setup (Configuration 3), which
uses $n_{\text{presmooth}}=1$ and $n_{\text{postsmooth}}=2$.

In the non-nested two-level method in Algorithm 2, the coarse level pressure
system in Eq. (28) can be solved in different ways. Since Algorithm 2 is used
as a preconditioner for a BiCGStab iteration and the coarse level solution
$\widetilde{\delta H}$ is an approximation to the solution of the residual
equation, it is usually sufficient to only solve it to some loose tolerance,
which will lead to better overall performance. To explore this, we consider
three different setups for the coarse level solve:

Configuration 3a (exact coarse solve):

Solve Eq. (28) to a very tight tolerance, reducing the relative residual norm
below $10^{-14}$.

Configuration 3b (approximate coarse solve):

Solve Eq. (28) to a loose tolerance, reducing the residual norm by two orders
of magnitude only.

Configuration 3c (pressure multigrid V-cycle coarse solve):

Apply one V-cycle of the pressure multigrid algorithm from Section 2.3 to
approximately solve Eq. (28) as cheaply as possible.

In Configurations 3a and 3b, the coarse level solver is BiCGStab,
preconditioned with the existing pressure multigrid algorithm from Section
2.3. Configuration 3a, which essentially solves the coarse level problem in
Eq. (28) to machine precision, is useful to debug any issues with Algorithm 2.
It is worth recalling that we used similar setups in [25] to verify that
solving the pressure correction there with a single multigrid V-cyle gives the
best overall performance.

For reference, Table 1 summarises the different solver configurations that
were used for the numerical experiments below.

Config | mixed solve+precon | inner solve+precon | coarse solve+precon
---|---|---|---
1 | $\textsc{BiCGStab}(10^{-6})+\text{approx Schur}$ | $\textsc{pre-only}+\textsc{MG}$ | —
2 | pre-only+hybridisation |  $\textsc{BiCGStab}(10^{-6})$ \+ LineSmooth | —
3a | pre-only+hybridisation |  $\textsc{BiCGStab}(10^{-6})$+non-nested |  $\textsc{BiCGStab}(10^{-14})$+MG
3b | pre-only+hybridisation |  $\textsc{BiCGStab}(10^{-6})$+non-nested |  $\textsc{BiCGStab}(10^{-2})$+MG
3c | pre-only+hybridisation |  $\textsc{BiCGStab}(10^{-6})$+non-nested |  pre-only+MG
Table 1: Summary of solver configurations. Using the “pre-only” solve
corresponds to a single application of the preconditioner. “LineSmooth” is the
block-Jacobi smoother in Algorithm 1, “non-nested” refers to Algorithm 2
whereas “MG” stands for the pressure multigrid algorithm described in Section
2.3.

###  3.3 Model setup

To obtain the results which follow, each configuration in Table 1 is run for
12 timesteps with a horizontal wave Courant number of approximately 12; this
mimics what would be used in operational settings. The pressure multigrid
preconditioner uses 4 multigrid levels and 1 pre-smoothing and 2 post-
smoothing steps with a relaxation factor of $\omega_{\text{pressure}}=0.9$;
this setup is used both in Configuration 1 and when pressure multigrid is used
as a coarse level solver in Configuration 3. The same number of pre- and post-
smoothing steps are also used for the non-nested two-level algorithm in
Algorithm 2. However, in this algorithm and for the smoother that
preconditions the single-level method in Configuration 2, we set $\omega=0.6$.
The dependence of solver performance on the different parameters is explored
in Appendix B, where it is demonstrated that the chosen values give good
performance in general. As discussed in [25] and initially observed in [45],
it is sufficient to use a small number of multigrid levels. This is because
for the chosen Courant numbers the condition number of Eq. (28) is reduced to
$\approx 1$ on the coarsest multigrid level, independent of the horizontal
resolution. We perform all experiments in parallel, using the Met Office Cray
XC40 supercomputer which uses the Aries interconnect. Each node comprises
dual-socket, 18-core Broadwell Intel Xeon processors, that is, 36 CPU cores
per node. The model was compiled with the Intel 17 Fortran compiler.

We simulate the full compressible Euler equations for a perfect gas, as set
out in Section 2.1. The equations are solved in a spherical domain without
orography, which is covered with a cubed sphere mesh that contains $6\times
96\times 96$ vertical columns (this mesh is referred to as “C96” below). Each
column is subdivided into 30 vertical layers with a grid spacing which
increases quadratically with height above ground. Overall, the mesh contains
$1,658,880$ grid cells and the average horizontal grid spacing is $\approx
96~{}\text{km}$. The model timestep size is set to $1~{}\text{hour}$,
resulting in a horizontal acoustic Courant number of $\approx 12$. The radius
of the earth is assumed to be $6371.229~{}\text{km}$; the vertical lid is
fixed at $30~{}\text{km}$ above ground. All calculations were carried out in
double precision (64bit floating point) arithmetic.

###  3.4 Single-node Performance

Two key metrics are considered to assess the performance of each solver
configuration:

  1. 1.

The average time per model timestep.

  2. 2.

The convergence history of the preconditioned iterative solvers (for Eq. (9)
in Configuration 1 and Eq. (22) in Configurations 2 and 3).

While the ultimate goal is to reduce the total runtime, the convergence
history provides valuable insights into the algorithmic performance of the
solvers and can highlight the potential for further optimisation. For
multigrid we typically expect the residual to be reduced by one order of
magnitude per iteration, while (as already demonstrated in [25]) convergence
is much slower for single-level preconditioners; this is confirmed by
comparing Configurations 2 and 3 in the following results.

Figure 2 shows the average time per timestep for all the considered solver
configurations, running the model on a single node (36 cores) of the Met
Office Cray XC40 supercomputer and averaging over all 12 model timesteps. The
pressure multigrid setup (Configuration 1) has been further separated into a
method which iterates over the $4\times 4$ system as in [25] and the new
method described here, which eliminates potential temperature and density as
described in Sections 2.1 and 2.2 to obtain a $2\times 2$ system in Eq. (9).
While the old method takes $7.60\text{s}$ per timestep, iterating over the
$2\times 2$ system reduces this to $6.83\text{s}$, which demonstrates that
eliminating potential density and temperature before the iteration improves
performance by around $10\%$. For the rest of this paper only the $2\times 2$
system will be considered for Configuration 1.

Figure 2: Average time per timestep for the full compressible Euler equations
in LFRic.

We observe that the non-hybridised pressure multigrid solver (Configuration 1)
gives the best performance (6.82 s), but is closely followed by our hybridised
multigrid solver (Configuration 3c), which is around $4\%$ slower (7.08 s).
Solving the coarse level system in the non-nested two-level solver to a
tighter tolerance increases the runtime to 22.63 s for Configuration 3a and
8.92 s for Configuration 3b. One striking feature of Figure 2 is that the
hybridised single-level solver (Configuration 2) gives the worst runtime
overall – 49.29 s, 44.64 s, and 56.06 s for $n_{\text{smooth}}=1$, 2, and 3
respectively. In other words, replacing the linesmoother preconditioner in the
hybridised solver with the non-nested multigrid algorithm reduces the runtime
by a factor of more than $6\times$ and only with this change does our
hybridised solver become competitive with the pressure multigrid setup in
[25].

Figure 3 shows the convergence histories for the relative residual
$||r||/||r_{0}||$.

Figure 3: Convergence history of the iterative solvers (left), with a zoom to
the first 20 iterations (right).

We observe that the convergence rates for the three variations of the non-
nested two-level preconditioner (Configuration 3) are comparable to that
observed for the pressure multigrid solver (Configuration 1). In contrast, the
convergence rate of the hybridised single-level solver (Configuration 2) is
prohibitively slow and this is reflected in the time to solution in Figure 2:
more than 400 iterations are required to reduce the residual by eight orders
of magnitude for $n_{\text{smooth}}=1$, and in fact Configuration 2 with
$n_{\text{smooth}}=1$ fails to reduce the residual by a relative factor
$~{}10^{-9}$. In contrast, to reduce the residual by twelve orders of
magnitude our hybridised multigrid solver (Configuration 3c) requires 16
iterations – which is textbook multigrid performance – and the pressure
multigrid solver (Configuration 1) takes slightly longer to converge with 21
iterations.

The convergence history also explains the poor performance of Configurations
3a and 3b in Figure 2: Figure 3 demonstrates that the multigrid convergence
rate is virtually independent of the coarse level solution method. This is
readily explained by the fact that the coarse level equation in Eq. (28) is an
approximation of Eq. (22), so there is not much point in solving it to high
accuracy. This mirrors the findings in [25], where we showed that the best
performance is achieved by approximately solving the pressure correction
equation with a single multigrid V-cycle.

The convergence history for the hybridised single-level solver (Configuration
2) shows large spikes, with the residual increasing by more than one order of
magnitude from one iteration to the next in some cases. In practice, this can
cause serious issues when determining whether the solver has converged. It
also potentially makes the runtime unpredictable, which is undesirable for
operational forecast models that must run in a fixed time window. More
importantly, the slow convergence of the single-level solver will likely cause
serious issues when running the model in single precision floating point
arithmetic – the Krylov subspace vectors lose orthogonality in reduced
precision calculations and the convergence will stall.

Comparing the different number of smoothing applications $n_{\text{smooth}}$
for the hybridised single-level solver (Configuration 2) illustrates that the
linesmoother is effective in reducing the number of iterations to convergence.
However, this of course results in a greater cost per BiCGStab iteration and
thus the hybridised single-level solver is (as expected) not competitive with
the pressure multigrid (Configuration 1) and hybridised multigrid
(Configuration 3) solvers. We also observe that there is a limit to what the
smoother can do: the convergence improves to a much lesser extent when
increasing from $n_{\text{smooth}}=2$ to $n_{\text{smooth}}=3$, which means
that the runtime using $n_{\text{smooth}}=3$ is greatly increased, and is in
fact very similar to the runtime when using $n_{\text{smooth}}=1$.

Interestingly, in contrast to what was observed for a simplified gravity wave
model (not shown here) and the Krylov subspace solvers with single-level
preconditioners in the ENDGame model (also not shown here), in Figure 3 the
residual for the single-level method (Configuration 2) does not drop rapidly
in the first few iterations. Hence, even if the residual is reduced by less
than six orders of magnitude, the hybridised multigrid preconditioner
(Configuration 3c) will likely retain its advantage over the single-level
method (Configuration 2).

###  3.5 Breakdown of runtimes

To gain further insight into the different configurations, we also perform
timings on key components of the solver. Only the timings from the pressure
multigrid setup (Configuration 1) and our hybridised multigrid solver with the
“precondition only” setup (Configuration 3c) are shown here due to the much
greater cost of the other solver configurations. The results for these two
configurations are presented in the bar chart in Figure 4, where the total
height of each bar represents the total time spent in the mixed solver; this
is the average time it takes to solve Eq. (9) with a suitably preconditioned
Krylov subspace method. The figure shows the total time spent in all four
linear solves that are required per timestep. Note that this is different to
the times presented in Figure 2 which are the times for executing an entire
timestep (including, for example, the construction of the operators and non-
linear residual calculation). Each bar is then broken down into key
components, for example the time spent in solving the pressure correction
equation. For Configuration 3c, “Pressure solve” refers to the time spent in
the solution of the coarse level pressure system, whereas “Hybrid solve”
measures the combined time spent in the application of the linesmoother
(Algorithm 1) and trace operator $S$ in Eq. (23). Together, these times give a
good estimate for the total time it takes to solve Eq. (22). For Configuration
1 “Mixed operator” is the time spent in the application of the $2\times 2$
operator in Eq. (9), whereas “Forward/backward substitution” measures the
construction of the right hand side for the pressure equation in Eq. (11) and
the recovery of the velocity field from the pressure solve in Eq. (13).

Figure 4: Bar chart illustrating the breakdown of time spent in the mixed
solve.

For the pressure multigrid setup (Configuration 1), the solution of the
pressure correction system accounts for about 16 % of the solve time, which is
dominated by other parts of the mixed solver, such as in the application of
the mixed operator and forward/backward substitution. The hybridised multigrid
solver (Configuration 3c) spends a significant fraction of the total time in
application of the hybridised operator and the linesmoother (namely in
performing the hybrid solve), while the pressure multigrid is also important
here.

As a comparison, it is worth noting that the hybridised single-level solver
with $n_{\text{smooth}}=1$ (Configuration 2a) spends around 46.8 s per
timestep in the mixed solver, 72 % of which is in the application of the
hybridised operator and the linesmoother, which are called frequently in the
preconditioned Krylov solve. As above, it is clear that the slow convergence
rate greatly increases the time spent in the mixed solver – the pressure
multigrid setup (Configuration 1) and the hybridised multigrid solver with
precondition only (Configuration 3c) spend 10$\times$ less time in solving the
linear system in Eq. (9) (around 4.7 s and 4.6 s respectively).

For reference, we also present a breakdown of the total time (integrated over
the entire run) and the time per call for key solver components in Table 2.

| Time [s]
---|---
| Config 1 | Config 2a | Config 3c
Component | Total | per call | Total | per call | Total | per call
pressure operator | $5.13$ | $0.00048$ | — | — | $3.97$ | $0.00041$
approximate Schur complement preconditioner | $8.90$ | $0.029$ | — | — | — | —
pressure multigrid V-cycle | $9.33$ | $0.010$ | — | — | $7.10$ | $0.087$
hybridised operator | — | — | $284.40$ | $0.0085$ | $13.83$ | $0.0084$
hybridised forward/backward substitution | — | — | $2.10$ | $0.044$ | $2.13$ | $0.044$
non-nested multigrid algorithm (inc. pressure solve) | — | — | — | — | $42.29$ | $0.052$
single-level preconditioner for hybridised system | — | — | $131.33$ | $0.0039$ | $25.61$ | $0.016$
Table 2: Summary of the total time and time spent per call in several key
components of the three solver configurations.

###  3.6 Parallel scaling

The metrics considered above give a useful insight into how each solver
performs per timestep, thus highlighting the algorithmic advantages of
different methods. However, a key factor in choosing solvers is how they scale
within massively parallel architectures.

To this end, we repeat the numerical experiments in Section 3.4 at scale to
investigate the parallel performance of the two most promising solver
configurations, namely the pressure multigrid setup (Configuration 1) and the
hybridised multigrid solver (Configuration 3c). In all cases the timestep size
is adjusted such that the horizontal acoustic Courant number is $\approx 12$
as for the experiments in Section 3.4. For the weak scaling analysis the model
is run on 6, 24, 96 and 384 nodes (recall that each node consists of 36
cores). The highest horizontal resolution that was considered is 12km on
13,824 cores. Details on the setup for these runs are collected in Table 3.

# nodes | # cores | mesh | total # grid cells | resolution $\Delta x$ | timestep $\Delta t$
---|---|---|---|---|---
6 | 216 | C96 | $6\times 96\times 96\times 30$ |  = | 1,658,880 | $96\text{km}$ | $3600\text{s}$
24 | 864 | C192 | $6\times 192\times 192\times 30$ |  = | 6,635,520 | $48\text{km}$ | $1800\text{s}$
96 | 3,456 | C384 | $6\times 384\times 384\times 30$ |  = | 26,542,080 | $24\text{km}$ | $900\text{s}$
384 | 13,824 | C768 | $6\times 768\times 768\times 30$ |  = | 106,168,320 | $12\text{km}$ | $450\text{s}$
Table 3: Setup of weak scaling runs.

All runs use pure MPI parallelisation with 36 MPI regions per node, giving a
local volume of $16^{2}$ columns (all with 30 vertical layers) per region for
the weak scaling runs.

The strong scaling analysis is carried out for 6, 24, 96 and 384 nodes on a
C384 mesh with $6\times 384\times 384\times 30=26,542,080$ grid cells and a
resolution of 24km, fixing the timestep size to $\Delta t=900\text{s}$. As the
number of processors increases, the local volume decreases from $64^{2}$
columns on 6 nodes to $8^{2}$ columns on 384 nodes; note that since four
multigrid levels are used, this corresponds to a local volume of just 1 column
(30 cells) per MPI region on the coarsest level of the multigrid hierarchy.
More details on the setup for the strong scaling runs can be found in Table 4.

| # columns per core | # grid cells per core
---|---|---
# nodes | # cores | finest | coarsest | finest | coarsest
6 | 216 | $64\times 64$ | $8\times 8$ | 122,880 | 1,920
24 | 864 | $32\times 32$ | $4\times 4$ | 30,720 | 480
96 | 3,456 | $16\times 16$ | $2\times 2$ | 7,680 | 120
384 | 13,824 | $8\times 8$ | $1\times 1$ | 1,920 | 30
Table 4: Setup of strong scaling runs on the C384 mesh. The number of columns
and grid cells per core are given both for the finest and for the coarsest
mesh in the multigrid hierarchy.

The performance of the solver under both weak- and strong scaling is shown in
Figure 5, in both cases we report the average time per timestep. Denoting this
time by $t_{\text{step}}(p)$ when running on $p$ nodes, the corresponding
parallel efficiencies relative to a 6-node run can be defined as

weak parallel efficiency
$\displaystyle=\frac{t_{\text{step}}(6)}{t_{\text{step}}(p)},$ strong parallel
efficiency $\displaystyle=\frac{6\cdot t_{\text{step}}(6)}{p\cdot
t_{\text{step}}(p)}.$ (32)

These efficiencies are shown in Figure 6. For the strong scaling experiment we
also plot the speedup $t_{\text{step}}(p)/t_{\text{step}}(6)$ relative to a
6-node run in Figure 7.

(a)  (b)

Figure 5: Time per timestep for the pressure multigrid setup (Configuration 1)
and the hybridised multigrid solver (Configuration 3c) under weak scaling
(left) and strong scaling (right).

(a)  (b)

Figure 6: Weak scaling (left) and strong scaling (right) efficiency as defined
in Eq. (32) for the pressure multigrid setup (Configuration 1) and the
hybridised multigrid solver (Configuration 3c). (a)  Figure 7: Parallel
speedup of the time per timestep for the pressure multigrid setup
(Configuration 1) and the hybridised multigrid solver (Configuration 3c).

At this point it should be stressed that the parallel infrastructure for the
LFRic model is under continuous development and optimising parallel
performance for the solver is beyond the scope of this paper. At all
resolutions the timings of the two considered setups are similar, but the
hybrid solver generally exhibits slightly better scaling and results in a
shorter runtime at the highest node counts. Interestingly, the strong scaling
efficiency for both solver configurations exceeds $100\%$ for 24 nodes but
drops for larger node counts, with the pressure multigrid solver exhibiting a
more significant drop in efficiency. This effect is more pronounced for the
strong scaling experiment, where it can possibly be explained by the fact that
improved cache efficiency due to the smaller local volume counteracts a drop
in performance due to communication overheads. The weak scaling efficiency of
the hybrid solver is almost $100\%$ on the largest core count, whereas it
drops to around $90\%$ for the pressure multigrid solver (Configuration 1).
The strong scaling efficiency is significantly reduced on 384 nodes. Recall
that at this point each core will only own the unknowns in a single vertical
column with 30 grid cells. On the finest grid, which accounts for a larger
fraction of the runtime, the number of unknowns is $8\times 8\times 30=1920$.
Common experience shows that codes usually stop to scale when the number of
local unknowns is at the order of several thousands, so our observations are
in line with this.

##  4 Conclusion

In this work we introduced a new hybridisable finite element discretisation
for the solution of the linear equation that arises in semi-implicit
timestepping in the LFRic finite element core. A non-nested multigrid method
is used to precondition the resulting linear system in trace space. Our
numerical experiments show that the new method is competitive with the
pressure multigrid method in [25]. Both methods show good weak- and strong
scaling on parallel architectures.

###  4.1 Discussion

One reason why the approximate Schur complement solver with multigrid
preconditioner (Configuration 1) is still competitive with our hybridised non-
nested multigrid solver (Configuration 3c) is that, at lowest order, the
number of unknowns and the sizes of the local stiffness matrices are
comparable. For example, bearing in mind that the velocity unknowns are shared
between neighbouring cells, there are four unknowns per cell for the $2\times
2$ system in Eq. (9) and the local stiffness matrices are of size $7\times 7$.
For the trace system in Eq. (22), there are three unknowns per cell (the trace
unknowns on the six facets are shared) and the local stiffness matrices are of
size $6\times 6$.

For general polynomial degree $p$, the number of pressure- and velocity
unknowns can be worked out using, for example [46]. For the mixed system in
Eq. (9), the average number of unknowns per cell grows as $\mathcal{O}(p^{3})$
and hence the local matrices have $\mathcal{O}(p^{6})$ entries. Since the
trace unknowns for the hybridised system are associased with the lower-
dimensional facets, the number of degrees of freedom per cell grows only as
$\mathcal{O}(p^{2})$ and the cell-local matrices have $\mathcal{O}(p^{4})$
entries111Construction of the right-hand side in Eq. (24) and recovery of the
velocity and pressure with Eq. (25) still requires multiplication with
matrices of size $\mathcal{O}(p^{6})$, but this is only necessary once per
solve.; explicit numbers for polynomial degrees up to $p=3$ are given in Table
5.

degree | number of unknowns per cell | cell-local matrix size
---|---|---
$p$ | $N_{p}$ | $+$ | $N_{u}$ | $=$ | $N_{\text{mixed}}$ | $N_{\lambda}$ | mixed | hybridised
$0$ | $1$ | $+$ | $3$ | $=$ | $4$ | $3$ | $7\times 7$ | $6\times 6$
$1$ | $8$ | $+$ | $24$ | $=$ | $32$ | $12$ | $44\times 44$ | $24\times 24$
$2$ | $27$ | $+$ | $81$ | $=$ | $108$ | $27$ | $135\times 135$ | $54\times 54$
$3$ | $64$ | $+$ | $192$ | $=$ | $256$ | $48$ | $304\times 304$ | $96\times 96$
Table 5: Number of unknowns per cell and local matrix sizes for different
polynomial degrees $p$, using $DG_{p}$ pressure elements and $RT_{p}$ velocity
elements [46]. The number of pressure unknowns is denoted with $N_{p}$, the
number of velocity unknowns is $N_{u}$ and the number of mixed unknowns is
$N_{\text{mixed}}=N_{p}+N_{u}$; $N_{\lambda}$ denotes the number of trace
unknowns. Shared unknowns associated with facets of the mesh are weighted with
a factor $1/2$ in $N_{u}$ and $N_{\lambda}$.

Assuming that the multigrid solver achieves $p$-independent convergence, we
therefore expect the hybridised solver to be more efficient, simply because
application of the matrix $S$ in Eq. (22) is cheaper than multiplication by
the $2\times 2$ matrix in Eq. (9).

###  4.2 Future work

There are several avenues for future research. Firstly, the discussion in the
previous section implies that the relative runtime of the three different
solver configurations studied here will likely change for higher order
discretisations. An increase in discretisation order will also open the door
for additional optimisations, such as the use of sum-factorisation techniques
in matrix-free implementations (see e.g. [47, 48]), which can have a
significant effect on solver performance. Secondly, the multigrid algorithm
for the second order Dirichlet problem in [32] uses a conforming, piecewise
linear P1 coarse level space whereas in this work we chose to use a piecewise
constant space for practical reasons. Extending our non-nested multigrid
algorithm to a P1 coarse level space could potentially be more efficient, but
requires substantial changes to the LFRic infrastructure. This change and a
theoretical analysis similar to [32] should be pursued in a subsequent
publication. Currently, the pressure operator in the non-nested multigrid
algorithm is constructed as an approximate Schur-complement of the original
saddle point system. An alternative approach, which was used in [34], would be
to form the Schur-complement of the continuous equations before spatial
discretisation and then re-discretise in pressure space. While previous work
[26] has shown that the bespoke multigrid algorithms for finite element
discretisations of atmospheric model problems are competitive with algebraic
multigrid (AMG), it is certainly worth re-visiting this in the context of the
hybridisable discretisations considered here. One could, for example, use AMG
as a coarse-level solver for Eq. (28) in the non-nested two-level method in
Algorithm 2 or directly as a solver for the system of Lagrange multipliers in
Eq. (22).

Given that time-to-solution is critical in operational forecast models, it is
also worth exploring further low-level performance improvements. As discussed
in Section 3.5, a significant fraction of the runtime is spent in tridiagonal
solves. When run in parallel, the Thomas algorithm is limited by memory
bandwidth, i.e. the time it takes to load the right-hand side and the three
matrix diagonals into memory and write the solution vector. It is likely that
the matrix columns change very little across the horizontal domain, and one
could therefore consider using the same, cached values over small patches of
the grid. Although this might lead to a slightly weaker multigrid smoother,
this is likely offset by a potentially twofold increase in computational
efficiency in the smoother due to reduced memory traffic. A similar
optimisation would be to divide all equations by their diagonal, which would
then only require loading the off-diagonal matrix entries in the smoother and
operator application. However, additional divisions and multiplications by the
diagonal will be required before and after intergrid transfer in the multigrid
algorithm.

## Acknowledgements

This work was funded by two EPSRC Impact Acceleration Awards (grant numbers
EP/R511547/1, EP/R51164X/1) and an EPSRC PostDoctoral Award (grant number
EP/W522491/1) of Matthew Griffith.

##  Appendix A Implementation in the LFRic Solver API

Figures 8, 9 and 10 show the LFRic Solver API implementation of the three
different solver configurations described in Section 3.2. Each iterative
solver instance (shown in green in the diagrams) contains operator objects
(blue) and preconditioner objects (red). Since the types for iterative
solvers, operators and preconditioners are derived from their respective
abstract base types, they conform to a well-defined interface and can be
plugged together in the code. In addition to common Krylov subspace methods
such as BiCGStab, it is also possible to only apply the preconditioner in a
specialised “pre-only” iterative solver type. This approach is used in the
hybridised solvers in Figures 9 and 10 since the pressure- and velocity
unknowns that are reconstructed from the trace field are the exact solutions
of the original, non-hybridised system in Eq. (9).

Figure 8: Solver API diagram of solver Configuration 1, the approximate Schur
complement solver described in Section 2.3. Figure 9: Solver API diagram of
solver Configuration 2, the hybridised solver described in Section 2.4 and
Section 2.5, using a few iterations of the linesmoother (Algorithm 1) to solve
the trace system in Eq. (23). Figure 10: Solver API diagram of solver
Configuration 3, the hybridised solver described in Section 2.4 and Section
2.5, using the non-nested multigrid method (Algorithm 2) to precondition the
trace system in Eq. (23). The coarse level pressure system is solved with the
multigrid method described in Section 2.3.

##  Appendix B Solver Parameters

The performance of the sophisticated solvers considered in this work depends
on a number of parameters such as the number of pre- and post- smoothing steps
and the relaxation factor $\omega$ in the Jacobi iteration in Algorithm 1. We
refer to the corresponding relaxation factor for the pressure multigrid
algorithm in [25] as $\omega_{\text{pressure}}$. To verify that the chosen
parameters are sensible, we investigated the runtime of the solver in
Configurations 1, 2 and 3c for a range of parameters; the code was run for the
same single node setup as used in Section 3.4. For Configuration 1 (pressure
multigrid) the number of pre- and post- smoothing iterations are set to 1 or 2
and the relaxation factor is varied between $\omega_{\text{pressure}}=0.3$ and
$\omega_{\text{pressure}}=1.0$. The time spent in the solver for each setting
is shown in Figure 11.

Figure 11: Average time spent in the solver for Configuration 1 (pressure
multigrid solver) as a function of the relaxation factor
$\omega_{\text{pressure}}$ for different numbers of pre- and post- smoothing
steps.

For $\omega_{\text{pressure}}\leq 0.9$ the measured times show very little
variation, indicating that performance depends only weakly on the relaxation
factor and the number of pre- and post- smoothing steps. The only exception is
the setup consisting of only a single pre- and post- smoothing iteration,
which is distinctly slower than the other configurations for smaller values of
$\omega_{\text{pressure}}$.

For Configuration 2 (hybridised single-level) the range of relaxation factors
$\omega$ which produce good performance is reduced to $0.4\leq\omega\leq 0.8$.
There is a tradeoff between the efficiency of the preconditioner, which
improves with the number of smoother iterations, and the cost per iteration in
the Krylov solver, which favours fewer smoothing steps. Two smoothing
iterations is generally the optimal choice. Figure 12 shows that the solver
cost gradually increases as more iterations are used.

Figure 12: Average time spent in the solver for Configuration 2 (hybridised
single-level solver) as a function of the relaxation factor $\omega$ for
different numbers of iterations of the Jacobi smoother that is used as a
preconditioner. For an even number of iterations the solver only converges for
$\omega\leq 0.8$; for three smoothing iterations the solver does not converge
for $\omega=0.3$.

Finally, the experiment is repeated for Configuration 3 (hybridised
multigrid). In this case the coarse level pressure multigrid preconditioner is
fixed to use 1 pre- and 2 post-smooth iterations with
$\omega_{\text{pressure}}=0.9$ and the number of smoothing iterations and the
relaxation factor in the nonested two-level preconditioner (Algorithm 2) is
varied. Figure 13 shows the time spent in the solver as a function of the
relaxation factor and for up to three pre-smoothing steps. As long as
$\omega\leq 0.7$ there is little difference in cost between the
configurations, with a slight improvement in performance as $\omega$
increases. Unless the number of post-smooth iterations is limited to one,
performance drops substantially for $\omega=0.8$ and the solver fails to
converge for larger relaxation factors. The optimal setup is generally to use
1 pre- and 2 post smooth iterations.

Figure 13: Average time spent in the solver for Configuration 3c (hybridised
multigrid solver) as a function of the relaxation factor $\omega$ for
different numbers of pre- and post- smoothing steps of the non-nested two-
level method. The solver does not converge for relaxation factors larger than
$\omega=0.8$ if the number of postsmoothing steps exceeds one.

## References

  * [1] Aimé Fournier, Mark A Taylor, and Joseph J Tribbia.  The spectral element atmosphere model (SEAM): High-resolution parallel computation and localized resolution of regional dynamics.  Monthly Weather Review, 132(3):726–748, 2004.
  * [2] Stephen J Thomas and Richard D Loft.  The NCAR spectral element climate dynamical core: semi-implicit Eulerian formulation.  Journal of Scientific Computing, 25(1):307–322, 2005.
  * [3] John M Dennis, Jim Edwards, Katherine J Evans, Oksana Guba, Peter H Lauritzen, Arthur A Mirin, Amik St-Cyr, Mark A Taylor, and Patrick H Worley.  CAM-SE: A scalable spectral element dynamical core for the Community Atmosphere Model.  The International Journal of High Performance Computing Applications, 26(1):74–89, 2012.
  * [4] James F Kelly and Francis X Giraldo.  Continuous and discontinuous Galerkin methods for a scalable three-dimensional nonhydrostatic atmospheric model: Limited-area mode.  Journal of Computational Physics, 231(24):7988–8008, 2012.
  * [5] Francis X Giraldo, James F Kelly, and Emil M Constantinescu.  Implicit-explicit formulations of a three-dimensional nonhydrostatic unified model of the atmosphere (NUMA).  SIAM Journal on Scientific Computing, 35(5):B1162–B1194, 2013.
  * [6] Lei Bao, Robert Klöfkorn, and Ramachandran D Nair.  Horizontally explicit and vertically implicit (HEVI) time discretization scheme for a discontinuous Galerkin nonhydrostatic model.  Monthly Weather Review, 143(3):972–990, 2015.
  * [7] Simone Marras, James F Kelly, Margarida Moragues, Andreas Müller, Michal A Kopera, Mariano Vázquez, Francis X Giraldo, Guillaume Houzeaux, and Oriol Jorba.  A review of element-based Galerkin methods for numerical weather prediction: Finite elements, spectral elements, and discontinuous Galerkin.  Archives of Computational Methods in Engineering, 23(4):673–722, 2016.
  * [8] Colin J. Cotter and Jemma Shipton.  Mixed finite elements for numerical weather prediction.  Journal of Computational Physics, 231(21):7076–7091, 2012.
  * [9] Andrew Staniforth, Thomas Melvin, and Colin Cotter.  Analysis of a mixed finite-element pair proposed for an atmospheric dynamical core.  Quarterly Journal of the Royal Meteorological Society, 139(674):1239–1254, 2013.
  * [10] Colin J. Cotter and John Thuburn.  A finite element exterior calculus framework for the rotating shallow-water equations.  Journal of Computational Physics, 257:1506–1526, 2014.
  * [11] Andrew TT McRae and Colin J Cotter.  Energy-and enstrophy-conserving schemes for the shallow-water equations, based on mimetic finite elements.  Quarterly Journal of the Royal Meteorological Society, 140(684):2223–2234, 2014.
  * [12] Andrea Natale, Jemma Shipton, and Colin J. Cotter.  Compatible finite element spaces for geophysical fluid dynamics.  Dynamics and Statistics of the Climate System, 1(1):dzw005, 2016\.
  * [13] Thomas H. Gibson, Andrew T.T. McRae, Colin J. Cotter, Lawrence Mitchell, and David A. Ham.  Compatible Finite Element Methods for Geophysical Flows: Automation and Implementation Using Firedrake.  Springer International Publishing, 2019.
  * [14] Douglas N Arnold, Richard S Falk, and Ragnar Winther.  Finite element exterior calculus, homological techniques, and applications.  Acta numerica, 15:1–155, 2006.
  * [15] Douglas Arnold, Richard Falk, and Ragnar Winther.  Finite element exterior calculus: from Hodge theory to numerical stability.  Bulletin of the American mathematical society, 47(2):281–354, 2010\.
  * [16] Douglas Arnold and Gerard Awanou.  Finite element differential forms on cubical meshes.  Mathematics of computation, 83(288):1551–1570, 2014.
  * [17] F. Brezzi, D. Boffi, L. Demkowicz, R. G. Durán, R. S. Falk, and M. Fortin.  Mixed finite elements, compatibility conditions, and applications.  Springer, 2008.
  * [18] Daniele Boffi, Franco Brezzi, Michel Fortin, et al.  Mixed finite element methods and applications, volume 44.  Springer, 2013.
  * [19] Akio Arakawa and Vivian R Lamb.  Computational design of the basic dynamical processes of the UCLA general circulation model.  Methods in computational physics, 17:173–265, 1977.
  * [20] Thomas Melvin, Tommaso Benacchio, Ben Shipway, Nigel Wood, John Thuburn, and Colin Cotter.  A mixed finite-element, finite-volume, semi-implicit discretization for atmospheric dynamics: Cartesian geometry.  Quarterly Journal of the Royal Meteorological Society, 145(724):2835–2853, 2019.
  * [21] Samantha V Adams, Rupert W Ford, M Hambley, JM Hobson, I Kavčič, Christopher M Maynard, Thomas Melvin, Eike Hermann Müller, S Mullerworth, AR Porter, et al.  LFRic: Meeting the challenges of scalability and performance portability in Weather and Climate models.  Journal of Parallel and Distributed Computing, 132:383–396, 2019\.
  * [22] Pierre-Arnaud Raviart and Jean-Marie Thomas.  A mixed finite element method for 2-nd order elliptic problems.  In Mathematical aspects of finite element methods, pages 292–315. Springer, 1977.
  * [23] D. N. Arnold, R. S. Falk, and R. Winther.  Multigrid in $H(\text{div})$ and $H(\text{curl})$.  Numerische Mathematik, 85:197–217, 2000.
  * [24] Ralf Hiptmair and Jinchao Xu.  Nodal auxiliary space preconditioning in $H(\text{curl})$ and $H(\text{div})$ spaces.  SIAM Journal on Numerical Analysis, 45(6):2483–2509, 2007.
  * [25] Christopher Maynard, Thomas Melvin, and Eike Hermann Müller.  Multigrid preconditioners for the mixed finite element dynamical core of the LFRic atmospheric model.  Quarterly Journal of the Royal Meteorological Society, 146(733):3917–3936, 2020.
  * [26] Lawrence Mitchell and Eike Hermann Müller.  High level implementation of geometric multigrid solvers for finite element problems: applications in atmospheric modelling.  Journal of Computational Physics, 327:1–18, 2016.
  * [27] Steffen Börm and Ralf Hiptmair.  Analysis of tensor product multigrid.  Numerical Algorithms, 26(3):219–234, 2001.
  * [28] Douglas N. Arnold and Franco Brezzi.  Mixed and nonconforming finite element methods: implementation, postprocessing and error estimates.  ESAIM: Mathematical Modelling and Numerical Analysis, 19(1):7–32, 1985.
  * [29] Bernardo Cockburn and Jayadeep Gopalakrishnan.  A characterization of hybridized mixed methods for second order elliptic problems.  SIAM Journal on Numerical Analysis, 42(1):283–301, 2004.
  * [30] James H Bramble, Joseph E Pasciak, and Jinchao Xu.  The analysis of multigrid algorithms with nonnested spaces or noninherited quadratic forms.  Mathematics of Computation, 56(193):1–34, 1991.
  * [31] Jayadeep Gopalakrishnan and Shuguang Tan.  A convergent multigrid cycle for the hybridized mixed method.  Numerical Linear Algebra with Applications, 16(9):689–714, 2009\.
  * [32] Bernardo Cockburn, Olivier Dubois, Jay Gopalakrishnan, and Shuguang Tan.  Multigrid for an HDG method.  IMA Journal of Numerical Analysis, 34(4):1386–1425, 2014.
  * [33] Shinhoo Kang, Francis X Giraldo, and Tan Bui-Thanh.  IMEX HDG-DG: A coupled implicit hybridized discontinuous Galerkin and explicit discontinuous Galerkin approach for shallow water systems.  Journal of Computational Physics, 401:109010, 2020.
  * [34] Jack Betteridge, Thomas H Gibson, Ivan G Graham, and Eike H Müller.  Multigrid preconditioners for the hybridised discontinuous galerkin discretisation of the shallow water equations.  Journal of Computational Physics, 426:109948, 2021.
  * [35] D. Pardo, J. Álvarez-Aramberri, M. Paszynski, L. Dalcin, and V.M. Calo.  Impact of element-level static condensation on iterative solver performance.  Computers & Mathematics with Applications, 70(10):2331–2341, nov 2015.
  * [36] Nigel Wood, Andrew Staniforth, Andy White, Thomas Allen, Michail Diamantakis, Markus Gross, Thomas Melvin, Chris Smith, Simon Vosper, Mohamed Zerroukat, et al.  An inherently mass-conserving semi-implicit semi-Lagrangian discretization of the deep-atmosphere global non-hydrostatic equations.  Quarterly Journal of the Royal Meteorological Society, 140(682):1505–1520, 2014.
  * [37] Pierre-Arnaud Raviart and Jean-Marie Thomas.  A mixed finite element method for 2-nd order elliptic problems.  In Mathematical aspects of finite element methods, pages 292–315. Springer, 1977.
  * [38] William H Press, Saul A Teukolsky, William T Vetterling, and Brian P Flannery.  Numerical recipes 3rd edition: The art of scientific computing.  Cambridge university press, 2007.
  * [39] Satish Balay, Shrirang Abhyankar, Mark F. Adams, Steven Benson, Jed Brown, Peter Brune, Kris Buschelman, Emil Constantinescu, Lisandro Dalcin, Alp Dener, Victor Eijkhout, William D. Gropp, Václav Hapla, Tobin Isaac, Pierre Jolivet, Dmitry Karpeev, Dinesh Kaushik, Matthew G. Knepley, Fande Kong, Scott Kruger, Dave A. May, Lois Curfman McInnes, Richard Tran Mills, Lawrence Mitchell, Todd Munson, Jose E. Roman, Karl Rupp, Patrick Sanan, Jason Sarich, Barry F. Smith, Stefano Zampini, Hong Zhang, Hong Zhang, and Junchao Zhang.  PETSc/TAO users manual.  Technical Report ANL-21/39 - Revision 3.16, Argonne National Laboratory, 2021.
  * [40] Satish Balay, Shrirang Abhyankar, Mark F. Adams, Steven Benson, Jed Brown, Peter Brune, Kris Buschelman, Emil M. Constantinescu, Lisandro Dalcin, Alp Dener, Victor Eijkhout, William D. Gropp, Václav Hapla, Tobin Isaac, Pierre Jolivet, Dmitry Karpeev, Dinesh Kaushik, Matthew G. Knepley, Fande Kong, Scott Kruger, Dave A. May, Lois Curfman McInnes, Richard Tran Mills, Lawrence Mitchell, Todd Munson, Jose E. Roman, Karl Rupp, Patrick Sanan, Jason Sarich, Barry F. Smith, Stefano Zampini, Hong Zhang, Hong Zhang, and Junchao Zhang.  PETSc Web page.  https://petsc.org/, 2021.
  * [41] Satish Balay, William D. Gropp, Lois Curfman McInnes, and Barry F. Smith.  Efficient management of parallelism in object oriented numerical software libraries.  In E. Arge, A. M. Bruaset, and H. P. Langtangen, editors, Modern Software Tools in Scientific Computing, pages 163–202. Birkhäuser Press, 1997.
  * [42] Markus Blatt and Peter Bastian.  The iterative solver template library.  In International Workshop on Applied Parallel Computing, pages 666–675. Springer, 2006.
  * [43] Henk A Van der Vorst.  Bi-CGSTAB: A fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems.  SIAM Journal on scientific and Statistical Computing, 13(2):631–644, 1992.
  * [44] Yousef Saad.  Iterative methods for sparse linear systems.  SIAM, 2003.
  * [45] Steven Sandbach, John Thuburn, Danail Vassilev, and Michael G Duda.  A semi-implicit version of the MPAS-atmosphere dynamical core.  Monthly Weather Review, 143(9):3838–3855, 2015.
  * [46] Douglas N Arnold and Anders Logg.  Periodic table of the finite elements.  Siam News, 47(9):212, 2014.
  * [47] Peter EJ Vos, Spencer J Sherwin, and Robert M Kirby.  From h to p efficiently: Implementing finite and spectral/hp element methods to achieve optimal performance for low-and high-order discretisations.  Journal of Computational Physics, 229(13):5161–5181, 2010.
  * [48] Martin Kronbichler and Katharina Kormann.  A generic interface for parallel cell-based finite element operator application.  Computers & Fluids, 63:135–147, 2012.