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

Abstract. 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 & 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₂ × P_-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.