In [3]:

    

from SimPEG import *
import simpegEM as EM
from simpegem1d import Utils1D
%pylab inline


    

    




Populating the interactive namespace from numpy and matplotlib

In [4]:

    

meshType = 'CYL'
cs, ncx, ncz, npad = 20., 25, 30, 12
hx = [(cs,ncx), (cs,npad,1.3)]
hz = [(cs,npad,-1.4), (cs,ncz), (cs,npad,1.4)]
mesh = Mesh.CylMesh([hx,1,hz], '00C')


    

In [5]:

    

active = mesh.vectorCCz<0.
layer1 = (mesh.vectorCCz<0.) & (mesh.vectorCCz>=-60.)
layer3 = (mesh.vectorCCz<-60) & (mesh.vectorCCz>=-160.)
actMap = Maps.ActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
mapping = Maps.ExpMap(mesh) * Maps.Vertical1DMap(mesh) * actMap


    

In [6]:

    

sig_half = 1e-3
sig_air = 1e-8
sig_layer1 = 1./300
sig_layer3 = 1./10
sigma = np.ones(mesh.nCz)*sig_air
sigma[active] = sig_half
sigma[layer1] = sig_layer1
sigma[layer3] = sig_layer3
mtrue = np.log(sigma[active])


    

In [7]:

    

fig, ax = plt.subplots(1,1, figsize = (3, 6))
Utils1D.plotLayer(np.exp(mtrue), mesh.vectorCCz[active], showlayers=True, ax = ax)
ax.set_ylim(-500., 0.)


    

    
Out[7]:





(-500.0, 0.0)






    

In [8]:

    

prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping, verbose=False)
prb.Solver = SolverLU
prb.timeSteps = [(1e-4/10, 15), (1e-3/10, 15), (1e-2/10, 15), (1e-1/10, 15)]


    

In [9]:

    

time = np.logspace(-4, -2, 31)
rx = EM.TDEM.RxTDEM(np.array([[rxoffset, 0., 0.]]), time, 'bz')
tx = EM.TDEM.TxTDEM(np.array([0., 0., 0.]), 'CircularLoop_MVP', [rx])
tx.radius = 250.
survey = EM.TDEM.SurveyTDEM([tx])
if prb.ispaired:
    prb.unpair()
if survey.ispaired:
    survey.unpair()
prb.pair(survey)


    

    




---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-9-e92f951b6f87> in <module>()
      1 time = np.logspace(-4, -2, 31)
----> 2 rx = EM.TDEM.RxTDEM(np.array([[rxoffset, 0., 0.]]), time, 'bz')
      3 tx = EM.TDEM.TxTDEM(np.array([0., 0., 0.]), 'CircularLoop_MVP', [rx])
      4 tx.radius = 250.
      5 survey = EM.TDEM.SurveyTDEM([tx])

NameError: name 'rxoffset' is not defined

In [10]:

    

np.save('Mopt1_realistic', Mopt)
np.save('Dest1_realistic', Dest)


    

    




---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-10-291581f4fcc9> in <module>()
----> 1 np.save('Mopt1_realistic', Mopt)
      2 np.save('Dest1_realistic', Dest)

NameError: name 'Mopt' is not defined

In [ ]:

    

Mopt = np.load('Mopt1_realistic.npy')
Dest = np.load('Dest1_realistic.npy')


    

In [ ]:

    

# fig, ax = plt.subplots(1,1, figsize = (3, 6))
# Utils1D.plotLayer(np.exp(mopt), mesh.vectorCCz[active], showlayers=True, ax = ax)
# ax.set_ylim(-300., 0.)


    

In [ ]:

    

SigMat = np.exp(np.vstack(Mopt))
DpreMat = np.vstack(Dest)


    

In [ ]:

    

plt.pcolor(X, Z, np.log10(SigMat))
plt.ylim(-700, 0.)


    

In [ ]:

    

DobsMat = Utils.mkvc(Dobs1[ind, :]).reshape((30, 31), order='F')


    

In [ ]:

    

fig, ax = plt.subplots(1,2, figsize = (12, 5))
vmin = np.log10(Utils.mkvc(SigMat).min())
vmax = np.log10(Utils.mkvc(SigMat).max())
ax[0].contourf(X, Z, np.log10(SigMat), 31, vmin = -3, vmax = -1)
mesh3D.plotSlice(np.log10(sigma3D), ind = 20, normal='Y', ax = ax[1], clim=(-3, -1))
for i in range(2):
    ax[i].set_ylim(-700., 0.)
    ax[i].set_xlim(-300., -10.)


    

In [ ]:

    

fig, ax = plt.subplots(1,2, figsize = (12, 5))
vmin = np.log10(Utils.mkvc(DpreMat).min())
vmax = np.log10(Utils.mkvc(DpreMat).max())
ax[1].contourf(Xtime, np.log10(Time), np.log10(DpreMat), 31, vmin = vmin, vmax = vmax)
ax[0].contourf(Xtime, np.log10(Time), np.log10(DobsMat), 31, vmin = vmin, vmax = vmax)