Scalar Wave Equation

This page derives the constant-density acoustic scalar wave equation that deepwave.scalar() and deepwave.scalar_born() solve, and shows the exact update that Deepwave applies at each time step, including the PML absorbing boundary. It is reference material for readers who want the underlying physics and numerics; if you only want to run a propagation, the example below and Quick Start are enough.

Minimal example

deepwave.scalar() takes a wavespeed model, a grid spacing, a time step, and a source, and returns the final wavefields and the recorded receiver data:

import torch
import deepwave

v = 1500 * torch.ones(40, 40)
out = deepwave.scalar(
    v, grid_spacing=5, dt=0.001,
    source_amplitudes=deepwave.wavelets.ricker(25, 200, 0.001, 0.05).reshape(1, 1, -1),
    source_locations=torch.tensor([[[5, 20]]]),
    receiver_locations=torch.tensor([[[35, 20]]]),
    pml_width=10,
    pml_freq=25,
)
out.receiver_amplitudes  # shape [1, 1, 200]

See Quick Start for a complete runnable example that also plots the result, and Usage for the full set of arguments.

Derivation

We derive the constant density acoustic wave equation (a scalar wave equation) by applying Newton’s second law (\(F=ma\)) and the principle of conservation of mass. For a constant density medium, these principles yield, to first order:

(1)\[-\nabla P + \mathbf{s} = \rho_0\ddot{\mathbf{u}}\]

and

(2)\[\dot{\rho} = -\rho_0\nabla\cdot\dot{\mathbf{u}},\]

where:

• \(P\) is pressure.

• \(\mathbf{s}\) is an external body force per unit volume (in both 2D and 3D: every term in (1) has units of force per unit volume). See Units and source conventions for what changes between 2D and 3D, which is the physical interpretation of a grid cell’s amplitude rather than its units.

• \(\mathbf{u}\) is particle displacement.

• \(\rho\) is the total density.

• \(\rho_0\) is the constant background density.

By applying the divergence operator to (1) and the time derivative to (2), we can combine these two equations to obtain:

(3)\[-\nabla^2 P + \nabla\cdot\mathbf{s} = -\ddot{\rho}.\]

Assuming pressure and density are proportional, with a constant of proportionality \(c^2\) (where \(c\) is the wave speed):

(4)\[P = c^2\rho\]

Substituting this relationship into (3) and rearranging the terms yields the scalar wave equation:

(5)\[\nabla^2 P -\frac{1}{c^2}\ddot{P} = \nabla\cdot\mathbf{s}.\]

Implementation

Deepwave implements the scalar wave equation in the form:

(6)\[\nabla^2 u - \frac{1}{c^2}\ddot{u} = f,\]

where \(u\) is a scalar field (e.g., pressure) that obeys the wave equation, \(c\) is the wave speed, and \(f\) is a scalar source term. Sources are applied at grid cell centres. For instance, if \(u\) represents pressure (as in (5)), an explosive point source centred on a grid cell corresponds to a single scalar source in (6). Note that \(f\) is the divergence \(\nabla\cdot\mathbf{s}\) of the body force of (1), so it has units of \([u]\) per unit area, and that the value supplied for a cell is a density: it is not divided by the cell volume, so the equivalent continuum point-source strength is the amplitude times the cell volume (see Units and source conventions).

Deepwave employs central finite differences for both time (second-order accurate) and space (user-specified order of accuracy) to numerically solve this equation within a defined time and space domain.

Perfectly Matched Layer (PML)

To prevent unwanted reflections from the spatial boundaries of the computational domain, Deepwave uses a Perfectly Matched Layer (PML). For the scalar wave equation, Deepwave adopts the method proposed by Pasalic and McGarry. In this approach, spatial derivatives are modified within the PML regions:

(7)\[\frac{\partial}{\partial \tilde{x}} = \frac{\partial}{\partial x} + \psi,\]

where \(\psi\) is an auxiliary operator defined at time step \(t\) as:

(8)\[\psi^t = a\psi^{t-1} + b\left(\frac{\partial}{\partial x}\right)_t,\]

with \(a\) and \(b\) being grid-cell-dependent values. Applying this modification to (6) (simplified to one spatial dimension for illustration):

(9)\[\begin{split}\begin{align} c^2\frac{\partial^2 u^t}{\partial \tilde{x}^2} - \frac{\partial^2 u^t}{\partial t^2} &= c^2 f^t \\ c^2\frac{\partial}{\partial \tilde{x}}\left(\frac{\partial u^t}{\partial x} + p^t\right) - \frac{\partial^2 u^t}{\partial t^2} &= c^2 f^t \\ c^2\left(\frac{\partial^2 u^t}{\partial x^2} + \frac{\partial p^t}{\partial x} + z^t\right) - \frac{\partial^2 u^t}{\partial t^2} &= c^2 f^t, \end{align}\end{split}\]

where \(p^t\) and \(z^t\) are new auxiliary variables defined as:

(10)\[p^t = \psi^t u^t\]

(11)\[z^t = \psi^t \left(\frac{\partial u^t}{\partial x} + p^t\right)\]

Using finite differences in time with a time step interval \(h_t\), the update equation for \(u\) is:

(12)\[u^{t+1} = c^2h_t^2\left(\frac{\partial^2 u^t}{\partial x^2} + \frac{\partial p^t}{\partial x} + z^t\right) + 2u^t - u^{t-1} - c^2h_t^2 f^t\]

To ensure computational efficiency, especially on GPUs where independent calculations are preferred, the auxiliary variable update equations are rearranged. Using (8) in (10) and (11), we can calculate \(p\) and \(z\) as:

(13)\[p^t = ap^{t-1} + b\frac{\partial u^t}{\partial_x}\]

(14)\[z^t = az^{t-1} + b\left(\frac{\partial^2 u^t}{\partial_x^2} + \frac{\partial p^t}{\partial x}\right)\]

The direct dependence of \(z^t\) on the spatial derivative of \(p^t\) is problematic for parallel computation. Substituting (13) into (14) yields an expression for \(z^t\) that is more amenable to parallelisation:

(15)\[z^t = az^{t-1} + b\left(\frac{\partial^2 u^t}{\partial_x^2} + \frac{\partial \left(ap^{t-1} + b\frac{\partial u^t}{\partial_x}\right)}{\partial x}\right)\]

Applying the same transformation to (12) results in:

(16)\[u^{t+1} = c^2h_t^2\left(\frac{\partial^2 u^t}{\partial x^2} + \frac{\partial \left(ap^{t-1} + b\frac{\partial u^t}{\partial_x}\right)}{\partial x} + \left(az^{t-1} + b\left(\frac{\partial^2 u^t}{\partial_x^2} + \frac{\partial \left(ap^{t-1} + b\frac{\partial u^t}{\partial_x}\right)}{\partial x}\right)\right)\right) + 2u^t - u^{t-1} - c^2h_t^2 f^t\]

Matrix Form of Time Stepping

Through further rearrangement, each time step can be expressed in a compact matrix form. This representation facilitates the internal computation of the wavefield and auxiliary variables:

(17)\[\begin{split}\begin{pmatrix} u^{t+1} \\ u^t \\ z^t \\ p^t \\ r^t \end{pmatrix} = \begin{pmatrix} c^2h_t^2(1+b)\left((1+b)\partial_x^2 +\partial_x(b)\partial_x\right) + 2 & -1 & c^2h_t^2a & c^2h_t^2(1+b)\left(\partial_x a\right) & -c^2h_t^2 \\ 1 & 0 & 0 & 0 & 0\\ b\left((1+b)\partial_{x}^2+\partial_x(b)\partial_x\right) & 0 & a & b\left(\partial_x a\right) & 0\\ b\partial_x & 0 & 0 & a & 0 \\ \delta_r & 0 & 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} u^t \\ u^{t-1} \\ z^{t-1} \\ p^{t-1} \\ f^t \end{pmatrix}\end{split}\]

Output Data

Time sample \(t\) of the output receiver data \(r\) is generated by recording \(u^t\) at the specified receiver locations \(\delta_r\). This ensures that the receiver data covers the same time range and has the same time sampling as the source data. It’s important to note that the wavefield \(u\) is not affected by time sample \(t\) of the source until time step \(t+1\) (as shown in (12), where \(f^t\) is added to \(u^{t+1}\)). Consequently, the first time sample of the receiver data will not be influenced by the source, and the last time sample of the source will not affect the receiver data.

The propagator returns the final wavefield states \(u^{T}\) and \(u^{T-1}\) (wavefield_0/wavefield_m1), the absorbing boundary’s memory variables \(p_y^{T-1}\), \(p_x^{T-1}\), \(z_y^{T-1}\), \(z_x^{T-1}\) carried opaquely as one auxiliary_state field, and the recorded receiver data \(r\), where \(T\) is the total number of time steps. Passing the whole output back in as initial_state continues propagation in a subsequent call; see Usage.

Specifying sources and receivers

The channel-keyed mapping form is the canonical way to specify a survey:

out = deepwave.scalar(
    ...,
    sources={"u": (amplitudes, source_locations)},
    receivers={"u": receiver_locations},
)

This propagator’s source channels are u and its receiver channels are u (the z channel is absent in 2-D).

The equivalent flat source_amplitudes_*/source_locations_*/ receiver_locations keywords remain permanently supported and produce bit-identical results; giving one channel in both forms raises. See Usage for the full vocabulary.