notebook.community

Edit and run 





In [1]:



    


from SimPEG import *
from scipy.constants import mu_0, epsilon_0
%pylab inline



    


















    






Populating the interactive namespace from numpy and matplotlib

















In [3]:



    


from SimPEG.Utils import mkvc, sdiag, sdInv, speye



    











In [5]:



    


cs = 1.0*1e-2
hx = np.ones(200)*cs
hy = np.ones(200)*cs
mesh = Mesh.TensorMesh([hx, hy], 'CC')



    











In [6]:



    


data[mesh, mesh] = 1.

2D Geo Penetrating Radar (GPR) modeling: $TM^z$ mode with Generalized PML (GPML) boundary

Generalized Perfectly Mached Layer for the Absorption of Propagating and Evanescent Waves in Losless and Lossy Media (Fand and Zhonghua, 1996)

Maxwell's equations in time domain:

$ \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})$

In time domain assuming $TE^z$ mode ($h_x$, $h_y$, $e_z$):

$-\frac{\partial hx^{n}}{\partial y} = s{y0}[\epsilon\triangle t^{-1}(e{zy}^{n+1/2}-e{zy}^{n-1/2})

                             + (\sigma_0+\sigma_y)e_{zy}^{n-1/2}
                             + \epsilon^{-1}\sigma_0\sigma_ye_{zy}^{I \ n}

]$

$\frac{\partial hy^{n}}{\partial x} = s{x0}[\epsilon\triangle t^{-1}(e{zx}^{n+1/2}-e{zx}^{n-1/2})

                             + (\sigma_0+\sigma_x)e_{zx}^{n-1/2}
                             + \epsilon^{-1}\sigma_0\sigma_xe_{zx}^{I \ n}

]$

where $e_{zx}^{I \ n} = e_{zx}^{I \ n-1}+\triangle t e_{zx}^{n-1/2}$ and $e_{zy}^{I \ n} = e_{zy}^{I \ n-1}+\triangle t e_{zy}^{n-1/2}$

$\mathbf{Curl}^{vec}\mathbf{h}^{n} = \triangle t^{-1}\mathbf{S}_{\epsilon}(\mathbf{e}_d^{n+1/2}-\mathbf{e}_d^{n-1/2}) +\mathbf{S}_{\sigma}\mathbf{e}_d^{n-1/2} + \mathbf{S}_{\epsilon^{-1}\sigma}\mathbf{e}_d^{I \ n} $

$ \mathbf{e}_d = \begin{bmatrix} \mathbf{e_{zy}} & \\[0.3em] \mathbf{e_{zx}} \end{bmatrix} $ , $ \mathbf{e} = [\mathbf{I}^{cc}, \mathbf{I}^{cc}] \begin{bmatrix} \mathbf{e_{zy}} & \\[0.3em] \mathbf{e_{zx}} \end{bmatrix} = \mathbf{e_{zy}}+\mathbf{e_{zx}} $ and $ \mathbf{S}_{\epsilon} = \begin{bmatrix} \mathbf{diag}(\epsilon s_{y0}) & 0 \\[0.3em] 0 & \mathbf{diag}(\epsilon 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}_{\epsilon^{-1}\sigma} = \begin{bmatrix} \mathbf{diag}(\epsilon^{-1}(\sigma_0\sigma_ys_{y0}) & 0 \\[0.3em] 0 & \mathbf{diag}(\epsilon^{-1}(\sigma_0\sigma_xs_{x0}) \end{bmatrix} $





In [62]:



    


# Seps = sp.block_diag([sdiag(epsilon*sy0), sdiag(epsilon*sx0)])
# SepsI = sp.block_diag([sdiag(1./(epsilon*sy0)), sdiag(1./(epsilon*sx0))])
# Sepsisig = sp.block_diag([sdiag(sig0*sigy*(1./epsilon)*sy0), sdiag(sig0*sigx*(1./epsilon)*sx0)])
# Ssig = sp.block_diag([sdiag((sigy+sig0)*sy0), sdiag((sigx+sig0)*sx0)])

$\frac{\partial ez^{n+1/2}}{\partial y} = -s{y0}[\mu\triangle t^{-1}(h{x}^{n+1}-h{x}^{n})

                             + (\sigma_0^{*}+\sigma_y^{*})h_{x}^{n-1/2}
                             + \mu^{-1}\sigma_0^{*}\sigma_y^{*}h_{x}^{I \ n+1/2}

]$

$-\frac{\partial ez^{n+1/2}}{\partial x} = -s{x0}[\mu\triangle t^{-1}(h{y}^{n+1}-h{y}^{n})

                             + (\sigma_0^{*}+\sigma_x^{*})h_{y}^{n-1/2}
                             + \mu^{-1}\sigma_0^{*}\sigma_x^{*}h_{y}^{I \ n+1/2}

]$

where $h_{y}^{I \ n+1/2} = h_{y}^{I \ n-1/2}+\triangle t h_{y}^{n}$ and $h_x^{I \ n+1/2} = h_{x}^{I \ n-1/2}+\triangle t h_{x}^{n}$

$\mathbf{Curl}^T\mathbf{e} = -\triangle t^{-1}\mathbf{M}^e_{s\mu}(\mathbf{h}^{n+1}-\mathbf{h}^{n}) -\mathbf{M}^e_{s\sigma^{*}}\mathbf{h}^{n} - \mathbf{M}^e_{s\mu^{-1}\sigma^{*}} \mathbf{h}^{I \ n+1/2} $

$ \mathbf{M}^e_{s\mu} = \mathbf{diag} (\mathbf{A}^{e \ T}_{c \ vec} \begin{bmatrix} \mu s_{y0} \\[0.3em] \mu 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\mu^{-1}\sigma^{*}} = \mathbf{diag} (\mathbf{A}^{e \ T}_{c \ vec} \begin{bmatrix} \mu^{-1} \sigma^{*}_0\sigma^{*}_y s_{y0} \\[0.3em] \mu^{-1} \sigma^{*}_0\sigma^{*}_x s_{x0} \end{bmatrix} ) $





In [63]:



    


# Mesmuisigs = sdiag(mesh.aveE2CCV.T*np.r_[sigs*sigy/epsilon*sy0, sigs*sigx/epsilon*sx0])
# Messigs = sdiag(mesh.aveE2CCV.T*np.r_[(sigs+sigy*mu/epsilon)*sy0, (sigs+sigx*mu/epsilon)*sx0])
# Mesmu = sdiag(mesh.aveE2CCV.T*np.r_[mu*sy0, mu*sx0])
# MesmuI = sdInv(Mesmu)
# Icc = sp.hstack((speye(mesh.nC), speye(mesh.nC)))
# curl = mesh.edgeCurl
# curlvec = sp.block_diag((curl[:,:mesh.nEx], curl[:,mesh.nEx:]))

Discretize the system:

$\mathbf{Curl}^{vec}\mathbf{h}^{n} = \triangle t^{-1}\mathbf{S}_{\epsilon}(\mathbf{e}_d^{n+1/2}-\mathbf{e}_d^{n-1/2}) +\mathbf{S}_{\sigma}\mathbf{e}_d^{n-1/2} + \mathbf{S}_{\epsilon^{-1}\sigma}\mathbf{e}_d^{I \ n} + \mathbf{j}_e $

$\triangle t^{-1}\mathbf{S}_{\epsilon}(\mathbf{e}_d^{n+1/2}-\mathbf{e}_d^{n-1/2}) = \mathbf{Curl}^{vec}\mathbf{h}^{n} -\mathbf{S}_{\sigma}\mathbf{e}_d^{n-1/2} - \mathbf{S}_{\epsilon^{-1}\sigma}\mathbf{e}_d^{I \ n} - \mathbf{j}_e $

$\mathbf{e}_d^{n+1/2} = \mathbf{e}_d^{n-1/2} + \triangle t\mathbf{S}^{-1}_{\epsilon}[\mathbf{Curl}^{vec}\mathbf{h}^{n} -\mathbf{S}_{\sigma}\mathbf{e}_d^{n-1/2} - \mathbf{S}_{\epsilon^{-1}\sigma}\mathbf{e}_d^{I \ n} - \mathbf{j}_e] $

$\mathbf{Curl}^T\mathbf{e} = -\triangle t^{-1}\mathbf{M}^e_{s\mu}(\mathbf{h}^{n+1}-\mathbf{h}^{n}) -\mathbf{M}^e_{s\sigma^{*}}\mathbf{h}^{n} - \mathbf{M}^e_{s\mu^{-1}\sigma^{*}} \mathbf{h}^{I \ n+1/2} - \mathbf{j}_m $

$-\triangle t^{-1}\mathbf{M}^e_{s\mu}(\mathbf{h}^{n+1}-\mathbf{h}^{n}) = \mathbf{Curl}^T\mathbf{e} + \mathbf{M}^e_{s\sigma^{*}}\mathbf{h}^{n} + \mathbf{M}^e_{s\mu^{-1}\sigma^{*}} \mathbf{h}^{I \ n+1/2} + \mathbf{j}_m $

$\mathbf{h}^{n+1} = \mathbf{h}^{n}-\triangle t(\mathbf{M}^e_{s\mu})^{-1}[\mathbf{Curl}^T\mathbf{e} + \mathbf{M}^e_{s\sigma^{*}}\mathbf{h}^{n} + \mathbf{M}^e_{s\mu^{-1}\sigma^{*}} \mathbf{h}^{I \ n+1/2} + \mathbf{j}_m ] $

Some useful properties:

$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})$

First compute:

$ \mathbf{e}_d^{I \ n} = \mathbf{e}_d^{I \ n-1} + \triangle t \mathbf{e}_d^{n-1/2} $

$\mathbf{e}_d^{n+1/2} = \mathbf{e}_d^{n-1/2} + \triangle t\mathbf{S}^{-1}_{\epsilon}[\mathbf{Curl}^{vec}\mathbf{h}^{n} -\mathbf{S}_{\sigma}\mathbf{e}_d^{n-1/2} - \mathbf{S}_{\epsilon^{-1}\sigma}\mathbf{e}_d^{I \ n} - \mathbf{j}_e] $

Then we compute:

$ \mathbf{h}^{I \ n+1/2} = \mathbf{h}^{I \ n-1/2} + \triangle t \mathbf{h}^{n} $

$\mathbf{h}^{n+1} = \mathbf{h}^{n}-\triangle t(\mathbf{M}^e_{s\mu})^{-1}[\mathbf{Curl}^T\mathbf{e} + \mathbf{M}^e_{s\sigma^{*}}\mathbf{h}^{n} + \mathbf{M}^e_{s\mu^{-1}\sigma^{*}} \mathbf{h}^{I \ n+1/2} + \mathbf{j}_m ] $





In [7]:



    


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 [8]:



    


Lx = mesh.hx[-npad:].sum()
Ly = mesh.hy[-npad:].sum()



    











In [41]:



    


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 [42]:



    


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)



    


















    






>> sm: 3.00e+00, lamda: 7.07e-02, sigm: 1.87e-01

















In [43]:



    


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 [44]:



    


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 [44]:



    


Seps = sp.block_diag([sdiag(epsilon*sy0), sdiag(epsilon*sx0)])
SepsI = sp.block_diag([sdiag(1./(epsilon*sy0)), sdiag(1./(epsilon*sx0))])
Sepsisig = sp.block_diag([sdiag(sig0*sigy*(1./epsilon)*sy0), sdiag(sig0*sigx*(1./epsilon)*sx0)])
Ssig = sp.block_diag([sdiag((sigy+sig0)*sy0), sdiag((sigx+sig0)*sx0)])
Mesmuisigs = sdiag(mesh.aveE2CCV.T*np.r_[sigs*sigy/epsilon*sy0, sigs*sigx/epsilon*sx0])
Messigs = sdiag(mesh.aveE2CCV.T*np.r_[(sigs+sigy*mu/epsilon)*sy0, (sigs+sigx*mu/epsilon)*sx0])
Mesmu = sdiag(mesh.aveE2CCV.T*np.r_[mu*sy0, mu*sx0])
MesmuI = sdInv(Mesmu)
Icc = sp.hstack((speye(mesh.nC), speye(mesh.nC)))
curl = mesh.edgeCurl
curlvec = sp.block_diag((curl[:,:mesh.nEx], curl[:,mesh.nEx:]))



    











In [45]:



    


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[45]:








<matplotlib.text.Text at 0x122def98>










    
























In [46]:



    


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 [47]:



    


dt = 1e-11
fmain = 3e9
time = np.arange(700)*dt
wave = ricker(fmain, time, 50*dt)



    











In [48]:



    


plt.plot(time, wave, '.-')



    


















    
Out[48]:








[<matplotlib.lines.Line2D at 0x118a05c0>]










    
























In [49]:



    


# 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')

First compute:

$ \mathbf{e}_d^{I \ n} = \mathbf{e}_d^{I \ n-1} + \triangle t \mathbf{e}_d^{n-1/2} $

$\mathbf{e}_d^{n+1/2} = \mathbf{e}_d^{n-1/2} + \triangle t\mathbf{S}^{-1}_{\epsilon}[\mathbf{Curl}^{vec}\mathbf{h}^{n} -\mathbf{S}_{\sigma}\mathbf{e}_d^{n-1/2} - \mathbf{S}_{\epsilon^{-1}\sigma}\mathbf{e}_d^{I \ n} - \mathbf{j}_e] $

Then compute:

$ \mathbf{h}^{I \ n+1/2} = \mathbf{h}^{I \ n-1/2} + \triangle t \mathbf{h}^{n} $

$\mathbf{h}^{n+1} = \mathbf{h}^{n}-\triangle t(\mathbf{M}^e_{s\mu})^{-1}[\mathbf{Curl}^T\mathbf{e} + \mathbf{M}^e_{s\sigma^{*}}\mathbf{h}^{n} + \mathbf{M}^e_{s\mu^{-1}\sigma^{*}} \mathbf{h}^{I \ n+1/2} + \mathbf{j}_m ] $





In [56]:



    


h0 = np.zeros(mesh.nE)
h1 = np.zeros(mesh.nE)
hI0 = np.zeros(mesh.nE)
hI1 = np.zeros(mesh.nE)
hx = np.zeros((mesh.nC, time.size))
hy = np.zeros((mesh.nC, time.size))

ed0 = np.zeros(2*mesh.nC)
ed1 = np.zeros(2*mesh.nC)
eId0 = np.zeros(2*mesh.nC)
eId1 = np.zeros(2*mesh.nC)
e = np.zeros((mesh.nC, time.size))

for i in range (time.size):
    
    je = np.r_[q, q]*0.5*wave[i] 
    eId0 = eId1.copy()
    eId1 = eId0 + dt*ed0
    ed1 = ed0 + SepsI*dt*(curlvec*(h1) - Ssig*ed0 - Sepsisig*eId1-je)
    ed0 = ed1.copy()
    e[:,i] = Icc*ed1

    hI0 = hI1.copy()
    hI1 = hI0 + dt*h0
    h1 = h0 - MesmuI*dt*(curl.T*(Icc*ed0)+Messigs*h0+Mesmuisigs*hI1)
#     hd1 = hd0 - SmuI*dt*(curlvec*e0+Ssig*hd0+Smuisig*hId1)    
    hx[:,i] = Px*h1
    hy[:,i] = Py*h1
    h0 = h1.copy()



    











In [65]:



    


icount = 500
extent = [mesh.vectorCCx.min(), mesh.vectorCCx.max(), mesh.vectorCCy.min(), mesh.vectorCCy.max()]
plt.imshow(np.flipud(e[:,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[65]:








<matplotlib.image.AxesImage at 0x1a5c2c88>










    
























In [61]:



    


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('$e_z$', fontsize = 20)
ax[1].set_title('$h_x$', fontsize = 20)
ax[2].set_title('$h_y$', fontsize = 20)
nskip = 40
def animate(i_id):
    icount = i_id*nskip
    frame1 = ax[0].imshow(np.flipud(e[:,icount].reshape((mesh.nCx, mesh.nCy), order = 'F').T), cmap = 'RdBu', extent=extent)    
    frame2 = ax[1].imshow(np.flipud(hx[:,icount].reshape((mesh.nCx, mesh.nCy), order = 'F').T), cmap = 'RdBu', extent=extent)
    frame3 = ax[2].imshow(np.flipud(hy[:,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[61]:











    
    

    
    

    
    
    
    
    
    
    
    
    
  

     Once 
     Loop 
     Reflect

Content source: simpeg/simpegGPR

Similar notebooks:

notebook.community | gallery | about