Scalable Multigrid Solver for the Helmholtz Equation: Real-Shifted Coarse Grid Correction

Abstract: We present a convergent and scalable multigrid solver for high-frequency Helmholtz equations. Standard multigrid methods do not converge for high-frequency Helmholtz problems, and a common cure is adding a complex shift and using the shifted operator as a preconditioner. Nevertheless, the complex shift prevents scalability. In this work we present a new method that achieves scalable convergence of a 3-level cycle without a complex shift. Our key idea is real-shifting the coarsest grid Galerkin operator, to correct the numerical dispersion between the grids. We show that this real-shifted coarse grid correction leads to a scalable 3-level method, for problems with 12 grid points per wavelength on the fine grid, and a convergent cycle with very few iterations for 11 grid points per wavelength, using standard point-smoothers. For problems with 10 grid points per wavelength, our method combined with a modest complex shift outperforms the standard complex shifted Laplacian method by an order of magnitude. We demonstrate wavenumber independent convergence for heterogeneous geophysical media in 2D and 3D.

• The paper presents innovative real-shifted coarse grid correction (RS-CGC) to effectively reduce inter-level phase misalignments in high-frequency Helmholtz problems.
• It employs fourth-order discretization and high-order grid-transfer operators to optimize phase alignment, achieving robust scalability across 2D and 3D benchmarks.
• Numerical results on geophysical models validate RS-CGC’s superior convergence and scalability against traditional complex-shifted methods, offering practical benefits for large-scale simulations.

Scalable Multigrid Solver for the Helmholtz Equation via Real-Shifted Coarse Grid Correction

Introduction

The paper "Scalable Multigrid Solver for the Helmholtz Equation: Real-Shifted Coarse Grid Correction" (2604.19501) introduces a new multigrid methodology addressing the principal difficulties in solving high-frequency Helmholtz problems. The complexity of these problems arises mainly from the indefiniteness of the discretized linear system and the necessity for fine meshes, often requiring at least 10 points per wavelength for accurate simulation. Achieving scalable and wavenumber-independent convergence in this regime presents substantial numerical challenges, primarily due to phase misalignments (numerical dispersion) across multigrid levels.

Standard multigrid frameworks—preconditioned with complex-shifted Laplacian (CSLP) approaches—suffer from poor scalability as the complex shift required increases with the wavenumber, and their effectiveness degrades for deep multilevel cycles or on large, heterogeneous domains. The central contribution of this paper is the development of a real-shifted coarse grid correction (RS-CGC), applied at the coarsest Galerkin grid, fundamentally aimed at reducing inter-level phase discrepancies. The approach surpasses CSLP and comparable re-discretization techniques, especially in the high-frequency regime and on large-scale 2D/3D geophysical benchmarks.

Methodology

Discretization and Multigrid Cycle Structuring

Discretization employs fourth-order compact finite difference schemes (9-point stencil in 2D, 19-point in 3D) to control truncation error and mitigate numerical dispersion. High-order interpolatory grid-transfer operators (cubic and level-dependent schemes) ensure optimal interaction between fine and coarse grids. The method is limited to 3-level multigrid cycles, balancing stencil compactness, computational feasibility, and Nyquist sampling constraints.

The real shift parameter $\alpha$ modifies the Helmholtz operator at the coarsest grid level. The fine-level operator $H_h = -\Delta_h - k^2 M_h$ is replaced on the coarsest grid with $H_h^\alpha = -\Delta_h - (\alpha k)^2 M_h$, with $\alpha$ chosen to optimally minimize inter-grid phase mismatch.

Real-Shifted Grid-to-Grid Dispersion Analysis

A pivotal aspect is the detailed grid-to-grid dispersion analysis, tailored to multigrid, contrasting with classical local Fourier analysis. The alignment of phase velocities and radii of discrete dispersion relations between fine and coarse grids is explicitly computed for each propagation direction using the derived stencils. The error metric is $$e_g(\alpha, \phi) = \frac{V_1(\phi)}{4 V_3(\alpha, \phi)} - 1$$, quantifying the relative phase error across levels, accounting for grid coarsening by a factor of 4.

