FDEM Inversion for Conductivity

Forward simulation and inversion for a block (fracture) on a whole space



In [15]:

    
from SimPEG import *
import simpegEM as EM
from pymatsolver import MumpsSolver
from scipy.constants import mu_0

Set up Mesh



In [2]:

    
cs = 5
ncore = 35
ncore_x = 15
pad = 8
padfactor = 1.4
hx = [(cs,pad,-padfactor),(cs,ncore_x),(cs,pad,padfactor)]
h = [(cs,pad,-padfactor),(cs,ncore),(cs,pad,padfactor)]
mesh = Mesh.TensorMesh([hx, h, h],'CCC')



In [3]:

    
sigwh = 1e-2 # conductivity of the wholespace background
sigma = 1e-2*np.ones(prod(mesh.nC))  # whole space
sigmaB = 1e-2*np.ones(prod(mesh.nC)) # with block 
sigmaB[(np.abs(mesh.gridCC[:,0]) < 10.) & (np.abs(mesh.gridCC[:,1]) < 25.) & (np.abs(mesh.gridCC[:,2]) < 25.)] = 1
mesh.plotSlice(sigmaB,grid=True) # todo: fix axes
mesh.plotImage(sigmaB)           #









    



/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/matplotlib/lines.py:503: RuntimeWarning: invalid value encountered in greater_equal
  return np.alltrue(x[1:] - x[0:-1] >= 0)






    Out[3]:





<matplotlib.collections.QuadMesh at 0x10a2c6790>

Comparing against analytic

Start by comparing the solution for the whole-space model to what we expect for a dipole in a whole-space



In [4]:

    
x = np.linspace(-50,50,10)
XYZ = Utils.ndgrid(np.r_[0],np.r_[50],x)

Rx0 = EM.FDEM.RxFDEM(XYZ, 'bzr')
Rx1 = EM.FDEM.RxFDEM(XYZ, 'bzi')

txList = []
freq = 500.
Tx = EM.FDEM.TxFDEM(np.r_[0.,0.,0.], 'VMD',freq,[Rx0,Rx1])
txList.append(Tx)

survey = EM.FDEM.SurveyFDEM(txList)
mapping = Maps.ExpMap(mesh)
prb = EM.FDEM.ProblemFDEM_b(mesh, mapping=mapping)
prb.pair(survey)
prb.Solver = MumpsSolver



In [5]:

    
# SimPEG soln
fTest = prb.fields(np.log(sigma))



In [61]:

    
Bx, By, Bz = EM.Analytics.FDEM.AnalyticMagDipoleWholeSpace(Rx0.locs,Tx.loc,sigwh,freq, m = 1.0, orientation='Z')
Bx0, By0, Bz0 = EM.Analytics.FDEM.AnalyticMagDipoleWholeSpace(Rx0.locs,Tx.loc,sigwh,0.0, m = 1.0, orientation='Z')
BxA = Bx - Bx0
ByA = By - By0
BzA = Bz - Bz0
# XYZ2 = Utils.ndgrid(np.r_[0],np.r_[100],x)
# Bx, By, Bz = EM.Analytics.FDEM.AnalyticMagDipoleWholeSpace(np.r_[0,100,0],Tx.loc,sigwh,100., m = 1.0, orientation='Z')



In [66]:

    
plot(x, np.c_[BzA.real/data[Tx,Rx0]])









    Out[66]:





[<matplotlib.lines.Line2D at 0x114c64f90>]



In [53]:

    
plot(x, np.c_[np.arctan2(Bz.imag,Bz.real)*180./np.pi, np.arctan2(data[Tx,Rx1],data[Tx,Rx0])*180./np.pi])
legend(('ana','data'))









    Out[53]:





<matplotlib.legend.Legend at 0x110a13a90>



In [9]:

    
data = survey.projectFields(fTest)



In [19]:

    
plot(x, np.c_[data[Tx,Rx0], Bz.real])
legend(('data','Ana'))









    Out[19]:





<matplotlib.legend.Legend at 0x10a9688d0>



In [48]:

    
fTest









    Out[48]:





<simpegEM.FDEM.SurveyFDEM.FieldsFDEM at 0x122bed68>



In [ ]:



In [16]:

    
x = np.linspace(-50,50,10)
XYZ = Utils.ndgrid(np.r_[0],np.r_[50],x)
Rx0 = EM.FDEM.RxFDEM(XYZ, 'bzr')
Rx1 = EM.FDEM.RxFDEM(XYZ, 'bzi')
txList = []
for z in np.linspace(-50,50,10):
    for freq in [300.,600.,900.]:
        Tx = EM.FDEM.TxFDEM(np.r_[0.,-50.,z], 'VMD', freq, [Rx0,Rx1])
        txList.append(Tx)
survey = EM.FDEM.SurveyFDEM(txList)
mapping = Maps.ExpMap(mesh)
prb = EM.FDEM.ProblemFDEM_b(mesh, mapping=mapping)
prb.pair(survey)
prb.Solver = MumpsSolver



In [41]:

    
fB = prb.fields(np.log(sigmaB))



In [ ]:

    
dB = survey.projectFields(fB)



In [71]:

    
f300 = dB[survey.getTransmitters(300.)[-1],Rx0]
f600 = dB[survey.getTransmitters(600.)[-1],Rx0]
f900 = dB[survey.getTransmitters(900.)[-1],Rx0]

semilogy(x, np.abs(np.c_[f300,f600,f900]))
legend(('3','6','9'))
# plot(x, np.c_[dpred[Tx0,Rx0], dB[Tx0,Rx0]])









    Out[71]:





<matplotlib.legend.Legend at 0x10da95950>



In [80]:

    
survey.dobs = survey.dpred(None, u=fB)
survey.std = 0.01



In [81]:

    
reg = Regularization.Tikhonov(mesh)
dmis = DataMisfit.l2_DataMisfit(survey)
opt = Optimization.InexactGaussNewton(maxIter=3)
invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
beta = Directives.BetaSchedule()
betaest = Directives.BetaEstimate_ByEig()
inv = Inversion.BaseInversion(invProb, directiveList=[beta, betaest])
m0 = np.ones_like(sigmaB)*1e-2



In [82]:

    
inv.run(np.log(m0))









    



SimPEG.InvProblem will set Regularization.mref to m0.
SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.
                    ***Done using same solver as the problem***
SimPEG.l2_DataMisfit is creating default weightings for Wd.
============================ Inexact Gauss Newton ============================
  #     beta     phi_d     phi_m       f      |proj(x-g)-x|  LS    Comment   
-----------------------------------------------------------------------------
   0  4.46e+03  6.92e+04  0.00e+00  6.92e+04    2.31e+04      0              
   1  4.46e+03  6.51e+03  2.09e+00  1.58e+04    3.71e+02      0              
   2  4.46e+03  6.15e+03  2.16e+00  1.58e+04    2.17e+01      0              
   3  5.57e+02  6.13e+03  2.16e+00  7.33e+03    3.68e+03      0   Skip BFGS  
------------------------- STOP! -------------------------
0 : |fc-fOld| = 8.4323e+03 <= tolF*(1+|f0|) = 6.9198e+03
1 : |xc-x_last| = 4.2277e-02 <= tolX*(1+|x0|) = 8.0690e+01
0 : |proj(x-g)-x|    = 3.6813e+03 <= tolG          = 1.0000e-01
0 : |proj(x-g)-x|    = 3.6813e+03 <= 1e3*eps       = 1.0000e-02
1 : maxIter   =       3    <= iter          =      3
------------------------- DONE! -------------------------






    Out[82]:





array([-4.62771738, -4.62789076, -4.62806113, ..., -4.62159749,
       -4.62166636, -4.62176178])



In [83]:

    
mopt = _



In [90]:

    
colorbar(mesh.plotSlice(mapping*mopt)[0])









    Out[90]:





<matplotlib.colorbar.Colorbar instance at 0x1141e51b8>



In [85]: