PaperlayerView exact revision on arXiv

Source-linked AI summary

A dispersion minimizing scheme for the 3-D Helmholtz equation based on ray theory

Christiaan C. Stolk

arXiv:1504.01609v2math.NAcs.NAphysics.comp-ph

TL;DR

Long-distance Helmholtz simulations face phase errors that grow with source distance, while coarse meshes are attractive for reducing computational demands. The paper develops a ray-theory-based compact finite-difference scheme and reports accurate solutions on coarse meshes, with still coarser discretizations usable in multigrid. Its 3-D tests also report favorable computation times compared with some literature results.

  • Problem

    For Helmholtz discretizations, phase errors grow with distance from the source, and solver cost is also important, particularly for coarse-mesh applications.

    Source: §1 · p. 1

  • Method

    The paper develops a compact dispersion-minimizing finite-difference scheme in 2-D and 3-D, using ray theory for Helmholtz-like difference equations and amplitude correction operators.

    Source: §1 · p. 2

  • Results

    Five to six points per wavelength yield accurate solutions in constant and smoothly varying media, while three points per wavelength can be used for multigrid coarse levels; 3-D computation times compare favorably with some literature results.

    Source: §1 · p. 3

  • Takeaways & Limitations

    The scheme is presented for both direct coarse-mesh simulation and coarse-level use in multigrid solvers.

    Source: §1 · p. 2

    The analysis assumes the discrete symbol is Helmholtz-like, including a specified zero-set parameterization and nonzero gradient on its zero set.

    Source: §2.2 · p. 5

Abstract

from arXiv · show

We develop a new dispersion minimizing compact finite difference scheme for the Helmholtz equation in 2 and 3 dimensions. The scheme is based on a newly developed ray theory for difference equations. A discrete Helmholtz operator and a discrete operator to be applied to the source and the wavefields are constructed. Their coefficients are piecewise polynomial functions of $hk$, chosen such that phase and amplitude errors are minimal. The phase errors of the scheme are very small, approximately as small as those of the 2-D quasi-stabilized FEM method and substantially smaller than those of alternatives in 3-D, assuming the same number of gridpoints per wavelength is used. In numerical experiments, accurate solutions are obtained in constant and smoothly varying media using meshes with only five to six points per wavelength and wave propagation over hundreds of wavelengths. When used as a coarse level discretization in a multigrid method the scheme can even be used with downto three points per wavelength. Tests on 3-D examples with up to $10^8$ degrees of freedom show that with a recently developed hybrid solver, the use of coarser meshes can lead to corresponding savings in computation time, resulting in good simulation times compared to the literature.

1. Introduction

The paper introduces a dispersion-minimizing compact finite-difference scheme for 2-D and 3-D Helmholtz problems, justified by ray theory for difference equations. It targets accurate solutions on coarse meshes and use as a coarse-level discretization in multigrid solvers.

1. Introduction

  • Dispersion errors cause phase errors to grow with source distance, motivating coarse-mesh schemes that control phase slowness while keeping operators inexpensive to apply or invert.Compact finite-difference methods use small 3 × 3 or 3 × 3 × 3 stencils; QS-FEM has the smallest phase slowness errors among the cited schemes but is available only in 2-D.Source: §1 · pp. 1–2
  • The new second-order scheme targets phase slowness errors comparable to 2-D QS-FEM, with amplitude correction operators and a ray-theory justification.The paper develops ray theory for Helmholtz-like difference equations to approximate the phase and amplitude functions of the continuous operator.Source: §1 · p. 2
  • The method extends earlier parameter tuning by making coefficients depend smoothly on hk through third-order Hermite interpolation, further reducing phase slowness errors.Earlier schemes used chosen parameters or piecewise constant or linear coefficient dependence [26].Source: §1 · p. 2
  • Five to six points per wavelength yield accurate solutions in constant and smoothly varying media, while three points per wavelength can be used for multigrid coarse levels.At three points per wavelength, the solutions may be too inaccurate for direct use but can be refined iteratively in multigrid.Source: §1 · p. 3
  • In 3-D tests with the SEG-EAGE Salt Model, cost per degree of freedom changes little as frequency increases, and computation times compare favorably with some literature results.The tests use a two-grid solver with an inexact coarse-level inverse and IOFD as the coarse-level discretization.Source: §1 · p. 3

2. Theory of discrete Helmholtz equations with constant k

The section develops a Fourier-symbol framework for analyzing constant-coefficient discrete Helmholtz operators, their phase slowness, and their Green’s functions. It then gives conditions under which modified discrete solutions match the continuous solution’s asymptotic phase and amplitude.

2. Theory of discrete Helmholtz equations with constant k

  • The analysis uses Fourier transforms to define the operator’s symbol and phase slownesses, then studies the discrete Green’s function and its outgoing asymptotics.A dimensionally reduced symbol is also introduced in the framework.Source: §2 · p. 4, §2.3 · p. 6

2.1. Symbol and phase slownesses.

  • Phase slowness error measures the relative discrepancy between discrete and exact wave vectors, and accumulated phase error grows in proportion to propagation distance measured in wavelengths.For a constant-coefficient difference operator, propagating plane waves correspond to zeros of its symbol.Source: §2.1 · pp. 4–5

2.2. Dimensional reduction and Helmholtz-like symbols.

  • A Helmholtz-like symbol is characterized by geometric and regularity conditions on its zero set, including a smooth normal map; the section shows that this property follows when P1 is Helmholtz-like.The definition also constrains the zero set to the frequency domain and requires P(0) to be negative.Source: §2.2 · p. 5

2.3. The discrete Green’s function.

  • The discrete Green’s function has non-unique distributional solutions, and its outgoing solution is selected by a condition on the Fourier-domain term B.Its spatial asymptotics are obtained by stationary phase; for the usual 3-D Helmholtz operator, the leading asymptotic expansion equals the outgoing Green’s function.Source: §2.3 · pp. 6, 8, 9

2.4. The modified discrete Green’s function.

  • Matching the zero sets of the discrete and continuous symbols matches phase, while a condition on the source and wavefield operators also matches asymptotic amplitude.The section states these two conditions together as sufficient for asymptotic agreement with the continuous Helmholtz solution.Source: §2.4 · p. 9

3. Theory of discrete Helmholtz equations with variable k

For variable k, the paper defines discrete operators and symbols and derives ray-theoretic equations for phase and amplitude. In symbol form, these equations match the continuous equations, and an additional operator condition is used to align discrete and continuous amplitudes.

3.2. Ray-theoretic equations for amplitude and phase.

  • The discrete ray theory yields the same eikonal and transport equations as the continuous theory when both are expressed using operator symbols.The phase satisfies a nonlinear first-order equation, while the amplitude satisfies a transport equation.Source: §3 · p. 10, §3.2 · p. 13

3.1. Symbols and operators for variable k.

  • For variable k, the operator coefficients depend on k at points determined by a fixed parameter t, and t = 1/2 gives a symmetric discretization when the coefficient functions satisfy the stated symmetry.The symbol’s quantization type corresponds to t; t = 1/2 is Weyl quantization.Source: §3 · p. 10, §3.1 · pp. 10–11

3.2. Ray-theoretic equations for amplitude and phase.

  • At t = 1/2, the amplitude transport equation conserves |A|^2.The ray-theoretic construction uses the high-frequency limit with ωh held constant and smoothly varying c(x).Source: §3.2 · pp. 11, 13

3.3. Amplitude correction.

  • For amplitude correction, matching the phase zero sets and imposing the stated operator conditions makes the discrete ray phase coincide with the reference phase.The construction uses symmetric discretization and equal source and wavefield correction operators; the associated transport operators are scaled versions of one another.Source: §3 · p. 10, §3.3 · p. 14

4. Phase slowness errors for existing discretizations

This section describes standard and compact finite differences and Lagrange finite elements on regular meshes, then compares their phase slowness errors across propagation directions and grid resolutions. The reported comparison identifies QS-FEM as the smallest-error method among those discussed in 2-D, while noting it is unavailable in 3-D.

4. Phase slowness errors for existing discretizations

  • The comparison includes standard finite differences, compact schemes, and regular-mesh Lagrange finite elements, with phase slowness errors evaluated across propagation directions and grid resolutions.The maximum relative error is taken over θ and plotted as a function of inverse points per wavelength, 1/G.Source: §4 · p. 15, §4.4 · pp. 18–19

4.1. Standard finite differences.

  • Standard finite differences approximate each one-dimensional second derivative with a central difference and discretize −k^2u pointwise.The section mentions central-difference orders 2, 4, 6, and 8.Source: §4.1 · p. 15

4.2. Compact finite difference discretizations.

  • Existing compact finite difference schemes use parameter choices that minimize maximum phase slowness error, while other approaches vary coefficients piecewise to reduce errors.For [26], the objective is instead to minimize differences from a fine-scale operator, so its phase errors are not directly comparable.Source: §4.2 · p. 17

4.3. Lagrange finite elements on regular meshes.

  • The Lagrange finite element description builds regular-mesh trial and test functions from one-dimensional Lagrange polynomials and their tensor products.Its discretization is derived from the weak form.Source: §4.2 · p. 17, §4.3 · p. 18

