A scalable, matrix-free multigrid preconditioner for finite element discretizations of heterogeneous Stokes flow

Abstract

International audienceIn this paper we describe a computational methodology that is specifically designed for studying three-dimensional geodynamic processes governed by heterogeneous visco-plastic Stokes flow. The method employs a hybrid spatial discretization consisting of a View the MathML sourceQ2-P1disc mixed finite element formulation for the Stokes problem, coupled to a material-point formulation which is used for representing material state and history-dependent variables. The applicability and practicality of this methodology is realized through the development of an efficient, scalable and robust variable viscosity Stokes preconditioner. In this work, these objectives are achieved through exploiting matrix-free operators and a geometric multigrid preconditioner employing hybrid coarse level operators, Chebyshev smoothers and hybrid Krylov coarse level solvers. The robustness and parallel efficiency of this strategy is demonstrated using an idealized geodynamic model. Lastly, we apply the new methodology to study geodynamic models of continental rifting and break-up in order to understand the diverse range of passive continental margins we observe on Earth today

Extreme-Scale Multigrid Components within PETSc

“…The primary concern is how to best represent the coefficient structure on the coarse grids. When the coefficients are highly variable, hybrid strategies which employ different techniques to define the coarselevel operators have proven to yield improved time-to-solution with minimal storage overhead [16,27].…”

“…The depth of the multigrid hierarchy might be shallower than the equivalent sequential hierarchy, and thus an iterative solver on the "coarsest" level should be used. In practice it is observed that an inexact coarse-level solve up to a specified tolerance can yield an optimal multigrid preconditioner [16,28]. However, the cost of obtaining a fixed-precision inexact solve could be n k log n k and given the shallow nature of the hierarchy, this cost may be sufficient to degrade the expected O(n) time-to-solution.…”

“…For a range of problems containing less than about 3M unknowns and scaling using up to 32 cores, a comparison of using a parallel LU factorization, an iterative method, and an agglomeration strategy together with algebraic multigrid has indicated that agglomeration is a superior approach [7]. For medium-sized problems containing approximately 25M unknowns on 4096 cores, the benefits of using agglomeration compared to an iterative coarselevel solver combined with smoothed aggregation algebraic multigrid on 4096 cores was demonstrated in [16], with an overall solver speed-up of 1.8 being reported. For large-scale computations using the UG framework [19] with problem sizes in excessive of 1B unknowns, agglomeration was found essential to obtain good scalability past 4096 cores.…”

“…[5,7,16,18,19,21]. In general, the agglomeration methods require a-priori specification of when agglomeration should occur and how aggressive it should be.…”

“…On large-scale computations employing more than 100k processing cores, the cost of the coarse-level solve represented ∼15% of the total compute time. Hierarchical Krylov methods [17] were utilized as the coarse level solver in [16] to reduce the number of global reductions required, and thereby providing a better balance between communication and computation.…”