A matrix-free parallel solution method for the three-dimensional heterogeneous Helmholtz equation

TL;DR

This work tackles the difficulty of solving large-scale 3D heterogeneous Helmholtz equations by introducing a matrix-free, parallel CSLP-preconditioned Krylov framework. The key idea is to approximately invert the Complex Shifted Laplace Preconditioner using a 3D geometric multigrid cycle, while avoiding explicit assembly of global matrices through matrix-free stencil operations. The approach supports multiple Krylov solvers (GMRES, Bi-CGSTAB, IDR(s)) and leverages a blockwise MPI domain decomposition to achieve strong and weak scalability on realistic models, including the SEG/EAGE salt structure. Numerical experiments demonstrate accurate solutions with minimal pollution error and favorable parallel performance, making the method suitable for very large 3D heterogeneous Helmholtz problems. The work also provides a POP-based performance analysis to pinpoint bottlenecks and guide further optimizations.

Abstract

The Helmholtz equation is related to seismic exploration, sonar, antennas, and medical imaging applications. It is one of the most challenging problems to solve in terms of accuracy and convergence due to the scalability issues of the numerical solvers. For 3D large-scale applications, high-performance parallel solvers are also needed. In this paper, a matrix-free parallel iterative solver is presented for the three-dimensional (3D) heterogeneous Helmholtz equation. We consider the preconditioned Krylov subspace methods for solving the linear system obtained from finite-difference discretization. The Complex Shifted Laplace Preconditioner (CSLP) is employed since it results in a linear increase in the number of iterations as a function of the wavenumber. The preconditioner is approximately inverted using one parallel 3D multigrid cycle. For parallel computing, the global domain is partitioned blockwise. The matrix-vector multiplication and preconditioning operator are implemented in a matrix-free way instead of constructing large, memory-consuming coefficient matrices. Numerical experiments of 3D model problems demonstrate the robustness and outstanding strong scaling of our matrix-free parallel solution method. Moreover, the weak parallel scalability indicates our approach is suitable for realistic 3D heterogeneous Helmholtz problems with minimized pollution error.

Figures (12)

• Figure 2.1: The velocity distribution of the 3D wedge problem.
• Figure 2.2: The velocity distribution of the 3D SEG/EAGE Salt Model. The velocity varies from $1500m\per s$ to $4482m\per s$.
• Figure 3.1: Vertex-centered 3D standard coarsening.
• Figure 3.2: The allocation map of interpolation operator.
• Figure 4.1: The schematic diagram of the MPI topology for npx0 $\times$ npy0 $\times$ npz0 = $3 \times 4 \times 5$.