In [1]:
from SimPEG import Mesh, Utils, Maps, Survey
from SimPEG.EM.Static import DC, IP
from pymatsolver import MumpsSolver
%pylab inline
In [80]:
csx, csy, csz = 50., 50., 50.
ncx, ncy, ncz = 30, 30, 15
npad = 2
hx = [(csx,npad, -1.3),(csx,ncx),(csx,npad, 1.3)]
hy = [(csy,npad, -1.3),(csy,ncy),(csy,npad, 1.3)]
hz = [(csz,npad, -1.3),(csz,ncz), (csz/2.,6)]
mesh = Mesh.TensorMesh([hx, hy, hz],x0="CCN")
# sigma = mesh.readModelUBC("VTKout_DC.dat")
sigma = np.ones(mesh.nC)*1e-2
In [236]:
def vizeta(eta, ind, normal="Z"):
# sigma = np.ones(mesh.nC)*np.nan
if normal == "Z":
figsize(5, 5)
else:
figsize(5, 2.5)
# print mesh.vectorCCz[ind]
temp = eta.copy()
temp[airind] = np.nan
mesh.plotSlice(eta, ind=ind, normal=normal, grid=True, clim=(temp[~airind].min(), temp[~airind].max()))
# plt.axis("equal")
if normal == "Z":
xlim(-600, 600)
ylim(-600, 600.)
else:
xlim(-600, 600)
ylim(-600, 0.)
In [82]:
def vizdata(data, src, rx, rxcomponent="X", clim=None):
figsize(5,5)
temp = data[src, rx]
if rxcomponent=="X":
X = Xx.copy()
Y = Yx.copy()
else:
X = Xy.copy()
Y = Yy.copy()
temp = temp.reshape(X.shape, order="F")
if clim is not None:
vmin, vmax = clim[0], clim[1]
plt.contourf(X, Y, temp, 20, clim=clim, vmin=vmin, vmax=vmax)
else:
plt.contourf(X, Y, temp, 20)
In [83]:
def gettopoCC(mesh, airind):
# def gettopoCC(mesh, airind):
"""
Get topography from active indices of mesh.
"""
mesh2D = Mesh.TensorMesh([mesh.hx, mesh.hy], mesh.x0[:2])
zc = mesh.gridCC[:,2]
AIRIND = airind.reshape((mesh.vnC[0]*mesh.vnC[1],mesh.vnC[2]), order='F')
ZC = zc.reshape((mesh.vnC[0]*mesh.vnC[1], mesh.vnC[2]), order='F')
topo = np.zeros(ZC.shape[0])
topoCC = np.zeros(ZC.shape[0])
for i in range(ZC.shape[0]):
ind = np.argmax(ZC[i,:][~AIRIND[i,:]])
topo[i] = ZC[i,:][~AIRIND[i,:]].max() + mesh.hz[~AIRIND[i,:]][ind]*0.5
topoCC[i] = ZC[i,:][~AIRIND[i,:]].max()
XY = Utils.ndgrid(mesh.vectorCCx, mesh.vectorCCy)
return mesh2D, topoCC
In [84]:
airind = sigma==1e-8
mesh2D, topoCC = gettopoCC(mesh, airind)
In [85]:
Aloc0 = np.r_[-700., 0, 0.]
Bloc0 = np.r_[700., 0, 0.]
Aloc1 = np.r_[-600., 0, 0.]
Bloc1 = np.r_[600., 0, 0.]
Aloc2 = np.r_[-350., 0, 0.]
Bloc2 = np.r_[350., 0, 0.]
In [86]:
x = mesh.vectorCCx[np.logical_and(mesh.vectorCCx>-300., mesh.vectorCCx<300.)]
y = mesh.vectorCCy[np.logical_and(mesh.vectorCCy>-300., mesh.vectorCCy<300.)]
In [87]:
Mx = Utils.ndgrid(x[:-1], y, np.r_[-12.5/2.])
Nx = Utils.ndgrid(x[1:], y, np.r_[-12.5/2.])
My = Utils.ndgrid(x, y[:-1], np.r_[-12.5/2.])
Ny = Utils.ndgrid(x, y[1:], np.r_[-12.5/2.])
In [88]:
inds_Mx = Utils.closestPoints(mesh2D, Mx[:,:2])
inds_Nx = Utils.closestPoints(mesh2D, Nx[:,:2])
inds_My = Utils.closestPoints(mesh2D, My[:,:2])
inds_Ny = Utils.closestPoints(mesh2D, Ny[:,:2])
In [89]:
Mx_dr = np.c_[Mx[:,0], Mx[:,1], topoCC[inds_Mx]]
Nx_dr = np.c_[Nx[:,0], Nx[:,1], topoCC[inds_Nx]]
My_dr = np.c_[My[:,0], My[:,1], topoCC[inds_My]]
Ny_dr = np.c_[Ny[:,0], Ny[:,1], topoCC[inds_Ny]]
In [90]:
fig, ax = plt.subplots(1,1, figsize=(5,4))
ax.plot(Aloc0[0], Aloc0[1], 'ro', ms=10)
ax.plot(Bloc0[0], Bloc0[1], 'bo', ms=10)
# ax.plot(Mx[:,0], Mx[:,1], 'r.', ms=2)
# ax.plot(Nx[:,0], Nx[:,1], 'b.')
ax.plot(My[:,0], My[:,1], 'r.', ms=2)
ax.plot(Ny[:,0], Ny[:,1], 'b.', ms=2)
dat = mesh2D.plotImage(topoCC, ax=ax, grid=True)
plt.colorbar(dat[0])
ax.set_xlim(-1000, 1000)
ax.set_ylim(-1000, 1000)
Out[90]:
In [91]:
rx_x = DC.Rx.Dipole(Mx_dr, Nx_dr)
rx_y = DC.Rx.Dipole(My_dr, Ny_dr)
src0 = DC.Src.Dipole([rx_x, rx_y], Aloc0, Bloc0)
src1 = DC.Src.Dipole([rx_x, rx_y], Aloc1, Bloc1)
src2 = DC.Src.Dipole([rx_x, rx_y], Aloc2, Bloc2)
In [92]:
survey = DC.Survey([src0, src1, src2])
problem = DC.Problem3D_CC(mesh)
problem.Solver = MumpsSolver
problem.pair(survey)
f = problem.fields(sigma)
dobs = survey.dpred(sigma, f=f)
In [93]:
Mx.shape
Out[93]:
In [94]:
### dcdata = Survey.Data(survey, v=dobs)
Xx = 0.5*(Mx[:,0]+Nx[:,0]).reshape((x.size-1, x.size), order="F")
Yx = Mx[:,1].reshape((x.size-1, x.size), order="F")
Xy = My[:,0].reshape((x.size, x.size-1), order="F")
Yy = 0.5*(My[:,1]+Ny[:,1]).reshape((x.size, x.size-1), order="F")
In [95]:
# sigopt = np.load("sigest.npy")
In [100]:
eta = np.zeros(mesh.nC)
blkind = Utils.ModelBuilder.getIndicesSphere(np.r_[0., 0., -200], 100., mesh.gridCC)
eta[blkind] = 0.05
In [101]:
mesh.plotSlice(eta, normal="X", grid=True)
Out[101]:
In [102]:
# eta = mesh.readModelUBC("VTKout_eta.dat")
actmapIP = Maps.InjectActiveCells(mesh, ~airind, 0.)
problemIP = IP.Problem3D_CC(mesh, rho=1./sigma, Ainv=problem.Ainv, f=f, mapping=actmapIP)
problemIP.Solver = MumpsSolver
surveyIP = IP.Survey([src0, src1, src2])
problemIP.pair(surveyIP)
dataIP = surveyIP.dpred(eta[~airind])
In [103]:
ipdata = Survey.Data(surveyIP, v=dataIP)
In [104]:
vizdata(ipdata, src2, rx_x)
In [105]:
vizdata(ipdata, src2, rx_y, rxcomponent="Y")
In [106]:
a = hist(np.log10(dataIP[dataIP>0.]), bins=100)
b = hist(np.log10(-dataIP[dataIP<0.]), bins=100, color='r', alpha=0.5)
In [239]:
depth = 1./(abs(mesh.gridCC[:,2]))**1.5
depth = depth/depth.max()
In [ ]:
M = m.reshape((self.regmesh.nC, self.nModels), order="F")
In [240]:
from SimPEG import DataMisfit, Regularization, Optimization, Directives, InvProblem, Inversion
std = 0.
eps = 10**(-6)
surveyIP.std = std
surveyIP.eps = eps
m0 = np.ones(mesh.nC)[~airind]*1e-20
regmap = Maps.IdentityMap(nP=m0.size)
#TODO put warning when dobs is not set!
surveyIP.dobs = dataIP
dmisfit = DataMisfit.l2_DataMisfit(surveyIP)
reg = Regularization.Simple(mesh, mapping=regmap, indActive=~airind)
reg.wght = depth[~airind]
# reg.wght = np.ones_like(depth[~airind])
opt = Optimization.ProjectedGNCG(maxIter = 20)
opt.lower = 0.
invProb = InvProblem.BaseInvProblem(dmisfit, reg, opt)
# Create an inversion object
beta = Directives.BetaSchedule(coolingFactor=8, coolingRate=3)
betaest = Directives.BetaEstimate_ByEig()
save = Directives.SaveOutputEveryIteration()
target = Directives.TargetMisfit()
beta.beta = 10.
inv = Inversion.BaseInversion(invProb, directiveList=[betaest, beta, save, target])
reg.alpha_s = 1e-2
reg.alpha_x = 1.
reg.alpha_y = 1.
reg.alpha_z = 1.
problemIP.counter = opt.counter = Utils.Counter()
opt.LSshorten = 0.5
opt.remember('xc')
mIPopt = inv.run(m0)
problemIP.Ainv.clean()
In [265]:
plt.plot(invProb.dpred)
plt.plot(surveyIP.dobs, '.')
Out[265]:
In [242]:
ippred = Survey.Data(surveyIP, v=invProb.dpred)
ippred_x = ippred[src0, rx_x]
ippred_y = ippred[src0, rx_y]
ip_x = ipdata[src0, rx_x]
ip_y = ipdata[src0, rx_y]
IPpredx = ippred_x.reshape((x.size-1, x.size), order="F")
IPpredy = ippred_y.reshape((x.size, x.size-1), order="F")
IPx = ip_x.reshape((x.size-1, x.size), order="F")
IPy = ip_y.reshape((x.size, x.size-1), order="F")
In [243]:
plt.figure(figsize=(10,5))
ax0 = plt.subplot(121)
ax0.contourf(Xx, Yx, IPx, 20, vmin=IPx.min(), vmax=IPx.max(), clim=(IPx.min(), IPx.max()))
ax1 = plt.subplot(122)
ax1.contourf(Xx, Yx, IPpredx, 20, vmin=IPx.min(), vmax=IPx.max(), clim=(IPx.min(), IPx.max()))
Out[243]:
In [244]:
plt.figure(figsize=(10,5))
ax0 = plt.subplot(121)
ax0.contourf(Xy, Yy, IPy, 20, vmin=IPy.min(), vmax=IPy.max(), clim=(IPy.min(), IPy.max()))
ax1 = plt.subplot(122)
ax1.contourf(Xy, Yy, IPpredy, 20, vmin=IPy.min(), vmax=IPy.max(), clim=(IPy.min(), IPy.max()))
Out[244]:
In [245]:
xc = opt.recall('xc')
In [246]:
temp = eta.copy()
temp[eta==0.05]=0.05
In [247]:
from ipywidgets import interact, IntSlider
In [248]:
interact(lambda ind: vizeta(temp, ind, normal="Y"), ind=IntSlider(min=0, max=40,step=1, value=16))
Out[248]:
In [266]:
interact(lambda ind: vizeta(actmapIP*mIPopt, ind, normal="Y"), ind=IntSlider(min=0, max=40,step=1, value=16))
Out[266]:
In [267]:
interact(lambda ind: vizeta(actmapIP*mIPopt, ind, normal="Z"), ind=IntSlider(min=0, max=40,step=1, value=16))
In [230]:
plt.figure(figsize=(5*2, 2.5*2))
ax = plt.subplot(111)
mesh.plotSlice(actmapIP*xc[8], normal="Y", streamOpts={'color':'w'}, gridOpts={'color':'k', 'alpha':0.8}, ind=16,ax=ax)
mesh.plotSlice(temp, normal="Y", streamOpts={'color':'w'}, gridOpts={'color':'k', 'alpha':0.1}, grid=True, ind=16, clim=(0, 0.05),ax=ax, pcolorOpts={'cmap':'binary', 'alpha':0.2})
# xlim(-600, 600)
# ylim(-600, 0.)
Out[230]:
In [231]:
plt.figure(figsize=(5*2, 5*2))
ax = plt.subplot(111)
indz=20
mesh.plotSlice(actmapIP*xc[6], normal="Z", streamOpts={'color':'w'}, gridOpts={'color':'k', 'alpha':0.1}, grid=True, ind=indz,ax=ax)
mesh.plotSlice(temp, normal="Z", streamOpts={'color':'w'}, gridOpts={'color':'k', 'alpha':0.1}, grid=True, ind=indz, clim=(0, 0.05),ax=ax, pcolorOpts={'cmap':'jet', 'alpha':0.8})
xlim(-600, 600)
ylim(-600, 600.)
Out[231]:
In [ ]: