Global convection in Earth's mantle : advanced numerical methods and extreme-scale simulations

The thermal convection of rock in Earth's mantle and associated plate tectonics are modeled by nonlinear incompressible Stokes and energy equations. This dissertation focuses on the development of advanced, scalable linear and nonlinear solvers for numerical simulations of realistic instantaneous mantle flow, where we must overcome several computational challenges. The most notable challenges are the severe nonlinearity, heterogeneity, and anisotropy due to the mantle's rheology as well as a wide range of spatial scales and highly localized features. Resolving the crucial small scale features efficiently necessitates adaptive methods, while computational results greatly benefit from a high accuracy per degree of freedom and local mass conservation. Consequently, the discretization of Earth's mantle is carried out by high-order finite elements on aggressively adaptively refined hexahedral meshes with a continuous, nodal velocity approximation and a discontinuous, modal pressure approximation. These velocity--pressure pairings yield optimal asymptotic convergence rates of the finite element approximation to the infinite-dimensional solution with decreasing mesh element size, are inf-sup stable on general, non-conforming hexahedral meshes with "hanging nodes,'' and have the advantage of preserving mass locally at the element level due to the discontinuous pressure. However, because of the difficulties cited above and the desired accuracy, the large implicit systems to be solved are extremely poorly conditioned and sophisticated linear and nonlinear solvers including powerful preconditioning techniques are required. The nonlinear Stokes system is solved using a grid continuation, inexact Newton--Krylov method. We measure the residual of the momentum equation in the H⁻¹-norm for backtracking line search to avoid overly conservative update steps that are significantly reduced from one. The Newton linearization is augmented by a perturbation of a highly nonlinear term in mantle's rheology, resulting in dramatically improved nonlinear convergence. We present a new Schur complement-based Stokes preconditioner, weighted BFBT, that exhibits robust fast convergence for Stokes problems with smooth but highly varying (up to 10 orders of magnitude) viscosities, optimal algorithmic scalability with respect to mesh refinement, and only a mild dependence on the polynomial order of high-order finite element discretizations. In addition, we derive theoretical eigenvalue bounds to prove spectral equivalence of our inverse Schur complement approximation. Finally, we present a parallel hybrid spectral--geometric--algebraic multigrid (HMG) to approximate the inverses of the Stokes system's viscous block and variable-coefficient pressure Poisson operators within weighted BFBT. Building on the parallel scalability of HMG, our Stokes solver demonstrates excellent parallel scalability to 1.6 million CPU cores without sacrificing algorithmic optimality.Computational Science, Engineering, and Mathematic

Implicit time integration for high-order compressible flow solvers

The application of high-order spectral/hp element discontinuous Galerkin (DG) methods to unsteady compressible flow simulations has gained increasing popularity. However, the time step is seriously restricted when high-order methods are applied to an explicit solver. To eliminate this restriction, an implicit high-order compressible flow solver is developed using DG methods for spatial discretization, diagonally implicit Runge-Kutta methods for temporal discretization, and the Jacobian-free Newton-Krylov method as its nonlinear solver. To accelerate convergence, a block relaxed Jacobi preconditioner is partially matrix-free implementation with a hybrid calculation of analytical and numerical Jacobian.The problem of too many user-defined parameters within the implicit solver is then studied. A systematic framework of adaptive strategies is designed to relax the difficulty of parameter choices. The adaptive time-stepping strategy is based on the observation that in a fixed mesh simulation, when the total error is dominated by the spatial error, further decreasing of temporal error through decreasing the time step cannot help increase accuracy but only slow down the solver. Based on a similar error analysis, an adaptive Newton tolerance is proposed based on the idea that the iterative error should be smaller than the temporal error to guarantee temporal accuracy. An adaptive strategy to update the preconditioner based on the Krylov solver’s convergence state is also discussed. Finally, an adaptive implicit solver is developed that eliminates the need for repeated tests of tunning parameters, whose accuracy and efficiency are verified in various steady/unsteady simulations. An improved shock-capturing strategy is also proposed when the implicit solver is applied to high-speed simulations. Through comparisons among the forms of three popular artificial viscosities, we identify the importance of the density term and add density-related terms on the original bulk-stress based artificial viscosity. To stabilize the simulations involving strong shear layers, we design an extra shearstress based artificial viscosity. The new shock-capturing strategy helps dissipate oscillations at shocks but has negligible dissipation in smooth regions.Open Acces

Advanced Newton Methods for Geodynamical Models of Stokes Flow with Viscoplastic Rheologies

Strain localization and resulting plasticity and failure play an important role in the evolution of the lithosphere. These phenomena are commonly modeled by Stokes flows with viscoplastic rheologies. The nonlinearities of these rheologies make the numerical solution of the resulting systems challenging, and iterative methods often converge slowly or not at all. Yet accurate solutions are critical for representing the physics. Moreover, for some rheology laws, aspects of solvability are still unknown. We study a basic but representative viscoplastic rheology law. The law involves a yield stress that is independent of the dynamic pressure, referred to as von Mises yield criterion. Two commonly used variants, perfect/ideal and composite viscoplasticity, are compared. We derive both variants from energy minimization principles, and we use this perspective to argue when solutions are unique. We propose a new stress-velocity Newton solution algorithm that treats the stress as an independent variable during the Newton linearization but requires solution only of Stokes systems that are of the usual velocity-pressure form. To study different solution algorithms, we implement 2D and 3D finite element discretizations, and we generate Stokes problems with up to 7 orders of magnitude viscosity contrasts, in which compression or tension results in significant nonlinear localization effects. Comparing the performance of the proposed Newton method with the standard Newton method and the Picard fixed-point method, we observe a significant reduction in the number of iterations and improved stability with respect to problem nonlinearity, mesh refinement, and the polynomial order of the discretization.Comment: To appear in Geochemistry, Geophysics, Geosystem

A matrix-free high-order discontinuous Galerkin compressible Navier-Stokes solver: A performance comparison of compressible and incompressible formulations for turbulent incompressible flows

Both compressible and incompressible Navier-Stokes solvers can be used and are used to solve incompressible turbulent flow problems. In the compressible case, the Mach number is then considered as a solver parameter that is set to a small value, $\mathrm{M}\approx 0.1$, in order to mimic incompressible flows. This strategy is widely used for high-order discontinuous Galerkin discretizations of the compressible Navier-Stokes equations. The present work raises the question regarding the computational efficiency of compressible DG solvers as compared to a genuinely incompressible formulation. Our contributions to the state-of-the-art are twofold: Firstly, we present a high-performance discontinuous Galerkin solver for the compressible Navier-Stokes equations based on a highly efficient matrix-free implementation that targets modern cache-based multicore architectures. The performance results presented in this work focus on the node-level performance and our results suggest that there is great potential for further performance improvements for current state-of-the-art discontinuous Galerkin implementations of the compressible Navier-Stokes equations. Secondly, this compressible Navier-Stokes solver is put into perspective by comparing it to an incompressible DG solver that uses the same matrix-free implementation. We discuss algorithmic differences between both solution strategies and present an in-depth numerical investigation of the performance. The considered benchmark test cases are the three-dimensional Taylor-Green vortex problem as a representative of transitional flows and the turbulent channel flow problem as a representative of wall-bounded turbulent flows

A fully-coupled discontinuous Galerkin method for two-phase flow in porous media with discontinuous capillary pressure

In this paper we formulate and test numerically a fully-coupled discontinuous Galerkin (DG) method for incompressible two-phase flow with discontinuous capillary pressure. The spatial discretization uses the symmetric interior penalty DG formulation with weighted averages and is based on a wetting-phase potential / capillary potential formulation of the two-phase flow system. After discretizing in time with diagonally implicit Runge-Kutta schemes the resulting systems of nonlinear algebraic equations are solved with Newton's method and the arising systems of linear equations are solved efficiently and in parallel with an algebraic multigrid method. The new scheme is investigated for various test problems from the literature and is also compared to a cell-centered finite volume scheme in terms of accuracy and time to solution. We find that the method is accurate, robust and efficient. In particular no post-processing of the DG velocity field is necessary in contrast to results reported by several authors for decoupled schemes. Moreover, the solver scales well in parallel and three-dimensional problems with up to nearly 100 million degrees of freedom per time step have been computed on 1000 processors

A Parallel Geometric Multigrid Method for Adaptive Finite Elements

Applications in a variety of scientific disciplines use systems of Partial Differential Equations (PDEs) to model physical phenomena. Numerical solutions to these models are often found using the Finite Element Method (FEM), where the problem is discretized and the solution of a large linear system is required, containing millions or even billions of unknowns. Often times, the domain of these solves will contain localized features that require very high resolution of the underlying finite element mesh to accurately solve, while a mesh with uniform resolution would require far too much computational time and memory overhead to be feasible on a modern machine. Therefore, techniques like adaptive mesh refinement, where one increases the resolution of the mesh only where it is necessary, must be used. Even with adaptive mesh refinement, these systems can still be on the order of much more than a million unknowns (large mantle convection applications like the ones in [90] show simulations on over 600 billion unknowns), and attempting to solve on a single processing unit is infeasible due to limited computational time and memory required. For this reason, any application code aimed at solving large problems must be built using a parallel framework, allowing the concurrent use of multiple processing units to solve a single problem, and the code must exhibit efficient scaling to large amounts of processing units. Multigrid methods are currently the only known optimal solvers for linear systems arising from discretizations of elliptic boundary valued problems. These methods can be represented as an iterative scheme with contraction number less than one, independent of the resolution of the discretization [24, 54, 25, 103], with optimal complexity in the number of unknowns in the system [29]. Geometric multigrid (GMG) methods, where the hierarchy of spaces are defined by linear systems of finite element discretizations on meshes of decreasing resolution, have been shown to be robust for many different problem formulations, giving mesh independent convergence for highly adaptive meshes [26, 61, 83, 18], but these methods require specific implementations for each type of equation, boundary condition, mesh, etc., required by the specific application. The implementation in a massively parallel environment is not obvious, and research into this topic is far from exhaustive. We present an implementation of a massively parallel, adaptive geometric multigrid (GMG) method used in the open-source finite element library deal.II [5], and perform extensive tests showing scaling of the v-cycle application on systems with up to 137 billion unknowns run on up to 65,536 processors, and demonstrating low communication overhead of the algorithms proposed. We then show the flexibility of the GMG by applying the method to four different PDE systems: the Poisson equation, linear elasticity, advection-diffusion, and the Stokes equations. For the Stokes equations, we implement a fully matrix-free, adaptive, GMG-based solver in the mantle convection code ASPECT [13], and give a comparison to the current matrix-based method used. We show improvements in robustness, parallel scaling, and memory consumption for simulations with up to 27 billion unknowns and 114,688 processors. Finally, we test the performance of IDR(s) methods compared to the FGMRES method currently used in ASPECT, showing the effects of the flexible preconditioning used for the Stokes solves in ASPECT, and the demonstrating the possible reduction in memory consumption for IDR(s) and the potential for solving large scale problems. Parts of the work in this thesis has been submitted to peer reviewed journals in the form of two publications ([36] and [34]), and the implementations discussed have been integrated into two open-source codes, deal.II and ASPECT. From the contributions to deal.II, including a full length tutorial program, Step-63 [35], the author is listed as a contributing author to the newest deal.II release (see [5]). The implementation into ASPECT is based on work from the author and Timo Heister. The goal for the work here is to enable the community of geoscientists using ASPECT to solve larger problems than currently possible. Over the course of this thesis, the author was partially funded by the NSF Award OAC-1835452 and by the Computational Infrastructure in Geodynamics initiative (CIG), through the NSF under Award EAR-0949446 and EAR-1550901 and The University of California -- Davis