the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
| 13 Jun 2024
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 (Q2 × Q1) or discontinuous (Q2 × 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.
- Preprint
(1547 KB) - Metadata XML
- BibTeX
- EndNote
Status: final response (author comments only)
-
RC1: 'Comment on egusphere-2024-1668', Anonymous Referee #1, 07 Aug 2024
This manuscript is a useful contribution for the developers of the finite element codes for geodynamic applications. Despite various pieces of a stable discretization puzzle were circulating in the community during decades, this work (together with its first part considering quadrilateral/hexahedral elements) represents the only attempt to bring them together in an integrated picture. Potential benefits for the community include increased general awareness about selection of an efficient/robust discretization, which can in turn lead to substantial time savings in geodynamic software development.
In general, I would certainly recommend publishing this paper after the following issues/suggestions are addressed/incorporated.
[1] Using an extensive set of benchmarks, the authors demonstrate that only sufficiently rich interpolation functions (e.g. at least quadratic for the velocity) are able to perform well in a geodynamically relevant context on the simplex geometries. This result is quite expected. It is surprising, however, that little difference was found between the elements with continuous and discontinuous pressure (according to the conclusion section). This clearly contradicts previous findings of Pelletier et al. (1989), who have explicitly advocated superiority of discontinuous pressure approximations. I believe it is at least worth citing this paper in the context of this study, and discussing the difference in the conclusions.
The conclusion “there is little difference between the Taylor-Hood variants” was most like made based on the results of SolVi benchmark alone. I think this is slightly biased, since all elements performed equally bad in this benchmark. On the other hand, in the SolCx benchmark the discontinuous pressure elements performed much better compared to the continuous counterparts, even if viscosity jumps were not aligned with the element edges. I think this also deserves to be mentioned in the conclusions.
Altogether, I actually believe that both these benchmarks are not really capable of demonstrating the difference between the continuous/discontinuous pressure elements. I suggest to complement this study with a Rayleigh-Taylor instability benchmark with a high density and high viscosity layer overlaying a low density and low viscosity layer. The viscosity contrast should be large enough (e.g. comparable with the one used in SolCx benchmark).
I also recommend authors not to term a conforming Crouzeix-Raviart element as a Taylor-Hood variant, especially in the conclusions section. Historically Taylor-Hood elements explicitly implied a continuous pressure interpolation in the geodynamic community. Significant confusion may result from the sentences: "Only Taylor-Hood elements are accurate and robust" and "There is little difference between the Taylor-Hood variants". Please rephrase them.
[2] To complete the overview presented in this study I believe it is necessary to include a discourse on the use of the so-called isoparametric elements. In geodynamic codes it is quite common to interpolate the velocities and coordinates with the same shape functions (e.g. Dabrowski et al., 2008). Stable elements which include mid-side and mid-face nodes, and use at least second order shape functions, can potentially allow elements with curvilinear edges and faces. This configuration can be achieved either when this geometrical flexibility is used on purpose to fit curvilinear boundaries, or in the context of an Arbitrary Lagrangian Eulerian (ALE) advection scheme, when nodal coordinates are updated using the calculated velocities.
At least for the mixed rectangular elements with discontinuous pressure the inf-sup stability estimates are not available for the curvilinear shapes, and special measures need to be taken to restore the convergence (e.g. Chilton, and Suri, 2000). The numerical experiments with more traditional displacements-based elements also suggest that curved edges should be avoided (e.g. Lee and Bathe, 1993). To my knowledge, there is currently a lack of similar studies in the available literature with a focus placed on simplicial element shapes (triangles and tetrahedra). It is also clear that this work is not affected by this issue, since only the elements with straight edges are considered.
Nevertheless, I believe that stability of the curvilinear elements should be also covered in this study, since potential side effects might be relevant for the community. If possible, a test should be designed to demonstrate the influence of element curvature on the convergence rates. In any case the text should contain a discussion on this topic with potential remedies explicitly indicated. The latter can include, for example, the use of sub-parameric elements, e.g. when only the linear shape functions and corner nodes are used to represent the element geometry, and the element edges and faces remain straight/planar.
[3] This manuscript only presents the two-dimensional tests performed with triangular elements. I disagree with the argumentation given by the authors to exclude the three-dimensional tests and tetrahedral elements. The convergence results cannot be trivially extrapolated from 2D to 3D. It cannot be a priory guaranteed that the same type of discretization exhibits the same convergence rate, is equivalently stable, or even exists both in 2D to 3D (e.g. wedge elements). I believe this study would still profit from including tetrahedral elements.
References
Chilton, L., Suri, M., 2000. On the construction of stable curvilinear p version elements for mixed formulations of elasticity and Stokes flow. Numerische Mathematik, 86, 29-48.
Dabrowski, M., Krotkiewski, M., Schmid, D., 2008. MILAMIN: Matlab based finite element solver for large problems. Geochem. Geophys. Geosyst., 9, Q04 030
Lee, N.-S., Bathe, K.-J., 1993. Effects of element distortions on the performance of isoparametric elements. International Journal for Numerical Methods in Engineering, 36, 3553-3576.
Pelletier, D., Fortin A., Camarero, R., 1989. Are FEM solutions of incompressible flows really incompressible? (or how simple flows can cause headaches!). International Journal for Numerical Methods In Fluids, 9, 99-112.
Citation: https://doi.org/10.5194/egusphere-2024-1668-RC1 -
RC2: 'Comment on egusphere-2024-1668', Albert de Montserrat Navarro, 22 Oct 2024
This article is complementary to recent work by the same authors that focused on different choices of hypercube (quadrilateral/hexahedral) finite elements for solving the Stokes equations in geodynamic modeling. This time, the authors focus on qualitatively assessing the accuracy of simplex (triangular/tetrahedral) finite elements. Interestingly, this work clearly shows that Taylor-Hood-type elements are the most accurate —and essentially the only viable option despite their higher computational cost— for solving Stokes equations with finite elements.Given the lack of any previous exhaustive analysis of this type, I believe this is a relevant study that provides valuable insights to the computational geodynamic community. The paper is also well-written and easy to follow, and the results are clearly presented. Therefore, I recommend this paper for publication in Solid Earth. Please find below some minor comments:- The authors discuss the performance of various elements. While it’s relatively straightforward to infer which element should be faster, I believe the paper would benefit from a figure comparing the performance of the assembly and linear solve for these elements at similar resolutions for one of the benchmarks. Metrics like DoFs/s could be useful for this comparison. Although the performance may not be entirely reliable since they are implemented in Python and likely lack many optimizations, they would still provide valuable qualitative insights into the relative speeds of the different elements.- It would be fantastic if you could upload the source files to a public repository. They would be highly educational and would greatly enhance reproducibility for anyone looking to compare future benchmarks against these results.L29: "were the" should be "are the."L215: I understand that expressing the resolution as h gives a better mental picture. However, since you are testing very different elements, perhaps using degrees of freedom instead of h is a more fair comparison?L215: "h is a same" should be "h is the same."L225: "python" should be "Python."L278: I assume the inclusion viscosity is larger than that of the matrix; could you slightly rephrase this to make it fully clear to the reader? The analytical solution of this problem also allows for different inclusion radii; could you add to the text what value you used? Additionally, could you specify the size of the domain to provide a better idea of the equivalence between the number of elements and h?L279: "the (structured) mesh" should be "the structured mesh."L295: Isn’t P2 × P-1 performing better than P2 × P1 for the unstructured mesh?Figure 6: P2 × P0 (red dots) seems to have some duplicated data points.L420-424: One could perhaps also argue that unstructured simplex meshes are easier to refine instead of using oct-/quad-trees with hanging nodes.L527: The DOI of your 2022 paper is incorrect; it leads to a different paper. The correct DOI should be https://doi.org/10.5194/se-13-229-2022.Albert de MontserratCitation: https://doi.org/
10.5194/egusphere-2024-1668-RC2 -
CC1: 'Comment on egusphere-2024-1668', Marcus Mohr, 08 Nov 2024
Dear authors,
I just finished reading your preprint "On the choice of finite element for applications in geodynamics. Part II: A comparison of simplex and hypercube elements". Like Part I, I found this really interesting and very much appreciate that you extended your testing to simplex elements.
In addition to the comments by the two reviewers, I'd like to ask some simple questions and make some small comments.
- Your conjecture on the performance on non-conforming elements would be supported (at least partially) by Terrel et al. (2012) [a chapter in "The FEniCS Book"]. There different simplicial element pairs are evaluated for a smooth 2D Stokes problem. A lowest order non-conforming Crouzeix-Raviart element for velocity with a P0 pressure is also included. It does not perform well, basically similarly bad for velocity than the mini-element and even worse for pressure. Excellent for divergence, though, as to be expected.
- On line 208 you state that "None of the computational experiments we perform herein presents the problem of a velocity nullspace". However, the SolCx benchmark in Part I is described as having "free-slip on all edges", same for the Sinking Block (line 313). Would this not result in a non-trivial kernel for the velocity block in the matrix and, hence, a non-trivial nullspace?
- In the SolVi case, I assume you have just used the actual viscosity at the quadrature points? Because there was no real superior case in the different averaging settings you compared in Part I, if I remember that correctly?
- The remark on line 129 that the bubble is a "quadratic polynomial" to me seemed a little misleading? As a product of the natural coordinates should the degree of the bubble not be (dim+1)?
- Line 313 says the domain is the unit "cube". Should that be unit "square" with the latter indicating (-1,1)x(-1,1) as otherwise part of the sinker would be outside the domain, if it was (0,1)x(0,1)?
- The DOI for the 1973 Crouzeix-Raviart paper would be https://doi.org/10.1051/m2an/197307R300331 (I looked that up. Actually the footnotes on the naming/history for TH and CR were very nice :)
- "Tayloor" on p.5 should be "Taylor" same in the references.
- W.r.t. the remark (line 150) on the P2 x P0 element being "the cheapest stable element with discontinuous pressure" my colleague Fabian Böhm reminded me that there exists the EG-P0 pair, where a P1 element for velocity is enriched by a single discontinuous vector-valued component. Hence, one has only one additional degree of freedom. See e.g.
- Si et al. (2022), SINUM, https://doi.org/10.1137/21M1391353
- Si et al. (2022), arXiv, https://doi.org/10.48550/arXiv.2110.05310
Fabians master thesis contains some tests for this pairing, e.g. a comparison of EG-P0 and P2-P1 for SolVi: https://www10.cs.fau.de/publications/theses/2023/Master_B%C3%B6hmFabian.pdf
Looking forward to the final version of the paper.
Citation: https://doi.org/10.5194/egusphere-2024-1668-CC1
Viewed
| HTML | XML | Total | BibTeX | EndNote | |
|---|---|---|---|---|---|
| 228 | 76 | 22 | 326 | 14 | 13 |
- HTML: 228
- PDF: 76
- XML: 22
- Total: 326
- BibTeX: 14
- EndNote: 13
Viewed (geographical distribution)
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1