The dynamics of geophysical fluids is characterised by phenomena which occur at very different time scales. While large-scale flow in the atmosphere is described by the slow evolution of quasi-geostrophic modes, fast acoustic waves play an important rôle in incorporating non-linear dynamics and rapid adjustment to geostrophic balance. However, acoustic waves carry very little energy and hence do not have to be represented to high accuracy in computer simulations. Efficient numerical models for weather- and climate-prediction often exploit this scale-separation by using semi-implicit time integrators. Treating the fast waves implicitly (and thereby removing any instabilities due to acoustic modes) avoids tight restrictions on the time step size. It is sufficient to chose the timestep such that it allows the stable and accurate representation of the slow modes, which are treated explicitly. Examples of popular time-integrators are semi-implicit semi-Lagrangian methods used by the UK Met Office  and implicit-explicit (IMEX) methods [2, 3, 4, 5] (see also  for applications in atmospheric modelling).
In a computer model the equations of fluid-dynamics have to be represented on a grid with finite resolution. While simple finite-difference- or finite-volumes schemes on structured meshes (such as latitude-longitude grids) work well for local area models, they suffer from the convergence of grid lines at the poles on global grids . This leads to artificially tight constraints on the timestep size due to the CFL condition because of the small grid cells at the poles. Furthermore, all-to-all communications limit scalability in massively parallel implementations. Finite element methods on quasi-uniform unstructured global grids overcome this problem, but care has to be taken when solving the resulting linear systems which arise in semi-implicit timestepping. The standard approach uses a Schur-complement reduction to a pressure system. Naively this results in a dense Schur-complement operator in pressure space. For conforming finite element methods this can be traced back to the fact that the mass-matrices are non-diagonal. In this case the problem is easily overcome by under-integrating and co-locating the quadrature points with nodal basis points, or forming an approximate Schur-complement .
Discontinuous Galerkin (DG) discretisations have recently received significant attention in atmospheric modelling 
. DG methods have several advantages, especially on modern chip architectures: since the unknowns are associated with cells of the mesh, data access can be coalesced, which is particularly important for chips with wide vector units. For sufficiently smooth solutions, higher order methods are numerically efficient since errors decrease exponentially as the polynomial degree of the DG approximation increases. Matrix-free implementations are particularly efficient on chips with limited memory bandwidth since they have a very high computational intensity; sum factorisation techniques can further reduce the computational complexity. Unfortunately, for DG methods the numerical flux introduces off-diagonal matrix elements which cause issues for semi-implicit discretisations. The flux generates artificial diffusion terms, which are numerical artifacts and can not be removed in an obvious way. This rules out the naive Schur-complement approach for solving the implicit linear system.
In  a solution to this problem is presented for the shallow water equations (SWEs), a popular two-dimensional model system for atmospheric dynamics. The authors use the hybridized Discontinuous Galerkin (HDG) method (see e.g. [12, 13, 14, 15]), which is adapted for the SWEs in [16, 17]. Introducing a set of flux-unknowns on the facets of the mesh allows the reduction to a sparse Schur-complement system in the flux space. In  this problem is assembled and solved with a direct method via LU factorisation. While this works for small problems, the cost of this direct solver grows rapidly as the problem size increases, even if regularities in the sparsity pattern are exploited (see  for a general overview of direct methods). This approach is therefore clearly unfeasible for large scale applications, especially if they are eventually to be used in time critical operational forecasting systems.
In this paper we present an alternative solution which overcomes those issues and we demonstrate the computational performance of our approach. We introduce a new non-nested multigrid method for solving the linear systems arising in the IMEX approach for the SWEs, when discretised using a hybridized (high order) DG method. The method is strongly influenced by the work of Gopalakrishnan and Tan  and Cockburn, Dubois, Gopalakrishnan and Tan . In particular  introduces a non-nested multigrid V-cycle for the hybridised mixed method for standard second order elliptic PDEs discretised using Raviart Thomas elements for the velocity approximation. This has been extended in  to the multigrid solution of the HDG discretisation for the same elliptic PDE. In this paper we adapt this approach to solve the linear and non-linear shallow water equations. The non-linear SWEs present additional challenges since they require the construction of a Lax-Friedrich flux, which results in a vector-valued flux space  for which the method in  is not directly applicable.
More specifically, the new contributions of this paper are:
The non-nested multigrid method of  is extended to the HDG solution of the linearised SWEs with IMEX methods.
To solve the non-linear SWEs by this approach we construct novel intergrid- and coarse-level operators for vector-valued flux spaces.
Based on this, the superior performance of the non-nested multigrid approach for both the linear- and non-linear SWEs is demonstrated for a parallel implementation.
Our approach differs from the more algebraic method presented in [23, 24], which constructs a multigrid hierarchy directly on the trace-space defined on the skeleton of the mesh; this is achieved by using -refinement and statically condensing out unknowns on fine edges which are not represented on the coarser mesh.
This paper is organised as follows: In Section 2 we briefly review the physics of the shallow water equations. The hybridized DG discretisation is described in Section 3 and the non-nested multigrid preconditioner for the resulting linear systems in IMEX timestepping is presented in Section 4. The specific integrators considered in this paper are introduced in Section 5 with focus on the linear solvers used for implicit methods. The implementation in the Firedrake library is described in Section 6 and results are presented in Section 7. We finally conclude and outline directions for future work in Section 8. Some more technical issues such as the explicit expressions for key bilinear forms are relegated to the appendices.
2 Shallow water equations
The shallow water equations (SWEs) describe the evolution of a layer of fluid under gravitational and inertial forces. In a rotating frame a Coriolis term has to be included as well. The SWEs can be written in conservative form for the geopotential height and the column-integrated momentum , which are defined for all points in a bounded domain . The time evolution of and is described by the following two equations:
Suitable boundary conditions have to be included at the boundary of . The geopotential height is defined as where denotes the perturbation of the fluid relative to its state at rest and is the gravitational acceleration. Furthermore where is the bathymetry of the sea-floor, measured from the surface at rest (see Fig. 1). The momentum is related to the column-averaged velocity of the fluid at each point by . is the latitude-dependent Coriolis parameter, and the vector is obtained by rotating by 90 degrees in the counter-clockwise direction.
Eq. (1) has been non-dimensionalised as follows: let be a reference height, and define . Further assume that any horizontal distances are measured in units of some reference length scale and times in units of a reference time scale . , are given in units of and the momentum in units of . The dimensionless reference gravity wave speed which multiplies the flux terms in Eq. (1) is then given by
For convenience, we introduce the variable , which allows us to write Eq. (1) in compact form as
The SWEs can be linearised about (this is a good approximation if at every point in the domain) to obtain
It is important to note that in a non-rotating frame () fast gravity waves are supported by Eq. (3). This can be seen by taking the divergence of the second equation and inserting it into the time derivative of the first, to obtain the wave equation
The (local) speed of the gravity waves is given by . More generally, in a rotating frame the fundamental solutions of the linear SWEs in Eq. (3) are inertia-gravity (Poincaré) waves. In this case the frequency is related to the wave vector by the dispersion relation . In analogy to Eq. (2), the linearised SWEs in Eq. (3) can be written as
with and the linear operator given by
3 IMEX HDG methods for the shallow water equations
3.1 Spatial discretisation
To discretise the SWEs with the discontinuous Galerkin (DG) method, we first construct a mesh of non-overlapping cells that partition the domain such that . The mesh spacing is defined as . Let be the set of facets of the mesh; is often called the “skeleton of ”. We further assume that the mesh consists of triangular or quadrilateral cells and has no hanging nodes. More specifically, each facet is a straight line which forms the side of exactly two triangles or quadrilaterals. For simplicity we use periodic boundary conditions in this paper, which implies that all facets are internal.
Since functions in the DG space which we introduce below are not continuous across cell-boundaries, care has to be taken when defining quantities on facets. For two neighbouring cells and that share a facet , let be the outward unit normal on the boundary of element , and the corresponding outward normal of cell (see Fig. 2), such that on the shared facet . On each facet , the following average , difference and jump operators can be defined for scalar () and vector valued () quantities:
At any point on the facet and for any scalar- or vector-valued piecewise polynomial functions , we define and . Note that the jump of a scalar quantity is vector-valued and vice versa.
Write for the space of polynomials of degree at most on an arbitrary domain . The (scalar-valued) DG space of order on the mesh is defined as
Case 1: Linear SWEs.
It is instructive to first discuss the DG discretisation of the linear SWEs in Eq. (4). This can be found by assuming that111In a slight abuse of notation, we use the same variable names for the continuous fields , , , etc. in Section 2 and for the DG fields discussed from this section onwards. , then multiplying by a test function and integrating by parts in each cell to obtain
Here is the outward normal on each facet of cell , with . Round and angular brackets denote scalar products in function spaces defined on the cell and its boundary respectively:
The numerical flux appearing in Eq. (6) (still to be defined) will depend on the values , of the field on both sides of the facets comprising and will couple the unknowns between neighbouring cells. For the linear problem both the local Lax-Friedrichs (Rusanov) flux  and the upwind flux (see e.g. ) are candidates for the numerical flux . As shown in [16, 17], those are defined as
where the matrix defined on is given by
If Eq. (6) is integrated in time with a semi-implicit method, a coupled sparse linear system for the fields has to be solved in each timestep. The standard approach in many atmospheric models is to eliminate the momentum unknowns from the discrete equations and solve the resulting elliptic Schur-complement system for the height perturbation , followed by reconstruction of the momentum field (see e.g.  for the equivalent procedure applied to the full Navier Stokes equations). Unfortunately, for the DG discretisation considered here, the numerical flux in Eq. (6) introduces artificial diffusion terms in the momentum equation. Although those terms disappear in the limit , they result in a dense Schur-complement operator which renders this naive Schur-complement approach unfeasible here. It is possible to address this issue by forming an approximate Schur-complement and using the resulting approximate solver as a preconditioner for a Krylov subspace method. A similar approach is pursued to deal with the non-diagonal velocity mass matrix in conforming mixed finite element discretisations in [8, 27]. However, this requires additional iterations over the large -system. As our numerical results in Section 7 demonstrate, this approximate Schur-complement method (which essentially drops the artificial diffusion terms in the preconditioner) is not competitive with the hybridized DG method pursued in this paper since the latter allows the formulation of an exact Schur-complement in a suitably defined finite element space on the facets.
3.2 Hybridised DG discretisation
The key idea of the hybridized DG method described in [16, 17] is to introduce a new variable , which is associated with the skeleton of the mesh, see Fig. 2 (right). As we will show below, it is possible to form an exact Schur-complement in . For this, recall that the facets are straight lines bounding the cells. On these we define the piecewise polynomial space
This allows the introduction of a flux which only depends on the value inside the cell and the field on its boundary. To see this, recall that for hyperbolic conservation laws the numerical upwind flux can be constructed by solving a Riemann problem on the interface between two cells (see e.g. section 2.4 of ). The solution contains both the field inside the cell and an intermediate state , which is usually eliminated with a jump condition. Since this jump condition involves the value of fields on both sides of the interface, the resulting flux in Eq. (7b) depends on and . As argued in [16, 17], the HDG method represents the intermediate state as in the space instead of eliminating ; the jump condition is enforced weakly.
replaces the numerical flux in Eq. (6) and is given by either (Lax-Friedrichs) or (upwind) defined on the boundary of each cell as
To close the system of equations, continuity of the flux is enforced weakly by requiring that on each facet
To simplify notation, introduce the following abbreviations:
As has been shown in , the solution of the original, non-hybridized system in Eq. (12) is identical to the solution of the HDG problem given in Eq. (13). In general, the additional field introduced in Eq. (13) on the facets has three components. As will be discussed in Section 4.1 below, can be eliminated for the upwind flux, whereas does not enter the equations if the Lax-Friedrichs method is used.
Case 2: non-linear SWEs.
Similarly, the weak form of the non-linear SWEs in Eq. (2) can be written for each cell as
where is the non-linear Lax-Friedrichs flux
with . Further define
By splitting the right hand side into two terms, we have isolated the fast dynamics due to gravity waves in , whereas any slower modes are described by . This observation is crucial for the construction of suitable semi-implicit time integrators, which avoid numerical instabilities due to the fast gravity waves.
3.3 Time discretisation
To integrate Eq. (18) (Eq. (13) can be dealt with in exactly the same way) in time, the term is integrated explicitly, whereas is treated implicitly. Denoting the solution at time by , a simple scheme would for example be the following “-method”:
with and subject to the flux condition Eq. (10). The equivalent expression for the linear SWEs in Eq. (13) can be obtained by replacing in Eq. (19). At each time step, calculating the field requires the solution of the linear system
Note that for a purely explicit method (), the system matrix arising from the left hand side of the system in Eq. (20) is simply the mass matrix (obviously, in this case Eq. (19) reduces to the explicit Euler integrator). On the other hand, -stage IMEX methods [2, 3, 4, 5, 6] are a generalisation of the method in Eq. (19). They require the construction of intermediate states which are obtained by solving the system of linear equations
The state at the next timestep satisfies
and is obtained by an additional mass-solve. Each particular IMEX method is determined by a choice of coefficients , , , which are commonly known as the Butcher tableau coefficients for the explicit and implicit terms. The equations in Eq. (21) can be solved stage-by-stage. At every stage a linear system of the abstract general form
with some positive constant has to be solved222
Note that the “-method” in Eq. (19) could be written as a 2-stage IMEX method with
and for all , .
4 Schur-complement multigrid solver
Constructing efficient solvers for Eq. (23) is the topic of this paper. In the following section we give explicit expressions for for the Lax-Friedrichs and upwind flux. We then describe a Schur complement approach which reduces the problem to solving an elliptic system in the flux variable on the facets. This flux system is solved with a new non-nested multigrid preconditioner based on ideas in .
4.1 Linear equation
Recall that in the domain, on the facets, and write the corresponding test-functions as , . Using the definitions in Eqs. (7a) and (7b), the bilinear form on the left hand side of Eq. (23) can be written down explicitly for the upwind- and Lax-Friedrichs flux. As shown in , the resulting bilinear form for the upwind flux does not depend on . For the Lax-Friedrichs flux the form does not depend on . Explicit expressions for those forms are given as follows:
The system in Eq. (23) can be solved with a Schur-complement approach by eliminating the field defined on the elements and leaving a problem for defined on the facets. For this write
where is the Riesz representer of . The scalar products in Eq. (26) are defined as
and the linear operators
4.2 Schur complement solve
Let the inverse of be defined in the obvious way as
It can now be shown that Eq. (23) can be solved in three steps.
- Step 1.
Given in Eq. (26), calculate the modified right hand side
- Step 2.
Solve the Schur-complement system
for the flux unknown .
- Step 3.
Reconstruct the solution as
The crucial observation is that both and only couple unknowns within one grid cell and therefore can both be represented by a set of cell-local matrices. Since further and couple the unknowns in one grid cell to the flux unknowns on its boundary, the matrix representation of the operator in Eq. (29) is sparse. Recall further that for the Lax-Friedrichs flux only and appear in Eq. (23) whereas for the upwind flux, the equation only contains and . Therefore Eq. (29) reduces to
4.3 Relationship to previous work in 
For the HDG upwind flux, the Schur-complement system in Eq. (29) is very similar to the elliptic problem written down in Theorem 2.1 of . To see this, observe that for constant , Eq. (2.5) in  can be written in our notation as
for some right hand side . The bilinear form is given by
where is a constant. If we assume spatially non-varying bathymetry () for the shallow water problem and set , , in Eq. (32), the resulting bilinear form differs from the one given in our Eq. (24) only by a zero order term, namely: