\[\renewcommand{\div}{\nabla\cdot\,} \newcommand{\grad}{\vec \nabla} \newcommand{\curl}{{\vec \nabla}\times\,} \newcommand {\J}{{\vec J}} \renewcommand{\H}{{\vec H}} \newcommand {\E}{{\vec E}} \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}}} \newcommand{\St}{{\mathbf \Sigma_\tau}} \newcommand{\T}{{\mathbf T}} \newcommand{\Tt}{{\mathbf T_\tau}} \newcommand{\diag}{\mathbf{diag}} \newcommand{\M}{{\mathbf M}} \newcommand{\MfMui}{{\M^f_{\mu^{-1}}}} \newcommand{\dMfMuI}{{d_m (\M^f_{\mu^{-1}})^{-1}}} \newcommand{\MeSig}{{\M^e_\sigma}} \newcommand{\MeSigInf}{{\M^e_{\sigma_\infty}}} \newcommand{\MeSigO}{{\M^e_{\sigma_0}}} \newcommand{\Me}{{\M^e}} \newcommand{\Mes}[1]{{\M^e_{#1}}} \newcommand{\Mee}{{\M^e_e}} \newcommand{\Mej}{{\M^e_j}} \newcommand{\BigO}[1]{\mathcal{O}\bigl(#1\bigr)} \newcommand{\bE}{\mathbf{E}} \newcommand{\bH}{\mathbf{H}} \newcommand{\B}{\vec{B}} \newcommand{\D}{\vec{D}} \renewcommand{\H}{\vec{H}} \newcommand{\s}{\vec{s}} \newcommand{\bfJ}{\bf{J}} \newcommand{\vecm}{\vec m} \renewcommand{\Re}{\mathsf{Re}} \renewcommand{\Im}{\mathsf{Im}} \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}}\]

Magnetics

The geomagnetic field can be ranked as the longest studied of all the geophysical properties of the earth. In addition, magnetic survey, has been used broadly in diverse realm e.g., mining, oil and gas industry and environmental engineering. Although, this geophysical application is quite common in geoscience; however, we do not have modular, well-documented and well-tested open-source codes, which perform forward and inverse problems of magnetic survey. Therefore, here we are going to build up magnetic forward and inverse modeling code based on two common methodologies for forward problem - differential equation and integral equation approaches.

First, we start with some backgrounds of magnetics, e.g., Maxwell’s equations. Based on that secondly, we use differential equation approach to solve forward problem with secondary field formulation. In order to discretzie our system here, we use finite volume approach with weak formulation. Third, we solve inverse problem through Gauss-Newton method.

Backgrounds

Maxwell’s equations for static case with out current source can be written as

\[ \begin{align}\begin{aligned}\begin{split}\nabla U = \frac{1}{\mu}\vec{B} \\\end{split}\\\nabla \cdot \vec{B} = 0\end{aligned}\end{align} \]

where \(\vec{B}\) is magnetic flux (\(T\)) and \(U\) is magnetic potential and \(\mu\) is permeability. Since we do not have any source term in above equations, boundary condition is going to be the driving force of our system as given below

\[(\vec{B}\cdot{\vec{n}})_{\partial\Omega} = B_{BC}\]

where \(\vec{n}\) means the unit normal vector on the boundary surface (\(\partial \Omega\)). By using seocondary field formulation we can rewrite above equations as

\[ \begin{align}\begin{aligned}\frac{1}{\mu}\vec{B}_s = (\frac{1}{\mu}_0-\frac{1}{\mu})\vec{B}_0+\nabla\phi_s\\\nabla \cdot \vec{B}_s = 0\\(\vec{B}_s\cdot{\vec{n}})_{\partial\Omega} = B_{sBC}\end{aligned}\end{align} \]

where \(\vec{B}_s\) is the secondary magnetic flux and \(\vec{B}_0\) is the background or primary magnetic flux. In practice, we consider our earth field, which we can get from International Geomagnetic Reference Field (IGRF) by specifying the time and location, as \(\vec{B}_0\). And based on this background fields, we compute secondary fields (\(\vec{B}_s\)). Now we introduce the susceptibility as

\[ \begin{align}\begin{aligned}\begin{split}\chi = \frac{\mu}{\mu_0} - 1 \\\end{split}\\\mu = \mu_0(1+\chi)\end{aligned}\end{align} \]

Since most materials in the earth have lower permeability than \(\mu_0\), usually \(\chi\) is greater than 0.

Note

Actually, this is an assumption, which means we are not sure exactly this is true, although we are sure, it is very rare that we can encounter those materials. Anyway, practical range of the susceptibility is \(0 < \chi < 1 \).

Since we compute secondary field based on the earth field, which can be different from different locations in the world, we can expect different anomalous responses in different locations in the earth. For instance, assume we have two susceptible spheres, which are exactly same. However, anomalous responses in Seoul and Vancouver are going to be different.

(Source code, png, hires.png, pdf)

../../_images/index-12.png

Since we can measure total fields ( \(\vec{B}\)), and usually have reasonably accurate earth field (\(\vec{B}_0\)), we can compute anomalous fields, \(\vec{B}_s\) from our observed data. If you want to download earth magnetic fields at specific location see this website (noaa).

What is our data?

In applied geophysics, which means in practice, it is common to refer to measurements as “the magnetic anomaly” and we can consider this as our observed data. For further descriptions in GPG materials for magnetic survey. Now we have the simple relation ship between “the magnetic anomaly” and the total field as

\[\triangle\vec{B} = |\hat{B}_o-\vec{B}_s|-|\hat{B}_o| \approx |\vec{B}_s|cos \theta\]

where \(\theta\) is the angle between total and anomalous fields, \(\hat{B}_o\) is the unit vector for \(\vec{B}_o\). Equivalently, we can use the vector dot product to show that the anomalous field is approximately equal to the projection of that field onto the direction of the inducing field. Using this approach we would write

\[\triangle\vec{B} = |\vec{B}_s|cos \theta = |\hat{B}_o||\vec{B}_s|cos \theta = \hat{B}_o \cdot \vec{B}_s\]

This is important because, in practice we usually use a total field magnetometer (like a proton precession or optically pumped sensor), which can measure only that part of the anomalous field which is in the direction of the earth’s main field.

Sphere in a whole space

Forward problem

Differential equation approach

\[ \begin{align}\begin{aligned}\mathbf{A}\mathbf{u} = \mathbf{rhs}\\\mathbf{A} = \Div(\MfMui)^{-1}\Div^{T}\\\mathbf{rhs} = \Div(\MfMui)^{-1}\mathbf{M}^f_{\mu_0^{-1}}\mathbf{B}_0 - \Div\mathbf{B}_0+\diag(v)\mathbf{D} \mathbf{P}_{out}^T \mathbf{B}_{sBC}\\\mathbf{B}_s = (\MfMui)^{-1}\mathbf{M}^f_{\mu_0^{-1}}\mathbf{B}_0-\mathbf{B}_0 -(\MfMui)^{-1}\Div^T \mathbf{u}\end{aligned}\end{align} \]

Mag Differential eq. approach

class SimPEG.PF.Magnetics.Problem3D_DiffSecondary(mesh, **kwargs)[source]

Bases: SimPEG.Problem.BaseProblem

Secondary field approach using differential equations!

Optional

Parameters:
  • 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
  • muMap (Mapping) – Mapping of Magnetic Permeability (H/m) to the inversion model., a SimPEG Map
  • mui (PhysicalProperty) – Inverse Magnetic Permeability (m/H), a physical property
  • muiMap (Mapping) – Mapping of Inverse Magnetic Permeability (m/H) to the inversion model., a SimPEG Map

Immutable

Attribute muDeriv:
 (Derivative) - Derivative of Magnetic Permeability (H/m) wrt the model.
Attribute muiDeriv:
 (Derivative) - Derivative of Inverse Magnetic Permeability (m/H) wrt the model.
surveyPair

alias of BaseMagSurvey

modelPair

alias of BaseMagMap

MfMuI
MfMui
MfMu0
makeMassMatrices(m)[source]
getB0(*args, **kwargs)[source]

To use getB0 method, SimPEG requires that the survey be specified.

getRHS(m)[source]
\[\mathbf{rhs} = \Div(\MfMui)^{-1}\mathbf{M}^f_{\mu_0^{-1}}\mathbf{B}_0 - \Div\mathbf{B}_0+\diag(v)\mathbf{D} \mathbf{P}_{out}^T \mathbf{B}_{sBC}\]
getA(m)[source]

GetA creates and returns the A matrix for the Magnetics problem

The A matrix has the form:

\[\mathbf{A} = \Div(\MfMui)^{-1}\Div^{T}\]
fields(m)[source]

Return magnetic potential (u) and flux (B) u: defined on the cell center [nC x 1] B: defined on the cell center [nF x 1]

After we compute u, then we update B.

\[\mathbf{B}_s = (\MfMui)^{-1}\mathbf{M}^f_{\mu_0^{-1}}\mathbf{B}_0-\mathbf{B}_0 -(\MfMui)^{-1}\Div^T \mathbf{u}\]
Jvec(*args, **kwargs)[source]

Computing Jacobian multiplied by vector

By setting our problem as

\[\mathbf{C}(\mathbf{m}, \mathbf{u}) = \mathbf{A}\mathbf{u} - \mathbf{rhs} = 0\]

And taking derivative w.r.t m

\[ \begin{align}\begin{aligned}\nabla \mathbf{C}(\mathbf{m}, \mathbf{u}) = \nabla_m \mathbf{C}(\mathbf{m}) \delta \mathbf{m} + \nabla_u \mathbf{C}(\mathbf{u}) \delta \mathbf{u} = 0\\\frac{\delta \mathbf{u}}{\delta \mathbf{m}} = - [\nabla_u \mathbf{C}(\mathbf{u})]^{-1}\nabla_m \mathbf{C}(\mathbf{m})\end{aligned}\end{align} \]

With some linear algebra we can have

\[ \begin{align}\begin{aligned}\nabla_u \mathbf{C}(\mathbf{u}) = \mathbf{A}\\\nabla_m \mathbf{C}(\mathbf{m}) = \frac{\partial \mathbf{A}}{\partial \mathbf{m}}(\mathbf{m})\mathbf{u} - \frac{\partial \mathbf{rhs}(\mathbf{m})}{\partial \mathbf{m}}\end{aligned}\end{align} \]
\[ \begin{align}\begin{aligned}\frac{\partial \mathbf{A}}{\partial \mathbf{m}}(\mathbf{m})\mathbf{u} = \frac{\partial \mathbf{\mu}}{\partial \mathbf{m}} \left[\Div \diag (\Div^T \mathbf{u}) \dMfMuI \right]\\\dMfMuI = \diag(\MfMui)^{-1}_{vec} \mathbf{Av}_{F2CC}^T\diag(\mathbf{v})\diag(\frac{1}{\mu^2})\\\frac{\partial \mathbf{rhs}(\mathbf{m})}{\partial \mathbf{m}} = \frac{\partial \mathbf{\mu}}{\partial \mathbf{m}} \left[ \Div \diag(\M^f_{\mu_{0}^{-1}}\mathbf{B}_0) \dMfMuI \right] - \diag(\mathbf{v})\mathbf{D} \mathbf{P}_{out}^T\frac{\partial B_{sBC}}{\partial \mathbf{m}}\end{aligned}\end{align} \]

In the end,

\[\frac{\delta \mathbf{u}}{\delta \mathbf{m}} = - [ \mathbf{A} ]^{-1}\left[ \frac{\partial \mathbf{A}}{\partial \mathbf{m}}(\mathbf{m})\mathbf{u} - \frac{\partial \mathbf{rhs}(\mathbf{m})}{\partial \mathbf{m}} \right]\]

A little tricky point here is we are not interested in potential (u), but interested in magnetic flux (B). Thus, we need sensitivity for B. Now we take derivative of B w.r.t m and have

\[ \begin{align}\begin{aligned}\frac{\delta \mathbf{B}} {\delta \mathbf{m}} = \frac{\partial \mathbf{\mu} } {\partial \mathbf{m} } \left[ \diag(\M^f_{\mu_{0}^{-1} } \mathbf{B}_0) \dMfMuI \ - \diag (\Div^T\mathbf{u})\dMfMuI \right ]\\ - (\MfMui)^{-1}\Div^T\frac{\delta\mathbf{u}}{\delta \mathbf{m}}\end{aligned}\end{align} \]

Finally we evaluate the above, but we should remember that

Note

We only want to evalute

\[\mathbf{J}\mathbf{v} = \frac{\delta \mathbf{P}\mathbf{B}} {\delta \mathbf{m}}\mathbf{v}\]

Since forming sensitivity matrix is very expensive in that this monster is “big” and “dense” matrix!!

Jtvec(*args, **kwargs)[source]
Computing Jacobian^T multiplied by vector.
\[ \begin{align}\begin{aligned}(\frac{\delta \mathbf{P}\mathbf{B}} {\delta \mathbf{m}})^{T} = \left[ \mathbf{P}_{deriv}\frac{\partial \mathbf{\mu} } {\partial \mathbf{m} } \left[ \diag(\M^f_{\mu_{0}^{-1} } \mathbf{B}_0) \dMfMuI \ - \diag (\Div^T\mathbf{u})\dMfMuI \right ]\right]^{T}\\ - \left[\mathbf{P}_{deriv}(\MfMui)^{-1}\Div^T\frac{\delta\mathbf{u}}{\delta \mathbf{m}} \right]^{T}\end{aligned}\end{align} \]

where

\[\mathbf{P}_{derv} = \frac{\partial \mathbf{P}}{\partial\mathbf{B}}\]

Note

Here we only want to compute

\[\mathbf{J}^{T}\mathbf{v} = (\frac{\delta \mathbf{P}\mathbf{B}} {\delta \mathbf{m}})^{T} \mathbf{v}\]
Jtvec_approx(m, v, f=None)

Approximate effect of transpose of J(m) on a vector v.

Parameters:
Return type:

numpy.array

Returns:

JTv

Jvec_approx(m, v, f=None)

Approximate effect of J(m) on a vector v

Parameters:
Return type:

numpy.array

Returns:

approxJv

class Solver(A, **kwargs)

Bases: object

clean()
Problem3D_DiffSecondary.counter = None
Problem3D_DiffSecondary.curModel

Setting the curModel is depreciated.

Use SimPEG.Problem.model instead.

Problem3D_DiffSecondary.deleteTheseOnModelUpdate = []
Problem3D_DiffSecondary.deserialize(value, trusted=False, verbose=True, **kwargs)

Creates new HasProperties instance from JSON dictionary

Problem3D_DiffSecondary.ispaired

True if the problem is paired to a survey.

Problem3D_DiffSecondary.mapPair

alias of IdentityMap

Problem3D_DiffSecondary.mapping

Setting an unnamed mapping has been depreciated in v0.4.0. Please see the release notes for more details.

Problem3D_DiffSecondary.mesh = None
Problem3D_DiffSecondary.model

Inversion model.

Problem3D_DiffSecondary.mu

Magnetic Permeability (H/m)

Problem3D_DiffSecondary.muDeriv

Derivative of Magnetic Permeability (H/m) wrt the model.

Problem3D_DiffSecondary.muMap

Mapping of Magnetic Permeability (H/m) to the inversion model.

Problem3D_DiffSecondary.mui

Inverse Magnetic Permeability (m/H)

Problem3D_DiffSecondary.muiDeriv

Derivative of Inverse Magnetic Permeability (m/H) wrt the model.

Problem3D_DiffSecondary.muiMap

Mapping of Inverse Magnetic Permeability (m/H) to the inversion model.

Problem3D_DiffSecondary.needs_model

True if a model is necessary

Problem3D_DiffSecondary.pair(d)

Bind a survey to this problem instance using pointers.

Problem3D_DiffSecondary.serialize(include_class=True, **kwargs)

Serializes a HasProperties instance to JSON

Problem3D_DiffSecondary.solverOpts = {'accuracyTol': 0.0001}
Problem3D_DiffSecondary.summary()
Problem3D_DiffSecondary.survey

The survey object for this problem.

Problem3D_DiffSecondary.unpair()

Unbind a survey from this problem instance.

Problem3D_DiffSecondary.validate()

Call all the registered ClassValidators

class SimPEG.PF.BaseMag.BaseMagSurvey(**kwargs)[source]

Bases: SimPEG.Survey.BaseSurvey

Base Magnetics Survey

rxLoc = None

receiver locations

rxType = None

receiver type

setBackgroundField(Inc, Dec, Btot)[source]
Qfx
Qfy
Qfz
projectFields(u)[source]

This function projects the fields onto the data space.

Especially, here for we use total magnetic intensity (TMI) data, which is common in practice.

First we project our B on to data location

\[\mathbf{B}_{rec} = \mathbf{P} \mathbf{B}\]

then we take the dot product between B and b_0

\[\text{TMI} = \vec{B}_s \cdot \hat{B}_0\]
projectFieldsDeriv(*args, **kwargs)[source]

This function projects the fields onto the data space.

\[\frac{\partial d_\text{pred}}{\partial \mathbf{B}} = \mathbf{P}\]

Especially, this function is for TMI data type

projectFieldsAsVector(B)[source]
counter = None
dobs = None
dpred(m, f=None)

Create the projected data from a model. The fields, f, (if provided) will be used for the predicted data instead of recalculating the fields (which may be expensive!).

\[d_\text{pred} = P(f(m))\]

Where P is a projection of the fields onto the data space.

Note

To use survey.dpred(), SimPEG requires that a problem be bound to the survey. If a problem has not been bound, an Exception will be raised. To bind a problem to the Data object:

survey.pair(myProblem)
dtrue = None
eps = None
eval(f)

This function projects the fields onto the data space.

\[d_\text{pred} = \mathbf{P} f(m)\]
evalDeriv(f)

This function s the derivative of projects the fields onto the data space.

\[\frac{\partial d_\text{pred}}{\partial u} = \mathbf{P}\]
getSourceIndex(sources)
isSynthetic

Check if the data is synthetic.

ispaired
makeSyntheticData(m, std=0.05, f=None, force=False)

Make synthetic data given a model, and a standard deviation.

Parameters:
  • m (numpy.array) – geophysical model
  • std (numpy.array) – standard deviation
  • u (numpy.array) – fields for the given model (if pre-calculated)
  • force (bool) – force overwriting of dobs
mesh

Mesh of the paired problem.

mtrue = None
nD

Number of data

nSrc

Number of Sources

pair(p)

Bind a problem to this survey instance using pointers

prob

The geophysical problem that explains this survey, use:

survey.pair(prob)
residual(m, f=None)
Parameters:
Return type:

numpy.array

Returns:

data residual

The data residual:

\[\mu_\text{data} = \mathbf{d}_\text{pred} - \mathbf{d}_\text{obs}\]
srcList

Source List

srcPair

alias of BaseSrc

std = None
unpair()

Unbind a problem from this survey instance

vnD

Vector number of data

Mag Integral eq. approach

Mag analytic solutions

SimPEG.PF.MagAnalytics.spheremodel(mesh, x0, y0, z0, r)[source]

Generate model indicies for sphere - (x0, y0, z0 ): is the center location of sphere - r: is the radius of the sphere - it returns logical indicies of cell-center model

SimPEG.PF.MagAnalytics.MagSphereAnaFun(x, y, z, R, x0, y0, z0, mu1, mu2, H0, flag='total')[source]

test Analytic function for Magnetics problem. The set up here is magnetic sphere in whole-space assuming that the inducing field is oriented in the x-direction.

  • (x0, y0, z0)
  • (x0, y0, z0 ): is the center location of sphere
  • r: is the radius of the sphere
\[\mathbf{H}_0 = H_0\hat{x}\]
SimPEG.PF.MagAnalytics.CongruousMagBC(mesh, Bo, chi)[source]

Computing boundary condition using Congrous sphere method. This is designed for secondary field formulation.

>> Input

  • mesh: Mesh class
  • Bo: np.array([Box, Boy, Boz]): Primary magnetic flux
  • chi: susceptibility at cell volume
\[\vec{B}(r) = \frac{\mu_0}{4\pi} \frac{m}{ \| \vec{r} - \vec{r}_0\|^3}[3\hat{m}\cdot\hat{r}-\hat{m}]\]
SimPEG.PF.MagAnalytics.MagSphereAnaFunA(x, y, z, R, xc, yc, zc, chi, Bo, flag)[source]

Computing boundary condition using Congrous sphere method. This is designed for secondary field formulation. >> Input mesh: Mesh class Bo: np.array([Box, Boy, Boz]): Primary magnetic flux Chi: susceptibility at cell volume

\[\vec{B}(r) = \frac{\mu_0}{4\pi}\frac{m}{\| \vec{r}-\vec{r}_0\|^3}[3\hat{m}\cdot\hat{r}-\hat{m}]\]
SimPEG.PF.MagAnalytics.IDTtoxyz(Inc, Dec, Btot)[source]

Convert from Inclination, Declination, Total intensity of earth field to x, y, z

SimPEG.PF.MagAnalytics.MagSphereFreeSpace(x, y, z, R, xc, yc, zc, chi, Bo)[source]

Computing the induced response of magnetic sphere in free-space.

>> Input x, y, z: Observation locations R: radius of the sphere xc, yc, zc: Location of the sphere chi: Susceptibility of sphere Bo: Inducing field components [bx, by, bz]*|H0|