Shifted HSS solvers for the indefinite Helmholtz equation

1. Problem Statement

The paper addresses the challenge of solving the indefinite Helmholtz equation at high wavenumbers ($k$) using finite element discretization. The equation is given by:
$$(-\delta + ik)^2 u - \nabla^2 u = f$$
where $\delta \ge 0$ is an absorption parameter and $k > 0$ is the wavenumber. The primary difficulty arises when $k \gg \delta$ (the high-frequency regime), particularly when $\delta = 0$ (purely indefinite).

Key Challenges:

• Pollution Effect: To maintain accuracy, the mesh size $h$ must scale such that $hk^{3/2}$ is small, leading to a massive number of degrees of freedom ($N$).
• Solver Robustness: Standard iterative solvers (like GMRES) and classical multigrid methods often fail to converge robustly as $k$ increases. Classical multigrid requires a shift $\hat{\delta} = O(k)$ to be robust, but preconditioning the original operator requires $\hat{\delta} = O(1)$ to keep the number of outer iterations independent of $k$.
• Scalability: Existing methods often struggle to scale efficiently on high-performance computing (HPC) systems for very large $k$.

2. Methodology

The authors propose a Shifted Hermitian Skew-Hermitian Splitting (HSS) iteration strategy that bridges the gap between the requirements for multigrid robustness and outer iteration robustness.

A. Core Concept: HSS Splitting

The method splits the operator into Hermitian and Skew-Hermitian parts. For a shifted operator with shift $\hat{\delta}$, the iteration is formulated as a two-step stationary method. By eliminating the intermediate step, the skew-Hermitian part appears only in its Cayley transform (a unitary operator), which allows for a contraction rate bound independent of the mesh size.

B. The Two-Stage Strategy

The authors distinguish between two cases based on the absorption parameter $\delta$:

1. Case $\delta > 0$ (Shifted System):

   • The operator is preconditioned with a shift $\hat{\delta} = O(k)$.
   • This shifted operator is suitable for standard multigrid solvers (robust convergence).
   • The HSS iteration is applied directly to this system.
   • Result: $O(k)$ HSS iterations are sufficient to reduce the error by a constant factor, independent of mesh resolution.

2. Case $\delta = 0$ (Indefinite System):

   • The original indefinite operator is preconditioned by a shifted operator with $\hat{\delta} = O(1)$.
   • The inverse of this shifted operator is approximated not by a direct solve, but by performing $O(k)$ "inner" HSS iterations.
   • These inner iterations solve a system with a larger shift ($\hat{\delta}_{inner} = O(k)$), which is efficiently solved via multigrid.
   • Result: The outer iteration (Richardson or FGMRES) becomes $k$-robust (number of iterations independent of $k$) provided the inner loop scales as $O(k)$.

C. Implementation Details

• Formulations: The method is developed for both Primal (standard $H^1$) and Mixed (first-order system) finite element formulations.
• Multigrid Integration: The inner HSS steps are solved approximately using standard Multigrid cycles (W-cycles) with Jacobi smoothers. The stencil widths and communication halos remain fixed regardless of $k$, ensuring excellent parallel scalability.
• Parallelization: The algorithm is designed for HPC using MPI (e.g., via PETSc), leveraging the fact that the whole iteration is built from standard components.

3. Key Contributions

1. Theoretical Proof of Robustness:
   • Proved that for the shifted system ($\delta > 0$), the HSS iteration converges in $O(k)$ steps with a contraction rate independent of the mesh size.
   • Proved that for the indefinite system ($\delta = 0$), using $O(k)$ inner HSS iterations to approximate the inverse of the shifted preconditioner results in an outer iteration count that is independent of both $k$ and the mesh resolution.
2. Bridging the $\hat{\delta}$ Dichotomy: The paper resolves the conflict between needing $\hat{\delta} = O(k)$ for multigrid efficiency and $\hat{\delta} = O(1)$ for preconditioner effectiveness by nesting the iterations.
3. Complexity Analysis: The authors argue that within the range of multigrid scalability, the total wallclock time complexity is $O(k)$. This is achieved because the total work scales as $O(k \cdot N)$ while the number of processors scales as $O(N)$, leading to linear time scaling with respect to the wavenumber.

4. Numerical Results

The authors validated their theory using 2D and 3D numerical experiments on the ARCHER2 supercomputer.

• Convergence Rates:
  • Inner HSS Iterations: The measured convergence rates matched the theoretical bounds tightly.
  • Outer Iterations: When using $O(k)$ inner iterations ($\theta=1$ in their notation), the number of outer Richardson or FGMRES iterations remained constant (or slightly decreasing) as $k$ increased from 16 to 1024 (2D) and 20 to 80 (3D).
  • Comparison: Using fewer inner iterations ($\theta < 1$) caused outer iterations to grow with $k$, while using more ($\theta > 1$) reduced outer iterations but increased total work unnecessarily.
• Scalability:
  • Weak Scaling: The algorithm demonstrated excellent weak scaling.
  • Wallclock Time: The time-to-solution scaled linearly with $k$ ($O(k)$), confirming the theoretical prediction.
  • FGMRES vs. Richardson: While Richardson iteration verified the theory, FGMRES significantly reduced the total number of iterations (often by more than half) in practice.
• Dimensions: Results were consistent across 2D and 3D, with the 3D case showing similar robustness despite the higher growth in degrees of freedom ($O(k^{4.5})$).

5. Significance

This work establishes a scalable, robust, and efficient solver for the indefinite Helmholtz equation at high frequencies.

• HPC Readiness: By relying on standard multigrid components with fixed stencil widths, the method avoids complex, non-standard transfer operators, making it highly suitable for modern parallel architectures.
• Overcoming Pollution: The $O(k)$ wallclock time complexity is crucial for problems where the mesh must be refined faster than $1/k$ to overcome the pollution effect, a regime where many existing solvers become prohibitively expensive.
• Generalizability: The framework is easily extendable to other wave equations (e.g., Maxwell's equations, elastodynamics) and higher-order discretizations.

In summary, the paper provides a rigorous mathematical framework and practical implementation for solving high-frequency Helmholtz problems, demonstrating that $O(k)$ wallclock time is achievable on large-scale supercomputers.