Fork me on GitHub

geoscience / scientific computing / visualization / web design / thoughts

simpeg / SimPEG Play - DC Resistivity Inversion

by Rowan Cockett
In [1]:
from SimPEG import *
Warning: mumps solver not available.

In [2]:
# Create the mesh
h1 = np.ones(100)
h2 = np.ones(21)
M = Mesh.TensorMesh([h1,h2])

# Create some parameters for the model
sig1 = np.log(0.001)
sig2 = np.log(1)

# Create a synthetic model from a block in a half-space
model = Model.LogModel(M)
mSynth = Utils.ModelBuilder.defineBlockConductivity(M.gridCC,[50, 10],[80, 16],[sig1, sig2])
plt.colorbar(M.plotImage(mSynth))
Tx = np.arange(10.5,90,5)
Rx = np.arange(9.5,90,3)
yLoc = 13.5
plot(Rx,np.ones_like(Rx)*yLoc,'r.',ms=20)
plot(Tx,np.ones_like(Tx)*yLoc,'c.',ms=20)
Out[2]:
[<matplotlib.lines.Line2D at 0x10f9e1510>]
In [3]:
inds = []
for ii in range(Rx.size):
    for jj in range(Rx.size - ii - 1):
        inds.append([ ii, ii + jj + 1 ])
posElectrodes = np.array([[Rx[ind[0]],yLoc] for ind in inds])
negElectrodes = np.array([[Rx[ind[1]],yLoc] for ind in inds])

Ppos = M.getInterpolationMat(posElectrodes, 'CC')
Pneg = M.getInterpolationMat(negElectrodes, 'CC')
P = Ppos - Pneg

inds = []
for ii in range(Tx.size):
    for jj in range(Tx.size - ii - 1):
        inds.append([ ii, ii + jj + 1 ])
posElectrodes = np.array([[Tx[ind[0]],yLoc] for ind in inds])
negElectrodes = np.array([[Tx[ind[1]],yLoc] for ind in inds])

Qpos = M.getInterpolationMat(posElectrodes, 'CC')
Qneg = M.getInterpolationMat(negElectrodes, 'CC')
Q = (Qpos - Qneg).T.toarray()
In [4]:
prob = Examples.DC.DCProblem(M, model)
# Create some data
data = prob.createSyntheticData(mSynth, std=0.05, P=P, RHS=Q)
u = prob.field(mSynth)
u = data.reshapeFields(u)
M.plotImage(u[:,2])
Out[4]:
<matplotlib.collections.QuadMesh at 0x10fc77310>
In [5]:
# Now set up the prob to do some minimization
# prob.dobs = dobs
# prob.std = dobs*0 + 0.05
m0 = np.zeros(M.nC) + sig2

reg = Regularization.Tikhonov(model)
beta = Parameters.BetaEstimate(beta0_ratio=0.5)
obj = ObjFunction.BaseObjFunction(data, reg, beta=beta)
opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, maxIterCG=3)
inv = Inversion.BaseInversion(obj, opt)

# Check Derivative
derChk = lambda m: [obj.dataObj(m), obj.dataObjDeriv(m)]
# Tests.checkDerivative(derChk, mSynth)
In [6]:
m = inv.run(m0)
Regularization has not set mref. SimPEG will set it to m0.
============================ Inexact Gauss Newton ============================
  #     beta     phi_d     phi_m       f      |proj(x-g)-x|  LS    Comment   
-----------------------------------------------------------------------------
   0  1.80e+02  1.25e+06  0.00e+00  1.25e+06    5.24e+03      0              
   1  1.80e+02  1.21e+06  6.15e+00  1.21e+06    6.17e+03      0              
   2  1.80e+02  1.14e+06  9.45e+01  1.16e+06    6.83e+03      0              
   3  1.80e+02  9.87e+05  1.53e+02  1.01e+06    1.69e+04      0              
   4  1.80e+02  8.33e+05  1.73e+02  8.64e+05    6.31e+04      2              
   5  1.80e+02  6.92e+05  2.17e+02  7.31e+05    6.27e+04      1              
   6  1.80e+02  5.64e+05  2.30e+02  6.05e+05    3.48e+04      0              
   7  1.80e+02  4.35e+05  2.70e+02  4.83e+05    1.48e+05      0              
   8  1.80e+02  3.55e+05  2.79e+02  4.05e+05    3.84e+04      0              
   9  1.80e+02  2.84e+05  2.82e+02  3.34e+05    5.42e+04      0              
  10  1.80e+02  2.01e+05  3.24e+02  2.60e+05    6.50e+04      0              
------------------------- STOP! -------------------------
1 : |fc-fOld| = 7.4684e+04 <= tolF*(1+|f0|) = 1.2546e+05
0 : |xc-x_last| = 6.7408e+00 <= tolX*(1+|x0|) = 1.0000e-01
0 : |proj(x-g)-x|    = 6.5036e+04 <= tolG          = 1.0000e-01
0 : |proj(x-g)-x|    = 6.5036e+04 <= 1e3*eps       = 1.0000e-02
1 : maxIter   =      10    <= iter          =     10
------------------------- DONE! -------------------------

In [7]:
plt.colorbar(M.plotImage(m))
Out[7]:
<matplotlib.colorbar.Colorbar instance at 0x111e14368>
In [7]: