pith. sign in

arxiv: 2604.19501 · v1 · submitted 2026-04-21 · 🧮 math.NA · cs.NA

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

Rachel Yovel , Eran Treister This is my paper

Pith reviewed 2026-05-10 02:02 UTC · model grok-4.3

classification 🧮 math.NA cs.NA
keywords Helmholtz equationmultigrid methodsnumerical dispersionreal shiftscalable solvershigh-frequency wavesgeophysical modelingcoarse-grid correction
0
0 comments X

The pith

Real-shifting only the coarsest grid operator yields scalable three-level multigrid convergence for high-frequency Helmholtz equations.

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The paper shows that standard multigrid fails for high-frequency Helmholtz problems because of dispersion errors that grow with wavenumber. By applying a real shift exclusively to the Galerkin operator on the coarsest grid, the method removes those inter-grid mismatches while preserving fine-grid accuracy and avoiding the scalability penalty of complex shifts. This produces wavenumber-independent convergence in three-level cycles with ordinary point smoothers when the fine grid has at least twelve points per wavelength, and rapid convergence at eleven points per wavelength. At ten points per wavelength the same real-shift idea combined with a modest complex shift beats the usual complex-shifted Laplacian preconditioner by roughly an order of magnitude. The resulting solver works for heterogeneous media in two and three dimensions, which matters for large-scale wave simulations where iteration counts must stay bounded as problem size increases.

Core claim

Real-shifting the coarsest-grid Galerkin operator corrects numerical dispersion between grids and thereby produces a convergent, scalable three-level multigrid cycle for the Helmholtz equation that uses standard point smoothers, achieves wavenumber-independent iteration counts for discretizations with twelve points per wavelength, converges in very few iterations at eleven points per wavelength, and, when paired with a small complex shift, outperforms the standard complex-shifted Laplacian method at ten points per wavelength, all while remaining effective on heterogeneous geophysical media in two and three dimensions.

What carries the argument

The real-shifted coarsest-grid Galerkin operator, which modifies only the coarsest level with a real shift chosen to compensate for dispersion mismatch between the fine-grid Helmholtz operator and its coarse representations.

If this is right

  • A three-level cycle becomes scalable for problems discretized with twelve grid points per wavelength.
  • Convergence occurs in only a few iterations at eleven grid points per wavelength using standard smoothers.
  • At ten grid points per wavelength the real-shift approach plus a modest complex shift reduces work by roughly an order of magnitude compared with the standard complex-shifted Laplacian.
  • Wavenumber-independent convergence holds for heterogeneous media in both two and three dimensions.

Where Pith is reading between the lines

These are editorial extensions of the paper, not claims the author makes directly.

  • The same real-shift idea might extend to four or more levels while preserving the observed scalability.
  • The shift magnitude could be chosen adaptively from local wavenumber estimates to handle strongly varying media.
  • Similar dispersion corrections might improve multigrid performance for other high-frequency wave equations that currently rely on complex shifts.

Load-bearing premise

That a real shift applied solely to the coarsest grid is sufficient to remove dispersion errors across all levels in heterogeneous media without introducing instabilities or losing fine-grid accuracy.

What would settle it

Numerical experiments on three-dimensional heterogeneous problems with exactly twelve points per wavelength showing that the number of iterations grows with wavenumber or that the method diverges would falsify the scalability claim.

read the original 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.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, simulated authors' rebuttal, and a circularity audit. Tearing a paper down is the easy half of reading it; the pith above is the substance, this is the friction.

Referee Report

3 major / 2 minor

Summary. The manuscript proposes a 3-level multigrid solver for the high-frequency Helmholtz equation in which a real shift is applied exclusively to the coarsest-grid Galerkin operator to correct inter-grid numerical dispersion. It claims that this yields a scalable cycle (wavenumber-independent iteration counts) for problems discretized at 12 points per wavelength and rapid convergence at 11 points per wavelength using standard point smoothers, with further improvement over complex-shifted Laplacian preconditioners at 10 points per wavelength when a modest complex shift is added. Numerical results are presented for 2D and 3D heterogeneous geophysical media demonstrating the claimed behavior.

