Acoustic Wave Equation

This page derives the variable-density acoustic wave equation that deepwave.acoustic() solves and describes the staggered-grid scheme used to discretise it. It is reference material for readers who already know the scalar propagator and need variable density, or who want the exact numerics; the example below is enough to get a first result running.

Minimal example

deepwave.acoustic() takes wavespeed and density models, a grid spacing, a time step, and a source, and returns the final wavefields, one opaque auxiliary_state holding the absorbing boundary’s memory variables, and the recorded receiver data. Passing the whole output back in as initial_state continues the propagation in a later call (see Usage):

import torch
import deepwave

v = 1500 * torch.ones(40, 40)
rho = 2200 * torch.ones(40, 40)
out = deepwave.acoustic(
    v, rho, grid_spacing=5, dt=0.001,
    source_amplitudes_volume_rate=deepwave.wavelets.ricker(25, 200, 0.001, 0.05).reshape(1, 1, -1),
    source_locations_volume_rate=torch.tensor([[[5, 20]]]),
    receiver_locations_p=torch.tensor([[[35, 20]]]),
    pml_width=10,
    pml_freq=25,
)
out.receiver_amplitudes_p  # the recorded pressure, shape [1, 1, 200]

See Usage for the full set of source/receiver components (pressure and each velocity component) and Units and source conventions for what the source amplitude physically represents.

Derivation

The acoustic wave equation in a heterogeneous medium with variable density \(\rho(\mathbf{x})\) and bulk modulus \(K(\mathbf{x})\) is derived from the linearised equations of momentum conservation and mass conservation (continuity equation).

Momentum Conservation (Equation of Motion):

(1)\[\rho \frac{\partial \mathbf{v}}{\partial t} = -\nabla p + \mathbf{f},\]

where:

  • \(\mathbf{v}\) is the particle velocity vector.

  • \(p\) is the acoustic pressure.

  • \(\mathbf{f}\) is the body force source term (force per unit volume).

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

Mass Conservation (Continuity Equation):

(2)\[\frac{1}{K} \frac{\partial p}{\partial t} = -\nabla \cdot \mathbf{v} + s,\]

where:

  • \(K = \rho c^2\) is the adiabatic bulk modulus, with \(c\) being the acoustic wave speed.

  • \(s\) is a source term (volume injection rate per unit volume, so it has units of 1/time; see Units and source conventions).

Implementation

The system is discretised using a finite-difference time-domain (FDTD) method on a staggered grid.

Staggered Grid Layout

For a 2D domain in the \((y, x)\) plane (using \(y\) for the first dimension and \(x\) for the second to match the array indexing convention):

  • Pressure \(p\) and bulk modulus \(K\) are defined at grid nodes \((y, x)\).

  • Vertical velocity \(v_y\) and density \(\rho_y\) are defined at half-integer steps in \(y\): \((y+1/2, x)\).

  • Horizontal velocity \(v_x\) and density \(\rho_x\) are defined at half-integer steps in \(x\): \((y, x+1/2)\).

Time Integration

Time stepping is performed using an explicit second-order leapfrog scheme. Pressure and particle velocities are staggered in time by half a time step \(h_t/2\).

  1. Update Particle Velocities:

    (3)\[\mathbf{v}^{t+1/2} = \mathbf{v}^{t-1/2} - h_t \frac{1}{\rho} \nabla p^t + h_t \frac{1}{\rho} \mathbf{f}^t.\]

  2. Update Pressure:

    (4)\[p^{t+1} = p^t - h_t K (\nabla \cdot \mathbf{v}^{t+1/2} - s^{t+1/2}).\]

Absorbing Boundaries (CPML)

To mitigate artificial reflections from the domain boundaries, a Convolutional Perfectly Matched Layer (CPML) is employed. This method introduces memory variables to implement the coordinate stretching.

For the velocity update, we replace the gradient \(\nabla p\) with \(\nabla p + \mathbf{\psi}\). The memory variable \(\mathbf{\psi}\) is updated as:

\[\psi_i^t = a_i \psi_i^{t-1} + b_i \partial_i p^t\]

where \(i \in \{y, x\}\).

For the pressure update, we replace the divergence \(\nabla \cdot \mathbf{v}\) with \(\nabla \cdot \mathbf{v} + \sum_i \phi_i\). The memory variable \(\mathbf{\phi}\) is updated as:

\[\phi_i^{t+1/2} = a_i \phi_i^{t-1/2} + b_i \partial_i v_i^{t+1/2}\]

The coefficients \(a\) and \(b\) are determined by the PML damping profile.

Matrix Form

We express the update steps in matrix form, which facilitates the derivation of the adjoint operators. We use \(B\) for buoyancy (\(1/\rho\)) and \(h_t\) for the time step \(h_t\). The notation \(y'\) and \(x'\) denotes locations shifted by half a grid cell (e.g., \(y+1/2\)).

1. Velocity Update

This step updates the particle velocities and the gradient memory variables.

\[\begin{split}\begin{pmatrix} v_y^{t+\frac{1}{2}} \\ v_x^{t+\frac{1}{2}} \\ \psi_y^t \\ \psi_x^t \end{pmatrix} = \begin{pmatrix} -h_t B^{y'x} (1+b^{y'}) \partial_{y'} & 1 & 0 & -h_t B^{y'x} a^{y'} & 0 & h_t B^{y'x} & 0 \\ -h_t B^{yx'} (1+b^{x'}) \partial_{x'} & 0 & 1 & 0 & -h_t B^{yx'} a^{x'} & 0 & h_t B^{yx'} \\ b^{y'} \partial_{y'} & 0 & 0 & a^{y'} & 0 & 0 & 0 \\ b^{x'} \partial_{x'} & 0 & 0 & 0 & a^{x'} & 0 & 0 \end{pmatrix} \begin{pmatrix} p^t \\ v_y^{t-\frac{1}{2}} \\ v_x^{t-\frac{1}{2}} \\ \psi_y^{t-1} \\ \psi_x^{t-1} \\ f_y^t \\ f_x^t \end{pmatrix} \label{acoustic_matrix_v}\end{split}\]

2. Pressure Update

This step updates the pressure and the divergence memory variables.

\[\begin{split}\begin{pmatrix} p^{t+1} \\ \phi_y^{t+\frac{1}{2}} \\ \phi_x^{t+\frac{1}{2}} \end{pmatrix} = \begin{pmatrix} -h_t K^{yx} (1+b^y) \partial_y & -h_t K^{yx} (1+b^x) \partial_x & 1 & -h_t K^{yx} a^y & -h_t K^{yx} a^x & h_t K^{yx} \\ b^y \partial_y & 0 & 0 & a^y & 0 & 0 \\ 0 & b^x \partial_x & 0 & 0 & a^x & 0 \end{pmatrix} \begin{pmatrix} v_y^{t+\frac{1}{2}} \\ v_x^{t+\frac{1}{2}} \\ p^t \\ \phi_y^{t-\frac{1}{2}} \\ \phi_x^{t-\frac{1}{2}} \\ s^{t+\frac{1}{2}} \end{pmatrix} \label{acoustic_matrix_p}\end{split}\]

Specifying sources and receivers

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

out = deepwave.acoustic(
    ...,
    sources={"volume_rate": (amplitudes, source_locations)},
    receivers={"p": receiver_locations},
)

This propagator’s source channels are volume_rate, z, y, x and its receiver channels are p, z, y, x (the z channel is absent in 2-D). "volume_rate" names the quantity the scalar source injects; see Units and source conventions. "p" is still accepted there, with a deepwave.DeepwaveFutureWarning, and is removed in 0.4.0. It is unchanged as a receiver channel name.

The equivalent flat source_amplitudes_*/source_locations_*/ receiver_locations_p keywords are equally permanent and produce bit-identical results; giving one channel in both forms raises. See Usage for the full vocabulary.