The optimal real shift $\alpha$ is determined by solving $\arg\min_\alpha \max_\phi |e_g(\alpha, \phi)|$ for each grid configuration. This approach ensures minimal maximal phase error over all propagation angles, a criterion more stringent and appropriate than $L^2$-averaged error minimization.

Figure 1: Orange: Dispersion relation of the coarse grid operator with the real shift for 3D Helmholtz, compared to the fine grid dispersion relation (blue, stretched by factor 4).

Algorithmic Implementation

The RS-CGC framework is compatible as both a standalone solver and as a Krylov preconditioner. Numerical linear algebra is implemented to full scalability with parallelized Jacobi smoothers and minimal pre/post-relaxations. Galerkin coarse grid operators exploit wide stencils, leveraging GPU-friendly computation despite increased stencil bandwidth.

Numerical Results

Homogeneous and Heterogeneous 2D Benchmarks

Extensive experiments on homogeneous and linearly varying velocity fields display wavenumber-independent iteration counts (6–8 up to $1024^2$ grids for $G=12$). For $H_h = -\Delta_h - k^2 M_h$0, scalability remains robust, although iteration counts increase modestly. For $H_h = -\Delta_h - k^2 M_h$1, combining RS-CGC with a small CSLP shift ($H_h = -\Delta_h - k^2 M_h$2) yields over an order of magnitude reduction in iteration relative to classical CSLP, demonstrating the strength of the hybridization strategy.

The Marmousi model, a standard for complexity and heterogeneity, is solved efficiently with both W- and V-cycle variants, with W-cycles offering superior scalability for the largest grids.

Figure 2: Marmousi velocity model used for heterogeneous 2D testing.

Three-Dimensional Scaling and Geophysical Models

On large 3D heterogeneous domains (Overthrust model), RS-CGC achieves constant iteration counts ($H_h = -\Delta_h - k^2 M_h$3) for both $H_h = -\Delta_h - k^2 M_h$4 and $H_h = -\Delta_h - k^2 M_h$5, and wall-clock time demonstrates linear scaling with problem size. The method is memory efficient owing to the level-dependent restriction and efficient implementation of wide coarse-grid stencils.

Figure 3: The Overthrust model, a challenging 3D geophysical benchmark with substantial heterogeneity.

Comparison With Alternative Methods

The method furnishes clear numerical superiority over re-discretization-based coarse grid corrections and CSLP (with either standard or high-order transfer), particularly as the grid is refined or as heterogeneity intensifies. Both convergence factor alignment and practical iteration counts validate the theoretical grid-to-grid dispersion analysis, with computed optimal shifts closely matching empirical optima.

Theoretical and Practical Implications

The RS-CGC method provides scalable, wavenumber-independent convergence for Helmholtz problems, obviating the need for complex shifts except in the most aggressive under-resolved regimes. The grid-to-grid dispersion analysis introduces a precise criterion for real-shift selection and a rule-of-thumb for maximal grid size versus allowable phase error, which offers practical guidance for deployment on very large problems.

From a theoretical perspective, the demonstrated alignment between fine and coarsened levels suggests new pathways for generalizing phase-corrected multigrid to broader classes of indefinite PDEs and for integration with physics-informed machine learning or adaptive discretization schemes. Practically, the method's success on geophysical benchmarks in both 2D and 3D implies substantial value for seismic imaging and wave-based inversion at scale.

Future Directions

Potential extensions include:

• Formalizing tighter a priori bounds on the pollution effect and accumulation of error,
• Extension to elastic Helmholtz and Maxwell formulations,
• Aggressive parallelization leveraging GPU-enabled wide-stencil operators,
• Adaptive or locally-tuned shift selection for strongly variable coefficients.

Conclusion

This work establishes RS-CGC as a scalable, practical, and theoretically motivated multigrid solver for indefinite Helmholtz equations, demonstrating strong convergence and computational efficiency across both homogeneous and strongly heterogeneous media in two and three dimensions. The synergy of real-shifted coarse correction and rigorous dispersion alignment analysis marks a significant advancement in the design of Helmholtz solvers for large-scale scientific and engineering applications.