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

Abstract In 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 Q 2 - P 1 disc 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.

[1]  S. Brune.  Evolution of stress and fault patterns in oblique rift systems: 3‐D numerical lithospheric‐scale experiments from rift to breakup , 2014 .

[2]  Laurent Jolivet,et al.  Kinematic interpretation of the 3D shapes of metamorphic core complexes , 2012 .

[3]  Patrick Amestoy,et al.  A Fully Asynchronous Multifrontal Solver Using Distributed Dynamic Scheduling , 2001, SIAM J. Matrix Anal. Appl..

[4]  Hans Muhlhaus,et al.  A Lagrangian integration point finite element method for large deformation modeling of viscoelastic geomaterials , 2003 .

[5]  Xiaoye S. Li,et al.  An overview of SuperLU: Algorithms, implementation, and user interface , 2003, TOMS.

[6]  F. A. Dahlen,et al.  CRITICAL TAPER MODEL OF FOLD-AND-THRUST BELTS AND ACCRETIONARY WEDGES , 1990 .

[7]  Alan George,et al.  Computer Solution of Large Sparse Positive Definite , 1981 .

[8]  B. Parsons,et al.  An analysis of the variation of ocean floor bathymetry and heat flow with age , 1977 .

[9]  N. Bellahsen,et al.  Buoyancy and localizing properties of continental mantle lithosphere: Insights from thermomechanical models of the eastern Gulf of Aden , 2013 .

[10]  P. Cobbold,et al.  Topography as a major factor in the development of arcuate thrust belts: insights from sandbox experiments , 2002 .

[11]  William Gropp,et al.  Performance Modeling and Tuning of an Unstructured Mesh CFD Application , 2000, ACM/IEEE SC 2000 Conference (SC'00).

[12]  C. Dohrmann,et al.  A stabilized finite element method for the Stokes problem based on polynomial pressure projections , 2004 .

[13]  Y. Podladchikov,et al.  Transition from passive to active rifting: Relative importance of asthenospheric doming and passive extension of the lithosphere , 2001 .

[14]  S. Brune,et al.  Oblique rifting of the Equatorial Atlantic: Why there is no Saharan Atlantic Ocean , 2014 .

[15]  Taras V. Gerya,et al.  Initiation of transform faults at rifted continental margins: 3D petrological-thermomechanical modeling and comparison to the Woodlark Basin , 2013, Petrology.

[16]  Maxim A. Olshanskii,et al.  An Iterative Method for the Stokes-Type Problem with Variable Viscosity , 2009, SIAM J. Sci. Comput..

[17]  Sascha Brune,et al.  The rift to break-up evolution of the Gulf of Aden: Insights from 3D numerical lithospheric-scale modelling , 2013 .

[18]  T. Gerya.  Three-dimensional thermomechanical modeling of oceanic spreading initiation and evolution , 2013 .

[19]  Roberto F. Weinberg,et al.  Polydiapirs: multiwavelength gravity structures , 1992 .

[20]  Homer F. Walker,et al.  Globally Convergent Inexact Newton Methods , 1994, SIAM J. Optim..

[21]  Robert B. Trumbull,et al.  Numerical models of crustal scale convection and partial melting beneath the Altiplano–Puna plateau , 2002 .

[22]  Homer F. Walker,et al.  Choosing the Forcing Terms in an Inexact Newton Method , 1996, SIAM J. Sci. Comput..

[23]  D. Keyes,et al.  Jacobian-free Newton-Krylov methods: a survey of approaches and applications , 2004 .

[24]  Daniele Boffi,et al.  On the quadrilateral Q2-P1 element for the Stokes problem , 2000 .

[25]  David A. Yuen,et al.  Characteristics-based marker-in-cell method with conservative finite-differences schemes for modeling geological flows with strongly variable transport properties , 2003 .

[26]  A. Poliakov,et al.  Diapirism and topography , 1992 .

[27]  Vaneeda Allken,et al.  Factors controlling the mode of rift interaction in brittle‐ductile coupled systems: A 3D numerical study , 2012 .

[28]  Hong Zhang,et al.  Hierarchical Krylov and nested Krylov methods for extreme-scale computing , 2014, Parallel Comput..

[29]  Andrew J. Wathen,et al.  Performance and analysis of saddle point preconditioners for the discrete steady-state Navier-Stokes equations , 2002, Numerische Mathematik.

[30]  C. Beaumont,et al.  Dynamical models of the role of crustal shear zones in asymmetric continental extension , 1989 .

[31]  Georg Stadler,et al.  Parallel scalable adjoint-based adaptive solution of variable-viscosity Stokes flow problems , 2009 .

[32]  S. Brune,et al.  Linking rift propagation barriers to excess magmatism at volcanic rifted margins , 2014 .

[33]  B. Huet,et al.  Links between long-term and short-term rheology of the lithosphere: Insights from strike-slip fault modelling , 2014 .

[34]  Gene H. Golub,et al.  A Note on Preconditioning for Indefinite Linear Systems , 1999, SIAM J. Sci. Comput..

[35]  J. Pasciak,et al.  Computer solution of large sparse positive definite systems , 1982 .

[36]  Jed Brown,et al.  Efficient Nonlinear Solvers for Nodal High-Order Finite Elements in 3D , 2010, J. Sci. Comput..

[37]  C. Thieulot,et al.  3D numerical modelling of graben interaction and linkage: a case study of the Canyonlands grabens, Utah , 2013 .

[38]  H. Muhlhaus,et al.  Shear banding analysis of plastic models formulated for incompressible viscous flows , 2008 .

[39]  Josef Stoer,et al.  Numerische Mathematik 1 , 1989 .

[40]  Louis Moresi,et al.  Computational approaches to studying non-linear dynamics of the crust and mantle , 2007 .

[41]  Ralph Müller,et al.  A scalable multi‐level preconditioner for matrix‐free µ‐finite element analysis of human bone structures , 2008 .

[42]  Michel Fortin,et al.  Mixed and Hybrid Finite Element Methods , 2011, Springer Series in Computational Mathematics.

[43]  Boris J. P. Kaus,et al.  A stabilization algorithm for geodynamic numerical simulations with a free surface , 2010 .

[44]  T. Gerya,et al.  Asymmetric three-dimensional topography over mantle plumes , 2014, Nature.

[45]  Hong Zhang,et al.  Sparse triangular solves for ILU revisited: data layout crucial to better performance , 2011, Int. J. High Perform. Comput. Appl..

[46]  Pavel B. Bochev,et al.  Analysis and Computation of Compatible Least-Squares Methods for div-curl Equations , 2011, SIAM J. Numer. Anal..

[47]  Gabriele Morra,et al.  A benchmark comparison of spontaneous subduction models – towards a free surface , 2008 .

[48]  Gabriel Wittum,et al.  Competence in High Performance Computing 2010 , 2012, Springer Berlin Heidelberg.

[49]  D. Sulsky,et al.  A particle method for history-dependent materials , 1993 .

[50]  T. Hughes,et al.  Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations , 1990 .

[51]  Jonathan J. Hu,et al.  Parallel multigrid smoothing: polynomial versus Gauss--Seidel , 2003 .

[52]  G. Schreurs,et al.  ANALOGUE AND NUMERICAL MODELLING OF CRUSTAL-SCALE PROCESSES , 2006 .

[53]  Marie-Pierre Doin,et al.  A comparison of methods for the modeling of thermochemical convection , 1997 .

[54]  Taras Gerya,et al.  Dynamical Instability Produces Transform Faults at Mid-Ocean Ridges , 2010, Science.

[55]  Cedric Thieulot,et al.  FANTOM: Two- and three-dimensional numerical modelling of creeping flows for the solution of geological problems , 2011 .

[56]  Clark R. Dohrmann,et al.  Stabilization of Low-order Mixed Finite Elements for the Stokes Equations , 2004, SIAM J. Numer. Anal..

[57]  Eunseo Choi,et al.  Thermomechanics of mid-ocean ridge segmentation , 2008 .

[58]  D. Blackman,et al.  Dynamics of continental rift propagation: the end-member modes , 2005 .

[59]  S. Eisenstat,et al.  Variational Iterative Methods for Nonsymmetric Systems of Linear Equations , 1983 .

[60]  Timothy A. Davis,et al.  Algorithm 832: UMFPACK V4.3---an unsymmetric-pattern multifrontal method , 2004, TOMS.

[61]  A. Poliakov,et al.  Erosion and rheology controls on synrift and postrift evolution: Verifying old and new ideas using a fully coupled numerical model , 2001 .

[62]  C. Thieulot,et al.  DOUAR: A new three-dimensional creeping flow numerical model for the solution of geological problems , 2008 .

[63]  A. Poliakov,et al.  Self-consistent rolling-hinge model for the evolution of large-offset low-angle normal faults , 1999 .

[64]  W. R. Buck,et al.  Modes of continental lithospheric extension , 1991 .

[65]  S. V. Sobolev,et al.  SLIM3D: A tool for three-dimensional thermomechanical modeling of lithospheric deformation with elasto-visco-plastic rheology , 2008 .

[66]  Stephan V. Sobolev,et al.  Modeling suggests that oblique extension facilitates rifting and continental break‐up , 2012 .

[67]  Vincent Courtillot,et al.  Propagating rifts and continental breakup , 1982 .

[68]  Allen C. Robinson,et al.  Toward an h-Independent Algebraic Multigrid Method for Maxwell's Equations , 2006, SIAM J. Sci. Comput..

[69]  Boris J. P. Kaus,et al.  Comparing thin-sheet models with 3-D multilayer models for continental collision , 2011 .

[70]  Taras Gerya,et al.  A free plate surface and weak oceanic crust produce single‐sided subduction on Earth , 2012 .

[71]  M. B. Van Gijzen,et al.  Comparison of the deflated preconditioned conjugate gradient method and algebraic multigrid for composite materials , 2012 .

[72]  Yousef Saad,et al.  A Flexible Inner-Outer Preconditioned GMRES Algorithm , 1993, SIAM J. Sci. Comput..

[73]  D. May,et al.  Numerical modelling of magma transport in dykes , 2012 .

[74]  G. Schreurs,et al.  The numerical sandbox: comparison of model results for a shortening and an extension experiment , 2006, Geological Society, London, Special Publications.

[75]  Todd Gamblin,et al.  Scaling Algebraic Multigrid Solvers: On the Road to Exascale , 2010, CHPC.

[76]  Jed Brown,et al.  pTatin3D: High-Performance Methods for Long-Term Lithospheric Dynamics , 2014, SC14: International Conference for High Performance Computing, Networking, Storage and Analysis.

[77]  Gunar Matthies,et al.  The Inf-Sup Condition for the Mapped Qk−Pk−1disc Element in Arbitrary Space Dimensions , 2002, Computing.

[78]  S. MacLachlan,et al.  Scalable robust solvers for unstructured FE geodynamic modeling applications: Solving the Stokes equation for models with large localized viscosity contrasts , 2009 .

[79]  David A. Yuen,et al.  Robust characteristics method for modelling multiphase visco-elasto-plastic thermo-mechanical problems , 2007 .

[80]  B. Hager,et al.  It takes three to tango: 2. Bubble dynamics in basaltic volcanoes and ramifications for modeling normal Strombolian activity , 2010 .

[81]  Stéphane Zaleski,et al.  Numerical simulation of Rayleigh-Taylor instability for single and multiple salt diapirs , 1992 .

[82]  P. Fullsack.  An arbitrary Lagrangian-Eulerian formulation for creeping flows and its application in tectonic models , 1995 .

[83]  D. May,et al.  Preconditioned iterative methods for Stokes flow problems arising in computational geodynamics , 2008 .

[84]  Susan Ellis,et al.  The role of boundary conditions in numerical models of subduction zone dynamics , 2011 .

[85]  F. Harlow,et al.  Numerical Calculation of Time‐Dependent Viscous Incompressible Flow of Fluid with Free Surface , 1965 .