# 2D DC inversion of Dipole Dipole arrayΒΆ

This is an example for 2D DC Inversion. The model consists of 2 cylinders, one conductive, the other one resistive compared to the background.

We restrain the inversion to the Core Mesh through the use an Active Cells mapping that we combine with an exponetial mapping to invert in log conductivity space. Here mapping, $$\mathcal{M}$$, indicates transformation of our model to a different space:

$\sigma = \mathcal{M}(\mathbf{m})$

Following example will show you how user can implement a 2D DC inversion.

from SimPEG import (
Mesh, Maps, Utils,
DataMisfit, Regularization, Optimization,
InvProblem, Directives, Inversion
)
from SimPEG.EM.Static import DC, Utils as DCUtils
import numpy as np
import matplotlib.pyplot as plt
try:
from pymatsolver import Pardiso as Solver
except ImportError:
from SimPEG import SolverLU as Solver

# Reproducible science
np.random.seed(12345)

# 2D Mesh
#########

# Cells size
csx, csz = 0.25, 0.25
# Number of core cells in each direction
ncx, ncz = 123, 41
# Vectors of cell lengthts in each direction
hz = [(csz, npad, -1.5), (csz, ncz)]
# Create mesh
mesh = Mesh.TensorMesh([hx, hz], x0="CN")
mesh.x0[1] = mesh.x0[1] + csz / 2.

# 2-cylinders Model Creation

# Spheres parameters
x0, z0, r0 = -6., -5., 3.
x1, z1, r1 = 6., -5., 3.

ln_sigback = -5.
ln_sigc = -3.
ln_sigr = -6.

mtrue = ln_sigback * np.ones(mesh.nC)

csph = (np.sqrt((mesh.gridCC[:, 1] - z0) **
2. + (mesh.gridCC[:, 0] - x0)**2.)) < r0
mtrue[csph] = ln_sigc * np.ones_like(mtrue[csph])

# Define the sphere limit
rsph = (np.sqrt((mesh.gridCC[:, 1] - z1) **
2. + (mesh.gridCC[:, 0] - x1)**2.)) < r1
mtrue[rsph] = ln_sigr * np.ones_like(mtrue[rsph])

mtrue = Utils.mkvc(mtrue)
xmin, xmax = -15., 15
ymin, ymax = -15., 0.
xyzlim = np.r_[[[xmin, xmax], [ymin, ymax]]]
actind, meshCore = Utils.meshutils.ExtractCoreMesh(xyzlim, mesh)

# Function to plot cylinder border
def getCylinderPoints(xc, zc, r):
xLocOrig1 = np.arange(-r, r + r / 10., r / 10.)
xLocOrig2 = np.arange(r, -r - r / 10., -r / 10.)
# Top half of cylinder
zLoc1 = np.sqrt(-xLocOrig1**2. + r**2.) + zc
# Bottom half of cylinder
zLoc2 = -np.sqrt(-xLocOrig2**2. + r**2.) + zc
# Shift from x = 0 to xc
xLoc1 = xLocOrig1 + xc * np.ones_like(xLocOrig1)
xLoc2 = xLocOrig2 + xc * np.ones_like(xLocOrig2)

topHalf = np.vstack([xLoc1, zLoc1]).T
topHalf = topHalf[0:-1, :]
bottomHalf = np.vstack([xLoc2, zLoc2]).T
bottomHalf = bottomHalf[0:-1, :]

cylinderPoints = np.vstack([topHalf, bottomHalf])
cylinderPoints = np.vstack([cylinderPoints, topHalf[0, :]])
return cylinderPoints

# Setup a Dipole-Dipole Survey
xmin, xmax = -15., 15.
ymin, ymax = 0., 0.
zmin, zmax = mesh.vectorCCy[-1], mesh.vectorCCy[-1]
endl = np.array([[xmin, ymin, zmin], [xmax, ymax, zmax]])
survey = DCUtils.gen_DCIPsurvey(endl, "dipole-dipole", dim=mesh.dim,
a=1, b=1, n=10, d2flag='2D')

# Setup Problem with exponential mapping and Active cells only in the core mesh
expmap = Maps.ExpMap(mesh)
mapactive = Maps.InjectActiveCells(mesh=mesh, indActive=actind,
valInactive=-5.)
mapping = expmap * mapactive
problem = DC.Problem3D_CC(mesh, sigmaMap=mapping)
problem.pair(survey)
problem.Solver = Solver

survey.dpred(mtrue[actind])
survey.makeSyntheticData(mtrue[actind], std=0.05, force=True)

# Tikhonov Inversion


Out:

SimPEG.Survey assigned new std of 5.00%

array([-13.47219189,  -5.47768312,  -3.22386406,  -2.07292378,
-1.23308792,  -0.85059664,  -0.5787703 ,  -0.36229954,
-0.2447519 ,  -0.19141787, -12.73444321,  -6.04147174,
-3.13779281,  -1.98502043,  -1.20462239,  -0.76521131,
-0.53951935,  -0.33037245,  -0.2263695 ,  -0.21597329,
-13.85673898,  -5.60392808,  -2.67590793,  -2.03228037,
-1.19056716,  -0.68006795,  -0.42362403,  -0.30014103,
-0.25029665,  -0.20369697, -14.09165599,  -6.04091436,
-3.51152453,  -1.94747423,  -1.05459764,  -0.58475779,
-0.42757887,  -0.30403622,  -0.26749497,  -0.21709242,
-15.36849737,  -5.85885279,  -3.22037025,  -1.6857607 ,
-0.96128554,  -0.55108455,  -0.38557776,  -0.31716502,
-0.24892682,  -0.20299512, -13.82008168,  -5.63581001,
-2.81195205,  -1.40411552,  -0.79153886,  -0.57119619,
-0.42185611,  -0.30313399,  -0.23859573,  -0.21806845,
-13.32936903,  -5.5922933 ,  -2.68480182,  -1.38209539,
-0.80971133,  -0.58205453,  -0.39792684,  -0.33836749,
-0.25193643,  -0.21891943, -14.91239957,  -5.24327076,
-2.59148002,  -1.42256517,  -0.85061003,  -0.60820563,
-0.43434002,  -0.33143187,  -0.26129075,  -0.21144138,
-13.94138031,  -5.90580389,  -2.68619417,  -1.56079442,
-0.94612095,  -0.60737406,  -0.49910506,  -0.34275962,
-0.29003082,  -0.23770572, -13.78637301,  -5.07199898,
-2.98580348,  -1.58087701,  -1.12804583,  -0.74286388,
-0.55807682,  -0.45214183,  -0.35452251,  -0.32039428,
-14.92500669,  -5.92994674,  -3.10539543,  -1.97766675,
-1.27813646,  -0.88120897,  -0.62911494,  -0.51917219,
-0.48955489,  -0.3863628 , -14.06169967,  -5.61277268,
-2.9863396 ,  -2.09595537,  -1.38177453,  -0.98377392,
-0.76421828,  -0.66926321,  -0.56021322,  -0.53800916,
-15.05670205,  -5.98647425,  -3.03834534,  -1.87289946,
-1.37042568,  -1.05072797,  -0.86934544,  -0.77708735,
-0.65523074,  -0.53195687, -14.20255208,  -6.05575287,
-2.7872264 ,  -1.83952827,  -1.3718127 ,  -1.09918271,
-0.89225787,  -0.82040559,  -0.66946691,  -0.52151808,
-13.79882196,  -5.100677  ,  -2.9202009 ,  -2.03247342,
-1.47363289,  -1.23680168,  -1.01483717,  -0.81388464,
-0.63146811,  -0.51017728, -13.62348916,  -5.19312665,
-3.05392374,  -2.22655011,  -1.55017882,  -1.36946746,
-1.03320042,  -0.76383402,  -0.60409095,  -0.52597475,
-13.26816975,  -5.38790974,  -3.37352274,  -1.94514285,
-1.72469688,  -1.200259  ,  -1.01664811,  -0.7482127 ,
-0.57588847,  -0.49784498, -11.45031383,  -5.58246984,
-3.09551495,  -2.36293048,  -1.79673014,  -1.32388012,
-1.04955382,  -0.8251507 ,  -0.62229811,  -0.5341308 ,
-13.19183625,  -5.59532669,  -3.45993713,  -2.09388654,
-1.58407705,  -1.19647411,  -1.04021809,  -0.81175295,
-0.64871704,  -0.53058539, -13.30784574,  -5.74482409,
-3.16122202,  -2.21926542,  -1.57737716,  -1.18232127,
-1.01487245,  -0.82503595,  -0.58641208, -12.66857223,
-5.17068845,  -3.281094  ,  -2.10533859,  -1.58899538,
-1.23682739,  -0.95137305,  -0.82930553, -12.90674906,
-5.57302385,  -3.07097984,  -1.86069478,  -1.39060005,
-1.10695264,  -0.79136702, -12.93612946,  -4.93162469,
-3.07711298,  -1.84698584,  -1.25672047,  -1.01303895,
-13.23062998,  -5.19324888,  -3.05686905,  -2.02248874,
-1.24703757, -12.3716122 ,  -5.06005593,  -2.9092524 ,
-1.87377494, -12.96908101,  -5.63768361,  -3.1949766 ,
-13.85008144,  -5.02995447, -13.54722926])

m0 = np.median(ln_sigback) * np.ones(mapping.nP)
dmis = DataMisfit.l2_DataMisfit(survey)
regT = Regularization.Simple(mesh, indActive=actind)

# Personal preference for this solver with a Jacobi preconditioner
opt = Optimization.ProjectedGNCG(maxIter=20, lower=-10, upper=10,
maxIterLS=20, maxIterCG=30, tolCG=1e-4)

opt.remember('xc')
invProb = InvProblem.BaseInvProblem(dmis, regT, opt)

beta = Directives.BetaEstimate_ByEig(beta0_ratio=1.)
Target = Directives.TargetMisfit()
betaSched = Directives.BetaSchedule(coolingFactor=5., coolingRate=2)
update_Jacobi = Directives.UpdatePreconditioner()

inv = Inversion.BaseInversion(invProb, directiveList=[beta, Target,
update_Jacobi])

minv = inv.run(m0)

# Final Plot
############

fig, ax = plt.subplots(1, 2, figsize=(12, 5))
ax = Utils.mkvc(ax)

cyl0v = getCylinderPoints(x0, z0, r0)
cyl1v = getCylinderPoints(x1, z1, r1)

clim = [(mtrue[actind]).min(), (mtrue[actind]).max()]

dat = meshCore.plotImage(((mtrue[actind])), ax=ax[0], clim=clim)
ax[0].set_title('Ground Truth')
ax[0].set_aspect('equal')

meshCore.plotImage((minv), ax=ax[1], clim=clim)
ax[1].set_aspect('equal')
ax[1].set_title('Inverted Model')

for i in range(2):
ax[i].plot(cyl0v[:, 0], cyl0v[:, 1], 'k--')
ax[i].plot(cyl1v[:, 0], cyl1v[:, 1], 'k--')

cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7])
cb = plt.colorbar(dat[0], ax=cbar_ax)
cb.set_label('ln conductivity')

cbar_ax.axis('off')

plt.show()


Out:

SimPEG.DataMisfit.l2_DataMisfit assigning default eps of 1e-5 * ||dobs||
SimPEG.InvProblem will set Regularization.mref to m0.

SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.
***Done using same Solver and solverOpts as the problem***
model has any nan: 0
=============================== Projected GNCG ===============================
#     beta     phi_d     phi_m       f      |proj(x-g)-x|  LS    Comment
-----------------------------------------------------------------------------
x0 has any nan: 0
0  1.93e+01  5.80e+03  0.00e+00  5.80e+03    5.06e+02      0
1  1.93e+01  2.22e+02  9.81e+00  4.12e+02    7.26e+01      0
------------------------- STOP! -------------------------
1 : |fc-fOld| = 0.0000e+00 <= tolF*(1+|f0|) = 5.7965e+02
1 : |xc-x_last| = 7.9396e+00 <= tolX*(1+|x0|) = 3.6689e+01
0 : |proj(x-g)-x|    = 7.2546e+01 <= tolG          = 1.0000e-01
0 : |proj(x-g)-x|    = 7.2546e+01 <= 1e3*eps       = 1.0000e-02
0 : maxIter   =      20    <= iter          =      2
------------------------- DONE! -------------------------
/Users/lindseyjh/git/simpeg/simpeg/examples/06-dc/plot_dipoledipole_2Dinversion_twocylinders.py:178: UserWarning: Matplotlib is currently using agg, which is a non-GUI backend, so cannot show the figure.
plt.show()


Total running time of the script: ( 0 minutes 14.550 seconds)

Gallery generated by Sphinx-Gallery