Slate: extending Firedrake's domain-specific abstraction to hybridized solvers for geoscience and beyond

Within the finite element community, discontinuous Galerkin (DG) and mixed finite element methods have become increasingly popular in simulating geophysical flows. However, robust and efficient solvers for the resulting saddle point and elliptic systems arising from these discretizations continue to be an ongoing challenge. One possible approach for addressing this issue is to employ a method known as hybridization, where the discrete equations are transformed such that classic static condensation and local post-processing methods can be employed. However, it is challenging to implement hybridization as performant parallel code within complex models whilst maintaining a separation of concerns between applications scientists and software experts. In this paper, we introduce a domain-specific abstraction within the Firedrake finite element library that permits the rapid execution of these hybridization techniques within a code-generating framework. The resulting framework composes naturally with Firedrake's solver environment, allowing for the implementation of hybridization and static condensation as runtime-configurable preconditioners via the Python interface to the Portable, Extensible Toolkit for Scientific Computation (PETSc), petsc4py. We provide examples derived from second-order elliptic problems and geophysical fluid dynamics. In addition, we demonstrate that hybridization shows great promise for improving the performance of solvers for mixed finite element discretizations of equations related to large-scale geophysical flows.

1 Introduction

The development of simulation software is an increasingly important aspect of modern scientific computing, in the geosciences in particular. Such software requires a vast range of knowledge spanning several disciplines, ranging from applications expertise to mathematical analysis to high-performance computing and low-level code optimization. Software projects developing automatic code generation systems have become quite popular in recent years, as such systems help create a separation of concerns which focuses on a particular complexity independent from the rest. This allows for agile collaboration between computer scientists with hardware and software expertise, computational scientists with numerical algorithm expertise, and domain scientists such as meteorologists, oceanographers and climate scientists. Examples of such projects in the domain of finite element methods include FreeFEM++ (Hecht, 2012), Sundance (Long et al., 2010), the FEniCS Project (Logg et al., 2012a), Feel++ (Prud'Homme et al., 2012), and Firedrake (Rathgeber et al., 2016).

The finite element method (FEM) is a mathematically robust framework for computing numerical solutions of partial differential equations (PDEs) that has become increasingly popular in fluids and solids models across the geosciences, with a formulation that is highly amenable to code generation techniques. A description of the weak formulation of the PDEs, together with appropriate discrete function spaces, is enough to characterize the finite element problem. Both the FEniCS and Firedrake projects employ the Unified Form Language (UFL) (Alnæs et al., 2014) to specify the finite element integral forms and discrete spaces necessary to properly define the finite element problem. UFL is a highly expressive domain-specific language (DSL) embedded in Python, which provides the necessary abstractions for code generation systems.

There are classes of finite element discretizations resulting in discrete systems that can be solved more efficiently by directly manipulating local tensors. For example, the static condensation technique for the reduction of global finite element systems (Guyan, 1965; Irons, 1965) produces smaller globally coupled linear systems by eliminating interior unknowns to arrive at an equation for the degrees of freedom defined on cell interfaces only. This procedure is analogous to the point-wise elimination of variables used in staggered finite difference codes, such as the ENDGame dynamical core (Melvin et al., 2010; Wood et al., 2014) of the UK Meteorological Office (Met Office) but requires the local inversion of finite element systems. For finite element discretizations of coupled PDEs, the hybridization technique provides a mechanism for enabling the static condensation of more complex linear systems. First introduced by Fraeijs de Veubeke (1965) and analyzed further by Brezzi and Fortin (1991), Cockburn et al. (2009a), and Boffi et al. (2013), the hybridization method introduces Lagrange multipliers enforcing certain continuity constraints. Local static condensation can then be applied to the augmented system to produce a reduced equation for the multipliers. Methods of this type are often accompanied by local post-processing techniques, which exploit the approximation properties of the Lagrange multipliers. This enables the manufacturing of fields exhibiting superconvergent phenomena or enhanced conservation properties (Arnold and Brezzi, 1985; Brezzi et al., 1985; Bramble and Xu, 1989; Stenberg, 1991; Cockburn et al., 2009b, 2010b). These procedures require invasive manual intervention during the equation assembly process in intricate numerical code.

In this paper, we provide a simple yet effective high-level abstraction for localized dense linear algebra on systems derived from finite element problems. Using embedded DSL technology, we provide a means to enable the rapid development of hybridization and static condensation techniques within an automatic code generation framework. In other words, the main contribution of this paper is in solving the problem of automatically translating from the mathematics of static condensation and hybridization to compiled code. This automated translation facilitates the separation of concerns between applications scientists and computational/computer scientists and facilitates the automated optimization of compiled code. This framework provides an environment for the development and testing of numerics relevant to the Gung-Ho Project, an initiative by the UK Met Office in designing the next-generation atmospheric dynamical core using mixed finite element methods (Melvin et al., 2019). Our work is implemented in the Firedrake finite element library and the PETSc solver library (Balay et al., 1997, 2019), accessed via the Python interface petsc4py (Dalcin et al., 2011).

The rest of the paper is organized as follows. We introduce common notation used throughout the paper in Sect. 1.1. The embedded DSL, called “Slate”, is introduced in Sect. 2, which allows concise expression of localized linear algebra operations on finite element tensors. We provide some contextual examples for static condensation and hybridization in Sect. 3, including a discussion on post-processing. We then outline in Sect. 4 how, by interpreting static condensation techniques as a preconditioner, we can go further and automate many of the symbolic manipulations necessary for hybridization and static condensation. We first demonstrate our implementation on a manufactured problem derived from a second-order elliptic equation, starting in Sect. 5. The first example compares a hybridizable discontinuous Galerkin (HDG) method with an optimized continuous Galerkin method. Section 5.2 illustrates the composability and relative performance of hybridization for compatible mixed methods applied to a semi-implicit discretization of the nonlinear rotating shallow water equations. Our final example in Sect. 5.3 demonstrates time-step robustness of a hybridizable solver for a compatible finite element discretization of a rotating linear Boussinesq model. Conclusions follow in Sect. 6.

1.1 Notation

We begin by establishing notation used throughout this paper. Let 𝒯_h denote a tessellation of Ω⊂ℝⁿ, the computational domain, consisting of polygonal elements K associated with a mesh size parameter h, and $\partial {\mathcal{T}}_h=\mathit{\{}e\in \partial K:K\in {\mathcal{T}}_h\mathit{\}}$ the set of facets of 𝒯_h. The set of facets interior to the domain Ω is denoted by ${\mathcal{E}}_h^{{}^{\circ }}:=\partial {\mathcal{T}}_h\setminus \partial \mathrm{\Omega }$. Similarly, we denote the set of exterior facets as ${\mathcal{E}}_h^{\partial }:=\partial {\mathcal{T}}_h\cap \partial \mathrm{\Omega }$. For brevity, we denote the finite element integral forms over 𝒯_h and any facet set $\mathrm{\Gamma }\subset \partial {\mathcal{T}}_h$ by

$$\begin{array}{}\text{(1)}& {\displaystyle }{\left(u,v\right)}_K& {\displaystyle }=\underset{K}{\int }u\cdot v\mathrm{d}x,& {\displaystyle }{〈u,v〉}_e& {\displaystyle }=\underset{e}{\int }u\cdot v\mathrm{d}s,\text{(2)}& {\displaystyle }{\left(u,v\right)}_{{\mathcal{T}}_h}& {\displaystyle }=\sum _{K\in {\mathcal{T}}_h}{\left(u,v\right)}_K,& {\displaystyle }{〈u,v〉}_{\mathrm{\Gamma }}& {\displaystyle }=\sum _{e\in \mathrm{\Gamma }}{〈u,v〉}_e,\end{array}$$

where dx and ds denote appropriate integration measures. The operation ⋅ should be interpreted as standard multiplication for scalar functions or a dot product for vector functions.

For any double-valued vector field w on a facet $e\in \partial {\mathcal{T}}_h$, we define the jump of its normal component across e by

$$\begin{array}{}\text{(3)}& {\left[\left[\mathbf{w}\right]\right]}_e=\left\{\begin{array}{ll}\mathbf{w}{|}_{e^{+}}\cdot {\mathbf{n}}_{e^{+}}+\mathbf{w}{|}_{e^{-}}\cdot {\mathbf{n}}_{e^{-}},& e\in {\mathcal{E}}_h^{{}^{\circ }}\\ \mathbf{w}{|}_e\cdot {\mathbf{n}}_e,& e\in {\mathcal{E}}_h^{\partial }\end{array}\right.,\end{array}$$

where + and − denote arbitrarily but globally defined sides of the facet. Here, ${\mathbf{n}}_{e^{+}}$ and ${\mathbf{n}}_{e^{-}}$ are the unit normal vectors with respect to the positive and negative sides of the facet e. Whenever the facet domain is clear by the context, we omit the subscripts for brevity and simply write [⋅].

2 A system for localized algebra on finite element tensors

We present an expressive language for dense linear algebra on the elemental matrix systems arising from finite element problems. The language, which we call Slate, provides typical mathematical operations performed on matrices and vectors; hence the input syntax is comparable to high-level linear algebra software such as MATLAB. The Slate language provides basic abstract building blocks which can be used by a specialized compiler for linear algebra to generate low-level code implementations.

Slate is heavily influenced by the Unified Form Language (UFL) (Alnæs et al., 2014; Logg et al., 2012a), a DSL embedded in Python which provides symbolic representations of finite element forms. The expressions can be compiled by a form compiler, which translates UFL into low-level code for the local assembly of a form over the cells and facets of a mesh. In a similar manner, Slate expressions are compiled to low-level code that performs the requested linear algebra element-wise on a mesh.

2.1 An overview of Slate

To clarify conventions and the scope of Slate, we start by establishing our notation for a general finite element form following the convention of Alnæs et al. (2014). We define a real-valued multi-linear form as an operator which maps a list of arguments $\mathbf{v}=(v_{\mathrm{0}},\mathrm{\cdots },v_{\mathit{\alpha }-\mathrm{1}})\in V_{\mathrm{0}}\times \mathrm{\cdots }\times V_{\mathit{\alpha }-\mathrm{1}}$ into ℝ:

$$\begin{array}{}\text{(4)}& a:V_{\mathrm{0}}\times \mathrm{\cdots }\times V_{\mathit{\alpha }-\mathrm{1}}\to \mathbb{R},a\mapsto a(v_{\mathrm{0}},\mathrm{\cdots },v_{\mathit{\alpha }-\mathrm{1}})=a\left(\mathbf{v}\right),\end{array}$$

where a is linear in each argument v_k. The arity of a form is α, an integer denoting the total number of form arguments. In traditional finite element nomenclature (for α≤2), V₀ is referred to as the space of test functions and V₁ as the space of trial functions. Each V_k are referred to as argument spaces. Forms with arity $\mathit{\alpha }=\mathrm{0},\mathrm{1}$, or 2 are best interpreted as the more familiar mathematical objects: scalars (0-forms), linear forms or functionals (1-forms), and bilinear forms (2-forms), respectively.

If a given form a is parameterized by one or more coefficients, say $\mathbf{c}=(c_{\mathrm{0}},\mathrm{\cdots },c_q)\in C_{\mathrm{0}}\times \mathrm{\cdots }\times C_q$, where $\mathit{\{}{C}_k{\mathit{\}}}_{k=\mathrm{0}}^q$ are coefficient spaces, then we write

$$\begin{array}{}\text{(5)}& \begin{array}{rl}& a:C_{\mathrm{0}}\times \mathrm{\cdots }\times C_q\times V_{\mathrm{0}}\times \mathrm{\cdots }\times V_{\mathit{\alpha }-\mathrm{1}}\to \mathbb{R},\\ & a\mapsto a(c_{\mathrm{0}},\mathrm{\cdots },c_q;v_{\mathrm{0}},\mathrm{\cdots },v_{\mathit{\alpha }-\mathrm{1}})=a(\mathbf{c};\mathbf{v}).\end{array}\end{array}$$

From here on, we shall work exclusively with forms that are linear in v and possibly nonlinear in the coefficients c. This is reasonable since nonlinear methods based on Newton iterations produce linear problems via Gâteaux differentiation of a nonlinear form corresponding to a PDE (also known as the form Jacobian). We refer the interested reader to Alnæs et al. (2014, Sect. 2.1.2) for more details. For clarity, we present examples of multi-linear forms of arity $\mathit{\alpha }=\mathrm{0},\mathrm{1}$, and 2 that frequently appear in finite element discretizations:

$$\begin{array}{}\text{(6)}& \begin{array}{rl}& a(\mathit{\kappa };v,u):={\left(\mathrm{\nabla }v,\mathit{\kappa }\mathrm{\nabla }u\right)}_{{\mathcal{T}}_h}\equiv \sum _{K\in {\mathcal{T}}_h}\underset{K}{\int }\mathrm{\nabla }v\cdot \left(\mathit{\kappa }\mathrm{\nabla }u\right)\mathrm{d}x,\\ & \mathit{\kappa }\in C_{\mathrm{0}},u\in V_{\mathrm{1}},v\in V_{\mathrm{0}},\mathit{\alpha }=\mathrm{2},q=\mathrm{1},\end{array}\end{array}$$

$$\begin{array}{}\text{(7)}& \begin{array}{rl}& a(f;v):={\left(v,f\right)}_{{\mathcal{T}}_h}\equiv \sum _{K\in {\mathcal{T}}_h}\underset{K}{\int }vf\mathrm{d}x,\\ & f\in C_{\mathrm{0}},v\in V_{\mathrm{0}},\mathit{\alpha }=\mathrm{1},q=\mathrm{1},\end{array}\end{array}$$

$$\begin{array}{}\text{(8)}& \begin{array}{rl}& a(f,g;):={\left(f-g,f-g\right)}_{{\mathcal{T}}_h}\equiv \sum _{K\in {\mathcal{T}}_h}\underset{K}{\int }|f-g{|}^{\mathrm{2}}\mathrm{d}x,\\ & g\in C_{\mathrm{1}},f\in C_{\mathrm{0}},\mathit{\alpha }=\mathrm{0},q=\mathrm{2},\end{array}\end{array}$$

$$\begin{array}{}\text{(9)}& \begin{array}{rl}& a(\mathit{\gamma },\mathit{\sigma }):={〈\mathit{\gamma },\left[\left[\mathit{\sigma }\right]\right]〉}_{\partial {\mathcal{T}}_h}\equiv \sum _{e\in {\mathcal{E}}_h^{{}^{\circ }}}\underset{e}{\int }\mathit{\gamma }\left[\left[\mathit{\sigma }\right]\right]\mathrm{d}s\\ & +\sum _{e\in {\mathcal{E}}_h^{\partial }}\underset{e}{\int }\mathit{\gamma }\mathit{\sigma }\cdot \mathbf{n}\mathrm{d}s,\mathit{\sigma }\in V_{\mathrm{1}},\mathit{\gamma }\in V_{\mathrm{0}},\mathit{\alpha }=\mathrm{2},q=\mathrm{0}.\end{array}\end{array}$$

In general, a finite element form will consist of integrals over various geometric domains: integration over cells 𝒯_h, interior facets ℰ_h^∘, and exterior facets ${\mathcal{E}}_h^{\partial }$. Therefore, we express a general multi-linear form in terms of integrals over each set of geometric entities:

$$\begin{array}{}\text{(10)}& \begin{array}{rl}& a(\mathbf{c};\mathbf{v})=\sum _{K\in {\mathcal{T}}_h}\underset{K}{\int }{\mathcal{I}}_K^{\mathcal{T}}(\mathbf{c};\mathbf{v})\mathrm{d}x+\sum _{e\in {\mathcal{E}}_h^{{}^{\circ }}}\underset{e}{\int }{\mathcal{I}}_e^{\mathcal{E}{,}^{\circ }}(\mathbf{c};\mathbf{v})\mathrm{d}s\\ & +\sum _{e\in {\mathcal{E}}_h^{\partial }}\underset{e}{\int }{\mathcal{I}}_e^{\mathcal{E},\partial }(\mathbf{c};\mathbf{v})\mathrm{d}s,\end{array}\end{array}$$

where ${\mathcal{I}}_K^{\mathcal{T}}$ denotes a cell integrand on K∈𝒯_h, ${\mathcal{I}}_e^{\mathcal{E},{}^{\circ }}$ is an integrand on the interior facet e∈ℰ_h^∘, and ${\mathcal{I}}_e^{\mathcal{E},\partial }$ is an integrand defined on the exterior facet $e\in {\mathcal{E}}_h^{\partial }$. The form a(c;v) describes a finite element form globally over the entire problem domain.

Here, we will consider the case where the interior facet integrands ${\mathcal{I}}_e^{\mathcal{E},{}^{\circ }}$ (c;v) can be decomposed into two independent parts on each interior facet e: one for the positive restriction (+) and the negative restriction (−). That is, for each e∈ℰ_h^∘, we may write ${\mathcal{I}}_e^{\mathcal{E},{}^{\circ }}(\mathbf{c};\mathbf{v})={\mathcal{I}}_{e^{+}}^{\mathcal{E},{}^{\circ }}(\mathbf{c};\mathbf{v})+{\mathcal{I}}_{e^{-}}^{\mathcal{E},{}^{\circ }}(\mathbf{c};\mathbf{v})$. This allows us to express the integral over an interior facet e connecting two adjacent elements, say K⁺ and K⁻, as the sum of integrals:

$$\begin{array}{}\text{(11)}& \begin{array}{rl}& \underset{e\subset \partial K^{+}\cup \partial K^{-}}{\int }{\mathcal{I}}_e^{\mathcal{E}{,}^{\circ }}(\mathbf{c};\mathbf{v})\mathrm{d}s=\underset{e\subset \partial K^{+}}{\int }{\mathcal{I}}_{e^{+}}^{\mathcal{E}{,}^{\circ }}(\mathbf{c};\mathbf{v})\mathrm{d}s\\ & +\underset{e\subset \partial K^{-}}{\int }{\mathcal{I}}_{e^{-}}^{\mathcal{E}{,}^{\circ }}(\mathbf{c};\mathbf{v})\mathrm{d}s.\end{array}\end{array}$$

The local contribution of Eq. (10) in each cell K, along with its associated facets $e\subset \partial K$, is then

$$\begin{array}{}\text{(12)}& \begin{array}{rl}& a_K(\mathbf{c};\mathbf{v})=\underset{K}{\int }{\mathcal{I}}_K^{\mathcal{T}}(\mathbf{c};\mathbf{v})\mathrm{d}x+\sum _{e\subset \partial K\partial \mathrm{\Omega }}\underset{e}{\int }{\mathcal{I}}_e^{\mathcal{E}{,}^{\circ }}(\mathbf{c};\mathbf{v})\mathrm{d}s\\ & +\sum _{e\subset \partial K\cap \partial \mathrm{\Omega }}\underset{e}{\int }{\mathcal{I}}_e^{\mathcal{E},\partial }(\mathbf{c};\mathbf{v})\mathrm{d}s.\end{array}\end{array}$$

We call Eq. (12) the cell-local contribution of a(c;v), with

$$\begin{array}{}\text{(13)}& a(\mathbf{c};\mathbf{v})=\sum _{K\in {\mathcal{T}}_h}{a}_K(\mathbf{c};\mathbf{v}).\end{array}$$

To make matters concrete, let us suppose a(c;v) is a bilinear form with arguments $\mathbf{v}=(v_{\mathrm{0}},v_{\mathrm{1}})\in V_{\mathrm{0}}\times V_{\mathrm{1}}$. Now let $\mathit{\{}{\mathrm{\Phi }}_i{\mathit{\}}}_{i=\mathrm{1}}^N$ and $\mathit{\{}{\mathrm{\Psi }}_i{\mathit{\}}}_{i=\mathrm{1}}^M$ denote bases for V₀ and V₁, respectively. Then the global N×M matrix A corresponding to $a(\mathbf{c};v_{\mathrm{0}},v_{\mathrm{1}})$ has its entries defined via

$$\begin{array}{}\text{(14)}& \begin{array}{rl}& {\mathbf{A}}_{ij}=a\left(\mathbf{c};{\mathrm{\Phi }}_i,{\mathrm{\Psi }}_j\right)=\sum _{K\in {\mathcal{T}}_h}{\mathbf{A}}_{K,ij},{\mathbf{A}}_{K,ij}\\ & =a_K\left(\mathbf{c};{\mathrm{\Phi }}_i,{\mathrm{\Psi }}_j\right).\end{array}\end{array}$$

By construction, ${\mathbf{A}}_{K,ij}\ne \mathbf{0}$ if and only if Φ_i and Ψ_j take non-zero values in K. Now we introduce the cell-node map $i=e(K,\widehat{i})$ as the mapping from the local node number $\widehat{i}$ in K to the global node number i. Suppose there are n and m nodes defining the degrees of freedom for V₀ and V₁, respectively, in K. Then all non-zero entries of A_K,ij arise from integrals involving basis functions with local indices corresponding to the global indices i,j:

$$\begin{array}{}\text{(15)}& \begin{array}{rl}& {\mathbf{A}}_{\widehat{i}\widehat{j}}^K:=a_K\left(\mathbf{c};{\mathrm{\Phi }}_{e(K,\widehat{i})},{\mathrm{\Psi }}_{e(K,\widehat{j})}\right),\widehat{i}\in \left[\mathrm{1},\mathrm{\cdots },n\right],\\ & \widehat{j}\in \left[\mathrm{1},\mathrm{\cdots },m\right].\end{array}\end{array}$$

These local contributions are collected in the n×m dense matrix A^K, which we call the element tensor. The global matrix A is assembled from the collection of element tensors: $\mathbf{A}\leftarrow \mathit{\{}{\mathbf{A}}^K{\mathit{\}}}_{K\in {\mathcal{T}}_h}$. For details on the general evaluation of finite element basis functions and multi-linear forms, we refer the reader to Kirby (2004), Kirby and Logg (2006), Logg et al. (2012b), and Homolya et al. (2018). Further details on the global assembly of finite element operators, with a particular focus on code generation, are summarized in the work of Logg and Wells (2010) and Markall et al. (2013).

In standard finite element software packages, the element tensor is mapped entry-wise into a global sparse array using the cell-node map $e(K,\cdot )$. Within Firedrake, this operation is handled by PyOP2 (Rathgeber et al., 2012) and serves as the main user-facing abstraction for global finite element assembly. For many applications, one may want to produce a new global operator by algebraically manipulating different element tensors. This is relatively invasive in numerical code, as it requires bypassing direct operator assembly to produce the new tensor. This is precisely the scope of Slate.

Like UFL, Slate relies on the grammar of the host-language: Python. The entire Slate language is implemented as a Python module which defines its types (classes) and operations on said types. Together, this forms a high-level language for expressing dense linear algebra on element tensors. The Slate language consists of two primary abstractions for linear algebra:

1. terminal element tensors corresponding to multi-linear integral forms (matrices, vectors, and scalars) or assembled data (for example, coefficient vectors of a finite element function) and

2. expressions consisting of algebraic operations on terminal tensors.

The composition of binary and unary operations on terminal tensors produces a Slate expression. Such expressions can be composed with other Slate objects in arbitrary ways, resulting in concise representations of complex algebraic operations on locally assembled arrays. We summarize all currently supported Slate abstractions here.

2.1.1 Terminal tensors

In Slate, one associates a tensor with data on a cell either by using a multi-linear form, or assembled coefficient data:

• Tensor(a(c;v))
  associates a form, expressed in UFL, with its local element tensor

  $$\begin{array}{}\text{(16)}& {\mathbf{A}}^K\leftarrow a_K(\mathbf{c};\mathbf{v}),\text{for all}K\in {\mathcal{T}}_h.\end{array}$$

  The form arity α of a_K(c;v) determines the rank of the corresponding Tensor; i.e., scalars, vectors, and matrices are produced from scalars, linear forms, and bilinear forms, respectively.¹ The shape of the element tensor is determined by both the number of arguments and total number of degrees of freedom local to the cell.

• AssembledVector(f),
  where f is some finite element function. The function f∈V is expressed in terms of the finite element basis of V: $f\left(x\right)={\sum }_{i=\mathrm{1}}^N{f}_i{\mathrm{\Phi }}_i\left(x\right)$. The result is the local coefficient vector of f on K:

  $$\begin{array}{}\text{(17)}& {\mathbf{F}}^K\leftarrow {\left\{f_{e(K,\widehat{i})}\right\}}_{\widehat{i}=\mathrm{1}}^n,\end{array}$$

  where $e(K,\widehat{i})$ is the local node numbering and n is the number of nodes local to the cell K.

2.1.2 Symbolic linear algebra

Slate supports typical binary and unary operations in linear algebra, with a high-level syntax close to mathematics. At the time of this paper, these include the following.

• A + B, the addition of two equal shaped tensors: A^K+B^K.

• A ∗ B, a contraction over the last index of A and the first index of B. This is the usual multiplicative operation on matrices, vectors, and scalars: A^KB^K.

• -A, the additive inverse (negation) of a tensor: −A^K.

• A.T, the transpose of a tensor: ${\left({\mathbf{A}}^K\right)}^T$.

• A.inv, the inverse of a square tensor: ${\left({\mathbf{A}}^K\right)}^{-\mathrm{1}}$.

• A.solve(B, decomposition=``...''), the result, X^K, of solving a local linear system A^KX^K=B^K, optionally specifying a factorization strategy.

• A.blocks[indices], where A is a tensor from a mixed finite element space. This allows for the extraction of sub-blocks, which are indexed by field (slices are allowed). For example, if a matrix A corresponds to the bilinear form $a:V\times W\to \mathbb{R}$, where $V=V_{\mathrm{0}}\times \mathrm{\cdots }\times V_n$ and $W=W_{\mathrm{0}}\times \mathrm{\cdots }\times W_m$ are product spaces consisting of finite element spaces $\mathit{\{}{V}_i{\mathit{\}}}_{i=\mathrm{0}}^n$, $\mathit{\{}{W}_i{\mathit{\}}}_{i=\mathrm{0}}^m$, then the element tensors have the form

  $$\begin{array}{}\text{(18)}& {\mathbf{A}}^K=\left[\begin{array}{llll}{\mathbf{A}}_{\mathrm{00}}^K& {\mathbf{A}}_{\mathrm{01}}^K& \mathrm{\cdots }& {\mathbf{A}}_{\mathrm{0}m}^K\\ {\mathbf{A}}_{\mathrm{10}}^K& {\mathbf{A}}_{\mathrm{11}}^K& \mathrm{\cdots }& {\mathbf{A}}_{\mathrm{1}m}^K\\ \mathrm{⋮}& \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\\ {\mathbf{A}}_{n\mathrm{0}}^K& {\mathbf{A}}_{n\mathrm{1}}^K& \mathrm{\cdots }& {\mathbf{A}}_{nm}^K\end{array}\right].\end{array}$$

  The associated submatrix of Eq. (18) with indices $\mathbf{i}=(\mathbf{p},\mathbf{q})$, $\mathbf{p}=\mathit{\{}{p}_{\mathrm{1}},\mathrm{\cdots },p_r\mathit{\}}$, $\mathbf{q}=\mathit{\{}{q}_{\mathrm{1}},\mathrm{\cdots },q_c\mathit{\}}$, is

  $$\begin{array}{}\text{(19)}& {\mathbf{A}}_{\mathbf{pq}}^K=\left[\begin{array}{lll}{\mathbf{A}}_{p_{\mathrm{1}}{q}_{\mathrm{1}}}^K& \mathrm{\cdots }& {\mathbf{A}}_{p_{\mathrm{1}}{q}_c}^K\\ \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\\ {\mathbf{A}}_{p_r{q}_{\mathrm{1}}}^K& \mathrm{\cdots }& {\mathbf{A}}_{p_r{q}_c}^K\end{array}\right]={\mathbf{A}}^K\mathtt{\text{.blocks}}[\mathbf{p},\mathbf{q}],\end{array}$$

  where $\mathbf{p}\subseteq \mathit{\{}\mathrm{0},\mathrm{\cdots },n\mathit{\}}$, $\mathbf{q}\subseteq \mathit{\{}\mathrm{0},\mathrm{\cdots },m\mathit{\}}$.

Each Tensor object knows all the information about the underlying UFL form that defines it, such as form arguments, coefficients, and the underlying finite element space(s) it operates on. This information is propagated through as unary or binary transformations are applied. The unary and binary operations shown here provide the necessary algebraic framework for a large class of problems, some of which we present in this paper.

In Firedrake, Slate expressions are transformed into low-level code by a linear algebra compiler. The compiler interprets Slate expressions as a syntax tree, where the tree is visited to identify what local arrays need to be assembled and the sequence of array operations. At the time of this work, our compiler generates C++ code, using the templated library Eigen (Guennebaud et al., 2015) for dense linear algebra. The translation from Slate to C++ is fairly straightforward, as all operations supported by Slate have a representation in Eigen.

The compiler pass will generate a single “macro” kernel, which performs the dense linear algebra operations represented in Slate. The resulting code will also include (often multiple) function calls to local assembly kernels generated by TSFC (Two Stage Form Compiler) (Homolya et al., 2018) to assemble all necessary sub-blocks of an element tensor. All code generated by the linear algebra compiler conforms to the application programming interface (API) of the PyOP2 framework, as detailed by Rathgeber et al. (2012, Sect. 3). Figure 1 provides an illustration of the complete tool chain.

Figure 1 The Slate language wraps UFL objects describing the finite element system. The resulting Slate expressions are passed to a specialized linear algebra compiler, which produces a single macro kernel assembling the local contributions and executes the dense linear algebra represented in Slate. The kernels are passed to the Firedrake's PyOP2 interface, which wraps the Slate kernel in a mesh-iteration kernel. Parallel scheduling, code generation, and compilation occurs after the PyOP2 layer.

Most optimization of the resulting dense linear algebra code is handled directly by Eigen. In the case of unary and binary operations such as A.inv and A.solve(B), stable default behaviors are applied by the linear algebra compiler. For example, A.solve(B) without a specified factorization strategy will default to using an in-place LU factorization with partial pivoting. For local matrices smaller than 5×5, the inverse is translated directly into Eigen's A.inverse(), which employs stable analytic formulas. For larger matrices, the linear algebra replaces A.inv with an LU factorization.² Currently, we only support direct matrix factorizations for solving local linear systems. However, it would not be difficult to extend Slate to support more general solution techniques like iterative methods.

3 Examples

We now present examples and discuss solution methods which require element-wise manipulations of finite element systems and their specification in Slate. We stress here that Slate is not limited to these model problems; rather these examples are chosen for clarity and to demonstrate key features of the Slate language. For our discussion, we use a model elliptic equation defined in a computational domain Ω. Consider the second-order PDE with both Dirichlet and Neumann boundary conditions:

$$\begin{array}{}\text{(20)}& {\displaystyle }-\mathrm{\nabla }\cdot \left(\mathit{\kappa }\mathrm{\nabla }p\right)+cp=f,\text{in}\mathrm{\Omega },\text{(21)}& {\displaystyle }p=p_{\mathrm{0}},\text{on}\partial {\mathrm{\Omega }}_D,\text{(22)}& {\displaystyle }-\mathit{\kappa }\mathrm{\nabla }p\cdot \mathbf{n}=g,\text{on}\partial {\mathrm{\Omega }}_N,\end{array}$$

where $\partial {\mathrm{\Omega }}_D\cup \partial {\mathrm{\Omega }}_N=\partial \mathrm{\Omega }$ and κ, $c:\mathrm{\Omega }\to {\mathbb{R}}^{+}$ are positive-valued coefficients. To obtain a mixed formulation of Eqs. (20)–(22), we introduce the auxiliary velocity variable $\mathbf{u}=-\mathit{\kappa }\mathrm{\nabla }p$. We then obtain the first-order system of PDEs:

$$\begin{array}{}\text{(23)}& {\displaystyle }\mathit{\mu }\mathbf{u}+\mathrm{\nabla }p=\mathrm{0},\text{in}\mathrm{\Omega },\text{(24)}& {\displaystyle }\mathrm{\nabla }\cdot \mathbf{u}+cp=f,\text{in}\mathrm{\Omega },\text{(25)}& {\displaystyle }p=p_{\mathrm{0}},\text{on}\partial {\mathrm{\Omega }}_D,\text{(26)}& {\displaystyle }\mathbf{u}\cdot \mathbf{n}=g,\text{on}\partial {\mathrm{\Omega }}_N,\end{array}$$

where $\mathit{\mu }={\mathit{\kappa }}^{-\mathrm{1}}$.

3.1 Hybridization of mixed methods

To motivate our discussion in this section, we start by recalling the mixed method for Eqs. (23)–(26). Methods of this type seek approximations (u_h,p_h) in finite-dimensional subspaces $U_h\times V_h\subset H(\mathrm{div};\mathrm{\Omega })\times L^{\mathrm{2}}\left(\mathrm{\Omega }\right)$, defined by

$$\begin{array}{}\text{(27)}& \begin{array}{rl}& U_h=\left\{\mathbf{w}\in H(\mathrm{div};\mathrm{\Omega }):\mathbf{w}{|}_K\in U\left(K\right),\forall K\in {\mathcal{T}}_h,\right.\\ & \left.\mathbf{w}\cdot \mathbf{n}=g\text{ on }\partial {\mathrm{\Omega }}_N\right\},\end{array}\end{array}$$

$$\begin{array}{}\text{(28)}& V_h=\left\{\mathit{\varphi }\in L^{\mathrm{2}}\left(\mathrm{\Omega }\right):\mathit{\varphi }{|}_K\in V\left(K\right),\forall K\in {\mathcal{T}}_h\right\}.\end{array}$$

The space U_h consists of H(div)-conforming piecewise vector polynomials, where choices of U(K) typically include the Raviart–Thomas (RT), Brezzi–Douglas–Marini (BDM), or Brezzi–Douglas–Fortin–Marini (BDFM) elements (Raviart and Thomas, 1977; Nédélec, 1980; Brezzi et al., 1985, 1987). The space V_h is the Lagrange family of discontinuous polynomials. These spaces are of particular interest when simulating geophysical flows, since choosing the right pairing results in stable discretizations with desirable conservation properties and avoids spurious computational modes. We refer the reader to Cotter and Shipton (2012), Cotter and Thuburn (2014), Natale et al. (2016), and Shipton et al. (2018) for a discussion of mixed methods relevant for geophysical fluid dynamics. Two examples of such discretizations are presented in Sect. 5.2.

The mixed formulation of Eqs. (23)–(26) is arrived at by multiplying Eqs. (23)–(24) by test functions and integrating by parts. The resulting finite element problem reads as follows: find $({\mathbf{u}}_h,p_h)\in U_h\times V_h$ satisfying

$$\begin{array}{}\text{(29)}& \begin{array}{rl}& {\left(\mathbf{w},\mathit{\mu }{\mathbf{u}}_h\right)}_{{\mathcal{T}}_h}-{\left(\mathrm{\nabla }\cdot \mathbf{w},p_h\right)}_{{\mathcal{T}}_h}=-{〈\mathbf{w}\cdot \mathbf{n},p_{\mathrm{0}}〉}_{\partial {\mathrm{\Omega }}_D},\\ & \forall \mathbf{w}\in U_{h,\mathrm{0}},\end{array}\end{array}$$

$$\begin{array}{}\text{(30)}& {\left(\mathit{\varphi },\mathrm{\nabla }\cdot {\mathbf{u}}_h\right)}_{{\mathcal{T}}_h}+{\left(\mathit{\varphi },c{p}_h\right)}_{{\mathcal{T}}_h}={\left(\mathit{\varphi },f\right)}_{{\mathcal{T}}_h},\forall \mathit{\varphi }\in V_h,\end{array}$$

where U_h,0 is the subspace of U_h with functions whose normal components vanish on ∂Ω_N. The discrete system is obtained by first expanding the solutions in terms of the finite element bases:

$$\begin{array}{}\text{(31)}& {\mathbf{u}}_h=\sum _{i=\mathrm{1}}^{N_{\mathbf{u}}}{U}_i{\mathbf{\Psi }}_i,p_h=\sum _{i=\mathrm{1}}^{N_p}{P}_i{\mathit{\xi }}_i,\end{array}$$

where $\mathit{\{}{\mathbf{\Psi }}_i{\mathit{\}}}_{i=\mathrm{1}}^{N_{\mathbf{u}}}$ and $\mathit{\{}{\mathit{\xi }}_i{\mathit{\}}}_{i=\mathrm{1}}^{N_p}$ are bases for U_h and V_h, respectively. Here, U_i and P_i are the coefficients to be determined. As per standard Galerkin-based finite element methods, taking w=Ψ_j, $j\in \mathit{\{}\mathrm{1},\mathrm{\cdots },N_{\mathbf{u}}\mathit{\}}$ and ϕ=ξ_j, $j\in \mathit{\{}\mathrm{1},\mathrm{\cdots },N_p\mathit{\}}$ in Eqs. (29)–(30) produces the following discrete saddle point system:

$$\begin{array}{}\text{(32)}& \left[\begin{array}{ll}\mathbf{A}& -{\mathbf{B}}^T\\ \mathbf{B}& \mathbf{D}\end{array}\right]\left\{\begin{array}{l}\mathbf{U}\\ \mathbf{P}\end{array}\right\}=\left\{\begin{array}{l}{\mathbf{F}}_{\mathrm{0}}\\ {\mathbf{F}}_{\mathrm{1}}\end{array}\right\}.\end{array}$$

where $\mathbf{U}=\mathit{\{}{U}_i{\mathit{\}}}_{i=\mathrm{1}}^{N_{\mathbf{u}}}$, $\mathbf{P}=\mathit{\{}{P}_i{\mathit{\}}}_{i=\mathrm{1}}^{N_p}$ are the coefficient vectors, and

$$\begin{array}{}\text{(33)}& {\displaystyle }{\mathbf{A}}_{ij}& {\displaystyle }={\left({\mathbf{\Psi }}_i,\mathit{\mu }{\mathbf{\Psi }}_j\right)}_{{\mathcal{T}}_h},\text{(34)}& {\displaystyle }{\mathbf{B}}_{ij}& {\displaystyle }={\left({\mathit{\xi }}_i,\mathrm{\nabla }\cdot {\mathbf{\Psi }}_j\right)}_{{\mathcal{T}}_h},\text{(35)}& {\displaystyle }{\mathbf{D}}_{ij}& {\displaystyle }={\left({\mathit{\xi }}_i,c{\mathit{\xi }}_j\right)}_{{\mathcal{T}}_h},\text{(36)}& {\displaystyle }{\mathbf{F}}_{\mathrm{0},j}& {\displaystyle }=-{〈{\mathbf{\Psi }}_j\cdot \mathbf{n},p_{\mathrm{0}}〉}_{\partial {\mathrm{\Omega }}_D},\text{(37)}& {\displaystyle }{\mathbf{F}}_{\mathrm{1},j}& {\displaystyle }={\left({\mathit{\xi }}_j,f\right)}_{{\mathcal{T}}_h}.\end{array}$$

Methods to efficiently invert such systems include H(div)-multigrid (Arnold et al., 2000) (requiring complex overlapping Schwarz smoothers), global Schur complement factorizations (which require an approximation to the inverse of the dense³ elliptic Schur complement $\mathbf{D}+{\mathbf{BA}}^{-\mathrm{1}}{\mathbf{B}}^T$), or auxiliary space multigrid (Hiptmair and Xu, 2007). Here, we focus on a solution approach using a hybridized mixed method (Arnold and Brezzi, 1985; Brezzi and Fortin, 1991; Boffi et al., 2013).

The hybridization technique replaces the original system with a discontinuous variant, decoupling the velocity degrees of freedom between cells. This is done by replacing the discrete solution space for u_h with the “broken” space $U_h^d$, defined as

$$\begin{array}{}\text{(38)}& U_h^d=\left\{\mathbf{w}\in {\left[L^{\mathrm{2}}\left(\mathrm{\Omega }\right)\right]}^n:\mathbf{w}{|}_K\in U\left(K\right),\forall K\in {\mathcal{T}}_h\right\}.\end{array}$$

The vector finite element space $U_h^d$ is a subspace of ${\left[L^{\mathrm{2}}\left(\mathrm{\Omega }\right)\right]}^n$ consisting of local H(div) functions, but normal components are no longer required to be continuous on ∂𝒯_h. The approximation space for p_h remains unchanged.

Next, Lagrange multipliers are introduced as an auxiliary variable in the space M_h, defined only on cell interfaces:

$$\begin{array}{}\text{(39)}& M_h=\mathit{\{}\mathit{\gamma }\in L^{\mathrm{2}}(\partial {\mathcal{T}}_h):\mathit{\gamma }{|}_e\in M(e),\forall e\in \partial {\mathcal{T}}_h\mathit{\}},\end{array}$$

where M(e) denotes a polynomial space defined on each facet. We call M_h the space of approximate traces. Functions in M_h are discontinuous across vertices in two dimensions and vertices or edges in three dimensions.

Deriving the hybridizable mixed system is accomplished through integration by parts over each element K. Testing with $\mathbf{w}\in U_h^d\left(K\right)$ and integrating Eq. (23) over the cell K produces

$$\begin{array}{}\text{(40)}& \begin{array}{rl}& {\left(\mathbf{w},\mathit{\mu }{\mathbf{u}}_h^d\right)}_K-{\left(\mathrm{\nabla }\cdot \mathbf{w},p_h\right)}_K+{〈\mathbf{w}\cdot \mathbf{n},{\mathit{\lambda }}_h〉}_{\partial K}\\ & =-{〈\mathbf{w}\cdot \mathbf{n},p_{\mathrm{0}}〉}_{\partial K\cap \partial {\mathrm{\Omega }}_D}.\end{array}\end{array}$$

The trace function λ_h is introduced in the surface integral as an approximation to $p{|}_{\partial K}$. An additional constraint equation, called the transmission condition, is added to close the system. The resulting hybridizable formulation reads as follows: find $({\mathbf{u}}_h^d,p_h,{\mathit{\lambda }}_h)\in U_h^d\times V_h\times M_h$ such that

$$\begin{array}{}\text{(41)}& \begin{array}{rl}& {\left(\mathbf{w},\mathit{\mu }{\mathbf{u}}_h^d\right)}_{{\mathcal{T}}_h}-{\left(\mathrm{\nabla }\cdot \mathbf{w},p_h\right)}_{{\mathcal{T}}_h}+{〈\left[\left[\mathbf{w}\right]\right],{\mathit{\lambda }}_h〉}_{\partial {\mathcal{T}}_h\setminus \partial {\mathrm{\Omega }}_D}\\ & =-{〈\mathbf{w}\cdot \mathbf{n},p_{\mathrm{0}}〉}_{\partial {\mathrm{\Omega }}_D},\forall \mathbf{w}\in U_h^d,\end{array}\end{array}$$

$$\begin{array}{}\text{(42)}& {\left(\mathit{\varphi },\mathrm{\nabla }\cdot {\mathbf{u}}_h^d\right)}_{{\mathcal{T}}_h}+{\left(\mathit{\varphi },c{p}_h\right)}_{{\mathcal{T}}_h}={\left(\mathit{\varphi },f\right)}_{{\mathcal{T}}_h},\forall \mathit{\varphi }\in V_h,\end{array}$$

$$\begin{array}{}\text{(43)}& {〈\mathit{\gamma },\left[\left[{\mathbf{u}}_h^d\right]\right]〉}_{\partial {\mathcal{T}}_h\setminus \partial {\mathrm{\Omega }}_D}={〈\mathit{\gamma },g〉}_{\partial {\mathrm{\Omega }}_N},\forall \mathit{\gamma }\in M_{h,\mathrm{0}},\end{array}$$

where M_h,0 denotes the space of traces vanishing on ∂Ω_D. The transmission condition Eq. (43) enforces both the continuity of ${\mathbf{u}}_h^d\cdot \mathbf{n}$ across element boundaries as well as the boundary condition ${\mathbf{u}}_h^d\cdot \mathbf{n}=g$ on ∂Ω_N. If the space of Lagrange multipliers M_h is chosen appropriately, then the broken velocity ${\mathbf{u}}_h^d$, albeit sought a priori in a discontinuous space, will coincide with its H(div)-conforming counterpart. Specifically, the formulations in Eqs. (41)–(42) and (29)–(30) are solving equivalent problems if the normal components of w∈U_h lie in the same polynomial space as the trace functions (Arnold and Brezzi, 1985).

The discrete matrix system arising from Eqs. (41)–(43) has the general form

$$\begin{array}{}\text{(44)}& \left[\begin{array}{lll}{\mathbf{A}}_{\mathrm{00}}& {\mathbf{A}}_{\mathrm{01}}& {\mathbf{A}}_{\mathrm{02}}\\ {\mathbf{A}}_{\mathrm{10}}& {\mathbf{A}}_{\mathrm{11}}& {\mathbf{A}}_{\mathrm{12}}\\ {\mathbf{A}}_{\mathrm{20}}& {\mathbf{A}}_{\mathrm{21}}& {\mathbf{A}}_{\mathrm{22}}\end{array}\right]\left\{\begin{array}{l}{\mathbf{U}}^d\\ \mathbf{P}\\ \mathbf{\Lambda }\end{array}\right\}=\left\{\begin{array}{l}{\mathbf{F}}_{\mathrm{0}}\\ {\mathbf{F}}_{\mathrm{1}}\\ {\mathbf{F}}_{\mathrm{2}}\end{array}\right\},\end{array}$$

where the discrete system is produced by expanding functions in terms of the finite element bases for $U_h^d$, V_h, and M_h like before. Upon initial inspection, it may not appear to be advantageous to replace our original formulation with this augmented equation set; the hybridizable system has substantially more total degrees of freedom. However, Eq. (44) has a considerable advantage over Eq. (32) in the following ways.

1. Since both $U_h^d$ and V_h are discontinuous spaces, U^d and P are coupled only within the cell. This allows us to simultaneously eliminate both unknowns via local static condensation to produce a significantly smaller global (hybridized) problem for the trace unknowns, Λ:

   $$\begin{array}{}\text{(45)}& \mathbf{S}\mathbf{\Lambda }=\mathbf{E},\end{array}$$

   where $\mathbf{S}\leftarrow \mathit{\{}{\mathbf{S}}^K{\mathit{\}}}_{K\in {\mathcal{T}}_h}$ and $\mathbf{E}\leftarrow \mathit{\{}{\mathbf{E}}^K{\mathit{\}}}_{K\in {\mathcal{T}}_h}$ are assembled via the local element tensors:

   $$\begin{array}{}\text{(46)}& {\mathbf{S}}^K={\mathbf{A}}_{\mathrm{22}}^K-\left[\begin{array}{ll}{\mathbf{A}}_{\mathrm{20}}^K& {\mathbf{A}}_{\mathrm{21}}^K\end{array}\right]{\left[\begin{array}{ll}{\mathbf{A}}_{\mathrm{00}}^K& {\mathbf{A}}_{\mathrm{01}}^K\\ {\mathbf{A}}_{\mathrm{10}}^K& {\mathbf{A}}_{\mathrm{11}}^K\end{array}\right]}^{-\mathrm{1}}\left[\begin{array}{l}{\mathbf{A}}_{\mathrm{02}}^K\\ {\mathbf{A}}_{\mathrm{12}}^K\end{array}\right],\end{array}$$
   $$\begin{array}{}\text{(47)}& {\mathbf{E}}^K={\mathbf{F}}_{\mathrm{2}}^K-\left[\begin{array}{ll}{\mathbf{A}}_{\mathrm{20}}^K& {\mathbf{A}}_{\mathrm{21}}^K\end{array}\right]{\left[\begin{array}{ll}{\mathbf{A}}_{\mathrm{00}}^K& {\mathbf{A}}_{\mathrm{01}}^K\\ {\mathbf{A}}_{\mathrm{10}}^K& {\mathbf{A}}_{\mathrm{11}}^K\end{array}\right]}^{-\mathrm{1}}\left\{\begin{array}{l}{\mathbf{F}}_{\mathrm{0}}^K\\ {\mathbf{F}}_{\mathrm{1}}^K\end{array}\right\}.\end{array}$$

   Note that the inverse of the block matrix in Eqs. (46) and (47) is never evaluated globally; the elimination can be performed locally by performing a sequence of Schur complement reductions within each cell.

2. The matrix S is sparse, symmetric, positive-definite, and spectrally equivalent to the dense Schur complement $\mathbf{D}+{\mathbf{BA}}^{-\mathrm{1}}{\mathbf{B}}^T$ from Eq. (32) of the original mixed formulation (Gopalakrishnan, 2003; Cockburn et al., 2009a).

3. Once Λ is computed, both U^d and P can be recovered locally in each element. This can be accomplished in a number of ways. One way is to compute P^K by solving

   $$\begin{array}{}\text{(48)}& \begin{array}{rl}& \left({\mathbf{A}}_{\mathrm{11}}^K-{\mathbf{A}}_{\mathrm{10}}^K{\left({\mathbf{A}}_{\mathrm{00}}^K\right)}^{-\mathrm{1}}{\mathbf{A}}_{\mathrm{01}}^K\right){\mathbf{P}}^K={\mathbf{F}}_{\mathrm{1}}^K-{\mathbf{A}}_{\mathrm{10}}^K\\ & {\left({\mathbf{A}}_{\mathrm{00}}^K\right)}^{-\mathrm{1}}{\mathbf{F}}_{\mathrm{0}}^K-\left({\mathbf{A}}_{\mathrm{12}}^K-{\mathbf{A}}_{\mathrm{10}}^K{\left({\mathbf{A}}_{\mathrm{00}}^K\right)}^{-\mathrm{1}}{\mathbf{A}}_{\mathrm{02}}^K\right){\mathbf{\Lambda }}^K,\end{array}\end{array}$$

   followed by solving for ${\left({\mathbf{U}}^d\right)}^K$:

   $$\begin{array}{}\text{(49)}& {\mathbf{A}}_{\mathrm{00}}^K{\left({\mathbf{U}}^d\right)}^K={\mathbf{F}}_{\mathrm{0}}^K-{\mathbf{A}}_{\mathrm{01}}^K{\mathbf{P}}^K-{\mathbf{A}}_{\mathrm{02}}^K{\mathbf{\Lambda }}^K.\end{array}$$

   Similarly, one could rearrange the order in which each variable is reconstructed.

4. If desired, the solutions can be improved further through local post-processing. We highlight two such procedures for U^d and P, respectively, in Sect. 3.3.

Figure 2 displays the corresponding Slate code for assembling the trace system, solving Eq. (45), and recovering the eliminated unknowns. For a complete reference on how to formulate the hybridized mixed system Eqs. (41)–(43) in UFL, we refer the reader to Alnæs et al. (2014). Complete Firedrake code using Slate to solve a hybridizable mixed system is also publicly available in Zenodo/Tabula-Rasa (2019, “Code verification”). We remark that, in the case of this hybridizable system, Eq. (44) contains zero-valued blocks which can simplify the resulting expressions in Eqs. (46)–(47) and (48)–(49). This is not true in general and therefore the expanded form using all sub-blocks of Eq. (44) is presented for completeness.

Figure 2 Firedrake code for solving Eq. (44) via static condensation and local recovery, given UFL expressions a, L for Eqs. (41)–(43). Arguments of the mixed space $U_h^d\times V_h\times M_h$ are indexed by 0, 1, and 2, respectively. Lines 8 and 9 are symbolic expressions for Eqs. (46) and (47), respectively. Any vanishing conditions on the trace variables should be provided as boundary conditions during operator assembly (line 10). Lines 26 and 28 are expressions for Eqs. (48) and (49) (using LU). Code generation occurs in lines 10, 11, 30, and 31. A global linear solver for the reduced system is created and used in line 15. Configuring the linear solver is done by providing an appropriate Python dictionary of solver options for the PETSc library.

3.2 Hybridization of discontinuous Galerkin methods

The HDG method is a natural extension of discontinuous Galerkin (DG) discretizations. Here, we consider a specific HDG discretization, namely the LDG-H method (Cockburn et al., 2010b). Other forms of HDG that involve local lifting operators can also be implemented in this software framework by the introduction of additional local (i.e., discontinuous) variables into the definition of the local solver.

Deriving the LDG-H formulation follows exactly from standard DG methods. All prognostic variables are sought in the discontinuous spaces $U_h\times V_h\subset {\left[L^{\mathrm{2}}\left(\mathrm{\Omega }\right)\right]}^n\times L^{\mathrm{2}}\left(\mathrm{\Omega }\right)$. Within a cell K, integration by parts yields

$$\begin{array}{}\text{(50)}& \begin{array}{rl}& {\left(\mathbf{w},\mathit{\mu }{\mathbf{u}}_h\right)}_K-{\left(\mathrm{\nabla }\cdot \mathbf{w},p_h\right)}_K+{〈\mathbf{w}\cdot \mathbf{n},\widehat{p}〉}_{\partial K}=\mathrm{0},\\ & \forall \mathbf{w}\in U\left(K\right),\end{array}\end{array}$$

$$\begin{array}{}\text{(51)}& \begin{array}{rl}& -{\left(\mathrm{\nabla }\mathit{\varphi },{\mathbf{u}}_h\right)}_K+{〈\mathit{\varphi },\widehat{\mathbf{u}}\cdot \mathbf{n}〉}_{\partial K}+{\left(\mathit{\varphi },c{p}_h\right)}_K={\left(\mathit{\varphi },f\right)}_K,\\ & \forall \mathit{\varphi }\in V\left(K\right),\end{array}\end{array}$$

where U(K) and V(K) are vector and scalar polynomial spaces, respectively. Now, we define the numerical fluxes $\widehat{p}$ and $\widehat{\mathbf{u}}$ to be functions of the trial unknowns and a new independent unknown in the trace space M_h:

$$\begin{array}{}\text{(52)}& {\displaystyle }& {\displaystyle }\widehat{\mathbf{u}}({\mathbf{u}}_h,p_h,{\mathit{\lambda }}_h;\mathit{\tau })={\mathbf{u}}_h+\mathit{\tau }\left(p_h-\widehat{p}\right)\mathbf{n},\text{(53)}& {\displaystyle }& {\displaystyle }\widehat{p}\left({\mathit{\lambda }}_h\right)={\mathit{\lambda }}_h,\end{array}$$

where λ_h∈M_h is a function approximating p on ∂𝒯_h and τ is a positive stabilization function that may vary on each facet $e\in \partial {\mathcal{T}}_h$. We further require that λ_h satisfies the Dirichlet condition for p on ∂Ω_D in an L²-projection sense. The full LDG-H formulation reads as follows. Find $({\mathbf{u}}_h,p_h,{\mathit{\lambda }}_h)\in U_h\times V_h\times M_h$ such that

$$\begin{array}{}\text{(54)}& \begin{array}{rl}& {\left(\mathbf{w},\mathit{\mu }{\mathbf{u}}_h\right)}_{{\mathcal{T}}_h}-{\left(\mathrm{\nabla }\cdot \mathbf{w},p_h\right)}_{{\mathcal{T}}_h}+{〈\left[\left[\mathbf{w}\right]\right],{\mathit{\lambda }}_h〉}_{\partial {\mathcal{T}}_h}=\mathrm{0},\\ & \forall \mathbf{w}\in U_h,\end{array}\end{array}$$

$$\begin{array}{}\text{(55)}& \begin{array}{rl}& -{\left(\mathrm{\nabla }\mathit{\varphi },{\mathbf{u}}_h\right)}_{{\mathcal{T}}_h}+{〈\mathit{\varphi },\left[\left[{\mathbf{u}}_h+\mathit{\tau }\left(p_h-{\mathit{\lambda }}_h\right)\mathbf{n}\right]\right]〉}_{\partial {\mathcal{T}}_h}\\ & +{\left(\mathit{\varphi },c{p}_h\right)}_{{\mathcal{T}}_h}={\left(\mathit{\varphi },f\right)}_{{\mathcal{T}}_h},\forall \mathit{\varphi }\in V_h,\end{array}\end{array}$$

$$\begin{array}{}\text{(56)}& \begin{array}{rl}& {〈\mathit{\gamma },\left[\left[{\mathbf{u}}_h+\mathit{\tau }\left(p_h-{\mathit{\lambda }}_h\right)\mathbf{n}\right]\right]〉}_{\partial {\mathcal{T}}_h\setminus \partial {\mathrm{\Omega }}_D}={〈\mathit{\gamma },g〉}_{\partial {\mathrm{\Omega }}_N},\\ & \forall \mathit{\gamma }\in M_h,\end{array}\end{array}$$

$$\begin{array}{}\text{(57)}& {〈\mathit{\gamma },{\mathit{\lambda }}_h〉}_{\partial {\mathrm{\Omega }}_D}={〈\mathit{\gamma },p_{\mathrm{0}}〉}_{\partial {\mathrm{\Omega }}_D},\forall \mathit{\gamma }\in M_h.\end{array}$$

Equation (56) is the transmission condition, which enforces the continuity of $\widehat{\mathbf{u}}\cdot \mathbf{n}$ on ∂𝒯_h and q. Equation (57) ensures λ_h satisfies the Dirichlet condition. This ensures that the numerical flux is single-valued on the facets. Hence, the LDG-H method defines a conservative DG method (Cockburn et al., 2010b). Note that the choice of τ has a significant influence on the expected convergence rates of the computed solutions.

The LDG-H method retains the advantages of standard DG methods while also enabling the assembly of reduced linear systems through static condensation. The matrix system arising from Eqs. (54)–(57) has the same general form as the hybridized mixed method in Eq. (44), except all sub-blocks are now populated with non-zero entries due to the coupling of trace functions with both p_h and u_h. However, all previous properties of the discrete matrix system from Sect. 3.1 still apply. The Slate expressions for the local elimination and reconstruction operations will be identical to those illustrated in Fig. 2. For the interested reader, a unified analysis of hybridization methods (both mixed and DG) for second-order elliptic equations is presented in Cockburn et al. (2009a) and Cockburn (2016).

3.3 Local post-processing

For both mixed (Arnold and Brezzi, 1985; Brezzi et al., 1985; Bramble and Xu, 1989; Stenberg, 1991) and discontinuous Galerkin methods (Cockburn et al., 2010b, 2009b), it is possible to locally post-process solutions to obtain superconvergent approximations (gaining 1 order of accuracy over the unprocessed solution). These methods can be expressed as local solves on each element and are straightforward to implement using Slate. In this section, we present two post-processing techniques: one for scalar fields and another for the vector unknown. The Slate code follows naturally from previous discussions in Sect. 3.1 and 3.2, using the standard set of operations on element tensors summarized in Sect. 2.1.

3.3.1 Post-processing of the scalar solution

Our first example is a modified version of the procedure presented by Stenberg (1991) for enhancing the accuracy of the scalar solution. This was also highlighted within the context of hybridizing eigenproblems by Cockburn et al. (2010a). This post-processing technique can be used for both the hybridizable mixed and LDG-H methods. We proceed by posing the finite element systems cell-wise.

Let 𝒫_k(K) denote a polynomial space of degree ≤k on a cell K∈𝒯_h. Then for a given pair of computed solutions u_h,p_h of the hybridized methods, we define the post-processed scalar $p_h^{\star }\in {\mathcal{P}}_{k+\mathrm{1}}\left(K\right)$ as the unique solution of the local problem:

$$\begin{array}{}\text{(58)}& {\displaystyle }& {\displaystyle }{\left(\mathrm{\nabla }w,\mathrm{\nabla }{p}_h^{\star }\right)}_K=-{\left(\mathrm{\nabla }w,{\mathit{\kappa }}^{-\mathrm{1}}{\mathbf{u}}_h\right)}_K,\forall w\in {\mathcal{P}}_{k+\mathrm{1}}^{⟂,l}\left(K\right),\text{(59)}& {\displaystyle }& {\displaystyle }{\left(v,p_h^{\star }\right)}_K={\left(v,p_h\right)}_K,\forall v\in {\mathcal{P}}_l\left(K\right),\end{array}$$

where $\mathrm{0}\le l\le k$. Here, the space ${\mathcal{P}}_{k+\mathrm{1}}^{⟂,l}\left(K\right)$ denotes the L²-orthogonal complement of 𝒫_l(K). This post-processing method directly uses the definition of the flux u_h, the approximation of −κ∇p. In practice, the space ${\mathcal{P}}_{k+\mathrm{1}}^{⟂,l}\left(K\right)$ may be constructed using an orthogonal hierarchical basis, and solving Eqs. (58)–(59) amounts to inverting a symmetric positive definite system in each cell of the mesh.

At the time of this work, Firedrake does not support the construction of such a finite element basis. However, we can introduce Lagrange multipliers to enforce the orthogonality constraint. The resulting local problem then becomes the following mixed system: find $(p_h^{\star },\mathit{\psi })\in {\mathcal{P}}_{k+\mathrm{1}}\left(K\right)\times {\mathcal{P}}_l\left(K\right)$ such that

$$\begin{array}{}\text{(60)}& \begin{array}{rlr}& {\left(\mathrm{\nabla }w,\mathrm{\nabla }{p}_h^{\star }\right)}_K+{\left(w,\mathit{\psi }\right)}_K& =-{\left(\mathrm{\nabla }w,{\mathit{\kappa }}^{-\mathrm{1}}{\mathbf{u}}_h\right)}_K,\\ & \forall w\in {\mathcal{P}}_{k+\mathrm{1}}\left(K\right),\end{array}\end{array}$$

$$\begin{array}{}\text{(61)}& {\left(\mathit{\varphi },p_h^{\star }\right)}_K={\left(\mathit{\varphi },p_h\right)}_K,\forall \mathit{\varphi }\in {\mathcal{P}}_l\left(K\right),\end{array}$$

where $\mathrm{0}\le l\le k$. The local problems Eqs. (60)–(61) and Eqs. (58)–(59) are equivalent, with the Lagrange multiplier ψ enforcing orthogonality of test functions in 𝒫_k+1(K) with functions in 𝒫_l(K).

This post-processing method produces a new approximation which superconverges at a rate of k+2 for hybridized mixed methods (Stenberg, 1991; Cockburn et al., 2010a). For the LDG-H method, k+2 superconvergence is achieved when τ=𝒪(1) and τ=𝒪(h), but only k+1 convergence is achieved when $\mathit{\tau }=\mathcal{O}(\mathrm{1}/h)$ (Cockburn et al., 2009b, 2010b). We demonstrate the increased accuracy in computed solutions in Sect. 5.1. An abridged example using Firedrake and Slate to solve the local linear systems is provided in Fig. 3.

Figure 3 Example of local post-processing using Firedrake and Slate. Here, we locally solve the mixed system defined in Eqs. (58)–(59). The corresponding symbolic local tensors are defined in lines 9 and 11. The Slate expression for directly inverting the local system is written in line 12. In line 16, a Slate-generated kernel is produced which solves the resulting linear system in each cell. Since we are not interested in the multiplier, we only return the block corresponding to the new pressure field.

3.3.2 Post-processing of the flux

Our second example illustrates a procedure that uses the numerical flux of an HDG discretization for Eqs. (23)–(26). Within the context of the LDG-H method, we can use the numerical trace in Eq. (52) to produce a vector field that is H(div)-conforming. The technique we outline here follows that of Cockburn et al. (2009b).

Let 𝒯_h be a mesh consisting of simplices. On each cell K∈𝒯_h, we define a new function ${\mathbf{u}}_h^{\star }$ to be the unique element of the local Raviart–Thomas space ${\left[{\mathcal{P}}_k\left(K\right)\right]}^n+\mathbf{x}{\mathcal{P}}_k\left(K\right)$ satisfying

$$\begin{array}{}\text{(62)}& {\displaystyle }& {\displaystyle }{\left(\mathbf{r},{\mathbf{u}}_h^{\star }\right)}_K={\left(\mathbf{r},{\mathbf{u}}_h\right)}_K,\forall \mathbf{r}\in {\left[{\mathcal{P}}_{k-\mathrm{1}}\left(K\right)\right]}^n,\text{(63)}& {\displaystyle }& {\displaystyle }{〈\mathit{\mu },{\mathbf{u}}_h^{\star }\cdot \mathbf{n}〉}_e={〈\mathit{\mu },\widehat{\mathbf{u}}\cdot \mathbf{n}〉}_e,\forall \mathit{\mu }\in {\mathcal{P}}_k\left(e\right)\end{array}$$

for all facets e on ∂K, where $\widehat{\mathbf{u}}$ is the numerical flux defined in Eq. (52). This local problem produces a new velocity ${\mathbf{u}}_h^{\star }$ with the following properties.

1. ${\mathbf{u}}_h^{\star }$ converges at the same rate as u_h for all choices of τ producing a solvable system for Eqs. (54)–(57). However,

2. ${\mathbf{u}}_h^{\star }\in H(\mathrm{div};\mathrm{\Omega })$. That is, ${\left[\left[{\mathbf{u}}_h^{\star }\right]\right]}_e=\mathrm{0},\forall e\in {\mathcal{E}}_h$^∘.

3. Additionally, the divergence of ${\mathbf{u}}_h^{\star }$ convergences at a rate of k+1.

The Firedrake implementation using Slate is similar to the scalar post-processing example (see Fig. 3); the cell-wise linear systems Eqs. (62)–(63) can be expressed in UFL, and therefore the necessary Slate expressions to invert the local systems follows naturally from the set of operations presented in Sect. 2.1. We use the very sensitive parameter dependency in the post-processing methods to validate our software implementation in Zenodo/Tabula-Rasa (2019, “Code-verification”).

4 Static condensation as a preconditioner

Slate enables static condensation approaches to be expressed very concisely. Nonetheless, the application of a particular approach to different variational problems using Slate still requires a certain amount of code repetition. By formulating each form of static condensation as a preconditioner, code can be written once and then applied to any mathematically suitable problem. Rather than writing the static condensation by hand, in many cases, it is sufficient to just select the appropriate, Slate-based, preconditioner.

For context, it is helpful to frame the problem in the particular context of the solver library: PETSc. Firedrake uses PETSc as its main solver abstraction framework and can provide operator-based preconditioners for solving linear systems as PC objects expressed in Python via petsc4py (Dalcin et al., 2011). For a comprehensive overview on solving linear systems using PETSc, we refer the interested reader to Balay et al. (2019, Sect. 4).

Suppose we wish to solve a linear system: Ax=b. We can think of (left) preconditioning the system in residual form:

$$\begin{array}{}\text{(64)}& \mathbf{r}=\mathbf{r}(\mathbf{A},\mathbf{b})\equiv \mathbf{b}-\mathbf{Ax}=\mathbf{0}\end{array}$$

by an operator 𝒫 (which may not necessarily be linear) as a transformation into an equivalent system of the form

$$\begin{array}{}\text{(65)}& \mathcal{P}\mathbf{r}=\mathcal{P}\mathbf{b}-\mathcal{P}\mathbf{Ax}=\mathbf{0}.\end{array}$$

Given a current iterate x_i the residual at the ith iteration is simply ${\mathbf{r}}_i\equiv \mathbf{b}-{\mathbf{Ax}}_i$, and 𝒫 acts on the residual to produce an approximation to the error ${\mathit{ϵ}}_i\equiv \mathbf{x}-{\mathbf{x}}_i$. If 𝒫 is an application of an exact inverse, the residual is converted into an exact (up to numerical roundoff) error.

We will denote the application of a particular Krylov subspace method (KSP) for the linear system Eq. (64) as 𝒦_x(r(A,b)). Upon preconditioning the system via 𝒫 as in Eq. (65), we write

$$\begin{array}{}\text{(66)}& {\mathcal{K}}_{\mathbf{x}}\left(\mathcal{P}\mathbf{r}\right(\mathbf{A},\mathbf{b}\left)\right).\end{array}$$

If Eq. (66) is solved directly via $\mathcal{P}={\mathbf{A}}^{-\mathrm{1}}$, then $\mathcal{P}\mathbf{r}(\mathbf{A},\mathbf{b})={\mathbf{A}}^{-\mathrm{1}}\mathbf{b}-\mathbf{x}$. So Eq. (66) then becomes ${\mathcal{K}}_x\left(\mathbf{r}\right(\mathbf{I},{\mathbf{A}}^{-\mathrm{1}}\mathbf{b}\left)\right)$, producing the exact solution of Eq. (64) in a single iteration of 𝒦. Having established notation, we now present our implementation of static condensation via Slate by defining the appropriate operator, 𝒫.

4.1 Interfacing with PETSc via custom preconditioners

The implementation of preconditioners for the systems considered in this paper requires the manipulation not of assembled matrices but rather their symbolic representation. To do this, we use the preconditioning infrastructure developed by Kirby and Mitchell (2018), which gives preconditioners written in Python access to the symbolic problem description. In Firedrake, this means all derived preconditioners have direct access to the UFL representation of the PDE system. From this mathematical specification, we manipulate this appropriately via Slate and provide operators assembled from Slate expressions to PETSc for further algebraic preconditioning. Using this approach, we have developed a static condensation interface for the hybridization of H(div)×L² mixed problems and a generic interface for statically condensing finite element systems. The advantage of writing even the latter as a preconditioner is the ability to switch out the solution scheme for the system, even when nested inside a larger set of coupled equations or nonlinear solver (Newton-based methods) at runtime.

4.1.1 A static condensation interface for hybridization

As discussed in Sect. 3.1 and 3.2, one of the main advantages of using a hybridizable variant of a DG or mixed method is that such systems permit the use of cell-wise static condensation. To facilitate this, we provide a PETSc PC static condensation interface: firedrake.SCPC. This preconditioner takes the discretized system as in Eq. (44) and performs the local elimination and recovery procedures. Slate expressions are generated from the underlying UFL problem description.

More precisely, the incoming system has the form

$$\begin{array}{}\text{(67)}& \left[\begin{array}{ll}{\mathbf{A}}_{e,e}& {\mathbf{A}}_{e,c}\\ {\mathbf{A}}_{c,e}& {\mathbf{A}}_{c,c}\end{array}\right]\left\{\begin{array}{l}{\mathbf{X}}_e\\ {\mathbf{X}}_c\end{array}\right\}=\left\{\begin{array}{l}{\mathbf{R}}_e\\ {\mathbf{R}}_c\end{array}\right\},\end{array}$$

where X_e is the vector of unknowns to be eliminated, X_c is the vector of unknowns for the condensed field, and R_e and R_c are the incoming right-hand sides. The partitioning in Eq. (67) is determined by the solver option: pc_sc_eliminate_fields. Field indices are provided in the same way one configures solver options to PETSc. These indices determine which field(s) to statically condense into. For example, on a three-field problem (with indices 0, 1, and 2), setting -pc_sc_eliminate_fields 0,1 will configure firedrake.SCPC to cell-wise eliminate fields 0 and 1; the resulting condensed system is associated with field 2.

The firedrake.SCPC preconditioner can be interpreted as a Schur complement method for Eq. (67) of the form

$$\begin{array}{}\text{(68)}& \mathcal{P}=\left[\begin{array}{ll}\mathbf{I}& -{\mathbf{A}}_{e,e}^{-\mathrm{1}}{\mathbf{A}}_{e,c}\\ \mathbf{0}& \mathbf{I}\end{array}\right]\left[\begin{array}{ll}{\mathbf{A}}_{e,e}^{-\mathrm{1}}& \mathbf{0}\\ \mathbf{0}& {\mathbf{S}}^{-\mathrm{1}}\end{array}\right]\left[\begin{array}{ll}\mathbf{I}& \mathbf{0}\\ -{\mathbf{A}}_{c,e}{\mathbf{A}}_{e,e}^{-\mathrm{1}}& \mathbf{I}\end{array}\right],\end{array}$$

where $\mathbf{S}={\mathbf{A}}_{c,c}-{\mathbf{A}}_{c,e}{\mathbf{A}}_{e,e}^{-\mathrm{1}}{\mathbf{A}}_{e,c}$ is the Schur complement operator for the X_c system. The distinction here from block preconditioners via the PETSc fieldsplit option (Brown et al., 2012), for example, is that 𝒫 does not require global actions; by design ${\mathbf{A}}_{e,e}^{-\mathrm{1}}$ can be inverted locally and S is sparse. As a result, S can be assembled or applied exactly, up to numerical roundoff, via Slate-generated kernels.

In practice, the only globally coupled system requiring iterative inversion is S:

$$\begin{array}{}\text{(69)}& {\mathcal{K}}_{{\mathbf{X}}_c}\left({\mathcal{P}}_{\mathrm{1}}r\left(\mathbf{S},{\mathbf{R}}_s\right)\right),\end{array}$$

where ${\mathbf{R}}_s={\mathbf{R}}_c-{\mathbf{A}}_{c,e}{\mathbf{A}}_{e,e}^{-\mathrm{1}}{\mathbf{R}}_e$ is the condensed right-hand side and 𝒫₁ is a preconditioner for S. Once X_c is computed, X_e is reconstructed by inverting the system ${\mathbf{X}}_e={\mathbf{A}}_{e,e}^{-\mathrm{1}}\left({\mathbf{R}}_c-{\mathbf{A}}_{e,c}{\mathbf{X}}_c\right)$ cell-wise.

By construction, this preconditioner is suitable for both hybridized mixed and HDG discretizations. It can also be used within other contexts, such as the static condensation of continuous Galerkin discretizations (Guyan, 1965; Irons, 1965) or primal-hybrid methods (Devloo et al., 2018). As with any PETSc preconditioner, solver options can be specified for inverting S via the appropriate options prefix (condensed_field). The resulting KSP for Eq. (69) is compatible with existing solvers and external packages provided through the PETSc library. This allows users to experiment with a direct method and then switch to a more parallel-efficient iterative solver without changing the core application code.

4.1.2 Preconditioning mixed methods via hybridization

The preconditioner firedrake.HybridizationPC expands on the previous one, this time taking an H(div)×L² system and automatically forming the hybridizable problem. This is accomplished through manipulating the UFL objects representing the discretized PDE. This includes replacing argument spaces with their discontinuous counterparts, introducing test functions on an appropriate trace space, and providing operators assembled from Slate expressions in a similar manner as described in Sect. 4.1.1.

More precisely, let AX=R be the incoming mixed saddle point problem, where $\mathbf{R}={\begin{array}{ll}{\mathbf{R}}_U& {\mathbf{R}}_P\end{array}}^T$, $\mathbf{X}={\begin{array}{ll}\mathbf{U}& \mathbf{P}\end{array}}^T$, and U and P are the velocity and scalar unknowns, respectively. Then this preconditioner replaces AX=R with the extended problem:

$$\begin{array}{}\text{(70)}& \left[\begin{array}{ll}\widehat{\mathbf{A}}& {\mathbf{C}}^T\\ \mathbf{C}& \mathbf{0}\end{array}\right]\left\{\begin{array}{l}\widehat{\mathbf{X}}\\ \mathbf{\Lambda }\end{array}\right\}=\left\{\begin{array}{l}\widehat{\mathbf{R}}\\ {\mathbf{R}}_g\end{array}\right\}\end{array}$$

where Λ are the Lagrange multipliers, $\widehat{\mathbf{R}}={\begin{array}{ll}{\widehat{\mathbf{R}}}_U& {\mathbf{R}}_P\end{array}}^T$, ${\widehat{\mathbf{R}}}_U$ and R_P are the right-hand sides for the flux and scalar equations, respectively, and $\widehat{\cdot }$ indicates modified matrices and co-vectors with discontinuous functions. Here, $\widehat{\mathbf{X}}={\begin{array}{ll}{\mathbf{U}}^d& \mathbf{P}\end{array}}^T$ are the hybridizable (discontinuous) unknowns to be determined and $\mathbf{C}\widehat{\mathbf{X}}={\mathbf{R}}_g$ is the matrix representation of the transmission condition for the hybridizable mixed method (see Eq. 43).

The application of firedrake.HybridizationPC can be interpreted as the Schur complement reduction of Eq. (70):

$$\begin{array}{}\text{(71)}& \widehat{\mathcal{P}}=\left[\begin{array}{ll}\mathbf{I}& -{\widehat{\mathbf{A}}}^{-\mathrm{1}}{\mathbf{C}}^T\\ \mathbf{0}& \mathbf{I}\end{array}\right]\left[\begin{array}{ll}{\widehat{\mathbf{A}}}^{-\mathrm{1}}& \mathbf{0}\\ \mathbf{0}& {\mathbf{S}}^{-\mathrm{1}}\end{array}\right]\left[\begin{array}{ll}\mathbf{I}& \mathbf{0}\\ -\mathbf{C}{\widehat{\mathbf{A}}}^{-\mathrm{1}}& \mathbf{I}\end{array}\right],\end{array}$$

where S is the Schur complement matrix $\mathbf{S}=-\mathbf{C}{\widehat{\mathbf{A}}}^{-\mathrm{1}}{\mathbf{C}}^T$. As before, a single global system for Λ can be assembled cell-wise using Slate-generated kernels. Configuring the solver for inverting S is done via the PETSc options prefix: -hybridization. The recovery of U^d and P happens in the same manner as firedrake.SCPC.

Since the hybridizable flux solution is constructed in the broken H(div) space $U_h^d$, we must project the computed solution into U_h⊂H(div). This can be done cheaply via local facet averaging. The resulting solution is then updated via U←Π_divU^d, where ${\mathrm{\Pi }}_{\mathrm{div}}:U_h^d\to U_h$ is a projection operator. This ensures that the residual for the original mixed problem is properly evaluated to test for solver convergence. With $\widehat{\mathcal{P}}$ as in Eq. (71), the preconditioning operator for the original system AX=R then has the form

$$\begin{array}{}\text{(72)}& \mathcal{P}=\mathbf{\Pi }\widehat{\mathcal{P}}{\mathbf{\Pi }}^T,\mathbf{\Pi }=\left[\begin{array}{lll}{\mathrm{\Pi }}_{\mathrm{div}}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \mathbf{I}& \mathbf{0}\end{array}\right].\end{array}$$

We note here that assembly of the right-hand side for the Λ system requires special attention. Firstly, when Neumann conditions are present, then R_g is not necessarily 0. Since the hybridization preconditioner has access to the entire Python context (which includes a list of boundary conditions and the spaces in which they are applied), surface integrals on the exterior boundary are added where appropriate and incorporated into the generated Slate expressions. A more subtle issue that requires extra care is the incoming right-hand side tested in the H(div) space U_h.

The situation we are given is that we have R_U=R_U(w) for w∈U_h but require ${\widehat{\mathbf{R}}}_U\left({\mathbf{w}}^d\right)$ for ${\mathbf{w}}^d\in U_h^d$. For consistency, we also require for any w∈U_h that

$$\begin{array}{}\text{(73)}& {\widehat{\mathbf{R}}}_U\left(\mathbf{w}\right)={\mathbf{R}}_U\left(\mathbf{w}\right).\end{array}$$

We can construct such an ${\widehat{\mathbf{R}}}_U$ satisfying Eq. (73) in the following way. By construction, we have for each basis function Ψ_i∈U_h

$$\begin{array}{}\text{(74)}& {\mathbf{\Psi }}_i=\left\{\begin{array}{ll}{\mathbf{\Psi }}_i^d& {\mathbf{\Psi }}_i\text{associated with an exterior},\\ & \text{facet node}\\ {\mathbf{\Psi }}_i^{d,+}+{\mathbf{\Psi }}_i^{d,-}& {\mathbf{\Psi }}_i\text{associated with an interior},\\ & \text{facet node}\\ {\mathbf{\Psi }}_i^d& {\mathbf{\Psi }}_i\text{associated with a cell},\\ & \text{interior node}\\ ,& \end{array}\right.\end{array}$$

where ${\mathbf{\Psi }}_i^d,{\mathbf{\Psi }}_i^{d,\pm }\in U_h^d$, and ${\mathbf{\Psi }}_i^{d,\pm }$ are basis functions corresponding to the positive and negative restrictions associated with the ith facet node.⁴ We then define our broken right-hand side via the local definition:

$$\begin{array}{}\text{(75)}& {\widehat{\mathbf{R}}}_U\left({\mathbf{\Psi }}_i^d\right)={\displaystyle \frac{{\mathbf{R}}_U\left({\mathbf{\Psi }}_i\right)}{N_i}},\end{array}$$

where N_i is the number of cells that the degree of freedom corresponding to Ψ_i∈U_h is topologically associated with. Using Eq. (74), Eq. (75), and the fact that R_U is linear in its argument, we can verify that our construction of ${\widehat{\mathbf{R}}}_U$ satisfies Eq. (73).

5 Numerical studies

We now present results utilizing the Slate DSL and our static condensation preconditioners for a set of test problems. Since we are using the interfaces outlined in Sect. 4, Slate is accessed indirectly and requires no manually written solver code for hybridization or static condensation or local recovery. All parallel results were obtained on a single fully loaded compute node of dual-socket Intel E5-2630v4 (Xeon) processors with 2×10 cores (2 threads per core) running at 2.2GHz. In order to avoid potential memory effects due to the operating system migrating processes between sockets, we pin MPI processes to cores.

The verification of the generated code is performed using parameter-sensitive convergence tests. The study consists of running a variety of discretizations spanning the methods outlined in Sect. 3. Details and numerical results are made public and can be viewed in Zenodo/Tabula-Rasa (2019) (see “Code and data availability”). All results are in full agreement with the theory.

5.1 HDG method for a three-dimensional elliptic equation

In this section, we take a closer look at the LDG-H method for the model elliptic equation (sign-definite Helmholtz):

$$\begin{array}{}\text{(76)}& {\displaystyle }& {\displaystyle }-\mathrm{\nabla }\cdot \mathrm{\nabla }p+p=f,\text{ in }\mathrm{\Omega }={\left[\mathrm{0},\mathrm{1}\right]}^{\mathrm{3}},\text{(77)}& {\displaystyle }& {\displaystyle }p=g,\text{ on }\partial \mathrm{\Omega },\end{array}$$

where f and g are chosen such that the analytic solution is p=exp {sin (πx)sin (πy)sin (πz)}. We use a regular mesh consisting of 6⋅N³ tetrahedral elements ($N\in \mathit{\{}\mathrm{4},\mathrm{8},\mathrm{16},\mathrm{32},\mathrm{64}\mathit{\}}$). First, we reformulate Eqs. (76)–(77) as the mixed problem:

$$\begin{array}{}\text{(78)}& {\displaystyle }& {\displaystyle }\mathbf{u}+\mathrm{\nabla }p=\mathrm{0},\text{(79)}& {\displaystyle }& {\displaystyle }\mathrm{\nabla }\cdot \mathbf{u}+p=f,\text{(80)}& {\displaystyle }& {\displaystyle }p=g,\text{ on }\partial \mathrm{\Omega }.\end{array}$$

We start with linear polynomial approximations, up to cubic, for the LDG-H discretization of Eqs. (78)–(80). Additionally, we compute a post-processed scalar approximation $p_h^{\star }$ of the HDG solution. This raises the approximation order of the computed solution by an additional degree. In all numerical studies here, we set the HDG parameter τ=1. All results were computed in parallel, utilizing a single compute node (described previously).

A continuous Galerkin (CG) discretization of the primal problem Eqs. (76)–(77) serves as a reference for this experiment. Due to the superconvergence in the post-processed solution for the HDG method, we use CG discretizations of polynomial orders 2, 3, and 4. This takes into account the enhanced accuracy of the HDG solution, despite being initially computed as a lower-order approximation. We therefore expect both methods to produce equally accurate solutions to the model problem.

Our aim here is not to compare the performance of HDG and CG, which has been investigated elsewhere (for example, see Kirby et al., 2012; Yakovlev et al., 2016). Instead, we provide a reference that the reader might be more familiar with in order to evaluate whether our software framework produces a sufficiently performant HDG implementation relative to what might be expected.

To invert the CG system, we use a conjugate gradient solver with Hypre's BoomerAMG implementation of algebraic multigrid (AMG) as a preconditioner (Falgout et al., 2006). For the HDG method, we use the preconditioner described in Sect. 4.1.1 and the same solver setup as the CG method for the trace system. While the trace operator is indeed symmetric and positive-definite, one should keep in mind that conclusions regarding the performance of off-the-shelf AMG packages on the HDG trace system are still relatively unclear. As a result, efforts on developing more efficient multigrid strategies are a topic of ongoing interest (Cockburn et al., 2014; Kronbichler and Wall, 2018).

To avoid over-solving, we iterate to a relative tolerance such that the discretization error is minimal for a given mesh. In other words, the solvers are configured to terminate when there is no further reduction in the L² error of the computed solution compared with the analytic solution. This means we are not iterating to a fixed solver tolerance across all mesh resolutions. Therefore, we can expect the total number of Krylov iterations (for both the CG and HDG methods) to increase as the mesh resolution becomes finer. The rationale behind this approach is to directly compare the execution time to solve for the best possible approximation to the solution given a fixed resolution.

5.1.1 Error versus execution time

The total execution time is recorded for the CG and HDG solvers, which includes the setup time for the AMG preconditioner, matrix assembly, and the time to solution for the Krylov method. In the HDG case, we include all setup costs, the time spent building the Schur complement for the traces, local recovery of the scalar and flux approximations, and post-processing. The L² error against execution time and Krylov iterations to reach discretization error for each mesh are summarized in Fig. 4.

Figure 4 Comparison of continuous Galerkin and LDG-H solvers for the model three-dimensional positive-definite Helmholtz equation. Panel (a): a log–log plot showing the error against execution time for the CG and HDG with post-processing (τ=1) methods. Panel (b): a log–linear plot showing Krylov iterations of the AMG-preconditioned conjugate gradient algorithm (to reach discretization error) against number of cells.

The HDG method of order k−1 (HDG_k−1) with post-processing, as expected, produces a solution which is as accurate as the CG method of order k (CG_k). While the full HDG system is never explicitly assembled, the larger execution time is a result of several factors. The primary factor is that the total number of trace unknowns for the HDG₁, HDG₂, and HDG₃ discretizations is roughly 4, 3, and 2 times larger, respectively, than the corresponding number of CG unknowns. Therefore, each iteration is more expensive. We also observe that the trace system requires more Krylov iterations to reach discretization error, which appears to improve relative to the CG method as the approximation order increases. Further analysis on a multigrid methods for HDG systems is required to draw further conclusions. The main computational bottleneck in HDG methods is the global linear solver. We therefore expect our implementation to be dominated by the cost associated with inverting the trace operator. If one considers just the time-to-solution, the CG method is clearly ahead of the HDG method. However, the superior scaling, local conservation, and stabilization properties of the HDG method make it a particularly appealing choice for fluid dynamics applications (Yakovlev et al., 2016; Kronbichler and Wall, 2018). Therefore, the development of good preconditioning strategies for the HDG method is critical for its competitive use.

5.1.2 Breakdown of solver time

The HDG method requires many more degrees of freedom than CG or primal DG methods. This is largely due to the fact that the HDG method simultaneously approximates the primal solution and its velocity. The global matrix for the traces is significantly larger than the one for the CG system at low polynomial order. The execution time for HDG is then compounded by a more expensive global solve.

Figure 5 Breakdown of the CG_k and HDG_k−1 execution times on a 6⋅64³ simplicial mesh.

Table 1 Breakdown of the raw timings for the HDG_k−1 (τ=1) and CG_k methods; k=2, 3, and 4. Each method corresponds to a mesh size N=64 on a fully loaded compute node.

Figure 5 displays a breakdown of total execution times on a simplicial mesh consisting of 1.5 million elements. The execution times have been normalized by the CG total time, showing that the HDG method is roughly 3 times more expensive than the CG method. This is expected given the larger degree-of-freedom count and expensive global solve. The raw numerical breakdown of the HDG and CG solvers are shown in Table 1. We isolate each component of the HDG method contributing to the total execution time. Local operations include static condensation (trace operator assembly), forward elimination (right-hand-side assembly for the trace system), backwards substitution to recover the scalar and velocity unknowns, and local post-processing of the scalar solution. For all k, our HDG implementation is solver-dominated as expected.

Both trace operator and right-hand-side assembly are dominated by the costs of inverting a local square mixed matrix coupling the scalar and velocity unknowns, which is performed directly via an LU factorization. This is also the case for backwards substitution. They should all therefore be of the same magnitude in time spent. We observe that this is the case across all degrees, with times ranging between approximately 6 % and 11 % of total execution time for assembling the condensed system. Back-substitution takes roughly the same time as the static condensation and forward elimination stages (approximately 12 % of execution time on average). Finally, the additional cost of post-processing accrues negligible time (roughly 2 % of execution time across all degrees). This is a small cost for an increase in order of accuracy.

We note that the caching of local tensors does not occur. Each pass to perform the local eliminations and backwards reconstructions rebuilds the local element tensors. It is not clear at this time whether the performance gained from avoiding rebuilding the local operators will offset the memory costs of storing the local matrices. Moreover, in time-dependent problems where the operators may contain state-dependent variables, rebuilding local matrices will be necessary in each time step regardless.

5.2 Hybridizable mixed methods for the shallow water equations

A primary motivator for our interest in hybridizable methods revolves around developing efficient solvers for problems in geophysical flows. In this section, we present some results integrating the nonlinear, rotating shallow water equations on the sphere using test case 5 (flow past an isolated mountain) from Williamson et al. (1992). For our discretization approach, we use the framework of compatible finite elements (Cotter and Shipton, 2012; Cotter and Thuburn, 2014).

The model equations we consider are the vector-invariant rotating nonlinear shallow water system defined on a two-dimensional spherical surface Ω embedded in ℝ³:

$$\begin{array}{}\text{(81)}& {\displaystyle \frac{\partial \mathbf{u}}{\partial t}}+\left({\mathrm{\nabla }}^{⟂}\cdot \mathbf{u}+f\right){\mathbf{u}}^{⟂}+\mathrm{\nabla }\left(g\left(D+b\right)+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}|\mathbf{u}{|}^{\mathrm{2}}\right)=\mathrm{0},\end{array}$$

$$\begin{array}{}\text{(82)}& {\displaystyle \frac{\partial D}{\partial t}}+\mathrm{\nabla }\cdot \left(\mathbf{u}D\right)=\mathrm{0},\end{array}$$

where u is the fluid velocity, D is the depth field, f is the Coriolis parameter, g is the acceleration due to gravity, b is the bottom topography, and $(\cdot {)}^{⟂}\equiv \widehat{\mathbf{z}}\times \cdot$, with $\widehat{\mathbf{z}}$ being the unit normal to the surface Ω. After discretizing in time and space using a semi-implicit scheme and Picard linearization, following Natale et al. (2016), we must solve a sequence of the saddle point system at each time step of the form

$$\begin{array}{}\text{(83)}& \left[\begin{array}{ll}\mathbf{A}& -g\frac{\mathrm{\Delta }t}{\mathrm{2}}{\mathbf{B}}^T\\ H\frac{\mathrm{\Delta }t}{\mathrm{2}}\mathbf{B}& \mathbf{M}\end{array}\right]\left\{\begin{array}{l}\mathrm{\Delta }\mathbf{U}\\ \mathrm{\Delta }\mathbf{D}\end{array}\right\}=\left\{\begin{array}{l}{\mathbf{R}}_{\mathbf{u}}\\ {\mathbf{R}}_D\end{array}\right\}.\end{array}$$

See Appendix A for a complete description of the entire discretization strategy. The system Eq. (83) is the matrix equation corresponding to the linearized equations in Eqs. (A3)–(A4).

The Picard updates ΔU and ΔD are sought in the mixed finite element spaces U_h⊂H(div) and V_h⊂L², respectively. Stable mixed finite element pairings correspond to the well-known RT and BDM mixed methods, such as ${\text{RT}}_k\times {\text{DG}}_{k-\mathrm{1}}$ or ${\text{BDM}}_k\times {\text{DG}}_{k-\mathrm{1}}$. These also fall within the set of compatible mixed spaces ideal for geophysical fluid dynamics (Cotter and Shipton, 2012; Natale et al., 2016; Melvin et al., 2019). In particular, the lowest-order RT method (RT₁×DG₀) on a structured quadrilateral grid (such as the latitude–longitude grid used by many operational dynamical cores) is analogous to the Arakawa C-grid finite difference discretization.

In staggered finite difference models, the standard approach for solving Eq. (83) is to neglect the Coriolis term and eliminate the velocity unknown ΔU to obtain a discrete elliptic equation for ΔD, where smoothers like Richardson iterations or relaxation methods are convergent. This is more problematic in the compatible finite element framework, since A has a dense inverse. Instead, we use the preconditioner described in Sect. 4.1.2 to form the equivalent hybridizable formulation, where both ΔU and ΔD are eliminated locally to produce a sparse elliptic equation for the Lagrange multipliers.

5.2.1 Atmospheric flow over a mountain

As a test problem, we solve test case 5 of Williamson et al. (1992), on the surface of a sphere with radius R=6371km. We refer the reader to Cotter and Shipton (2012) and Shipton et al. (2018) for a more comprehensive study on mixed finite elements for shallow water systems of this type. We use the mixed finite element pairs (RT₁,DG₀) (lowest-order RT method) and (BDM₂,DG₁) (next-to-lowest-order BDM method) for the velocity and depth spaces. A mesh of the sphere is generated from seven refinements of an icosahedron, resulting in a triangulation 𝒯_h consisting of 327 680 elements in total. The grid information for both mixed methods is summarized in Table 2.

Table 2 The number of unknowns to be determined are summarized for each compatible finite element method. Resolution is the same for both methods.

We run for a total of 25 time steps, with a fixed number of four Picard iterations in each time step. We compare the overall simulation time using two different solver configurations for the implicit linear system. First, we use a flexible variant of the generalized minimal residual method (GMRES) ⁵ acting on the system Eq. (83) with an approximate Schur complement preconditioner:

$$\begin{array}{}\text{(84)}& {\mathcal{P}}_{\text{SC}}=\left[\begin{array}{ll}\mathbf{I}& g\frac{\mathrm{\Delta }t}{\mathrm{2}}{\mathbf{A}}^{-\mathrm{1}}{\mathbf{B}}^T\\ \mathbf{0}& \mathbf{I}\end{array}\right]\left[\begin{array}{ll}{\mathbf{A}}^{-\mathrm{1}}& \mathbf{0}\\ \mathbf{0}& {\stackrel{\mathrm{̃}}{\mathbf{S}}}^{-\mathrm{1}}\end{array}\right]\end{array}$$

$$\left[\begin{array}{ll}\mathbf{I}& \mathbf{0}\\ -H\frac{\mathrm{\Delta }t}{\mathrm{2}}{\mathbf{BA}}^{-\mathrm{1}}& \mathbf{I}\end{array}\right],$$

where $\stackrel{\mathrm{̃}}{\mathbf{S}}=\mathbf{M}+gH\frac{\mathrm{\Delta }{t}^{\mathrm{2}}}{\mathrm{4}}\mathbf{B}\mathrm{diag}(\mathbf{A}{)}^{-\mathrm{1}}{\mathbf{B}}^T$ and diag(A) is a diagonal approximation to the velocity mass matrix (plus the addition of a Coriolis matrix). The Schur complement system is inverted via GMRES due to the asymmetry from the Coriolis term, with the inverse of $\stackrel{\mathrm{̃}}{\mathbf{S}}$ as the preconditioning operator. The sparse approximation $\stackrel{\mathrm{̃}}{\mathbf{S}}$ is inverted using PETSc's smoothed aggregation multigrid (GAMG). The Krylov method is set to terminate once the preconditioned residual norm is reduced by a factor of 10⁸. A⁻¹ is computed approximately using a single application of incomplete LU (ILU) with zero fill-in.

Next, we use only the application of our hybridization preconditioner (no outer Krylov method), which replaces the original linearized mixed system with its hybridizable equivalent. After hybridization, we have the following extended problem for the Picard updates: find $(\mathrm{\Delta }{\mathbf{u}}_h^d,\mathrm{\Delta }{D}_h,{\mathit{\lambda }}_h)\in U_h^d\times V_h\times M_h$ satisfying

$$\begin{array}{}\text{(85)}& \begin{array}{rl}& {\left(\mathbf{w},\mathrm{\Delta }{\mathbf{u}}_h^d\right)}_{{\mathcal{T}}_h}+{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{2}}}{\left(\mathbf{w},f{\left(\mathrm{\Delta }{\mathbf{u}}_h^d\right)}^{⟂}\right)}_{{\mathcal{T}}_h}\\ & -{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{2}}}{\left(\mathrm{\nabla }\cdot \mathbf{w},g\mathrm{\Delta }{D}_h\right)}_{{\mathcal{T}}_h}+{〈\left[\left[\mathbf{w}\right]\right],{\mathit{\lambda }}_h〉}_{\partial {\mathcal{T}}_h}\\ & ={\widehat{R}}_{\mathbf{u}},\forall \mathbf{w}\in U_h^d,\end{array}\end{array}$$

$$\begin{array}{}\text{(86)}& {\left(\mathit{\varphi },\mathrm{\Delta }{D}_h\right)}_{{\mathcal{T}}_h}+{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{2}}}{\left(\mathit{\varphi },H\mathrm{\nabla }\cdot \mathrm{\Delta }{\mathbf{u}}_h^d\right)}_{{\mathcal{T}}_h}=R_D,\forall \mathit{\varphi }\in V_h,\end{array}$$

$$\begin{array}{}\text{(87)}& {〈\mathit{\gamma },\left[\left[\mathrm{\Delta }{\mathbf{u}}_h^d\right]\right]〉}_{\partial {\mathcal{T}}_h}=\mathrm{0},\forall \mathit{\gamma }\in M_h.\end{array}$$

Note that the space M_h is chosen such that the trace functions, when restricted to a facet $e\in \partial {\mathcal{T}}_h$, are in the same polynomial space as $\mathrm{\Delta }{\mathbf{u}}_h\cdot \mathbf{n}{|}_e$. Moreover, it can be shown that the Lagrange multiplier λ_h is an approximation to the depth unknown ΔtgΔD∕2 restricted to ∂𝒯_h.

The resulting three-field problem in Eqs. (85)–(87) produces the following matrix equation:

$$\begin{array}{}\text{(88)}& \left[\begin{array}{ll}\widehat{\mathcal{A}}& {\mathbf{C}}^T\\ \mathbf{C}& \mathbf{0}\end{array}\right]\left\{\begin{array}{l}\mathrm{\Delta }\mathbf{X}\\ \mathbf{\Lambda }\end{array}\right\}=\left\{\begin{array}{l}{\widehat{\mathbf{R}}}_{\mathrm{\Delta }X}\\ \mathbf{0}\end{array}\right\},\end{array}$$

where $\widehat{\mathcal{A}}$ is the discontinuous operator coupling $\mathrm{\Delta }\mathbf{X}={\begin{array}{ll}\mathrm{\Delta }{\mathbf{U}}^d& \mathrm{\Delta }\mathbf{D}\end{array}}^T$ and ${\mathbf{R}}_{\mathrm{\Delta }X}={\begin{array}{ll}{\widehat{\mathbf{R}}}_{\mathbf{u}}& {\mathbf{R}}_D\end{array}}^T$ are the problem residuals. An exact Schur complement factorization is performed on Eq. (88), using Slate to generate the local elimination kernels. We use the same set of solver options for the inversion of $\stackrel{\mathrm{̃}}{\mathbf{S}}$ in Eq. (84) to invert the Lagrange multiplier system. The increments ΔU^d and ΔD are recovered locally, using Slate-generated kernels. Once recovery is complete, ΔU^d is projected back into the conforming H(div) finite element space via ΔU←Π_divΔU^d. Based on the discussion in Sect. 4.1.2, we apply

$$\begin{array}{}\text{(89)}& {\mathcal{P}}_{\text{hybrid}}=\mathbf{\Pi }\left(\left[\begin{array}{ll}\mathbf{I}& -{\widehat{\mathcal{A}}}^{-\mathrm{1}}{\mathbf{C}}^T\\ \mathbf{0}& \mathbf{I}\end{array}\right]\left[\begin{array}{ll}{\widehat{\mathcal{A}}}^{-\mathrm{1}}& \mathbf{0}\\ \mathbf{0}& {\mathbf{S}}^{-\mathrm{1}}\end{array}\right]\right.\end{array}$$

$$\left.\left[\begin{array}{ll}\mathbf{I}& \mathbf{0}\\ -\mathbf{C}{\widehat{\mathcal{A}}}^{-\mathrm{1}}& \mathbf{I}\end{array}\right]\right){\mathbf{\Pi }}^T.$$

Table 3 Preconditioner solve times for a 25-step run with Δt=100s. These are cumulative times in each stage of the two preconditioners throughout the entire profile run. We display the average iteration count (rounded to the nearest integer) for both the outer and the inner Krylov solvers. The significant speedup when using hybridization is a direct result of eliminating the outermost solver.

Table 4 Breakdown of the cost (average) of a single application of the preconditioned flexible GMRES method and hybridization preconditioner. Hybridization takes approximately the same time per iteration.

Table 3 displays a summary of our findings. The advantages of a hybridizable method versus a mixed method are more clearly realized in this experiment. When using hybridization, we observe a significant reduction in time spent in the implicit solver compared to the approximate Schur complement approach. This is primarily because we have reduced the number of “outer” iterations to zero; the hybridization preconditioner is performing an exact factorization of the global hybridizable system. This is empirically supported when considering per-application solve times. The values reported in Table 4 show the average cost of a single outer GMRES iteration (which includes the application of 𝒫_SC) and a single application of 𝒫_hybrid. Hybridization and the approximate Schur complement preconditioner are comparable in terms of average execution time, with hybridization being slightly faster. This further demonstrates that the primary cause for the longer execution time of the latter is directly related to the additional outer iterations induced from using an approximate factorization. In terms of over all time-to-solution, the hybridizable methods are clearly ahead of the original mixed methods.

We also measure the relative reductions in the problem residual of the linear system Eq. (83). Our hybridization preconditioner reduces the residual by a factor of 10⁸ on average, which coincides with the specified relative tolerance for the Krylov method on the trace system. In other words, the reduction in the residual for the trace system translates into an overall reduction in the residual for the mixed system by the same factor.

The test case was run up to day 15 on a coarser resolution (20 480 simplicial cells with Δx≈210km) and a time-step size Δt=500 s. Snapshots of the entire simulation are provided in Fig. 6 using the semi-implicit scheme described in Appendix A. The results we have obtained for days 5, 10, and 15 are comparable to the corresponding results of Nair et al. (2005), Ullrich et al. (2010), and Kang et al. (2020). We refer the reader to Shipton et al. (2018) for further demonstrations of shallow water test cases featuring the use of the hybridization preconditioner described in Sect. 4.1.2.

5.3 Hybridizable methods for a linear Boussinesq model

As a final example, we consider the simplified atmospheric model obtained from a linearization of the compressible Boussinesq equations in a rotating domain:

$$\begin{array}{}\text{(90)}& {\displaystyle }& {\displaystyle \frac{\partial \mathbf{u}}{\partial t}}+\mathrm{2}\mathbf{\Omega }\times \mathbf{u}=-\mathrm{\nabla }p+b\widehat{\mathbf{z}},\text{(91)}& {\displaystyle }& {\displaystyle \frac{\partial p}{\partial t}}=-c^{\mathrm{2}}\mathrm{\nabla }\cdot \mathbf{u},\text{(92)}& {\displaystyle }& {\displaystyle \frac{\partial b}{\partial t}}=-N^{\mathrm{2}}\mathbf{u}\cdot \widehat{\mathbf{z}},\end{array}$$

where u is the fluid velocity, p the pressure, b the buoyancy, Ω the planetary angular rotation vector, c the speed of sound ($\approx \mathrm{343}\text{m}{\text{s}}^{-\mathrm{1}}$), and N the buoyancy frequency ($\approx \mathrm{0.01}{\text{s}}^{-\mathrm{1}}$). Equations (90)–(92) permit fast-moving acoustic waves driven by perturbations in b. This is the model presented in Skamarock and Klemp (1994), which uses a quadratic equation of state to avoid some of the complications of the full compressible Euler equations (the hybridization of which we shall address in future work). We solve these equations subject to the rigid-lid condition $\mathbf{u}\cdot \mathbf{n}=\mathrm{0}$ on all boundaries.

Our domain consists of a spherical annulus, with the mesh constructed from a horizontal “base” mesh of the surface of a sphere of radius R, extruded upwards by a height H_Ω. The vertical discretization is a structured one-dimensional grid, which facilitates the staggering of thermodynamic variables, such as b. We consider two kinds of meshes: one obtained by extruding an icosahedral sphere mesh and another from a cubed sphere.

Since our mesh has a natural tensor product structure, we construct suitable finite element spaces constructed by taking the tensor product of a horizontal space with a vertical space. To ensure our discretization is “compatible,” we use the one- and two-dimensional finite element de Rham complexes: $V_h^{\mathrm{0}}\stackrel{{\partial }_z}{\to }{V}_h^{\mathrm{1}}$ and $U_h^{\mathrm{0}}\stackrel{{\mathrm{\nabla }}^{⟂}}{\to }{U}_h^{\mathrm{1}}\stackrel{\mathrm{\nabla }\cdot }{\to }{U}_h^{\mathrm{2}}$. We can then construct the three-dimensional complex: $W_h^{\mathrm{0}}\stackrel{\nabla }{\to }{W}_h^{\mathrm{1}}\stackrel{\mathrm{\nabla }\times }{\to }{W}_h^{\mathrm{2}}\stackrel{\mathrm{\nabla }\cdot }{\to }{W}_h^{\mathrm{3}}$, where

$$\begin{array}{}\text{(93)}& W_h^{\mathrm{0}}=U_h^{\mathrm{0}}\otimes V_h^{\mathrm{0}},\end{array}$$

$$\begin{array}{}\text{(94)}& \begin{array}{rl}& W_h^{\mathrm{1}}=\mathtt{\text{HCurl}}(U_h^{\mathrm{1}}\otimes V_h^{\mathrm{0}})\oplus \mathtt{\text{HCurl}}(U_h^{\mathrm{0}}\otimes V_h^{\mathrm{1}})=:\\ & W_h^{\mathrm{1},h}\oplus W_h^{\mathrm{1},v},\end{array}\end{array}$$

$$\begin{array}{}\text{(95)}& \begin{array}{rl}& W_h^{\mathrm{2}}=\mathtt{\text{HDiv}}(U_h^{\mathrm{1}}\otimes V_h^{\mathrm{1}})\oplus \mathtt{\text{HDiv}}(U_h^{\mathrm{2}}\otimes V_h^{\mathrm{0}})=:\\ & W_h^{\mathrm{2},h}\oplus W_h^{\mathrm{2},v},\end{array}\end{array}$$

$$\begin{array}{}\text{(96)}& W_h^{\mathrm{3}}=U_h^{\mathrm{2}}\otimes V_h^{\mathrm{1}}.\end{array}$$

Here, HCurl and HDiv denote operators which ensure that the correct Piola transformations are applied when mapping from physical to reference element. We refer the reader to McRae et al. (2016) for an overview of constructing tensor product finite element spaces in Firedrake. For the analysis of compatible finite element discretizations and their relation to the complex Eqs. (93)–(96), we refer the reader to Natale et al. (2016). Each discretization used in this section is constructed from more familiar finite element families, shown in Table 5.

Figure 6 Snapshots (view from the northern pole) from the isolated mountain test case. The surface height (m) at days 5, 10, and 15. The snapshots were generated on a mesh with 20 480 simplicial cells, a BDM₂×DG₁ discretization, and Δt=500 s. The linear system during each Picard cycle was solved using the hybridization preconditioner.

Table 5 Vertical and horizontal spaces for the three-dimensional compatible finite element discretization of the linear Boussinesq model. The RT_k and BDFM_k+1 methods are constructed on triangular prism elements, while the RTCF_k method is defined on extruded quadrilateral elements.

5.3.1 Compatible finite element discretization

A compatible finite element discretization of Eqs. (90)–(92) constructs solutions in the following finite element spaces:

$$\begin{array}{}\text{(97)}& {\mathbf{u}}_h\in {\stackrel{\mathrm{˚}}{W}}_h^{\mathrm{2}},p_h\in W_h^{\mathrm{3}},b_h\in W_h^b,\end{array}$$

where ${\stackrel{\mathrm{˚}}{W}}_h^{\mathrm{2}}$ is the subspace of $W_h^{\mathrm{2}}\subset H\left(\mathrm{div}\right)$ whose functions w satisfy $\mathbf{w}\cdot \mathbf{n}=\mathrm{0}$ on ∂Ω, $W_h^{\mathrm{3}}\subset L^{\mathrm{2}}$, and $W_h^b\equiv U_h^{\mathrm{2}}\otimes V_h^{\mathrm{0}}$. Note that $W_h^b$ is just the scalar version of the vertical velocity space.⁶ That is, $W_h^b$ and $W_h^{\mathrm{2},v}$ have the same number of degrees of freedom but differ in how they are pulled back to the reference element.

To obtain the discrete system, we simply multiply Eqs. (90)–(92) by test functions $\mathbf{w}\in {\stackrel{\mathrm{˚}}{W}}_h^{\mathrm{2}}$, $\mathit{\varphi }\in W_h^{\mathrm{3}}$, and $\mathit{\eta }\in W_h^b$ and integrate by parts. We introduce the increments $\mathit{\delta }{\mathbf{u}}_h\equiv {\mathbf{u}}_h^{n+\mathrm{1}}-{\mathbf{u}}_h^n$, and set ${\mathbf{u}}_{\mathrm{0}}\equiv {\mathbf{u}}_h^n$ (similarly for δp_h, p₀, δb_h, and b₀). Using an implicit midpoint rule discretization, we need to solve the following mixed problem at each time step: find $\mathit{\delta }{\mathbf{u}}_h\in {\stackrel{\mathrm{˚}}{W}}_h^{\mathrm{2}}$, $\mathit{\delta }{p}_h\in W_h^{\mathrm{3}}$ and $\mathit{\delta }{b}_h\in W_h^b$ such that

$$\begin{array}{}\text{(98)}& \begin{array}{rl}& {\left(\mathbf{w},\mathit{\delta }{\mathbf{u}}_h\right)}_{{\mathcal{T}}_h}+{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{2}}}{\left(\mathbf{w},\mathrm{2}\mathbf{\Omega }\times \mathit{\delta }{\mathbf{u}}_h\right)}_{{\mathcal{T}}_h}-{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{2}}}\\ & {\left(\mathrm{\nabla }\cdot \mathbf{w},\mathit{\delta }{p}_h\right)}_{{\mathcal{T}}_h}-{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{2}}}{\left(\mathbf{w},\mathit{\delta }{b}_h\widehat{\mathbf{z}}\right)}_{{\mathcal{T}}_h}=r_u,\forall \mathbf{w}\in {\stackrel{\mathrm{˚}}{W}}_h^{\mathrm{2}}\end{array}\end{array}$$

$$\begin{array}{}\text{(99)}& {\left(\mathit{\varphi },\mathit{\delta }{p}_h\right)}_{{\mathcal{T}}_h}+{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{2}}}{c}^{\mathrm{2}}{\left(\mathit{\varphi },\mathrm{\nabla }\cdot \mathit{\delta }{\mathbf{u}}_h\right)}_{{\mathcal{T}}_h}=r_p,\forall \mathit{\varphi }\in W_h^{\mathrm{3}},\end{array}$$

$$\begin{array}{}\text{(100)}& {\left(\mathit{\eta },\mathit{\delta }{b}_h\right)}_{{\mathcal{T}}_h}+{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{2}}}{N}^{\mathrm{2}}{\left(\mathit{\eta },\mathit{\delta }{\mathbf{u}}_h\cdot \widehat{\mathbf{z}}\right)}_{{\mathcal{T}}_h}=r_b,\forall \mathit{\eta }\in W_h^b,\end{array}$$

where the residuals are $r_u=-\mathrm{\Delta }t{\left(\mathbf{w},\mathrm{2}\mathbf{\Omega }\times {\mathbf{u}}_{\mathrm{0}}\right)}_{{\mathcal{T}}_h}$, $r_p=-c^{\mathrm{2}}\mathrm{\Delta }t{\left(\mathit{\varphi },\mathrm{\nabla }\cdot {\mathbf{u}}_{\mathrm{0}}\right)}_{{\mathcal{T}}_h}$, and $r_b=-N^{\mathrm{2}}\mathrm{\Delta }t{\left(\mathit{\eta },{\mathbf{u}}_{\mathrm{0}}\cdot \widehat{\mathbf{z}}\right)}_{{\mathcal{T}}_h}$.

The resulting matrix equations have the form

$$\begin{array}{}\text{(101)}& \left[\begin{array}{lll}{\mathbf{A}}_{\mathbf{u}}& -\frac{\mathrm{\Delta }t}{\mathrm{2}}{\mathbf{D}}^T& -\frac{\mathrm{\Delta }t}{\mathrm{2}}{\mathbf{Q}}^T\\ \frac{\mathrm{\Delta }t}{\mathrm{2}}{c}^{\mathrm{2}}\mathbf{D}& {\mathbf{M}}_p& \mathbf{0}\\ \frac{\mathrm{\Delta }t}{\mathrm{2}}{N}^{\mathrm{2}}\mathbf{Q}& \mathbf{0}& {\mathbf{M}}_b\end{array}\right]\left\{\begin{array}{l}\mathbf{U}\\ \mathbf{P}\\ \mathbf{B}\end{array}\right\}=\left\{\begin{array}{l}{\mathbf{R}}_{\mathbf{u}}\\ {\mathbf{R}}_p\\ {\mathbf{R}}_b\end{array}\right\},\end{array}$$

where ${\mathbf{A}}_{\mathbf{u}}={\mathbf{M}}_{\mathbf{u}}+\frac{\mathrm{\Delta }t}{\mathrm{2}}{\mathbf{C}}_{\mathbf{\Omega }}$, C_Ω is the asymmetric matrix associated with the Coriolis term, M_u, M_p, and M_b are mass matrices, D is the weak divergence term, and Q is an operator containing the vertical components of δu_h. In the absence of orography, we can use the point-wise expression for the buoyancy,

$$\begin{array}{}\text{(102)}& \mathit{\delta }{b}_h=r_b-{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{2}}}{N}^{\mathrm{2}}\mathit{\delta }{\mathbf{u}}_h\cdot \widehat{\mathbf{z}},\end{array}$$

and eliminate δb_h from Eq. (101) by substituting Eq. (102) into Eq. (98). This produces the following mixed velocity–pressure system:

$$\begin{array}{}\text{(103)}& \mathcal{A}\left\{\begin{array}{l}\mathbf{U}\\ \mathbf{P}\end{array}\right\}=\left[\begin{array}{ll}{\stackrel{\mathrm{̃}}{\mathbf{A}}}_{\mathbf{u}}& -\frac{\mathrm{\Delta }t}{\mathrm{2}}{\mathbf{D}}^T\\ c^{\mathrm{2}}\frac{\mathrm{\Delta }t}{\mathrm{2}}\mathbf{D}& {\mathbf{M}}_p\end{array}\right]\left\{\begin{array}{l}\mathbf{U}\\ \mathbf{P}\end{array}\right\}=\left\{\begin{array}{l}{\stackrel{\mathrm{̃}}{\mathbf{R}}}_{\mathbf{u}}\\ {\mathbf{R}}_p\end{array}\right\},\end{array}$$

where ${\stackrel{\mathrm{̃}}{\mathbf{A}}}_{\mathbf{u}}={\mathbf{A}}_{\mathbf{u}}+\frac{\mathrm{\Delta }{t}^{\mathrm{2}}}{\mathrm{4}}{N}^{\mathrm{2}}{\mathbf{Q}}^T{\mathbf{M}}_b^{-\mathrm{1}}\mathbf{Q}$ and ${\stackrel{\mathrm{̃}}{\mathbf{R}}}_{\mathbf{u}}={\mathbf{R}}_{\mathbf{u}}+\frac{\mathrm{\Delta }t}{\mathrm{2}}{\mathbf{Q}}^T{\mathbf{M}}_b^{-\mathrm{1}}{\mathbf{R}}_b$ are the modified velocity operator and right-hand side, respectively. Note that in our elimination strategy, ${\stackrel{\mathrm{̃}}{\mathbf{A}}}_{\mathbf{u}}$ corresponds to the bilinear form obtained after eliminating the buoyancy at the equation level:

$$\begin{array}{}\text{(104)}& \begin{array}{rl}& {\stackrel{\mathrm{̃}}{\mathbf{A}}}_{\mathbf{u}}\leftarrow {\left(\mathbf{w},\mathit{\delta }\mathbf{u}\right)}_{{\mathcal{T}}_h}+{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{2}}}{\left(\mathbf{w},\mathrm{2}\mathbf{\Omega }\times \mathit{\delta }\mathbf{u}\right)}_{{\mathcal{T}}_h}+{\displaystyle \frac{\mathrm{\Delta }{t}^{\mathrm{2}}}{\mathrm{4}}}\\ & N^{\mathrm{2}}{\left(\mathbf{w}\cdot \widehat{\mathbf{z}},\mathit{\delta }\mathbf{u}\cdot \widehat{\mathbf{z}}\right)}_{{\mathcal{T}}_h}.\end{array}\end{array}$$

A similar construction holds for ${\stackrel{\mathrm{̃}}{\mathbf{R}}}_{\mathbf{u}}$. Once Eq. (103) is solved, δb_h is reconstructed by solving

$$\begin{array}{}\text{(105)}& {\mathbf{M}}_b\mathbf{B}={\mathbf{R}}_b-{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{2}}}{N}^{\mathrm{2}}\mathbf{QU}.\end{array}$$

Equation (105) can be efficiently inverted using the conjugate gradient method.

5.3.2 Preconditioning the mixed velocity pressure system

The primary difficulty is finding efficient solvers for Eq. (103). This was studied by Mitchell and Müller (2016) within the context of developing a preconditioner which is robust against parameters, like Δt and mesh resolution. However, the implicit treatment of the Coriolis term was not taken into account. We consider two preconditioning strategies.

One strategy proposed by Mitchell and Müller (2016) is to build a preconditioner based on the Schur complement factorization of 𝒜 in Eq. (103):

$$\begin{array}{}\text{(106)}& {\mathcal{A}}^{-\mathrm{1}}=\left[\begin{array}{ll}\mathbf{I}& \frac{\mathrm{\Delta }t}{\mathrm{2}}{\stackrel{\mathrm{̃}}{\mathbf{A}}}_{\mathbf{u}}^{-\mathrm{1}}{\mathbf{D}}^T\\ \mathbf{0}& \mathbf{I}\end{array}\right]\left[\begin{array}{ll}{\stackrel{\mathrm{̃}}{\mathbf{A}}}_{\mathbf{u}}^{-\mathrm{1}}& \mathbf{0}\\ \mathbf{0}& {\mathbf{H}}^{-\mathrm{1}}\end{array}\right]\end{array}$$

$$\left[\begin{array}{ll}\mathbf{I}& \mathbf{0}\\ -c^{\mathrm{2}}\frac{\mathrm{\Delta }t}{\mathrm{2}}\mathbf{D}{\stackrel{\mathrm{̃}}{\mathbf{A}}}_{\mathbf{u}}^{-\mathrm{1}}& \mathbf{I}\end{array}\right],$$

where $\mathbf{H}={\mathbf{M}}_p+c^{\mathrm{2}}\frac{\mathrm{\Delta }{t}^{\mathrm{2}}}{\mathrm{4}}\mathbf{D}{\stackrel{\mathrm{̃}}{\mathbf{A}}}_{\mathbf{u}}^{-\mathrm{1}}{\mathbf{D}}^T$ is the dense pressure Helmholtz operator. Because we have chosen to include the Coriolis term, the operator H is nonsymmetric and has the form

$$\begin{array}{}\text{(107)}& \mathbf{H}={\mathbf{M}}_p+c^{\mathrm{2}}{\displaystyle \frac{\mathrm{\Delta }{t}^{\mathrm{2}}}{\mathrm{4}}}\mathbf{D}{\left({\stackrel{\mathrm{̃}}{\mathbf{M}}}_{\mathbf{u}}+{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{2}}}{\mathbf{C}}_{\mathbf{\Omega }}\right)}^{-\mathrm{1}}{\mathbf{D}}^T,\end{array}$$

where ${\stackrel{\mathrm{̃}}{\mathbf{M}}}_{\mathbf{u}}$ is a modified mass matrix. As Δt increases, the contribution of C_Ω becomes more prominent in H, making sparse approximations of H more challenging. We shall elaborate on this further below when we present the results of our second solver strategy.

Our preferred strategy solves the hybridizable formulation of the system Eq. (103). Let $W_h^{\mathrm{2},d}$ denote the broken version of $W_h^{\mathrm{2}}$ and M_h the space of Lagrange multipliers. Then the hybridizable formulation for the velocity–pressure system reads as follows: find $\mathit{\delta }{\mathbf{u}}_h^d\in W_h^{\mathrm{2},d}$, $\mathit{\delta }{p}_h\in W_h^{\mathrm{3}}$, and λ_h∈M_h such that

$$\begin{array}{}\text{(108)}& \begin{array}{rl}& {\left(\mathbf{w},\mathit{\delta }{\mathbf{u}}_h^d\right)}_{{\mathcal{T}}_h}+{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{2}}}{\left(\mathbf{w},\mathrm{2}\mathbf{\Omega }\times \mathit{\delta }{\mathbf{u}}_h^d\right)}_{{\mathcal{T}}_h}\\ & +{\displaystyle \frac{\mathrm{\Delta }{t}^{\mathrm{2}}}{\mathrm{4}}}{N}^{\mathrm{2}}{\left(\mathbf{w}\cdot \widehat{\mathbf{z}},\mathit{\delta }{\mathbf{u}}_h^d\cdot \widehat{\mathbf{z}}\right)}_{{\mathcal{T}}_h}-{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{2}}}{\left(\mathrm{\nabla }\cdot \mathbf{w},\mathit{\delta }{p}_h\right)}_{{\mathcal{T}}_h}\\ & +{〈\left[\left[\mathbf{w}\right]\right],{\mathit{\lambda }}_h〉}_{\partial {\mathcal{T}}_h}={\stackrel{\mathrm{̃}}{r}}_u,\forall \mathbf{w}\in W_h^{\mathrm{2},d}\end{array}\end{array}$$

$$\begin{array}{}\text{(109)}& {\left(\mathit{\varphi },\mathit{\delta }{p}_h\right)}_{{\mathcal{T}}_h}+{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{2}}}{c}^{\mathrm{2}}{\left(\mathit{\varphi },\mathrm{\nabla }\cdot \mathit{\delta }{\mathbf{u}}_h\right)}_{{\mathcal{T}}_h}=r_p,\forall \mathit{\varphi }\in W_h^{\mathrm{3}},\end{array}$$

$$\begin{array}{}\text{(110)}& {〈\mathit{\gamma },\left[\left[\mathit{\delta }{\mathbf{u}}_h^d\right]\right]〉}_{\partial {\mathcal{T}}_h}=\mathrm{0},\forall \mathit{\gamma }\in M_h.\end{array}$$

The system Eqs. (108)–(110) is automatically formed by the Firedrake preconditioner: firedrake.HybridizationPC. We then locally eliminate the velocity and pressure after hybridization, producing the following condensed problem:

$$\begin{array}{}\text{(111)}& {\mathbf{H}}_{\partial }\mathbf{\Lambda }=\mathbf{E},{\mathbf{H}}_{\partial }=\left[\begin{array}{ll}\mathbf{C}& \mathbf{0}\end{array}\right]{\left[\begin{array}{ll}{\widehat{\stackrel{\mathrm{̃}}{\mathbf{A}}}}_{\mathbf{u}}& -\frac{\mathrm{\Delta }t}{\mathrm{2}}{\widehat{\mathbf{D}}}^T\\ c^{\mathrm{2}}\frac{\mathrm{\Delta }t}{\mathrm{2}}\widehat{\mathbf{D}}& {\mathbf{M}}_p\end{array}\right]}^{-\mathrm{1}}\end{array}$$

$$\left[\begin{array}{l}{\mathbf{C}}^T\\ \mathbf{0}\end{array}\right],E=\left[\begin{array}{ll}\mathbf{C}& \mathbf{0}\end{array}\right]{\left[\begin{array}{ll}{\widehat{\stackrel{\mathrm{̃}}{\mathbf{A}}}}_{\mathbf{u}}& -\frac{\mathrm{\Delta }t}{\mathrm{2}}{\widehat{\mathbf{D}}}^T\\ c^{\mathrm{2}}\frac{\mathrm{\Delta }t}{\mathrm{2}}\widehat{\mathbf{D}}& {\mathbf{M}}_p\end{array}\right]}^{-\mathrm{1}}\left\{\begin{array}{l}{\widehat{R}}_{\mathbf{u}}\\ R_p\end{array}\right\},$$

where $\widehat{\cdot }$ denotes matrices or vectors with discontinuous test and trial functions. The nonsymmetric operator H_∂ is inverted using a preconditioned generalized conjugate residual (GCR) Krylov method, as suggested in Thomas et al. (2003). For our choice of preconditioner, we follow strategies outlined in Elman et al. (2001) and employ an algebraic multigrid method (V cycle) with GMRES (five iterations) smoothers on the coarse levels. The GMRES smoothers are preconditioned with block ILU on each level. For the finest level, block ILU produces a line smoother (necessary for efficient solution on thin domains) when the trace variable nodes are numbered in vertical lines, as is the case in our Firedrake implementation. On the coarser levels, less is known about the properties of ILU under the AMG coarsening strategies, but as we shall see, we observe performance that suggests ILU is still behaving as a line smoother. More discussion on multigrid for nonsymmetric problems can be found in Bramble et al. (1994, 1988) and Mandel (1986). A gravity wave test using our solution strategy and hybridization preconditioner is illustrated in Fig. 7 for a problem on a condensed Earth (radius scaled down by a factor of 125) and 10 km lid.

Figure 7 Buoyancy perturbation (y−z cross section) at t=3600 s from a simple gravity wave test (Δt=100 s). The initial conditions (in lat–long coordinates) for the velocity are a simple solid-body rotation: u=20e_λ, where e_λ is the unit vector pointing in the direction of decreasing longitude. A buoyancy anomaly is defined via $b=\frac{d^{\mathrm{2}}}{d^{\mathrm{2}}+q^{\mathrm{2}}}\sin (\mathit{\pi }z/\mathrm{10}\mathrm{000})$, where $q=R{\cos }^{-\mathrm{1}}\left(\cos \right(\mathit{\varphi }\left)\cos \right(\mathit{\lambda }-{\mathit{\lambda }}_{\mathit{\varphi }}\left)\right)$, d=5000 m, $R=\mathrm{6371}\text{km}/\mathrm{125}$ is the planet radius, and ${\mathit{\lambda }}_{\mathit{\varphi }}=\mathrm{2}/\mathrm{3}$. The equations are discretized using the lowest-order method RTCF₁, with 24 576 quadrilateral cells in the horizontal and 64 extrusion levels. The velocity–pressure system is solved using hybridization.

5.3.3 Robustness against acoustic Courant number with implicit Coriolis

In this final experiment, we repeat a similar study to that presented in Mitchell and Müller (2016). Our setup closely resembles the gravity wave test of Skamarock and Klemp (1994) extended to a spherical annulus. We initialize the velocity in a simple solid-body rotation and introduce a localized buoyancy anomaly. A Coriolis term is introduced as a function of the Cartesian coordinate z and is constant along lines of latitude (f plane): $\mathrm{2}\mathbf{\Omega }=\mathrm{2}{\mathrm{\Omega }}_r\frac{z}{R}\widehat{\mathbf{z}}$, with angular rotation rate ${\mathrm{\Omega }}_r=\mathrm{7.292}\times {\mathrm{10}}^{-\mathrm{5}}{\text{s}}^{-\mathrm{1}}$. We fix the resolution of the problem and run the solver over a range of Δt. We measure this by adjusting the horizontal acoustic Courant number ${\mathit{\lambda }}_C=c\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}$, where c is the speed of sound and Δx is the horizontal resolution.

Note that the range of Courant numbers used in this paper exceeds what is typical in operational forecast settings (typically between 𝒪(2)–𝒪(10)). The grid setup mirrors that of actual global weather models; we extrude a spherical mesh of the Earth upwards to a height of 85 km. The setup for the different discretizations (including degrees of freedom for the velocity–pressure and hybridized systems) is presented in Table 6.

Table 6 Grid setup and discretizations for the acoustic Courant number study. The total number of degrees of freedom (dofs) for the mixed (velocity and pressure) and hybridizable (velocity, pressure, and trace) discretizations are shown in the last two columns (millions). The vertical resolution is fixed across all discretizations.

It was shown in Mitchell and Müller (2016) that using a sparse approximation of the pressure Schur complement of the form

$$\begin{array}{}\text{(112)}& \stackrel{\mathrm{̃}}{\mathbf{H}}={\mathbf{M}}_p+c^{\mathrm{2}}{\displaystyle \frac{\mathrm{\Delta }{t}^{\mathrm{2}}}{\mathrm{4}}}\mathbf{D}{\left(\text{Diag}\left({\stackrel{\mathrm{̃}}{\mathbf{A}}}_{\mathbf{u}}\right)\right)}^{-\mathrm{1}}{\mathbf{D}}^T\end{array}$$

served as a good preconditioner, leading to a system that was amenable to multigrid methods and resulted in a Courant-number-independent solver. However, when the Coriolis term is included, this is no longer the case: the diagonal approximation $\text{Diag}\left({\stackrel{\mathrm{̃}}{\mathbf{A}}}_{\mathbf{u}}\right)$ becomes worse with increasing λ_C. To demonstrate this, we solve the gravity wave system on a low-resolution grid (10 km lid, 10 vertical levels, maintaining the same cell aspect ratio as in Table 6) using the Schur complement factorization Eq. (106). LU factorizations are applied to invert both ${\stackrel{\mathrm{̃}}{\mathbf{A}}}_{\mathbf{u}}^{-\mathrm{1}}$ and ${\stackrel{\mathrm{̃}}{\mathbf{H}}}^{-\mathrm{1}}$. Inverting the Schur complement H⁻¹ is done using preconditioned GMRES iterations, and a flexible-GMRES algorithm is used on the full velocity–pressure system. If ${\stackrel{\mathrm{̃}}{\mathbf{H}}}^{-\mathrm{1}}$ is a good approximation to H⁻¹, then we should see low iteration counts in the Schur complement solve. Figure 8 shows the results of this study for a range of Courant numbers.

Figure 8 Number of Krylov iterations to invert the Helmholtz system using ${\stackrel{\mathrm{̃}}{H}}^{-\mathrm{1}}$ as a preconditioner. The preconditioner is applied using a direct LU factorization within a GMRES method on the entire pressure equation. While the lowest-order methods grow slowly over the Courant number range, the higher-order (by only 1 approximation order) methods quickly degrade and diverge after the critical range λ_C=𝒪(2)–𝒪(10). At λ_C>32, the solvers take over 150 iterations.

Figure 9 Courant number parameter test run on a fully loaded compute node. Both figures display the hybridized solver for each discretization, described in Table 6. Panel (a) displays to total iteration count (preconditioned GCR) to solve the trace system to a relative tolerance of 10⁻⁵. Panel (b) displays the relative work of each solver, which takes into account the time required to forward eliminate and locally recover the velocity and pressure.

For the lower-order methods, the number of iterations to invert H grows slowly but remains under control. Increasing the approximation degree by 1 results in degraded performance. As Δt increases, the number of Krylov iterations needed to invert the system to a relative tolerance of 10⁻⁵ grows rapidly. It is clear that this sparse approximation is not robust against Courant number. This can be explained by the fact that diagonalizing the velocity operator fails to take into account the effects of the Coriolis term (which appear in off-diagonal positions in the operator). Even if one were to use traditional mass lumping (row-wise sums), the Coriolis effects are effectively canceled out due to asymmetry.

Hybridization avoids this problem entirely: we always construct an exact Schur complement and only have to worry about solving the trace system Eq. (111). We now show that this approach (described in Sect. 5.3.2) is much more robust to changes in Δt. We use the same workstation as for the three-dimensional CG/HDG problem in Sect. 5.1 (executed with a total of 40 MPI processes). Figure 9 shows the parameter test for all the discretizations described in Table 6. We see that, in terms of total number of GCR iterations needed to invert the trace system, hybridization is far more controlled as Courant number increases. They largely remain constant throughout the entire parameter range, only varying by an iteration or two. It is not until after λ_C>32 that we begin to see a larger jump in the number of GCR iterations. This is expected, since the Coriolis operator causes the problem to become singular for very large Courant numbers. However, unlike with the approximate Schur complement solver, iteration counts are still under control. In particular, each method (lowest and higher order) remains constant throughout the critical range (shaded in gray in Fig. 9).

In Fig. 9b, we display the ratio of execution time and the time-to-solution at the lowest Courant number of 2. We perform this normalization to better compare the lower and higher-order methods (and discretizations on triangular prisms vs. extruded quadrilaterals). The calculation of the ratios includes the time needed to eliminate and reconstruct the hybridized velocity–pressure variables. The fact that the hybridization solver remains close to 1 demonstrates that the entire solution procedure is largely λ_C-independent until around λ_C=32. The overall trend is largely the same as the number of Krylov iterations to reach solver convergence. This is due to our hybridization approach being solver dominated, with local operations like forward elimination together with local recovery taking approximately one-third of the total execution time for each method. The percentage breakdown of the hybridization solver is similar to what is already presented in Sect. 5.1.2.

Implicitly treating the Coriolis term has been discussed for semi-implicit discretizations of large-scale geophysical flows (Temperton, 1997; Cullen, 2001; Nechaev and Yaremchuk, 2004). The addition of the Coriolis term presents a particular challenge to the solution finite element discretizations of these equations since it increases the difficulty of finding a good sparse approximation of the nonsymmetric elliptic operator. Hybridization shows promise here, as it allows for the assembly of the elliptic equation that both captures the effects of rotation and results in a sparse linear system.

6 Conclusions

We have presented Slate, and shown how this language can be used to create concise mathematical representations of localized linear algebra on the tensors corresponding to finite element forms. We have shown how this DSL can be used in tandem with UFL in Firedrake to implement solution approaches making use of automated code generation for static condensation, hybridization, and localized post-processing. Setup and configuration are done at runtime, allowing one to switch in different discretizations at will. In particular, this framework alleviates much of the difficulty in implementing such methods within intricate numerical code and paves the way for future low-level optimizations. In this way, the framework in this paper can be used to help enable the rapid development and exploration of new hybridization and static condensation techniques for a wide class of problems. We remark here that the reduction of global matrices via element-wise algebraic static condensation, as described in Guyan (1965) and Irons (1965) is also possible using Slate, including other more general static condensation procedures outside the context of hybridization.

Our approach to preconditioner design revolves around its composable nature, in that these Slate-based implementations can be seamlessly incorporated into complicated solution schemes. In particular, there is current research in the design of dynamical cores for numerical weather prediction using implementations of hybridization and static condensation with Slate (Bauer and Cotter, 2018; Shipton et al., 2018). The performance of such methods for geophysical flows are a subject of ongoing investigation.

In this paper, we have provided some examples of hybridization procedures for compatible finite element discretizations of geophysical flows. These approaches avoid the difficulty in constructing sparse approximations of dense elliptic operators. Static condensation arising from hybridizable formulations can best be interpreted as producing an exact Schur complement factorization on the global hybridizable system. This eliminates the need for outer iterations from a suitable Krylov method to solve the full mixed system and replaces the original global mixed equations with a condensed elliptic system. More extensive performance benchmarks, which require detailed analysis of the resulting operator systems arising from hybridization, are a necessary next step to determine whether hybridization provides a scalable solution strategy for compatible finite elements in operational settings.

Appendix A: Semi-implicit method for the shallow water system

For some tessellation, 𝒯_h, our semi-discrete mixed method for Eqs. (81)–(82) seeks approximations $({\mathbf{u}}_h,D_h)\in U_h\times V_h\subset H\left(\mathrm{div}\right)\times L^{\mathrm{2}}$ satisfying

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}& {\left(\mathbf{w},{\displaystyle \frac{\partial {\mathbf{u}}_h}{\partial t}}\right)}_{{\mathcal{T}}_h}-{\left({\mathrm{\nabla }}^{⟂}\left(\mathbf{w}\cdot {\mathbf{u}}_h^{⟂}\right),{\mathbf{u}}_h^{⟂}\right)}_{{\mathcal{T}}_h}+{\left(\mathbf{w},f{\mathbf{u}}_h^{⟂}\right)}_{{\mathcal{T}}_h}\\ & +{〈\left[\left[{\mathbf{n}}^{⟂}\mathbf{w}\cdot {\mathbf{u}}_h^{⟂}\right]\right],{\stackrel{\mathrm{̃}}{\mathbf{u}}}_h^{⟂}〉}_{\partial {\mathcal{T}}_h}\end{array}\end{array}$$

$$-{\left(\mathrm{\nabla }\cdot \mathbf{w},g\left(D_h+b\right)+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}|{\mathbf{u}}_h{|}^{\mathrm{2}}\right)}_{{\mathcal{T}}_h}=\mathrm{0},\forall \mathbf{w}\in U_h,$$

$$\begin{array}{}\text{(A2)}& \begin{array}{rl}& {\left(\mathit{\varphi },{\displaystyle \frac{\partial D_h}{\partial t}}\right)}_{{\mathcal{T}}_h}-{\left(\mathrm{\nabla }\mathit{\varphi },{\mathbf{u}}_h{D}_h\right)}_{{\mathcal{T}}_h}+{〈\left[\left[\mathit{\varphi }{\mathbf{u}}_h\right]\right],{\stackrel{\mathrm{̃}}{D}}_h〉}_{\partial {\mathcal{T}}_h},\\ & =\mathrm{0},\forall \mathit{\varphi }\in V_h,\end{array}\end{array}$$

where $\stackrel{\mathrm{̃}}{\cdot }$ indicates that the value of the function should be taken from the upwind side of each facet. The discretization of the velocity advection operator is an extension of the energy-conserving scheme of Natale and Cotter (2017) to the shallow water equations.

The time-stepping scheme follows a Picard iteration semi-implicit approach, where predictive values of the relevant fields are determined via an explicit step of the advection equations and corrective updates are generated by solving an implicit linear system (linearized about a state of rest) for $(\mathrm{\Delta }{\mathbf{u}}_h,\mathrm{\Delta }{D}_h)\in {\mathbf{U}}_h\times V_h$, given by

$$\begin{array}{}\text{(A3)}& \begin{array}{rl}& {\left(\mathbf{w},\mathrm{\Delta }{\mathbf{u}}_h\right)}_{{\mathcal{T}}_h}+{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{2}}}{\left(\mathbf{w},f\mathrm{\Delta }{\mathbf{u}}_h^{⟂}\right)}_{{\mathcal{T}}_h}-{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{2}}}{\left(\mathrm{\nabla }\cdot \mathbf{w},g\mathrm{\Delta }{D}_h\right)}_{{\mathcal{T}}_h}\\ & =-R_{\mathbf{u}}[{\mathbf{u}}_h^{n+\mathrm{1}},D_h^{n+\mathrm{1}};\mathbf{w}],\forall \mathbf{w}\in U_h,\end{array}\end{array}$$

$$\begin{array}{}\text{(A4)}& \begin{array}{rl}& {\left(\mathit{\varphi },\mathrm{\Delta }{D}_h\right)}_{{\mathcal{T}}_h}+{\displaystyle \frac{H\mathrm{\Delta }t}{\mathrm{2}}}{\left(\mathit{\varphi },\mathrm{\nabla }\cdot \mathrm{\Delta }{\mathbf{u}}_h\right)}_{{\mathcal{T}}_h}\\ & =-R_D[{\mathbf{u}}_h^{n+\mathrm{1}},D_h^{n+\mathrm{1}};\mathit{\varphi }],\forall \mathit{\varphi }\in V_h,\end{array}\end{array}$$

where H is the mean layer depth and R_u and R_D are residual linear forms that vanish when ${\mathbf{u}}_h^{n+\mathrm{1}}$ and $D_h^{n+\mathrm{1}}$ are solutions to the implicit midpoint rule time discretization of Eqs. (A1)–(A2). The residuals are evaluated using the predictive values of ${\mathbf{u}}_h^{n+\mathrm{1}}$ and $D_h^{n+\mathrm{1}}$.

The implicit midpoint rule time discretization of the nonlinear rotating shallow water Eqs. (A1)–(A2) is

$$\begin{array}{}\text{(A5)}& \begin{array}{rl}& {\left(\mathbf{w},{\mathbf{u}}_h^{n+\mathrm{1}}-{\mathbf{u}}_h^n\right)}_{{\mathcal{T}}_h}-\mathrm{\Delta }t{\left({\mathrm{\nabla }}^{⟂}\left(\mathbf{w}\cdot {{\mathbf{u}}_h^{*}}^{⟂}\right),{{\mathbf{u}}_h^{*}}^{⟂}\right)}_{{\mathcal{T}}_h}+\\ & \mathrm{\Delta }t{\left(\mathbf{w},f{{\mathbf{u}}_h^{*}}^{⟂}\right)}_{{\mathcal{T}}_h}+\mathrm{\Delta }t{〈\left[\left[{\mathbf{n}}^{⟂}\mathbf{w}\cdot {{\mathbf{u}}_h^{*}}^{⟂}\right]\right],{{\stackrel{\mathrm{̃}}{\mathbf{u}}}_h^{*}}^{⟂}〉}_{\partial {\mathcal{T}}_h}\\ & -\mathrm{\Delta }t{\left(\mathrm{\nabla }\cdot \mathbf{w},g\left(D_h^{*}+b\right)+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}|{\mathbf{u}}_h^{*}{|}^{\mathrm{2}}\right)}_{{\mathcal{T}}_h}=\mathrm{0},\forall \mathbf{w}\in U_h,\end{array}\end{array}$$

$$\begin{array}{}\text{(A6)}& \begin{array}{rl}& {\left(\mathit{\varphi },D_h^{n+\mathrm{1}}-D_h^n\right)}_{{\mathcal{T}}_h}-\mathrm{\Delta }t{\left(\mathrm{\nabla }\mathit{\varphi },{\mathbf{u}}_h^{*}{D}_h^{*}\right)}_{{\mathcal{T}}_h}\\ & +\mathrm{\Delta }t{〈\left[\left[\mathit{\varphi }{\mathbf{u}}_h^{*}\right]\right],{\stackrel{\mathrm{̃}}{D}}_h^{*}〉}_{\partial {\mathcal{T}}_h},=\mathrm{0},\forall \mathit{\varphi }\in V_h,\end{array}\end{array}$$

where ${\mathbf{u}}_h^{*}=({\mathbf{u}}_h^{n+\mathrm{1}}+{\mathbf{u}}_h^n)/\mathrm{2}$ and $D_h^{*}=(D_h^{n+\mathrm{1}}+D_h^n)/\mathrm{2}$.

One approach to construct the residual functionals R_u and R_D would be to simply define these from Eqs. (A5)–(A6). However, this leads to a small critical time step for the stability of the scheme. To make the numerical scheme more stable, we define residuals as follows. For R_u, we first solve for v_h∈U_h such that

$$\begin{array}{}\text{(A7)}& \begin{array}{rl}& {\left(\mathbf{w},{\mathbf{v}}_h-{\mathbf{u}}_h^n\right)}_{{\mathcal{T}}_h}-\mathrm{\Delta }t{\left({\mathrm{\nabla }}^{⟂}\left(\mathbf{w}\cdot {{\mathbf{u}}_h^{*}}^{⟂}\right),{{\mathbf{v}}_h^{\mathrm{♯}}}^{⟂}\right)}_{{\mathcal{T}}_h}\\ & +\mathrm{\Delta }t{\left(\mathbf{w},f{{\mathbf{v}}_h^{\mathrm{♯}}}^{⟂}\right)}_{{\mathcal{T}}_h}+\mathrm{\Delta }t{〈\left[\left[{\mathbf{n}}^{⟂}\mathbf{w}\cdot {{\mathbf{u}}_h^{*}}^{⟂}\right]\right],{{\stackrel{\mathrm{̃}}{\mathbf{v}}}_h^{\mathrm{♯}}}^{⟂}〉}_{\partial {\mathcal{T}}_h}\\ & -\mathrm{\Delta }t{\left(\mathrm{\nabla }\cdot \mathbf{w},g\left(D_h^{*}+b\right)+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}|{\mathbf{u}}_h^{*}{|}^{\mathrm{2}}\right)}_{{\mathcal{T}}_h}=\mathrm{0},\forall \mathbf{w}\in U_h,\end{array}\end{array}$$

where ${\mathbf{v}}_h^{\mathrm{♯}}=({\mathbf{v}}_h+{\mathbf{u}}_h^n)/\mathrm{2}$. This is a linear variational problem. Then,

$$\begin{array}{}\text{(A8)}& R_{\mathbf{u}}[{\mathbf{u}}_h^{n+\mathrm{1}},D_h^{n+\mathrm{1}};\mathbf{w}]={\left(\mathbf{w},{\mathbf{v}}_h-{\mathbf{u}}_h^{n+\mathrm{1}}\right)}_{{\mathcal{T}}_h}.\end{array}$$

Similarly, for R_D we first solve for E_h∈V_h such that

$$\begin{array}{}\text{(A9)}& \begin{array}{rl}& {\left(\mathit{\varphi },E_h-D_h^n\right)}_{{\mathcal{T}}_h}-\mathrm{\Delta }t{\left(\mathrm{\nabla }\mathit{\varphi },{\mathbf{u}}_h^{*}{E}_h^{\mathrm{♯}}\right)}_{{\mathcal{T}}_h}\\ & +\mathrm{\Delta }t{〈\left[\left[\mathit{\varphi }{\mathbf{u}}_h^{*}\right]\right],{\stackrel{\mathrm{̃}}{E}}_h^{\mathrm{♯}}〉}_{\partial {\mathcal{T}}_h}=\mathrm{0},\forall \mathit{\varphi }\in V_h,\end{array}\end{array}$$

where $E_h^{\mathrm{♯}}=(E_h+D_h^n)/\mathrm{2}$. This is also a linear problem. Then,

$$\begin{array}{}\text{(A10)}& R_D[{\mathbf{u}}_h^{n+\mathrm{1}},D_h^{n+\mathrm{1}};\mathit{\varphi }]={\left(\mathit{\varphi },E_h-D_h^{n+\mathrm{1}}\right)}_{{\mathcal{T}}_h}.\end{array}$$

This process can be thought of as iteratively solving for the average velocity and depth that satisfies the implicit midpoint rule discretization. Both Eqs. (A7) and (A9) can be solved separately, since there is no coupling between them. The fields v_h and E_h are then used to construct the right-hand side for the implicit linearized system in Eqs. (A3)–(A4). Once the system is solved, the solution (Δu_h,ΔD_h) is then used to update the iterative values of ${\mathbf{u}}_h^{n+\mathrm{1}}$ and $D_h^{n+\mathrm{1}}$ according to $({\mathbf{u}}_h^{n+\mathrm{1}},D_h^{n+\mathrm{1}})\leftarrow ({\mathbf{u}}_h^{n+\mathrm{1}}+\mathrm{\Delta }{\mathbf{u}}_h,D_h^{n+\mathrm{1}}+\mathrm{\Delta }{D}_h)$, having initially chosen $({\mathbf{u}}_h^{n+\mathrm{1}},D_h^{n+\mathrm{1}})=({\mathbf{u}}_h^n,D_h^n)$.