4.4. Phase slowness errors.

  • The paper plots phase slowness errors for some 2-D schemes in Figure 1 and some 3-D schemes in Figure 2, each against the inverse number of points per wavelength, 1/G.Source: Figure 1 · p. 19, Figure 2 · p. 19

5. A dispersion minimizing scheme with amplitude corrections

The section introduces a compact Helmholtz discretization whose operators are chosen to minimize phase and amplitude errors using ray theory for difference equations. It reports roughly 20-fold lower phase errors than compact sixth-order schemes and accurate long-range results at five points per wavelength for IOFD and QS-FEM.

5. A dispersion minimizing scheme with amplitude corrections

  • The IOFD scheme solves a discrete Helmholtz system with compact operators P and Q selected to minimize phase and amplitude discrepancies in ray-theoretic solutions.The design seeks to match the phase and amplitude behavior of the continuous Helmholtz operator.Source: §5 · pp. 19–20

5.1. IOFD discretization of the Helmholtz operator.

  • The Helmholtz-operator coefficients α_j are piecewise-polynomial functions of 1/G, specified by control-node values and first derivatives, and fitted to reduce integrated phase slowness error.The objective also includes a regularization term involving derivatives of α_j.Source: §5.1 · p. 20

5.3. Comparison of phase slowness errors.

  • Roughly 20-fold lower phase errors are reported for IOFD and QS-FEM than for compact sixth-order schemes SUT and CHO6.The comparison includes IOFD in two and three dimensions and QS-FEM in two dimensions.Source: Figure 3 · p. 22, §5.2 · p. 22, §5.3 · p. 23

5.2. Amplitude correction operators.

  • The amplitude-correction operator Q uses coefficients β_j defined by Hermite interpolation and chosen through a linearly constrained least-squares problem.The control values are chosen to minimize a discrete approximation of an integral objective.Source: §5.2 · p. 22

5.3. Comparison of phase slowness errors.

  • At five points per wavelength over hundreds of wavelengths, only QS-FEM and IOFD meet the stated δph ≲ 0.01λL criterion.At six points per wavelength, CHO6 is near the bound and FE8 qualifies; at eight points per wavelength, FE N=8 has relative phase slowness errors roughly equal to QS-FEM and IOFD.Source: §5.3 · p. 23

6. Numerical examples

The numerical examples test IOFD on constant and smoothly varying media, including wave propagation over hundreds of wavelengths. They report small phase shifts at coarse resolutions and low amplitude errors in a smoothed Marmousi model.

6. Numerical examples

  • A 45-degree annular segment makes phase shifts visible by comparing the numerical wavefield with the exact solution’s real-part minimum line.Cubic spline interpolation maps Cartesian-grid values to polar coordinates for the plot.Source: §6 · p. 24, Figure 4 · p. 24
  • At 500 wavelengths, IOFD’s maximum phase shifts were 0.0065 radians at 6 ppw, 0.020 at 5 ppw, and 0.089 at 4 ppw, versus CHO6’s 0.27, 0.84, and π radians.At 3 ppw the IOFD phase shift was clearly visible; at 2.5 ppw it became large, with the field plotted at 100 wavelengths.Source: §6 · p. 25, Figure 5 · p. 25
  • For IOFD, amplitude errors decreased as resolution rose above 4 ppw; the section also presents amplitude-error plots at 3 and 4 ppw.The plotted errors are relative to the exact solution.Source: §6 · p. 25, Figure 6 · p. 26
  • In the smoothed Marmousi model, relative amplitude errors were below 0.01 over most of the domain, with isolated 0.05 or 0.08 errors where absolute amplitude was small.The examples used IOFD with at least six points per wavelength and compared it with a fourth-order finite-element solution using twice as many grid points per direction.Source: §6 · pp. 25–26

7. Application in multigrid based solvers

The multigrid discussion motivates using IOFD at both fine and coarse levels because prior coarse-grid optimization matched the phase slowness of a less accurate fine-grid operator. Earlier work reported good convergence with optimized finite differences down to three points per wavelength.

7. Application in multigrid based solvers

  • The section situates this approach among two-grid methods with inexact coarse solvers, including a shifted-Laplacian-preconditioned Krylov solver and a double-sweep domain-decomposition preconditioner.These methods are discussed in and [25], respectively.Source: §7 · p. 26
  • Optimized finite differences at the coarse level gave good convergence down to three points per wavelength, whereas standard coarse discretizations needed about 10.The prior results were reported in [26].Source: §7 · p. 27
  • Using IOFD at both levels avoids designing the coarse operator to match the relatively large phase-slowness errors of standard second-order finite differences.In [26], those standard differences were used at the fine level, so the coarse-level method was tuned to match them rather than the true phase slowness.Source: §7 · p. 27

