On the choice of finite element for applications in geodynamics – Part 2: A comparison of simplex and hypercube elements

Many geodynamical models are formulated in terms of the Stokes equations that are then coupled to other equations. For the numerical solution of the Stokes equations, geodynamics codes over the past decades have used essentially every finite element that has ever been proposed for the solution of this equation, on both triangular/tetrahedral (“simplex”) and quadrilaterals/hexahedral (“hypercube”) meshes. However, in many and perhaps most cases, the specific choice of element does not seem to have been the result of careful benchmarking efforts but based on implementation efficiency or the implementers' background.

In a first part of this paper (Thieulot and Bangerth, 2022), we have provided a comprehensive comparison of the accuracy and efficiency of the most widely used hypercube elements for the Stokes equations. We have done so using a number of benchmarks that illustrate “typical” geodynamic situations, specifically taking into account spatially variable viscosities. Our findings there showed that only Taylor–Hood-type elements with either continuous (Q₂×Q₁) or discontinuous ($Q_{\mathrm{2}}\times P_{-\mathrm{1}}$) pressure are able to adequately and efficiently approximate the solution of the Stokes equations.

In this, the second part of this work, we extend the comparison to simplex meshes. In particular, we compare triangular Taylor–Hood elements against the MINI element and one often referred to as the “Crouzeix–Raviart” element. We compare these choices against the accuracy obtained on hypercube Taylor–Hood elements with approximately the same computational cost. Our results show that, like on hypercubes, the Taylor–Hood element is substantially more accurate and efficient than the other choices. Our results also indicate that hypercube meshes yield slightly more accurate results than simplex meshes but that the difference is relatively small and likely unimportant given that hypercube meshes often lead to slightly denser (and consequently more expensive) matrices.

1 Introduction

Over the past decades, a large number of geodynamics simulation codes have been built on the finite element method. In the specific context of mantle convection and long-term dynamics simulators, the key component of many models that needs to be solved is the Stokes equations for which finite element methods are well suited but which leaves many choices still to be made: (1) should the element choice be of Taylor–Hood type (where the polynomial degree used for the velocity is chosen one higher than that for the pressure) or a stabilized equal-order element combination, or any of the nonconforming elements? (2) Should the reference cell for the underlying mesh be simplices (triangles or tetrahedra) or hypercubes (quadrilaterals or hexahedra)?

In an earlier part of this work (see Thieulot and Bangerth, 2022), we have extensively compared different hypercube choices for the finite element combination used to discretize the Stokes equations in the context of models that are relevant to geodynamic applications. Our conclusions there were that the lowest-order Taylor–Hood-type element (denoted by Q₂×Q₁ on quadrilaterals or $Q_{\mathrm{2}}\times P_{-\mathrm{1}}$ if one uses a discontinuous pressure element that then leads to local mass conservation) is the only one that produces accurate results in all circumstances. This conclusion is notwithstanding the fact that these elements are not cheap, owing to their higher-order shape functions and the consequent large number of entries in the system matrix. Yet, all other choices we have compared there – specifically, the stabilized Q₁×Q₁ and the unstable Q₁×P₀ combinations – were too inaccurate, were too unstable, and had too much difficulty representing the hydrostatic pressure component to be competitive. Of course these other choices are widely used in many existing codes, making studies such as Thieulot and Bangerth (2022) useful to inform what the next generation of codes should build on.

At the same time, in Thieulot and Bangerth (2022) we did not investigate whether simplex or hypercube meshes are better suited to the task. Historically, geodynamics has largely settled on the use of quadrilateral or hexahedral (“hypercube”) elements – somewhat separate from the rest of the finite element world that has traditionally predominantly used triangular or tetrahedral (“simplex”) meshes. The reasons for this deviation are likely rooted in the fact that the geometries of the domains used in geodynamics are largely rather simple: rectangles and boxes, along with circles, spheres, and shells. These geometries present no difficulties in meshing with hypercube cells, whereas the complex geometries frequently used in solid and fluid mechanics can often only reasonably be meshed using mesh generators that create simplex meshes. Still, one could of course also use simplex meshes in geodynamics, and in fact many codes have done so over the past decades; see, for example, Barr and Houseman (1996) (BASIL code, P₂×P₁), Dabrowski et al. (2008) (MILAMIN code, $P_{\mathrm{2}}^{+}\times P_{-\mathrm{1}}$), Tommasi et al. (2009) (FORGE2005 software, $P_{\mathrm{1}}^{+}\times P_{\mathrm{1}}$), Davies et al. (2011) (Fluidity code, P₂×P₁), Chertova et al. (2014) (SEPRAN, P₂×P₁), Paczkowski et al. (2014) (COMSOL, P₂×P₁), de Montserrat et al. (2019) (LaCoDe, $P_{\mathrm{2}}^{+}\times P_{-\mathrm{1}}$), Jones et al. (2021) (FEniCS project, P₂×P₁), and Ilangovan et al. (2024) (HyTeG framework, P₂×P₁). It is, therefore, a reasonable question of whether it would result in more accurate simulations for the same computational cost or less costly simulations at the same accuracy.

We are not aware of systematic comparisons between the two choices of reference cell – simplex or hypercube – in the geodynamics literature. Perhaps surprisingly, there is also not a large body of literature on the topic in other disciplines, nor is there a strong oral “lore” in the scientific computing community about which of the two approaches is better. In our search for past work, we have found a modestly informative recent publication that clearly illustrates the benefits of quadratic over linear elements but only a weak preference for triangles/tetrahedra over quadrilaterals/hexahedra (Schneider et al., 2022). That publication also contains references to other, earlier studies in the same direction; it is worth also pointing out Terrel et al. (2012) as another example of a study that compares different elements, though not in great depth and not for applications relevant in geodynamics. On the other hand, while discretization accuracy matters, so does solver speed. In this regard, modern solver techniques intended to better utilize the power of CPUs over the limitations of memory latency, specifically matrix-free approaches, heavily build on the fact that shape functions on hypercube cells can be written as a tensor product of 1D functions and so are naturally more suited to hypercube cells (Kronbichler and Kormann, 2019; Munch et al., 2023).

Regardless of which reference element is better suited, using simplex meshes also opens up a number of other possibilities. Specifically, the number of stable Stokes element combinations for triangles and tetrahedra is substantially larger than it is for quadrilaterals and hexahedra, owing to decades of research on “nonconforming” elements – i.e., finite element choices whose basis functions are not continuous but have a sufficient number of structural properties (such as being continuous at edge midpoints) that it is not necessary to add specific stabilization terms to the weak formulation of the Stokes equations.¹ Examples of such nonconforming elements include the Brezzi–Douglas–Marini (BDM) element (Brezzi et al., 1985) and the Crouzeix–Raviart nonconforming P₁ element² (Braess, 2007). Nonconforming elements have of course also been developed for hypercube meshes – see, e.g., the Rannacher–Turek element (Rannacher and Turek, 1992) or the DSSY element (Douglas et al., 1999). However, in contrast to their relatives defined on simplex meshes, they have not found widespread use and are less often implemented in widely used finite element libraries. Thus, nonconforming elements are generally only considered viable choices on simplex meshes.

In practice, however, nonconforming elements have never found much use in the geodynamics community. As a consequence, while we consider them viable alternatives (and potential targets for future studies), we will not include any nonconforming elements in this work (see also Sect. 3.3).

We end this overview by mentioning that while we have not found much literature that quantitatively compares reference cell and element choices, the book by Gresho and Sani (2000) contains an extensive and excellent overview of the many available choices for the Stokes equations in Sects. 3.13.2 to 3.13.6, covering more than 150 pages. Tables 3.13-1 to 3.13-4, along with lengthy comments throughout the section, provide arguments that lead the authors of that book to favor hypercube cells over simplex cells and to favor Taylor–Hood-type elements over others, based on qualitative arguments and references to the literature from the 1970s, 1980s, and 1990s – opinions that we share, based on the results of this paper and of Thieulot and Bangerth (2022). Yet, the authors also state that the question is not at all trivial, is not settled, and is in need of systematic quantitative comparisons. In any case, both the book and the literature cited therein exclusively consider the isoviscous Stokes equations which typically has a much smoother solution than the ones we find in geodynamic applications with their highly variable viscosity coefficient. As a consequence, we believe that the comparisons we provide here are useful not only because they are quantitative, but also because they are specific to the kinds of applications we typically encounter in our discipline. We should also mention the paper by Pelletier et al. (1989) that contains qualitative comparisons between Taylor–Hood elements with continuous and discontinuous pressures and that advocates for the use of discontinuous pressure elements. However, we consider its interpretation difficult in today's context because the meshes used there are quite coarse, and we think it likely that their recommendations are perhaps no longer as applicable to today's much finer meshes as they were in the historical context nearly 40 years ago when the paper was written.

Goals of this paper

Given the setting described above, our goal in this contribution is to compare finite element choices on simplex and hypercube meshes both qualitatively and quantitatively. For hypercube meshes, our previous work in Thieulot and Bangerth (2022) has already indicated that only the Taylor–Hood variants Q₂×Q₁ or $Q_{\mathrm{2}}\times P_{-\mathrm{1}}$ are reasonable choices, whereas equal-order elements are not. Based on this observation, we then pose the following two questions for the current work:

1. What is the best choice of finite element on simplex meshes?

2. How does the best choice of finite element on simplex meshes compare to the choice of the Q₂×Q₁ or $Q_{\mathrm{2}}\times P_{-\mathrm{1}}$ element on hypercube meshes?

For our numerical comparisons, we will consider both the accuracy and computational cost of a finite element as a function of the mesh size (or number of unknowns) as a criterion. The elements we will consider for simplex meshes include the P₂×P₁ Taylor–Hood-type element, the P₂×P₀ element (the cheapest stable element with discontinuous pressure), the MINI element $P_{\mathrm{1}}^{+}\times P_{\mathrm{1}}$, and the “Crouzeix–Raviart” element $P_{\mathrm{2}}^{+}\times P_{-\mathrm{1}}$.

Outline of the paper

In the remainder of this paper, we will first briefly state the equations we seek to solve (Sect. 2). In Sect. 3 we will then discuss the finite elements one can choose on simplex meshes for the discretization of the Stokes equations, along with a description of the elements we do and do not compare in this work. Section 4 then provides a numerical comparison of these elements using a series of benchmarks that illustrate how solutions of geodynamic models often behave. We conclude in Sect. 5.

2 The governing equations

As in the first part of this work, we will here be concerned with the accurate numerical solution of the incompressible Stokes equations:

$$\begin{array}{}\text{(1)}& {\displaystyle }-\mathrm{\nabla }\cdot \left[\mathrm{2}\mathit{\eta }\mathit{\epsilon }\right(\mathit{u}\left)\right]+\mathrm{\nabla }p& {\displaystyle }=\mathit{\rho }\mathit{g},\text{(2)}& {\displaystyle }-\mathrm{\nabla }\cdot \mathit{u}& {\displaystyle }=\mathrm{0},\end{array}$$

where η is the viscosity, ρ the density, and g the gravity vector, and we will denote by ε(⋅) the symmetric gradient operator defined by $\mathit{\epsilon }\left(\mathit{u}\right)=\frac{\mathrm{1}}{\mathrm{2}}(\mathrm{\nabla }\mathit{u}+\mathrm{\nabla }{\mathit{u}}^T)$. Ω⊂ℝ^d, d = 2 or 3 is the domain of interest. Both the viscosity η and the density ρ will, in typical applications, be spatially variable; the variability is often introduced through nonlinear dependencies on the strain rate ε(u) and/or the pressure p, but the exact reasons are not of relevance to us here: the important point is that these coefficients may vary strongly and on short length scales.

In actual applications, the equations above will be completed by appropriate boundary conditions and will be augmented by additional and often time-dependent equations, such as ones that describe the evolution of the temperature field or of the composition of rocks (see, for example, Schubert et al., 2001; Turcotte and Schubert, 2012). This coupling is also not of interest to us here, nor is the fact the “true” equations in geodynamics are often compressible – in most cases, the equations above will have to be solved as a “sub-problem” to what one really wants to do, and the efficiency of a discretization of these equations then translates to a lower bound for the efficiency of solving the outer problem.

3 Discretization

3.1 Elements and element combinations

The finite element discretization of the Stokes equations is complicated by the fact that one cannot choose the piecewise polynomial spaces for velocity and pressure independently. Rather, to obtain a stable discretization, the pair of spaces needs to satisfy a compatibility condition known as the Ladyzhenskaya–Babuška–Brezzi (LBB) or inf-sup condition (Braess, 2007; John, 2016); the condition, in essence, states that the velocity space must be sufficiently large compared to the pressure space. A common, stable choice is the “Taylor–Hood” space (Taylor and Hood, 1973) that uses piecewise quadratic elements for the velocity and piecewise linear elements for the pressure.³

Yet, there are many more combinations than just the Taylor–Hood choice one could consider (and that are used in practical applications). Specifically, among the conforming velocity elements⁴ on simplex meshes, we can consider the following choices:

• P₁, the space of piecewise linear, continuous elements;

• $P_{\mathrm{1}}^{+}$, the same space as above but with the addition of a “bubble function” on each cell that is a polynomial of order d+1 (where d is the number of space dimensions) and is zero on the faces of the cell (see for example Sect. 3.6.1 of John, 2016);

• P₂, the space of piecewise quadratic, continuous elements;

• $P_{\mathrm{2}}^{+}$, the same space as before but enriched with cubic bubble functions on each cell.

For the pressure, common choices that match those for the velocity above are

• P₀, the space of piecewise constant and consequently discontinuous elements;

• P₁, the space of piecewise linear, continuous elements;

• P₋₁, the space of piecewise linear but discontinuous elements.

