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_realistic.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 [7]:
frequency = np.r_[1., 10., 100.]
ind = np.argwhere(xyz2[:,1]==0.).flatten()
In [8]:
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]
DobsLine = Utils.mkvc(Dobs[ind,:,:,:])
np.save('bzobs_FD_realistic_line', DobsLine)
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 [11]:
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 [12]:
xcirc1, ycirc1 = circfun(-150., 0., 250., 60)
xcirc2, ycirc2 = circfun(150., 0., 250., 60)
In [14]:
DobsLine = Dobs[ind,:,:,:]
In [15]:
plt.plot(DobsLine[0,0,:,1], 'b.-')
plt.plot(DobsLine[0,1,:,1], 'k.-')
Out[15]:
In [16]:
print Dobs.shape
print DobsLine.shape
In [17]:
# Dobs = dobs.reshape((900, 2, frequency.size, 2), order='F')
# Dobsr1 = Dobs[:,0,:,0]
In [18]:
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[18]:
In [38]:
sig_test = (sigma[active])
In [39]:
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 [24]:
z = mesh.vectorCCz[active]
Frequency, X = np.meshgrid(frequency, x)
In [25]:
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[25]:
In [26]:
prb = EM.FDEM.ProblemFDEM_b(mesh, mapping=mapping, verbose=False)
In [27]:
prb.Solver = SolverLU
In [31]:
Mopt = []
Dest=[]
for i in range(r2[ind].size):
# for i in range(1):
# i = 29
txList = []
for freq in frequency:
rxoffset=r2[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,:,:,1])
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 [29]:
# 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.)
Out[29]:
In [30]:
# 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[30]:
In [34]:
np.save('Mopt2_realistic_FD', Mopt)
np.save('Dest2_realistic_FD', Dest)
In [35]:
Mopt = np.load('Mopt2_realistic_FD.npy')
Dest = np.load('Dest2_realistic_FD.npy')
In [44]:
SigMat = np.exp(np.vstack(Mopt))
DpreMat = np.vstack(Dest)
DpreMatR = DpreMat[:,[0,2,4]]
DpreMatI = DpreMat[:,[1,3,5]]
In [45]:
plt.pcolor(X, Z, np.log10(SigMat))
plt.ylim(-700, 0.)
Out[45]:
In [47]:
plt.semilogy(xyz2[ind,0], -DpreMatR[:,0], 'k')
plt.semilogy(xyz2[ind,0], -DpreMatR[:,1], 'b')
plt.semilogy(xyz2[ind,0],-DpreMatR[:,2], 'r')
plt.semilogy(xyz2[ind,0], -DobsLine[:,0,0,1], 'ko')
plt.semilogy(xyz2[ind,0], -DobsLine[:,0,1,1], 'bo')
plt.semilogy(xyz2[ind,0],-DobsLine[:,0,2,1], 'ro')
Out[47]:
In [53]:
plt.semilogy(xyz2[ind,0], -DpreMatI[:,0], 'k')
plt.semilogy(xyz2[ind,0], -DpreMatI[:,1], 'b')
plt.semilogy(xyz2[ind,0],-DpreMatI[:,2], 'r')
plt.semilogy(xyz2[ind,0], -DobsLine[:,1,0,1], 'ko')
plt.semilogy(xyz2[ind,0], -DobsLine[:,1,1,1], 'bo')
plt.semilogy(xyz2[ind,0],-DobsLine[:,1,2,1], 'ro')
Out[53]:
In [49]:
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(10., 300.)