In [62]:
import simpegDarcy as Darcy
import simpegSP as SP
from SimPEG import Mesh, Maps, Utils
from pymatsolver import PardisoSolver
%pylab inline
In [63]:
# Generate 3D mesh
csx, csy, csz = 5., 5., 2.5
ncx, ncy, ncz = 30, 30, 20
hx = [(csx, 5, -1.3), (csx, ncx), (csx, 5, 1.3)]
hy = [(csy, 5, -1.3), (csy, ncy), (csy, 5, -1.3),]
hz = [(csz, 5, -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 [64]:
mesh.plotGrid()
In [65]:
# Run simulation
h = Darcyprb.fields(K)
In [66]:
gradh = Darcyprb.gradh(h)
vel = Darcyprb.vel(h)
In [115]:
out = mesh.plotSlice(h, grid=True, ind= 15, 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[115]:
In [116]:
from SimPEG import Survey
In [117]:
out = mesh.plotSlice(np.log10(K), grid=True, ind= 15, normal="Z", clim=(-9, -5))
plt.colorbar(out[0])
Out[117]:
In [118]:
out = mesh.plotSlice(np.log10(K), grid=True, normal="Y", clim=(-9, -5))
plt.colorbar(out[0])
Out[118]:
In [119]:
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[119]:
In [72]:
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[72]:
In [73]:
fxm, fxp, fym, fyp, fzm, fzp = mesh.faceBoundaryInd
find = np.r_[fxm+fxp, fym+fyp, fzm+fzp]
vel[find] = 0.
In [74]:
out = mesh.plotSlice(vel, normal="Z", view="vec", vType="F", streamOpts={"color":"w"}, ind=15)
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[74]:
In [75]:
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[75]:
In [76]:
L0 = 1e-5
# jsCC = np.r_[Qv, Qv, Qv]*(mesh.aveF2CCV*vel)
In [77]:
# 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 [78]:
plt.plot(xyzM[:,0], xyzM[:,1], 'bo')
plt.plot(xyzN[:,0], xyzN[:,1], 'r.')
Out[78]:
In [79]:
#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 [80]:
print mesh
In [81]:
## current density J calculate...
In [82]:
prb = SP.Problem_CC(mesh, sigma=sigma, jsMap=Maps.IdentityMap(nP=mesh.nC*3), Solver=PardisoSolver)
rx = SP.Rx.Dipole(xyzN, xyzM)
src = SP.Src.StreamingCurrents([rx], L=np.ones(mesh.nC)*L0, mesh=mesh, modelType="CurrentDensity")
survey = SP.Survey([src])
survey.pair(prb)
In [83]:
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[83]:
In [84]:
# # 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 [85]:
out = Utils.plot2Ddata(xyzN, dobs*1e3)
plt.colorbar(out[0])
Out[85]:
In [86]:
depthweight = 1./ ((abs(mesh.gridCC[:,2]-mesh.vectorCCz[-1])+0.5)**0.8)
depthweight /= depthweight.max()
In [87]:
out = mesh.plotSlice(np.log10(depthweight), grid=True, normal="Y", ind=8)
plt.colorbar(out[0])
Out[87]:
In [88]:
out = hist(np.log10(abs(dobs)), bins=100)
In [89]:
from SimPEG.EM.Static.SIP.Regularization import MultiRegularization
In [93]:
from SimPEG import (Mesh, Maps, Utils, DataMisfit, Regularization,
Optimization, Inversion, InvProblem, Directives)
survey.std = 0.
survey.eps = abs(dobs).max() * 0.03
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/depthweight
reg.wght = depthweight
reg.alpha_s = 1.
reg.alpha_x = 1.
reg.alpha_y = 1.
reg.alpha_z = 1.
opt = Optimization.ProjectedGNCG(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-8
# 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 [94]:
xc = opt.recall("xc")
In [95]:
plt.plot(dobs)
plt.plot(invProb.dpred)
Out[95]:
In [96]:
# 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 [97]:
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 [98]:
0.00001
Out[98]:
In [99]:
# 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 [100]:
print mopt.min()
print mopt.max()
In [101]:
mesh.vectorCCz
Out[101]:
In [102]:
out = mesh.plotSlice(mopt, view='vec', vType="CCv", ind=10)
plt.colorbar(out[0])
plt.xlim(-80, 80)
plt.ylim(-80, 80)
mesh.plotSlice(mopt, view='vec', vType="CCv", normal="Y", clim=out[0].get_clim(), ind=5)
plt.colorbar(out[0])
plt.xlim(-80, 80)
plt.ylim(0, 44)
Out[102]:
In [103]:
# plt.plot(survey.dobs)
# plt.plot(invProb.dpred)
In [106]:
out = Utils.plot2Ddata(xyzN, invProb.dpred*1e3)
plt.colorbar(out[0])
Out[106]:
In [ ]:
In [ ]: