PaperlayerView exact revision on arXiv

Source-linked AI summary

A Fast Nonuniform Solver for the Poisson Equation over a Disk

Charlie Pyle, Prabir Daripa

arXiv:2609.40344v1math.NAcs.NA

TL;DR

Uniform-grid FFTRR is efficient but requires evenly spaced angular samples, limiting its fit to nonuniform meshes. NUFFTRR combines nonuniform Fourier solvers with Green’s-function radial recurrences; tailored radial nodes improve accuracy by up to three orders of magnitude in the reported setting.

  • Problem

    Uniform-grid FFT Poisson solvers require evenly spaced angular samples, while interpolation of irregular samples can introduce dominant error on coarse or strongly distorted grids.

    Source: §1 · p. 2

  • Method

    NUFFTRR combines nonuniform Fourier analysis and synthesis with vectorized Green’s-function radial recurrences to support nonuniform angular and radial grids.

    Source: §3.1 · p. 4, §8 · p. 29

  • Results

    Tailoring radial nodes to resolve origin peaks or steep boundary layers improves accuracy by up to three orders of magnitude with negligible computational overhead.

    Source: §8 · p. 29

  • Takeaways & Limitations

    The radial recurrences decouple across Fourier modes and use localized one-dimensional quadrature, allowing radial nodes to be tailored without a global radial solve.

    Source: §8 · p. 29

    Large angular gaps elevate conditioning and limit stable frequency recovery relative to a well-resolved uniform grid.

    Source: §6.3 · p. 15

Abstract

from arXiv · show

We study fast numerical methods for the Poisson equation on a disk within the FFTRR (Fast Fourier Transform Radial Recurrence) framework, which is built on Green's function representations. Classical FFTRR schemes apply FFTs in the azimuthal variable and evaluate mode-by-mode radial recurrences, achieving a complexity of \(O(MN\log N)\) on an \(N\times M\) uniform grid. However, they require a uniformly spaced azimuthal mesh of \(N\) points. In this work we develop a Nonuniform (NUFFTRR) solver that admits nonuniform grids in both the radial and azimuthal directions while retaining the favorable structure of the original FFTRR formulation. The azimuthal analysis-synthesis step is implemented using either a dense NUDFT least-squares solver or one of two NUFFT-based iterative schemes: a Toeplitz solver using circulant-preconditioned conjugate gradients (PCG), and a preconditioned conjugate gradient for least squares (PCGLS) solver with Pipe--Menon density compensation. These yield azimuthal complexities \(O(N^3 + MN^2)\) for the NUDFT variant and \(O(K_{\mathrm{iter}} MN\log N)\) for the NUFFT-based variants on an \(N\times M\) grid and Krylov iteration count \(K_{\mathrm{iter}}\). Numerical experiments on nonuniform meshes demonstrate that the proposed method is robust, spectrally accurate in the azimuthal variable, and competitive in runtime with existing fast Poisson solvers. We make use of vectorization and batched BLAS/GPU-accelerated linear algebra operations to eliminate explicit loops while also formulating azimuthal and radial steps entirely in terms of dense array operations, FFTs, and NUFFTs, allowing for straightforward GPU acceleration. The implementation is released as an open-source Python package at https://github.com/CharliePyle4/NUFFTRR_Poisson, and its methodology can be extended directly to related elliptic problems such as the Helmholtz equation.

1 Introduction

Poisson solvers are needed across scientific applications, and practical meshes may be nonuniform. NUFFTRR extends the Green’s-function-based FFTRR framework to nonuniform sampling while retaining its Fourier and radial-recurrence structure.

1 Introduction

  • Practical measurement and simulation meshes may be nonuniform, motivating Poisson solvers that accommodate irregular sampling directly.Repeated large Poisson solves also arise in applications including fluid dynamics and computational geometry.Source: §1 · pp. 1–2
  • Uniform-grid FFTRR combines angular FFTs with Green’s-function radial recurrences, retaining O(MN log N) scaling while avoiding a radial grid shift at the origin.Its radial integration rule can be replaced by higher-order quadrature when additional radial accuracy is required.Source: §1 · p. 2
  • Uniform angular FFT solvers require evenly spaced samples, while interpolation of irregular data onto a uniform grid can introduce error that dominates total error.The passage identifies coarse or strongly distorted angular grids as settings where interpolation error may be especially consequential.Source: §1 · p. 2
  • NUFFTRR incorporates dense NUDFT least squares and two NUFFT-based iterative formulations into FFTRR for nonuniform angular data.The Toeplitz-PCG method uses circulant preconditioning; PCGLS uses Pipe–Menon-style density compensation, with solver choice reflecting conditioning and hardware trade-offs.Source: §1 · pp. 2–3

