Derivation
To derive the constant density acoustic wave equation (which is a scalar wave equation), we use Newton’s second law () and conservation of mass. Applied to a constant density medium, these yield, to first order,
(1)
and
(2)
where is pressure, is an external force per unit area (in 2D, or per unit volume in 3D), is particle displacement, is total density, and 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)
Pressure and density are proportional, and we will call the constant of proportionality ,
(4)
Substituting this into (3) and rearranging yields the scalar wave equation
(5)
The dimensions of each term in 2D are .
Implementation
The scalar wave equation implemented in Deepwave is
(6)
where is a scalar field that obeys the wave equation, such as pressure, is the wave speed, and is a scalar source term. Sources are added to grid cell centres. If 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)
where is an operator defined at time step as
(8)
and and 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)
where we define new auxiliary variables
(10)
and
(11)
Using finite differences in time, with time step interval , the update equation for is thus
(12)
Using (8) in (10) and (11), we can calculate and using
(13)
and
(14)
We see that depends on the spatial derivative of . 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)
Applying the same change to (12), we get
(16)
With some rearrangement, we can express each timestep in matrix form,
(17)
Time sample of the output receiver data is produced by recording at the specified receiver locations . This means that the receiver data covers the same time range as the source data. Since the wavefield isn’t affected by time sample of the source until time step (see (12), where is added to ), 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 , , , , , , and , where is the number of time steps.