In [1]:
from SimPEG import Mesh, Utils, Maps, Survey
from SimPEG.EM.Static import DC, IP
from pymatsolver import MumpsSolver
%pylab inline
In [2]:
csx, csy, csz = 25., 25., 25.
ncx, ncy, ncz = 48, 48, 20
npad = 7
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")
In [11]:
import matplotlib
matplotlib.rcParams['font.size'] = 14
In [12]:
sigma = mesh.readModelUBC("VTKout_DC.dat")
In [13]:
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=(eta[~airind].min(), eta[~airind].max()))
# plt.axis("equal")
if normal == "Z":
xlim(-600, 600)
ylim(-600, 600.)
else:
xlim(-600, 600)
ylim(-600, 0.)
In [ ]:
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]
dat = plt.contourf(X, Y, temp, 20, clim=clim, vmin=vmin, vmax=vmax)
else:
dat = plt.contourf(X, Y, temp, 20)
plt.colorbar(dat)
In [15]:
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 [16]:
airind = sigma==1e-8
mesh2D, topoCC = gettopoCC(mesh, airind)
In [117]:
Aloc1_x = np.r_[-600., 0, 0.]
Bloc1_x = np.r_[600., 0, 0.]
Aloc2_x = np.r_[-350., 0, 0.]
Bloc2_x = np.r_[350., 0, 0.]
Aloc1_y = np.r_[0, -600., 0.]
Bloc1_y = np.r_[0, 600. , 0.]
Aloc2_y = np.r_[0, -350., 0.]
Bloc2_y = np.r_[0, 350. , 0.]
In [118]:
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 [119]:
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 [120]:
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 [121]:
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 [122]:
fig, ax = plt.subplots(1,1, figsize=(5,4))
ax.plot(Aloc1_x[0], Aloc1_x[1], 'ro', ms=5)
ax.plot(Bloc1_x[0], Bloc1_x[1], 'bo', ms=5)
ax.plot(Aloc2_x[0], Aloc2_x[1], 'ro', ms=5)
ax.plot(Bloc2_x[0], Bloc2_x[1], 'bo', ms=5)
ax.plot(Aloc1_y[0], Aloc1_y[1], 'ro', ms=5)
ax.plot(Bloc1_y[0], Bloc1_y[1], 'bo', ms=5)
ax.plot(Aloc2_y[0], Aloc2_y[1], 'ro', ms=5)
ax.plot(Bloc2_y[0], Bloc2_y[1], 'bo', ms=5)
ax.plot(My[:,0], My[:,1], 'r.', ms=2)
ax.plot(Ny[:,0], Ny[:,1], 'b.', ms=2)
# dat = mesh2D.plotImage(topoCC, ax=ax, pcolorOpts={'alpha':0.2})
# plt.colorbar(dat[0])
ax.set_xlim(-800, 800)
ax.set_ylim(-800, 800)
Out[122]:
In [132]:
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_x, Bloc0_x)
src1 = DC.Src.Dipole([rx_x, rx_y], Aloc1_x, Bloc1_x)
src2 = DC.Src.Dipole([rx_x, rx_y], Aloc2_x, Bloc2_x)
src3 = DC.Src.Dipole([rx_x, rx_y], Aloc1_y, Bloc1_y)
src4 = DC.Src.Dipole([rx_x, rx_y], Aloc2_y, Bloc2_y)
In [134]:
survey = DC.Survey([src1, src2, src3, src4])
problem = DC.Problem3D_CC(mesh)
problem.Solver = MumpsSolver
problem.pair(survey)
f = problem.fields(sigma)
dobs = survey.dpred(sigma, f=f)
In [135]:
dcdata = Survey.Data(survey, v=dobs)
Xx = 0.5*(Mx[:,0]+Nx[:,0]).reshape((23, 24), order="F")
Yx = Mx[:,1].reshape((23, 24), order="F")
Xy = My[:,0].reshape((24, 23), order="F")
Yy = 0.5*(My[:,1]+Ny[:,1]).reshape((24, 23), order="F")
In [136]:
sigopt = np.load("sigest.npy")
In [137]:
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([src1, src2, src3, src4])
problemIP.pair(surveyIP)
dataIP = surveyIP.dpred(eta[~airind])
In [138]:
# survey = DC.Survey([src0, src1, src2, src3, src4, src5])
# problem = DC.Problem3D_CC(mesh)
# problem.Solver = MumpsSolver
# problem.pair(survey)
# f = problem.fields(sigopt)
# problemIP = IP.Problem3D_CC(mesh, rho=1./sigopt, Ainv=problem.Ainv, f=f, mapping=actmapIP)
# problemIP.Solver = MumpsSolver
# surveyIP = IP.Survey([src0, src1, src2])
# problemIP.pair(surveyIP)
In [139]:
# interact(lambda ind: viz(eta, ind, normal="Z"), ind=IntSlider(min=0, max=32,step=1, value=30))
In [140]:
ipdata = Survey.Data(surveyIP, v=dataIP)
In [141]:
abs(dobs).min()
Out[141]:
In [142]:
vizdata(ipdata, src4, rx_y)
In [143]:
vizdata(ipdata, src4, rx_x)
In [146]:
a = hist(np.log10(dataIP[dataIP>0.]), bins=100)
b = hist(np.log10(-dataIP[dataIP<0.]), bins=100, color='r', alpha=0.5)
In [147]:
depth = 1./(abs(mesh.gridCC[:,2]))**1.5
depth = depth/depth.max()
In [148]:
from SimPEG import DataMisfit, Regularization, Optimization, Directives, InvProblem, Inversion
std = 0.
eps = 10**(-5.5)
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]
opt = Optimization.ProjectedGNCG(maxIter = 10)
opt.lower = 0.
invProb = InvProblem.BaseInvProblem(dmisfit, reg, opt)
# Create an inversion object
beta = Directives.BetaSchedule(coolingFactor=5, coolingRate=3)
betaest = Directives.BetaEstimate_ByEig(beta0_ratio=1.)
save = Directives.SaveOutputEveryIteration()
target = Directives.TargetMisfit()
# beta.beta = 10.
inv = Inversion.BaseInversion(invProb, directiveList=[betaest, beta, save, target])
reg.alpha_s = 1e-1
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 [151]:
plt.plot(invProb.dpred, '.', ms=1)
plt.plot(surveyIP.dobs)
Out[151]:
In [153]:
ippred = Survey.Data(surveyIP, v=invProb.dpred)
ippred_x = ippred[src1, rx_x]
ippred_y = ippred[src1, rx_y]
ip_x = ipdata[src1, rx_x]
ip_y = ipdata[src1, rx_y]
IPpredx = ippred_x.reshape((23, 24), order="F")
IPpredy = ippred_y.reshape((24, 23), order="F")
IPx = ip_x.reshape((23, 24), order="F")
IPy = ip_y.reshape((24, 23), order="F")
In [154]:
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[154]:
In [155]:
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[155]:
In [156]:
xc = opt.recall('xc')
In [157]:
temp = eta.copy()
temp[eta==0.05]=0.05
In [158]:
from ipywidgets import interact, IntSlider
In [159]:
interact(lambda ind: vizeta(temp, ind, normal="Y"), ind=IntSlider(min=0, max=40,step=1, value=34))
In [160]:
interact(lambda ind: vizeta(actmapIP*xc[6], ind, normal="Y"), ind=IntSlider(min=0, max=40,step=1, value=34))
Out[160]:
In [161]:
interact(lambda ind: vizeta(temp, ind, normal="Z"), ind=IntSlider(min=0, max=40,step=1, value=26))
Out[161]:
In [162]:
interact(lambda ind: vizeta(actmapIP*xc[9], ind, normal="Z"), ind=IntSlider(min=0, max=40,step=1, value=26))
Out[162]:
In [62]:
plt.figure(figsize=(5*2, 2.5*2))
ax = plt.subplot(111)
mesh.plotSlice(actmapIP*xc[4], normal="Y", streamOpts={'color':'w'}, gridOpts={'color':'k', 'alpha':0.1}, grid=True, ind=33,ax=ax)
mesh.plotSlice(temp, normal="Y", streamOpts={'color':'w'}, gridOpts={'color':'k', 'alpha':0.1}, grid=True, ind=33, clim=(0, 0.05),ax=ax, pcolorOpts={'cmap':'jet', 'alpha':0.3})
xlim(-600, 600)
ylim(-600, 0.)
Out[62]:
In [63]:
plt.figure(figsize=(5*2, 5*2))
ax = plt.subplot(111)
indz=25
mesh.plotSlice(actmapIP*xc[9], 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[63]:
In [ ]: