

In [1]:

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









    



Populating the interactive namespace from numpy and matplotlib



In [2]:

    
from SimPEG.Utils import mkvc, sdiag, sdInv



In [3]:

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

2D Geo Penetrating Radar (GPR) modeling: $TM^z$ mode

Maxwell's equations in time domain:

$ \nabla \times \vec{e} = -\mu \frac{\partial\vec{h}}{\partial t}$

$ \nabla \times \vec{h} -\sigma \vec{e} - \epsilon \frac{\partial\vec{e}}{\partial t} = \vec{j}_s $

Assuming TE mode ($h_x$, $h_y$, $e_z$):

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

$\mathbf{Curl}\mathbf{h}^{n+1/2} - \mathbf{diag}(\sigma)\mathbf{e}^{n} - \mathbf{diag}(\epsilon)\triangle t^{-1}(\mathbf{e}^{n+1}-\mathbf{e}^{n}) = \mathbf{j}^n_s$



In [4]:

    
sighalf = 1e-4
mu = mu_0*np.ones(mesh.nC)
epsilon = epsilon_0*np.ones(mesh.nC)*25
sigma = sighalf*np.ones(mesh.nC)
Memu = sdiag(mesh.aveE2CC.T*mu)
MemuI = sdInv(Memu)
MepsI = sdiag(1./epsilon)
Msig = sdiag(sigma)
curl = mesh.edgeCurl
V = sdiag(mesh.vol)

First compute:

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

Then compute:

$\mathbf{e}^{n+1} = \mathbf{e}^{n} - (\mathbf{diag}(\epsilon))^{-1}\triangle t(-\mathbf{Curl}\mathbf{h}^{n+1/2}+\mathbf{j}^n_s+\mathbf{diag}(\sigma)\mathbf{e}^{n})$



In [5]:

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

    
dt = 2e-10
time = np.arange(500)*dt
wave = ricker(200*1e6, time, 35*dt)



In [7]:

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









    Out[7]:





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



In [8]:

    
Px = mesh.getInterpolationMat(mesh.gridCC, 'Ex')
Py = mesh.getInterpolationMat(mesh.gridCC, 'Ey')



In [9]:

    
h0 = np.zeros(mesh.nE)
h1 = np.zeros(mesh.nE)
e0 = np.zeros(mesh.nC)
e1 = np.zeros(mesh.nC)
txind = Utils.closestPoints(mesh, [0., 0.], gridLoc='CC')
q = np.zeros(mesh.nC)
q[txind] = 1.
q = Utils.sdiag(1/mesh.vol)*q
e = np.zeros((mesh.nC, time.size))
hx = np.zeros((mesh.nC, time.size))
hy = np.zeros((mesh.nC, time.size))

for i in range (time.size):
    
    js = q*wave[i] 
    e1 = e0 - MepsI*dt*(-curl*h1+js+Msig*e0)
    e[:,i] = e1    
    e0 = e1.copy()
    h1 = h0 - MemuI*dt*curl.T*e0
    h0 = h1.copy()
    hx[:,i] = Px*h1
    hy[:,i] = Py*h1



In [10]:

    
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('Easting (m)', fontsize = 16)
    ax[i].set_ylabel('Depth (m)', fontsize = 16)
    ax[i].set_xlim(extent[:2])
    ax[i].set_ylim(extent[2:])

nskip = 30
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)    
    return frame1, frame2, frame3
animation.FuncAnimation(fig, animate, frames=14, interval=40, blit=True)