Not all combinations of these are stable (that is, satisfy the LBB condition). Table 1 illustrates which combinations are stable and can consequently be used. At the same time, not all of the combinations are useful; for example, it makes perhaps little sense to use high-order polynomials for the velocity when using P₀ for the pressure because the latter might limit the convergence order of the former. As a consequence, we will here only consider a subset of the combinations.

Table 1 A summary of which simplex element combination is stable. A dash indicates that the element combination is not stable, whereas a check mark indicates that it is (John, 2016; Arnold et al., 1984; Reddy and Gartling, 2010). Circles indicate element combinations we consider in this study.

Figure 1 A graphical representation of the elements and their degrees of freedom we consider herein. Filled dots indicate locations where velocity degrees are defined, whereas open circles indicate where pressure degrees of freedom are defined. The figure does not reflect whether the shape function associated with a degree of freedom is continuous across cell boundaries.

3.2 Element combinations used in this study

Concretely, we will show results for the following.

• $P_{\mathrm{1}}^{+}\times P_{\mathrm{1}}$. This element is also often called the “MINI” element. It has not been widely used in the geodynamics community, with the noticeable exception of Zlotnik et al. (2007) in 2D and Tommasi et al. (2009) in 3D.

  Because the bubble degrees of freedom only couple to the other degrees on one cell, it can easily be eliminated from the overall linear system through “static elimination” or “static condensation”, making the element as cheap as P₁×P₁ (but stable!) (see for example Braess, 2007). Since we are mostly interested in questions of accuracy, we will use this element combination but not make use of static elimination in our implementation.

• P₂×P₀. We could not find any example of this element's use in the (geodynamical) literature. Nonetheless, we decided to include it in this study to document its performance (see previous section) as it is a cheap and stable element with discontinuous pressure and easily constructible from typical building blocks available in many finite element codes.

• P₂×P₁. This element is commonly called the “Taylor–Hood” element (see also footnote ³ above). It is used in geodynamics in, for example, the Fluidity (Davies et al., 2011) and TerraFerma⁵ codes (Wilson et al., 2017). It is also used in Schubert and Anderson (1985) and Cuffaro et al. (2020).

  This element corresponds to the widely used Q₂×Q₁ space on hypercube cells that is used, for example, in the Aspect code (Kronbichler et al., 2012; Heister et al., 2017).

• $P_{\mathrm{2}}^{+}\times P_{-\mathrm{1}}$. This element is often referenced as the “Crouzeix–Raviart” element (Dabrowski et al., 2008; Gresho and Sani, 2000).⁶ It has a discontinuous pressure, leading to local mass conservation. The relatively large pressure space requires the augmentation of the P₂ velocity space by bubble functions to guarantee stability, but – just like in the case of the $P_{\mathrm{1}}^{+}\times P_{\mathrm{1}}$ space above – the bubble degrees of freedom can be removed by static elimination.

  This element is used in geodynamics in, for example, Poliakov and Podlachikov (1992) to study the deformation of the surface above a rising diapir. It is also used in the MILAMIN code (Dabrowski et al., 2008) and in LaCoDe (de Montserrat et al., 2019).

  The closest analog to this element on hypercube elements is $Q_{\mathrm{2}}\times P_{-\mathrm{1}}$ used for example in May et al. (2015).

All of these choices are represented graphically in Fig. 1. In our numerical results below, we will compare these choices against the Q₂×Q₁ and $Q_{\mathrm{2}}\times P_{-\mathrm{1}}$ elements on hypercube cells, as we have found these to be the best choice in the first part of this study (Thieulot and Bangerth, 2022).

3.3 Alternative elements and element combinations and alternative mappings

There are many more choices one could consider beyond the ones discussed in the previous section. For example, the following come to mind:

• Nearly all of the elements listed above have analogues with higher polynomial degrees. For example, the Taylor–Hood element P₂×P₁ can be generalized to $P_{k+\mathrm{1}}\times P_k$ with k>1; all of these combinations are known to be stable and at least theoretically result in higher convergence rates. At the same time – see the discussion in Sect. 3.2 of Thieulot and Bangerth (2022) – the lack of regularity of solutions in typical geodynamic applications makes these choices unattractive: they are more expensive without delivering higher accuracy because the solution is not smooth enough to actually allow for higher convergence orders. As a consequence, we will not consider higher polynomial degrees herein than the ones mentioned previously.

• There are variations in the spaces above in which a P₁ pressure space is enriched by piecewise constants, yielding the P₁+P₀ space. The resulting element, when paired with a sufficiently large velocity space, is then mass conserving.

• Another variation is to replace a P₂ velocity space by a P₁ space on a once-refined mesh. This is commonly referred to as the “P₁isoP₂” space (Bercovier and Pironneau, 1979). The original intent in developing this element was to re-use parts of existing implementations, as well as the robustness of linear elements (for example the fact that they always attain their minima and maxima at node points, unlike higher-order shape functions).

• There are also numerous nonconforming velocity spaces in which the velocity is not continuous and which can then be made convergent either through penalty terms or by requiring structural properties such as that the velocity is at least continuous at face midpoints (see Gresho and Sani, 2000, or John, 2016, and references therein).

While perhaps useful, these alternatives are not widely used in geodynamics, and we will consequently not consider them herein. Given the conclusions we will come to in Sect. 5, one can also (retroactively) speculate that at least the nonconforming elements will not be competitive with the best elements we will find in the numerical results in Sect. 4. This is because most of the nonconforming elements were developed with the specific purpose to scope out how small one can make elements (in terms of degrees of freedom), dating back to a time when that was a prime consideration given how small computer memory was at the time rather than with the purpose of coming up with accurate and universally robust elements. Indeed, the “small” elements we consider herein will prove to be unacceptable for some reason or other below. Of course, whether the speculation that nonconforming elements are not competitive is in fact true would make for an interesting topic for follow-up work.

A separate issue we defer to a later study is how the choice of mapping from the reference cell (triangles or quadrilaterals in our 2D examples herein) affects the accuracy. For elements with quadratic velocity shape functions, it would not be unreasonable to use quadratic (i.e., “isoparametric”) mappings resulting in curved edges. In contrast, herein we only ever use straight-edged elements. Indeed, one could completely separate the polynomial degree of the mapping from that of the finite element – for example, Aspect by default uses quartic (i.e., “supraparametric”) mappings to accurately resolve curved surfaces (Kronbichler et al., 2012). It is not entirely clear how the choice of higher-order mappings would affect accuracy: practical experience shows that using higher-order mappings results in smaller errors when the geometry is curved (and one might conjecture that this would also be the case if one could resolve internal boundaries). At the same time, there remain open theoretical questions about the stability of the usual Stokes elements when curved boundaries are used (see for example Chilton and Suri, 2000). In the end, we believe that the choice of mapping is orthogonal to the choice of element, as we see no reason that an element that is not competitive with the best elements we identify here when using straight edges should become competitive when using curved edges. Rather, we consider the current study as a “filter” that allows us to identify which elements are competitive and which are not; we will then leave it to a later study to determine how they can be used with higher-order mappings.

3.4 Computational setup

For the numerical simulations shown in the following sections, we use the elements mentioned in Sect. 3.2. Because there is no reason to believe that the elements we choose perform differently in three space dimensions, we restrict our computations to 2D benchmarks because (i) these computations are substantially cheaper and (ii) it is far easier to observe convergence rates accurately in 2D: it is possible to reach much higher mesh resolutions and, consequently, get more data points in the asymptotic range where errors strictly follow 𝒪(h^α) rates, where h is the mesh size and α describes the convergence rate.

The finite element method requires the computations of integrals, for which we will use quadrature with a number of points that guarantees exact integration as long as coefficients are constant. For example, when using the Taylor–Hood element with quadratic shape functions, we use a quadrature formula with six quadrature points per triangle, arranged in the usual fashion of Gauss-type schemes. The result of the finite element integration is then a matrix for the Stokes system that is passed to a linear solver. Although advanced linear solvers are usually preferable for geodynamical codes (e.g., Kronbichler et al., 2012; May et al., 2015; Clevenger et al., 2020; Clevenger and Heister, 2021), we here resort to building the whole Stokes matrix as a sparse array and use a direct solver provided via the SciPy package.⁷ None of the computational experiments we perform herein presents the problem of a velocity nullspace; consequently, after solving the linear system, we normalize the pressure by subtracting a constant so that the average pressure is zero. All of these steps were implemented in a Python code written for the purposes of this study.

From the velocity and pressure fields computed via the procedure described above, we can then compute errors by subtracting the exact solutions (where known) and applying appropriate norms. In order to make results comparable to those in Thieulot and Bangerth (2022), we show these errors as a function of the “mesh size” h. For a specific cell K, we define its size as $h_K=\sqrt{\mathrm{2}\text{ area}\left(K\right)}$ for triangles and $h_K=\sqrt{\text{area}\left(K\right)}$ for quadrilaterals. As a consequence, the cell sizes h_K are the same for a volume subdivided into quadrilaterals and for one in which every quadrilateral is then further subdivided into two triangles; when using corresponding polynomial degrees, these two meshes generally have the same (or approximately the same) number of degrees of freedom, and the definitions of h_K above then guarantee comparability of results. In practice, on structured meshes, all cells have the same h_K; on the unstructured meshes we will use, they are approximately equal. We therefore only report results as a function of h, which we define as the average value of the h_K values.

Figure 2 Meshes used in this benchmark:(a) an example of a structured quadrilateral mesh. (b) An example of a structured triangular mesh. (c) An example of an unstructured triangular mesh. (d) An example of the unstructured triangular mesh used for the SolVi benchmark of Sect. 4.3; note the nodes (and joining edges) aligned on the quarter circle at the bottom left highlighted in red.

Clearly, different elements will have different costs for the same value of h. However, in practice, the actual cost depends on many specifics of the element choice as well as the linear solver used; consequently, we present our comparisons primarily in terms of the mesh size h rather than the number of unknowns N or any other measure primarily because h is what we used before and because theoretical results about convergence rates that readers are likely familiar with are typically shown in terms of powers of h. However, Sect. 4.1 also contains a brief discussion of run times for the various element choices.

In this study we present results obtained on structured and unstructured meshes. Structured meshes are obtained by tessellating the domain with N_x×N_y quadrilaterals (and in practice setting N_y=N_x for simplicity) as shown in Fig. 2a. For simplex meshes, these quadrilaterals are then cut along a diagonal. In order to avoid very anisotropic meshes and potentially problematic cases where three vertices of a triangle would be on the boundary (see Boffi et al., 2012, or Cioncolini and Boffi, 2019, for reasons to avoid this situation), simplex meshes are built so that we vary the direction of splitting quadrilaterals as shown in Fig. 2b.

We create unstructured simplex meshes by creating meshes via the Triangle module, which is a Python wrapper around Jonathan Richard Shewchuk's 2D quality mesh generator and Delaunay triangulator library (Shewchuk, 1996, 2014) based on a target mesh size (see Fig. 2c). Sequences of unstructured meshes are always created de novo rather than by refinement of the previous mesh, since successive refinement results in block structured meshes. For the particular case of the SolVi benchmark in Sect. 4.3, we instruct Triangle to place a number of nodes along a quarter circle to match the discontinuity in the coefficients of the benchmark (see Fig. 2d).

Figure 3 Donea and Huerta benchmark: velocity (a, c) and pressure (b, d) errors as a function of (average) mesh size for structured (a, b) and unstructured meshes (c, d).

4 Numerical results

Having so set the scene, let us now turn to quantitative evaluations of the performance of the elements discussed in the previous section. Specifically, in the current section, we will use carefully selected benchmarks that are widely used in the literature to assess geodynamics software and that we have already used (at least in parts) in the first part of this paper. More complete descriptions of the Donea–Huerta, SolCx, and SolVi benchmarks shown in Sects. 4.1–4.4, along with visual depictions of their solutions, can be found in the first part of this work (Thieulot and Bangerth, 2022). We refer the reader there for more details rather than repeating them here. The Rayleigh–Taylor benchmark of Sect. 4.5 is outlined in its own section. We end with a comparison of convergence rates between the different elements in Sect. 4.6.

4.1 The Donea and Huerta manufactured solution benchmark

The setup for this benchmark – originally described in Donea and Huerta (2003) – considers a situation where the solution is described by smooth polynomials and where the coefficients in the Stokes equations are all constant. The solution is driven by a (nonphysical) gravity field. Given the smooth solution, the different elements ought to all reach their theoretically optimal convergence rate. We use this benchmark, among other reasons, to verify the correctness of our implementations.

We show results in Fig. 3 that illustrate the accuracy with which the various discretizations approximate the exact solution. It shows that – on both structured and unstructured meshes – the discretizations that use piecewise quadratic polynomials for the velocity and linear polynomials for the pressure reach their expected velocity error of $\Vert \mathit{u}-{\mathit{u}}_h{\Vert }_{L_{\mathrm{2}}}=\mathcal{O}\left(h^{\mathrm{3}}\right)$, whereas the others only achieve $\Vert \mathit{u}-{\mathit{u}}_h{\Vert }_{L_{\mathrm{2}}}=\mathcal{O}\left(h^{\mathrm{2}}\right)$. The latter category includes the P₂×P₀ discretization that uses a large number of degrees of freedom for the velocity but achieves an error only smaller than the $P_{\mathrm{1}}^{+}\times P_{\mathrm{1}}$ by a factor despite having a number of degrees of freedom roughly 4 times higher (in 2D); conversely, it has approximately the same number of degrees of freedom as the P₂×P₁ elements but an error about 2 orders of magnitude larger on the finest meshes. The figure also shows that the P₂×P₀ element fares even worse in approximating the pressure, being substantially less accurate than the far cheaper $P_{\mathrm{1}}^{+}\times P_{\mathrm{1}}$ element (at least on structured meshes).

The figure also shows that, at least in terms of accuracy as a function of mesh size h (and consequently number of unknowns), the best-performing discretization is the P₂×P₁ element and that it produces errors quite close to the Q₂×Q₁ and $Q_{\mathrm{2}}\times P_{-\mathrm{1}}$ elements we have found to be best on quadrilaterals. Finally, the figure shows that at least for some of the elements, unstructured meshes can lead to errors nearly an order of magnitude worse than structured meshes with the same mesh size.

In practice, of course, accuracy is only one indicator of performance. A different way to measure how well the different elements perform is to measure the time required to solve a given problem. As a consequence, let us also report run times for the different benchmarks as a function of the number of degrees of freedom using a laptop with an Intel Core i7-7820HQ CPU run at 2.90 GHz with 32 GB of memory. We caution that unlike the results obtained in the first part of the paper (where we used the Aspect code), the run times shown here are obtained with a test code written in Python rather than a production code written in C++.

Figure 4 Donea and Huerta benchmark: run time to assemble the linear system (a) and to solve the linear system (b) for the different element combinations as a function of the number of degrees of freedom on a sequence of meshes.

Figure 5 Donea and Huerta benchmark: run time to assemble the linear system (a) and to solve the linear system (b) for the different element combinations as a function of the mesh size (which is inversely proportional to the square root of the number of cells) on a sequence of meshes.

Figures 4 and 5 show run times for the two dominant operations – the assembly of the linear system and the solution of the linear system – as a function of the number of unknowns and the mesh size h, respectively. The figures show that the cost of both of these operations is, in essence, a function of the number of degrees of freedom of an element combination rather than the specific details of how an element's shape functions are defined. As a consequence, the costs of different elements only differ by a (modestly sized) constant factor, rather than leading to different rates. What this also implies is that if we had shown the results of, say, Fig. 3 as a function of run time instead of mesh size h, curves would only have been moved up and down, but they would have retained their relative convergence rates, and elements shown there with a higher convergence rate will also have a higher convergence rate as a function of run time. A secondary observation is that the run times for the various Taylor–Hood-type elements are not substantially different from each other; as a consequence, run time is not a criterion that will help us choose one of these variants over the others.

Similar observations will hold for the benchmarks in the following sub-sections, and we will consequently omit comparisons of run times there.

4.2 SolCx

SolCx is a substantially more difficult benchmark to solve since it involves a viscosity that jumps by a factor of 10⁶ along the vertical mid-line of the domain. This results in a nearly discontinuous pressure as well as a kink in the velocity along this line. With properly aligned meshes, some elements can resolve these singularities, though this of course makes the benchmark not representative of real-world situations where the locations and directions of jumps in the viscosity of geodynamic models can typically not be predicted a priori and may change with time. Elements using continuous pressures consequently exhibit poor convergence. This benchmark is widely used in many geodynamical papers (e.g., Zhong, 1996; Duretz et al., 2011; Kronbichler et al., 2012; Thielmann et al., 2014; de Montserrat et al., 2019; Thieulot and Bangerth, 2022).

Figure 6 SolCx: velocity (a, c) and pressure (b, d) errors as a function of (average) mesh size for structured (a, b) and unstructured meshes (c, d).

Figure 6 shows the approximation errors we obtain for this benchmark. It illustrates the difficulties elements with continuous pressure (such as Q₂×Q₁ and P₂×P₁, specifically as opposed to $Q_{\mathrm{2}}\times P_{-\mathrm{1}}$ and $P_{\mathrm{2}}^{+}\times P_{-\mathrm{1}}$) have with this benchmark: they all only achieve a convergence rate of 𝒪(h^0.5), reflecting the lack of regularity in the exact pressure. Notably, the convergence rate for elements using a continuous piecewise linear pressure is even worse than for the P₂×P₀ element that uses (discontinuous) piecewise constant pressures.

4.3 SolVi

Of course, the difficulties of the SolCx benchmark of the previous section are somewhat artificial given that the discontinuity in the viscosity is along a vertical line that is easily matched by the mesh (if desired). In other words, while the struggles of elements with continuous pressure are real, the fact that elements with discontinuous pressures work well on such meshes could be considered a lucky break because the jump in viscosity is aligned with the discontinuity of pressures along cell interfaces – at least on structured meshes with an even number of cells per coordinate direction, as we use here.

At the same time, this is perhaps not so. Figure 6c and d already suggest that elements with discontinuous pressure spaces can adequately resolve the discontinuous pressure even on unstructured meshes, where the jump in viscosity in the SolCx benchmark is no longer aligned with cell interfaces. The SolVi benchmark we consider in this section illustrates this in more detail. It models a situation where the viscosity inside an inclusion is 1000 times larger than outside the inclusion and where we make no attempt at resolving this boundary with the structured mesh – similar to realistic situations of slab subduction or other cases of large and perhaps dynamically changing viscosity jumps that cannot practically be resolved using the meshes in use. The setup is identical to the one in Thieulot and Bangerth (2022), although here we only model one quadrant of the problem: the domain is the unit square (0,1)², the inclusion is centered on the origin, and the analytical velocity is prescribed on all sides. Schmid and Podlachikov (2003) derived a simple analytical solution for the pressure and velocity fields for this case, which was subsequently used in many other publications (Deubelbeiss and Kaus, 2008; Suckale et al., 2010; Duretz et al., 2011; Kronbichler et al., 2012; Gerya et al., 2013; Thielmann et al., 2014; de Montserrat et al., 2019).

Figure 7 SolVi: velocity (a, c) and pressure (b, d) L₂ errors as a function of (average) mesh size for structured (a, b) and unstructured meshes (c, d).

We show results for the errors in velocity and pressure in Fig. 7. In the case of structured meshes (Fig. 7a and b), the figures show that the lack of regularity in the solution, coupled with the fact that the line where this singularity occurs is not captured by the mesh, leads to a situation where all elements only obtain the convergence rate allowed by the solution rather than based on their polynomial degrees. Indeed, the quality of approximation is largely determined simply by the number of degrees of freedom an element can offer for a given mesh size h.

For unstructured meshes, we use the modified procedure shown in Fig. 2d to obtain a triangular mesh whose edges are aligned with the discontinuity of the viscosity – an approach that is admittedly artificial and would not be possible in “real” applications. The corresponding results are shown in Fig. 7c and d. They show that the alignment of cell edges to the discontinuity can recover one order of convergence for the velocity (from 𝒪(h) to 𝒪(h²) for all of the elements we investigated) and up to one order of convergence for the pressure if one uses discontinuous pressure elements. Yet, even with these aligned meshes, none of the elements achieves its optimal convergence rate.

A comparison of the curves for the structured meshes shows that, for this complex situation, the Taylor–Hood elements Q₂×Q₁ and P₂×P₁ fare the best, at least as far as “error for a given mesh size” is concerned. For the very specific discontinuity-aligned unstructured meshes, the $P_{\mathrm{2}}^{+}\times P_{-\mathrm{1}}$ pair emerges as the overall best with a quadratic pressure error convergence. The observations from these experiments also support the assertion in Sect. 3.3 (as well as the conclusions of Thieulot and Bangerth, 2022) that higher-order Taylor–Hood elements (i.e., $Q_k\times Q_{k-\mathrm{1}}$ or $Q_k\times P_{-(k-\mathrm{1})}$ on hypercubes and $P_k\times P_{k-\mathrm{1}}$ or $P_k^{+}\times P_{-k}$ on simplices, in both cases with k>2) would not yield better convergence orders despite their additional cost and are therefore not worth investigating further for geodynamic applications. This justifies why we do not consider them for this study.

4.4 The sinking block

In the SolCx and SolVi cases, the difficulty is driven by a discontinuous coefficient (the viscosity) in the differential operator of the Stokes equations (Eqs. 1 and 2). In contrast, for the sinking block benchmark, one considers a situation where a square part of the domain differs not only in viscosity, but also in density from the surrounding material – that is, in the right hand side of the equation. This results in singularities in the solution at the edges of the inclusion that have a qualitatively different behavior than that one observes in the SolCx and SolVi benchmarks. Similar or identical benchmarks can be found, for example, in May and Moresi (2008), Gerya (2019), Thieulot (2011), Mishin et al. (2022), and Schuh-Senlis et al. (2020). The current benchmark also involves having to deal with buoyancy forces (that is, a non-trivial hydrostatic pressure) that are of course the driving force for many effects in geodynamics and whose resolution is therefore important; we have found in the first part of this paper that dealing with buoyancy presented substantial problems to the stabilized Q₁×Q₁ element.

In the current benchmark, we consider a “sinker” inclusion that has a density ${\mathit{\rho }}_{\text{sinker}}={\mathit{\rho }}_{\text{fluid}}+\mathit{\delta }\mathit{\rho }$ and viscosity ${\mathit{\eta }}_{\text{sinker}}={\mathit{\eta }}^{\ast }{\mathit{\eta }}_{\text{fluid}}$. Boundary conditions are free slip on all sides, and gravity is given by $\mathit{g}=-{\mathit{e}}_y$. The domain is the unit square, and we set ρ_fluid=1 and η_fluid=1. The sinker is a square of size 0.25 × 0.25 centered at $(x_s,y_s)=(\mathrm{0},\mathrm{0.75})$. As explained in Thieulot and Bangerth (2022), “in a geodynamical context, this setup could be interpreted as a detached slab (δρ>0) or a plume head (δρ<0). As such its viscosity and density can vary (a cold slab has a higher effective viscosity than the surrounding mantle, while it is the other way around for a plume head).”

We consider two cases: (1) the fluid and the sinker densities are as described above (the “full density” case) and (2) the fluid has zero density and the density of the block is set to ρ_sinker=δρ (the “reduced density” case). The two cases of course lead to the same exact velocity field but differ in the fact that the pressure field contains a “hydrostatic” (or, in the current context, “lithostatic”) component only in the first case, whereas the background fluid (having zero density ρ_fluid) does not contribute to the pressure field in the second case. Even though the difference between the two cases is only the addition of a pressure that grows linearly with depth, the discretized equations may show an element-dependent behavior. For example, it is clear that resolving a linear pressure with an element that uses piecewise constant pressures (such as the P₂×P₀ element) will incur a substantial accuracy penalty; likewise, as shown in Thieulot and Bangerth (2022), stabilized elements yield different solutions based on whether or not the hydrostatic pressure is included.

In order to evaluate the accuracy of different elements for this benchmark, we will make use of the observation shown in Appendix A.2 of Thieulot (2011): while one can independently vary η_fluid, ρ_sinker, and η_sinker and measure $\left|v_y\right|$ in the middle of the sinker for each combination, the quantity $v^{\ast }=\left|v_y\right|{\mathit{\eta }}_{\text{fluid}}/\mathit{\delta }\mathit{\rho }$ is found to be a function of only the ratio ${\mathit{\eta }}^{\ast }={\mathit{\eta }}_{\text{sinker}}/{\mathit{\eta }}_{\text{fluid}}$. At high-enough mesh resolution, all data points then collapse onto a single line (but this may not be the case on coarse meshes: different values of the material constants may correspond to the same η^∗ but numerically result in different values of v^∗). Similarly, the normalized pressure $p^{\ast }=p/\mathit{\delta }\mathit{\rho }g{L}_b$ measured in the middle of the block is, on sufficient fine meshes, a function of η^∗ only.

We will therefore show figures that report the computed values of v^∗ and p^∗ as a function of η^∗ for all six elements. For each η^∗, we show data for $\mathit{\delta }\mathit{\rho }/{\mathit{\rho }}_{\text{fluid}}\in \mathit{\{}\mathrm{0.25}\mathit{\%},\mathrm{1}\mathit{\%},\mathrm{40}\mathit{\%}\mathit{\}}$. As mentioned, the values of v^∗ and p^∗ obtained with these three density ratios should be the same but are not the same on coarse meshes; however, this is only visible in the figures for the P₂×P₀ element where for each η^∗ up to three different values of v^∗ (one for each value of $\mathit{\delta }\mathit{\rho }/\mathit{\rho }$ considered) are apparent. We here restrict ourselves to structured meshes with resolutions of 16², 32², 64², and 128² so that element edges align with the boundary of the block.

Figure 8 The sinking block benchmark with full densities: normalized velocity v^∗ in the middle of the block (obtained for three density ratios $\mathit{\delta }\mathit{\rho }/{\mathit{\rho }}_{\text{fluid}}$) as a function of viscosity ratio η^∗. Each panel corresponds to a different mesh resolution. For the P₂×P₀ element, some of the data points fall outside of the range of the plots. (See the main text for an explanation of the scattered red dots for the P₂×P₀ element.) For reference, we also show results obtained with Aspect on a 256 × 256 mesh.

4.4.1 Full density

Figure 8 shows results for all elements and four different mesh resolutions for the case where we include the lithostatic pressure in the model. We find that, as we increase the mesh resolution, all elements but the P₂×P₀ converge to reference results obtained with the Aspect code at 256 × 256 with the $Q_{\mathrm{2}}\times P_{-\mathrm{1}}$ element. Because the overall pressure is dominated by the lithostatic component that grows linearly with depth, it is not surprising that the P₂×P₀ has a hard time approximating the pressure well; the figures show that this also translates to a poor approximation of the normalized velocity v^∗. This error becomes smaller the larger η^∗ becomes since η^∗ is a measure of the ratio of the dynamic to the lithostatic pressure.

Figure 9 The sinking block benchmark with reduced densities: normalized velocity v^∗ as a function of viscosity ratio η^∗ for various resolutions.

Figure 10 The sinking block benchmark with reduced densities: normalized pressure p^∗ in the middle of the block as a function of viscosity ratio η^∗ for various resolutions. For the P₂×P₀, $P_{\mathrm{2}}^{+}\times P_{-\mathrm{1}}$, and $Q_{\mathrm{2}}\times P_{-\mathrm{1}}$ elements with their discontinuous pressure spaces, we show p^∗ at several slightly displaced points $(x_s\pm \mathit{\delta }x,y_s\pm \mathit{\delta }y$). For the $P_{\mathrm{2}}^{+}\times P_{-\mathrm{1}}$ and $Q_{\mathrm{2}}\times P_{-\mathrm{1}}$ elements the difference is not visible at high resolution, and values for the P₂×P₀ element (red dots) fall outside the range shown here at low resolution and still show substantial differences at high resolution.

4.4.2 Reduced density

In the second case, where the density outside the inclusion is zero, the lithostatic pressure is absent and we can investigate both the dimensionless velocity (Fig. 9) and pressure (Fig. 10) in the middle of the block.

While the figure shows that the P₂×P₀ element has recovered some of its accuracy in approximating the velocity, it is unable to provide an accurate approximation of the pressure. A comparison of the convergence behavior (going from coarse to fine meshes) shows that the $P_{\mathrm{1}}^{+}\times P_{\mathrm{1}}$ element also behaves pretty poorly. The remaining elements are all of Taylor–Hood type; of these, the P₂×P₁ element with continuous pressure is substantially more accurate than the $P_{\mathrm{2}}^{+}\times P_{-\mathrm{1}}$ element with discontinuous pressure.

4.5 Rayleigh–Taylor wave benchmark

As a final comparison, we have also carried out the buoyancy-driven Rayleigh–Taylor wave instability benchmark found, for example, in Gerya (2019), Deubelbeiss and Kaus (2008), and Thieulot (2011) and for which an analytical solution of the initial growth rate can be found in Ramberg (1968). The benchmark consists of a two-layer system in a box of size L_x×L_y driven by gravity. A layer of fluid 1 (with viscosity and density η₁ and ρ₁ and thickness $h_{\mathrm{1}}=L_y/\mathrm{2}$) overlies a layer of fluid 2 (with viscosity and density η₂ and ρ₂ and thickness $h_{\mathrm{2}}=L_y/\mathrm{2}$). The interface between the two layers is disturbed by a sinusoidal displacement characterized by its amplitude Δ=0.01 and wavelength $\mathit{\lambda }=L_x/\mathrm{2}$. No-slip boundary conditions are imposed on the top and the bottom of the domain, while free slip is imposed on the sides. Gravity is set to $\mathit{g}=-{\mathit{e}}_y$. We use a mesh that is slightly distorted so as to accommodate the sinusoidal interface; however, we use straight element edges in keeping with the other benchmarks solved in this contribution. In our experiments, we specifically choose $L_x=L_y=\mathrm{1}$.

The non-horizontal interface in the setup leads to diapiric growth (illustrated in Fig. 11) whose initial vertical velocity v, at points where it is maximal, can be shown to satisfy the analytic relationship $v=-\mathrm{\Delta }K\frac{{\mathit{\rho }}_{\mathrm{1}}-{\mathit{\rho }}_{\mathrm{2}}}{\mathrm{2}{\mathit{\eta }}_{\mathrm{2}}}{h}_{\mathrm{2}}g$, with K being a dimensionless growth factor that depends on ϕ₁, ϕ₂, η₁, and η₂ (see Gerya, 2019). Instead of targeting a specific node in the domain, we evaluate v by taking the maximum vertical velocity $\left|v_y\right|$ in the domain.

Figure 11 Rayleigh–Taylor wave benchmark. The figure shows the two layers and the velocity field that results from their unequal densities, as obtained on a 64 × 64 mesh using Q₂×Q₁ elements, with η₂=10².

We then solve this benchmark with η₁=1, ρ₁=1.1, and ρ₂=1, and we vary η₂ between 10⁻² and 10² and computationally determine the vertical growth velocity at the initial time for all six element pairs and for various resolutions.

Figure 12 Rayleigh–Taylor wave benchmark. For each of the six element combinations considered in this paper, we show the maximum (absolute) vertical velocity $\left|v_y\right|$ as a function of viscosity η₂ for several different mesh resolutions.

We show results in Fig. 12. We find that all elements but the P₂×P₀ perform as expected for the range of explored viscosity values η₂: the obtained velocities fall on the analytical dashed line, with fairly little variation between element combinations that is only visible on the coarsest mesh. On the other hand, the results obtained with the P₂×P₀ combination are far from the exact values; we find that with increasing resolution, obtained velocities get closer to the analytical values (especially for η₂>1), but even at a resolution of 256 × 256 elements, v is more than a factor 2 off for η₂=0.01. This behavior is consistent with what we have found for the sinking block benchmark in Sect. 4.4.

In order to elucidate the underlying reasons, we have re-run this experiment with reduced densities (see Sect. 4.4) where we choose ρ₁=0.1 and ρ₂=0; in other words, we have subtracted a constant from the densities in both layers since flow is only driven by density differences, not densities themselves. This results in a different pressure field but the same velocity. We find that in this case the obtained velocities for the P₂×P₀ converge much faster to the analytical ones over the entire range of viscosities, as shown in Fig. 13. This is again in line with the observations for the sinking block benchmark and also matches our findings in Thieulot and Bangerth (2022) that constant-pressure elements perform poorly in buoyancy-driven flow experiments where the lithostatic pressure is dominant.

Figure 13 Rayleigh–Taylor wave benchmark. We show the maximum (absolute) vertical velocity $\left|v_y\right|$ as a function of viscosity η₂ for several different mesh resolutions for the P₂×P₀ element. Compared with the bottom-left panel in Fig. 12, we here use reduced densities, as explained in the main text.

4.6 A quantitative comparison of convergence rates

For three of the benchmarks shown in the previous sub-sections, an analytic solution is available that allowed us to compute errors. For these cases, we can also compute error rates in the L₂ norm, namely $\Vert u-u_h{\Vert }_{L_{\mathrm{2}}}\propto h^{\mathit{\alpha }}$ and $\Vert p-p_h{\Vert }_{L_{\mathrm{2}}}\propto h^{\mathit{\beta }}$. Generally, for Taylor–Hood-type elements with polynomial degree k for the velocity, one would expect $\mathit{\alpha }=k+\mathrm{1}$ and β=k if the solution is smooth, but not all elements always achieve this rate and the rate is also limited by the smoothness of the solution – see the discussion in Sects. 3.1 and 3.2 of Thieulot and Bangerth (2022).

Table 2 Observed convergence rates for the three benchmarks for which an analytic solution is available, along with the theoretically predicted optimal convergence rate for each of the elements assuming a sufficiently smooth solution. Each entry in the table consists of a pair α and β of convergence rates for the L₂ norms of the error in the velocity and pressure. “struct.”: structured meshes; “unstruct.”: unstructured meshes (for simplex elements only). Note that the optimal pressure convergence rate for the MINI element $P_{\mathrm{1}}^{+}\times P_{\mathrm{1}}$ depends on the type of mesh; on general meshes, standard finite element theory predicts it to be 1 but in certain conditions can be up to 1.5 as observed for the Donea–Huerta benchmark (see Cioncolini and Boffi, 2019, and John, 2016, p.157, for experimental evidence and Eichel et al., 2011, for an earlier theoretical investigation).

We summarize the rates we observe in our computations in Table 2, along with the optimal rate one would expect theoretically for each of these elements. The table illustrates that in cases where the solution is smooth, the Taylor–Hood-type elements achieve a higher order of convergence and, consequently, will be asymptotically more efficient than the other elements. (In practice, the results of the previous sections as well as the first part of this paper show that the Taylor–Hood-type elements are already more efficient for rather coarse meshes.) This observation will apply to the large parts of the domain in geodynamics simulations where the viscosity varies smoothly. The second observation one can draw from the table is that for cases where the solution is not smooth because the viscosity or density is discontinuous, all discretizations take a hit (unless the mesh is aligned with the discontinuity) and convergence rates are limited by the regularity of the solution.

We end this section by noting that we also computed solutions to the SolKz benchmark (Zhong, 1996) that, like the Donea–Huerta benchmark, has a smooth solution and that has been widely used in the community for similar purposes (Duretz et al., 2011; Kronbichler et al., 2012; Gerya et al., 2013; de Montserrat et al., 2019). The results are very similar to those of the Donea–Huerta case. We have also run the benchmark described in John (2016, p. 752), with results matching those provided there. In both cases, the results confirm the correctness of our implementation but do not provide any insight not already available from the benchmarks shown above; we have consequently chosen not to show these results in this contribution.

Finally, in the first part of this paper, we followed our results on benchmarks with a more concrete geodynamic application. The observations there reinforced the conclusions we had drawn based on benchmarks. Based on the results of the current paper, we see no reason to believe that solving the same application again with triangular meshes would result in any different outcomes than reported there, and we consequently omit it here.

5 Conclusions

Historically as well as recently, geodynamics codes that solve the Stokes equations have based their numerical methods on a wide variety of finite element discretizations – nearly every element ever invented has been used in some geodynamics code or other. This diversity of approaches may not always have been motivated by careful considerations of what the best method is but also by human elements such as what the implementer was familiar with or felt feasible to implement. At the same time, today's finite element discretization libraries upon which most new codes are built support a broad range of elements, both low and high order, and as a consequence, evidence-based decisions about which element to use are now both possible and called for. As a consequence, comparative studies such as the current one for simplex elements and the first part of our work in Thieulot and Bangerth (2022) for hypercube elements are both useful and necessary.

Having compared a number of possible finite element choices for the Stokes equations using a carefully selected set of benchmarks, we can summarize our findings as follows.

• The $P_{\mathrm{1}}^{+}\times P_{\mathrm{1}}$ element is not accurate enough. Although appealing on paper because of its stability and small number of unknowns, the $P_{\mathrm{1}}^{+}\times P_{\mathrm{1}}$ element is also the least accurate one in most benchmarks.

• The P₂×P₀ element can not accurately represent the lithostatic pressure. Similarly, the P₂×P₀ element is appealing because of its small pressure space and the fact that it is mass conservative due to the discontinuous pressure. At the same time, the low-order pressure does not allow the velocity to reach the optimal convergence rate, and using a piecewise constant pressure simply does not result in sufficient accuracy for applications in which an accurate representation of the lithostatic pressure field is important – e.g., for problems with pressure-dependent rheologies.

• Only Taylor–Hood elements are accurate and robust. As a consequence of these considerations, only the Taylor–Hood-type elements P₂×P₁ and $P_{\mathrm{2}}^{+}\times P_{-\mathrm{1}}$ are truly competitive across all applications we have considered. This is in line with the conclusions of the first part of this work (Thieulot and Bangerth, 2022) where we have found that on hypercube cells, only the Q₂×Q₁ and $Q_{\mathrm{2}}\times P_{-\mathrm{1}}$ elements are consistently able to provide sufficient accuracy across benchmarks. This is despite the non-trivial costs of these elements due to their large number of velocity degrees of freedom, in particular in 3D, and consequent large number of nonzero entries in system matrices – apparently this is a price one needs to pay for consistently high accuracy.

• It is not obvious which of the Taylor–Hood variants is better. Comparing between the two Taylor–Hood-like elements on triangles, the P₂×P₁ element provides a substantially better pressure approximation than the $P_{\mathrm{2}}^{+}\times P_{-\mathrm{1}}$ element for smooth solutions; in other cases, the difference is marginal, and in yet other cases the discontinuous pressure elements are substantially better. In essence, the difference is not universally large enough either way to recommend one over the other based on accuracy alone. The same is true when considering run times: all of the Taylor–Hood-type elements take about the same time in the assembly of the matrix and in the solution of the linear system; run time is consequently not a criterion to choose one Taylor–Hood variant over another. On the other hand, if local mass conservation is important, or if one wanted to use a linear solver that can exploit the block diagonal structure of the pressure mass matrix of the $P_{\mathrm{2}}^{+}\times P_{-\mathrm{1}}$ combination, then this element may have a benefit over the P₂×P₁ element.

• Per degree of freedom, hypercube elements are slightly more accurate than the corresponding simplex elements. Comparing between the Taylor–Hood-type elements on simplex and hypercube meshes, the P₂×P₁ element is typically less accurate than its counterpart Q₂×Q₁. Likewise, the $P_{\mathrm{2}}^{+}\times P_{-\mathrm{1}}$ element is typically less accurate than its counterpart $Q_{\mathrm{2}}\times P_{-\mathrm{1}}$ for smooth solutions. In neither case are the differences very large, however.

These conclusions conform with the results of the first part of this study: at the end of the day, only Taylor–Hood-type elements are consistently able to provide reliable and robust accuracy in geodynamic applications, not because they are inherently superior but because all of the other choices fail on one benchmark or other in a way that make them unsuitable for the task. It is reassuring that this conclusion is the same for simplex and hypercube elements as this hints at the universality of the properties of finite element families, regardless of the choice of reference cell.

The comparisons we have made also support another conclusion: while triangular and tetrahedral meshes have rightfully been dominant in engineering applications for their ability to mesh complex geometries (and perhaps situations in which coefficients jump at predictable locations), they are generally slightly less accurate than the corresponding finite element on quadrilateral and hexahedral cells. Taking into account that they typically lead to matrices with fewer entries, one can speculate that per unit computational cost, their performance in terms of error as a function of computational work is roughly comparable to that of hypercube cells. But, given that geodynamic applications oftentimes do not need complex geometries, this also implies that simplex meshes and elements offer no specific benefit over hypercubes and that there is no reason to abandon the common practice in the field to build codes based on hypercube cells.