Significance. If the reported numerical performance holds under broader testing, the method would represent a meaningful advance for scalable Helmholtz solvers in applications such as seismic imaging, where complex shifts are known to degrade scalability. The restriction of the shift to the coarsest level preserves real arithmetic on finer grids and avoids the parameter-tuning overhead of complex shifts, provided the uniform real shift remains effective across velocity contrasts.

major comments (3)
  1. [§4, Table 2] §4 (Numerical Experiments), Table 2 and Figure 5: iteration counts are shown for heterogeneous 2D/3D cases at 11-12 ppw, but no systematic variation of the real-shift magnitude with contrast ratio or local wavenumber is reported; a single global shift is used throughout, leaving open whether residual phase errors accumulate in high-contrast regions as hypothesized in the stress-test note.
  2. [§3.1] §3.1 (Coarse-grid correction): the real-shifted Galerkin operator is defined with a constant shift parameter whose selection is not derived from an a priori dispersion analysis; the manuscript provides no proof or bound showing that this constant suffices to cancel position-dependent dispersion errors induced by spatially varying coefficients, which is load-bearing for the scalability claim.
  3. [§5] §5 (Comparison with complex-shifted Laplacian): at 10 ppw the hybrid method is reported to require an order-of-magnitude fewer iterations, yet the baseline complex-shift parameter and its scaling with wavenumber are not tabulated, preventing direct assessment of whether the improvement stems from the real-shift component or from the specific choice of the modest complex shift.
minor comments (2)
  1. [Abstract and §2] Notation for the real-shift parameter should be introduced once and used consistently; its symbol appears without prior definition in the abstract and early sections.
  2. [Figure 6] Figure captions for the 3D heterogeneous examples should explicitly state the velocity contrast ratio and the number of grid points per wavelength to allow readers to judge the regime.

Simulated Author's Rebuttal

3 responses · 1 unresolved

We thank the referee for the thorough review and valuable feedback on our manuscript. We address each of the major comments point by point below, indicating the revisions we plan to make to strengthen the paper.

read point-by-point responses

  1. Referee: [§4, Table 2] §4 (Numerical Experiments), Table 2 and Figure 5: iteration counts are shown for heterogeneous 2D/3D cases at 11-12 ppw, but no systematic variation of the real-shift magnitude with contrast ratio or local wavenumber is reported; a single global shift is used throughout, leaving open whether residual phase errors accumulate in high-contrast regions as hypothesized in the stress-test note.

    Authors: We agree that a systematic study of the real-shift parameter's dependence on contrast ratio and local wavenumber would provide additional insight. In our experiments, a single global shift was chosen based on tests that achieved the reported wavenumber-independent convergence for the presented heterogeneous media. To address this concern, we will add new numerical results in a revised §4 or an appendix, systematically varying the shift parameter in high-contrast test cases to demonstrate its robustness and show that residual phase errors do not accumulate significantly. revision: yes

  2. Referee: [§3.1] §3.1 (Coarse-grid correction): the real-shifted Galerkin operator is defined with a constant shift parameter whose selection is not derived from an a priori dispersion analysis; the manuscript provides no proof or bound showing that this constant suffices to cancel position-dependent dispersion errors induced by spatially varying coefficients, which is load-bearing for the scalability claim.

    Authors: The constant shift is motivated by dispersion analysis for the constant-coefficient case, where it corrects the numerical dispersion on the coarse grid to align with the fine-grid operator. For heterogeneous problems, we rely on numerical evidence that this fixed shift yields scalable convergence. We acknowledge the absence of a rigorous bound for variable coefficients. In the revision, we will expand §3.1 to include the dispersion analysis for the homogeneous case and discuss the empirical extension to heterogeneous media, while noting that a full theoretical proof remains an open question for future work. revision: partial

  3. Referee: [§5] §5 (Comparison with complex-shifted Laplacian): at 10 ppw the hybrid method is reported to require an order-of-magnitude fewer iterations, yet the baseline complex-shift parameter and its scaling with wavenumber are not tabulated, preventing direct assessment of whether the improvement stems from the real-shift component or from the specific choice of the modest complex shift.

    Authors: We appreciate this observation. The complex-shift parameters for the baseline method were selected following standard practices in the literature (e.g., proportional to the square of the wavenumber), but their specific values were not explicitly listed. In the revised manuscript, we will add a table in §5 detailing the complex-shift magnitudes used in the comparisons, along with their scaling, to facilitate direct assessment of the results. revision: yes

standing simulated objections not resolved
  • A rigorous mathematical proof or bound proving that a constant real shift cancels position-dependent dispersion errors in spatially varying coefficient Helmholtz problems.

Circularity Check

0 steps flagged

No circularity: real-shifted coarse correction is a proposed ansatz validated by experiments

full rationale

The derivation chain consists of proposing a real shift applied only to the coarsest Galerkin operator to correct inter-grid dispersion, then reporting numerical results showing wavenumber-independent 3-level convergence at 11-12 points per wavelength in heterogeneous media. No equations, parameters, or uniqueness theorems are presented that reduce the claimed scalability to a fit, self-definition, or self-citation chain; the central result is an empirical demonstration of the new operator choice rather than a tautological renaming or prediction of its own inputs.

Axiom & Free-Parameter Ledger

1 free parameters · 1 axioms · 0 invented entities

Only the abstract is available, so the ledger is necessarily incomplete. The central claim rests on the unproven premise that a real shift applied solely at the coarsest level corrects dispersion without side effects.

free parameters (1)
  • real shift magnitude
    The amount of real shift applied to the coarsest Galerkin operator is not quantified in the abstract and is presumably chosen or derived to match dispersion error.
axioms (1)
  • domain assumption Real shift on coarsest grid alone corrects numerical dispersion between levels
    Invoked to justify why the method works without complex shifts on finer levels.

pith-pipeline@v0.9.0 · 5482 in / 1311 out tokens · 50317 ms · 2026-05-10T02:02:51.457876+00:00 · methodology

discussion (0)

Sign in with ORCID, Apple, or X to comment. Anyone can read and Pith papers without signing in.

Reference graph

Works this paper leans on

42 extracted references · 42 canonical work pages

  1. [1]

    Aminzadeh, B

    F. Aminzadeh, B. Jean, and T. Kunz , 3-D salt and overthrust models , Society of Exploration Geophysicists, 1997

    work page 1997

  2. [2]

    Azulay and E

    Y. Azulay and E. Treister , Multigrid-augmented deep learning preconditioners for the Helmholtz equation , SIAM Journal on Scientific Computing, (2022)

    work page 2022

  3. [3]

    Bayliss, C

    A. Bayliss, C. I. Goldstein, and E. Turkel , On accuracy conditions for the numerical computation of waves , J. Comput. Phys., 59 (1985), pp. 396--404

    work page 1985

  4. [4]

    Benamou and B

    J.-D. Benamou and B. Despr \`e s , A domain decomposition method for the Helmholtz equation and related optimal control problems , J. Comput. Phys., 136 (1997), pp. 68--82

    work page 1997

  5. [5]

    Berenger , A perfectly matched layer for the absorption of electromagnetic waves , Journal of Computational Physics, 114 (1994), pp

    J.-P. Berenger , A perfectly matched layer for the absorption of electromagnetic waves , Journal of Computational Physics, 114 (1994), pp. 185--200

    work page 1994

  6. [6]

    Bezanson, A

    J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah , Julia : A fresh approach to numerical computing , SIAM Review, 59 (2017), pp. 65--98

    work page 2017

  7. [7]

    Brandt , Multi-level adaptive solutions to boundary-value problems , Mathematics of computation, 31 (1977), pp

    A. Brandt , Multi-level adaptive solutions to boundary-value problems , Mathematics of computation, 31 (1977), pp. 333--390

    work page 1977

  8. [8]

    W. L. Briggs, V. E. Henson, and S. F. McCormick , A multigrid tutorial , SIAM, second ed., 2000

    work page 2000

  9. [9]

    Brougois, M

    A. Brougois, M. Bourget, P. Lailly, M. Poulet, P. Ricarte, and R. Versteeg , Marmousi, model and data , in EAEG workshop-practical aspects of seismic data inversion, European Association of Geoscientists & Engineers, 1990, pp. cp--108

    work page 1990

  10. [10]

    J. Chen, V. Dwarka, and C. Vuik , A matrix-free parallel two-level deflation preconditioner for two-dimensional heterogeneous Helmholtz problems , Journal of Computational Physics, 514 (2024), p. 113264

    work page 2024

  11. [11]

    W. Chen, Y. Liu, and X. Xu , A robust domain decomposition method for the Helmholtz equation with high wave number , ESAIM: Mathematical Modelling and Numerical Analysis, 50 (2016), pp. 921--944

    work page 2016

  12. [12]

    Cocquet and M

    P.-H. Cocquet and M. J. Gander , How large a shift is needed in the shifted Helmholtz preconditioner for its effective inversion by multigrid? , SIAM Journal on Scientific Computing, 39 (2017), pp. A438--A478

    work page 2017

  13. [13]

    Cocquet and M

    P.-H. Cocquet and M. J. Gander , Asymptotic dispersion correction in general finite difference schemes for helmholtz problems , SIAM Journal on Scientific Computing, 46 (2024), pp. A670--A696

    work page 2024

  14. [14]

    Cocquet, M

    P.-H. Cocquet, M. J. Gander, and X. Xiang , A finite difference method with optimized dispersion correction for the helmholtz equation , in International Conference on Domain Decomposition Methods, Springer, 2017, pp. 205--213

    work page 2017

  15. [15]

    Cocquet, M

    P.-H. Cocquet, M. J. Gander, and X. Xiang , Closed form dispersion corrections including a real shifted wavenumber for finite difference discretizations of 2d constant coefficient helmholtz problems , SIAM Journal on Scientific Computing, 43 (2021), pp. A278--A308

    work page 2021

  16. [16]

    C. Cui, K. Jiang, and S. Shu , A neural multigrid solver for helmholtz equations with high wavenumber and heterogeneous media , SIAM Journal on Scientific Computing, 47 (2025), pp. C655--C679

    work page 2025

  17. [17]

    Dwarka and C

    V. Dwarka and C. Vuik , Scalable convergence using two-level deflation preconditioning for the Helmholtz equation , SIAM Journal on Scientific Computing, 42 (2020), pp. A901--A928

    work page 2020

  18. [18]

    H. C. Elman, O. G. Ernst, and D. P. O'leary , A multigrid method enhanced by Krylov subspace iteration for discrete Helmholtz equations , SIAM Journal on scientific computing, 23 (2001), pp. 1291--1315

    work page 2001

  19. [19]

    Engquist and A

    B. Engquist and A. Majda , Absorbing boundary conditions for numerical simulation of waves , Proceedings of the National Academy of Sciences, 74 (1977), pp. 1765--1766

    work page 1977

  20. [20]

    Y. A. Erlangga, C. W. Oosterlee, and C. Vuik , A novel multigrid based preconditioner for heterogeneous Helmholtz problems , SIAM J. Sci. Comput., 27 (2006), pp. 1471--1492

    work page 2006

  21. [21]

    O. G. Ernst and M. J. Gander , Multigrid methods for helmholtz problems: A convergent scheme in 1d using standard components , Direct and Inverse Problems in Wave Propagation and Applications. De Gruyer, (2013), pp. 135--186

    work page 2013

  22. [22]

    R. D. Falgout, M. Lecouvez, P. Ramet, and C. Richefort , Toward an algebraic multigrid method for the indefinite helmholtz equation , SIAM Journal on Scientific Computing, (2025), pp. S285--S310

    work page 2025

  23. [23]

    M. J. Gander, I. G. Graham, and E. A. Spence , Applying GMRES to the Helmholtz equation with shifted Laplacian preconditioning: what is the largest shift for which wavenumber-independent convergence is guaranteed? , Numerische Mathematik, 131 (2015), pp. 567--614

    work page 2015

  24. [24]

    M. J. Gander and H. Zhang , Domain decomposition methods for the Helmholtz equation: a numerical investigation , in Domain Decomposition Methods in Science and Engineering XX, Springer, 2013, pp. 215--222

    work page 2013

  25. [25]

    Haber and S

    E. Haber and S. MacLachlan , A fast method for the solution of the Helmholtz equation , J. Comput. Phys., 230 (2011), pp. 4403--4418

    work page 2011

  26. [26]

    Harari, M

    I. Harari, M. Slavutin, and E. Turkel , Analytical and numerical studies of a finite element PML for the Helmholtz equation , Journal of Computational Acoustics, 8 (2000), pp. 121--137

    work page 2000

  27. [27]

    C.-H. Jo, C. Shin, and J. H. Suh , An optimal 9-point, finite-difference, frequency-space, 2-D scalar wave extrapolator , Geophysics, 61 (1996), pp. 529--537

    work page 1996

  28. [28]

    Livshits , A scalable multigrid method for solving indefinite Helmholtz equations with constant wave numbers , Numerical Linear Algebra with Applications, 21 (2014), pp

    I. Livshits , A scalable multigrid method for solving indefinite Helmholtz equations with constant wave numbers , Numerical Linear Algebra with Applications, 21 (2014), pp. 177--193

    work page 2014

  29. [29]

    L. N. Olson and J. B. Schroder , Smoothed aggregation for Helmholtz problems , Numerical Linear Algebra with Applications, 17 (2010), pp. 361--386

    work page 2010

  30. [30]

    Rabinovich, D

    D. Rabinovich, D. Givoli, and E. B \'e cache , Comparison of high-order absorbing boundary conditions and perfectly matched layers in the frequency domain , International Journal for Numerical Methods in Biomedical Engineering, 26 (2010), pp. 1351--1369

    work page 2010

  31. [31]

    Singer and E

    I. Singer and E. Turkel , High-order finite difference methods for the Helmholtz equation , Computer Methods in Applied Mechanics and Engineering, 163 (1998), pp. 343--358

    work page 1998

  32. [32]

    C. C. Stolk , A rapidly converging domain decomposition method for the Helmholtz equation , J. Comput. Phys., 241 (2013), pp. 240--252

    work page 2013

  33. [33]

    C. C. Stolk , A two-grid method with dispersion matching for finite-element helmholtz problems: C. stolk , Advances in Computational Mathematics, 51 (2025), p. 43

    work page 2025

  34. [34]

    C. C. Stolk, M. Ahmed, and S. K. Bhowmik , A multigrid method for the Helmholtz equation with optimized coarse grid corrections , SIAM Journal on Scientific Computing, 36 (2014), pp. A2819--A2841

    work page 2014

  35. [35]

    Treister and E

    E. Treister and E. Haber , A multigrid solver to the helmholtz equation with a point source based on travel time and amplitude , Numerical Linear Algebra with Applications, 26 (2019), p. e2206

    work page 2019

  36. [36]

    Trottenberg, C

    U. Trottenberg, C. Oosterlee, and A. Sch\" u ller , Multigrid , Academic Press, London and San Diego, 2001

    work page 2001

  37. [37]

    Tsuji and R

    P. Tsuji and R. Tuminaro , Augmented AMG -shifted Laplacian preconditioners for indefinite Helmholtz problems , Numerical Linear Algebra with Applications, 22 (2015), pp. 1077--1101

    work page 2015

  38. [38]

    Turkel, D

    E. Turkel, D. Gordon, R. Gordon, and S. Tsynkov , Compact 2D and 3D sixth order schemes for the Helmholtz equation with variable wave number , J. Comput. Phys., 232 (2013), pp. 272--287

    work page 2013

  39. [39]

    Umetani, S

    N. Umetani, S. P. MacLachlan, and C. W. Oosterlee , A multigrid-based shifted Laplacian preconditioner for a fourth-order Helmholtz discretization , Numerical Linear Algebra with Applications, 16 (2009), pp. 603--626

    work page 2009

  40. [40]

    Yovel, Y

    R. Yovel, Y. He, and E. Treister , Vanka-smoothed shifted laplacian multigrid preconditioners for the helmholtz equations , arXiv preprint arXiv:2511.16808, (2025)

    work page arXiv 2025

  41. [41]

    Yovel and E

    R. Yovel and E. Treister , A block-acoustic preconditioner for the elastic helmholtz equation , arXiv preprint arXiv:2411.15897, (2024)

    work page arXiv 2024

  42. [42]

    Yovel and E

    R. Yovel and E. Treister , LFA -tuned matrix-free multigrid method for the elastic Helmholtz equation , SIAM Journal on Scientific Computing, (2024), pp. S1--S21

    work page 2024