2 Problem Statement and Method Derivation

The disk Dirichlet problem is solved by splitting the solution into a source-driven particular solution and a homogeneous boundary correction. A Green’s-function representation and Fourier decomposition turn the computation into radial integrals evaluated by recurrences.

2 Problem Statement and Method Derivation

  • The problem is the Poisson equation on a disk with prescribed Dirichlet boundary values.The domain is the open disk B(0, R).Source: §2 · p. 3
  • The solution is decomposed as u = v + w, where v accounts for the source and w solves the homogeneous equation with boundary data g − v.This makes the correction enforce the required boundary values after the particular solution is computed.Source: §2 · p. 3
  • The particular solution v is represented by integrating the source against the two-dimensional free-space Green’s function.The Green’s function describes the influence at an observation point of a unit source at another point.Source: §2 · p. 3
  • Fourier decomposition separates the computation into azimuthal analysis and synthesis plus one-dimensional radial integrals evaluated using closed-form recurrences.The naive quadrature evaluation of the original integral would incur O(N^4) complexity on the grid.Source: §2 · p. 4

3 Azimuthal Component

The azimuthal stage maps samples to Fourier modes and back, using FFTs on uniform grids and inverse-transform solvers on nonuniform grids. Three nonuniform solver options balance computational cost and robustness to mesh distortion.

3.1 Matrix Formulation and the Analysis Problem

  • The Fourier matrix maps modal coefficients to samples, while analysis on irregular nodes requires solving an inverse problem because the matrix is no longer unitary and its adjoint is not its inverse.The implementation uses a symmetric N-mode range and the Moore–Penrose pseudoinverse for the matrix formulation.Source: §3.1 · pp. 4–5

3.2 Azimuthal Mesh Algorithms

  • On uniform angular grids, batched FFTs analyze the source and boundary data and synthesize the solution at total azimuthal cost O(MN log N).The same Fourier operator is shared across radial columns when the angular nodes are shared.Source: §3.1 · p. 4, §3.2 · p. 5
  • Dense NUDFT solves the shared nonuniform analysis system by batched QR or SVD, with complexity O(N^3 + MN^2).Optimized dense BLAS operations have negligible constant overhead and are reported as exceptionally fast in practice for N ≤ 128.Source: §3.2 · p. 5
  • For larger nonuniform grids, Toeplitz-PCG uses circulant preconditioning on mildly perturbed meshes, while PCGLS with Pipe–Menon density compensation targets more strongly distorted grids.The method selection depends on angular resolution, geometric distortion, and computational considerations.Source: §3.2 · pp. 5–6
  • The Toeplitz normal operator can be applied using 2N-point FFT convolutions, and its circulant preconditioner is most effective for mildly deformed meshes.Angular density weights compensate for clustering to restore approximate discrete orthogonality.Source: §3.2 · p. 6
  • Strong clustering, large gaps, or systematic distortions can degrade circulant preconditioning and slow PCG, while the lumped normal operator can amplify roundoff errors.PCGLS instead applies paired forward and adjoint NUFFTs directly and tracks residuals in physical space.Source: §3.2 · pp. 6–7
  • PCGLS uses Pipe–Menon fixed-point updates to build density-compensation weights, and two iterations are reported as sufficient for stable weights in practice.Weights damp oversampled clusters and boost sparse regions; the solve is performed across radial columns in batch.Source: §3.2 · p. 7
  • After radial recurrences, a batched type-2 NUFFT synthesizes physical solution values on the nonuniform angular grid.The shared node set is used across the radial columns.Source: §3.2 · pp. 5, 7

4 Radial Recurrences

Fourier-mode Green’s-function representations reduce the disk problem to inner and outer radial integrals. The radial recurrences are vectorized with prefix operations, and the Neumann case modifies the boundary correction while requiring compatibility and a reference value.

4 Radial Recurrences

  • The modal Green’s-function representation splits each radial integral at the observation radius into inner and outer intervals, enabling stepwise radial recurrences.The representation gives the Dirichlet solution modes using boundary coefficients and the particular solution modes.Source: §4 · pp. 7–8

4.1 Vectorized 1D Radial Recurrences

  • The original sequential radial sweeps are rewritten as first-order recurrences and closed forms so they can be evaluated with prefix products and sums rather than explicit radial loops.The recurrence multipliers have magnitude at most one in the stated formulation, limiting geometric amplification of accumulated perturbations.Source: §4.1 · pp. 8–9
  • The vectorized outward sweep starts at the first positive radial node to avoid evaluating the coordinate-singular ratio at the origin.The zero-mode value at the origin is enforced analytically.Source: §4.1 · p. 9
  • For real-valued source data, positive-mode particular solutions are populated by conjugate symmetry after computing the nonpositive modes.The implementation evaluates radial recurrences only for n ≤ 0.Source: §4.1 · p. 9

4.2 Neumann Boundary Conditions

  • For Neumann conditions, the particular solution and its radial recurrences remain unchanged; only the harmonic boundary correction is modified to impose the flux.The implementation computes the nonzero-mode boundary correction across modes and radii by broadcasting in O(MN) operations.Source: §4.2 · pp. 9–10
  • A Neumann solution requires global compatibility of the data, and the additive constant is fixed by prescribing the outer-boundary reference value u0(R).The implementation fixes this gauge directly rather than rejecting or projecting flux perturbed by quadrature approximation.Source: §4.2 · p. 10

5 Quadrature Schemes

The quadrature schemes improve radial accuracy while retaining local recurrences on the original radial grid. Their two-step Simpson updates can be factorized for vectorized evaluation, and nonuniform spacing is handled through local weights.

5 Quadrature Schemes

  • Direct two-step Simpson increments use existing neighboring circles, avoiding auxiliary radii while improving on the trapezoidal rule’s second-order accuracy.The modified recurrences link consecutive grid points using integrals evaluated on the original mesh.Source: §5 · pp. 10–11
  • Because Simpson updates link only indices of the same parity, even and odd radial subsequences can be solved independently and factorized with prefix products and sums.The same closed-form factorization applies to both outward and inward recurrence branches.Source: §5 · pp. 11–12
  • Nonuniform Simpson weights are computed from adjacent interval lengths, and the resulting local integrals enter the recurrence framework without changing its factorization.For 1/2 < h1/h0 < 2, the endpoint weights are strictly positive; outside this range some local weights may be negative although the rule remains exact for quadratic polynomials.Source: §5 · p. 12
  • On uniform meshes, the three-point rule is fourth-order; on smoothly graded meshes it retains fourth-order global convergence, while arbitrary unequal spacing is generally third-order.The smoothly graded condition is h1 − h0 = O(H^2); arbitrary unequal spacing has h1 − h0 = O(H).Source: §5 · p. 12
  • The decoupled radial formulation permits node concentration near boundary layers or origin singularities without a global radial solve or increased asymptotic complexity.The same recurrence framework accommodates arbitrary radial spacing.Source: §5 · p. 12
  • The integral formulation handles the origin analytically, enabling systematic higher-order quadrature without the stability issues associated with polar finite-difference singularities.The formulation evaluates boundary conditions at r1 = 0 smoothly and robustly.Source: §5 · p. 12

6 NUFFTRR Algorithm

NUFFTRR combines azimuthal Fourier analysis and synthesis with radial recurrences to solve the disk Poisson problem on uniform or nonuniform grids. Its solver variants trade moderate-resolution efficiency against large-grid scaling, while accuracy and stability depend on sampling and radial quadrature choices.

6.1 Algorithm

  • The algorithm takes source and boundary data on radial nodes with shared uniform or irregular azimuthal nodes, then returns solution values across the disk.Boundary input is Dirichlet data or Neumann flux with a reference value.Source: §6.1 · p. 13
  • The real-valued implementation symmetrizes modes on nonuniform even-N grids, a small approximation when the highest retained modes are sufficiently resolved.General complex-valued data require independent treatment of positive and negative modes.Source: §6.1 · p. 13
  • Azimuthal analysis uses FFTs on uniform meshes, dense NUDFT least squares or iterative NUFFT solvers on nonuniform meshes, followed by corresponding FFT or NUFFT synthesis.The nonuniform iterative options are Toeplitz PCG and PCGLS with Pipe–Menon density weights.Source: §6.1 · pp. 13–14

6.2 Computational Complexity

  • The radial operations in Steps 2–6 require O(MN) work, below the angular transform stages’ asymptotic costs.These steps include quadrature increments, recurrence sweeps, mode combination, and boundary corrections.Source: §6.2 · p. 14
  • Dense NUDFT is optimal for moderate resolutions (N ≤128), while the iterative NUFFT solvers scale with near-linear complexity on large grids.The dense solver’s moderate-resolution advantage is attributed to minimal BLAS overhead.Source: §6.2 · p. 14

6.3 Accuracy

  • On uniform meshes, smooth data can reach machine-precision angular resolution with modest N, but unresolved oscillations alias; radial error instead depends on quadrature order.Trapezoidal quadrature gives O(δr^2); Simpson is nominally O(δr^4), but can be O(δr^3) for n = 0, ±1 because of the origin startup interval.Source: §6.3 · p. 15
  • Nonuniform angular nodes worsen conditioning, and gaps exceeding the nominal Nyquist interval degrade stable frequency recovery regardless of total N.Toeplitz PCG is fastest among iterative variants but its preconditioner works best on mildly deformed meshes; PCGLS is a robust alternative for large, severely distorted grids.Source: §6.3 · p. 15
  • Radial node clustering can restore optimal fourth-order convergence and reduce error by up to three orders of magnitude without increasing the O(M) radial operation count.The stated gains apply to solutions with localized steep gradients when nodes are adapted to those regions.Source: §6.3 · p. 15

7 Experiments

The experiments examine accuracy and runtime across uniform and nonuniform test problems, including how angular and radial resolution, mesh adaptation, and solver choice affect results. They also assess CPU/GPU performance tradeoffs and how smoothness, domain radius, and boundary discontinuities relate to accuracy.

Uniform Problem 1

  • For the smooth uniform problem, error is essentially independent of N and converges quadratically with M using trapezoidal quadrature; Simpson quadrature gives cubic convergence.With trapezoidal quadrature, Neumann errors are about an order of magnitude higher than Dirichlet errors, though both retain second-order convergence; within Simpson, the two boundary conditions yield nearly identical accuracy.Source: Table 2 · p. 16, §7 · pp. 16–17, Table 3 · p. 18

Uniform Problem 2

  • For α = 20, angular resolution N ≥ 256 is needed to resolve oscillations and recover second-order radial convergence; for α = 5, even N = 64 suffices.At α = 5, doubling M reduces error by a factor of four; at α = 20, refining M alone cannot correct the large errors on coarse angular grids.Source: Table 4 · p. 18, §7 · pp. 18–19

Nonuniform Problem 1 – Jittered Angular Grid

  • On jittered grids, NUFFT and NUDFT match accuracy and third-order radial convergence, while spline-resampled FFT errors plateau through N = 512.The interpolation gap subsides at N = 1024, where the FFT approaches the direct solvers’ accuracy; spline interpolation otherwise obstructs radial convergence.Source: §7 · pp. 19, 21, Table 5 · p. 20, Figure 4 · p. 20
  • At N = 1024, M = 1024, NUFFT Toeplitz runs in 0.491 s versus NUDFT’s 0.869 s, while spline-resampled FFT is fastest overall.At N = 128, M = 1024, NUFFT and NUDFT achieve relative L2 error 2.23 × 10^-6 in 0.280 s and 0.217 s, respectively, while the FFT takes 0.022 s but has relative error 6.54 × 10^-2; matching their accuracy with the FFT requires N = 1024 and takes 0.136 s.Source: Figure 5 · p. 21, §7 · p. 21

Nonuniform Problem 2 – Structured Grid

  • For the structured distorted grid, NUFFT and NUDFT match accuracy, while spline-resampled FFT errors remain limited even at N = 256.At N = 64, the direct solvers’ error is governed by radial quadrature; at N ≤ 128 the FFT error stalls across M. At N = 64, M = 128, NUDFT is faster than PCGLS on both CPU (0.0125 s vs. 0.1076 s) and GPU (0.0395 s vs. 1.3938 s); matching direct-solver accuracy requires the FFT to use N = 1024.Source: Table 6 · p. 22, §7 · pp. 22–23, Table 7 · p. 23

Nonuniform Problem 3 – Radial Concentration (Origin)

  • At N = 32, squared radial clustering reaches L∞ ≈ 4.59 × 10^-8 at M = 257, nearly 300× more accurate than uniform spacing.Both nonuniform meshes sustain fourth-order Simpson convergence, versus approximately third-order convergence for the uniform mesh; radial adaptation adds no runtime penalty, and CPU execution is faster at this small scale.Source: §7 · pp. 23, 25, Table 8 · p. 24, Table 9 · p. 24

Nonuniform Problem 4 – Radial Concentration (Boundary)

  • At fixed N = 32, tuning γ to the boundary-layer thickness lets the sinh mesh attain highest accuracy across tested radial resolutions, with L∞ ≈ 7.63 × 10^-8 at M = 257 and up to three orders of magnitude better than uniform spacing.Source: §7 · p. 25

Nonuniform Problem 4

  • At M = 33, the sinh-adapted radial mesh removes the uniform mesh's localized error peak near r ≈ 0.95 and maintains up to two orders of magnitude lower error across the domain.The uniform mesh degrades to approximately O(M^-2.8) convergence because fixed steps fail to resolve the steep boundary layer as r approaches 1.Source: Figure 10 · p. 26, Table 10 · p. 26, §7 · p. 27

CPU and GPU timings

  • On CPU, Uniform FFT is fastest overall; NUDFT wins at low N, but NUFFT variants overtake it at large N, with Toeplitz faster than PCGLS.For low fixed N, direct NUDFT multiplication avoids NUFFT setup costs; for large N, NUFFT scaling dominates, and Toeplitz avoids PCGLS's repeated iterative evaluations.Source: Figure 11 · p. 27, §7 · p. 28
  • On GPU, Uniform FFT remains fastest overall; NUDFT is competitive at small N, while Toeplitz is the fastest nonuniform solver at moderate-to-large N.Batched NUDFT across radial rings can keep runtime nearly flat as M increases, while PCGLS is hindered by host-device synchronization and repeated kernel launches.Source: §7 · pp. 28–29, Figure 12 · p. 28

Other testing

  • Other tests found expected trapezoidal and Simpson radial convergence despite mild smoothness loss, while accuracy on varying-radius disks depended on effective radial spacing and joint M,N selection.The varying-radius tests found better accuracy for smaller disks and worse accuracy for larger disks, even when N was doubled.Source: §7 · p. 29
  • For a sharply peaked solution, under-resolved angular sampling caused localized aliasing that persisted despite radial refinement until N resolved the peak, after which convergence became global.Source: §7 · p. 29
  • With discontinuous Dirichlet data, pointwise error remained near angular discontinuities as N grew, while error away from them continued to converge.Source: §7 · p. 29

8 Conclusion

NUFFTRR combines nonuniform Fourier analysis with vectorized radial recurrences to solve disk Poisson problems on nonuniform radial and azimuthal grids. The conclusion reports accuracy gains from tailored radial nodes and resolution-dependent trade-offs among azimuthal solvers.

8 Conclusion

  • Tailoring radial nodes improves accuracy by up to three orders of magnitude for origin peaks or steep boundary layers, with negligible overhead and no global radial solve.The radial recurrences decouple across Fourier modes and use localized one-dimensional quadrature.Source: §8 · p. 29
  • Dense NUDFT is fastest on coarse-to-moderate grids (N ≤128), while iterative NUFFT methods surpass it at larger resolutions (N ≥256).Toeplitz PCG performs rapid circular convolution on jittered meshes; PCGLS supports convergence under more structured clustering.Source: §8 · p. 29
  • NUFFTRR supports nonuniform radial and azimuthal discretizations, using nonuniform Fourier analysis and closed-form, vectorized radial recurrences that eliminate explicit mode-by-mode loops.The array operations map to batched GPU execution.Source: §8 · p. 29
  • The operator-decoupling structure extends the framework to annular and exterior domains and related elliptic systems such as the Helmholtz equation.Source: §8 · p. 29
Loading 2609.40344v1…