diff --git a/simpegPF/Dev/MAG/Intgrl_MAG_Inv_ActiveCells.py b/simpegPF/Dev/MAG/Intgrl_MAG_Inv_ActiveCells.py index 489fcb39..ddd3df14 100644 --- a/simpegPF/Dev/MAG/Intgrl_MAG_Inv_ActiveCells.py +++ b/simpegPF/Dev/MAG/Intgrl_MAG_Inv_ActiveCells.py @@ -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') diff --git a/simpegPF/Dev/MAG/Intgrl_MAG_Inv_Driver.py b/simpegPF/Dev/MAG/Intgrl_MAG_Inv_Driver.py index a81d8b27..f3ea912b 100644 --- a/simpegPF/Dev/MAG/Intgrl_MAG_Inv_Driver.py +++ b/simpegPF/Dev/MAG/Intgrl_MAG_Inv_Driver.py @@ -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