A fully adaptive, high-order, fast Poisson solver for complex two-dimensional geometries

TL;DR

The paper addresses fast, high-order solution of inhomogeneous elliptic PDEs (Poisson) on complex 2D domains. It introduces a two-region decomposition with a bulk box-code in a truncated domain and a boundary-conforming spectral strip, patched together by Laplace layer potentials to produce a global particular solution. Key contributions include an adaptive, width-tuned fictitious boundary tilde{Γ}, a quadtree-based adaptive assembly of the bulk source, a multi-domain spectral strip solver, and a layer-potential patching strategy that preserves smoothness and achieves $O(N)$ complexity. Numerical experiments on diverse multiscale geometries demonstrate 10-digit accuracy and substantial efficiency advantages over FFT-based solvers in highly multiscale regimes.

Abstract

We present a new framework for the fast solution of inhomogeneous elliptic boundary value problems in domains with smooth boundaries. High-order solvers based on adaptive box codes or the fast Fourier transform can efficiently treat the volumetric inhomogeneity, but require care to be taken near the boundary to ensure that the volume data is globally smooth. We avoid function extension or cut-cell quadratures near the boundary by dividing the domain into two regions: a bulk region away from the boundary that is efficiently treated with a truncated free-space box code, and a variable-width boundary-conforming strip region that is treated with a spectral collocation method and accompanying fast direct solver. Particular solutions in each region are then combined with Laplace layer potentials to yield the global solution. The resulting solver has an optimal computational complexity of $O(N)$ for an adaptive discretization with $N$ degrees of freedom. With an efficient two-dimensional (2D) implementation we demonstrate adaptive resolution of volumetric data, boundary data, and geometric features across a wide range of length scales, to typically 10-digit accuracy. The cost of all boundary corrections remains small relative to that of the bulk box code. The extension to 3D is expected to be straightforward in many cases because the strip ``thickens'' an existing boundary quadrature.

Figures (9)

• Figure 1: Overview of our adaptive Poisson solver for a domain $\Omega$. The left three terms comprise the particular solution $v = v_\textup{bulk} + v_\textup{strip} + v_\textup{glue}$, while $w$ is the homogeneous solution. In the left two panels, blue indicates where the particular solution is evaluated. The grey region ($\Omega_B$, which includes ${\tilde{\Omega}}$) created by well-separated box truncation generates the volume source for $v_\textup{bulk}$. In the right two panels, layer densities are shown in red, and evaluated throughout $\Omega$. See \ref{['alg:solver']}.
• Figure 1: Definition of the fictitious curve ${\tilde{\Gamma}}$ by inward normal extension. (Left) On a multiscale geometry, using a uniform strip width based on the smallest length scale results in an unnecessarily thin strip where the panel size is large. (Center) Using a width based on the largest length scale leads to self-intersection at the smallest length scales. (Right) Using a width function $h(t)$ that smoothly adapts to local panel size gives a strip that correctly resolves all geometric features.
• Figure 1: (Left) The raindrop shape is adaptively panelized, with small panels clustering in the rounded cusp. Black circles correspond to panel endpoints. The strip region (plotted as a dashed line) conforms to this adaptive panelization. The test solution is also shown and consists of a series of Gaussians with decreasing variances clustering inside the cusp. (Center) The inhomogeneity induced by the given solution is adaptively resolved on a 16th-order background quadtree with 594944.0 unknowns and truncated according to the strip refinement criterion outlined in \ref{['sec:quadtree']}. (Right) The maximum absolute error is around $10^{-10}$ over the whole domain.
• Figure 2: Schematic for creating a smooth fictitious curve that adapts to local panel size. (a) The input panelization, given as a set of high-order Gauss--Legendre nodes sampled from the original curve $\Gamma$. (b) A crude piecewise linear width function $h(t)$ is constructed from average local panel size (black circles). A rounded approximation to each kink is created, with the amount of smoothing commensurate with the local panel size (colored curves). (c) The rounded approximations are blended together into a globally smooth width function $h(t)$ using a partition of unity. (d) The original panelization is perturbed in the normal direction by $h(t)$, yielding a smooth fictitious curve ${\tilde{\Gamma}}$ that adapts to local panel size.
• Figure 2: Performance comparison between the 2D FFT (red circles) and our adaptive solver (blue squares). We benchmark our adaptive solver on a series of teardrop problems with decreasing length scale $\eta$ and a requested error tolerance of $10^{-10}$, and record runtime and memory consumption of the entire solver. For comparison, we benchmark the 2D FFT on successively refined uniform grids using MATLAB's fft2, using ten gridpoints to resolve $\eta$; see \ref{['sec:fft']}. Dashed lines indicate predicted values for the FFT. Note that our adaptive solver is solving the full Poisson problem, while the 2D FFT would only solve the bulk problem.

Theorems & Definitions (8)

• Remark 2.1
• Remark 4.1
• Remark 4.2
• Theorem 1
• Remark 4.3
• Proof 1
• Remark 4.4: Boundary regularity of classical Newton potential
• Remark 4.5