Algebraic Multigrid (AMG) for Saddle Point Systems

We introduce an algebraic multigrid method for the solution of matrices with saddle point structure. Such matrices e.g. arise after discretization of a second order partial differential equation (PDE) subject to linear constraints. Algebraic multigrid (AMG) methods provide optimal linear solvers for many applications in science, engineering or economics. The strength of AMG is the automatic construction of a multigrid hierarchy adapted to the linear system to be solved. However, the scope of AMG is mainly limited to symmetric positive definite matrices. An essential feature of these matrices is that they define an inner product and a norm. In AMG, matrix-dependent norms play an important role to investigate the action of the smoother, to verify approximation properties for the interpolation operator and to show convergence for the overall multigrid cycle. Furthermore, the non-singularity of all coarse grid operators in a AMG hierarchy is ensured by the positive definiteness of the initial fine level matrix. Saddle point matrices have positive and negative eigenvalues and hence are indefinite. In consequence, if conventional AMG is applied to these matrices, the method will not always converge or may even break down if a singular coarse grid operator is computed. In this thesis, we describe how to circumvent these difficulties and to build a stable saddle point AMG hierarchy. We restrict ourselves to the class of Stokes-like problems, i.e. saddle point matrices which contain a symmetric positive definite submatrix that arises from the discretization of a second order PDE. Our approach is purely algebraic, i.e. it does not require any information not contained in the matrix itself. We identify the variables associated to the positive definite submatrix block (the so-called velocity components) and compute an inexact symmetric positive Schur complement matrix for the remaining degrees of freedom (in the following called pressure components). Then, we employ classical AMG methods for these definite operators individually and obtain an interpolation operator for the velocity components and an interpolation operator for the pressure matrix. The key idea of our method is to not just merge these interpolation matrices into a single prolongation operator for the overall system, but to introduce additional couplings between velocity and pressure. The coarse level operator is computed using this "stabilized" interpolation operator. We present three different interpolation stabilization techniques, for which we show that they resulting coarse grid operator is non-singular. For one of these methods, we can prove two-grid convergence. The numerical results obtained from finite difference and finite element discretizations of saddle point PDEs demonstrate the practical applicability of our approach

Why diffusion-based preconditioning of Richards equation works: spectral analysis and computational experiments at very large scale

We consider here a cell-centered finite difference approximation of the Richards equation in three dimensions, averaging for interface values the hydraulic conductivity $K=K(p)$, a highly nonlinear function, by arithmetic, upstream, and harmonic means. The nonlinearities in the equation can lead to changes in soil conductivity over several orders of magnitude and discretizations with respect to space variables often produce stiff systems of differential equations. A fully implicit time discretization is provided by \emph{backward Euler} one-step formula; the resulting nonlinear algebraic system is solved by an inexact Newton Armijo-Goldstein algorithm, requiring the solution of a sequence of linear systems involving Jacobian matrices. We prove some new results concerning the distribution of the Jacobians eigenvalues and the explicit expression of their entries. Moreover, we explore some connections between the saturation of the soil and the ill-conditioning of the Jacobians. The information on eigenvalues justifies the effectiveness of some preconditioner approaches which are widely used in the solution of Richards equation. We also propose a new software framework to experiment with scalable and robust preconditioners suitable for efficient parallel simulations at very large scales. Performance results on a literature test case show that our framework is very promising in the advance towards realistic simulations at extreme scale

Multigrid methods for complex engineering geometries and unstructured meshes

The convergence of standard multigrid methods decays significantly if locally poor quality cells are present, and it is found that the poor convergence is due to the local failure of the smoothing property. The high frequency error localised in regions of low quality cells is not eliminated by standard multigrid smoothers and persists through multigrid cycles. We propose a global–local combined smoother for the geometric multigrid to deal with engineering meshes with a small number of poor quality cells, which includes two steps: a global smoother on the whole domain, followed by a local correction on the subdomains with low quality cells. The high frequency error remaining in the low quality regions can be damped out completely by the local correction. The idea is extended to the algebraic multigrid (AMG), including both classical AMG and smoothed aggregation AMG. It is suggested that the high frequency error produced by the smoother propagates outward the low quality region on the fine grid to the neighbouring areas on the coarse grid. An algorithm to track low quality regions on the abstract coarse grid of AMG has been developed based on the information transfer between grid levels via the transfer operators. With the local correction applied on low quality regions tracked on the abstract coarse grid, the high frequency error due to low grid quality can be removed. In the smoothed aggregation AMG, the construction of the smoothed prolongation operator depends on the spectral radius of the system. However, regions of low quality cells in a mesh increase the largest eigenvalue of the linear system. We propose a shifted largest eigenvalue strategy to approximate a reasonable spectral radius to construct the smoothed prolongation. Two and three dimensional numerical experiments, from illustrate to complicated, are demonstrated to validate the proposed smoother. Elliptic type PDEs, including Poisson and elasticity problems, are solved. For each example, the performance of multigrid on a high quality mesh is also presented as a reference case, and it is shown that the poor convergence of multigrid for low quality meshes can be recovered to the reference case by the proposed smoother. A realistic thermomechanical simulation of turbomachinery problem has also been successfully solved

HAZniCS -- Software Components for Multiphysics Problems

We introduce the software toolbox HAZniCS for solving interface-coupled multiphysics problems. HAZniCS is a suite of modules that combines the well-known FEniCS framework for finite element discretization with solver and graph library HAZmath. The focus of the paper is on the design and implementation of a pool of robust and efficient solver algorithms which tackle issues related to the complex interfacial coupling of the physical problems often encountered in applications in brain biomechanics. The robustness and efficiency of the numerical algorithms and methods is shown in several numerical examples, namely the Darcy-Stokes equations that model flow of cerebrospinal fluid in the human brain and the mixed-dimensional model of electrodiffusion in the brain tissue