Journal cover Journal topic
Solid Earth An interactive open-access journal of the European Geosciences Union
Journal topic

Journal metrics

IF value: 2.921
IF 5-year value: 3.087
IF 5-year
CiteScore value: 4.8
SNIP value: 1.314
IPP value: 2.87
SJR value: 0.993
Scimago H <br class='widget-line-break'>index value: 38
Scimago H
h5-index value: 36
© Author(s) 2020. This work is distributed under
the Creative Commons Attribution 4.0 License.
© Author(s) 2020. This work is distributed under
the Creative Commons Attribution 4.0 License.

  19 May 2020

19 May 2020

Review status
A revised version of this preprint was accepted for the journal SE.

Pragmatic Solvers for 3D Stokes and Elasticity Problems with Heterogeneous Coefficients: Evaluating Modern Incomplete LDLT Preconditioners

Patrick Sanan1, Dave A. May2, Matthias Böllhofer3, and Olaf Schenk4 Patrick Sanan et al.
  • 1Department of Earth Sciences, ETH Zurich, Sonneggstrasse 5, Zürich 8092, Switzerland
  • 2Department of Earth Sciences, University of Oxford, South Parks Road, Oxford OX1 3AN, United Kingdom
  • 3Institute for Computational Mathematics, TU Braunschweig, D-38106 Braunschweig, Germany
  • 4Advanced Computing Laboratory, Università della Svizzera italiana, Via Buffi 6, Lugano 6900, Switzerland

Abstract. The need to solve large saddle point systems within computational Earth sciences is ubiquitous. Physical processes giving rise to these systems include porous flow (the Darcy equations), poroelasticity, elastostatics, and highly viscous flows (the Stokes equations). The numerical solution of saddle point systems is non-trivial since the operators are indefinite.

Primary tools to solve such systems are direct solution methods (exact triangular factorization) or Approximate Block Factorization (ABF) preconditioners. While ABF solvers have emerged as the state-of-the-art scalable option, they are invasive solvers requiring splitting of pressure and velocity degrees of freedom, a multigrid hierarchy with tuned transfer operators and smoothers, machinery to construct complex Schur complement preconditioners, and the expertise to select appropriate parameters for a given coefficient regime – they are far from being "black box" solvers. Modern direct solvers, which robustly produce solutions to almost any system, do so at the cost of rapidly growing time and memory requirements for large problems, especially in 3D. Incomplete LDL (ILDL) factorizations, with symmetric maximum weighted matching preprocessing, used as preconditioners for Krylov (iterative) methods, have emerged as an efficient means to solve indefinite systems. These methods have been developed within the numerical linear algebra community but have yet to become widely used in non-trivial applications, despite their practical potential; they can be used whenever a direct solver can, only requiring an assembled operator, yet can offer comparable or superior to performance, with the added benefit of having a much lower memory footprint. In comparison to ABF solvers, they only require the specification of a drop tolerance and thus provide an easy-to-use addition to the solver toolkit for practitioners.

Here, we present solver experiments employing incomplete LDL factorization with symmetric maximum weighted matching preprocessing to precondition operators, and compare these to direct solvers and ABF-preconditioned iterative solves. To ensure the comparison study is meaningful for Earth scientists, we utilize matrices arising from two prototypical problems, namely Stokes flow and quasi-static (linear) elasticity, discretized using standard mixed finite element spaces. Our test suite targets problems with large coefficient discontinuities across non-grid-aligned interfaces, which represent a common, challenging-for-solvers, scenario in Earth science applications. Our results show: (i) as the coefficient structure is made increasingly challenging (high contrast, complex topology), the ABF solver can break down, becoming less efficient than the ILDL solver before breaking down entirely; (ii) ILDL is robust, with a time-to-solution that is largely independent of the coefficient topology and mildly dependent on the coefficient contrast; (iii) the time-to-solution obtained using ILDL is typically faster than that obtained from a direct solve, beyond 10^5 unknowns; (iv) ILDL always uses less memory than a direct solve.

Patrick Sanan et al.

Interactive discussion

Status: final response (author comments only)
Status: final response (author comments only)
AC: Author comment | RC: Referee comment | SC: Short comment | EC: Editor comment
[Login for Authors/Topical Editors] [Subscribe to comment alert] Printer-friendly Version - Printer-friendly version Supplement - Supplement

Patrick Sanan et al.

Patrick Sanan et al.


Total article views: 370 (including HTML, PDF, and XML)
HTML PDF XML Total Supplement BibTeX EndNote
206 121 43 370 32 38 41
  • HTML: 206
  • PDF: 121
  • XML: 43
  • Total: 370
  • Supplement: 32
  • BibTeX: 38
  • EndNote: 41
Views and downloads (calculated since 19 May 2020)
Cumulative views and downloads (calculated since 19 May 2020)

Viewed (geographical distribution)

Total article views: 171 (including HTML, PDF, and XML) Thereof 171 with geography defined and 0 with unknown origin.
Country # Views %
  • 1



No saved metrics found.


No discussed metrics found.
Latest update: 20 Sep 2020
Publications Copernicus
Short summary
Mantle and lithospheric dynamics, elasticity, subsurface flow, and other fields involve solving indefinite linear systems. Tools include direct solvers (robust, easy to use, expensive) or advanced iterative solvers (complex, problem-sensitive). We show that a third option, ILDL preconditioners, requires less memory than direct solvers but is easy to use, as applied to 3D problems with parameters jumps. With included software, we hope to allow researchers to solve previously infeasible problems.
Mantle and lithospheric dynamics, elasticity, subsurface flow, and other fields involve solving...