What Deepwave Calculates

Deepwave provides forward modelling and backpropagation of wave equations using the finite-difference method within the PyTorch framework. This page explains, at a conceptual level, how sources and receivers work, how internal time stepping and gradients relate to the dt and model you supply, and what Deepwave guarantees about the gradients it returns; it is for readers who already have a propagation running (see Quick Start) and want to understand what is happening underneath it, or who need the exact gradient guarantees for validation.

Sources

Source amplitudes are scalar functions of time. In scalar propagators, sources are located at grid cell centres. In the elastic propagator, sources are positioned at the \(v_y\) and \(v_x\) component locations. Creating an effective source with an arbitrary direction, or at a fractional grid cell location, may require distributing the source across an array of multiple sources over neighbouring cells. The Hicks method can assist with this.

Receivers

Receiver data is generated by recording the relevant wavefield at specified receiver locations. The recorded data spans the same time range and has the same time sampling as the input source data. Similar to sources, approximating directed receivers or those at fractional grid cell locations can be achieved by combining data recorded from multiple neighbouring grid cells. The Hicks method is also useful here.

Internal time stepping and the CFL condition

To satisfy the Courant-Friedrichs-Lewy (CFL) condition, Deepwave may internally use a smaller time step than that of the input source and output receiver data. This internal time step is always a divisor of the source/receiver time sampling interval (calculated using cfl_condition). Source amplitudes are internally upsampled to this finer rate, and conversely, receiver data is downsampled before being returned to the user (via upsample and downsample).

When the internal step is finer than your dt (step_ratio > 1), compiled Deepwave kernels normally sample the model-parameter gradient at the coarser input/output rate. This saves wavefield-history memory but approximates the full discrete adjoint sum. The effective stride is step_ratio * model_gradient_sampling_interval. Source and initial-state gradients still follow every internal step. See What Deepwave guarantees about gradients for the other conditions and the backend distinction.

Wavefield outputs

Deepwave propagators return the final wavefields, one opaque auxiliary_state (the absorbing boundary’s memory variables and any other auxiliary state, see Usage), and the recorded receiver data. These wavefields represent the complete wavefields used internally by Deepwave, including padding for the Perfectly Matched Layer (PML). Because of this, the whole output can be passed directly as initial_state to another propagator instance, to continue a propagation.

The sections below state Deepwave’s gradient guarantees precisely: exactly which gradients are exact, which are approximations, and under what conditions each holds. This is reference material for validating gradients rigorously (for example, against finite differences) or for understanding an edge case. If you do not need that level of detail, the per-propagator pages linked at the bottom of this page are a better next stop.

What Deepwave guarantees about gradients

Deepwave’s backpropagation is the discrete adjoint of the finite-difference forward recurrence, not a separately discretised continuous adjoint. This section states precisely what that does and does not guarantee. (For the forward scheme these gradients are adjoint to, see Discretisation: the exact scheme.)

Source and initial-wavefield gradients follow every internal step. The gradients of the outputs with respect to the source amplitudes and initial wavefields (wavefield_0/wavefield_m1 and the absorbing boundary’s memory variables, carried opaquely in auxiliary_state) are the exact discrete adjoint (transpose) of the forward recurrence, with the numerical setup held fixed. The symbolic tests in tests/test_adjoint_symbolic.py independently differentiate manually transcribed interior updates for selected scalar and elastic-family terms. They do not read the kernel sources or cover every family, PML recurrence or boundary operation. Whole-operator transpose tests in tests/test_adjoint_dot_product.py check the \(\langle Fs, r\rangle = \langle s, F^{*}r\rangle\) identity for their selected configurations. These are complementary checks, rather than a symbolic proof of all compiled kernels. The practical meaning: for these inputs the returned gradient differentiates the fixed discrete loss, subject to floating-point roundoff. Its relative numerical error also depends on conditioning; a machine-precision relative bound is not promised. It is not a discretised continuous adjoint that happens to be close; it is the exact derivative of the discrete operator when discretisation and inferred parameters are fixed. A torch.autograd.gradcheck.gradcheck() comparison also has finite-difference truncation and roundoff error; agreement to machine precision is not a general expectation.

Why sampling at the data interval is useful. A fixed, linear, time-independent medium preserves the temporal frequencies of its forced response, including under spatial differentiation. If both correlation factors are band-limited below \(1/(2\,dt)\), sinc-basis orthogonality gives

\[\int_{-\infty}^{\infty} a(t)\overline{b(t)}\,dt = dt\sum_{n=-\infty}^{\infty} a(n\,dt)\overline{b(n\,dt)}.\]

The product has twice the individual bandwidth, but this identity needs its integral (the zero-frequency component), not reconstruction of the entire product. Thus the ideal gradient integral can be obtained at the data interval; it does not require every CFL substep. With internal interval \(h=dt/r\), storing every internal step would require roughly \(r\) times the wavefield-history storage and history traffic. Fixed model and working-wavefield allocations do not grow by that factor. This is why the sampled model-gradient default is retained.

Both correlation factors must meet the bandwidth condition, including the adjoint forced through the transpose of receiver resampling. Velocity and stress live at staggered times; their correlation uses values at the proper half-step phase, rather than pretending they share the same timestamp. The finite-record DFT alias identity and its exact error term are given in Forward operators and their derivatives. At a fixed duration the memory comparison is M_fixed + r*M_history versus M_fixed + M_history for the default sampled history. Both choices still advance the live fields at every CFL step.

These assumptions concern the actual forward and adjoint fields. Sampling source and receiver arrays at dt alone does not prove strict band limitation: finite-record endpoints, initial/terminal transients and arbitrary initial wavefields or wavefield objectives can violate it. Further sampling with model_gradient_sampling_interval > 1 requires a correspondingly lower bandwidth. For the sampling basis, see Prandoni and Vetterli, Signal Processing for Communications; Nyquist-rate wavefield storage is also used in time-subsampled FWI.

The honest caveat: the model gradient at step_ratio > 1. The gradient with respect to the model parameters (e.g. the wavespeed v) is a time integral of the forward and adjoint wavefields. To save compute and memory, Deepwave samples this integral only every step_ratio * model_gradient_sampling_interval internal steps and multiplies by that stride — a strided rectangle-rule quadrature of the exact adjoint integral (with the default model_gradient_sampling_interval = 1 the stride is exactly step_ratio). In the C source this is visible as the forward pass storing the wavefield only when (t % stride) == 0 and the backward pass accumulating wfc_t[i] * w_store_1_t[i] * stride. Consequently:

• At step_ratio = 1 and model_gradient_sampling_interval = 1 the compiled model-gradient sum uses every internal step. It is then the derivative of the fixed discrete recurrence, subject to the conditions below and lossless storage.

• At step_ratio > 1 the sparse sum need not equal the finite discrete adjoint sum. Error depends on the correlation bandwidth and record endpoints; it need not decrease monotonically as the stride changes. The quadrature stride the kernel actually uses is step_ratio * model_gradient_sampling_interval, so exact differentiation of the fixed discrete recurrence is guaranteed by using both equal to 1, subject to the conditions below. Supply a dt at or below the CFL limit (so step_ratio = 1) and leave model_gradient_sampling_interval = 1. Increasing model_gradient_sampling_interval coarsens the quadrature further to save memory, with potentially greater error. It must divide nt whenever a model requires a gradient: the current history layout requires complete sampling blocks. Observation callbacks have their own independent schedule. Deepwave raises a ValueError naming the nearest valid nt values for an incomplete block. When no model requires a gradient there is nothing to sample and the interval is ignored.

If you are validating model gradients against finite differences, do so at step_ratio = 1; tests/test_scalar.py includes an absolute finite-difference gradcheck that measures quadrature error for selected step_ratio values and its disappearance at 1.

Finite-record sampling measurements

tests/test_gradient_sampling_characterisation.py compares the default compiled model gradient with full eager differentiation on fixed 16-by-16 float64 models, with 5 m spacing, 1 ms user sampling, internal stride two, fixed max_vel=3000 and pml_freq=25, and a tapered receiver-energy objective. All nine families have nonzero sensitivities. The entries below are relative gradient-norm errors, in percent, measured on the Linux CPU candidate. They describe this experiment rather than a universal error bound.

Finite-record quadrature error (%)

Family	25 Hz, 64 ms	25 Hz, 128 ms	180 Hz, 128 ms
scalar	0.00412	0.0000191	0.0000582
scalar_born	0.000778	0.0000246	0.0000373
acoustic	0.000496	0.00000451	0.0000482
acoustic_born	0.0105	0.0000149	0.0000570
elastic	0.0277	0.0000954	0.000564
elastic_born	0.00620	0.0000594	0.000346
elastic_vti	0.0303	0.000104	0.000679
acoustic_vti	0.000255	0.00000701	0.000541
viscoacoustic	0.00120	0.0000101	0.000592

The longer low-frequency record reduces endpoint error here. At 128 ms the higher-frequency cases have greater sampling error; their short wavelengths also make them deliberately underresolved spatially. Sampling error is measured against the same discrete spatial experiment, so this comparison does not establish physical waveform accuracy at 180 Hz. The scaled directional error abs(dot(g_sampled-g_full,d))/(norm(g_full)*norm(d)) remains below 5.4e-6 across these 36 cases. The independent DFT test includes an aliasing counterexample, and the initial-wavefield test in test_acoustic_vti.py checks sampled weight and phase against separate per-step differentiated leaves. Arbitrary endpoint or wavefield objectives can produce larger errors.

Further gradient conditions

The maximum-wave-speed caveat. The PML absorption profile’s amplitude and the internal time step are both built from max_vel. When max_vel is given they are constants of the call. When it is omitted, the selected speed bound may depend on the model through several physical parameters, especially for anisotropic or attenuating materials. Tied or switching maxima are nonsmooth. The returned derivative holds this resolved bound fixed; it does not follow its dependence on those model entries. Fixing max_vel removes that setup dependence, while exact differentiation of the fixed discrete recurrence also requires effective gradient stride one and lossless storage. The same fixed-setup condition applies to Born and forward-mode derivatives.

out.info.step_ratio is an integer function of the same maximum, so it too is constant in the model except at the isolated models where it steps. pml_freq, when it is not given, is inferred from the peak of the source amplitudes’ spectrum, so the source amplitudes then act on the profile as well; gradients with respect to them do not follow that path either. For a model-gradient check, fix max_vel and pml_freq, use model_gradient_sampling_interval=1, choose dt so out.info.step_ratio == 1, and use lossless storage (storage_compression=False). Keep accuracy, grid spacing, PML width and boundary fixed across perturbations. This verifies the derivative of that particular discrete call. The eager autograd backends have their own capability restrictions; inspect Capability matrix before extrapolating this compiled-kernel recipe to them.

If your scientific data has a coarser sample interval, the following scalar recipe makes a finer external call and preserves the default source and receiver resampling around it:

import torch
import deepwave
from deepwave.discretisation import cfl_condition_n, upsample, downsample

dx, coarse_dt, max_vel, pml_freq = 4.0, 0.004, 1500.0, 25.0
v = torch.full((32, 32), 1500.0, dtype=torch.float64).requires_grad_()
src = deepwave.wavelets.ricker(
    25.0, 24, coarse_dt, 0.04, dtype=torch.float64
).reshape(1, 1, -1)
inner_dt, ratio = cfl_condition_n([dx, dx], coarse_dt, max_vel)
fine_src = upsample(src, ratio)
out = deepwave.scalar(
    v, dx, inner_dt,
    source_amplitudes=fine_src,
    source_locations=torch.tensor([[[16, 16]]]),
    receiver_locations=torch.tensor([[[17, 17]]]),
    nt=fine_src.shape[-1], accuracy=4, pml_width=4,
    max_vel=max_vel, pml_freq=pml_freq,
    model_gradient_sampling_interval=1,
    storage_compression=False,
)
assert out.info.step_ratio == 1
predicted_at_coarse_rate = downsample(out.receiver_amplitudes, ratio)
predicted_at_coarse_rate.square().sum().backward()
assert v.grad is not None

For nondefault freq_taper_frac, time_pad_frac or time_taper, pass the same values to both resampling functions. Other propagators may use a receiver time shift; preserve it when downsampling. Keep the numerical and PML settings fixed across perturbations. This checks the derivative of the finer composite call, which is a different call from the original coarse-dt call. It does not validate the latter’s approximate model gradient.

tests/test_gradient_exactness.py measures all of this: a directional Taylor test along a direction that leaves the maximum alone, and cell-by-cell comparisons against finite differences at the fastest cell, a corner, an edge, a source cell and a receiver cell.

Determinism and reproducibility. On the compiled C/CUDA backend (backend="compiled"), Deepwave uses no CUDA atomics anywhere (there is no atomicAdd in any .cu kernel); gradient accumulation on the GPU is done with a fixed, data-independent traversal. The same care is taken in the PyTorch operations around the kernels: every propagator pads its models out to the PML before stepping, and that padding is built so that its backward sums each edge cell’s contributions in a fixed order rather than scattering them back with atomic adds. For an identical build of libdeepwave_C.so, an identical call configuration, and the same device, repeated runs on the compiled backend give bitwise-reproducible GPU gradients. This guarantee is narrow: it is not a claim that two different devices, two different builds (different CUDA toolkit/driver versions or GPU architectures), or the eager/torch_compile Python backends agree with each other, or with the compiled backend, bitwise. Floating-point summation order, kernel scheduling, and, for torch.compile in particular, graph capture and recompilation can each change a result at the level of the last bit, so an apparent mismatch across builds, devices or backends is not necessarily a bug. Where exact bitwise agreement matters, verify it directly for your own build, configuration and device rather than assuming it, and use float64 for the tightest comparison against finite differences.

On the CPU, every supported propagator has a compiled C/CUDA backend and parallelises a single shot’s interior stencil across threads (intra-shot spatial OpenMP): each internal time step opens an OpenMP parallel region that splits the spatial grid disjointly across threads, while the outer loop over shots is sequential, not parallelised. Because each grid cell’s model gradient is written by exactly one thread per time step — there is no per-thread scratch buffer and no cross-thread summation into the same address — the sequence of updates to any given cell is the same no matter how many threads happen to be available.

That disjoint split is necessary but not sufficient: the thread count also decides where the compiler’s vectorised loop body ends and its scalar remainder begins, so a build allowed to reassociate floating-point arithmetic can evaluate the same cell’s stencil in two different orders depending on how the grid was divided. The kernels are therefore built with reassociation disabled (the loops are elementwise rather than reductions, so this costs little). That removes the main cause, but the flag does not by itself force the vectorised and scalar paths to perform the same operations in the same order, so the outcome is checked rather than assumed: tests/test_determinism.py compares every compiled propagator’s forward outputs and gradients across a range of thread counts. Compiled-backend CPU results — forward outputs and gradients alike — are thread-count invariant: for a given build and input they are bitwise identical regardless of torch.get_num_threads(). This is a tested property of these kernel bodies rather than a guarantee the compiler flag provides, so if you build them yourself with a different compiler or different flags, run those tests.

PML damping is normalised per face. Each face’s physical width sets its own damping amplitude, with continuous normal-incidence design target \(R=0.001\). This target is not a guarantee of measured reflection after discretisation. Unlike the former global-longest-width rule, widening one face does not weaken another. See Boundary conditions for the formula and limits.

Gradient checkpointing and structured Jacobians. torch.utils.checkpoint composes correctly with the propagators after the issue #120 fix: the scalar propagator’s single-backward gradient is verified identical to the non-checkpointed gradient under both use_reentrant=True and use_reentrant=False (see tests/test_scalar.py::test_checkpoint_gradient). Double backprop through a checkpoint requires use_reentrant=False (reentrant checkpoint does not support double backward). The recommended, first-class memory lever is the built-in storage_mode ("device" / "cpu" / "disk"); for long propagations, segmenting the run and wrapping each segment in torch.utils.checkpoint.checkpoint (reentrant) is also available – see Gradient checkpointing in Usage. For structured Jacobians on the fast (compiled) backend, forward-mode jvp / Born gives Jacobian columns and reverse-mode vjp / backward gives Jacobian rows; vmap over rows is currently eager-only (backend="eager"). See the torch.func example (Using torch.func features) for the recipes.

CPU threading and performance

Deepwave does not choose its own CPU thread count: every propagator uses torch.get_num_threads() to decide how many OpenMP threads to launch for a call, so use torch.set_num_threads() to control it. Because CPU parallelism is spatial (within a single shot’s stencil, as described above, rather than a loop over shots), even a single-shot problem can make use of every thread you make available — it is not necessary to batch many shots together to benefit from extra cores.

Benchmark your own workload and hardware, but the general shape to expect when scaling the thread count up on one machine is: speedup rises with the thread count while threads still map to distinct physical cores; once the thread count exceeds the number of physical cores and additional threads start sharing cores (e.g. via hyperthreading/SMT), speedup typically plateaus or dips before partially recovering at higher logical-core counts. Requesting far more threads than physical cores rarely helps and can measurably hurt, so setting the thread count to (at most) your machine’s physical core count is a reasonable starting point if you have not benchmarked otherwise.

If Deepwave shares the machine with other CPU-heavy processes (rather than being the only workload running), consider setting the environment variable OMP_WAIT_POLICY=passive before your process starts. This is a standard OpenMP runtime environment variable (see, for example, GNU libgomp’s OMP_WAIT_POLICY documentation): with an unspecified policy, an OpenMP runtime commonly waits actively (spinning) for a short period between parallel regions before falling back to passive waiting, so that idle threads are ready to resume work with minimal latency; forcing passive instead makes idle threads yield the CPU immediately. This trades a little peak throughput (the latency of waking a sleeping thread for the next parallel region) for lower CPU usage and less contention with other processes. It is worth benchmarking on a shared or multi-tenant machine; it is unlikely to help on a machine dedicated to Deepwave.

More detailed information specific to each propagator – its equation, its models, its source and receiver channels, and a minimal example – is on that propagator’s own page, under “Propagator reference” in the sidebar: Scalar Wave Equation, Born Modelling with the Scalar Wave Equation, Acoustic Wave Equation, Born Modelling with the Variable-Density Acoustic Wave Equation, Viscoacoustic (attenuating) wave equation, Elastic wave equation, Born Modelling with the Elastic Wave Equation, Elastic VTI (anisotropic) wave equation and Pseudo-acoustic VTI wave equation.