from MT1D import MT1DProblem, MT1DSurvey, MT1DSrc, ZxyRx, Survey, AppResPhaRx
from SimPEG import Maps, DataMisfit
from scipy.constants import mu_0
import numpy as np
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.
For more details, see the SimPEG docs
Mesh
: mesh geometry and differential operatorsProblem
: physics engine. contains the machinery to construct and solve the PDESurvey
: sources and receiversSrc
: sourcesRx
: receiversProblem
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
.
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.
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}$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
We perform both order and adjoint tests shown in Appendix_A_MT1D_Sensitivity using 1D MT problem in SimPEG
's framework.
Jvec
¶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.
True
Jtvec
¶# 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||
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.
True
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
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.