We show a derivation for the time update expression used for the constant density acoustic solver. You can compare the end result of this derivation in the last equation line below with lines 58-59 in the file examples/seismic/acoustic/operators.py.
| Symbol | Description | Dimensionality |
|---|---|---|
| $\delta t$ | Temporal sampling interval | constant |
| $m(x,y,z)$ | slowness squared | function of space |
| $\eta(x,y,z)$ | Damping coefficient | function of space |
| $u(t,x,y,z)$ | Pressure wavefield | function of time and space |
| $q(t,x,y,z)$ | Source term | function of time, localized in space |
| $\partial_{t}$ | first derivative wrt $t$ | time |
| $\partial_{tt}$ | second derivative wrt $t$ | time |
| $\nabla^2$ | Laplacian operator | space |
To implement the Devito modeling operator we define the equation used to update the pressure wavefield as a function of time. What follows is a bit of algebra using the wave equation and finite difference approximations to time derivatives to express the pressure wavefield forward in time $u(t+\delta t)$ as a function of the current $u(t)$ and previous $u(t-\delta t)$ pressure wavefields.
The first order accurate forward approximation to the first time derivative involves two wavefields: $u(t-\delta t)$, and $u(t)$. We can use this expression as is.
$$ \partial_{t}\ u(t) = \frac{u(t+\delta t) - u(t)}{\delta t} $$The second order accurate centered approximation to the second time derivative involves three wavefields: $u(t-\delta t)$, $u(t)$, and $u(t+\delta t)$.
$$ \partial_{tt}\ u(t) = \frac{u(t+\delta t) - 2\ u(t) + u(t-\delta t)}{\delta t^2} $$In order to advance our finite difference solution in time, we solve for $u(t+\delta t)$.
$$ u(t+\delta t) = \delta t^2\ \partial_{tt}\ u(t) + 2\ u(t) - u(t-\delta t) $$Our Cerjan (reference below) damped wave equation, which we solve for $\partial_{tt}$:
$$ \begin{aligned} m\ \partial_{tt}\ u(t) + \eta\ \partial_{t}\ u(t) &= \nabla^2 u(t) + q(t) \\[10pt] \partial_{tt}\ u(t) &= \frac{1}{m} \Bigr[ \nabla^2 u(t) + q(t) - \eta\ \partial_{t}\ u(t) \Bigr] \end{aligned} $$Next we plug in the expression for $\partial_{tt}\ u$ (from the wave equation) and $\partial_{t}\ u$ (from the numerical derivative) into the the time update expression for $u(t+\delta t)$ from step 3.
$$ \begin{aligned} u(t+\delta t) &= \frac{\delta t^2}{m} \Bigr[ \nabla^2 u(t) + q(t) - \frac{\eta}{\delta t} \bigr\{ u(t+\delta t) - u(t) \bigr\} \Bigr] + 2\ u(t) - u(t-\delta t) \end{aligned} $$Finally we simplify this expression to the form used in the Devito Operator.
Please compare the last equation above with lines 58-59 in examples/seismic/acoustic/operators.py
eq_time = ((lap + q) * s**2 + s * damp * field +
m * (2 * field - prev))/(s * damp + m)
In [ ]: