Superstep wavefield propagation

TL;DR

This work tackles the inefficiency of traditional one-step FD wavefield propagation by introducing a two-stage superstep: precompute propagator matrices (PMs) per grid location to encode physics, then apply these PMs to advance $k$ time steps at once. For first-order systems, time evolution is written as $(oldsymbol{u}_{n+k},oldsymbol{v}_{n+k})^T = S(k)igl(oldsymbol{u}_n,oldsymbol{v}_nigr)^T$ with $S(k)=S_{11}(k)S_{12}(k)S_{21}(k)S_{22}(k)$, while second-order systems use a companion propagator $G_k$ that satisfies a recurrence $G_{k+1}=G_1G_k-G_{k-1}$; in both cases $k$-step evolution is captured by these precomputed operators. PMs, assembled per grid location, allow a separation between physics and compute, enabling domain-specific language (DSL) mappings and scalable cloud-based implementations. Synthetic Marmousi tests validate that PMs can encode boundary reflections and that superstep wavefields reproduce traditional results with high fidelity, demonstrating a practical pathway to accelerate large-scale linear wave propagation, gradient computation for inversion, and migration. Overall, the framework shifts the computational bottleneck from physics integration to memory-aware data management and hardware mapping, with broad implications for high-performance full-waveform inversion and seismic imaging.

Abstract

This paper describes how to propagate wavefields for arbitrary numbers of traditional time steps in a single step, called a superstep. We show how to construct operators that accomplish this task for finite-difference time domain schemes, including temporal first-order schemes in isotropic, anisotropic and elastic media, as well as temporal second-order schemes for acoustic media. This task is achieved by implementing a computational tradeoff differing from traditional single step wavefield propagators by precomputing propagator matrices for each model location for k timesteps (a superstep) and using these propagator matrices to advance the wavefield k time steps at once. This tradeoff separates the physics of the propagator matrix computation from the computer science of wavefield propagation and allows each discipline to provide their optimal modular solutions.

Figures (8)

• Figure 1: Sketch illustrating the propagating matrices (PM) in their physical domain. 2 PMs $G(A)$ and $G(B)$ shown corresponding to 30 timesteps of propagation from locations $A$ and $B$. The compact support for $G(A)$ and $G(B)$ are shown as gray squares. The wavefield $u_{60}$ at time step 60 is indicated in light gray shade and it is computed for the whole physical domain.
• Figure 2: Illustration of the dot product calculation between a propagating matrix and the corresponding wavefield. Top left: PM after 30 time steps. Top right: Wavefield after 300 timesteps in the physical 2D domain. Bottom left: Wavefield segment corresponding to the physical area of the PM. Bottom right: Element-wise multiplication of the PM and the wavefield segment.
• Figure 3: Comparison between the numerical and physical domains of influence. The same 30-step propagator matrix is plotted with "Accent" color scale in (a) and with gray in (b). The yellow diamond in (a) is the computational domain while the formed ring of wave in (b) is the physical domain.
• Figure 4: A propagator matrix near a reflecting boundary (a) without compensation for the boundary layer; (b) compensated for the boundary layer and (c) their difference.
• Figure 5: Wavefield after 420 traditional time steps. Top left: Wavefield after 420 steps calculated by traditional method; Top right: Wavefield after 420 steps calculated by superstep method. Initial wavefield at 30 steps was advanced by 30 steps at once using propagator matrices; Bottom left: Homogeneous velocity model overlaid by propagator matrices at proportional size and the source location (white dot); Bottom right: Difference between the 1-step and superstep methods at the time scale of the wave field.