Brendan Smithyman | bsmithym@uwo.ca | January, 2015
This section of the IPython notebook imports the necessary dependencies for the demo. First, we import NumPy, which is the de facto standard numerical computing library for Python. In the next line, we import the SeisFDFDKernel class from zephyr.Kernel, which is the Zephyr submodule that contains the numerical solver.
Second, we import matplotlib, which is one of the most commonly-used plotting libraries for Python, and was originally designed to be somewhat similar to MatLab's plotting interface. We enable inline plotting, and adjust some defaults such that the rendering engine generates 300 DPI PNG images for us to look at.
In [1]:
import numpy as np
from zephyr.Kernel import SeisFDFDKernel
In [2]:
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import matplotlib
%matplotlib inline
from IPython.display import set_matplotlib_formats
set_matplotlib_formats('png')
matplotlib.rcParams['savefig.dpi'] = 300 # Change this to adjust figure size
# Plotting options
font = {
'family': 'Bitstream Vera Sans',
'weight': 'normal',
'size': 8,
}
matplotlib.rc('font', **font)
This section configures the base system for the numerical Helmholtz solver. We configure:
cellsize, nx, nz)c, Q, density)srcs and recs, which go into geom)freeSurf, nPML)These other parameters are of particular interest here, since the parameter freeSurf and the parameter nPML determine how waves interact with the boundaries of the model. In typical usage, freeSurf might be set to [False, False, False, False], meaning that free-surface conditions aren't applied on any of the boundaries. However, it might be convenient to enable free-surface conditions on the top boundary in diving-wave examples with buried sources. The nPML parameter controls the size of the PML (Perfectly Matched Layer) region on the sides of the model. Too small, and waves heading towards the boundaries may not be absorbed. Too large, and the computational costs may be prohibitive.
In [3]:
# ------------------------------------------------
# Model geometry
cellSize = 1 # m
freq = 2e2 # Hz
nx = 164 # count
nz = 264 # count
dims = (nz, nx)
# ------------------------------------------------
# Model properties
velocity = 2500 # m/s
vanom = 500 # m/s
density = 2700 # units of density
Q = np.inf # can be inf
cPert = np.zeros(dims)
cPert[(nz/2)-20:(nz/2)+20,(nx/2)-20:(nx/2)+20] = vanom
c = np.fliplr(np.ones(dims) * velocity)
c += np.fliplr(cPert)
rho = np.fliplr(np.ones((nz,nx)) * density)
# ------------------------------------------------
# Survey geometry
srcs = np.array([np.ones(101)*32, np.zeros(101), np.linspace(32, 232, 101)]).T
recs = np.array([np.ones(101)*132, np.zeros(101), np.linspace(32, 232, 101)]).T
nsrc = len(srcs)
nrec = len(recs)
recmode = 'fixed'
geom = {
'src': srcs,
'rec': recs,
'mode': 'fixed',
}
# ------------------------------------------------
# Other parameters
ky = 0
freeSurf = [False, False, False, False] # t r b l
nPML = 32
# Base configuration for all subproblems
systemConfig = {
'dx': cellSize, # m
'dz': cellSize, # m
'c': c, # m/s
'rho': rho, # density
'Q': Q, # can be inf
'nx': nx, # count
'nz': nz, # count
'freeSurf': freeSurf, # t r b l
'nPML': nPML,
'geom': geom,
'freq': freq,
'ky': 0,
}
In [4]:
sms = 4
rms = 0.5
def plotField(u):
clip = 0.1*abs(u).max()
plt.imshow(u.real, cmap=cm.bwr, vmin=-clip, vmax=clip)
def plotModel(v):
lclip = 2000
hclip = 3000
plt.imshow(v.real, cmap=cm.jet, vmin=lclip, vmax=hclip)
def plotGeometry():
srcpos = srcs[activeSource][::2]
recpos = recs[:,::2]
axistemp = plt.axis()
plt.plot(srcpos[0], srcpos[1], 'kx', markersize=sms)
plt.plot(recpos[:,0], recpos[:,1], 'kv', markersize=rms)
plt.axis(axistemp)
These examples model wavefields and data for a single source located near the top-left of the model.
In [10]:
activeSource = 0
In [11]:
systemConfig_PML = systemConfig.copy()
sp_PML = SeisFDFDKernel(systemConfig_PML)
u_PML, d_PML = sp_PML.forward(activeSource, False)
u_PML.shape = (nz+1, nx+1)
In [12]:
systemConfig_FIXED = systemConfig.copy()
systemConfig_FIXED.update({'nPML': 2})
sp_FIXED = SeisFDFDKernel(systemConfig_FIXED)
u_FIXED, d_FIXED = sp_FIXED.forward(activeSource, False)
u_FIXED.shape = (nz+1, nx+1)
In [13]:
fig = plt.figure()
ax1 = fig.add_subplot(1,3,1)
plotModel(c)
plotGeometry()
ax1.set_title('Velocity Model')
ax1.set_xlabel('X')
ax1.set_ylabel('Z')
ax2 = fig.add_subplot(1,3,2)
plotField(u_PML)
ax2.set_title('PML: %d cells'%systemConfig_PML['nPML'])
ax2.set_xlabel('X')
ax2.set_ylabel('Z')
ax3 = fig.add_subplot(1,3,3)
plotField(u_FIXED)
ax3.set_title('PML: %d cells'%systemConfig_FIXED['nPML'])
ax3.set_xlabel('X')
ax3.set_ylabel('Z')
fig.tight_layout()
In [14]:
norm = abs(d_PML).max()
fig = plt.figure()
ax1 = fig.add_axes([0.05, 0.05, 0.20, 0.90])
plt.plot(d_PML.real/norm, np.arange(nrec), 'b-', label='PML')
plt.plot(abs(d_PML)/norm, np.arange(nrec), 'y-')
plt.plot(-abs(d_PML)/norm, np.arange(nrec), 'y-')
ax1.set_ylabel('Receiver Number')
ax1.set_title('PML case')
ax2 = fig.add_axes([0.28, 0.05, 0.62, 0.90])
plt.plot(d_PML.real/norm, np.arange(nrec), 'b-', label='PML')
plt.plot(d_FIXED.real/norm, np.arange(nrec), 'r-', label='Fixed')
plt.tick_params(\
axis='y', # changes apply to the x-axis
labelleft='off') # labels along the bottom edge are off
ax2.set_title('Data comparison: PML case vs. Fixed case')
plt.legend()
Out[14]:
In [9]: