In [1]:
from SimPEG import *
from scipy.constants import mu_0, epsilon_0
%pylab inline
In [42]:
from SimPEG.Utils import mkvc, sdiag, sdInv, speye
In [43]:
cs = 1.0*1e-2
hx = np.ones(200)*cs
hy = np.ones(200)*cs
mesh = Mesh.TensorMesh([hx, hy], 'CC')
In [46]:
data[mesh, mesh] = 1.
$ \nabla_s \times \vec{E} = -\imath\omega\mu^{'}\vec{H}$
$ \nabla_s \times \vec{H} -\imath\omega\epsilon^{'}\vec{E} = 0 $
where
$\mu^{'} = \mu + \frac{\sigma_0^{*}}{\imath\omega}$
$\epsilon^{'} = \epsilon + \frac{\sigma_0}{\imath\omega}$
$\nabla_s = \hat{u}_x\frac{1}{s_x}\frac{\partial}{\partial x}+\hat{u}_y\frac{1}{s_y}\frac{\partial}{\partial y} +\hat{u}_z\frac{1}{s_z}\frac{\partial}{\partial z}$
where
$s_x (x) = s_{x0}(x)(1+\frac{\sigma_x(x)}{\imath\omega\epsilon})$
$s_y (y) = s_{y0}(y)(1+\frac{\sigma_y(y)}{\imath\omega\epsilon})$
$s_z (z) = s_{z0}(z)(1+\frac{\sigma_z(z)}{\imath\omega\epsilon})$
$-\frac{\partial ex^{n}}{\partial y} = -s{y0}[\mu\triangle t^{-1}(h{zy}^{n+1/2}-h{zy}^{n-1/2})
+ (\sigma^{*}_0+\sigma^{*}_y)h_{zy}^{n-1/2}
+ \mu^{-1}\sigma^{*}_0\sigma^{*}_yh_{zy}^{I \ n}
]$
$\frac{\partial ey^{n}}{\partial x} = -s{x0}[\mu\triangle t^{-1}(h{zx}^{n+1/2}-h{zx}^{n-1/2})
+ (\sigma^{*}_0+\sigma^{*}_x)h_{zx}^{n-1/2}
+ \mu^{-1}\sigma_0^{*}\sigma^{*}_xh_{zx}^{I \ n}
]$
where $h_{zx}^{I \ n} = h_{zx}^{I \ n-1}+\triangle t h_{zx}^{n-1/2}$ and $h_{zy}^{I \ n} = h_{zy}^{I \ n-1}+\triangle t h_{zy}^{n-1/2}$
$-\frac{\partial hz^{n+1/2}}{\partial x} = s{x0}[\epsilon\triangle t^{-1}(e{y}^{n+1}-e{y}^{n})
+ (\sigma_0+\sigma_x)e_{y}^{n-1/2}
+ \epsilon^{-1}\sigma_0\sigma_xe_{y}^{I \ n+1/2}
]$
$\frac{\partial hz^{n+1/2}}{\partial y} = s{y0}[\epsilon\triangle t^{-1}(e{x}^{n+1}-e{x}^{n})
+ (\sigma_0+\sigma_y)e_{x}^{n-1/2}
+ \epsilon^{-1}\sigma_0\sigma_ye_{x}^{I \ n+1/2}
]$
where $e_y^{I \ n+1/2} = e_{y}^{I \ n-1/2}+\triangle t e_{y}^{n}$ and $e_x^{I \ n+1/2} = e_{x}^{I \ n-1/2}+\triangle t e_{x}^{n}$
$\mathbf{Curl}^{vec}\mathbf{e}^{n} = -\triangle t^{-1}\mathbf{S}_{\mu}(\mathbf{h}_d^{n+1/2}-\mathbf{h}_d^{n-1/2}) -\mathbf{S}_{\sigma^{*}}\mathbf{h}_d^{n-1/2} - \mathbf{S}_{\mu^{-1}\sigma^{*}}\mathbf{h}_d^{I \ n} $
$\mathbf{h}_d^{n+1/2} = \mathbf{h}_d^{n-1/2} - \triangle t\mathbf{S}^{-1}_{\mu}[\mathbf{Curl}^{vec}\mathbf{e}^{n} +\mathbf{S}_{\sigma^{*}}\mathbf{h}_d^{n-1/2} + \mathbf{S}_{\mu^{-1}\sigma^{*}}\mathbf{h}_d^{I \ n}] $
where
$ \mathbf{h}_d = \begin{bmatrix} \mathbf{h_{zy}} & \\[0.3em] \mathbf{h_{zx}} \end{bmatrix} $ , $ \mathbf{S}_{\mu} = \begin{bmatrix} \mathbf{diag}(\mu s_{y0}) & 0 \\[0.3em] 0 & \mathbf{diag}(\mu s_{x0}) \end{bmatrix} $
$ \mathbf{S}_{\sigma^{*}} = \begin{bmatrix} \mathbf{diag}((\sigma_0^{*}+\sigma_y^{*}) s_{y0}) & 0 \\[0.3em] 0 & \mathbf{diag}((\sigma_0^{*}+\sigma_x^{*})s_{x0}) \end{bmatrix} $ , and $ \mathbf{S}_{\mu^{-1}\sigma^{*}} = \begin{bmatrix} \mathbf{diag}(\mu^{-1}(\sigma_0^{*}\sigma_y^{*}s_{y0}) & 0 \\[0.3em] 0 & \mathbf{diag}(\mu^{-1}(\sigma_0^{*}\sigma_x^{*}s_{x0}) \end{bmatrix} $
$\mathbf{Curl}^T\mathbf{h} = \triangle t^{-1} \mathbf{M}^e_{s\epsilon}(\mathbf{e}^{n+1}-\mathbf{e}^{n}) +\mathbf{M}^e_{s\sigma}\mathbf{e}^{n} + \mathbf{M}^e_{s\epsilon^{-1}\sigma} \mathbf{e}^{I \ n+1/2} $
$ \mathbf{e}^{n+1} = \mathbf{e}^{n} + \triangle t (\mathbf{M}^e_{s\epsilon})^{-1}[\mathbf{Curl}^T\mathbf{h}
- \mathbf{M}^e_{s\sigma}\mathbf{e}^{n} - \mathbf{M}^e_{s\epsilon^{-1}\sigma} \mathbf{e}^{I \ n+1/2}]
$
where
$ \mathbf{M}^e_{s\epsilon} = \mathbf{diag} (\mathbf{A}^{e \ T}_{c \ vec} \begin{bmatrix} \epsilon s_{y0} \\[0.3em] \epsilon s_{x0} \end{bmatrix} ) $ , $ \mathbf{M}^e_{s\sigma} = \mathbf{diag} (\mathbf{A}^{e \ T}_{c \ vec} \begin{bmatrix} (\sigma_0+\sigma_y) s_{y0} \\[0.3em] (\sigma_0+\sigma_x) s_{x0} \end{bmatrix} ) $
$ \mathbf{M}^e_{s\epsilon^{-1}\sigma} = \mathbf{diag} (\mathbf{A}^{e \ T}_{c \ vec} \begin{bmatrix} \epsilon^{-1} \sigma_0\sigma_y s_{y0} \\[0.3em] \epsilon^{-1} \sigma_0\sigma_x s_{x0} \end{bmatrix} ) $ and $ \mathbf{h} = [\mathbf{I}^{cc}, \mathbf{I}^{cc}] \begin{bmatrix} \mathbf{h_{zy}} & \\[0.3em] \mathbf{h_{zx}} \end{bmatrix} = \mathbf{h_{zy}}+\mathbf{h_{zx}} $
$s_{x0}(x) = 1+s_m(\frac{x}{\delta_x})^2$
$\sigma_x(x) = \sigma_m sin^2(\frac{\pi x}{2\delta_x})$
$\sigma_m = \frac{\epsilon c \delta^{-1}}{1+s_m(1/3+2/\pi^2)}log(R_{th})$
$ \mathbf{h}_d^{I \ n} = \mathbf{h}_d^{I \ n-1} + \triangle t \mathbf{h}_d^{n-1/2} $
$\mathbf{h}_d^{n+1/2} = \mathbf{h}_d^{n-1/2} - \triangle t\mathbf{S}^{-1}_{\mu}[\mathbf{Curl}^{vec}\mathbf{e}^{n} +\mathbf{S}_{\sigma^{*}}\mathbf{h}_d^{n-1/2} + \mathbf{S}_{\mu^{-1}\sigma^{*}}\mathbf{h}_d^{I \ n}] $
$ \mathbf{e}^{I \ n+1/2} = \mathbf{e}^{I \ n-1/2} + \triangle t \mathbf{e}^{n} $
$ \mathbf{e}^{n+1} = \mathbf{e}^{n} + \triangle t (\mathbf{M}^e_{s\epsilon})^{-1}[\mathbf{Curl}^T\mathbf{h}
- \mathbf{M}^e_{s\sigma}\mathbf{e}^{n} - \mathbf{M}^e_{s\epsilon^{-1}\sigma} \mathbf{e}^{I \ n+1/2}]
$
In [4]:
npad = 20
ax = mesh.vectorCCx[-npad]
ay = mesh.vectorCCy[-npad]
indy = np.logical_or(mesh.gridCC[:,1]<=-ay, mesh.gridCC[:,1]>=ay)
indx = np.logical_or(mesh.gridCC[:,0]<=-ax, mesh.gridCC[:,0]>=ax)
tempx = zeros_like(mesh.gridCC[:,0])
tempx[indx] = (abs(mesh.gridCC[:,0][indx])-ax)**2
tempx[indx] = tempx[indx]-tempx[indx].min()
tempx[indx] = tempx[indx]/tempx[indx].max()
tempy = zeros_like(mesh.gridCC[:,1])
tempy[indy] = (abs(mesh.gridCC[:,1][indy])-ay)**2
tempy[indy] = tempy[indy]-tempy[indy].min()
tempy[indy] = tempy[indy]/tempy[indy].max()
tempx[tempx>1.] = 1.
tempy[tempy>1.] = 1.
In [5]:
Lx = mesh.hx[-npad:].sum()
Ly = mesh.hy[-npad:].sum()
In [6]:
sighalf = 1e-3
mu = mu_0*np.ones(mesh.nC)
epsilon = epsilon_0*np.ones(mesh.nC)*1.
epsilon[mesh.gridCC[:,1]<0.5] = epsilon_0*2.
sig0 = sighalf*np.ones(mesh.nC)
sigs = 0.
c = 1/np.sqrt(mu*epsilon)
In [14]:
fmain = 3e9
lamda = c.min()/fmain
sm = lamda/(2*hx.min())-1
sm = 3.
Rth = 1e-8
sigm = -(epsilon.max()*c.max()/(0.5*(Lx+Ly)))/(1.+sm*(1./3+2./(np.pi**2)))*np.log(Rth)
print ('>> sm: %5.2e, lamda: %5.2e, sigm: %5.2e') % (sm, lamda, sigm)
In [15]:
sigx = sigm*np.sin(0.5*np.pi*np.sqrt(tempx))**2
sigy = sigm*np.sin(0.5*np.pi*np.sqrt(tempy))**2
sx0 = 1.+sm*tempx
sy0 = 1.+sm*tempy
In [16]:
plt.plot(mesh.vectorCCx, sigx.reshape((mesh.nCx, mesh.nCy), order='F')[:,0])
plt.plot(mesh.vectorCCx, sx0.reshape((mesh.nCx, mesh.nCy), order='F')[:,0])
plt.plot(mesh.vectorCCx, (sigx*sx0).reshape((mesh.nCx, mesh.nCy), order='F')[:,0])
plt.legend(('$\sigma_x$', '$s_{x0}$', '$\sigma_xs_{x0}$'), fontsize = 16)
plt.xlabel('x')
plt.grid(True)
In [17]:
Smu = sp.block_diag([sdiag(mu*sy0), sdiag(mu*sx0)])
SmuI = sp.block_diag([sdiag(1./(mu*sy0)), sdiag(1./(mu*sx0))])
Smuisig = sp.block_diag([sdiag(sigs*sigy*(1./epsilon)*sy0), sdiag(sigs*sigx*(1./epsilon)*sx0)])
Ssig = sp.block_diag([sdiag((sigy*mu/epsilon+sigs)*sy0), sdiag((sigx*mu/epsilon+sigs)*sx0)])
Mesepsisig = sdiag(mesh.aveE2CCV.T*np.r_[1./epsilon*sig0*sigy*sy0, 1./epsilon*sig0*sigx*sx0])
Messig = sdiag(mesh.aveE2CCV.T*np.r_[(sig0+sigy)*sy0, (sig0+sigx)*sx0])
Meseps = sdiag(mesh.aveE2CCV.T*np.r_[epsilon*sy0, epsilon*sx0])
MesepsI = sdInv(Meseps)
Icc = sp.hstack((speye(mesh.nC), speye(mesh.nC)))
curl = mesh.edgeCurl
curlvec = sp.block_diag((curl[:,:mesh.nEx], curl[:,mesh.nEx:]))
In [18]:
fig, ax = plt.subplots(1,3, figsize = (15, 6))
dat = mesh.plotImage(sigx, ax = ax[0])
plt.colorbar(dat[0], ax = ax[0], orientation='horizontal')
dat = mesh.plotImage(sigy, ax = ax[1])
plt.colorbar(dat[0], ax = ax[1], orientation='horizontal')
dat = mesh.plotImage(epsilon, ax = ax[2])
plt.colorbar(dat[0], ax = ax[2], orientation='horizontal')
ax[0].set_title("$\sigma_x$", fontsize = 18)
ax[1].set_title("$\sigma_y$", fontsize = 18)
ax[2].set_title("$\epsilon$", fontsize = 18)
Out[18]:
In [19]:
def ricker(fpeak, t, tlag):
"""
Generating Ricker Wavelet
.. math ::
"""
# return (1-2*np.pi**2*fpeak**2*(t-tlag)**2)*np.exp(-np.pi**2*fpeak**2*(t-tlag)**2)
return np.exp(-2*fpeak**2*(t-tlag)**2)*np.cos(np.pi*fpeak*(t-tlag))
In [20]:
dt = 1e-11
fmain = 3e9
time = np.arange(700)*dt
wave = ricker(fmain, time, 50*dt)
In [21]:
plt.plot(time, wave, '.-')
Out[21]:
In [26]:
# Put Jx source
# txind = Utils.closestPoints(mesh, [0., 0.], gridLoc='Ex')
# q = np.zeros(mesh.nE)
# q[txind] = 1.
# q = Utils.sdiag(1/mesh.edge)*q
# Put Mz source
txind = Utils.closestPoints(mesh, [0., 0.], gridLoc='CC')
q = np.zeros(mesh.nC)
q[txind] = 1.
Px = mesh.getInterpolationMat(mesh.gridCC, 'Ex')
Py = mesh.getInterpolationMat(mesh.gridCC, 'Ey')
$ \mathbf{h}_d^{I \ n} = \mathbf{h}_d^{I \ n-1} + \triangle t \mathbf{h}_d^{n-1/2} $
$\mathbf{h}_d^{n+1/2} = \mathbf{h}_d^{n-1/2} - \triangle t\mathbf{S}^{-1}_{\mu}[\mathbf{Curl}^{vec}\mathbf{e}^{n} +\mathbf{S}_{\sigma^{*}}\mathbf{h}_d^{n-1/2} + \mathbf{S}_{\mu^{-1}\sigma^{*}}\mathbf{h}_d^{I \ n}] $
$ \mathbf{e}^{I \ n+1/2} = \mathbf{e}^{I \ n-1/2} + \triangle t \mathbf{e}^{n} $
$ \mathbf{e}^{n+1} = \mathbf{e}^{n} + \triangle t (\mathbf{M}^e_{s\epsilon})^{-1}[\mathbf{Curl}^T\mathbf{h}
- \mathbf{M}^e_{s\sigma}\mathbf{e}^{n} - \mathbf{M}^e_{s\epsilon^{-1}\sigma} \mathbf{e}^{I \ n+1/2} + \mathbf{j}_s]
$
In [34]:
hd0 = np.zeros(2*mesh.nC)
hd1 = np.zeros(2*mesh.nC)
hId0 = np.zeros(2*mesh.nC)
hId1 = np.zeros(2*mesh.nC)
e0 = np.zeros(mesh.nE)
e1 = np.zeros(mesh.nE)
eI0 = np.zeros(mesh.nE)
eI1 = np.zeros(mesh.nE)
h = np.zeros((mesh.nC, time.size))
ex = np.zeros((mesh.nC, time.size))
ey = np.zeros((mesh.nC, time.size))
for i in range (time.size):
js = q*wave[i]
eI0 = eI1.copy()
eI1 = eI0 + dt*e0
# e1 = e0 + MesepsI*dt*(curl.T*(Icc*hd1) - Messig*e0 - Mesepsisig*eI1-js)
e1 = e0 + MesepsI*dt*(curl.T*(Icc*hd1) - Messig*e0 - Mesepsisig*eI1)
e0 = e1.copy()
ex[:,i] = Px*e1
ey[:,i] = Py*e1
hId0 = hId1.copy()
hId1 = hId0 + dt*hd0
hd1 = hd0 - SmuI*dt*(curlvec*e0+Ssig*hd0+Smuisig*hId1+np.r_[js, js]*0.5)
# hd1 = hd0 - SmuI*dt*(curlvec*e0+Ssig*hd0+Smuisig*hId1)
hd0 = hd1.copy()
h[:,i] = Icc*hd1
In [40]:
icount = 600
extent = [mesh.vectorCCx.min(), mesh.vectorCCx.max(), mesh.vectorCCy.min(), mesh.vectorCCy.max()]
plt.imshow(np.flipud(h[:,icount].reshape((mesh.nCx, mesh.nCy), order = 'F').T), cmap = 'RdBu', extent=extent)
# plt.imshow(np.flipud(c.reshape((mesh.nCx, mesh.nCy), order = 'F').T), cmap = 'jet', extent=extent, alpha=0.2)
Out[40]:
In [33]:
from JSAnimation import IPython_display
from matplotlib import animation
extent = [mesh.vectorCCx.min(), mesh.vectorCCx.max(), mesh.vectorCCy.min(), mesh.vectorCCy.max()]
fig, ax = plt.subplots(1,3, figsize = (24, 8))
for i in range(3):
ax[i].set_xlabel('x (m)', fontsize = 16)
ax[i].set_ylabel('y (m)', fontsize = 16)
ax[i].set_xlim(extent[:2])
ax[i].set_ylim(extent[2:])
ax[0].set_title('$h_z$', fontsize = 20)
ax[1].set_title('$e_x$', fontsize = 20)
ax[2].set_title('$e_y$', fontsize = 20)
nskip = 40
def animate(i_id):
icount = i_id*nskip
frame1 = ax[0].imshow(np.flipud(h[:,icount].reshape((mesh.nCx, mesh.nCy), order = 'F').T), cmap = 'RdBu', extent=extent)
frame2 = ax[1].imshow(np.flipud(ex[:,icount].reshape((mesh.nCx, mesh.nCy), order = 'F').T), cmap = 'RdBu', extent=extent)
frame3 = ax[2].imshow(np.flipud(ey[:,icount].reshape((mesh.nCx, mesh.nCy), order = 'F').T), cmap = 'RdBu', extent=extent)
for i in range(3):
ax[i].imshow(np.flipud(c.reshape((mesh.nCx, mesh.nCy), order = 'F').T), cmap = 'jet', extent=extent, alpha=0.2)
return frame1, frame2, frame3
animation.FuncAnimation(fig, animate, frames=16, interval=40, blit=True)
Out[33]:
In [ ]: