In [2]:
from SimPEG import *
import simpegEM as EM
from simpegem1d import Utils1D
%pylab inline
In [3]:
mesh3D = Utils.meshutils.readUBCTensorMesh('mesh.msh')
sigma3D = Utils.meshutils.readUBCTensorModel('sigma_simple.con', mesh3D)
In [4]:
x2 = np.arange(30)*10 + 10.
y2 = np.arange(30)*10 - 150.
xyz2 = Utils.ndgrid(x2, y2, np.r_[0.])
xc2 = 150
yc2 = 0.
r2 = np.sqrt((xyz2[:,0]-xc2)**2+(xyz2[:,1]-yc2)**2)
In [5]:
dobs = np.load('bzobs_simple.npy')
Dobs = dobs.reshape((900, 31, 2), order='F')
Dobs1 = Dobs[:,:,0]
Dobs2 = Dobs[:,:,1]
In [6]:
fig, ax = plt.subplots(1,2, figsize = (12, 5))
ax[0].contourf(Dobs1)
ax[1].contourf(Dobs2)
Out[6]:
In [7]:
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 [8]:
active = mesh.vectorCCz<0.
layer1 = (mesh.vectorCCz<0.) & (mesh.vectorCCz>=-60.)
layer2 = (mesh.vectorCCz<-60) & (mesh.vectorCCz>=-100.)
layer3 = (mesh.vectorCCz<-100) & (mesh.vectorCCz>=-200.)
actMap = Maps.ActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
mapping = Maps.ExpMap(mesh) * Maps.Vertical1DMap(mesh) * actMap
In [9]:
sig_half = 1e-3
sig_air = 1e-8
sig_layer1 = 1./300
sig_layer2 = 1./100
sig_layer3 = 1./10
sigma = np.ones(mesh.nCz)*sig_air
sigma[active] = sig_half
sigma[layer1] = sig_layer1
sigma[layer2] = sig_layer2
sigma[layer3] = sig_layer3
mtrue = np.log(sigma[active])
In [10]:
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 [11]:
xcirc1, ycirc1 = circfun(-150., 0., 250., 60)
xcirc2, ycirc2 = circfun(150., 0., 250., 60)
In [12]:
ind = np.argwhere(xyz2[:,1]==0)
In [14]:
fig, ax = plt.subplots(1,1, figsize=(7,3))
indz = 20
print mesh.vectorCCz[indz]
mesh3D.plotSlice(np.log10(sigma3D), ind = indz, ax = ax, clim=(-3, -0.5))
ax.plot(xyz2[:,0], xyz2[:,1], 'b.')
ax.plot(xyz2[ind,0], xyz2[ind,1], 'k.', ms=10)
# ax.plot(xyz2[:,0], xyz2[:,1], 'b.')
# ax.plot(xcirc1, ycirc1, 'r-')
ax.plot(xcirc2, ycirc2, 'b-')
ax.set_xlim(-500, 500)
ax.set_ylim(-300, 300)
Out[14]:
In [140]:
sig_test = (sigma[active])
In [141]:
Sig_test = (sig_test.reshape([1,-1])).repeat(10, axis=0)
x = xyz2[ind,0]
z = mesh.vectorCCz[active]
Z, X = np.meshgrid(z, x)
In [142]:
z = mesh.vectorCCz[active]
Time, Xtime = np.meshgrid(time, x)
In [143]:
# plt.pcolor(X, Z, np.log10(Sig_test))
# plt.ylim()
In [144]:
fig, ax = plt.subplots(1,1, figsize = (3, 6))
Utils1D.plotLayer(sig_test, mesh.vectorCCz[active], showlayers=True, ax = ax)
ax.set_ylim(-300., 0.)
Out[144]:
In [145]:
prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping, verbose=False)
In [146]:
prb.Solver = SolverLU
prb.timeSteps = [(1e-4/10, 15), (1e-3/10, 15), (1e-2/10, 15), (1e-1/10, 15)]
In [147]:
# Mopt = []
# Dest=[]
# for i in range(Dobs2[ind,:].shape[0]):
# # for i in range(1):
# rxoffset=r2[ind][i]
# 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)
# std = 0.2
# survey.dobs = Utils.mkvc(Dobs2[ind,:][i,:])
# survey.std = survey.dobs*0 + std
# dmisfit = DataMisfit.l2_DataMisfit(survey)
# dmisfit.Wd = 1/(abs(survey.dobs)*std)
# regMesh = Mesh.TensorMesh([mesh.hz[mapping.maps[-1].indActive]])
# reg = Regularization.Tikhonov(regMesh)
# opt = Optimization.InexactGaussNewton(maxIter = 5)
# invProb = InvProblem.BaseInvProblem(dmisfit, reg, opt)
# # Create an inversion object
# beta = Directives.BetaSchedule(coolingFactor=5, coolingRate=2)
# betaest = Directives.BetaEstimate_ByEig(beta0_ratio=1e0)
# inv = Inversion.BaseInversion(invProb, directiveList=[beta,betaest])
# m0 = np.log(np.ones(mtrue.size)*2e-3)
# reg.alpha_s = 1e-2
# reg.alpha_x = 1.
# prb.counter = opt.counter = Utils.Counter()
# opt.LSshorten = 0.5
# opt.remember('xc')
# mopt = inv.run(m0)
# Mopt.append(mopt)
# Dest.append(invProb.dpred)
In [148]:
# np.save('Mopt2', Mopt)
# np.save('Dest2', Dest)
In [149]:
Mopt = np.load('Mopt2.npy')
Dest = np.load('Dest2.npy')
In [150]:
# 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 [151]:
SigMat = np.exp(np.vstack(Mopt))
DpreMat = np.vstack(Dest)
In [156]:
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 = -0.5)
mesh3D.plotSlice(np.log10(sigma3D), ind = 20, normal='Y', ax = ax[1], clim=(-3, -0.5))
for i in range(2):
ax[i].set_ylim(-700., 0.)
ax[i].set_xlim(10., 300.)
In [153]:
DobsMat = Utils.mkvc(Dobs2[ind, :]).reshape((30, 31), order='F')
In [154]:
fig, ax = plt.subplots(1,2, figsize = (12, 5))
vmin = np.log10(Utils.mkvc(DpreMat).min())
vmax = np.log10(Utils.mkvc(DpreMat).max())
ax[0].pcolor(Xtime, np.log10(Time), np.log10(DpreMat), vmin = vmin, vmax = vmax)
ax[1].pcolor(Xtime, np.log10(Time), np.log10(DobsMat), vmin = vmin, vmax = vmax)
Out[154]: