Scalar wave equation

Derivation

To derive the constant density acoustic wave equation (which is a scalar wave equation), we use Newton’s second law (\(F=ma\)) and conservation of mass. Applied to a constant density medium, these yield, to first order,

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

and

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

where \(P\) is pressure, \(\mathbf{f}\) is an external force per unit area (in 2D, or per unit volume in 3D), \(\mathbf{u}\) is particle displacement, \(\rho\) is total density, and \(\rho_0\) is the constant background density. Applying the divergence operator to (1) and the time derivative to (2), we can combine these two equations to yield

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

Pressure and density are proportional, and we will call the constant of proportionality \(c^2\),

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

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

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

The dimensions of each term in 2D are \([ML^{-2}T^{-2}]\).

Implementation

The scalar wave equation implemented in Deepwave is

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

where \(u\) is a scalar field that obeys the wave equation, such as pressure, \(c\) is the wave speed, and \(f\) is a scalar source term. Sources are added to grid cell centres. If \(u\) is pressure (as in (5)) then an explosive point source centred on a grid cell corresponds to a single scalar source in (6).

Deepwave uses central finite differences in time (with second-order accuracy) and space (with a user-specified order of accuracy) to approximately solve this equation on a given time and space domain.

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

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

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

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

and \(a\) and \(b\) are values that are determined for each grid cell depending on its location.

Applying this to (6), and, for simplicity, just using one spatial dimension,

(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 we define new auxiliary variables

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

and

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

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

(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\]

Using (8) in (10) and (11), we can calculate \(p\) and \(z\) using

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

and

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

We see that \(z^t\) depends on the spatial derivative of \(p^t\). This is problematic to calculate, especially on a GPU (where it is preferable for each element’s calculation to be independent from those of its neighbours). We therefore substitute (13) into (14) to get

(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 change to (12), we get

(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\]

With some rearrangement, we can express each timestep in matrix form,

(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}\]

Time sample \(t\) of the output receiver data \(r\) is produced by recording \(u^t\) at the specified receiver locations \(\delta_r\). This means that the receiver data covers the same time range as the source data. Since the wavefield \(u\) isn’t affected by time sample \(t\) of the source until time step \(t+1\) (see (12), where \(f^t\) is added to \(u^{t+1}\)), the first time sample of the receiver data will not be affected by the source, and the last time sample of the source does not affect the receiver data.

The propagator returns \(u^{T}\), \(u^{T-1}\), \(p_y^{T-1}\), \(p_x^{T-1}\), \(z_y^{T-1}\), \(z_x^{T-1}\), and \(r\), where \(T\) is the number of time steps.