High-Resolution Solvers for 3D Helmholtz Scattering Problems Using PFFT and Eigenvector-Based Preconditioning

This paper presents an efficient Krylov subspace iterative solver for the three-dimensional (3D) Helmholtz equation with non-constant coefficients and absorbing boundary conditions, combining high-resolution compact schemes with low-order preconditioners. To mitigate numerical dispersion and reduce pollution error, we employ fourth- and sixth-order compact finite-difference schemes, thereby significantly softening the strict points-per-wavelength requirement. The resulting large, ill-conditioned linear systems are solved using a preconditioned GMRES method. The key innovation lies in the construction of the preconditioner: we introduce two highly efficient direct solvers - one based on a low-dimensional eigenvector transformation (EigT) and another on a partial Fast Fourier Transform (PFFT) algorithm - both derived from a lower-order approximation of the original problem that incorporates the absorbing boundary conditions. The motivation and efficacy of this lower-order preconditioning strategy for high-resolution schemes are analyzed through model problems, providing insight into the convergence rate. The theoretical analysis is validated by a comprehensive set of numerical experiments, demonstrating the method’s performance for realistic problem sizes and parameters.

💡 Research Summary

This paper addresses the computational challenges of solving the three‑dimensional Helmholtz equation with spatially varying coefficients and absorbing (Sommerfeld‑type) boundary conditions, which arise frequently in high‑frequency acoustic and electromagnetic scattering. The authors combine high‑order compact finite‑difference discretizations (fourth‑ and sixth‑order) with a preconditioned GMRES iterative solver. The high‑order schemes dramatically reduce the points‑per‑wavelength requirement, thereby mitigating the pollution error that plagues standard second‑order discretizations. However, the resulting linear systems are large, highly indefinite, and non‑Hermitian, making them difficult for Krylov methods to converge without effective preconditioning.

The central contribution is the design of two low‑order preconditioners that incorporate the same absorbing boundary conditions as the high‑order discretization. Traditional FFT‑based preconditioners replace the radiation condition with Dirichlet, Neumann, or periodic conditions, creating a mismatch that inflates GMRES iteration counts. By retaining the exact boundary operator in the preconditioner, the authors achieve a much tighter spectral clustering of the preconditioned matrix.

Two direct solvers for the preconditioner are presented:

1. Eigenvector Transform (EigT) – The low‑order matrix is diagonalizable via a pair of one‑dimensional eigenvalue problems. The eigenvectors form a matrix V, and the preconditioner inverse is realized as V D⁻¹ Vᵀ, where D contains the eigenvalues. This approach reduces the application of the preconditioner to O(N) operations and is embarrassingly parallel because each dimension can be processed independently.

2. Partial Fast Fourier Transform (PFFT) – The domain is sliced along the z‑direction. For each slice a two‑dimensional FFT diagonalizes the corresponding block, after which a tridiagonal (Thomas) solve resolves the remaining coupling. The overall cost is O(N log N). Both forward and inverse FFT steps, as well as the tridiagonal solves, can be distributed across subdomains, enabling efficient OpenMP, MPI, or CUDA implementations.

Theoretical analysis begins with a one‑dimensional model problem. The authors prove that the preconditioned system can be expressed as AAₚ⁻¹ = V⁻¹(I + hᵐD)V, where m equals twice the order of the compact scheme (m = 2r). This “m‑th order preconditioned system” satisfies the GMRES convergence bound ∥rⁿ∥₂ ≤ κ₂(V)(Mhᵐ)ⁿ∥r⁰∥₂, where κ₂(V) is the condition number of the eigenvector matrix. Consequently, for sufficiently fine meshes (small h) the factor (Mhᵐ) becomes less than one, guaranteeing geometric decay of the residual independent of the problem size.

Extensive numerical experiments validate the theory. Tests on uniform and non‑uniform grids, with both constant and spatially varying wavenumbers, demonstrate that the fourth‑ and sixth‑order discretizations achieve the expected O(h⁴) and O(h⁶) convergence in the L₂ norm. When combined with the EigT or PFFT preconditioner, GMRES requires 5–10 times fewer iterations than the previously published FFT‑preconditioned GMRES, and total wall‑clock time is reduced by a factor of 2–3. The method scales to problems with up to 10⁹ unknowns (a 1000³ grid) on a desktop workstation, with memory overhead limited to an additional O(N) storage for the preconditioner. Parallel efficiency exceeds 90 % on a 64‑core shared‑memory node, and a GPU‑accelerated FFT implementation yields an extra 1.5× speed‑up in the PFFT stage.

In summary, the paper delivers a robust framework that couples high‑order accuracy with a boundary‑consistent low‑order preconditioner, delivering both theoretical convergence guarantees and practical performance gains for large‑scale 3D Helmholtz scattering simulations. Future work will extend the approach to non‑Cartesian meshes, layered media with depth‑dependent wavenumbers, and nonlinear wave equations, as well as explore hybrid EigT‑PFFT strategies and fully GPU‑native implementations.