# Required imports:
import numpy as np
import sympy as sp
from devito import *
from examples.seismic.source import RickerSource, TimeAxis
from examples.seismic import ModelViscoelastic, plot_image
09 - Viscoelastic wave equation implementation on a staggered grid
This is a first attempt at implementing the viscoelastic wave equation as described in [1]. See also the FDELMODC implementation by Jan Thorbecke [2].
In the following example, a three dimensional toy problem will be introduced consisting of a single Ricker source located at (100, 50, 35) in a 200 m \(\times\) 100 m \(\times\) 100 m domain.
The model domain is now constructed. It consists of an upper layer of water, 50 m in depth, and a lower rock layer separated by a 4 m thick sediment layer.
# Domain size:
= (200., 100., 100.) # 200 x 100 x 100 m domain
extent = 1.0 # Desired grid spacing
h = (int(extent[0]/h+1), int(extent[1]/h+1), int(extent[2]/h+1))
shape
# Model physical parameters:
= np.zeros(shape)
vp = np.zeros(shape)
qp = np.zeros(shape)
vs = np.zeros(shape)
qs = np.zeros(shape)
rho
# Set up three horizontally separated layers:
int(0.5*shape[2])+1] = 1.52
vp[:,:,:int(0.5*shape[2])+1] = 10000.
qp[:,:,:int(0.5*shape[2])+1] = 0.
vs[:,:,:int(0.5*shape[2])+1] = 0.
qs[:,:,:int(0.5*shape[2])+1] = 1.05
rho[:,:,:
int(0.5*shape[2])+1:int(0.5*shape[2])+1+int(4/h)] = 1.6
vp[:,:,int(0.5*shape[2])+1:int(0.5*shape[2])+1+int(4/h)] = 40.
qp[:,:,int(0.5*shape[2])+1:int(0.5*shape[2])+1+int(4/h)] = 0.4
vs[:,:,int(0.5*shape[2])+1:int(0.5*shape[2])+1+int(4/h)] = 30.
qs[:,:,int(0.5*shape[2])+1:int(0.5*shape[2])+1+int(4/h)] = 1.3
rho[:,:,
int(0.5*shape[2])+1+int(4/h):] = 2.2
vp[:,:,int(0.5*shape[2])+1+int(4/h):] = 100.
qp[:,:,int(0.5*shape[2])+1+int(4/h):] = 1.2
vs[:,:,int(0.5*shape[2])+1+int(4/h):] = 70.
qs[:,:,int(0.5*shape[2])+1+int(4/h):] = 2. rho[:,:,
Now create a Devito vsicoelastic model generating an appropriate computational grid along with absorbing boundary layers:
# Create model
= (0, 0, 0)
origin = (h, h, h)
spacing = 4 # FD space order (Note that the time order is by default 1).
so = 20 # Number of absorbing boundary layers cells
nbl = ModelViscoelastic(space_order=so, vp=vp, qp=qp, vs=vs, qs=qs,
model =1/rho, origin=origin, shape=shape, spacing=spacing,
b=nbl) nbl
Operator `initdamp` ran in 0.01 s
# As pointed out in Thorbecke's implementation and documentation, the viscoelastic wave euqation is
# not always stable with the standard elastic CFL condition. We enforce a smaller critical dt here
# to ensure the stability.
= .9 model.dt_scale
The source frequency is now set along with the required model parameters:
# Source freq. in MHz (note that the source is defined below):
= 0.12
f0
# Thorbecke's parameter notation
= model.lam
l = model.mu
mu = model.b
ro
= 1.0/(l + 2*mu)
k = l + 2*mu
pi
= (sp.sqrt(1.+1./model.qp**2)-1./model.qp)/f0
t_s = 1./(f0**2*t_s)
t_ep = (1.+f0*model.qs*t_s)/(f0*model.qs-f0**2*t_s) t_es
# Time step in ms and time range:
= 0., 30.
t0, tn = model.critical_dt
dt = TimeAxis(start=t0, stop=tn, step=dt) time_range
Generate Devito time functions for the velocity, stress and memory variables appearing in the viscoelastic model equations. By default, the initial data of each field will be set to zero.
# PDE fn's:
= model.grid.dimensions
x, y, z = model.damp
damp
# Staggered grid setup:
# Velocity:
= VectorTimeFunction(name="v", grid=model.grid, time_order=1, space_order=so)
v
# Stress:
= TensorTimeFunction(name='t', grid=model.grid, space_order=so, time_order=1)
tau
# Memory variable:
= TensorTimeFunction(name='r', grid=model.grid, space_order=so, time_order=1)
r
= model.grid.stepping_dim.spacing # Symbolic representation of the model grid spacing s
And now the source and PDE’s are constructed:
# Source
= RickerSource(name='src', grid=model.grid, f0=f0, time_range=time_range)
src = np.array([100., 50., 35.])
src.coordinates.data[:]
# The source injection term
= src.inject(field=tau[0, 0].forward, expr=src*s)
src_xx = src.inject(field=tau[1, 1].forward, expr=src*s)
src_yy = src.inject(field=tau[2, 2].forward, expr=src*s)
src_zz
# Particle velocity
= v.dt - ro * div(tau)
pde_v = Eq(v.forward, model.damp * solve(pde_v, v.forward))
u_v # Strain
= grad(v.forward) + grad(v.forward).transpose(inner=False)
e
# Stress equations
= tau.dt - r.forward - l * t_ep / t_s * diag(div(v.forward)) - mu * t_es / t_s * e
pde_tau = Eq(tau.forward, model.damp * solve(pde_tau, tau.forward))
u_t
# Memory variable equations:
= r.dt + 1 / t_s * (r + l * (t_ep/t_s-1) * diag(div(v.forward)) + mu * (t_es/t_s-1) * e)
pde_r = Eq(r.forward, damp * solve(pde_r, r.forward)) u_r
We now create and then run the operator:
# Create the operator:
= Operator([u_v, u_r, u_t] + src_xx + src_yy + src_zz,
op =model.spacing_map) subs
#NBVAL_IGNORE_OUTPUT
# Execute the operator:
=dt) op(dt
Operator `Kernel` ran in 9.82 s
PerformanceSummary([(PerfKey(name='section0', rank=None),
PerfEntry(time=0.003796, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),
(PerfKey(name='section1', rank=None),
PerfEntry(time=9.803400999999997, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),
(PerfKey(name='section2', rank=None),
PerfEntry(time=0.007545000000000001, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])
Before plotting some results, let us first look at the shape of the data stored in one of our time functions:
0].data.shape v[
(2, 241, 141, 141)
Since our functions are first order in time, the time dimension is of length 2. The spatial extent of the data includes the absorbing boundary layers in each dimension (i.e. each spatial dimension is padded by 20 grid points to the left and to the right).
The total number of instances in time considered is obtained from:
time_range.num
158
Hence 223 time steps were executed. Thus the final time step will be stored in index given by:
2) np.mod(time_range.num,
0
Now, let us plot some 2D slices of the fields vx
and szz
at the final time step:
#NBVAL_SKIP
# Mid-points:
= int(0.5*(v[0].data.shape[1]-1))+1
mid_x = int(0.5*(v[0].data.shape[2]-1))+1
mid_y
# Plot some selected results:
0].data[1, :, mid_y, :], cmap="seismic")
plot_image(v[0].data[1, mid_x, :, :], cmap="seismic")
plot_image(v[
2, 2].data[1, :, mid_y, :], cmap="seismic")
plot_image(tau[2, 2].data[1, mid_x, :, :], cmap="seismic") plot_image(tau[
#NBVAL_IGNORE_OUTPUT
assert np.isclose(norm(v[0]), 0.102959, atol=1e-4, rtol=0)
References
[1] Johan O. A. Roberston, et.al. (1994). “Viscoelatic finite-difference modeling” GEOPHYSICS, 59(9), 1444-1456.
[2] https://janth.home.xs4all.nl/Software/fdelmodcManual.pdf