In [1]:
import sys
sys.path.append('../')
import numpy as np
from zephyr.backend import Eurus, StackedSimpleSource, KaiserSource, AnisotropicKaiserSource, MultiFreq
In [2]:
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import matplotlib
import matplotlib.animation as animation
import math
%matplotlib inline
from IPython.display import set_matplotlib_formats
set_matplotlib_formats('png')
matplotlib.rcParams['savefig.dpi'] = 150 # Change this to adjust figure size
In [3]:
class Source(object):
def __init__(self, sc):
self._x, self._z = np.mgrid[
0:sc['dx']*sc['nx']:sc['dx'],
0:sc['dz']*sc['nz']:sc['dz']
]
def __call__(self, x, z):
dist = np.sqrt((self._x - x)**2 + (self._z - z)**2)
srcterm = 1.*(dist == dist.min())
nz, nx = self._x.shape
return np.hstack([srcterm.T.ravel() / srcterm.sum(), np.zeros(nx*nz, dtype=np.complex128)])
In [4]:
# Geometry parameters
dx = 5.
dz = 5.
nx = 50
nz = 50
# Bulk parameters
velReal = 1000. * np.ones((nz,nx))
density = 1. * np.ones((nz,nx))
Q = np.inf
velImag = ((-1 * velReal) * 2)/Q
velocity = velReal + (1j * velImag)
# Anisotropy parameters
theta = 0 * np.ones((nz,nx))
epsilon = 0.2 * np.ones((nz,nx))
delta = 0.1 * np.ones((nz,nx))
# Other parameters
freq = 4.
freeSurf = [False, False, False, False]
nPML = 5
fDamp = 0.4
# Source and Modelling Parameters
sourceFreq = 32.
# determine the maximum length of the Model
xMax = (nx-1) * dx
zMax = (nz-1) * dz
LMax = np.max([xMax,zMax])
#calculate the lowest velocity
cMin = velocity.min()
#calculate the smallest wavelength, in meters
#lambdaMin = cMin / fMax
# Set the sampling interval and number of time domain samples
dt = 0.002
nTSamp = 256
#nt = tMax / dt
# Set the maximum modelled time and the damping factor
tMax = dt * nTSamp
tau = fDamp * tMax
# determine the frequency interval, max frequency and frequency range from the maximum modelled time
df= 1 / tMax
fMax = 1 / dt
nf = int(fMax / df)
fMin = fMax - ((nf - 1) * df)
freqRange = np.arange(fMin,fMax + df,df)
# Pack values into systemConfig dictionary
systemConfig = {
'nx': nx,
'nz': nz,
'dx': dx,
'dz': dz,
'c': velocity,
'rho': density,
'theta': theta,
'eps': epsilon,
'delta': delta,
'freq': freqRange[0],
'freqs': freqRange,
'freeSurf': freeSurf,
'nPML': nPML,
'cPML': 1e3,
'tau': tau,
'f_max': fMax,
'source_freq': sourceFreq,
'nt': nTSamp,
'nf': nf,
'df': df,
'nf': nf,
'disc': Eurus, # discretization
'parallel': True,
}
In [5]:
#keuper wavelet
excursions=1
m = (excursions + 2) / excursions
nsrc = int((1. / sourceFreq) / dt)
delta = excursions * np.pi * sourceFreq
tsrc = np.zeros(nsrc)
tempSource = np.zeros(nsrc)
for i in range(1,nsrc+1):
tsrc[i-1] = (i - 1.) * dt
tempSource[i-1] = (delta * np.cos(delta *tsrc[i-1])) - (delta *np.cos(m * delta *tsrc[i-1]))
timeSource = np.zeros(nTSamp)
timeSource[:len(tempSource)] = tempSource[:]
# plot the source
plt.plot(np.diff(timeSource), markersize=0.1, marker=',')
# perform forward fft to generate frequency domain source wavefields
freqSource=np.fft.fft(np.asarray(timeSource))
freqSource[-((nTSamp / 2) + 1):]=0
#plt.plot(abs(freqSource), markersize=0.1, marker=',')
In [6]:
uPSimple = np.zeros((nf,nz,nx),dtype=np.complex128)
uPKaiser = np.zeros((nf,nz,nx),dtype=np.complex128)
uPAniKaiser = np.zeros((nf,nz,nx),dtype=np.complex128)
#Ainv = Eurus(systemConfig)
for i in range(1,len(freqRange)):
#if (freqRange[i-1] <= sourceFreq):
if (freqRange[i-1] <= (sourceFreq * 4.)):
systemConfig['freq'] = freqRange[i-1]
Ainv = Eurus(systemConfig)
srcSimple = StackedSimpleSource(systemConfig)
srcAniKaiser = AnisotropicKaiserSource(systemConfig)
srcKaiser = KaiserSource(systemConfig)
qSimple = srcSimple(np.array([[125,125]]))
qKaiser = srcKaiser(np.array([[125,125]]))
qAniKaiser = srcAniKaiser(np.array([[125,125]]))
uPSimple[i-1] = (Ainv* (qSimple * (-1 * freqSource[i])))[:nx*nz].reshape((nz,nx))
uPKaiser[i-1] = (Ainv* (qKaiser * (-1 * freqSource[i])))[:nx*nz].reshape((nz,nx))
uPAniKaiser[i-1] = (Ainv* (qAniKaiser * (-1 * freqSource[i])))[:nx*nz].reshape((nz,nx))
else:
break
In [7]:
#Convert frequency domain wavefields to time domain using ifft
tempFreqSimple = np.zeros((nTSamp,nz,nx),dtype=np.complex128)
tempFreqKaiser = np.zeros((nTSamp,nz,nx),dtype=np.complex128)
tempFreqAniKaiser = np.zeros((nTSamp,nz,nx),dtype=np.complex128)
istart = 1
iend = math.ceil(nTSamp/2)+1
tempFreqSimple[istart:iend] = uPSimple[:iend-1]
tempFreqKaiser[istart:iend] = uPKaiser[:iend-1]
tempFreqAniKaiser[istart:iend] = uPAniKaiser[:iend-1]
tempFreqSimple[iend:] =np.flipud(uPSimple[:iend-2]).conjugate()
tempFreqKaiser[iend:] =np.flipud(uPKaiser[:iend-2]).conjugate()
tempFreqAniKaiser[iend:] =np.flipud(uPAniKaiser[:iend-2]).conjugate()
uTSimple = np.zeros((nTSamp,nz,nx),dtype=np.complex128)
uTSimple = np.fft.ifft(tempFreqSimple.conjugate(),axis=0)
uTKaiser = np.zeros((nTSamp,nz,nx),dtype=np.complex128)
uTKaiser = np.fft.ifft(tempFreqKaiser.conjugate(),axis=0)
uTAniKaiser = np.zeros((nTSamp,nz,nx),dtype=np.complex128)
uTAniKaiser = np.fft.ifft(tempFreqAniKaiser.conjugate(),axis=0)
In [8]:
clip = 50
tRef = 40
plotopts = {
'vmin': -clip,
'vmax': clip,
'extent': [0., dx * nx, dz * nz, 0.],
'cmap': cm.bwr,
}
fig = plt.figure()
# ax1 = fig.add_subplot(1,2,1)
# plt.imshow(uTSimple[tRef].reshape((nz,nx)).real, **plotopts)
# plt.title('Simple Source')
# plt.xlabel('X')
# plt.ylabel('Z')
# plt.colorbar()
ax1 = fig.add_subplot(1,2,1)
plt.imshow(uTKaiser[tRef].reshape((nz,nx)).real, **plotopts)
plt.title('Kaiser Source')
plt.xlabel('X')
plt.ylabel('Z')
#plt.colorbar()
clip = 50
plotopts.update({
'vmin': -clip,
'vmax': clip,
})
ax2 = fig.add_subplot(1,2,2)
plt.imshow(uTAniKaiser[tRef].reshape((nz,nx)).real, **plotopts)
plt.title('Anisotropic Kaiser Source')
plt.xlabel('X')
plt.ylabel('Z')
#plt.colorbar()
fig.tight_layout()
In [9]:
#plt.figure()
#plt.plot(uT[20].reshape((nz,nx)).real[:,50])
In [10]:
clip = 300
fRef = 30
plotopts = {
'vmin': -clip,
'vmax': clip,
'extent': [0., dx * nx, dz * nz, 0.],
'cmap': cm.bwr,
}
fig = plt.figure()
ax1 = fig.add_subplot(1,2,1)
plt.imshow(uPSimple[fRef].reshape((nz,nx)).real, **plotopts)
plt.title('Simple Source')
plt.xlabel('X')
plt.ylabel('Z')
plt.colorbar()
#ax2 = fig.add_subplot(1,3,2)
#plt.imshow(uPKaiser[tRef].reshape((nz,nx)).real, **plotopts)
#plt.title('Kaiser Source')
#plt.xlabel('X')
#plt.ylabel('Z')
#plt.colorbar()
clip = 50
plotopts.update({
'vmin': -clip,
'vmax': clip,
})
ax2 = fig.add_subplot(1,2,2)
plt.imshow(uPAniKaiser[fRef].reshape((nz,nx)).real, **plotopts)
plt.title('Anisotropic Kaiser Source')
plt.xlabel('X')
plt.ylabel('Z')
plt.colorbar()
fig.tight_layout()
In [11]:
#clip = 2e-2
#fig = plt.figure()
#ax = fig.add_subplot(1,1,1, aspect=1)
#plotopts = {
# 'extent': [0, nx*dx, nz*dz, 0],
# 'cmap': cm.bwr,
# 'vmin': -clip,
# 'vmax': clip,
#}
#im = plt.imshow(uT[0].reshape((nz,nx)).real, **plotopts)
#im = plt.imshow(uT[0].reshape((nz,nx)).real)
#def updatefig(j):
# im.set_array(uT[j].reshape((nz,nx)).real)
# return im,
#ani = animation.FuncAnimation(fig, updatefig)
#fig = plt.figure(figsize=(8,6), dpi=160)
#ani.save()
#plt.imshow(uT[200].reshape((nz,nx)).real, **plotopts)
#plt.title('Real Wavefield')
#plt.xlabel('X')
#plt.ylabel('Z')
#plt.show()
In [12]:
print nf
print df
print fMax
print sourceFreq
print dt
In [13]:
wx = (1. + (2*epsilon) +np.sqrt(1+(2*delta)))/(1 + epsilon + np.sqrt(1+(2*delta)))
wz = (1. + np.sqrt(1+(2*delta)))/(1 + epsilon + np.sqrt(1+(2*delta)))
print wx.min()
print wz.min()
In [14]:
print Ainv.K
print Ainv.Kavg