In [1]:
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.
In [2]:
nx = 250 # space discretisation
dt = 0.1 # time discretisation
nt = 1000 # number of timesteps
grid = Grid(shape=(nx, nx), extent=(1000., 1000.), dtype=np.double)
Equation¶
The PDE we want to solve is $$\frac{d^2u}{dt^2} = c\nabla^2 u + f$$ but where $c$ is a time-dependant function. We show here an example of a cooling medium, where the propagation speed decreases following the equation $$c(x, y, t) = c_i \times (300-30\sqrt{2t})$$ up to $t=40$, and then $c(x, y, t) = c_i \times (300-30\sqrt{80})$.
$c_i$ is an initial field composed of four layers. For the sake of comparison, we show the evolution of the same perturbation in a field where $c(x, y, t)$ is constant and follows the law $$c(x, y, t) = c_i \times 300$$
In [3]:
t = np.array([np.full((nx, nx), i)*dt for i in range(nt)])
c = Function(name="c", grid=grid, space_order=2)
c.data[:nx//4] = 1
c.data[nx//4:nx//2] = 0.5
c.data[nx//2:3*nx//4] = 1
c.data[3*nx//4:] = 1.5
c_dec = TimeFunction(name="c_dec", grid=grid, space_order=2, save=nt)
c_dec.data[:] = 300-30*((t*2)**0.5)
c_dec.data[2*nt//5:] = c_dec.data[2*nt//5, 0, 0]
c_cons = TimeFunction(name="c_dec", grid=grid, save=nt)
c_cons.data[:] = 300
u_dec = TimeFunction(name="u_d", grid=grid, space_order=2, time_order=2, save=nt)
u_cons = TimeFunction(name="u_c", grid=grid, space_order=2, time_order=2, save=nt)
f_dec = TimeFunction(name="f_d", grid=grid, space_order=2, time_order=2)
f_dec.data[0:5, nx//2-5:nx//2+5, nx//2-5:nx//2+5] = 1
f_dec.data[5:10, nx//2-5:nx//2+5, nx//2-5:nx//2+5] = 0
f_cons = TimeFunction(name="f_c", grid=grid, space_order=2, time_order=2)
f_cons.data[0:5, nx//2-5:nx//2+5, nx//2-5:nx//2+5] = 1
f_cons.data[5:10, nx//2-5:nx//2+5, nx//2-5:nx//2+5] = 0
eq_dec = Eq(u_dec.dt2, f_dec + c_dec*c*u_dec.laplace)
eq_cons = Eq(u_cons.dt2, f_cons + c_cons*c*u_cons.laplace)
stencil_dec = solve(eq_dec, u_dec.forward)
stencil_cons = solve(eq_cons, u_cons.forward)
op_dec = Operator([Eq(u_dec.forward, stencil_dec)], opt='noop')
op_cons = Operator([Eq(u_cons.forward, stencil_cons)], opt='noop')
op_dec = op_dec.apply(dt=dt)
op_cons = op_cons.apply(dt=dt)
The factor $a$ is decreasing through time:
In [4]:
plt.plot([i[0, 0] for i in t], [i[0, 0] for i in c_dec.data])
plt.xlabel("time")
plt.ylabel("a")
plt.title("Multiplication factor through time")
plt.grid()
plt.show()
We plot the propagation of the two waves, depending of the propagation field evolution:
In [5]:
fig, axes = plt.subplots(2, 6, figsize=(30, 10))
for i in range(5):
axes[0, i].imshow(u_cons.data[i*200], vmin=np.min(u_cons.data), vmax=np.max(u_cons.data), cmap=get_black_cmap())
axes[0, i].imshow(c.data*c_cons.data[i*200], cmap="gray", vmin=0, vmax=450, alpha=0.4)
axes[0, i].set_title("t = " + str(i*200*dt))
axes[1, i].imshow(u_dec.data[i*200], vmin=np.min(u_cons.data), vmax=np.max(u_cons.data), cmap=get_black_cmap())
axes[1, i].imshow(c.data*c_dec.data[i*200], cmap="gray", vmin=0, vmax=450, alpha=0.4)
axes[0, 5].plot([i[0, 0] for i in t], [i[0, 0] for i in c_cons.data])
axes[0, 5].set_xlabel("time")
axes[0, 5].set_ylabel("a")
axes[0, 5].set_title("a multiplication factor")
axes[0, 5].grid()
axes[1, 5].plot([i[0, 0] for i in t], [i[0, 0] for i in c_dec.data])
axes[1, 5].set_xlabel("time")
axes[1, 5].set_ylabel("a")
axes[1, 5].grid()
plt.tight_layout()
plt.savefig("varying_c.jpg", dpi=250)
plt.show()
We generate a video containing the simulation:
In [10]:
ddt = 5
img = u_cons.data[::ddt]
img1 = u_dec.data[::ddt]
name = "varying_c"
for i in tqdm(range(len(img))):
fig, ax = plt.subplots(2, 2, figsize=(10, 10))
im = ax[0, 0].imshow(img[i], cmap=get_black_cmap(), vmin=np.min(u_cons.data), vmax=np.max(u_cons.data))
ax[0, 0].imshow(c.data*c_cons.data[::ddt][i], cmap="gray", vmin=0, vmax=450, alpha=0.4)
ax[0, 0].axis('off')
plt.colorbar(im, shrink=0.75, ax=ax[0, 0])
ax[0, 0].set_title("t = "+str(i*ddt*dt))
ax[0, 1].plot([i[0, 0] for i in t], [i[0, 0] for i in c_cons.data])
ax[0, 1].set_ylim([-50, 350])
ax[0, 1].scatter([i[0, 0] for i in t][i*ddt], [i[0, 0] for i in c_cons.data][i*ddt], color="red")
ax[0, 1].set_xlabel("time")
ax[0, 1].set_ylabel("a")
ax[0, 1].set_title("a multiplication factor")
ax[0, 1].grid()
im1 = ax[1, 0].imshow(img1[i], cmap=get_black_cmap(), vmin=np.min(u_cons.data), vmax=np.max(u_cons.data))
ax[1, 0].imshow(c.data*c_dec.data[::ddt][i], cmap="gray", vmin=0, vmax=450, alpha=0.4)
ax[1, 0].axis('off')
plt.colorbar(im1, shrink=0.75, ax=ax[1, 0])
ax[1, 1].plot([i[0, 0] for i in t], [i[0, 0] for i in c_dec.data])
ax[1, 1].set_ylim([-50, 350])
ax[1, 1].scatter([i[0, 0] for i in t][i*ddt], [i[0, 0] for i in c_dec.data][i*ddt], color="red")
ax[1, 1].set_xlabel("time")
ax[1, 1].set_ylabel("a")
ax[1, 1].grid()
plt.savefig(name + "%02d.png" % i, dpi=100)
plt.close()
call([
'ffmpeg', '-loglevel', 'panic', '-framerate', str(int(ddt / dt)), '-i', name + '%02d.png', '-r', '32', '-pix_fmt',
'yuv420p', name + ".mp4", '-y'
])
for file_name in glob(name + "*.png"):
remove(file_name)
0%| | 0/200 [00:00<?, ?it/s]
In [12]:
%%HTML
<video width="810" controls>
<source src="varying_c.mp4" type="video/mp4">
</video>
In [ ]: