notebook.community

Edit and run 





In [20]:



    


from SimPEG import *
%pylab inline



    


















    






Populating the interactive namespace from numpy and matplotlib

















In [70]:



    


import matplotlib
from IPython.display import set_matplotlib_formats
set_matplotlib_formats('png')
matplotlib.rcParams['savefig.dpi'] = 100



    











In [71]:



    


cs, ncx, ncy, ncz, npad = 20., 30, 20, 30, 12
hx = [(cs,npad,-1.4), (cs,ncx), (cs,npad,1.4)]
hy = [(cs,npad,-1.4), (cs,ncy), (cs,npad,1.4)]
hz = [(cs,npad,-1.4), (cs,ncz), (cs,npad,1.4)]
mesh = Mesh.TensorMesh([hx,hy,hz], 'CCC')
print ("Padding distance x: %10.5f m") % (np.sum(mesh.hx[:npad]))
print ("Padding distance z: %10.5f m") % (np.sum(mesh.hz[:npad]))
print ("Min dx: %10.5f m") % (mesh.hx.min())
print ("Min dz: %10.5f m") % (mesh.hz.min())



    


















    






Padding distance x: 3898.57387 m
Padding distance z: 3898.57387 m
Min dx:   20.00000 m
Min dz:   20.00000 m

















In [72]:



    


print mesh



    


















    






  ---- 3-D TensorMesh ----  
   x0: -4198.57
   y0: -4098.57
   z0: -4198.57
  nCx: 54
  nCy: 44
  nCz: 54
   hx: 1133.88, 809.91, 578.51, 413.22, 295.16, 210.83, 150.59, 107.56, 76.83, 54.88, 39.20, 28.00, 30*20.00, 28.00, 39.20, 54.88, 76.83, 107.56, 150.59, 210.83, 295.16, 413.22, 578.51, 809.91, 1133.88
   hy: 1133.88, 809.91, 578.51, 413.22, 295.16, 210.83, 150.59, 107.56, 76.83, 54.88, 39.20, 28.00, 20*20.00, 28.00, 39.20, 54.88, 76.83, 107.56, 150.59, 210.83, 295.16, 413.22, 578.51, 809.91, 1133.88
   hz: 1133.88, 809.91, 578.51, 413.22, 295.16, 210.83, 150.59, 107.56, 76.83, 54.88, 39.20, 28.00, 30*20.00, 28.00, 39.20, 54.88, 76.83, 107.56, 150.59, 210.83, 295.16, 413.22, 578.51, 809.91, 1133.88

















In [73]:



    


def circfun(xc, yc, r, npoint):
    theta = np.linspace(np.pi, -np.pi, npoint)
    x = r*np.cos(theta)
    y = r*np.sin(theta)
    return x+xc, y+yc



    











In [74]:



    


xcirc1, ycirc1 = circfun(-150., 0., 250., 60)
xcirc2, ycirc2 = circfun(150., 0., 250., 60)



    











In [75]:



    


Utils.meshutils.writeUBCTensorMesh(mesh, 'mesh.msh')



    











In [76]:



    


sigma = Utils.meshutils.readUBCTensorModel('sigma_realistic.con', mesh)



    











In [77]:



    


mopt = np.load('./inv3D/model_10.npy')
active = mesh.gridCC[:,2] < 0.
actMap = Maps.ActiveCells(mesh, active, np.log(1e-8), nC=mesh.nC)
mapping = Maps.ExpMap(mesh) * actMap



    











In [78]:



    


x1 = np.arange(30)*10 - 300.
y1 = np.arange(30)*10 - 150.
xyz1 = Utils.ndgrid(x1, y1, np.r_[0.])
x2 = np.arange(30)*10 + 10.
y2 = np.arange(30)*10 - 150.
xyz2 = Utils.ndgrid(x2, y2, np.r_[0.])



    











In [97]:



    


fig, ax = plt.subplots(1,1, figsize=(9,6))
indz1 = 20
indz2 = indz1
indy=20
print mesh.vectorCCz[indz1]
mesh.plotSlice(np.log10(sigma)*np.nan, ind = indz1, ax = ax, grid=True, pcolorOpts={'cmap':'binary'})
ind = np.argwhere(xyz1[:,1]==0.)
ax.plot(xyz1[:,0], xyz1[:,1], 'r.')
ax.plot(xyz2[:,0], xyz2[:,1], 'b.')
ax.plot(xcirc1, ycirc1, 'r-')
ax.plot(xcirc2, ycirc2, 'b-')
ax.plot(xyz1[ind,0], xyz1[ind,1], 'k.', ms=10)
ax.plot(xyz2[ind,0], xyz2[ind,1], 'k.', ms=10)
ax.set_xlim(-400, 400)
ax.set_ylim(-300, 300)
ax.set_xlabel('Easting (m)', fontsize = 16)
ax.set_ylabel('Northing (m)', fontsize = 16)
ax.set_title('Plan view', fontsize = 16)
ax.text(-300., 250, 'Tx Loop#1', color='r', fontsize = 16)
ax.text(0, 250, 'Tx Loop#2', color='b', fontsize = 16)
ax.text(-300., 160, 'Rx group#1', color='r', fontsize = 16)
ax.text(0, 160, 'Rx group#2', color='b', fontsize = 16)
fig.savefig('./figures/geometry.png', dpi = 200)



    


















    






-130.0










    
























In [82]:



    


# indy=20
# fig, ax = plt.subplots(1,2, figsize=(16,4))
# print mesh.vectorCCy[indy]
# dat0 = mesh.plotSlice(np.log10(sigma), normal='Y', ind = indy, ax = ax[0], clim=(-3, -1))
# dat1 = mesh.plotSlice(np.log10(mapping*mopt), normal='Y', ind = indy, ax = ax[1], clim=(-3, -1))
# for i in range(2):
#     plt.colorbar(dat0[0], ax = ax[i])
#     ax[i].set_xlim(-300, 300)
#     ax[i].set_ylim(-700, 0.)



    











In [20]:



    


# indy=23
# fig, ax = plt.subplots(1,2, figsize=(16,4))
# print mesh.vectorCCy[indy]
# dat0 = mesh.plotSlice(np.log10(sigma), normal='Y', ind = indy, ax = ax[0], clim=(-3, -0.5))
# dat1 = mesh.plotSlice(np.log10(mapping*mopt), normal='Y', ind = indy, ax = ax[1], clim=(-3, -0.5))
# for i in range(2):
#     plt.colorbar(dat0[0], ax = ax[i])
#     ax[i].set_xlim(-300, 300)
#     ax[i].set_ylim(-300, 0.)



    











In [21]:



    


import simpegEM as EM
mapping = Maps.IdentityMap(mesh)



    











In [22]:



    


ntx = 2
nrx1 = xyz1.shape[0]
time = np.logspace(-4, -2, 31)



    











In [23]:



    


print time



    


















    






[ 0.0001      0.00011659  0.00013594  0.00015849  0.00018478  0.00021544
  0.00025119  0.00029286  0.00034145  0.00039811  0.00046416  0.00054117
  0.00063096  0.00073564  0.0008577   0.001       0.00116591  0.00135936
  0.00158489  0.00184785  0.00215443  0.00251189  0.00292864  0.00341455
  0.00398107  0.00464159  0.0054117   0.00630957  0.00735642  0.00857696
  0.01      ]

















In [24]:



    


rx1 = EM.TDEM.RxTDEM(xyz1, time, 'bz')
tx1 = EM.TDEM.TxTDEM(np.array([0., -150., 0.]), 'CircularLoop_MVP', [rx1])
tx1.radius = 250.
rx2 = EM.TDEM.RxTDEM(xyz2, time, 'bz')
tx2 = EM.TDEM.TxTDEM(np.array([0.,  150., 0.]), 'CircularLoop_MVP', [rx2])
tx2.radius = 250.



    











In [25]:



    


# survey = EM.TDEM.SurveyTDEM([tx1, tx2])
survey = EM.TDEM.SurveyTDEM([tx1])
prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping, verbose=True)
# prb.solver = MumpsSolver
# prb.solverOpts = {"symmetric":True}
# prb.timeSteps = [(1e-4/15, 10), (1e-3/15, 10), (1e-2/15, 5)]
prb.timeSteps = [(1e-4/15, 10)]
if prb.ispaired:
    prb.unpair()
if survey.ispaired:
    survey.unpair()
prb.pair(survey)



    











In [26]:



    


dpred = np.load('../codes/bz.npy')



    











In [27]:



    


noise = abs(dpred)*np.random.randn(dpred.size)*0.05 
dobs = dpred+noise



    











In [28]:



    


dest = np.load('inv3D/dpred_15.npy')



    











In [29]:



    


Dpred = dobs.reshape((900, 31, 2), order='F')
Dpred1 = Dpred[:,:,0]
Dpred2 = Dpred[:,:,1]
Dest = dest.reshape((900, 31, 2), order='F')
Dest1 = Dest[:,:,0]
Dest2 = Dest[:,:,1]



    











In [30]:



    


itime = 16
fig, ax = plt.subplots(1,2, figsize = (12, 5))
for i in range(2):
    ax[i].set_xlabel('x (m)', fontsize = 16)
    ax[i].set_ylabel('y (m)', fontsize = 16)
vmin = Utils.mkvc(Dpred[:,itime,0]).min()
vmax = Utils.mkvc(Dpred[:,itime,0]).max()
dat1 = ax[0].contourf(xyz1[:,0].reshape((30, 30), order='F'), xyz1[:,1].reshape((30, 30), order='F'),  Dpred[:,itime,0].reshape((30, 30), order='F'), 30, vmin=vmin, vmax=vmax)
dat1 = ax[1].contourf(xyz2[:,0].reshape((30, 30), order='F'), xyz2[:,1].reshape((30, 30), order='F'),  Dest[:,itime,0].reshape((30, 30), order='F'), 30, vmin=vmin, vmax=vmax)



    


















    
























In [31]:



    


itime = 9
fig, ax = plt.subplots(1,2, figsize = (12, 5))
for i in range(2):
    ax[i].set_xlabel('x (m)', fontsize = 16)
    ax[i].set_ylabel('y (m)', fontsize = 16)
vmin = Utils.mkvc(Dpred[:,itime,1]).min()
vmax = Utils.mkvc(Dpred[:,itime,1]).max()    
dat1 = ax[0].contourf(xyz1[:,0].reshape((30, 30), order='F'), xyz1[:,1].reshape((30, 30), order='F'),  Dpred[:,itime,1].reshape((30, 30), order='F'), 30)
dat1 = ax[1].contourf(xyz2[:,0].reshape((30, 30), order='F'), xyz2[:,1].reshape((30, 30), order='F'),  Dest[:,itime,1].reshape((30, 30), order='F'), 30)



    


















    
























In [32]:



    


rxind = 20
plt.loglog(time, Dpred[rxind,:,0])
plt.loglog(time, Dest[rxind,:,0])



    


















    
Out[32]:








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










    
























In [33]:



    


itime = 0
fig, ax = plt.subplots(1,2, figsize = (10, 4))
ax[0].hist(Dpred1[:,itime])
ax[1].hist(Dpred2[:,itime])



    


















    
Out[33]:








(array([   5.,   13.,   65.,   98.,  157.,  165.,  194.,  139.,   52.,   12.]),
 array([  5.41459373e-10,   6.02450111e-10,   6.63440848e-10,
          7.24431586e-10,   7.85422324e-10,   8.46413061e-10,
          9.07403799e-10,   9.68394537e-10,   1.02938527e-09,
          1.09037601e-09,   1.15136675e-09]),
 <a list of 10 Patch objects>)










    
























In [34]:



    


from JSAnimation import IPython_display
from matplotlib import animation
fig, ax = plt.subplots(1,2, figsize = (12, 5))
for i in range(2):
    ax[i].set_xlabel('x (m)', fontsize = 16)
    ax[i].set_ylabel('y (m)', fontsize = 16)

def animate(itime):
    frame1 = ax[0].contourf(xyz1[:,0].reshape((30, 30), order='F'), xyz1[:,1].reshape((30, 30), order='F'),  Dpred[:,itime,0].reshape((30, 30), order='F'), 30)
    frame2 = ax[1].contourf(xyz2[:,0].reshape((30, 30), order='F'), xyz2[:,1].reshape((30, 30), order='F'),  Dpred[:,itime,1].reshape((30, 30), order='F'), 30)
#     cb1 = plt.colorbar(frame1, ax = ax[0])
#     cb2 = plt.colorbar(frame2, ax = ax[1])    
    return frame1, frame2
animation.FuncAnimation(fig, animate, frames=31, interval=40, blit=True)



    


















    
Out[34]:











    
    

    
    

    
    
    
    
    
    
    
    
    
  

     Once 
     Loop 
     Reflect 
  























In [ ]:



    


# np.save('../bzobs', dobs)
# dmis = DataMisfit.l2_DataMisfit(survey)

Content source: sgkang/AGU2014MovingDimensionsinEM

Similar notebooks:

notebook.community | gallery | about