Discretisation: the exact scheme¶
This page derives the shared interior discretisation and gives its exact coefficients: the time-stepping scheme, the spatial finite-difference stencils, the stability (CFL) limit and where it comes from, the grid sampling needed to control dispersion, and the source/receiver and resampling conventions. Everything below is the scheme Deepwave actually runs; the coefficients are transcribed from src/deepwave/regular_grid.h and src/deepwave/staggered_grid.h, and the stability numbers are re-derived from those same coefficients in the appendix.
Time stepping¶
The non-attenuating interior wave equations use an explicit, second-order-accurate leapfrog (centred-in-time) update. Constitutive attenuation and absorbing-boundary memory updates are described on the Viscoacoustic (attenuating) wave equation and Boundary conditions pages. Time is the only dimension whose accuracy order is fixed; the spatial order is the user-selectable accuracy.
Regular-grid scalar propagator¶
The scalar and scalar_born propagators solve
where \(u\) is the modelled field, \(v\) the wavespeed, and \(f\) the source term. On a single regular grid the leapfrog update is (ignoring the PML memory variables, which add absorbing terms but not time levels)
with \(\nabla^2_h\) the discrete Laplacian below and \(f\) the public source density defined in Units and source conventions. The homogeneous update is the following line in scalar.c; the separately injected source has already been scaled by \(-v^2\Delta t^2\) in scalar.py:
wfp_t[i] = v_shot[i] * v_shot[i] * dt2 * w_sum + 2 * wfc_t[i] - wfp_t[i];
where w_sum is \((\nabla^2_h u^n + \text{PML terms})\), wfc is \(u^n\), and wfp holds \(u^{n-1}\) on entry and \(u^{n+1}\) on exit (the two-time-level ring buffer, see below).
Staggered-grid acoustic and elastic propagators¶
The acoustic and elastic propagators use the Virieux (1986) velocity–stress staggered grid (extended to higher spatial order in the manner of Levander, 1988): the pressure/stress and the particle-velocity components live on grids offset by half a cell in space and are advanced on time levels offset by half a step. Schematically, for the acoustic system (pressure \(p\), particle velocity \(\mathbf w\), buoyancy \(b=1/\rho\), modulus \(\kappa=\rho v^2\)),
where \(\nabla_h^{\pm}\) are the staggered first-derivative operators below (forward/backward half-cell). The elastic propagator advances the velocity–stress system of the same family. Because velocity and stress are each updated from the other, this is still a single second-order leapfrog in time; the half-step time stagger is the standard velocity–stress arrangement, not an extra time level.
Two-time-level state and quiescent initial conditions¶
The scalar recurrence needs two field levels, \(u^n\) and \(u^{n-1}\) (wavefield_0 and wavefield_m1). The staggered first-order systems instead retain one level of each component: stress at \(n\Delta t\) and velocity at \((n-1/2)\Delta t\). Absorbing-boundary and constitutive memory is carried in auxiliary_state. Together these form the complete continuation state. Resuming with initial_state preserves the recurrence when the grid, internal step and coefficients are the same. Whole-record FFT resampling must also be preserved; independently resampling source segments need not reproduce one full call (see Forward operators and their derivatives).
When you do not supply initial wavefields they default to zero. Zero in both stored levels is a consistent quiescent past: it represents a medium that was exactly at rest for all \(t<0\) (\(u\equiv 0\), \(\partial_t u\equiv 0\)), not merely a single zeroed snapshot. The propagator’s boundary bookkeeping relies on this — the C sources document that “the first and last accuracy/2 elements in each spatial dimension are zero in all wavefields” and that the auxiliary PML fields start zero — so the default run begins from a physically consistent, source-driven initial-value problem.
Spatial finite-difference stencils¶
Deepwave uses central Taylor (maximally accurate) finite-difference coefficients — the classical weights of Fornberg (1988). The order is set by accuracy \(M\in\{2,4,6,8\}\), and the radius of the stencil (the number of guard cells, FD_PAD) is \(M/2\). All coefficients below are exact rationals as they appear in the source (stored as double literals such as 4.0/12.0).
Regular-grid second derivative¶
The scalar Laplacian uses the central second-derivative stencil
order | \(c_0\) | \(c_1\) | \(c_2\) | \(c_3\) | \(c_4\) |
|---|---|---|---|---|---|
2 | \(-2\) | \(1\) | |||
4 | \(-\tfrac{5}{2}\) | \(\tfrac{4}{3}\) | \(-\tfrac{1}{12}\) | ||
6 | \(-\tfrac{49}{18}\) | \(\tfrac{3}{2}\) | \(-\tfrac{3}{20}\) | \(\tfrac{1}{90}\) | |
8 | \(-\tfrac{205}{72}\) | \(\tfrac{8}{5}\) | \(-\tfrac{1}{5}\) | \(\tfrac{8}{315}\) | \(-\tfrac{1}{560}\) |
The first-derivative terms that appear only inside the PML (DIFFZ1 etc.) use the matching central first-derivative weights \(\{\tfrac12\}\), \(\{\tfrac{2}{3},-\tfrac{1}{12}\}\), \(\{\tfrac34,-\tfrac{3}{20},\tfrac{1}{60}\}\), \(\{\tfrac45,-\tfrac15,\tfrac{4}{105},-\tfrac{1}{280}\}\) for orders 2/4/6/8.
Staggered-grid first derivative¶
The acoustic/elastic propagators evaluate first derivatives between grid points (half-cell staggered). The derivative at the point half a cell above \(i\) is
i.e. \(w_m\) weights the difference of the pair of samples \((2m-1)\) half-cells either side of the derivative point.
order | \(w_1\) | \(w_2\) | \(w_3\) | \(w_4\) |
|---|---|---|---|---|
2 | \(1\) | |||
4 | \(\tfrac{9}{8}\) | \(-\tfrac{1}{24}\) | ||
6 | \(\tfrac{75}{64}\) | \(-\tfrac{25}{384}\) | \(\tfrac{3}{640}\) | |
8 | \(\tfrac{1225}{1024}\) | \(-\tfrac{245}{3072}\) | \(\tfrac{49}{5120}\) | \(-\tfrac{5}{7168}\) |
Stability: the CFL condition¶
An explicit leapfrog is only stable if the time step is small enough. A von Neumann analysis of the leapfrog wave update (substitute a plane wave \(u^n_{\mathbf i}=g^n e^{i\mathbf k\cdot\mathbf x_{\mathbf i}}\) and require \(|g|\le1\)) gives the classical Courant–Friedrichs–Lewy limit
where the dimensionless critical Courant number \(C_{\text{crit}}\) depends only on the spatial stencil (order and regular/staggered scheme). It is re-derived exactly from the stencil coefficients in the appendix; the results are
order | regular grid | staggered grid |
|---|---|---|
2 | \(1.000\) | \(1.000\) |
4 | \(0.866\) | \(0.857\) |
6 | \(0.813\) | \(0.805\) |
8 | \(0.784\) | \(0.777\) |
What Deepwave actually uses. By default Deepwave does not run at these limits. deepwave.discretisation.cfl_condition_n() picks the internal time step from a single conservative constant Courant number of 0.6, the same for every accuracy order and every number of dimensions:
If the requested dt is larger, the propagator subdivides it into the smallest integer number of internal steps (step_ratio) that satisfies this bound. One constant covers all dimensionalities because the number of dimensions enters the CFL limit only through the \(\sqrt{\sum_i 1/\Delta x_i^2}\) term — the critical Courant number itself is a property of the 1-D stencil, unchanged by adding dimensions. Using 0.6 (rather than, say, 0.78–1.0) provides a margin below the homogeneous stencil limit. A pointwise wavespeed alone does not bound arbitrary heterogeneous velocity/stress operators: a high-order derivative can couple a stiff cell to a distant velocity node with much lower density.
The CFL calculation uses dimensionless ratios, so equivalent length and time units give the same subdivision without squaring very large or very small speeds. A zero speed or zero dt needs no subdivision. If the required internal time step cannot be represented, the helper raises an error.
The staggered propagators therefore include a density-weighted operator bound in the conservative speed reported as RunInfo.max_vel. Acoustic and viscoacoustic use the weighted derivative stencil; elastic families bound the constitutive matrix and staggered buoyancy. Spatially constant elastic materials retain their directional Christoffel speed bound. This can increase automatic subdivision at sharp material contrasts, even if the pointwise physical speed is constant. An explicit max_vel must cover this bound as well as the physical speed. Use deepwave.plan() to inspect it before a large run. The bound concerns the interior energy pairing; physical boundary closures and PML require their separate validation.
The margin between 0.6 and the true \(C_{\text{crit}}\) is real headroom. Because the number of internal steps scales with \(1/C\), a low-order run that is safe at \(C=0.9\) instead of \(0.6\) uses \(0.6/0.9\approx0.67\times\) as many steps — a ~1.5x speed-up (up to ~1.7x at \(C_{\text{crit}}=1.0\)). You can opt in to a larger, still-stable Courant number via the courant_number argument of the propagators (see Accuracy and grid design); it is validated against the \(C_{\text{crit}}\) table above and errors if you request an unstable value. The default (None) is equivalent to explicitly selecting 0.6. The heterogeneous operator correction can change subdivision relative to earlier releases.
Grid sampling and numerical dispersion¶
Even when stable, a finite-difference wave travels at a slightly wrong, frequency-dependent speed (numerical dispersion). The spatial part depends on how many grid points sample a wavelength. For a plane wave sampled at \(G\) points per wavelength (so \(\theta = k\,\Delta x = 2\pi/G\)), the numerical phase velocity of the regular-grid scheme is
(axial propagation), with the analogous expression \(\tilde k\,\Delta x = 2\sum_m w_m\sin\!\bigl((2m-1)\theta/2\bigr)\) for the staggered first-derivative scheme.
Criterion. For each order we report the smallest points-per-wavelength \(G\) at which the axial phase-velocity error \(|v_{\text{num}}/v - 1|\) first falls to \(\le 1\%\) and to \(\le 0.1\%\). This is a local error; over a propagation distance \(D\) the accumulated phase error grows linearly, \(\Delta\phi = 2\pi\,(D/\lambda)\,|v_{\text{num}}/v-1|\), so a 1% velocity error slips a full cycle after ~100 wavelengths of propagation. Use the 1% column for short offsets and the 0.1% column for long models.
order | regular, 1% | regular, 0.1% | staggered, 1% | staggered, 0.1% |
|---|---|---|---|---|
2 | 12.8 | 40.6 | 12.8 | 40.6 |
4 | 5.3 | 9.6 | 5.1 | 9.2 |
6 | 3.9 | 6.0 | 3.8 | 5.8 |
8 | 3.4 | 4.8 | 3.3 | 4.6 |
Deepwave’s built-in “at least 6 points per wavelength” warning (see Accuracy and grid design) is calibrated for the default accuracy=4: at 6 ppw the order-4 axial phase error is ~0.6%, order-6 ~0.1%, order-8 ~0.02%, but order-2 is still ~4.5% — a 2nd-order run needs a much finer grid (or a higher order) for the same fidelity.
Letting Deepwave do the arithmetic. deepwave.suggest_discretisation() puts the table above and the CFL condition together: give it the slowest wavespeed in your model and the highest frequency you need to propagate, and it returns a cell spacing, a time step, and a sentence saying which constraint was binding:
>>> import deepwave
>>> suggestion = deepwave.suggest_discretisation(1500, 25, 2, max_vel=4500)
>>> suggestion.grid_spacing, round(suggestion.dt, 6)
((10.0, 10.0), 0.000943)
It applies both criteria and takes whichever asks for more cells: the dispersion table above, and the six-cells-per-wavelength floor of Deepwave’s own runtime check. Those two disagree — the 1% column asks for only 5.3, 3.9 and 3.4 cells at orders 4, 6 and 8 — and the floor is what binds there, because the plane-wave dispersion criterion is a statement about a single Fourier component on an infinite uniform grid, while the runtime check is a blunter guard on everything that analysis omits (how a point source, a sharp contrast and the PML are sampled). A grid the suggestion returns therefore never provokes Deepwave’s own warning. The price is that at the default 1% target orders 4, 6 and 8 all land on the same six-cell grid; pass min_points_per_wavelength=0 to follow the dispersion table alone, and expect the warning when you then call a propagator.
The table is axial, and axial is the conservative direction (see grid anisotropy immediately below), so a diagonally travelling wave is slightly better resolved than the suggestion promises.
Grid anisotropy. Numerical dispersion is not isotropic: a wave travelling along a grid axis and one travelling diagonally are dispersed by different amounts, because the diagonal direction samples the stencil differently. For the regular 2-D order-2 scheme at 5 ppw, for instance, the axial phase velocity is \(0.936\,v\) while the diagonal is \(0.967\,v\) — the diagonal is faster (less dispersed). The anisotropy shrinks with more points per wavelength and with higher order, but never vanishes; keep it in mind when comparing arrivals at different azimuths.
Temporal dispersion is separate. The \(O(\Delta t^2)\) time-stepping error is independent of the spatial accuracy order. Raising accuracy reduces spatial dispersion but does nothing to the temporal dispersion — that is controlled only by dt (hence by the Courant number). A dedicated time-dispersion-correction utility, deepwave.time_dispersion (a forward/inverse Koene-Robertsson transform of the wavelet and data), is available for callers who want to reduce temporal dispersion when using a larger dt — see Time-dispersion correction and Time-dispersion correction.
Source and receiver convention¶
A point source in Deepwave is injected as a cell delta: the source amplitude at time step \(n\) is added to the single grid cell containing the source, with no \(1/(\Delta x_1\cdots\Delta x_n)\) cell-volume normalisation. Concretely, in the scalar propagator the amplitude is pre-scaled by \(-v^2\,\Delta t^2\) (the -v**2 * dt**2 factor applied in scalar.py before stepping) and added to one cell — it is not divided by the cell area or volume.
For a point source of strength \(S(t)\), the source density supplied at a cell is \(S(t)/\Delta V\), where \(\Delta V\) is the product of spatial grid spacings. Equivalently, a fixed supplied amplitude represents a strength proportional to \(\Delta V\). In reduced dimensions that is strength per unit omitted length or area. See Units and source conventions for the physical quantity and sign of each channel.
The transpose relation is between geometrical scatter and gather on the same component grid. Adding an unscaled amplitude to a cell is the transpose of reading that cell under the ordinary tensor dot product. Physical source injection additionally includes channel-dependent time-step and material factors; different physical source and receiver channels are not automatically transposes of one another. The full derivative reverses every operation in the forward composition, including those factors and resampling; see Forward operators and their derivatives and What Deepwave guarantees about gradients.
Temporal resampling¶
When the internal step is finer than the user’s dt (step_ratio > 1), Deepwave uses Fourier interpolation of a finite periodic record, with a band-limiting projection. It is a periodic sinc interpolation, not an assumption that the physical record outside the interval is known. deepwave.discretisation.upsample() takes the real FFT of the source wavelet, zero-pads the spectrum up to the finer Nyquist (zeroing the original Nyquist bin so no energy is created above the input band), and inverse transforms; deepwave.discretisation.downsample() truncates the spectrum to the coarser Nyquist (zeroing that bin for even output lengths) before the inverse transform. Spectral padding and truncation retain only existing Fourier bins; optional cosine frequency tapering and time padding are available to suppress ringing and wrap-around. A positive taper fraction does not change the Nyquist-bin rule. Downsampling requires the input time length to be divisible by step_ratio; pad or trim explicitly when a different output length is required. Odd output lengths are supported and have no Nyquist bin to remove. Time-domain windows and finite-record padding can broaden spectra; the bandlimit statement concerns the Fourier interpolation step itself. Its normalisation and transpose are derived in Forward operators and their derivatives.
Appendix: re-deriving the critical Courant numbers¶
The \(C_{\text{crit}}\) table above is derived here directly from the stencil coefficients, so nothing is taken on faith.
Regular grid¶
Insert a plane wave into the leapfrog update. The discrete second-derivative operator has Fourier symbol (per dimension, in units of \(1/\Delta x^2\))
which is most negative at \(\theta=k\Delta x=\pi\) (the grid Nyquist). Writing \(S_M \equiv -\widehat{D^2}(\pi) = -\bigl(c_0 + 2\sum_m c_m(-1)^m\bigr)\) for its maximum magnitude, the leapfrog amplification factor stays on the unit circle (\(|g|=1\)) provided
Evaluating \(S_M\) from the second-derivative weights:
order | \(S_M\) | \(C_{\text{crit}}=2/\sqrt{S_M}\) |
|---|---|---|
2 | \(4\) | \(1.0000\) |
4 | \(16/3\) | \(\sqrt3/2 = 0.8660\) |
6 | \(272/45\) | \(0.8135\) |
8 | \(2048/315\) | \(0.7844\) |
Worked example, order 4: \(S_M = -\bigl(-\tfrac52 + 2(-\tfrac43 - \tfrac1{12})\bigr) = \tfrac{16}{3}\), so \(C_{\text{crit}} = 2/\sqrt{16/3} = \sqrt3/2 \approx 0.866\).
Staggered grid¶
The staggered first-derivative operator has the purely imaginary symbol \(i\,\tilde k(\theta)\Delta x\) with \(\tilde k(\theta)\Delta x = 2\sum_m w_m\sin\!\bigl((2m-1)\theta/2\bigr)\). Its magnitude is maximised at \(\theta=\pi\), where every sine term reaches \(\pm1\) and aligns with the weight’s sign, giving
The velocity–stress leapfrog composes two such first derivatives, so its stability condition is the same form as above with \(\sqrt{S_M}\) replaced by \(K_M\):
order | \(K_M = 2\sum_m|w_m|\) | \(C_{\text{crit}}=2/K_M\) |
|---|---|---|
2 | \(2\) | \(1.0000\) |
4 | \(7/3\) | \(6/7 = 0.8571\) |
6 | \(149/60\) | \(0.8054\) |
8 | \(2161/840\) | \(0.7774\) |
Worked example, order 4: \(K_M = 2(\tfrac98 + \tfrac1{24}) = \tfrac73\), so \(C_{\text{crit}} = 2/(7/3) = 6/7 \approx 0.857\).
Both derived tables match the values baked into deepwave.discretisation.critical_courant_number and quoted at the top of this page.
References¶
J. Virieux (1986), “P-SV wave propagation in heterogeneous media: Velocity- stress finite-difference method”, Geophysics 51(4), 889–901.
A. R. Levander (1988), “Fourth-order finite-difference P-SV seismograms”, Geophysics 53(11), 1425–1436.
B. Fornberg (1988), “Generation of finite difference formulas on arbitrarily spaced grids”, Mathematics of Computation 51(184), 699–706.