In [1]:
from SimPEG import *
import simpegEM as EM
from simpegem1d import Utils1D
%pylab inline
In [2]:
mesh3D = Utils.meshutils.readUBCTensorMesh('mesh.msh')
sigma3D = Utils.meshutils.readUBCTensorModel('sigma_simple.con', mesh3D)
In [3]:
x1 = np.arange(30)*10 - 300.
y1 = np.arange(30)*10 - 150.
xyz1 = Utils.ndgrid(x1, y1, np.r_[0.])
xc1 = -150
yc1 = 0.
r1 = np.sqrt((xyz1[:,0]-xc1)**2+(xyz1[:,1]-yc1)**2)
In [4]:
dobs = np.load('bzobs_simple.npy')
Dobs = dobs.reshape((900, 31, 2), order='F')
Dobs1 = Dobs[:,:,0]
Dobs2 = Dobs[:,:,1]
In [5]:
fig, ax = plt.subplots(1,2, figsize = (12, 5))
ax[0].contourf(Dobs1)
ax[1].contourf(Dobs2)
Out[5]:
In [6]:
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 [7]:
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 [8]:
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 [9]:
xc = -100
yc = 100.
dind = np.argmin(abs( xyz1[:,0]-xc)+abs( xyz1[:,1]-yc))
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(xyz1[:,1]==0.)
In [13]:
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(xyz1[:,0], xyz1[:,1], 'r.')
ax.plot(xyz1[ind,0], xyz1[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[13]:
In [14]:
sig_test = (sigma[active])
In [15]:
Sig_test = (sig_test.reshape([1,-1])).repeat(10, axis=0)
x = xyz1[ind,0]
z = mesh.vectorCCz[active]
Z, X = np.meshgrid(z, x)
In [16]:
time = np.logspace(-4, -2, 31)
z = mesh.vectorCCz[active]
Time, Xtime = np.meshgrid(time, x)
In [17]:
# plt.pcolor(X, Z, np.log10(Sig_test))
# plt.ylim()
In [18]:
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[18]:
In [19]:
prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping, verbose=False)
In [20]:
prb.Solver = SolverLU
prb.timeSteps = [(1e-4/10, 15), (1e-3/10, 15), (1e-2/10, 15), (1e-1/10, 15)]
In [29]:
Mopt = []
Dest=[]
# for i in range(Dobs1[ind,:].shape[0]):
for i in range(1):
rxoffset=r1[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(Dobs1[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)
reg.mref = np.log(sig_test)
opt = Optimization.InexactGaussNewton(maxIter = 3)
invProb = InvProblem.BaseInvProblem(dmisfit, reg, opt)
# Create an inversion object
beta = Directives.BetaSchedule(coolingFactor=5, coolingRate=2)
invProb.beta = 1e0
betaest = Directives.BetaEstimate_ByEig(beta0_ratio=1e2)
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 [27]:
print m0
print reg.mref
In [24]:
Directives.BetaSchedule??
In [23]:
# np.save('Mopt1', Mopt)
# np.save('Dest1', Dest)
In [24]:
Mopt = np.load('Mopt1.npy')
Dest = np.load('Dest1.npy')
In [25]:
# 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 [26]:
SigMat = np.exp(np.vstack(Mopt))
DpreMat = np.vstack(Dest)
In [27]:
plt.pcolor(X, Z, np.log10(SigMat))
plt.ylim(-700, 0.)
Out[27]:
In [28]:
DobsMat = Utils.mkvc(Dobs1[ind, :]).reshape((30, 31), order='F')
In [69]:
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 [54]:
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)
Out[54]:
In [ ]: