

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: $TE^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 ($e_x$, $e_y$, $h_z$):

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

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



In [4]:

    
sighalf = 0.
mu = mu_0*np.ones(mesh.nC)
epsilon = epsilon_0*np.ones(mesh.nC)*25
sigma = sighalf*np.ones(mesh.nC)
Mmui = sdiag(1./mu)
Meeps = sdiag(mesh.aveE2CC.T*epsilon)
MeepsI = sdInv(Meeps)
Mesig = sdiag(mesh.aveE2CC.T*sigma)
curl = mesh.edgeCurl

First compute:

$\mathbf{e}^{n+1/2} = \mathbf{e}^{n-1/2} - (\mathbf{M}^{e}_{\epsilon})^{-1}\triangle t(-\mathbf{Curl}^{T}\mathbf{h}^{n}+\mathbf{M}^{e}_{\sigma}\mathbf{e}^{n}+\mathbf{j}^{n-1/2}_s)$

Then compute:

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



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 0xae7a400>]



In [8]:

    
txind = Utils.closestPoints(mesh, [0., 0.], gridLoc='Ex')
q = np.zeros(mesh.nE)
q[txind] = 1.
q = Utils.sdiag(1/mesh.edge)*q



In [9]:

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



In [10]:

    
h0 = np.zeros(mesh.nC)
h1 = np.zeros(mesh.nC)
e0 = np.zeros(mesh.nE)
e1 = 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] 
    e1 = e0 - MeepsI*dt*(-curl.T*h1+Mesig*e0+js)
    e0 = e1.copy()
    ex[:,i] = Px*e0
    ey[:,i] = Py*e0
    h1 = h0 - Mmui*dt*curl*e0
    h0 = h1.copy()
    h[:,i] = h1



In [11]:

    
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:])

nskip = 30
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, vmin = ex[:,icount].min()*0.05, vmax = ex[:,icount].max()*0.05)
    frame3 = ax[2].imshow(np.flipud(ey[:,icount].reshape((mesh.nCx, mesh.nCy), order = 'F').T), cmap = 'RdBu', extent=extent, vmin = ey[:,icount].min()*0.05, vmax = ey[:,icount].max()*0.05)    
    return frame1, frame2, frame3
animation.FuncAnimation(fig, animate, frames=14, interval=40, blit=True)