Bringing Trimmed Serendipity Methods to Computational Practice in Firedrake

Over the last 15 years, conforming finite element methods for Hodge-Laplace type problems on simplicial and cubical meshes have been analyzed and categorized using the mathematical theory of Finite Element Exterior Calculus (FEEC) (Arnold et al., 2006, 2010, 2015). While this effort was initially focused on placing related theoretical results under a single, unified mathematical framework, it also exposed the frontier of knowledge and spawned a wealth of spin-off projects that improved awareness, understanding, and implementation of the methods described. Among these follow-on projects was a revisiting of the notion of “serendipity” finite elements, a particular type of finite element method dating back to the 1970s. Named for its seemingly too-good-to-be-true computational benefit, a serendipity finite element method converged to the correct solution of a partial differential equation (PDE) at an equivalent rate but with fewer degrees of freedom than the corresponding tensor product method. The simplest and most well-known serendipity method replaces the quadratic 9-node square element with a “serendipity” 8-node element that has no interior degree of freedom, but still converges at a quadratic rate in the appropriate sense.

The core idea of serendipity elements, i.e. reducing degrees of freedom from tensor product methods without reducing the accuracy of the method, floated around computational engineering communities until it found a resurgence under the FEEC framework. In a seminal work by Arnold and Awanou (Arnold and Awanou, 2011) the scalar-valued serendipity elements were formalized as providing approximation on a space of polynomials $\mathcal{S}_{r}$ that nests between maximum-degree $r$ polynomials $\mathcal{P}_{r}$ and tensor product degree $r$ polynomials $\mathcal{Q}_{r}$. In a subsequent work (Arnold and Awanou, 2014), this notion was extended using FEEC to introduce “serendipity vector elements” to the literature – analogues of the famous Nédélec elements (Nédélec, 1980, 1986) – which led to the widely circulated Periodic Table of Finite Elements (Arnold and Logg, 2014).

Despite excitement at these developments, implementing serendipity elements and realizing any potential computational benefits proved to be a significant challenge. While creating data structures for “all polynomials up to degree $r$” is straightforward, simply trying to write down the approximation spaces used by serendipity elements – especially the vector-valued elements – requires a substantial amount of mathematical notation and explanation. As a consequence, the serendipity elements have never been implemented or rigorously studied beyond a few special use cases.

Recent work on a closely related family of methods – called trimmed serendipity elements (Gillette and Kloefkorn, 2019) – has opened the door to potential widespread usage and benefit from the notions of serendipity theory just described. The trimmed serendipity spaces are identical to the serendipity spaces in the scalar-valued cases but are distinct in the vector-valued cases where, notably, they have fewer degrees of freedom than the serendipity elements of the same order. Since the computational motivation for serendipity methods is entirely about reducing degree of freedom count while preserving approximation power, the smaller dimensionality of the trimmed serendipity spaces obviates the need to implement the vector-valued “non-trimmed” serendipity elements.

Implementation of trimmed serendipity spaces is feasible by merging two independent efforts: the systematic definition of computational basis functions for these methods from work by Gillette, Kloefkorn and Sanders (Gillette et al., 2019) and the easily extensible open-source finite element software package Firedrake  (Rathgeber et al., 2016). The Unified Form Language (Logg et al., 2012; Alnæs et al., 2014), inspired by FEEC, provides a common backbone for translating the basis functions from their formal statement into a code structure using a well-established, high level interface. We explain the idea of the implementation by walking through the process of discretizing a PDE and selecting a finite element method for approximating its solution.

Among its many advantages, Finite Element Exterior Calculus (FEEC) (Arnold et al., 2006, 2010, 2015) gives an easy, unified way to notate different element types. The four best known families of elements are denoted $\mathcal{P}^{-}_{r}\Lambda^{k}$, $\mathcal{P}_{r}\Lambda^{k}$, $\mathcal{Q}^{-}_{r}\Lambda^{k}$, and $\mathcal{S}_{r}\Lambda^{k}$, which are, respectively, the trimmed polynomial, polynomial, tensor product, and serendipity elements of order $r$ using $k$-forms. The $\mathcal{P}$ and $\mathcal{P}^{-}$ spaces are defined over meshes of simplices (triangles, tetrahedra, etc) while the $\mathcal{Q}^{-}$ and $\mathcal{S}$ spaces are defined over meshes of hypercubes (squares, cubes, etc). The optional notation $(\square_{n})$ specifies that the space is constructed over the $n$-dimensional cube in $\mathbb{R}^{n}$, but we will frequently omit this addition if $n$ is clear from context.

The mathematical results from FEEC regarding these four families synthesize decades of research into what constitutes a stable finite element, i.e. a numerical method that can be proven to converge to the correct solution of certain PDEs in certain norms at a prescribed rate, indicated by the subscript $r$. In any dimension, the $0$-form spaces provide scalar-valued, $H^{1}$-conforming elements. In 2D, $1$-forms can represent both ${H}(\operatorname{curl})$ and ${H}(\operatorname{div})$ elements, depending on the orientation of the DOFs defined on the edges of a mesh. In 3D, $1$-forms correspond to ${H}(\operatorname{curl})$ elements while $2$-forms correspond to ${H}(\operatorname{div})$ elements. Notably, the serendipity family $\mathcal{S}_{r}\Lambda^{k}$ is the youngest, least implemented, and hence least understood among all these families. In particular, the 1-form and 2-form (regular) serendipity elements in 3D were only characterized in 2014 by Arnold and Awanou (2014), whereas the equivalent tensor product elements were first described more than 30 years earlier in Nédélec (1980, 1986).

Firedrake uses FIAT (Kirby, 2004, 2012) to provide finite element basis functions on reference elements. To implement a new element in FIAT, we must provide both rules to tabulate basis functions and their derivatives at reference element points and a data structure that assigns basis functions to particular reference element facets. Said element is then made available in Firedrake by providing a symbolic name (in UFL (Alnæs et al., 2014)) and a translation from symbolic name to concrete implementation in the form compiler TSFC (Homolya et al., 2018). While FIAT initially considered a very wide range of finite elements (Kirby et al., 2012), it would seek to express their bases as linear combinations of orthogonal polynomials. However, for some elements, it is easier to directly describe the basis functions. We follow the construction of Gillette et al. (2019) which provides explicit formulae for the basis functions for each of the trimmed serendipity elements and directly implement tabulation of the basis functions. To provide tabulations of derivatives, we implement the basis functions symbolically with SymPy (Meurer et al., 2017) and compute derivatives symbolically.

We use the decompositions from Gillette et al. (2019) to group the basis functions according to the geometric portions of a reference mesh element (vertices, edges, faces, or cell interiors). For instance, a basis for $\mathcal{S}^{-}_{r}\Lambda^{1}(\square_{3})$ – the trimmed serendipity ${H}(\operatorname{curl})$-conforming element in 3D – can be decomposed as

Subsets in these decompositions denoted with an $E$, $F$, or $I$ are defined on edges, faces, and interior, respectively, of the cubical cell in 3D. The subsets $\tilde{F}$ and $\tilde{I}$ are extra sets of basis functions defined on the faces or the interior that follow a different pattern in their definitions than those of the functions in $E$, $F$, or $I$.

To see how this plays out in the Firedrake implementation, consider the ${H}(\operatorname{curl})$ elements for trimmed serendipity at order $r=2$ in $n=3$ dimensions, indicating the space $S_{2}^{-}\Lambda^{1}$ in 3D. The space $\tilde{I}_{r}\Lambda^{1}$ is defined to be empty for $r<4$, so there are only two sets of functions to include in this case: one set associated to edges of the reference element (the $E_{i}\Lambda^{1}$ sum) and one set for the faces of the reference element (the $\tilde{F}_{2}\Lambda^{1}$ functions). According to Eq. 2, the basis functions can be decomposed as

We then implement these in Firedrake as follows. The first step is to determine the number of DOFs we will need on the reference element where we will define the basis functions. To do this, we count the DOFs for each mesh entity on the reference element (vertex, edge, face, and interior). For example, $\mathcal{S}_{2}^{-}\Lambda^{1}(\square_{3})$ should have no DOFs at the vertices (these only come into play for the $H^{1}$-conforming elements), two DOFs on each edge, two DOFs on each face, and no DOFs on the interior. This agrees with Eq. 3, where $E_{i}\Lambda^{1}$ supplies one DOF to each edge for each $i=0,1$, and then $\tilde{F}_{2}\Lambda^{1}$ supplies two DOFs to each face. An illustration of this can be seen in the second column, second row of Table 1.

With the number of DOFs assigned to each mesh entity in the reference element, we can then define the basis functions. The order of definition is important so that basis functions are matched to the proper mesh entity. Note that Eq. 3 doesn’t explicitly give a way to order the basis functions. Instead, we need to use the properties of the basis functions to determine the correct ordering. This is best illustrated by an example. Two of the “edge” basis functions that are contained in the sum of the $E_{i}\Lambda^{1}$ sets are $(y+1)(z+1)dx$ and $x(y+1)(z+1)dx$. Notice that these functions have no $dy$ or $dz$ portions. Therefore, these function vanish on any edge not parallel to the $x$ axis. Further, the polynomial coefficients of these forms indicate that they also vanish on the planes $y=-1$ and $z=-1$. Thus, the only edge of the cube on which these functions do not vanish is the edge contained in the line $\{y=1\}\cap\{z=1\}$. This edge is shown in blue and labeled with an “e” in Fig. 3.

This process is what determines the ordering for the basis functions. If, in the first step described above where we are assigning DOFs to mesh entities, we assign the first two DOFs to be on the edge of the reference element where $y=1$ and $z=1$, then we must define the basis functions $(y+1)(z+1)dx$ and $x(y+1)(z+1)dx$ as our first two basis functions.

While the formulas above are taken from Gillette et al. (2019), these monomials have poor conditioning at higher orders. Therefore, we use Legendre polynomials obtained symbolically from SymPy via the legendre function denoted leg. Then we write the differential forms in vector notation. For our examples of $(y+1)(z+1)dx$ and $x(y+1)(z+1)dx$, we get tuple([(leg(j,x_mid)*dz[1]*dy[1],0,0)]), where leg(j, x_mid) is computed at the midpoint between two vertices and used for the $1$ or $x$ coefficient, and the dy[1] and dz[1] are used for the values of $(y+1)$ and $(z+1)$, respectively. After repeating this process for each of the edges and associated basis functions, we then also do a similar process for the face functions in the set $\tilde{F}_{2}\Lambda^{1}$.

This process changes only slightly if we instead had considered the the $2$-forms $\mathcal{S}^{-}_{2}\Lambda^{2}(\square_{3})$. In this scenario, we would have the basis functions given by the equation

The newly supported elements, mapping FEEC spaces onto names in UFL are shown in the lower half of Table 2. Modifying the code from LABEL:lst:mixedPoissonCode to utilize trimmed serendipity spaces rather than tensor product spaces is then simply a case of replacing the element names in the FunctionSpace definitions with appropriate trimmed space names. Concretely, the new FunctionSpace definitions are shown in LABEL:lst:pde_using_trimmed_serendipity, the rest of the code remains unchanged.

The following experiments show the benefits and costs of using trimmed serendipity elements in comparison to traditional tensor product elements. We first present a basic projection example as a means of confirming approximation properties of our elements. Next, we present results on a primal Poisson problem (to test $H^{1}$ elements), a mixed Poisson problem (to test ${H}(\operatorname{div})$ and $L^{2}$ elements), and a cavity resonator problem (to test ${H}(\operatorname{curl})$ elements). Since the ${H}(\operatorname{curl})$ and ${H}(\operatorname{div})$ elements are only a rotation of the DOFs on the edges of elements in 2D, we do not give an experiment using ${H}(\operatorname{curl})$ elements in 2D.

The experiments were performed either on a single cluster compute node with 2x AMD EPYC 7642 48-core (Rome) processors (2.4GHz) and 512GB of memory running CentOS 7 or a similar node with 3TB of memory. The 2D experiments were all done using the 512GB node, while in 3D, the 4th order experiments were run on the 3TB node. Each job was run by submitting a SLURM script that requested one node in isolation to ensure no other jobs were running at the same time. We utilized on-node parallelism by requesting a full node and executing the jobs with mpirun set to use 24 processes, with a few exceptions for smaller cases that will be pointed out later. Timing data was collected after first performing a dry run of the code to warm the cache and then taking the minimum time over three subsequent runs. The timing results presented here depend upon the solver choice, and changing that may give different results illustrating the relative efficiency between $\mathcal{Q}^{-}$ and $\mathcal{S}^{-}$.

For simplicity, our numerical experiments all use a sparse direct solver. We expect the multigrid theory in Arnold et al. (2000, 2006) to carry over from existing $\mathcal{Q}^{-}$ spaces to $\mathcal{S}^{-}$ spaces. Optimal smoothers require aggregating degrees of freedom for vertex patches, and we anticipate that the reduction in local dimensionality that trimmed serendipity spaces offer will be beneficial in these contexts as well.