Arbitrary high order shaped stencils for time domain finite difference schemes in seismic wave propagation

Technical Summary: Arbitrary High Order Shaped Stencils for Time Domain Finite Difference Schemes in Seismic Wave Propagation

Problem Statement
Seismic imaging relies heavily on the numerical solution of wave equations, where high-fidelity solutions are required to suppress spurious numerical waves and accurately preserve physical dispersion relations. While Finite Difference (FD) schemes are the primary choice for seismic imaging due to their computational efficiency, they suffer from inherent numerical dispersion. This error causes the phase velocity to depend on grid spacing and discretization direction (anisotropy), often manifesting as a slowing of high-frequency components. Although classical methods based on Taylor series expansions are dominant, they frequently yield inadequate representations of wave dispersion. The paper addresses the question of whether effective gains in computational savings and better representation of realistic velocity models can be achieved through optimized FD stencils of arbitrary shapes and orders, moving beyond the standard cross-stencil geometry.

Methodology
The authors propose a general mathematical framework for deriving high-order time-domain FD schemes for the acoustic wave equation with arbitrary stencil shapes. The study focuses on the isotropic acoustic case, implemented within the Devito symbolic computational framework, which allows for the generation of highly optimized, stencil-based code.

The methodology encompasses three primary components:

1. Stencil Geometries: The paper investigates various stencil shapes, including:

   • Cross Stencils: Classical line stencils along the axes.
   • Rhombus and Cross-Rhombus Stencils: Adding "inner" points off the axes.
   • Square and Cross-Square Stencils: Filling the stencil with a square-like geometry.
     The study defines the geometric parameters ($M$, $N$, $P$, $Q$) that control the number of grid points and the symmetry of the weights.

2. Derivation Approaches: The authors unify and extend several derivation strategies to determine the stencil weights:

   • SpatTE: Spatial Taylor Expansions (classical approach).
   • SpecTE: Spectral Taylor Expansions of the spatial operator.
   • SpecLS: Least Squares approximation of the Spectral Taylor series.
   • DispTE: Dispersion-based Taylor truncation, trading temporal derivatives for spatial derivatives to reduce dispersion errors.
   • DispLS: Least Squares approximation of the full time-space spectral operator to minimize dispersion error.

3. Evaluation Metrics: Beyond classical $L_1$, $L_2$, and $L_\infty$ error norms, the authors employ Discrete Fourier Transform (DFT) and Continuous Wavelet Transform (CWT) analyses. These tools quantify dispersion effects by measuring weighted phase differences and amplitude deviations between numerical and reference solutions across frequency and time-frequency domains.

Key Contributions

• Unified Framework: The paper provides a general mathematical formulation for deriving arbitrary-order FD schemes with arbitrary stencil shapes, unifying existing literature on Taylor, Spectral, and Least Squares approaches.
• Systematic Comparison: It conducts a comprehensive numerical comparison of different stencil geometries (Cross, Rhombus, Square) combined with various derivation methods (SpatTE, SpecTE, DispTE, SpecLS, DispLS).
• Implementation in Devito: The methods are implemented in Devito, demonstrating how high-level symbolic definitions can be compiled into efficient, low-level code suitable for full-scale acoustic seismic inversion problems.
• Analysis of Non-Cross Stencils: The study rigorously tests the hypothesis that adding non-cross (inner) points to stencils improves accuracy.

Results
Numerical experiments were conducted on four velocity models: Homogeneous, Two-Layer, 2D SEG/EAGE Salt, and Marmousi.

• Homogeneous and Simple Heterogeneous Models: In constant or simple two-layer velocity models, optimized schemes (particularly DispTE-CrossRb and DispLS-CrossRb) generally outperform classical SpatTE schemes in reducing dispersion, especially at higher orders ($M \geq 4$).
• Complex Heterogeneous Models (Salt/Marmousi): In models with significant velocity contrasts, the performance of velocity-dependent optimized schemes (like DispTE and DispLS) degrades when a single stencil is applied across the entire grid based on the maximum velocity.
• Geometry Findings:
  • Cross Stencils vs. Complex Shapes: Contrary to the intuition that more points yield better accuracy, the results indicate that Cross stencils generally provide the best trade-off between accuracy and computational cost.
  • Diminishing Returns: Adding non-cross points (Rhombus or Square shapes) does not necessarily provide added accuracy or dispersion reduction. In many cases, particularly with Square stencils, the increased number of weights leads to ill-posed linear systems, numerical instabilities, or negligible improvements compared to the increased computational cost.
  • Robustness of Spatial-Only Methods: In highly heterogeneous media, schemes that do not rely on velocity-dependent weights for optimization (specifically SpecLS-CrossRb and SpatTE-Cross) demonstrated greater robustness and lower dispersion artifacts than velocity-dependent schemes.
• Trade-offs: For high-order schemes ($M \geq 4$) in homogeneous regions, optimized Cross stencils offer significant dispersion reduction. However, in complex media, the sensitivity of optimized schemes to velocity variations suggests that simpler, velocity-independent schemes or adaptive strategies may be preferable.

Significance and Claims
The paper claims that while the literature often explores complex stencil geometries to improve accuracy, non-cross stencil shapes (rhombus and square) do not necessarily provide added accuracy or dispersion reduction in general, despite their increased computational cost.

The primary significance of the work lies in:

1. Validating Cross Stencils: Confirming that dispersion-optimized Cross stencils offer adequate accuracy with reduced computational cost compared to classic approaches and more complex shapes.
2. Highlighting Velocity Sensitivity: Demonstrating that the benefits of dispersion-optimized schemes are highly dependent on the velocity model. In complex, heterogeneous media, schemes relying on velocity-dependent weight computation can underperform compared to spatially optimized or classical schemes.
3. Practical Implementation: Providing a ready-to-use implementation in Devito, facilitating the application of these methods to full-scale acoustic seismic inversion (FWI) and Reverse Time Migration (RTM) problems.
4. Future Adaptation: Noting that the framework can be adapted to anisotropic models (VTI/TTI) and that the pre-computation of stencils (as in SpecLS) offers efficiency advantages for iterative inversion workflows.

The authors conclude that for practical applications in complex media, the choice of temporal and spatial parameters, as well as the velocity model, critically influences solution quality, and that the "Cross" geometry remains the most efficient and robust choice for the tested scenarios.