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)P+f=ρ0u¨

and

(2)ρ˙=ρ0u˙,

where P is pressure, f is an external force per unit area (in 2D, or per unit volume in 3D), u is particle displacement, ρ is total density, and ρ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)2P+f=ρ¨.

Pressure and density are proportional, and we will call the constant of proportionality c2,

(4)P=c2ρ

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

(5)2P1c2P¨=f.

The dimensions of each term in 2D are [ML2T2].

Implementation

The scalar wave equation implemented in Deepwave is

(6)2u1c2u¨=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)x~=x+ψ,

where ψ is an operator defined at time step t as

(8)ψt=aψt1+b(x)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)c22utx~22utt2=c2ftc2x~(utx+pt)2utt2=c2ftc2(2utx2+ptx+zt)2utt2=c2ft,

where we define new auxiliary variables

(10)pt=ψtut

and

(11)zt=ψt(utx+pt)

Using finite differences in time, with time step interval ht, the update equation for u is thus

(12)ut+1=c2ht2(2utx2+ptx+zt)+2utut1c2ht2ft

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

(13)pt=apt1+butx

and

(14)zt=azt1+b(2utx2+ptx)

We see that zt depends on the spatial derivative of pt. 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)zt=azt1+b(2utx2+(apt1+butx)x)

Applying the same change to (12), we get

(16)ut+1=c2ht2(2utx2+(apt1+butx)x+(azt1+b(2utx2+(apt1+butx)x)))+2utut1c2ht2ft

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

(17)(ut+1utztptrt)=(c2ht2(1+b)((1+b)x2+x(b)x)+21c2ht2ac2ht2(1+b)(xa)c2ht210000b((1+b)x2+x(b)x)0ab(xa)0bx00a0δr0000)(utut1zt1pt1ft)

Time sample t of the output receiver data r is produced by recording ut at the specified receiver locations δ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 ft is added to ut+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 uT, uT1, pyT1, pxT1, zyT1, zxT1, and r, where T is the number of time steps.