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

A fluctuating boundary integral method for Brownian suspensions

We present a fluctuating boundary integral method (FBIM) for overdamped Brownian Dynamics (BD) of two-dimensional periodic suspensions of rigid particles of complex shape immersed in a Stokes fluid. We develop a novel approach for generating Brownian displacements that arise in response to the thermal fluctuations in the fluid. Our approach relies on a first-kind boundary integral formulation of a mobility problem in which a random surface velocity is prescribed on the particle surface, with zero mean and covariance proportional to the Green's function for Stokes flow (Stokeslet). This approach yields an algorithm that scales linearly in the number of particles for both deterministic and stochastic dynamics, handles particles of complex shape, achieves high order of accuracy, and can be generalized to three dimensions and other boundary conditions. We show that Brownian displacements generated by our method obey the discrete fluctuation-dissipation balance relation (DFDB). Based on a recently-developed Positively Split Ewald method [A. M. Fiore, F. Balboa Usabiaga, A. Donev and J. W. Swan, J. Chem. Phys., 146, 124116, 2017], near-field contributions to the Brownian displacements are efficiently approximated by iterative methods in real space, while far-field contributions are rapidly generated by fast Fourier-space methods based on fluctuating hydrodynamics. FBIM provides the key ingredient for time integration of the overdamped Langevin equations for Brownian suspensions of rigid particles. We demonstrate that FBIM obeys DFDB by performing equilibrium BD simulations of suspensions of starfish-shaped bodies using a random finite difference temporal integrator.Comment: Submitted to J. Comp. Phy

Robust and stable discrete adjoint solver development for shape optimisation of incompressible flows with industrial applications

PhD, 156ppThis thesis investigates stabilisation of the SIMPLE-family discretisations for incompressible flow and their discrete adjoint counterparts. The SIMPLE method is presented from typical \prediction-correction" point of view, but also using a pressure Schur complement approach, which leads to a wider class of schemes. A novel semicoupled implicit solver with velocity coupling is proposed to improve stability. Skewness correction methods are applied to enhance solver accuracy on non-orthogonal grids. An algebraic multi grid linear solver from the HYPRE library is linked to flow and discrete adjoint solvers to further stabilise the computation and improve the convergence rate. With the improved implementation, both of flow and discrete adjoint solvers can be applied to a wide range of 2D and 3D test cases. Results show that the semi-coupled implicit solver is more robust compared to the standard SIMPLE solver. A shape optimisation of a S-bend air flow duct from a VW Golf vehicle is studied using a CAD-based parametrisation for two Reynolds numbers. The optimised shapes and their flows are analysed to con rm the physical nature of the improvement. A first application of the new stabilised discrete adjoint method to a reverse osmosis (RO) membrane channel flow is presented. A CFD model of the RO membrane process with a membrane boundary condition is added. Two objective functions, pressure drop and permeate flux, are evaluated for various spacer geometries such as open channel, cavity, submerged and zigzag spacer arrangements. The flow and the surface sensitivity of these two objective functions is computed and analysed for these geometries. An optimisation with a node-base parametrisation approach is carried out for the zigzag con guration channel flow in order to reduce the pressure drop. Results indicate that the pressure loss can be reduced by 24% with a slight reduction in permeate flux by 0.43%.Queen Mary-China Scholarship Council Co-funded Scholarship No. 201206280018