Q. The operator P is used directly as coarse level discretization and at the fine level we are only interested in solving equation (111).

The section tests IOFD in multigrid, including a 3-D SEG-EAGE salt-model example with an inexact coarse solver. Results report good convergence at three points per wavelength and compare computation time with prior work, while noting limits at coarser or higher-frequency settings.

Q. The operator P is used directly as coarse level discretization and at the fine level we are only interested in solving equation (111).

  • For the two-grid method, IOFD behaved well down to three ppw; at 2.5 ppw it still converged, but iterations rose substantially and became more problem-size sensitive.The iteration study covered constant, Marmousi, and a slice of the 3-D SEG-EAGE salt model.Source: Q. The operator P is used directly as… · p. 28, Table 5 · p. 29
  • At three points per wavelength, IOFD is intended as a coarse-grid discretization; direct fine-level use generally produces errors that are too large.The coarse-level use can also support other fine-level discretizations, provided they use a regular mesh.Source: Q. The operator P is used directly as… · p. 28
  • At higher frequencies, convergence cannot be assumed from prior tests, since many solvers were tested at 10 or more mesh points per wavelength and multigrid convergence depends on frequency.The frequency dependence of the double-sweep domain-decomposition method was described as unclear.Source: Q. The operator P is used directly as… · p. 29
  • The inexact-coarse-solver test modifies the two-grid method of [25] to use IOFD at both levels and tests the SEG-EAGE Salt Model at a minimum of six ppw.The example increases frequency by a factor of 5/3; the solver uses a double-sweep domain-decomposition preconditioner for its inexact coarse inverse.Source: Q. The operator P is used directly as… · p. 29
  • Using coarser meshes can yield corresponding computation-time savings, although the 27-point compact operator costs more to apply than second-order finite differences.The Table 6 discussion reports little cost increase compared with [25] despite the frequency increase by a factor of 5/3.Source: Q. The operator P is used directly as… · p. 29
  • At 9.91 Hz, the SEG-EAGE salt-model solve took 45 seconds per right-hand side on 128 cores, versus 270 seconds at 10 Hz on 256 cores in.The reported residual reduction was 10^6 for this solve and 10^5 in.Source: Table 6 · p. 30, Q. The operator P is used directly as… · p. 30

8. Discussion

For smooth k, IOFD supports coarse meshes of five or six points per wavelength, with smaller phase errors than earlier methods for long-distance propagation. Strong gradients and inversion use impose practical constraints, and a comparison with method is incomplete.

8. Discussion

  • Method uses less memory and has been applied to larger examples, but no full accuracy comparison was made.The authors therefore do not provide a complete accuracy comparison with that method.Source: §8 · p. 30
  • IOFD in 2-D and 3-D, and QS-FEM in 2-D, can make phase errors much smaller than methods and for coarse-mesh, long-distance simulations.The discussion identifies five or six points per wavelength as a minimum for smooth k.Source: §8 · pp. 30–31
  • Where k(x) has strong gradients, locally finer meshes or different discretizations may be needed for accurate solutions.The discussion notes that multigrid may still be useful in this setting.Source: §8 · p. 31
  • For inversion, IOFD and QS-FEM make derivatives of the finite-difference operator with respect to medium coefficients more complicated to compute than methods and.Hermite interpolation makes these derivatives continuous for IOFD.Source: §8 · p. 31

Appendix A. Phase slowness computations for finite element methods

For periodically repeating discretizations, the phase-slowness computation uses an operator matrix and selects zero eigenvalues, with a plane-wave association requiring additional assumptions. The numerical selection method uses the eigenvector most correlated with the constant function.

Appendix A. Phase slowness computations for finite element methods

  • For periodically repeating operators, the symmetry property applies only when p, q, and r are divisible by N, and the operator acts as an N^3 × N^3 matrix on periodic components.The matrix is used to find ξ values with a zero eigenvalue.Source: Appendix A · p. 32
  • A zero eigenvalue does not automatically identify a plane wave because the periodic component can represent a linear combination of wave vectors.The associated wave vector can be determined modulo N when the zero eigenvalue is simple and its eigenfunction is close to a plane wave.Source: Appendix A · p. 32
  • In the computations, the selected eigenvalue is the one whose eigenvector has the largest absolute inner product with the constant function.A phase-velocity vector is identified when this eigenvalue is zero.Source: Appendix A · pp. 32–33
  • The selection approach has limitations, but appeared sufficient for standard finite elements when the mesh had more than four points per wavelength.Source: Appendix A · p. 33
Loading 1504.01609v2…