from devito import *
from examples.seismic.source import RickerSource, Receiver, TimeAxis
from examples.seismic import plot_image, demo_model
import numpy as np
import matplotlib.pyplot as plt
from sympy import init_printing, latex
init_printing(use_latex='mathjax')
# Some ploting setup
plt.rc('text', usetex=True)
plt.rc('font', family='serif')
plt.rc('xtick', labelsize=20)
plt.rc('ytick', labelsize=20)06 - Elastic wave equation implementation on a staggered grid
This second elastic tutorial extends the previous constant parameter implementation to varying parameters (Lame parameters) and takes advantage of the Tensorial capabilities of Devito to write the elastic wave equation following its mathematical definition. The staggering is automated via the TensorFunction API.
Explosive source
We will first attempt to replicate the explosive source test case described in [1], Figure 4. We start by defining the source signature \(g(t)\), the derivative of a Gaussian pulse, given by Eq 4:
\[g(t) = -2 \alpha(t - t_0)e^{-\alpha(t-t_0)^2}\]
#NBVAL_IGNORE_OUTPUT
# Initial grid: 3km x 3km, with spacing 10m
nlayers = 5
so = 8
model = demo_model(preset='layers-elastic', nlayers=nlayers, shape=(301, 301), spacing=(10., 10.),
space_order=so)Operator `initdamp` ran in 0.01 s#NBVAL_SKIP
aspect_ratio = model.shape[0]/model.shape[1]
plt_options_model = {'cmap': 'jet', 'extent': [model.origin[0], model.origin[0] + model.domain_size[0],
model.origin[1] + model.domain_size[1], model.origin[1]]}
fig, ax = plt.subplots(nrows=3, ncols=1, figsize=(15, 15))
slices = [slice(model.nbl, -model.nbl), slice(model.nbl, -model.nbl)]
img1 = ax[0].imshow(np.transpose(model.lam.data[slices]), vmin=1.5**2, vmax=4.0**2, **plt_options_model)
fig.colorbar(img1, ax=ax[0])
ax[0].set_title(r"First Lam\'e parameter $\lambda$", fontsize=20)
ax[0].set_xlabel('X (m)', fontsize=20)
ax[0].set_ylabel('Depth (m)', fontsize=20)
ax[0].set_aspect('auto')
img2 = ax[1].imshow(np.transpose(model.mu.data[slices]), vmin=0, vmax=15, **plt_options_model)
fig.colorbar(img2, ax=ax[1])
ax[1].set_title(r"Shear modulus $\mu$", fontsize=20)
ax[1].set_xlabel('X (m)', fontsize=20)
ax[1].set_ylabel('Depth (m)', fontsize=20)
ax[1].set_aspect('auto')
img3 = ax[2].imshow(1/np.transpose(model.b.data[slices]), vmin=1.0, vmax=3.0, **plt_options_model)
fig.colorbar(img3, ax=ax[2])
ax[2].set_title(r"Density $\rho$", fontsize=20)
ax[2].set_xlabel('X (m)', fontsize=20)
ax[2].set_ylabel('Depth (m)', fontsize=20)
ax[2].set_aspect('auto')
plt.tight_layout()
# Timestep size from Eq. 7 with V_p=6000. and dx=100
t0, tn = 0., 2000.
dt = model.critical_dt
time_range = TimeAxis(start=t0, stop=tn, step=dt)
src = RickerSource(name='src', grid=model.grid, f0=0.015, time_range=time_range)
src.coordinates.data[:] = [1500., 10.]#NBVAL_SKIP
src.show()
Vectorial form
While conventional litterature writes the elastic wave-equation as a set of scalar PDEs, the higher level representation comes from Hooke’s law and the equation of motion and writes as:
\[\begin{cases} &\frac{dv}{dt} = \nabla . \tau \\ &\frac{d \tau}{dt} = \lambda tr(\nabla v) \mathbf{I} + \mu (\nabla v + (\nabla v)^T) \end{cases}\]and as \(tr(\nabla v)\) is the divergence of \(v\) we can reqrite it as
\[\begin{cases} &\frac{dv}{dt} = \nabla . \tau \\ &\frac{d \tau}{dt} = \lambda \text{diag}(\nabla . v) + \mu (\nabla v + (\nabla v)^T) \end{cases}\]where \(v\) is a vector valued function:
\(v(t, x, y) = (v_x(t, x, y), v_y(t, x, y)\)
and the stress \(\tau\) is a symmetric tensor valued function:
\(\tau(t, x, y) = \begin{bmatrix}\tau_{xx}(t, x, y) & \tau_{xy}(t, x, y)\\\tau_{xy}t, x, y) & \tau_{yy}(t, x, y)\end{bmatrix}\)
We show in the following how to setup the elastic wave-equation form Devito’s high-level tensorial types.
# Now we create the velocity and pressure fields
x, z = model.grid.dimensions
t = model.grid.stepping_dim
time = model.grid.time_dim
s = time.spacing
v = VectorTimeFunction(name='v', grid=model.grid, space_order=so, time_order=1)
tau = TensorTimeFunction(name='t', grid=model.grid, space_order=so, time_order=1)# The source injection term
src_xx = src.inject(field=tau.forward[0, 0], expr=s*src)
src_zz = src.inject(field=tau.forward[1, 1], expr=s*src)
# The receiver
nrec = 301
rec = Receiver(name="rec", grid=model.grid, npoint=nrec, time_range=time_range)
rec.coordinates.data[:, 0] = np.linspace(0., model.domain_size[0], num=nrec)
rec.coordinates.data[:, -1] = 5.
rec2 = Receiver(name="rec2", grid=model.grid, npoint=nrec, time_range=time_range)
rec2.coordinates.data[:, 0] = np.linspace(0., model.domain_size[0], num=nrec)
rec2.coordinates.data[:, -1] = 3000.0/nlayers
rec3 = Receiver(name="rec3", grid=model.grid, npoint=nrec, time_range=time_range)
rec3.coordinates.data[:, 0] = np.linspace(0., model.domain_size[0], num=nrec)
rec3.coordinates.data[:, -1] = 3000.0/nlayers
rec_term = rec.interpolate(expr=tau[0, 0] + tau[1, 1])
rec_term += rec2.interpolate(expr=v[1])
rec_term += rec3.interpolate(expr=v[0])#NBVAL_SKIP
from examples.seismic import plot_velocity
plot_velocity(model, source=src.coordinates.data,
receiver=rec.coordinates.data[::10, :])
plot_velocity(model, source=src.coordinates.data,
receiver=rec2.coordinates.data[::10, :])
# Now let's try and create the staggered updates
# Lame parameters
l, mu, ro = model.lam, model.mu, model.b
# First order elastic wave equation
pde_v = v.dt - ro * div(tau)
pde_tau = tau.dt - l * diag(div(v.forward)) - mu * (grad(v.forward) + grad(v.forward).transpose(inner=False))
# Time update
u_v = Eq(v.forward, model.damp * solve(pde_v, v.forward))
u_t = Eq(tau.forward, model.damp * solve(pde_tau, tau.forward))
op = Operator([u_v] + [u_t] + src_xx + src_zz + rec_term)v.time_order\(\displaystyle 1\)
ro._eval_at(v[0]).evaluate\(\displaystyle 0.5 b(x, y) + 0.5 b(x + h_x, y)\)
We can now see that both the particle velocities and stress equations are vectorial and tensorial equations. Devito takes care of the discretization and staggered grids automatically for these types of object.
u_v\(\displaystyle \left[\begin{matrix}v_x(t + dt, x + h_x/2, y)\\v_y(t + dt, x, y + h_y/2)\end{matrix}\right] = \left[\begin{matrix}dt \left(\left(\frac{\partial}{\partial x} t_xx(t, x, y) + \frac{\partial}{\partial y} t_xy(t, x + h_x/2, y + h_y/2)\right) b(x, y) + \frac{v_x(t, x + h_x/2, y)}{dt}\right) damp(x, y)\\dt \left(\left(\frac{\partial}{\partial x} t_xy(t, x + h_x/2, y + h_y/2) + \frac{\partial}{\partial y} t_yy(t, x, y)\right) b(x, y) + \frac{v_y(t, x, y + h_y/2)}{dt}\right) damp(x, y)\end{matrix}\right]\)
u_t\(\displaystyle \left[\begin{matrix}t_xx(t + dt, x, y) & t_xy(t + dt, x + h_x/2, y + h_y/2)\\t_xy(t + dt, x + h_x/2, y + h_y/2) & t_yy(t + dt, x, y)\end{matrix}\right] = \left[\begin{matrix}dt \left(\left(\frac{\partial}{\partial x} v_x(t + dt, x + h_x/2, y) + \frac{\partial}{\partial y} v_y(t + dt, x, y + h_y/2)\right) lam(x, y) + 2 mu(x, y) \frac{\partial}{\partial x} v_x(t + dt, x + h_x/2, y) + \frac{t_xx(t, x, y)}{dt}\right) damp(x, y) & dt \left(\left(\frac{\partial}{\partial y} v_x(t + dt, x + h_x/2, y) + \frac{\partial}{\partial x} v_y(t + dt, x, y + h_y/2)\right) mu(x, y) + \frac{t_xy(t, x + h_x/2, y + h_y/2)}{dt}\right) damp(x, y)\\dt \left(\left(\frac{\partial}{\partial y} v_x(t + dt, x + h_x/2, y) + \frac{\partial}{\partial x} v_y(t + dt, x, y + h_y/2)\right) mu(x, y) + \frac{t_xy(t, x + h_x/2, y + h_y/2)}{dt}\right) damp(x, y) & dt \left(\left(\frac{\partial}{\partial x} v_x(t + dt, x + h_x/2, y) + \frac{\partial}{\partial y} v_y(t + dt, x, y + h_y/2)\right) lam(x, y) + 2 mu(x, y) \frac{\partial}{\partial y} v_y(t + dt, x, y + h_y/2) + \frac{t_yy(t, x, y)}{dt}\right) damp(x, y)\end{matrix}\right]\)
#NBVAL_IGNORE_OUTPUT
# Partial ru for 1.2sec to plot the wavefield
op(dt=model.critical_dt, time_M=int(1000/model.critical_dt))Operator `Kernel` ran in 0.22 sPerformanceSummary([(PerfKey(name='section0', rank=None),
PerfEntry(time=0.19680499999999987, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),
(PerfKey(name='section1', rank=None),
PerfEntry(time=0.005080999999999991, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),
(PerfKey(name='section2', rank=None),
PerfEntry(time=0.005789000000000039, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),
(PerfKey(name='section3', rank=None),
PerfEntry(time=0.005510000000000015, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),
(PerfKey(name='section4', rank=None),
PerfEntry(time=0.005339000000000027, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])#NBVAL_SKIP
scale = .5*1e-3
plt_options_model = {'extent': [model.origin[0] , model.origin[0] + model.domain_size[0],
model.origin[1] + model.domain_size[1], model.origin[1]]}
fig, ax = plt.subplots(nrows=2, ncols=2, figsize=(15, 15))
ax[0, 0].imshow(np.transpose(v[0].data[0][slices]), vmin=-scale, vmax=scale, cmap="RdGy", **plt_options_model)
ax[0, 0].imshow(np.transpose(model.lam.data[slices]), vmin=2.5, vmax=15.0, cmap="jet", alpha=.5, **plt_options_model)
ax[0, 0].set_aspect('auto')
ax[0, 0].set_xlabel('X (m)', fontsize=20)
ax[0, 0].set_ylabel('Depth (m)', fontsize=20)
ax[0, 0].set_title(r"$v_{x}$", fontsize=20)
ax[0, 1].imshow(np.transpose(v[1].data[0][slices]), vmin=-scale, vmax=scale, cmap="RdGy", **plt_options_model)
ax[0, 1].imshow(np.transpose(model.lam.data[slices]), vmin=2.5, vmax=15.0, cmap="jet", alpha=.5, **plt_options_model)
ax[0, 1].set_aspect('auto')
ax[0, 1].set_title(r"$v_{z}$", fontsize=20)
ax[0, 1].set_xlabel('X (m)', fontsize=20)
ax[0, 1].set_ylabel('Depth (m)', fontsize=20)
ax[1, 0].imshow(np.transpose(tau[0,0].data[0][slices]+tau[1,1].data[0][slices]),
vmin=-10*scale, vmax=10*scale, cmap="RdGy", **plt_options_model)
ax[1, 0].imshow(np.transpose(model.lam.data[slices]), vmin=2.5, vmax=15.0, cmap="jet",
alpha=.5, **plt_options_model)
ax[1, 0].set_aspect('auto')
ax[1, 0].set_title(r"$\tau_{xx} + \tau_{zz}$", fontsize=20)
ax[1, 0].set_xlabel('X (m)', fontsize=20)
ax[1, 0].set_ylabel('Depth (m)', fontsize=20)
ax[1, 1].imshow(np.transpose(tau[0,1].data[0][slices]), vmin=-scale, vmax=scale, cmap="RdGy", **plt_options_model)
ax[1, 1].imshow(np.transpose(model.lam.data[slices]), vmin=2.5, vmax=15.0, cmap="jet", alpha=.5, **plt_options_model)
ax[1, 1].set_aspect('auto')
ax[1, 1].set_title(r"$\tau_{xy}$", fontsize=20)
ax[1, 1].set_xlabel('X (m)', fontsize=20)
ax[1, 1].set_ylabel('Depth (m)', fontsize=20)
plt.tight_layout()
model._physical_parameters{'b', 'damp', 'lam', 'mu'}#NBVAL_IGNORE_OUTPUT
# Full run for the data
op(dt=model.critical_dt, time_m=int(1000/model.critical_dt))Operator `Kernel` ran in 0.24 sPerformanceSummary([(PerfKey(name='section0', rank=None),
PerfEntry(time=0.20575800000000005, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),
(PerfKey(name='section1', rank=None),
PerfEntry(time=0.006976999999999987, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),
(PerfKey(name='section2', rank=None),
PerfEntry(time=0.007339000000000036, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),
(PerfKey(name='section3', rank=None),
PerfEntry(time=0.0069080000000000166, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),
(PerfKey(name='section4', rank=None),
PerfEntry(time=0.006813000000000026, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])# Data on a standard 2ms tim axis
rec_plot = rec.resample(num=1001)
rec2_plot = rec2.resample(num=1001)
rec3_plot = rec3.resample(num=1001)scale_for_plot = np.diag(np.linspace(1.0, 2.5, 1001)**2.0)#NBVAL_SKIP
# Pressure (txx + tzz) data at sea surface
extent = [rec_plot.coordinates.data[0, 0], rec_plot.coordinates.data[-1, 0], 1e-3*tn, t0]
aspect = rec_plot.coordinates.data[-1, 0]/(1e-3*tn)/.5
plt.figure(figsize=(15, 15))
plt.imshow(np.dot(scale_for_plot, rec_plot.data), vmin=-.01, vmax=.01, cmap="seismic",
interpolation='lanczos', extent=extent, aspect=aspect)
plt.ylabel("Time (s)", fontsize=20)
plt.xlabel("Receiver position (m)", fontsize=20)Text(0.5, 0, 'Receiver position (m)')
#NBVAL_SKIP
# OBC data of vx/vz
plt.figure(figsize=(15, 15))
plt.subplot(121)
plt.imshow(rec2_plot.data, vmin=-1e-3, vmax=1e-3, cmap="seismic",
interpolation='lanczos', extent=extent, aspect=aspect)
plt.ylabel("Time (s)", fontsize=20)
plt.xlabel("Receiver position (m)", fontsize=20)
plt.subplot(122)
plt.imshow(rec3_plot.data, vmin=-1e-3, vmax=1e-3, cmap="seismic",
interpolation='lanczos', extent=extent, aspect=aspect)
plt.ylabel("Time (s)", fontsize=20)
plt.xlabel("Receiver position (m)", fontsize=20)Text(0.5, 0, 'Receiver position (m)')
assert np.isclose(norm(rec), 22.68, atol=0, rtol=1e-3)
assert np.isclose(norm(rec2), 1.8299, atol=0, rtol=1e-3)
assert np.isclose(norm(rec3), 3.8671, atol=0, rtol=1e-3)Second order in time
Now that looks pretty! But let’s do it again with a 2nd order in time
so = 8
v2 = VectorTimeFunction(name='v2', grid=model.grid, space_order=so, time_order=2)
tau0 = TensorFunction(name='t0', grid=model.grid, space_order=so)
# The source injection term
src_xx = src.inject(field=tau0[0, 0], expr=src.dt)
src_zz = src.inject(field=tau0[1, 1], expr=src.dt)
s = model.grid.time_dim.spacing
# Second order elastic wave equation
pde_v2 = v2.dt2 - ro * div(tau0) + (1 - model.damp) * v2.dt
# Time update
u_v = Eq(v2.forward, solve(pde_v2, v2.forward))
# The stress equation isn't time dependent so we don't need solve.
u_t = Eq(tau0, model.damp * (l * diag(div(v2.forward)) + mu * (grad(v2.forward) + grad(v2.forward).transpose(inner=False))))
rec_term2 = rec2.interpolate(expr=v2[0])
rec_term2 += rec3.interpolate(expr=v2[1])
op = Operator([u_v] + [u_t] + src_xx + src_zz + rec_term2)#NBVAL_IGNORE_OUTPUT
# Partial ru for 1.2sec to plot the wavefield
op(dt=model.critical_dt, time_M=int(1000/model.critical_dt))Operator `Kernel` ran in 0.21 sPerformanceSummary([(PerfKey(name='section0', rank=None),
PerfEntry(time=0.1906600000000002, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),
(PerfKey(name='section1', rank=None),
PerfEntry(time=0.005354999999999988, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),
(PerfKey(name='section2', rank=None),
PerfEntry(time=0.005681000000000019, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),
(PerfKey(name='section3', rank=None),
PerfEntry(time=0.005547000000000023, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])#NBVAL_SKIP
scale = 1e-4
plt_options_model = {'extent': [model.origin[0] , model.origin[0] + model.domain_size[0],
model.origin[1] + model.domain_size[1], model.origin[1]]}
fig, ax = plt.subplots(nrows=2, ncols=2, figsize=(15, 15))
ax[0, 0].imshow(np.transpose(v2[0].data[0][slices]), vmin=-scale, vmax=scale, cmap="RdGy", **plt_options_model)
ax[0, 0].imshow(np.transpose(model.lam.data[slices]), vmin=2.5, vmax=15.0, cmap="jet", alpha=.5, **plt_options_model)
ax[0, 0].set_aspect('auto')
ax[0, 0].set_xlabel('X (m)', fontsize=20)
ax[0, 0].set_ylabel('Depth (m)', fontsize=20)
ax[0, 0].set_title(r"$v_{x}$", fontsize=20)
ax[0, 1].imshow(np.transpose(v2[1].data[0][slices]), vmin=-scale, vmax=scale, cmap="RdGy", **plt_options_model)
ax[0, 1].imshow(np.transpose(model.lam.data[slices]), vmin=2.5, vmax=15.0, cmap="jet", alpha=.5, **plt_options_model)
ax[0, 1].set_aspect('auto')
ax[0, 1].set_title(r"$v_{z}$", fontsize=20)
ax[0, 1].set_xlabel('X (m)', fontsize=20)
ax[0, 1].set_ylabel('Depth (m)', fontsize=20)
ax[1, 0].imshow(np.transpose(tau0[0,0].data[slices]+tau0[1,1].data[slices]),
vmin=-10*scale, vmax=10*scale, cmap="RdGy", **plt_options_model)
ax[1, 0].imshow(np.transpose(model.lam.data[slices]), vmin=2.5, vmax=15.0, cmap="jet",
alpha=.5, **plt_options_model)
ax[1, 0].set_aspect('auto')
ax[1, 0].set_title(r"$\tau_{xx} + \tau_{zz}$", fontsize=20)
ax[1, 0].set_xlabel('X (m)', fontsize=20)
ax[1, 0].set_ylabel('Depth (m)', fontsize=20)
ax[1, 1].imshow(np.transpose(tau0[0,1].data[slices]), vmin=-scale, vmax=scale, cmap="RdGy", **plt_options_model)
ax[1, 1].imshow(np.transpose(model.lam.data[slices]), vmin=2.5, vmax=15.0, cmap="jet", alpha=.5, **plt_options_model)
ax[1, 1].set_aspect('auto')
ax[1, 1].set_title(r"$\tau_{xy}$", fontsize=20)
ax[1, 1].set_xlabel('X (m)', fontsize=20)
ax[1, 1].set_ylabel('Depth (m)', fontsize=20)
plt.tight_layout()
#NBVAL_IGNORE_OUTPUT
op(dt=model.critical_dt, time_m=int(1000/model.critical_dt))Operator `Kernel` ran in 0.27 sPerformanceSummary([(PerfKey(name='section0', rank=None),
PerfEntry(time=0.22429000000000004, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),
(PerfKey(name='section1', rank=None),
PerfEntry(time=0.013601000000000016, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),
(PerfKey(name='section2', rank=None),
PerfEntry(time=0.012073000000000028, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),
(PerfKey(name='section3', rank=None),
PerfEntry(time=0.012027000000000043, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])rec2_plot2 = rec2.resample(num=1001)
rec3_plot2 = rec3.resample(num=1001)#NBVAL_SKIP
# OBC data of vx/vz
plt.figure(figsize=(15, 15))
plt.subplot(121)
plt.imshow(rec2_plot2.data, vmin=-1e-3, vmax=1e-3, cmap="seismic",
interpolation='lanczos', extent=extent, aspect=aspect)
plt.ylabel("Time (s)", fontsize=20)
plt.xlabel("Receiver position (m)", fontsize=20)
plt.subplot(122)
plt.imshow(rec3_plot2.data, vmin=-1e-3, vmax=1e-3, cmap="seismic",
interpolation='lanczos', extent=extent, aspect=aspect)
plt.ylabel("Time (s)", fontsize=20)
plt.xlabel("Receiver position (m)", fontsize=20)Text(0.5, 0, 'Receiver position (m)')
assert np.isclose(norm(rec2), .407183, atol=0, rtol=1e-3)
assert np.isclose(norm(rec3), .209852, atol=0, rtol=1e-3)Rotated staggered grid
Now let’s use the rotated staggered grid to avoid dispersion and instabilities. This method is usually used for TTI elastic modeling but we show it here for a simplified isotropic case
Simple example of RSFD
We first show the stencil generated by the RSFD method compared to the standard staggered grid method
x = model.grid.dimensions[0]
f = Function(name='f', grid=model.grid, space_order=2, staggered=NODE)
dfdx, dfdxrsfd = f.dx, f.dx45dfdx.evaluate\(\displaystyle \frac{f(x, y)}{h_{x}} - \frac{f(x - h_x, y)}{h_{x}}\)
dfdxrsfd.evaluate\(\displaystyle \left(- \frac{f(x - h_x, y - h_y)}{4 h_{x}} + \frac{f(x + h_x, y + h_y)}{4 h_{x}}\right) + \left(- \frac{f(x - h_x, y + h_y)}{4 h_{x}} + \frac{f(x + h_x, y - h_y)}{4 h_{x}}\right)\)
Elastic modeling
from devito import div45, grad45
all_node = [[NODE for _ in range(model.grid.dim)] for _ in range(model.grid.dim)]
all_vert = [model.grid.dimensions for _ in range(model.grid.dim)]
so = 8
v_rsfd = VectorTimeFunction(name='vr', grid=model.grid, space_order=so, time_order=1, staggered=all_vert)
tau_rsfd = TensorTimeFunction(name='tr', grid=model.grid, space_order=so, time_order=1, staggered=all_node)
# The source injection term
src_xx = src.inject(field=v_rsfd.forward.diagonal(), expr=s*src)
# First order elastic wave equation
pde_v = v_rsfd.dt - ro * div45(tau_rsfd)
pde_tau = tau_rsfd.dt - l * diag(div45(v_rsfd.forward)) - mu * (grad45(v_rsfd.forward) + grad45(v_rsfd.forward).transpose(inner=False))
# Time update
u_v = Eq(v_rsfd.forward, model.damp * solve(pde_v, v_rsfd.forward))
u_t = Eq(tau_rsfd.forward, model.damp * solve(pde_tau, tau_rsfd.forward))
# Receiver
rec_term = rec.interpolate(expr=tau_rsfd[0, 0] + tau_rsfd[1, 1])
rec_term += rec2.interpolate(expr=v_rsfd[1])
rec_term += rec3.interpolate(expr=v_rsfd[0])
op = Operator([u_v] + [u_t] + src_xx + rec_term)#NBVAL_IGNORE_OUTPUT
op(dt=model.critical_dt)Operator `Kernel` ran in 0.57 sPerformanceSummary([(PerfKey(name='section0', rank=None),
PerfEntry(time=0.5206050000000007, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),
(PerfKey(name='section1', rank=None),
PerfEntry(time=0.010398000000000041, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),
(PerfKey(name='section2', rank=None),
PerfEntry(time=0.012117999999999964, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),
(PerfKey(name='section3', rank=None),
PerfEntry(time=0.011993000000000005, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),
(PerfKey(name='section4', rank=None),
PerfEntry(time=0.011515999999999974, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])rec_plot2 = rec.resample(num=1001)
rec2_plot2 = rec2.resample(num=1001)
rec3_plot2 = rec3.resample(num=1001)#NBVAL_SKIP
# Pressure (txx + tzz) data at sea surface
extent = [rec_plot.coordinates.data[0, 0], rec_plot.coordinates.data[-1, 0], 1e-3*tn, t0]
aspect = rec_plot.coordinates.data[-1, 0]/(1e-3*tn)/.5
plt.figure(figsize=(15, 15))
plt.imshow(np.dot(scale_for_plot, rec_plot.data), vmin=-.01, vmax=.01, cmap="seismic",
interpolation='lanczos', extent=extent, aspect=aspect)
plt.ylabel("Time (s)", fontsize=20)
plt.xlabel("Receiver position (m)", fontsize=20)Text(0.5, 0, 'Receiver position (m)')
#NBVAL_SKIP
# OBC data of vx/vz
plt.figure(figsize=(15, 15))
plt.subplot(121)
plt.imshow(rec2_plot.data, vmin=-1e-3, vmax=1e-3, cmap="seismic",
interpolation='lanczos', extent=extent, aspect=aspect)
plt.ylabel("Time (s)", fontsize=20)
plt.xlabel("Receiver position (m)", fontsize=20)
plt.subplot(122)
plt.imshow(rec3_plot.data, vmin=-1e-3, vmax=1e-3, cmap="seismic",
interpolation='lanczos', extent=extent, aspect=aspect)
plt.ylabel("Time (s)", fontsize=20)
plt.xlabel("Receiver position (m)", fontsize=20)Text(0.5, 0, 'Receiver position (m)')
assert np.isclose(norm(rec), 29.83, atol=0, rtol=1e-3)
assert np.isclose(norm(rec2), 3.4437, atol=0, rtol=1e-3)
assert np.isclose(norm(rec3), 4.5632, atol=0, rtol=1e-3)