Libraries import

import numpy as np
from devito import *
import matplotlib.pyplot as plt
from pyawd import Marmousi
from pyawd.utils import *
from pyawd.GenerateVideo import generate_quiver_video
from tqdm.notebook import tqdm
from glob import glob
from subprocess import call
from os import remove

Grid definition

We generate a Grid in 2D which is $1000$ m large and long.

nx = 32 # space discretisation
dt = 0.01 # time discretisation
nt = 1000 # number of timesteps

grid = Grid(shape=(nx, nx), extent=(1000., 1000.), dtype=np.double)

External force

s_x, s_y = create_explosive_source(nx)

s_t = np.exp(-(dt)*(np.arange(nt)-(nt//10))**2)

fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(10, 5))

ax[0].quiver(np.arange(0, nx, 1), np.arange(0, nx, 1), s_x, s_y, scale=10)
ax[0].set_title("Example if external force field")
ax[0].set_xlabel("x")
ax[0].set_ylabel("y")

ax[1].plot(np.arange(nt)*dt, s_t)
ax[1].set_title("External force scaling factor evolution")
ax[1].set_xlabel("Time")
ax[1].set_ylabel("Scaling factor")

plt.tight_layout()
plt.savefig("external_force.jpg")
plt.show()

s_x_t = np.array([s_x*s_t[t] for t in range(len(s_t))])
s_y_t = np.array([s_y*s_t[t] for t in range(len(s_t))])

We define the external force evolution:

f = VectorTimeFunction(name='f', grid=grid, space_order=1, save=nt, time_order=1)
f[0].data[:] = s_x_t
f[1].data[:] = s_y_t

Equation

The PDE we want to solve is $$\frac{d^2u}{dt^2} = c\nabla^2 u + f$$

u = VectorTimeFunction(name='u', grid=grid, space_order=2, save=nt, time_order=2)

c = Function(name="c", grid=grid)
c.data[:] = Marmousi(nx).get_data()
c.data[:] *= (100/np.max(c.data[:]))

eq = Eq(u.dt2, f + (c**2)*u.laplace)

stencil = solve(eq, u.forward)

op = Operator([Eq(u.forward, stencil)], opt='noop')

op.apply(dt=dt)

PerformanceSummary([(PerfKey(name='section0', rank=None),
                     PerfEntry(time=0.012278999999999981, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])

eq

$\displaystyle \left[\begin{matrix}\frac{\partial^{2}}{\partial time^{2}} u_x(time, x + h_x/2, y)\\\frac{\partial^{2}}{\partial time^{2}} u_y(time, x, y + h_y/2)\end{matrix}\right] = \left[\begin{matrix}\left(\frac{\partial^{2}}{\partial x^{2}} u_x(time, x + h_x/2, y) + \frac{\partial^{2}}{\partial y^{2}} u_x(time, x + h_x/2, y)\right) c(x, y) + f_x(time, x + h_x/2, y)\\\left(\frac{\partial^{2}}{\partial x^{2}} u_y(time, x, y + h_y/2) + \frac{\partial^{2}}{\partial y^{2}} u_y(time, x, y + h_y/2)\right) c(x, y) + f_y(time, x, y + h_y/2)\end{matrix}\right]$

from sympy import latex
print(latex(eq))

\left[\begin{matrix}\frac{\partial^{2}}{\partial time^{2}} u_x(time, x + h_x/2, y)\\\frac{\partial^{2}}{\partial time^{2}} u_y(time, x, y + h_y/2)\end{matrix}\right] = \left[\begin{matrix}\left(\frac{\partial^{2}}{\partial x^{2}} u_x(time, x + h_x/2, y) + \frac{\partial^{2}}{\partial y^{2}} u_x(time, x + h_x/2, y)\right) c(x, y) + f_x(time, x + h_x/2, y)\\\left(\frac{\partial^{2}}{\partial x^{2}} u_y(time, x, y + h_y/2) + \frac{\partial^{2}}{\partial y^{2}} u_y(time, x, y + h_y/2)\right) c(x, y) + f_y(time, x, y + h_y/2)\end{matrix}\right]