Xiangmin Jiao
A Multigrid Generalized Finite Difference (GFD) Method

Dept of Applied Math & Stat
Stony Brook University
Stony Brook
NY 11794
xiangmin.jiao@stonybrook.edu
Cao Lu

Finite differences are well-established numerical methods for solving differential equations, but they were limited to structured meshes. The generalized finite differences (GFDs) extend the classical finite differences to unstructured meshes with high-order accuracy. In this paper, we propose a GFD method for elliptic PDEs, and tightly integrate it with a multigrid solver. We refer to this integrated approach as multigrid GFD. Specifically, our method discretizes the PDE and the boundary conditions using GFD to obtain a sparse linear system

$$\vec{A}\vec{u}=\vec{b},$$ (1)

where $\vec{A}$ is in general asymmetric and may be singular for PDEs with Neumann or periodic boundary conditions. We solve the resulting linear system using a hybrid geometric+algebraic method, or HyGA. The core of HyGA is a geometric multigrid method, which rediscretizes the PDE using GFD over a series of hierarchical meshes, combined with an algebraic multigrid at the coarsest level.

We develop GFD methods based on a weighted least squares (WLS) formulation locally over a weighted stencil at each point, which generalizes the classical interpolation-based finite differences. We have shown recently that the WLS can deliver the same order of accuracy as the interpolation-based approximations, while delivering superior stability, flexibility, and robustness for multivariate approximations over irregular unstructured meshes. The GFD method for PDEs is formulated as follows. Consider a general form of linear, time-independent partial differential equations

$$\mathcal{P} u(\vec{x})=f(\vec{x})$$ (2)

with some Dirichlet or Neumman boundary conditions, where $\mathcal{P}$ is a linear differential operator (such as the Laplacian operator) and $\vec{x}\in\mathbb{R}^{d}$ for $d\ge1$ . Let $\Psi(\vec{x})=\{\psi_{j}(\vec{x})\}$ denote the set of Dirac delta functions at the vertices of a mesh. We obtain a linear equation for each $\psi_{j}$ as

$$\int_{\Omega}\mathcal{P} u(\vec{x})\psi_{j}(\vec{x}) d\vec{x}=\int_{\Omega}f(\vec{x})\psi_{j}(\vec{x}) d\vec{x},$$ (3)

and then the boundary conditions may modify the linear system. We then introduce a set of basis functions $\Phi(\vec{x})=\{\phi_{i}(\vec{x})\}$ to approximate $\mathcal{P}u$ as $\mathcal{P} u(\vec{x}_{0})=\vec{u}^{T}\mathcal{P} \Phi(\vec{x}_{0})\approx\vec{u}^{T}\vec{d}$ , where $\vec{d}$ is obtained from weighted least squares over the stencil around point $\vec{x}_{0}$ . Suppose $\Phi$ and $\Psi$ are both stored in column vectors. Then in ()

$$\vec{A}=\int_{\Omega}\Psi(\vec{x})\left(\mathcal{P} \Phi(\vec{x})\right)^{T} d\vec{x}$$ (4)

and $\vec{b}=\int_{\Omega}\Psi(\vec{x})f(\vec{x}) d\vec{x}$ .

To construct the geometric multigrid GFD, which is the core of our multigrid GFD approach, we derive the restriction and prolongation operators rigorously for hierarchical meshes. Specifically, let $\Phi^{(k)}$ denote the set of basis functions on the $k$ th finest mesh, and let $\Psi^{(k)}$ denote the Dirac delta functions on the $k$ th finest mesh. Let $\vec{A}^{(k)}$ denote the coefficient matrix at the $k$ th level, i.e., $\vec{A}^{(k)}=\int_{\Omega}\Psi^{(k)}(\vec{x})\left(\mathcal{P} \Phi^{(k)}(\vec{x})\right)^{T} d\vec{x}.$ In a two-level setting, let $\vec{R}_{\Phi}^{(1,2)}$ denote a restriction matrix of the functional space such that $\Phi^{(2)}=\vec{R}_{\Phi}^{(1,2)}\Phi^{(1)}$ , where $\vec{R}_{\Phi}^{(1,2)}\in\mathbb{R}^{n_{2}\times n_{1}}$ with $n_{k}=\vert\Phi^{(k)}\vert$ . Similarly, $\Psi^{(2)}=\vec{R}_{\Psi}^{(1,2)}\Psi^{(1)}$ for $\vec{R}_{\Psi}^{(1,2)}\in\mathbb{R}^{n_{2}\times n_{1}}$ . We show that the restriction matrix $\vec{R}$ and the prolongation matrix $\vec{P}$ in a two-grid method are

$$\vec{R}=\vec{R}_{\Psi}^{(1,2)} \vec{P}=\left(\vec{R}_{\Phi}^{(1,2)}\right)^{T},$$ (5)

so that $\vec{A}^{(2)}=\vec{R} \vec{A}^{(1)}\vec{P}$ . In other words, $\vec{R}$ corresponds to an $n_{2}\times n_{1}$ permutation matrix (i.e., an injection operator) that selects the entries in the residual vector corresponding to the coarse-level nodes in the fine-level mesh. $\vec{P}$ is the transpose of of the restriction operator from the basis functions $\Phi^{(1)}$ to $\Phi^{(2)}$ , and it would correspond to the interpolation of nodal values if $\Phi^{(k)}$ were composed of Lagrange basis functions. Applying these operators to adjacent grids in a multilevel method, we then obtain a geometric multigrid for GFD, where the coefficient matrices at all levels are obtained directly from GFD discretizations without performing sparse matrix-matrix multiplications. By coupling this geometric multigrid with algebraic multigrid at the coarsest level, our proposed HyGA method then combines the rigor, high accuracy and runtime-and-memory efficiency of geometric multigrid with the flexibility of algebraic multigrid, and at the same time it is relatively easy to implement.

In this work, we also investigate the convergence properties of the multigrid solver for the GFDs of elliptic PDEs. We show the multigrid solver may fail to converge for high-order GFD discretizations when using Jacobi or Gauss-Seidel iterations as smoothers. To overcome this issue, we propose new weighting schemes for WLS to improve the diagonal dominance of the linear system from the GFD, and also adopt a damped Gauss-Seidel iteration to enable more robust convergence of multigrid methods. We present the derivations and the numerical experimentations of our method for a number of test problems with various boundary conditions.