Fast Solvers for Elliptic Partial Differential Equations Based on Spectral and High-Order Methods

Abstract

This thesis develops numerical methods for the solution of elliptic partial differential equations (PDEs), all focused on delivering high-order accuracy with low computational complexity. Three algorithms for solving elliptic PDEs are presented: a fast direct solver for a global spectral method on simple two- and three-dimensional domains; a fast direct solver for a spectral element method in two dimensions; and a multigrid solver for a discontinuous Galerkin method in two and three dimensions. In Chapter 2, an optimal complexity, spectrally-accurate method is presented for the solution of Poisson's equation on a square, solid cylinder, solid sphere, and cube. The method employs a carefully chosen spectral basis to discretize Poisson's equation as a Sylvester equation with pentadiagonal matrices, which is solved using the alternating direction implicit (ADI) method. A separated spectra property in our discretizations is exploited to allow the ADI method to be used as a direct solver. In Chapter 3, we introduce a sparse spectral element method in two dimensions based on the ultraspherical spectral method, and a corresponding fast direct solver based on the hierarchical Poincaré–Steklov scheme for solving second-order linear PDEs on polygonal domains with unstructured quadrilateral or triangular meshes. The spectral element method achieves an overall computational complexity of O(p^4/h^3) for mesh size h and polynomial order p in 2D, enabling hp-adaptivity to be efficiently performed. We develop an open-source software system, ultraSEM, for flexible, user-friendly spectral element computations in MATLAB. In Chapter 4, an efficient hp-multigrid scheme is presented for local discontinuous Galerkin (LDG) discretizations of elliptic problems, formulated around the idea of separately coarsening the underlying discrete gradient and divergence operators. We show that traditional multigrid coarsening of the primal formulation leads to poor and suboptimal multigrid performance, whereas coarsening of the flux formulation leads to essentially optimal convergence and is equivalent to a purely geometric multigrid method. We show that good multigrid convergence rates are achieved in a variety of numerical tests on 2D and 3D uniform and adaptive Cartesian grids, as well as for curved domains using implicitly defined meshes and for multi-phase elliptic interface problems with complex geometry.