Time-dispersion correction

This page explains temporal dispersion – the error a leapfrog time step introduces independently of the spatial grid – and how to remove it after the fact approximately with the optional deepwave.time_dispersion module. It is for anyone pushing dt close to the stability limit to save compute, and wanting the accuracy of a smaller dt without paying its cost.

Every Deepwave propagator advances the wavefield with an explicit, second-order-accurate leapfrog step (see Discretisation: the exact scheme). That scheme is only second order in time, so it introduces temporal dispersion: a frequency-dependent error in the wave speed that, unlike spatial dispersion, is not reduced by raising the accuracy order. It is controlled solely by the time step. Pushing dt toward the CFL stability limit to save compute therefore trades speed for temporal-dispersion error — unless you correct for it.

The deepwave.time_dispersion module implements the correction of Koene, Robertsson, Broggini & Andersson (“Eliminating time dispersion from seismic wave modeling”, Geophysical Journal International 213, 2018). It is opt-in, pure PyTorch, differentiable, and model-independent.

The current transforms evaluate finite-record spectra with dense frequency by time matrices. Their temporary memory grows quadratically with record length: a 10,000-sample record uses roughly 50 million complex entries per matrix. Use them for moderate records, and budget memory accordingly.

The dispersion relation

Substituting a single Fourier mode \(u_n = e^{i\omega n\Delta t}\) into the leapfrog second difference shows that a mode oscillating in the simulation at temporal angular frequency \(\omega\) carries the physics of the lower frequency

\[\Omega(\omega) = \frac{2}{\Delta t}\,\sin\!\Bigl(\frac{\omega\,\Delta t}{2}\Bigr) \;\le\; \omega .\]

Equivalently, physics that should oscillate at \(\Omega\) is shifted up in the simulation to

\[\omega(\Omega) = \frac{2}{\Delta t}\,\arcsin\!\Bigl(\frac{\Omega\,\Delta t}{2}\Bigr) \;\ge\; \Omega , \qquad \Omega < \frac{2}{\Delta t}.\]

These two maps are exact inverses on the simulation band \(|\omega| < \pi/\Delta t\), which maps to the physical band \(|\Omega| < 2/\Delta t\).

The scalar identity for a single oscillator is that the leapfrog recursion,

\[u_{n+1} - 2u_n + u_{n-1} = -\Delta t^2\,\Omega_0^2\,u_n ,\]

is solved exactly by \(u_n = \sin(\omega_{\mathrm{num}}\,n\,\Delta t)\) with \(\omega_{\mathrm{num}} = (2/\Delta t)\,\arcsin(\Omega_0\,\Delta t/2)\). This identity describes the scalar oscillator only. It does not establish exact correction of finite signals or of other propagator families.

The two transforms

The correction is a pair of frequency-domain remappings applied outside the propagator:

  • correct_source() pre-distorts the source wavelet, applying the inverse map so that content the caller places at physical frequency \(\Omega\) is injected at the simulation frequency \(\omega(\Omega)\). After the simulation shifts it back down, it lands where intended.

  • correct_receivers() post-processes the recorded data, moving content from the simulation frequency \(\omega\) back down to the physical frequency \(\Omega(\omega)\).

The maps describe the scalar leapfrog frequency relation. Their use on finite records is approximate: only the retained frequency band is remapped, and content outside that band is lost. Finite-record spectral leakage, padding and windowing also affect the result. Exact evaluation of the finite-record DTFT does not make a remap round trip exact. This scalar relation does not establish exact correction of attenuating, staggered-grid or PML dynamics, or remove spatial discretisation error.

Which time step to pass

The dispersion is governed by the leapfrog step the propagator actually takes — its internal time step internal_dt = dt / step_ratio — not by the user-facing dt. Deepwave subdivides dt into step_ratio internal steps to satisfy the CFL condition. Both transforms therefore take internal_dt explicitly, and treat the passed tensors’ time axis as sampled at that step.

The simplest recipe is to pick a dt that is already CFL-stable, so the propagator runs with step_ratio == 1 and internal_dt == dt; then correct the source and the receivers at dt. cfl_internal_dt() returns the exact internal_dt the propagator would use for a given grid spacing, dt, maximum velocity and courant_number. Pass the grid spacing as the sequence of per-dimension spacings you give the propagator ([dx, dx]), or as a single value together with ndim (cfl_internal_dt(dx, dt, max_vel, ndim=2)), since the CFL limit depends on the number of dimensions. When step_ratio > 1, resample your source/receiver tensors to internal_dt (with deepwave.discretisation.upsample() / deepwave.discretisation.downsample()) before correcting, or force step_ratio == 1.

See Time-dispersion correction for a runnable end-to-end example.

API

correct_source Function

deepwave.time_dispersion.correct_source(source_amplitudes: Tensor, internal_dt: float, nt_out: int | None = None) → Tensor

Pre-correct a source wavelet for leapfrog time dispersion.

Applies the Koene-Robertsson source pre-transform: it pre-distorts the wavelet so that, after the second-order leapfrog time stepping shifts every temporal frequency up to omega(Omega) = (2/dt) arcsin(Omega dt/2), the modelled physics follows the scalar leapfrog frequency map. Feed the returned wavelet to the propagator in place of the original; pair it with correct_receivers() on the recorded data to reduce temporal dispersion.

The transform is the inverse frequency map: the corrected spectrum at simulation angular frequency \(\omega\) equals the input spectrum at the lower physical frequency \(\Omega(\omega) = (2/dt)\sin(\omega dt/2)\). Because the pre-warp pushes energy toward higher frequencies, the corrected wavelet is more oscillatory and may need a longer record to avoid time wrap-around; extend it with nt_out.

internal_dt must be the propagator’s internal (leapfrog) time step dt / step_ratio, and source_amplitudes are treated as sampled at that step. See the module docstring and cfl_internal_dt() for how to obtain it; the simplest usage is a CFL-stable dt giving step_ratio == 1 so internal_dt == dt.

Parameters:

  • source_amplitudes – Real tensor of shape [..., nt] (for example the usual [shot, source, nt]) holding the source wavelet(s). The time axis is the last dimension. Any number of leading batch dimensions is supported.

  • internal_dt – The internal (leapfrog) time step governing the temporal dispersion, in seconds. Must be positive.

  • nt_out – Optional length of the returned time axis, in time samples. Defaults to the input nt. Use a larger value to give the (broader-band) corrected wavelet room and avoid wrap-around.

Returns:

The pre-corrected source amplitudes, shape [..., nt_out], with the same dtype and device as source_amplitudes.

Raises:

  • TypeError – If source_amplitudes is not a tensor, internal_dt is not a float, or nt_out is not an int.

  • ValueError – If source_amplitudes is not a floating tensor with at least one dimension, internal_dt is not positive, or nt_out is not positive.

correct_receivers Function

deepwave.time_dispersion.correct_receivers(receiver_amplitudes: Tensor, internal_dt: float, nt_out: int | None = None) → Tensor

Post-correct recorded data for leapfrog time dispersion.

Applies the Koene-Robertsson receiver post-transform: it moves recorded spectral content from the simulation frequency \(\omega\) back down to the physical frequency \(\Omega(\omega) = (2/dt)\sin(\omega dt/2)\) at which it should sit, reducing temporal dispersion described by the scalar leapfrog frequency map. Apply it to the receiver output of the propagator; pair it with correct_source() on the input wavelet.

Equivalently, the corrected trace’s spectrum at physical angular frequency \(\Omega\) is read from the recorded spectrum at the higher simulation frequency \(\omega(\Omega) = (2/dt)\arcsin(\Omega dt/2)\). Physical frequencies above \(2/dt\) cannot be represented by the leapfrog and are set to zero.

internal_dt must be the propagator’s internal (leapfrog) time step dt / step_ratio, and receiver_amplitudes are treated as sampled at that step. See the module docstring and cfl_internal_dt(); the simplest usage is a CFL-stable dt giving step_ratio == 1 so internal_dt == dt.

Parameters:

  • receiver_amplitudes – Real tensor of shape [..., nt] (for example the usual [shot, receiver, nt]) holding the recorded traces. The time axis is the last dimension. Any number of leading batch dimensions is supported.

  • internal_dt – The internal (leapfrog) time step governing the temporal dispersion, in seconds. Must be positive.

  • nt_out – Optional length of the returned time axis, in time samples. Defaults to the input nt. The corrected traces are narrower band (energy moves to lower frequency), so the default length is normally sufficient.

Returns:

The dispersion-corrected receiver amplitudes, shape [..., nt_out], with the same dtype and device as receiver_amplitudes.

Raises:

  • TypeError – If receiver_amplitudes is not a tensor, internal_dt is not a float, or nt_out is not an int.

  • ValueError – If receiver_amplitudes is not a floating tensor with at least one dimension, internal_dt is not positive, or nt_out is not positive.

cfl_internal_dt Function

deepwave.time_dispersion.cfl_internal_dt(grid_spacing: float | Sequence[float], dt: float, max_vel: float, courant_number: float | None = None, ndim: int | None = None) → float

Return the internal (leapfrog) time step a propagator would use.

Convenience wrapper around deepwave.discretisation.cfl_condition_n() that returns just the internal time step internal_dt = dt / step_ratio – the dt that governs temporal dispersion and that correct_source() / correct_receivers() require. Passing the same grid_spacing, dt, max_vel (and optional courant_number) you give the propagator guarantees you correct with exactly the step the propagator runs.

The CFL limit depends on the number of spatial dimensions (each dimension contributes to the bound), and a propagator applies a single grid_spacing value to every dimension of its model. A single value therefore does not on its own say what the propagator would do, so ndim must be given with it; a sequence of per-dimension spacings already carries that information and needs no ndim.

Parameters:

  • grid_spacing – The grid cell spacing in each spatial dimension (the same value or sequence passed to the propagator). If a single value (applied to every dimension, as the propagators do), ndim must also be given. Deepwave is unit-agnostic: use whatever consistent system of units grid_spacing, dt, and max_vel are expressed in (see Units and source conventions).

  • dt – The user-facing time step, in the time units of that same consistent system.

  • max_vel – The maximum absolute wavespeed in the model, in grid_spacing units per dt unit; the propagator derives step_ratio from this via the CFL condition.

  • courant_number – Optional Courant number, matching the propagator’s courant_number argument. If None (default), the conservative value the propagator uses by default (0.6) is applied, reproducing the historical internal step exactly.

  • ndim – The number of spatial dimensions of the propagation. Required when grid_spacing is a single value; optional (but checked against the sequence’s length) otherwise.

Returns:

The internal time step internal_dt, in the time units of dt.

Raises:

  • TypeError – If ndim is not an int.

  • ValueError – If grid_spacing is a single value and ndim was not given, or if ndim is not positive.

  • RuntimeError – If grid_spacing is a sequence whose length is not ndim.