notebook.community

Edit and run



In [1]:

    
from SimPEG import *
%pylab inline









    



Populating the interactive namespace from numpy and matplotlib



In [2]:

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

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

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

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



In [6]:

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



In [7]:

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



In [48]:

    
mopt = np.load('./inv2D_FD_realistic_line_1/model_14.npy')



In [49]:

    
mesh2D = Mesh.TensorMesh([hx,hz], 'CC')
active = mesh2D.gridCC[:,1] < 0.
actMap = Maps.ActiveCells(mesh2D, active, np.log(1e-8), nC=mesh2D.nC)
map2to3 = Maps.Map2Dto3D(mesh, normal = 'Y')
mapping =  Maps.ExpMap(mesh) * map2to3 * actMap



In [50]:

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

    
fig, ax = plt.subplots(1,2, figsize=(18,4))
indz = 20
print mesh.vectorCCz[indz]
mesh.plotSlice(np.log10(sigma), ind = indz, ax = ax[0], clim=(-3, -1))
mesh.plotSlice(np.log10(mapping*mopt), ind = indz, ax = ax[1], clim=(-3, -1))
# mesh.plotSlice(np.log10(mapping*mopt), ind = indz, ax = ax[1])
for i in range(2):
    ax[i].plot(xyz1[:,0], xyz1[:,1], 'r.')
    ax[i].plot(xyz2[:,0], xyz2[:,1], 'b.')
#     ax[i].plot(xcirc1, ycirc1, 'r-')
#     ax[i].plot(xcirc2, ycirc2, 'b-')
    ax[i].set_xlim(-320, 320)
    ax[i].set_ylim(-150, 150)



In [52]:

    
temp = mapping*mopt
np.save('./inv2D_FD_realistic_line/invTEM2D', temp)



In [53]:

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

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



In [27]:

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



In [28]:

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

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

    
# 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 [31]:

    
dobs = np.load('bzobs_FD_realistic_line.npy')



In [42]:

    
dest = np.load('inv2D_FD_realistic_line/dpred_14.npy')



In [43]:

    
frequency = np.r_[1, 10., 100.]



In [44]:

    
dest.shape









    Out[44]:





(360L,)



In [45]:

    
Dpred = abs(dobs.reshape((30, 2, frequency.size, 2), order='F'))
Dest = abs(dest.reshape((30, 2, frequency.size, 2), order='F'))



In [46]:

    
ind = xyz1[:,1] == 0.
absFD = lambda x, y: np.sqrt(x**2+y**2)
mradFD = lambda x, y: np.angle(x+1j*y)*1e3
for itime in range(3):
    abs1 = absFD(Utils.mkvc(Dpred[:,0,itime,0]), Utils.mkvc(Dpred[:,1,itime,1]))
    abs2 = absFD(Utils.mkvc(Dpred[:,0,itime,0]), Utils.mkvc(Dpred[:,1,itime,1]))
    plt.semilogy(np.r_[xyz1[ind,0], xyz2[ind,0]], np.r_[abs1, abs2] )
    abs1 = absFD(Utils.mkvc(Dest[:,0,itime,0]), Utils.mkvc(Dest[:,1,itime,1]))
    abs2 = absFD(Utils.mkvc(Dest[:,0,itime,0]), Utils.mkvc(Dest[:,1,itime,1]))    
    plt.semilogy(np.r_[xyz1[ind,0], xyz2[ind,0]], np.r_[abs1, abs2])



In [47]:

    
for itime in range(3):
    phase1 = mradFD(Utils.mkvc(Dpred[:,0,itime,0]), Utils.mkvc(Dpred[:,1,itime,1]))
    phase2 = mradFD(Utils.mkvc(Dpred[:,0,itime,0]), Utils.mkvc(Dpred[:,1,itime,1]))
    plt.plot(np.r_[xyz1[ind,0], xyz2[ind,0]], np.r_[phase1, phase2] )
    phase1 = mradFD(Utils.mkvc(Dest[:,0,itime,0]), Utils.mkvc(Dest[:,1,itime,1]))
    phase2 = mradFD(Utils.mkvc(Dest[:,0,itime,0]), Utils.mkvc(Dest[:,1,itime,1]))    
    plt.plot(np.r_[xyz1[ind,0], xyz2[ind,0]], np.r_[phase1, phase2])



In [33]:

    
# for itime in range(3):
#     plt.semilogy(np.r_[xyz1[ind,0], xyz2[ind,0]], np.r_[Utils.mkvc(Dpred[:,0,itime,0]), Utils.mkvc(Dpred[:,0,itime,1])])
#     plt.semilogy(np.r_[xyz1[ind,0], xyz2[ind,0]], np.r_[Utils.mkvc(Dest[:,0,itime,0]), Utils.mkvc(Dest[:,0,itime,1])])



In [175]: