\[\renewcommand{\div}{\nabla\cdot\,} \newcommand{\grad}{\vec \nabla} \newcommand{\curl}{{\vec \nabla}\times\,} \newcommand{\dcurl}{{\mathbf C}} \newcommand{\dgrad}{{\mathbf G}} \newcommand{\Acf}{{\mathbf A_c^f}} \newcommand{\Ace}{{\mathbf A_c^e}} \renewcommand{\S}{{\mathbf \Sigma}} \renewcommand{\Div}{{\mathbf {Div}}} \renewcommand{\Grad}{{\mathbf {Grad}}} \newcommand{\St}{{\mathbf \Sigma_\tau}} \newcommand{\diag}{\mathbf{diag}} \newcommand{\M}{{\mathbf M}} \newcommand{\Me}{{\M^e}} \newcommand{\Mes}[1]{{\M^e_{#1}}} \newcommand{\be}{\mathbf{e}} \newcommand{\bj}{\mathbf{j}} \newcommand{\bphi}{\mathbf{\phi}} \newcommand{\bq}{\mathbf{q}} \newcommand{\bJ}{\mathbf{J}} \newcommand{\bG}{\mathbf{G}} \newcommand{\bP}{\mathbf{P}} \newcommand{\bA}{\mathbf{A}} \newcommand{\bm}{\mathbf{m}} \newcommand{\B}{\vec{B}} \newcommand{\D}{\vec{D}} \renewcommand{\H}{\vec{H}} \renewcommand {\j} { {\vec j} } \newcommand {\h} { {\vec h} } \renewcommand {\b} { {\vec b} } \newcommand {\e} { {\vec e} } \newcommand {\c} { {\vec c} } \renewcommand {\d} { {\vec d} } \renewcommand {\u} { {\vec u} } \newcommand{\I}{\vec{I}}\]

Direct Current Resistivity

SimPEG.EM.Static.DC and SimPEG.EM.Static.IP uses SimPEG as the framework for the forward and inverse direct current (DC) resistivity and induced polarization (IP) geophysical problems.

DC resistivity survey

Electrical resistivity of subsurface materials is measured by causing an electrical current to flow in the earth between one pair of electrodes while the voltage across a second pair of electrodes is measured. The result is an “apparent” resistivity which is a value representing the weighted average resistivity over a volume of the earth. Variations in this measurement are caused by variations in the soil, rock, and pore fluid electrical resistivity. Surveys require contact with the ground, so they can be labour intensive. Results are sometimes interpreted directly, but more commonly, 1D, 2D or 3D models are estimated using inversion procedures (GPG).

Background

As direct current (DC) implies, in DC resistivity survey, we assume steady-state. We consider Maxwell’s equations in steady state as

\[ \begin{align}\begin{aligned}\begin{split}\curl \frac{1}{\mu} \vec{b} - \j = \j_s \\\end{split}\\\curl \e = 0\end{aligned}\end{align} \]

Then by taking \(\div\) of the first equation, we have

\[\begin{split}- \div\j = q \\\end{split}\]

where

\[\div \j_s = q = I(\delta(\vec{r}-\vec{r}_{s+})-\delta(\vec{r}-\vec{r}_{s-}))\]

Since \(\curl \e = 0\), we have

\[\e = \grad \phi\]

And by Ohm’s law, we have

\[\j = \sigma \grad \phi\]

Finally, we can compute the solution of the system:

\[ \begin{align}\begin{aligned}- \div\j = q\\\j = \sigma \grad \phi\\\frac{\partial \phi}{\partial r}\Big|_{\partial \Omega_{BC}} = 0\end{aligned}\end{align} \]

Discretization

By using finite volume method (FVM), we discretize our system as

\[ \begin{align}\begin{aligned}-\Div \bj = \bq\\\diag(\Acf^{T}\sigma^{-1}) \bj = \Grad \bphi\end{aligned}\end{align} \]

Here boundary condtions are embedded in the discrete differential operators. With some linear algebra we have

\[\bA\bphi = -\bq\]

where

\[\bA = \Div (\diag(\Acf^{T}\sigma^{-1}))^{-1} \Grad\]

By solving this linear equation, we can compute the solution of \(\phi\). Based on this discretization, we derive sensitivity in discretized space. Sensitivity matrix can be in general can be written as

\[\bJ = -\bP\bA^{-1}\bG\]

where

\[ \begin{align}\begin{aligned}\bP: \text{Projection}\\\bJ = \bP\frac{\partial \phi}{\partial \bm}\end{aligned}\end{align} \]

Here \(\bm\) indicates model parameters in discretized space.

../../_images/sphx_glr_plot_analytic_dipole_001.png

API for DC codes

Problem

class SimPEG.EM.Static.DC.ProblemDC.BaseDCProblem(mesh, **kwargs)[source]

Bases: SimPEG.EM.Base.BaseEMProblem

Base DC Problem

Optional Properties:

  • model (Model): Inversion model., a numpy array of <type ‘float’>, <type ‘int’> with shape (*)
  • mu (PhysicalProperty): Magnetic Permeability (H/m), a physical property, Default: 1.25663706144e-06
  • mui (PhysicalProperty): Inverse Magnetic Permeability (m/H), a physical property
  • rho (PhysicalProperty): Electrical resistivity (Ohm m), a physical property
  • rhoMap (Mapping): Mapping of Electrical resistivity (Ohm m) to the inversion model., a SimPEG Map
  • sigma (PhysicalProperty): Electrical conductivity (S/m), a physical property
  • sigmaMap (Mapping): Mapping of Electrical conductivity (S/m) to the inversion model., a SimPEG Map

Other Properties:

  • rhoDeriv (Derivative): Derivative of Electrical resistivity (Ohm m) wrt the model.
  • sigmaDeriv (Derivative): Derivative of Electrical conductivity (S/m) wrt the model.
surveyPair

alias of Survey

fieldsPair

alias of FieldsDC

Ainv = None
fields(m=None)[source]
Jvec(m, v, f=None)[source]
Jtvec(m, v, f=None)[source]
getSourceTerm()[source]

Evaluates the sources, and puts them in matrix form

Return type:tuple
Returns:q (nC or nN, nSrc)
class SimPEG.EM.Static.DC.ProblemDC.Problem3D_CC(mesh, **kwargs)[source]

Bases: SimPEG.EM.Static.DC.ProblemDC.BaseDCProblem

3D cell centered DC problem

Optional Properties:

  • model (Model): Inversion model., a numpy array of <type ‘float’>, <type ‘int’> with shape (*)
  • mu (PhysicalProperty): Magnetic Permeability (H/m), a physical property, Default: 1.25663706144e-06
  • mui (PhysicalProperty): Inverse Magnetic Permeability (m/H), a physical property
  • rho (PhysicalProperty): Electrical resistivity (Ohm m), a physical property
  • rhoMap (Mapping): Mapping of Electrical resistivity (Ohm m) to the inversion model., a SimPEG Map
  • sigma (PhysicalProperty): Electrical conductivity (S/m), a physical property
  • sigmaMap (Mapping): Mapping of Electrical conductivity (S/m) to the inversion model., a SimPEG Map

Other Properties:

  • rhoDeriv (Derivative): Derivative of Electrical resistivity (Ohm m) wrt the model.
  • sigmaDeriv (Derivative): Derivative of Electrical conductivity (S/m) wrt the model.
fieldsPair

alias of Fields_CC

bc_type = u'Neumann'
getA()[source]

Make the A matrix for the cell centered DC resistivity problem

A = D MfRhoI G

getADeriv(u, v, adjoint=False)[source]
getRHS()[source]

RHS for the DC problem

q

getRHSDeriv(src, v, adjoint=False)[source]

Derivative of the right hand side with respect to the model

setBC()[source]
class SimPEG.EM.Static.DC.ProblemDC.Problem3D_N(mesh, **kwargs)[source]

Bases: SimPEG.EM.Static.DC.ProblemDC.BaseDCProblem

3D nodal DC problem

Optional Properties:

  • model (Model): Inversion model., a numpy array of <type ‘float’>, <type ‘int’> with shape (*)
  • mu (PhysicalProperty): Magnetic Permeability (H/m), a physical property, Default: 1.25663706144e-06
  • mui (PhysicalProperty): Inverse Magnetic Permeability (m/H), a physical property
  • rho (PhysicalProperty): Electrical resistivity (Ohm m), a physical property
  • rhoMap (Mapping): Mapping of Electrical resistivity (Ohm m) to the inversion model., a SimPEG Map
  • sigma (PhysicalProperty): Electrical conductivity (S/m), a physical property
  • sigmaMap (Mapping): Mapping of Electrical conductivity (S/m) to the inversion model., a SimPEG Map

Other Properties:

  • rhoDeriv (Derivative): Derivative of Electrical resistivity (Ohm m) wrt the model.
  • sigmaDeriv (Derivative): Derivative of Electrical conductivity (S/m) wrt the model.
fieldsPair

alias of Fields_N

getA()[source]

Make the A matrix for the cell centered DC resistivity problem

A = G.T MeSigma G

getADeriv(u, v, adjoint=False)[source]

Product of the derivative of our system matrix with respect to the model and a vector

getRHS()[source]

RHS for the DC problem

q

getRHSDeriv(src, v, adjoint=False)[source]

Derivative of the right hand side with respect to the model

class SimPEG.EM.Static.DC.ProblemDC_2D.BaseDCProblem_2D(mesh, **kwargs)[source]

Bases: SimPEG.EM.Base.BaseEMProblem

Base 2.5D DC problem

Optional Properties:

  • model (Model): Inversion model., a numpy array of <type ‘float’>, <type ‘int’> with shape (*)
  • mu (PhysicalProperty): Magnetic Permeability (H/m), a physical property, Default: 1.25663706144e-06
  • mui (PhysicalProperty): Inverse Magnetic Permeability (m/H), a physical property
  • rho (PhysicalProperty): Electrical resistivity (Ohm m), a physical property
  • rhoMap (Mapping): Mapping of Electrical resistivity (Ohm m) to the inversion model., a SimPEG Map
  • sigma (PhysicalProperty): Electrical conductivity (S/m), a physical property
  • sigmaMap (Mapping): Mapping of Electrical conductivity (S/m) to the inversion model., a SimPEG Map

Other Properties:

  • rhoDeriv (Derivative): Derivative of Electrical resistivity (Ohm m) wrt the model.
  • sigmaDeriv (Derivative): Derivative of Electrical conductivity (S/m) wrt the model.
surveyPair

alias of Survey_ky

fieldsPair

alias of Fields_ky

nky = 15
kys = array([ 1.00000000e-04, 2.27584593e-04, 5.17947468e-04, 1.17876863e-03, 2.68269580e-03, 6.10540230e-03, 1.38949549e-02, 3.16227766e-02, 7.19685673e-02, 1.63789371e-01, 3.72759372e-01, 8.48342898e-01, 1.93069773e+00, 4.39397056e+00, 1.00000000e+01])
Ainv = [None, None, None, None, None, None, None, None, None, None, None, None, None, None, None]
nT = 15
fields(m)[source]
Jvec(m, v, f=None)[source]
Jtvec(m, v, f=None)[source]
getSourceTerm(ky)[source]

takes concept of source and turns it into a matrix

i = 14
class SimPEG.EM.Static.DC.ProblemDC_2D.Problem2D_CC(mesh, **kwargs)[source]

Bases: SimPEG.EM.Static.DC.ProblemDC_2D.BaseDCProblem_2D

2.5D cell centered DC problem

Optional Properties:

  • model (Model): Inversion model., a numpy array of <type ‘float’>, <type ‘int’> with shape (*)
  • mu (PhysicalProperty): Magnetic Permeability (H/m), a physical property, Default: 1.25663706144e-06
  • mui (PhysicalProperty): Inverse Magnetic Permeability (m/H), a physical property
  • rho (PhysicalProperty): Electrical resistivity (Ohm m), a physical property
  • rhoMap (Mapping): Mapping of Electrical resistivity (Ohm m) to the inversion model., a SimPEG Map
  • sigma (PhysicalProperty): Electrical conductivity (S/m), a physical property
  • sigmaMap (Mapping): Mapping of Electrical conductivity (S/m) to the inversion model., a SimPEG Map

Other Properties:

  • rhoDeriv (Derivative): Derivative of Electrical resistivity (Ohm m) wrt the model.
  • sigmaDeriv (Derivative): Derivative of Electrical conductivity (S/m) wrt the model.
fieldsPair

alias of Fields_ky_CC

getA(ky)[source]

Make the A matrix for the cell centered DC resistivity problem

A = D MfRhoI G

getADeriv(ky, u, v, adjoint=False)[source]
getRHS(ky)[source]

RHS for the DC problem

q

getRHSDeriv(ky, src, v, adjoint=False)[source]

Derivative of the right hand side with respect to the model

setBC()[source]
class SimPEG.EM.Static.DC.ProblemDC_2D.Problem2D_N(mesh, **kwargs)[source]

Bases: SimPEG.EM.Static.DC.ProblemDC_2D.BaseDCProblem_2D

2.5D nodal DC problem

Optional Properties:

  • model (Model): Inversion model., a numpy array of <type ‘float’>, <type ‘int’> with shape (*)
  • mu (PhysicalProperty): Magnetic Permeability (H/m), a physical property, Default: 1.25663706144e-06
  • mui (PhysicalProperty): Inverse Magnetic Permeability (m/H), a physical property
  • rho (PhysicalProperty): Electrical resistivity (Ohm m), a physical property
  • rhoMap (Mapping): Mapping of Electrical resistivity (Ohm m) to the inversion model., a SimPEG Map
  • sigma (PhysicalProperty): Electrical conductivity (S/m), a physical property
  • sigmaMap (Mapping): Mapping of Electrical conductivity (S/m) to the inversion model., a SimPEG Map

Other Properties:

  • rhoDeriv (Derivative): Derivative of Electrical resistivity (Ohm m) wrt the model.
  • sigmaDeriv (Derivative): Derivative of Electrical conductivity (S/m) wrt the model.
fieldsPair

alias of Fields_ky_N

MnSigma

Node inner product matrix for (sigma). Used in the E-B formulation

MnSigmaDeriv(u)[source]

Derivative of MnSigma with respect to the model

getA(ky)[source]

Make the A matrix for the cell centered DC resistivity problem

A = D MfRhoI G

getADeriv(ky, u, v, adjoint=False)[source]
getRHS(ky)[source]

RHS for the DC problem

q

getRHSDeriv(ky, src, v, adjoint=False)[source]

Derivative of the right hand side with respect to the model

Survey

class SimPEG.EM.Static.DC.SurveyDC.Survey(srcList, **kwargs)[source]

Bases: SimPEG.EM.Base.BaseEMSurvey

Base DC survey

rxPair

alias of BaseRx

srcPair

alias of BaseSrc

class SimPEG.EM.Static.DC.SurveyDC.Survey_ky(srcList, **kwargs)[source]

Bases: SimPEG.EM.Base.BaseEMSurvey

2.5D survey

rxPair

alias of BaseRx

srcPair

alias of BaseSrc

eval(f)[source]

Project fields to receiver locations :param Fields u: fields object :rtype: numpy.ndarray :return: data

class SimPEG.EM.Static.DC.SrcDC.BaseSrc(rxList, **kwargs)[source]

Bases: SimPEG.Survey.BaseSrc

Base DC source

current = 1.0
loc = None
eval(prob)[source]
evalDeriv(prob)[source]
class SimPEG.EM.Static.DC.SrcDC.Dipole(rxList, locA, locB, **kwargs)[source]

Bases: SimPEG.EM.Static.DC.SrcDC.BaseSrc

Dipole source

eval(prob)[source]
class SimPEG.EM.Static.DC.SrcDC.Pole(rxList, loc, **kwargs)[source]

Bases: SimPEG.EM.Static.DC.SrcDC.BaseSrc

eval(prob)[source]
class SimPEG.EM.Static.DC.RxDC.BaseRx(locs, rxType, **kwargs)[source]

Bases: SimPEG.Survey.BaseRx

Base DC receiver

locs = None
rxType = None
knownRxTypes = {u'phi': [u'phi', None], u'jx': [u'j', u'x'], u'jy': [u'j', u'y'], u'jz': [u'j', u'z'], u'ey': [u'e', u'y'], u'ex': [u'e', u'x'], u'ez': [u'e', u'z']}
projField

Field Type projection (e.g. e b ...)

projGLoc(f)[source]

Grid Location projection (e.g. Ex Fy ...)

eval(src, mesh, f)[source]
evalDeriv(src, mesh, f, v, adjoint=False)[source]
class SimPEG.EM.Static.DC.RxDC.Dipole(locsM, locsN, rxType=u'phi', **kwargs)[source]

Bases: SimPEG.EM.Static.DC.RxDC.BaseRx

Dipole receiver

nD

Number of data in the receiver.

getP(mesh, Gloc)[source]
class SimPEG.EM.Static.DC.RxDC.Dipole_ky(locsM, locsN, rxType=u'phi', **kwargs)[source]

Bases: SimPEG.EM.Static.DC.RxDC.BaseRx

Dipole receiver for 2.5D simulations

nD

Number of data in the receiver.

getP(mesh, Gloc)[source]
eval(kys, src, mesh, f)[source]
evalDeriv(ky, src, mesh, f, v, adjoint=False)[source]
IntTrapezoidal(kys, Pf, y=0.0)[source]
class SimPEG.EM.Static.DC.RxDC.Pole(locsM, rxType=u'phi', **kwargs)[source]

Bases: SimPEG.EM.Static.DC.RxDC.BaseRx

Pole receiver

nD

Number of data in the receiver.

getP(mesh, Gloc)[source]
class SimPEG.EM.Static.DC.RxDC.Pole_ky(locsM, rxType=u'phi', **kwargs)[source]

Bases: SimPEG.EM.Static.DC.RxDC.BaseRx

Pole receiver for 2.5D simulations

nD

Number of data in the receiver.

getP(mesh, Gloc)[source]
eval(kys, src, mesh, f)[source]
evalDeriv(ky, src, mesh, f, v, adjoint=False)[source]
IntTrapezoidal(kys, Pf, y=0.0)[source]

Fields

class SimPEG.EM.Static.DC.FieldsDC.FieldsDC(mesh, survey, **kwargs)[source]

Bases: SimPEG.Fields.Fields

knownFields = {}
dtype

alias of float

class SimPEG.EM.Static.DC.FieldsDC.Fields_CC(mesh, survey, **kwargs)[source]

Bases: SimPEG.EM.Static.DC.FieldsDC.FieldsDC

knownFields = {u'phiSolution': u'CC'}
aliasFields = {u'phi': [u'phiSolution', u'CC', u'_phi'], u'j': [u'phiSolution', u'F', u'_j'], u'e': [u'phiSolution', u'F', u'_e'], u'charge': [u'phiSolution', u'CC', u'_charge']}
startup()[source]
class SimPEG.EM.Static.DC.FieldsDC.Fields_N(mesh, survey, **kwargs)[source]

Bases: SimPEG.EM.Static.DC.FieldsDC.FieldsDC

knownFields = {u'phiSolution': u'N'}
aliasFields = {u'phi': [u'phiSolution', u'N', u'_phi'], u'j': [u'phiSolution', u'E', u'_j'], u'e': [u'phiSolution', u'E', u'_e'], u'charge': [u'phiSolution', u'N', u'_charge']}
startup()[source]
class SimPEG.EM.Static.DC.FieldsDC_2D.Fields_ky(mesh, survey, **kwargs)[source]

Bases: SimPEG.Fields.TimeFields

Fancy Field Storage for a 2.5D code.

u[:,’phi’, kyInd] = phi print(u[src0,’phi’])

Only one field type is stored for each problem, the rest are computed. The fields obejct acts like an array and is indexed by

f = problem.fields(m)
e = f[srcList,'e']
j = f[srcList,'j']

If accessing all sources for a given field, use the

f = problem.fields(m)
phi = f[:,'phi']
e = f[:,'e']
b = f[:,'b']

The array returned will be size (nE or nF, nSrcs \(\times\) nFrequencies)

knownFields = {}
dtype

alias of float

class SimPEG.EM.Static.DC.FieldsDC_2D.Fields_ky_CC(mesh, survey, **kwargs)[source]

Bases: SimPEG.EM.Static.DC.FieldsDC_2D.Fields_ky

Fancy Field Storage for a 2.5D cell centered code.

knownFields = {u'phiSolution': u'CC'}
aliasFields = {u'phi': [u'phiSolution', u'CC', u'_phi'], u'j': [u'phiSolution', u'F', u'_j'], u'e': [u'phiSolution', u'F', u'_e']}
startup()[source]
class SimPEG.EM.Static.DC.FieldsDC_2D.Fields_ky_N(mesh, survey, **kwargs)[source]

Bases: SimPEG.EM.Static.DC.FieldsDC_2D.Fields_ky

Fancy Field Storage for a 2.5D nodal code.

knownFields = {u'phiSolution': u'N'}
aliasFields = {u'phi': [u'phiSolution', u'N', u'_phi'], u'j': [u'phiSolution', u'E', u'_j'], u'e': [u'phiSolution', u'E', u'_e']}
startup()[source]

Utils

SimPEG.EM.Static.DC.Utils.WennerSrcList(nElecs, aSpacing, in2D=False, plotIt=False)[source]

Source list for a Wenner Array

SimPEG.EM.Static.DC.BoundaryUtils.getxBCyBC_CC(mesh, alpha, beta, gamma)[source]

This is a subfunction generating mixed-boundary condition:

\[ \begin{align}\begin{aligned}\nabla \cdot \vec{j} = -\nabla \cdot \vec{j}_s = q\\\rho \vec{j} = -\nabla \phi\\\alpha \phi + \beta \frac{\partial \phi}{\partial r} = \gamma \quad \vec{r} \subset \partial \Omega\\xBC = f_1(\alpha, \beta, \gamma)\\yBC = f(\alpha, \beta, \gamma)\end{aligned}\end{align} \]

Computes xBC and yBC for cell-centered discretizations