# Scaling the wave equation

Solving elastic wave equation in FP16 precision requires scaling the wave equation. In this notebook, we derive the scaling approach explained in the article, and verify it numerically. 

## Derivation

The elastic wave equation is:
$$
\begin{align}
&\partial_t v_i- \frac{1}{\rho}\partial_j \sigma_{ij} =f_i, \\
&\partial_t \sigma_{ij}- (M+2\mu) \partial_k v_k \delta_{ij} - \mu( \partial_j v_i + \partial_i v_j) =s_{ij}.
\end{align}
$$

Multiplying the stress equation by $2^{e_v}$:
$$
\begin{align}
&\partial_t v_i- \frac{1}{\rho}2^{-e_v}2^{e_v}\partial_j \sigma_{ij} =f_i, \\
&2^{e_v}\partial_t \sigma_{ij}- 2^{e_v}(M+2\mu) \partial_k v_k \delta_{ij} - 2^{e_v}\mu( \partial_j v_i + \partial_i v_j) =2^{e_v}s_{ij}.
\end{align}
$$

Defining $\sigma'=2^{e_v}\sigma$, $M'=2^{e_v}M$, $\mu'=2^{e_v}\mu$, $1/\rho'=1/\rho2^{-e_v}$:

$$
\begin{align}
&\partial_t v_i- \frac{1}{\rho'}\partial_j \sigma'_{ij} =f_i, \\
&\partial_t \sigma'_{ij}- (M'+2\mu') \partial_k v_k \delta_{ij} - \mu'( \partial_j v_i + \partial_i v_j) =2^{e_v}s_{ij}.
\end{align}
$$

Multiplying all equations by $2^{-e_v}2^{e_s}$:

$$
\begin{align}
&\partial_t 2^{-e_v}2^{e_s}v_i- \frac{1}{\rho'}\partial_j 2^{-e_v}2^{e_s}\sigma'_{ij} =2^{-e_v}2^{e_s}f_i, \\
&\partial_t 2^{-e_v}2^{e_s}\sigma'_{ij}- (M'+2\mu') \partial_k 2^{-e_v}2^{e_s}v_k \delta_{ij} - \mu'( \partial_j 2^{-e_v}2^{e_s}v_i + \partial_i 2^{-e_v}2^{e_s}v_j) =2^{e_s}s_{ij}.
\end{align}
$$

Defining $\sigma''=2^{-e_v}2^{e_s}\sigma'$, $v''=2^{-e_v}2^{e_s}v$, $s'' = 2^{e_s}s$, $f_i'' = 2^{-e_v}2^{e_s}f_i$:

$$
\begin{align}
&\partial_t v''_i- \frac{1}{\rho'}\partial_j \sigma''_{ij} =f_i'', \\
&\partial_t \sigma''_{ij}- (M'+2\mu') \partial_k v''_k \delta_{ij} - \mu'( \partial_j v''_i + \partial_iv''_j) =s''_{ij}.
\end{align}
$$

The last equation can be solved numericaly. The solution to the initial unscaled PDE can be recovered by the following transformations:

$$
\begin{align}
\sigma_{ij} &= 2^{-e_v}\sigma_{ij}' = 2^{-e_v}2^{e_v}2^{-e_s}\sigma_{ij}'' \\
&=  2^{-e_s}\sigma_{ij}''\\
v_i &= 2^{e_v}2^{-e_s} v_i''
\end{align}
$$

## Numerical verification

We first define the fourth order finite difference operators:


```python
import numpy as np
import copy
def Dp(u):
    return 1.1382*(u[2:-2] - u[1:-3])-0.046414*(u[3:-1] - u[0:-4])
def Dm(u):
    return 1.1382*(u[3:-1] - u[2:-2])-0.046414*(u[4:] - u[1:-3])
```

The kernel for updating the velocity is:


```python
def updatev(v, sigma, M, rho, dt):
    v[2:-2] = v[2:-2] + dt/rho*Dp(sigma)
    return v, sigma
```

and the kernel to update the stresses is given by:


```python
def updates(v, sigma, M, rho, dt):
    sigma[2:-2] =sigma[2:-2] + dt*M*Dm(v)
    return v, sigma
```

We then define random initial values for the velocities, the stresses and the material parameters.


```python
sigma = np.zeros(100) #Stress
v = np.zeros(100) #Velocity
M = 1000*np.random.random(96) #P-wave modulus
rho = 1000*np.random.random(96) #Density
s = 10**-4*np.random.random(96) #Source for stress
f =  10**-1*np.random.random(96) #Source for velocity
dt = np.random.random() #Time step
```

We can first propagate without any scaling for 100 steps.


```python
sigma_1 = copy.copy(sigma)
v_1 = copy.copy(v)
sigma_1[2:-2] += s
v_1[2:-2] += f
for ii in range(100):
    v_1, sigma_1 = updatev(v_1, sigma_1, M, rho, dt)
    v_1, sigma_1 = updates(v_1, sigma_1, M, rho, dt)
```

We then compute the scaling factors.


```python
ev = -np.log2(np.max(M))
es = -np.log2(np.max(s))
```

We perfom the scaling.


```python
Mp = 2**ev * M
rhop = 2**ev * rho
sp = 2**es * s
fp = 2**-ev * 2**es * f
sigma_2 = copy.copy(sigma)
v_2 = copy.copy(v)
```

Propagating with the scaled variables.


```python
sigma_2[2:-2] += sp
v_2[2:-2] += fp
for ii in range(100):
    v_2, sigma_2 = updatev(v_2, sigma_2, Mp, rhop, dt)
    v_2, sigma_2 = updates(v_2, sigma_2, Mp, rhop, dt)
```

Let's compute the difference between the unscaled equation and the scaled equation.


```python
print(np.max(np.abs(sigma_1 - sigma_2 * 2**-es))/np.max(np.abs(sigma_1)))
print(np.max(np.abs(v_1 - v_2 * 2**-es * 2**ev))/np.max(np.abs(v_1)))
```

    7.721878362239146e-15
    7.665104868863169e-15


This difference is the order of the machine precision, which shows that the scaled equation leads to the same solution.

50701bdd0e7e16e88ce18e6109c508e5d9bf5233
ipynb
Jupyter Notebook
Scaling_the_wave_equation.ipynb
GeoCode-polymtl/Seis_float16
5f9660cbdc37e5ab7f6054f7547df2ffb661a81d
MIT
Scaling_the_wave_equation.ipynb
GeoCode-polymtl/Seis_float16
5f9660cbdc37e5ab7f6054f7547df2ffb661a81d
MIT
2020-01-28T22:17:04.000Z
2022-02-09T23:33:07.000Z
Scaling_the_wave_equation.ipynb
GeoCode-polymtl/Seis_float16
5f9660cbdc37e5ab7f6054f7547df2ffb661a81d
MIT
2019-11-27T06:06:04.000Z
2020-06-05T17:18:15.000Z
Qwen/Qwen-72B
1. YES
2. YES