\[\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}}\]


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.


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.


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)


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!


  • 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


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.

alias of BaseMagSurvey


alias of BaseMagMap

getB0(*args, **kwargs)[source]

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

\[\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 creates and returns the A matrix for the Magnetics problem

The A matrix has the form:

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

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


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} \]


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


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.

Return type:




Jvec_approx(m, v, f=None)

Approximate effect of J(m) on a vector v

Return type:




class Solver(A, **kwargs)

Bases: object

Problem3D_DiffSecondary.counter = None

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


True if the problem is paired to a survey.


alias of IdentityMap


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

Problem3D_DiffSecondary.mesh = None

Inversion model.


Magnetic Permeability (H/m)


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


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


Inverse Magnetic Permeability (m/H)


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


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


True if a model is necessary


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}

The survey object for this problem.


Unbind a survey from this problem instance.


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]

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

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.


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:

dtrue = None
eps = None

This function projects the fields onto the data space.

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

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

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

Check if the data is synthetic.

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

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

  • 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 of the paired problem.

mtrue = None

Number of data


Number of Sources


Bind a problem to this survey instance using pointers


The geophysical problem that explains this survey, use:

residual(m, f=None)
Return type:



data residual

The data residual:

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

Source List


alias of BaseSrc

std = None

Unbind a problem from this survey instance


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|