A high order multigrid-preconditioned immersed interface solver for the Poisson equation with boundary and interface conditions

TL;DR

This work addresses the challenge of solving the Poisson equation on complex 3D domains with immersed boundaries and interfaces while achieving high-order accuracy. It introduces a matrix-free GMRES solver for a high-order immersed discretization, preconditioned by a low-order Shortley-Weller multigrid that remains robust on arbitrary geometries and multiresolution grids. The approach delivers up to sixth-order spatial accuracy, demonstrates spectral properties favorable to preconditioning, and scales to large 3D problems on thousands of cores, with verified convergence on 3D manufactured solutions and two substantial applications. The results indicate that high-order immersed methods, coupled with SW-MG preconditioning, provide a practical, scalable pathway for complex elliptic PDEs and offer potential extensions to elasticity, incompressible flows, and fluid-structure interactions in 3D settings.

Abstract

This work presents a multigrid preconditioned high order immersed finite difference solver to accurately and efficiently solve the Poisson equation on complex 2D and 3D domains. The solver employs a low order Shortley-Weller multigrid method to precondition a high order matrix-free Krylov subspace solver. The matrix-free approach enables full compatibility with high order IIM discretizations of boundary and interface conditions, as well as high order wavelet-adapted multiresolution grids. Through verification and analysis on 2D domains, we demonstrate the ability of the algorithm to provide high order accurate results to Laplace and Poisson problems with Dirichlet, Neumann, and/or interface jump boundary conditions, all effectively preconditioned using the multigrid method. We further show that the proposed method is able to efficiently solve high order discretizations of Laplace and Poisson problems on complex 3D domains using thousands of compute cores and on multiresolution grids. To our knowledge, this work presents the largest problem sizes tackled with high order immersed methods applied to elliptic partial differential equations, and the first high order results on 3D multiresolution adaptive grids. Together, this work paves the way for employing high order immersed methods to a variety of 3D partial differential equations with boundary or inter-face conditions, including linear and non-linear elasticity problems, the incompressible Navier-Stokes equations, and fluid-structure interactions.

Figures (17)

• Figure 1: Each crossing between a grid line and the boundary ($x_c$) is used to construct ghosts points for the affected grid points (light grey) using a multidimensional interpolant constructed from a half-elliptical region of grid points (half-ellipsoidal in 3D). For immersed boundaries (a), the interpolant is constructed using imposed boundary conditions (Dirichlet, Neumann). For immersed interfaces (b), polynomials from both sides are constructed using imposed jump conditions.
• Figure 2: Truncation error of Poisson discretizations with (a) fourth order and (b) sixth order interior schemes as a function of the grid spacing $\Delta x$ for both Dirichlet and Neumann boundary conditions. In all cases the truncation error of each $(n, k)$ discretization is of order $\Delta x^{k - 2}$, indicated by dashed lines, with the Neumann cases showing a slightly lower error magnitude than their Dirichlet counterparts.
• Figure 3: Solution error (4, k) and (6, k) Poisson discretizations with both Dirichlet and Neumann boundary conditions. The black dashed lines in each plot indicate slopes of $\order{\Delta x^{n - 1}}$, $\order{\Delta x^{n}}$, and $\order{\Delta x^{n+1}}$, from top to bottom.
• Figure 4: $L_\infty$ norm convergences of the unknown boundary quantity for (4, 5) and (6, 7) discretizations with Dirichlet or Neumann boundary conditions. For Dirichlet discretizations the unknown is the boundary normal gradient, while for Neumann discretizations it is the boundary solution values.
• Figure 5: $L_\infty$ norm convergences of the solution gradient for (4, 5) and (6, 7) discretizations with Dirichlet or Neumann boundary conditions. The solution gradient is computed using the same immersed interface discretization as used for the Laplacian, using a centered finite difference stencil.