In [1]:
from MT1D import MT1DProblem, MT1DSurvey, MT1DSrc, ZxyRx, Survey, AppResPhaRx
from SimPEG import Maps, DataMisfit
from scipy.constants import mu_0
import numpy as np

Put all things together to SimPEG Problem

In the 1_MT1D_NumericalSetup notebook, we coded up a 1D MT simulation function that computes complex impedances at multiple frequencies for a given conductivity model. To invert MT data in 1D using gradient-based optimization, sensitivity functions were derived, coded up, and tested in Appendix_A_MT1D_Sensitivity notebook. Basically those are three funcionalities:

  • dpred(m): takes model and predict MT data for a given model, m

  • Jvec(m,v): takes model and a vector v, then computes sensitivity-vector product.

  • Jtvec(m,v): takes model and a vector v, then computes adjoint sensitivity-vector product.

Now we put all these things to SimPEG's framework.

SimPEG framework

Our setup of the simluation follows the SimPEG framework.

plane_wave

For more details, see the SimPEG docs

  • Mesh: mesh geometry and differential operators
  • Problem: physics engine. contains the machinery to construct and solve the PDE
  • Survey: sources and receivers
  • Src: sources
  • Rx: receivers

Problem class in SimPEG is a physics. It discretizes the earth and computes fields in a discretized domain, and here we use Mesh class. For an MT problem, Maxwell's equations are used, and electromagnetic fields are computed. To do this, we need to know details about MT survey (e.g. frequencies and location of source and receivers), and they are defined in Survey class. Survey class takes Src class, and Src class takes Rx class. Therefore, both information in sources and receivers can be aware in the Survey.

Mapping

For simulation, we do need to input discretized physical property. For instance in 1D MT problem values of conductivity at each cell is required. However, when consdiering inversion, inversion model can be different from physical property (e.g. conductivity). For instance, often in electromagnetic (EM) inversions logarithmic conductivity was used as an inversion model: $\mathbf{m} = log(\boldsymbol{\sigma})$. While inversion process, still simulation is required to predict data hence updated model, $\mathbf{m}$ has to be converted to $\boldsymbol{\sigma}$, which is a discretized physical property. Here we define a mapping:

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

mapping is a function which takes an inversion model, and outputs physical property used in simulation. For more details see: http://docs.simpeg.xyz/content/api_core/api_Maps.html.

Codes:

Taking snippets of codes written in 1_MT1D_NumericalSetup and Appendix_A_MT1D_Sensitivity, we generated the Problem, Survey, Src, and Rx classes for 1D MT problem. And here are lists of classes:

  • MT1DProblem: This includes functions: dpred(m), Jvec(m,v), Jtvec(m,v)
  • MT1DSurvey: Need to input MT1DSrc
  • MT1DSrc: Need to input ZxyRx
  • ZxyRx: Receiver class / Real and Imaginary part of $Z_{xy}$
In [2]:
rxloc = np.r_[0.]
srcloc = np.r_[0.]
frequency = np.logspace(-3, 2, 25)

# Rx
rx = ZxyRx(rxloc, component="both", frequency=frequency)
rxList = [rx]

# Src
src = MT1DSrc(rxList, loc=srcloc)

# Survey
survey = MT1DSurvey([src])
mesh = survey.setMesh(
    sigma=0.01, max_depth_core=5000., 
    ncell_per_skind=10, n_skind=2, 
    core_meshType="log", max_hz_core=1000.
)

# Physical property: Conductivity
sigma = np.ones(mesh.nC) * 0.01

# Problem
prob = MT1DProblem(mesh, sigmaMap=Maps.ExpMap(mesh), verbose=False)
prob.pair(survey)
>> Smallest cell size = 50 m
>> Padding distance = 316227 m
>> # of padding cells 17
>> # of core cells cells 16

Order and Adjoint tests

We perform both order and adjoint tests shown in Appendix_A_MT1D_Sensitivity using 1D MT problem in SimPEG's framework.

Order test: Jvec

In [3]:
from SimPEG import Tests
def derChk(m):
    return [survey.dpred(m), lambda mx: prob.Jvec(m, mx)]
Tests.checkDerivative(derChk, np.log(sigma), plotIt=False, num=3, eps=1e-20, dx=np.log(sigma)*2)
==================== checkDerivative ====================
iter    h         |ft-f0|   |ft-f0-h*J0*dx|  Order
---------------------------------------------------------
 0   1.00e-01    2.669e-01     5.676e-02      nan
 1   1.00e-02    2.151e-02     4.914e-04      2.063
 2   1.00e-03    2.106e-03     4.846e-06      2.006
========================= PASS! =========================
You are awesome.

Out[3]:
True

Order test: Jtvec

In [4]:
# Evaluate predicted data
pred = survey.dpred(np.log(sigma))
survey.dobs = pred
dmis = DataMisfit.l2_DataMisfit(survey)
m0 = np.log(sigma)*2.
SimPEG.DataMisfit.l2_DataMisfit assigning default eps of 1e-5 * ||dobs||
In [5]:
Tests.checkDerivative(
    lambda m: [dmis(m), dmis.deriv(m)],
    m0,
    plotIt=False,
    num=4,
    dx = m0*1.
)
==================== checkDerivative ====================
iter    h         |ft-f0|   |ft-f0-h*J0*dx|  Order
---------------------------------------------------------
 0   1.00e-01    1.074e+06     4.010e+05      nan
 1   1.00e-02    7.038e+04     3.059e+03      2.118
 2   1.00e-03    6.762e+03     2.982e+01      2.011
 3   1.00e-04    6.735e+02     2.975e-01      2.001
========================= PASS! =========================
You are awesome.

Out[5]:
True

Adjoint test

In [6]:
v = np.random.rand(mesh.nC)
w = np.random.rand(pred.shape[0])
wtJv = w.dot(prob.Jvec(m0, v))
vtJtw = v.dot(prob.Jtvec(m0, w))
passed = np.abs(wtJv - vtJtw) < 1e-10
print('Adjoint Test', np.abs(wtJv - vtJtw), passed)
Adjoint Test 8.881784197001252e-16 True

Ready to run an inversion

Once you have created the Survey and Problem classes so that you can compute predicted data (dpred), forward and adjoint sensitivity-vector products (Jvec and Jtvec), you are all set to run inversion using SimPEG's inverersion framework. Check out the examples in 3_MT1D_5layer_inversion and 1_MT1D_NumericalSetup.