wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-07-08 15:12:41 +08:00
Code Issues Packages Projects Releases Wiki Activity

Small changes to Example

Browse Source
This commit is contained in:
D Fournier
2016-05-29 14:14:00 -07:00
parent 90fb125521
commit f1eaac7cd0
2 changed files with 26 additions and 24 deletions
Download Patch File Download Diff File Expand all files Collapse all files
simpegPF/Dev/MAG/Intgrl_MAG_Inv_ActiveCells.py
+23 -21
View File
@@ -34,6 +34,8 @@ midx = int(mesh.nCx/2)
midy = int(mesh.nCy/2)+1
midz = int(mesh.nCz/2)
vmin = 0
vmax = 1e-3
#%% Run inversion
prob = PF.Magnetics.Problem3D_Integral(mesh, mapping=idenMap, actInd=actv)
prob.solverOpts['accuracyTol'] = 1e-4
@@ -43,7 +45,7 @@ survey.pair(prob)
pred = prob.fields(driver.m0)
PF.Magnetics.writeUBCobs('Pred.dat', survey, pred)
wr = np.sum(prob.G**2.,axis=0)**0.5 / mesh.vol[actv]
wr = np.sum(prob.G**2.,axis=0)**0.5
wr = ( wr/np.max(wr) )
wr_out = actvMap * wr
@@ -56,11 +58,11 @@ plt.gca().set_aspect('equal', adjustable='box')
reg = Regularization.Simple(mesh, indActive=actv, mapping=idenMap)
reg.mref = driver.mref
reg.wght = wr
reg.cell_weights = wr
dmis = DataMisfit.l2_DataMisfit(survey)
dmis.Wd = 1/wd
opt = Optimization.ProjectedGNCG(maxIter=10,lower=0.,upper=1., maxIterCG= 20, tolCG=1e-3)
opt = Optimization.ProjectedGNCG(maxIter=10,lower=0.,upper=1., maxIterCG= 10, tolCG=1e-3)
invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
@@ -93,14 +95,14 @@ m_out[m_out==-100] = np.nan
plt.figure()
ax = plt.subplot(221)
mesh.plotSlice(m_out, ax = ax, normal = 'Z', ind=-5, clim = (mrec.min(), mrec.max()))
mesh.plotSlice(m_out, ax = ax, normal = 'Z', ind=-5, clim = (vmin,vmax))
plt.plot(np.array([mesh.vectorCCx[0],mesh.vectorCCx[-1]]), np.array([mesh.vectorCCy[yslice],mesh.vectorCCy[yslice]]),c='w',linestyle = '--')
plt.title('Z: ' + str(mesh.vectorCCz[-5]) + ' m')
plt.xlabel('x');plt.ylabel('z')
plt.gca().set_aspect('equal', adjustable='box')
ax = plt.subplot(222)
mesh.plotSlice(m_out, ax = ax, normal = 'Z', ind=-8, clim = (mrec.min(), mrec.max()))
mesh.plotSlice(m_out, ax = ax, normal = 'Z', ind=-8, clim = (vmin,vmax))
plt.plot(np.array([mesh.vectorCCx[0],mesh.vectorCCx[-1]]), np.array([mesh.vectorCCy[yslice],mesh.vectorCCy[yslice]]),c='w',linestyle = '--')
plt.title('Z: ' + str(mesh.vectorCCz[-8]) + ' m')
plt.xlabel('x');plt.ylabel('z')
@@ -108,15 +110,15 @@ plt.gca().set_aspect('equal', adjustable='box')
ax = plt.subplot(212)
mesh.plotSlice(m_out, ax = ax, normal = 'Y', ind=yslice, clim = (mrec.min(), mrec.max()))
mesh.plotSlice(m_out, ax = ax, normal = 'Y', ind=yslice, clim = (vmin,vmax))
plt.title('Smooth Unconstrained')
plt.xlabel('x');plt.ylabel('z')
plt.gca().set_aspect('equal', adjustable='box')
#%% Re-run inversion using a starting model with
# static cells
m0 = np.ones(mesh.nC)*1e-3
val = 0.005
m0 = np.ones(mesh.nC)*1e-4
val = 0.002
# Reshape the model in order to create a static block
m0 = np.reshape(m0,(mesh.nCx,mesh.nCy,mesh.nCz), order = 'F')
@@ -140,11 +142,11 @@ PF.Magnetics.writeUBCobs('Pred.dat', survey, pred)
reg = Regularization.Simple(mesh, indActive=actv, mapping=actvCells)
reg.mref = driver.mref[ind_act]
reg.wght = wr
reg.cell_weights = wr
dmis = DataMisfit.l2_DataMisfit(survey)
dmis.Wd = 1/wd
opt = Optimization.ProjectedGNCG(maxIter=10,lower=0.,upper=1., maxIterCG= 20, tolCG=1e-3)
opt = Optimization.ProjectedGNCG(maxIter=10,lower=0.,upper=1., maxIterCG= 10, tolCG=1e-3)
invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
@@ -177,21 +179,21 @@ m_out[m_out==-100] = np.nan
plt.figure()
ax = plt.subplot(221)
mesh.plotSlice(m_out, ax = ax, normal = 'Z', ind=-5, clim = (mrec.min(), mrec.max()))
mesh.plotSlice(m_out, ax = ax, normal = 'Z', ind=-5, clim = (vmin,vmax))
plt.plot(np.array([mesh.vectorCCx[0],mesh.vectorCCx[-1]]), np.array([mesh.vectorCCy[yslice],mesh.vectorCCy[yslice]]),c='w',linestyle = '--')
plt.title('Z: ' + str(mesh.vectorCCz[-5]) + ' m')
plt.xlabel('x');plt.ylabel('z')
plt.gca().set_aspect('equal', adjustable='box')
ax = plt.subplot(222)
mesh.plotSlice(m_out, ax = ax, normal = 'Z', ind=-8, clim = (mrec.min(), mrec.max()))
mesh.plotSlice(m_out, ax = ax, normal = 'Z', ind=-8, clim = (vmin,vmax))
plt.plot(np.array([mesh.vectorCCx[0],mesh.vectorCCx[-1]]), np.array([mesh.vectorCCy[yslice],mesh.vectorCCy[yslice]]),c='w',linestyle = '--')
plt.title('Z: ' + str(mesh.vectorCCz[-8]) + ' m')
plt.xlabel('x');plt.ylabel('z')
plt.gca().set_aspect('equal', adjustable='box')
ax = plt.subplot(212)
mesh.plotSlice(m_out, ax = ax, normal = 'Y', ind=yslice, clim = (mrec.min(), mrec.max()))
mesh.plotSlice(m_out, ax = ax, normal = 'Y', ind=yslice, clim = (vmin,vmax))
plt.title('Smooth Constrained')
plt.xlabel('x');plt.ylabel('z')
plt.gca().set_aspect('equal', adjustable='box')
@@ -202,17 +204,17 @@ phid = invProb.phi_d
# Set parameters for sparsity
reg = Regularization.Sparse(mesh, indActive = actv, mapping=actvCells)
reg.recModel = mrec
reg.curModel = mrec
reg.mref = driver.mref[ind_act]
reg.wght = wr
reg.cell_weights = wr
reg.eps_p = eps_p
reg.eps_q = eps_q
reg.norms = driver.lpnorms
dmis = DataMisfit.l2_DataMisfit(survey)
dmis.Wd = wd
opt = Optimization.ProjectedGNCG(maxIter=10 , lower=0.,upper=1., maxIterCG= 10, tolCG = 1e-4)
dmis.Wd = 1./wd
opt = Optimization.ProjectedGNCG(maxIter=10 , lower=0.,upper=1., maxIterCG= 20, tolCG = 1e-4)
invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta = invProb.beta)
@@ -220,7 +222,7 @@ invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta = invProb.beta)
beta = Directives.BetaSchedule(coolingFactor=1, coolingRate=1)
update_beta = Directives.Scale_Beta(tol = 0.05)
target = Directives.TargetMisfit()
IRLS =Directives.Update_IRLS( phi_m_last = phim, phi_d_last = phid )
IRLS =Directives.Update_IRLS( phi_m_last = phim )
update_Jacobi = Directives.Update_lin_PreCond(onlyOnStart=False)
save_log = Directives.SaveOutputEveryIteration()
save_log.fileName = 'LogName_blabla'
@@ -251,21 +253,21 @@ m_out[m_out==-100] = np.nan
plt.figure()
ax = plt.subplot(221)
mesh.plotSlice(m_out, ax = ax, normal = 'Z', ind=-5, clim = (mrec.min(), mrec.max()))
mesh.plotSlice(m_out, ax = ax, normal = 'Z', ind=-5, clim = (mrec.min(),mrec.max()))
plt.plot(np.array([mesh.vectorCCx[0],mesh.vectorCCx[-1]]), np.array([mesh.vectorCCy[yslice],mesh.vectorCCy[yslice]]),c='w',linestyle = '--')
plt.title('Z: ' + str(mesh.vectorCCz[-5]) + ' m')
plt.xlabel('x');plt.ylabel('z')
plt.gca().set_aspect('equal', adjustable='box')
ax = plt.subplot(222)
mesh.plotSlice(m_out, ax = ax, normal = 'Z', ind=-8, clim = (mrec.min(), mrec.max()))
mesh.plotSlice(m_out, ax = ax, normal = 'Z', ind=-8, clim = (mrec.min(),mrec.max()))
plt.plot(np.array([mesh.vectorCCx[0],mesh.vectorCCx[-1]]), np.array([mesh.vectorCCy[yslice],mesh.vectorCCy[yslice]]),c='w',linestyle = '--')
plt.title('Z: ' + str(mesh.vectorCCz[-8]) + ' m')
plt.xlabel('x');plt.ylabel('z')
plt.gca().set_aspect('equal', adjustable='box')
ax = plt.subplot(212)
mesh.plotSlice(m_out, ax = ax, normal = 'Y', ind=yslice, clim = (mrec.min(), mrec.max()))
mesh.plotSlice(m_out, ax = ax, normal = 'Y', ind=yslice, clim = (mrec.min(),mrec.max()))
plt.title('Compact Constrained')
plt.xlabel('x');plt.ylabel('z')
plt.gca().set_aspect('equal', adjustable='box')
simpegPF/Dev/MAG/Intgrl_MAG_Inv_Driver.py
+3 -3
View File
@@ -60,9 +60,9 @@ reg.mref = driver.mref
reg.cell_weights = wr
#reg.mref = np.zeros(mesh.nC)
eps_p = 2e-4
eps_p = 5e-5
eps_q = 5e-5
norms = [0., 2., 2., 2.]
norms = [0., 1., 1., 1.]
opt = Optimization.ProjectedGNCG(maxIter=100 ,lower=-2.,upper=2., maxIterLS = 20, maxIterCG= 10, tolCG = 1e-3)
invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
@@ -70,7 +70,7 @@ invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
#update_beta = Directives.Scale_Beta(tol = 0.05, coolingRate=5)
betaest = Directives.BetaEstimate_ByEig()
IRLS = Directives.Update_IRLS( norms=norms, eps_p=eps_p, eps_q=eps_q, f_min_change = 1e-2)
update_Jacobi = Directives.Update_lin_PreCond(onlyOnStart=True)
update_Jacobi = Directives.Update_lin_PreCond()
inv = Inversion.BaseInversion(invProb, directiveList=[IRLS,betaest,update_Jacobi])
# Run inversion