#!/usr/bin/env python
# coding: utf-8

# # Time update for *Cerjan* damped constant density acoustics with forward approximation for $\partial_t$
# 
# 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```.

# ## Table of symbols
# 
# | 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 |

# ## A word about notation 
# 
# For clarity in the following derivation we will drop the space notatation for certain variables:
# - $m(x,y,z) \rightarrow m$
# - $\eta(x,y,z) \rightarrow \eta$
# - $u(t,x,y,z) \rightarrow u(t)$
# - $q(t,x,y,z) \rightarrow q(t)$

# ## The time update equation
# 
# 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.
# 
# #### 1. First order numerical derivative (forward):
# 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}
# $$
# <br>
# 
# #### 2. Second order numerical derivative:
# 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}
# $$
# <br>
# 
# #### 3. Second order time update:
# 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)
# $$
# <br>
# 
# #### 4. Damped wave equation:
# 
# 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}
# $$
# 
# #### 5. Time update:
# 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}
# $$
# 
# #### 6. Simplify:
# 
# Finally we simplify this expression to the form used in the Devito ```Operator```.
# 
# $$
# \begin{aligned}
# \left(1 + \frac{\delta t\ \eta}{m}\right) u(t+\delta t) &= 
#         \frac{\delta t^2}{m} \Bigr\{ \nabla^2 u(t) + q(t) \Bigr\}
#         + \frac{\delta t\ \eta}{m}\ u(t) + 2\ u(t) - u(t-\delta t) \\[15pt]
# u(t+\delta t) &=
#     \left( \frac{1}{m+\delta t\ \eta} \right) \Bigr[
#         \delta t^2 \Bigr\{ \nabla^2 u(t) + q(t) \Bigr\}
#         + \delta t\ \eta\ u(t) + m\ \left[2\ u(t) - u(t-\delta t) \right]
#     \Bigr]
# \end{aligned}
# $$
# 
# 
# #### 7. Compare:
# 
# Please compare the last equation above with [lines 58-59 in examples/seismic/acoustic/operators.py](https://github.com/devitocodes/devito/blob/master/examples/seismic/acoustic/operators.py#L58-L59)
# 
# ```
# eq_time = ((lap + q) * s**2 + s * damp * field +
#            m * (2 * field - prev))/(s * damp + m)```

# ## References
# 
# - **A nonreflecting boundary condition for discrete acoustic and elastic wave equations** (1985)
# <br>Charles Cerjan, Dan Kosloft. Ronnie Kosloff, and Moshe Resheq
# <br> Geophysics, Vol. 50, No. 4
# <br>https://library.seg.org/doi/pdfplus/10.1190/segam2016-13878451.1

# In[ ]: