In [1]:
import simpegDarcy as Darcy
import simpegSP as SP
from SimPEG import Mesh, Maps, Utils
from pymatsolver import PardisoSolver
%pylab inline
In [2]:
!pwd
In [3]:
# Generate 3D mesh
csx, csy, csz = 5., 5., 2.5
ncx, ncy, ncz = 30, 30, 20
npad = 8
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-4), (csz/2., 4)]
mesh = Mesh.TensorMesh([hx, hy, hz], "CC0")
mesh._x0 = np.r_[mesh.x0[0], mesh.x0[1], -mesh.hz[:5].sum()]
# Generate Darcy problem
Darcyprb = Darcy.Problem_CC(mesh, KMap=Maps.IdentityMap(mesh))
# Set boundary condition
bc = [["dirichlet", "dirichlet"], ["neumann", "neumann"], ["neumann", "neumann"]]
hbc = [[50., 30.] ,[0., 0.], [0., 0.]]
Darcyprb.setBC(bc, hbc)
# Set Darcy survey
locs = Utils.ndgrid(mesh.vectorCCx, mesh.vectorCCy, np.r_[mesh.vectorCCz[-1]])
rx = Darcy.DarcyRx(locs)
Darcysurvey = Darcy.DarcySurvey([rx])
# Make hydraulic conductivity model (m/s)
K = np.ones(mesh.nC)*1e-7
layerind1 = np.logical_and(mesh.gridCC[:,2]>=20., mesh.gridCC[:,2]<30.)
blkind1 = np.logical_and(mesh.gridCC[:,0]>-30., mesh.gridCC[:,0]<30.) & np.logical_or(mesh.gridCC[:,1]<-30., mesh.gridCC[:,1]>30.)
K[layerind1] = 1e-5
K[blkind1 & layerind1] = 1e-9
# Pair the survey to the problem
Darcysurvey.pair(Darcyprb)
Darcyprb.Solver = PardisoSolver
In [4]:
# Run simulation
h = Darcyprb.fields(K)
In [5]:
gradh = Darcyprb.gradh(h)
vel = Darcyprb.vel(h)
In [6]:
out = mesh.plotSlice(h, grid=True, ind= 11, normal="Z")
plt.colorbar(out[0])
# out = mesh.plotSlice(gradh, grid=False, normal="Z", view="vec", vType="F", ind=15)
# plt.colorbar(out[0])
# out = mesh.plotSlice(gradh, grid=False, normal="Y", view="vec", vType="F", clim=out[0].get_clim(), ind =8)
# plt.colorbar(out[0])
Out[6]:
In [7]:
from SimPEG import Survey
In [8]:
out = mesh.plotSlice(np.log10(K), grid=True, ind= 11, normal="Z", clim=(-9, -5))
plt.colorbar(out[0])
Out[8]:
In [9]:
out = mesh.plotSlice(np.log10(K), grid=True, normal="Y", clim=(-9, -5))
plt.colorbar(out[0])
Out[9]:
In [10]:
p = Darcyprb.p(h)
p[p<0.] = np.nan
out = mesh.plotSlice(p, grid=True, normal="Y", clim=(0, 50))
plt.colorbar(out[0])
Out[10]:
In [11]:
out = mesh.plotSlice(gradh, grid=False, normal="Z", view="vec", vType="F", ind=11)
plt.colorbar(out[0])
out = mesh.plotSlice(gradh, grid=False, normal="Y", view="vec", vType="F", clim=out[0].get_clim(), ind =8)
plt.colorbar(out[0])
Out[11]:
In [12]:
fxm, fxp, fym, fyp, fzm, fzp = mesh.faceBoundaryInd
find = np.r_[fxm+fxp, fym+fyp, fzm+fzp]
vel[find] = 0.
In [13]:
mesh.vectorCCz[11]
Out[13]:
In [14]:
out = mesh.plotSlice(vel, normal="Z", view="vec", vType="F", streamOpts={"color":"w"}, ind=11)
plt.colorbar(out[0])
mesh.plotSlice(vel, normal="Y", view="vec", vType="F", clim=out[0].get_clim(), streamOpts={"color":"w"})
plt.colorbar(out[0])
Out[14]:
In [15]:
out = mesh.plotSlice(Darcyprb.divgradh(h), grid=True, normal="Z", ind=15)
plt.colorbar(out[0])
mesh.plotSlice(Darcyprb.divgradh(h), grid=True, normal="Y", clim=out[0].get_clim())
plt.colorbar(out[0])
Out[15]:
In [16]:
L0 = 1e-5
# jsCC = np.r_[Qv, Qv, Qv]*(mesh.aveF2CCV*vel)
In [17]:
# Make SP survey
dx = 5.
x = mesh.vectorCCx[np.logical_and(mesh.vectorCCx<80., mesh.vectorCCx>-80.)]
y = mesh.vectorCCy[np.logical_and(mesh.vectorCCy<80., mesh.vectorCCy>-80.)]
xyzM = Utils.ndgrid(x*0., y*0., np.r_[mesh.vectorCCz[-1]])
xyzN = Utils.ndgrid(x, y, np.r_[mesh.vectorCCz[-1]])
In [18]:
plt.plot(xyzM[:,0], xyzM[:,1], 'bo')
plt.plot(xyzN[:,0], xyzN[:,1], 'r.')
Out[18]:
In [19]:
#SP problem with Head map
sigma = np.ones(mesh.nC)*1e-3
sigma[layerind1] = 1e-3
sigma[blkind1 & layerind1] = 1e-1
prb = SP.Problem_CC(mesh, sigma=sigma, hMap=Maps.IdentityMap(mesh), Solver=PardisoSolver)
rx = SP.Rx.Dipole(xyzN, xyzM)
src = SP.Src.StreamingCurrents([rx], L=np.ones(mesh.nC)*L0, mesh=mesh, modelType="Head")
survey = SP.Survey([src])
survey.pair(prb)
dobs = survey.dpred(h)
q = src.eval(prb)
In [20]:
#print mesh
In [21]:
wires = Maps.Wires(('jsx', mesh.nC), ('jsy', mesh.nC), ('jsz', mesh.nC))
prb = SP.Problem_CC(mesh, sigma=sigma, jsxMap=wires.jsx, jsyMap=wires.jsy, jszMap=wires.jsz, Solver=PardisoSolver)
rx = SP.Rx.Dipole(xyzN, xyzM)
src = SP.Src.StreamingCurrents([rx], L=np.ones(mesh.nC), mesh=mesh, modelType="CurrentDensity")
survey = SP.Survey([src])
survey.pair(prb)
In [22]:
out = mesh.plotSlice(q, grid=True, normal="Z", ind=15)
plt.colorbar(out[0])
mesh.plotSlice(q, grid=True, normal="Y", clim=out[0].get_clim(), ind=5)
plt.colorbar(out[0])
Out[22]:
In [23]:
# # Generate Full sensitivity
# I = np.diag(np.ones_like(dobs))
# J = np.zeros((dobs.size, mesh.nC))
# for i in range(dobs.size):
# J[i,:] = prb.Jtvec(sigma, I[:,i])
# JtJ = (J**2).sum(axis=0)
# JtJ /= JtJ.max()
# prb.G = J
In [24]:
out = Utils.plot2Ddata(xyzN, dobs*1e3)
plt.colorbar(out[0])
Out[24]:
In [25]:
depthweight = 1./ ((abs(mesh.gridCC[:,2]-mesh.vectorCCz[-1])+0.5)**1.)
depthweight /= depthweight.max()
In [26]:
out = mesh.plotSlice(np.log10(depthweight), grid=True, normal="Y", ind=8)
plt.colorbar(out[0])
Out[26]:
In [27]:
out = hist(np.log10(abs(dobs)), bins=100)
In [28]:
from SimPEG.EM.Static.SIP.Regularization import MultiRegularization
In [29]:
from SimPEG import (Mesh, Maps, Utils, DataMisfit, Regularization,
Optimization, Inversion, InvProblem, Directives)
survey.std = 0.
survey.eps = abs(dobs).max() * 0.05
survey.dobs = dobs
dmisfit = DataMisfit.l2_DataMisfit(survey)
regmap = Maps.IdentityMap(nP = mesh.nC*3)
reg = MultiRegularization(mesh, mapping=regmap ,nModels=3)
reg.cell_weights = depthweight
reg.alpha_s = 1.
reg.alpha_x = 0.
reg.alpha_y = 0.
reg.alpha_z = 0.
opt = Optimization.InexactGaussNewton(maxIter=100, tolX=1e-20, tolF=1e-20)
opt.maxIterLS = 20
IRLS = Directives.Update_IRLS(norms=([0.,1.,1.,1.]),
eps=None, f_min_change=1e-3,
minGNiter=3)
# senseweight = Directives.Update_Wj()
invProb = InvProblem.BaseInvProblem(dmisfit, reg, opt)
target = Directives.TargetMisfit()
# Create an inversion object
beta = Directives.BetaSchedule(coolingFactor=5, coolingRate=3)
betaest = Directives.BetaEstimate_ByEig()
betaest.beta0_ratio = 1e-10
# updateprecond = Directives.Update_lin_PreCond()
# inv = Inversion.BaseInversion(invProb, directiveList=[betaest, updateprecond, IRLS])
inv = Inversion.BaseInversion(invProb, directiveList=[beta, betaest, target])
prb.counter = opt.counter = Utils.Counter()
opt.LSshorten = 0.5
opt.remember('xc')
m0 = np.ones(mesh.nC*3)*0.
reg.mref = m0
mopt = inv.run(m0)
In [30]:
xc = opt.recall("xc")
In [31]:
plt.plot(dobs)
plt.plot(invProb.dpred)
Out[31]:
In [32]:
# fig = plt.figure(figsize = (12, 4))
# ax = plt.subplot(111)
# ax_1 = ax.twinx()
# out = ax.hist(abs(reg.l2model), bins=100)
# ax.set_xscale("linear")
# ax.set_yscale("log")
# temp = np.sort(abs(reg.l2model))
# ax_1.plot(temp, temp / (temp**2 + (reg.eps_p)**2))
# plt.xlim(0, 0.0006)
In [33]:
fig = plt.figure(figsize = (12, 4))
ax = plt.subplot(111)
ax_1 = ax.twinx()
out = ax.hist(abs(mopt), bins=100)
ax.set_xscale("linear")
ax.set_yscale("log")
temp = np.sort(abs(mopt))
# ax_1.plot(temp, temp / (temp**2 + (reg.eps_p)**2))
# plt.xlim(0, 0.0006)
In [34]:
0.00001
Out[34]:
In [35]:
# out = mesh.plotSlice(reg.l2model, grid=True, normal="Z", ind=15, clim=(-0.0001, 0.0001))
# plt.colorbar(out[0])
# mesh.plotSlice(reg.l2model, grid=True, normal="Y", clim=out[0].get_clim(), ind=5)
# plt.colorbar(out[0])
In [36]:
print mopt.min()
print mopt.max()
In [37]:
mesh.vectorCCz
Out[37]:
In [38]:
mesh.vnC[2]-1
Out[38]:
In [39]:
out = mesh.plotSlice(mopt, view='vec', vType="CCv", ind=11)
plt.colorbar(out[0])
plt.xlim(-80, 80)
plt.ylim(-80, 80)
mesh.plotSlice(mopt, view='vec', vType="CCv", normal="Y", ind=15)
plt.colorbar(out[0])
plt.xlim(-80, 80)
plt.ylim(0, 55)
mesh.plotSlice(mopt, view='vec', vType="CCv", normal="X", ind=15)
plt.colorbar(out[0])
plt.xlim(-80, 80)
plt.ylim(0, 55)
Out[39]:
In [40]:
# plt.plot(survey.dobs)
# plt.plot(invProb.dpred)
In [41]:
out = Utils.plot2Ddata(xyzN, invProb.dpred*1e3)
plt.colorbar(out[0])
Out[41]: