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_realistic.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 [7]:
frequency = np.r_[1., 10., 100.]
In [9]:
dobs = np.load('bzobs_FD_realistic.npy')
Dobs = dobs.reshape((900, 2, frequency.size, 2), order='F')
Dobsr1 = Dobs[:,0,:,0]
Dobsi1 = Dobs[:,1,:,0]
Dobsr2 = Dobs[:,0,:,1]
Dobsi2 = Dobs[:,1,:,1]
In [10]:
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 [11]:
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 [12]:
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 [13]:
xc = -100
yc = 100.
dind = np.argmin(abs( xyz1[:,0]-xc)+abs( xyz1[:,1]-yc))
In [14]:
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 [15]:
xcirc1, ycirc1 = circfun(-150., 0., 250., 60)
xcirc2, ycirc2 = circfun(150., 0., 250., 60)
In [31]:
ind = np.argwhere(xyz1[:,1]==0.).flatten()
In [32]:
DobsLine = Dobs[ind,:,:,:]
In [77]:
plt.plot(DobsLine[0,0,:,1], 'b.-')
plt.plot(DobsLine[0,1,:,1], 'k.-')
Out[77]:
In [34]:
print Dobs.shape
print DobsLine.shape
In [ ]:
# Dobs = dobs.reshape((900, 2, frequency.size, 2), order='F')
# Dobsr1 = Dobs[:,0,:,0]
In [17]:
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[17]:
In [18]:
sig_test = (sigma[active])
In [19]:
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 [20]:
z = mesh.vectorCCz[active]
Time, Xtime = np.meshgrid(time, x)
In [21]:
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[21]:
In [44]:
test = survey.dpred(np.log(sigma[active]))
In [22]:
prb = EM.FDEM.ProblemFDEM_b(mesh, mapping=mapping, verbose=False)
In [20]:
prb.Solver = SolverLU
In [137]:
# Mopt = []
# Dest=[]
# for i in range(r1[ind].size):
# # for i in range(1):
# # i = 29
# txList = []
# for freq in frequency:
# rxoffset=r1[ind][i]
# rxr = EM.FDEM.RxFDEM(np.array([[rxoffset, 0., 0.]]),'bzr')
# rxi = EM.FDEM.RxFDEM(np.array([[rxoffset, 0., 0.]]),'bzi')
# tx = EM.FDEM.TxFDEM(np.array([0., 0., 0.]), 'CircularLoop', freq, [rxr, rxi])
# tx.radius = 250.
# txList.append(tx)
# survey = EM.FDEM.SurveyFDEM(txList)
# if prb.ispaired:
# prb.unpair()
# if survey.ispaired:
# survey.unpair()
# prb.pair(survey)
# std = 0.2
# survey.dobs = Utils.mkvc(DobsLine[i,:,:,0])
# survey.std = survey.dobs*0 + std
# dmisfit = DataMisfit.l2_DataMisfit(survey)
# dmisfit.Wd = 1/(abs(survey.dobs)*std+1e-12)
# 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 [ ]:
# fig, ax = plt.subplots(1,1, figsize = (3, 6))
# Utils1D.plotLayer(np.exp(mopt), mesh.vectorCCz[active], showlayers=True, ax = ax)
# ax.set_ylim(-500., 0.)
In [99]:
# plt.semilogy(abs(invProb.dpred.reshape((2,3), order='F')[0,:]), 'b.-')
# plt.semilogy(abs(invProb.dpred.reshape((2,3), order='F')[1,:]), 'r.-')
# plt.semilogy(abs(survey.dobs.reshape((2,3), order='F')[0,:]), 'bo')
# plt.semilogy(abs(survey.dobs.reshape((2,3), order='F')[1,:]), 'ro')
Out[99]:
In [138]:
# np.save('Mopt1_realistic_FD', Mopt)
# np.save('Dest1_realistic_FD', Dest)
In [115]:
Mopt = np.load('Mopt1_realistic_FD.npy')
Dest = np.load('Dest1_realistic_FD.npy')
In [123]:
SigMat = np.exp(np.vstack(Mopt))
DpreMat = np.vstack(Dest)
DpreMatR = DpreMat[:,[0,2,4]]
DpreMatI = DpreMat[:,[1,3,5]]
In [117]:
plt.pcolor(X, Z, np.log10(SigMat))
plt.ylim(-700, 0.)
Out[117]:
In [131]:
plt.semilogy(xyz1[ind,0], -DpreMatR[:,0], 'k')
plt.semilogy(xyz1[ind,0], -DpreMatR[:,1], 'b')
plt.semilogy(xyz1[ind,0],-DpreMatR[:,2], 'r')
plt.semilogy(xyz1[ind,0], -DobsLine[:,0,0,0], 'ko')
plt.semilogy(xyz1[ind,0], -DobsLine[:,0,1,0], 'bo')
plt.semilogy(xyz1[ind,0],-DobsLine[:,0,2,0], 'ro')
Out[131]:
In [132]:
plt.semilogy(xyz1[ind,0], -DpreMatI[:,0], 'k')
plt.semilogy(xyz1[ind,0], -DpreMatI[:,1], 'b')
plt.semilogy(xyz1[ind,0],-DpreMatI[:,2], 'r')
plt.semilogy(xyz1[ind,0], -DobsLine[:,1,0,0], 'ko')
plt.semilogy(xyz1[ind,0], -DobsLine[:,1,1,0], 'bo')
plt.semilogy(xyz1[ind,0],-DobsLine[:,1,2,0], 'ro')
Out[132]:
In [136]:
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].pcolor(X, Z, np.log10(SigMat), 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., -10.)
ax[i].set_xlim(-300., -10.)
In [ ]: