mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-13 17:45:30 +08:00
Documentation updates.
This commit is contained in:
rowanc1
@@ -107,4 +107,3 @@ FDEM Survey
|
||||
:show-inheritance:
|
||||
:members:
|
||||
:undoc-members:
|
||||
|
||||
|
||||
+2
-3
@@ -55,13 +55,12 @@ TDEM - B formulation
|
||||
:show-inheritance:
|
||||
:members:
|
||||
:undoc-members:
|
||||
:inherited-members:
|
||||
|
||||
|
||||
Field Storage
|
||||
=============
|
||||
|
||||
.. automodule:: simpegEM.TDEM.FieldsTDEM
|
||||
.. autoclass:: simpegEM.TDEM.SurveyTDEM.FieldsTDEM
|
||||
:show-inheritance:
|
||||
:members:
|
||||
:undoc-members:
|
||||
@@ -71,7 +70,7 @@ Field Storage
|
||||
TDEM Survey Classes
|
||||
===================
|
||||
|
||||
.. automodule:: simpegEM.TDEM.SurveyTDEM
|
||||
.. autoclass:: simpegEM.TDEM.SurveyTDEM.SurveyTDEM
|
||||
:show-inheritance:
|
||||
:members:
|
||||
:undoc-members:
|
||||
|
||||
@@ -135,129 +135,3 @@ class SurveyTDEM(Survey.BaseSurvey):
|
||||
return f
|
||||
|
||||
|
||||
|
||||
# class SurveyTDEM1D(BaseSurvey):
|
||||
# """
|
||||
# docstring for SurveyTDEM1D
|
||||
# """
|
||||
|
||||
# txLoc = None #: txLoc
|
||||
# txType = None #: txType
|
||||
# rxLoc = None #: rxLoc
|
||||
# rxType = None #: rxType
|
||||
# timeCh = None #: timeCh
|
||||
# nTx = 1 #: Number of transmitters
|
||||
|
||||
# @property
|
||||
# def nTimeCh(self):
|
||||
# """Number of time channels"""
|
||||
# return self.timeCh.size
|
||||
|
||||
# def __init__(self, **kwargs):
|
||||
# BaseSurvey.__init__(self, **kwargs)
|
||||
# Utils.setKwargs(self, **kwargs)
|
||||
|
||||
# def projectFields(self, u):
|
||||
# #TODO: this is hardcoded to 1Tx
|
||||
# return self.Qrx.dot(u.b[:,:,0].T).T
|
||||
|
||||
# def projectFieldsAdjoint(self, d):
|
||||
# # TODO: make the following self.nTimeCh
|
||||
# d = d.reshape((self.prob.nT, self.nTx), order='F')
|
||||
# #TODO: *Qtime.T need to multiply by a time projection. (outside for loop??)
|
||||
# ii = 0
|
||||
# F = FieldsTDEM(self.prob.mesh, self.nTx, self.prob.nT, 'b')
|
||||
# for ii in range(self.prob.nT):
|
||||
# b = self.Qrx.T*d[ii,:]
|
||||
# F.set_b(b, ii)
|
||||
# F.set_e(np.zeros((self.prob.mesh.nE,self.nTx)), ii)
|
||||
# return F
|
||||
|
||||
# ####################################################
|
||||
# # Interpolation Matrices
|
||||
# ####################################################
|
||||
|
||||
# @property
|
||||
# def Qrx(self):
|
||||
# if self._Qrx is None:
|
||||
# if self.rxType == 'bz':
|
||||
# locType = 'Fz'
|
||||
# self._Qrx = self.prob.mesh.getInterpolationMat(self.rxLoc, locType=locType)
|
||||
# return self._Qrx
|
||||
# _Qrx = None
|
||||
|
||||
|
||||
# class FieldsTDEM_OLD(object):
|
||||
# """docstring for FieldsTDEM"""
|
||||
|
||||
# phi0 = None #: Initial electric potential
|
||||
# A0 = None #: Initial magnetic vector potential
|
||||
# e0 = None #: Initial electric field
|
||||
# b0 = None #: Initial magnetic flux density
|
||||
# j0 = None #: Initial current density
|
||||
# h0 = None #: Initial magnetic field
|
||||
|
||||
# phi = None #: Electric potential
|
||||
# A = None #: Magnetic vector potential
|
||||
# e = None #: Electric field
|
||||
# b = None #: Magnetic flux density
|
||||
# j = None #: Current density
|
||||
# h = None #: Magnetic field
|
||||
|
||||
# def __init__(self, mesh, nTx, nT, store='b'):
|
||||
|
||||
# self.nT = nT #: Number of times
|
||||
# self.nTx = nTx #: Number of transmitters
|
||||
# self.mesh = mesh
|
||||
|
||||
# def update(self, newFields, tInd):
|
||||
# self.set_b(newFields['b'], tInd)
|
||||
# self.set_e(newFields['e'], tInd)
|
||||
|
||||
# def fieldVec(self):
|
||||
# u = np.ndarray((0, self.nTx))
|
||||
# for i in range(self.nT):
|
||||
# u = np.r_[u, self.get_b(i), self.get_e(i)]
|
||||
# if self.nTx == 1:
|
||||
# u = u.flatten()
|
||||
# return u
|
||||
|
||||
# ####################################################
|
||||
# # Get Methods
|
||||
# ####################################################
|
||||
|
||||
# def get_b(self, ind):
|
||||
# if ind == -1:
|
||||
# return self.b0
|
||||
# else:
|
||||
# return self.b[ind,:,:]
|
||||
|
||||
# def get_e(self, ind):
|
||||
# if ind == -1:
|
||||
# return self.e0
|
||||
# else:
|
||||
# return self.e[ind,:,:]
|
||||
|
||||
# ####################################################
|
||||
# # Set Methods
|
||||
# ####################################################
|
||||
|
||||
# def set_b(self, b, ind):
|
||||
# if self.b is None:
|
||||
# self.b = np.zeros((self.nT, np.sum(self.mesh.nF), self.nTx))
|
||||
# self.b[:] = np.nan
|
||||
# if len(b.shape) == 1:
|
||||
# b = b[:, np.newaxis]
|
||||
# self.b[ind,:,:] = b
|
||||
|
||||
# def set_e(self, e, ind):
|
||||
# if self.e is None:
|
||||
# self.e = np.zeros((self.nT, np.sum(self.mesh.nE), self.nTx))
|
||||
# self.e[:] = np.nan
|
||||
# if len(e.shape) == 1:
|
||||
# e = e[:, np.newaxis]
|
||||
# self.e[ind,:,:] = e
|
||||
|
||||
|
||||
# def __contains__(self, key):
|
||||
# return key in self.children
|
||||
|
||||
+99
-11
@@ -18,7 +18,7 @@ class ProblemTDEM_b(BaseTDEMProblem):
|
||||
def __init__(self, mesh, mapping=None, **kwargs):
|
||||
BaseTDEMProblem.__init__(self, mesh, mapping=mapping, **kwargs)
|
||||
|
||||
solType = 'b'
|
||||
solType = 'b' #: Type of the solution, in this case the 'b' field
|
||||
|
||||
surveyPair = SurveyTDEM
|
||||
|
||||
@@ -47,16 +47,42 @@ class ProblemTDEM_b(BaseTDEMProblem):
|
||||
####################################################
|
||||
|
||||
def Jvec(self, m, v, u=None):
|
||||
if u is None:
|
||||
u = self.fields(m)
|
||||
"""
|
||||
:param numpy.array m: Conductivity model
|
||||
:param numpy.ndarray v: vector (model object)
|
||||
:param simpegEM.TDEM.FieldsTDEM u: Fields resulting from m
|
||||
:rtype: numpy.ndarray
|
||||
:return: w (data object)
|
||||
|
||||
Multiplying \\\(\\\mathbf{J}\\\) onto a vector can be broken into three steps
|
||||
|
||||
* Compute \\\(\\\\vec{p} = \\\mathbf{G}v\\\)
|
||||
* Solve \\\(\\\hat{\\\mathbf{A}} \\\\vec{y} = \\\\vec{p}\\\)
|
||||
* Compute \\\(\\\\vec{w} = -\\\mathbf{Q} \\\\vec{y}\\\)
|
||||
|
||||
"""
|
||||
u = u or self.fields(m)
|
||||
p = self.Gvec(m, v, u)
|
||||
y = self.solveAh(m, p)
|
||||
Jv = self.survey.projectFieldsDeriv(u, v=y)
|
||||
return - mkvc(Jv)
|
||||
|
||||
def Jtvec(self, m, v, u=None):
|
||||
if u is None:
|
||||
u = self.fields(m)
|
||||
"""
|
||||
:param numpy.array m: Conductivity model
|
||||
:param numpy.ndarray,SimPEG.Survey.Data v: vector (data object)
|
||||
:param simpegEM.TDEM.FieldsTDEM u: Fields resulting from m
|
||||
:rtype: numpy.ndarray
|
||||
:return: w (model object)
|
||||
|
||||
Multiplying \\\(\\\mathbf{J}^\\\\top\\\) onto a vector can be broken into three steps
|
||||
|
||||
* Compute \\\(\\\\vec{p} = \\\mathbf{Q}^\\\\top \\\\vec{v}\\\)
|
||||
* Solve \\\(\\\hat{\\\mathbf{A}}^\\\\top \\\\vec{y} = \\\\vec{p}\\\)
|
||||
* Compute \\\(\\\\vec{w} = -\\\mathbf{G}^\\\\top y\\\)
|
||||
|
||||
"""
|
||||
u = u or self.fields(m)
|
||||
|
||||
if not isinstance(v, self.dataPair):
|
||||
v = self.dataPair(self.survey, v)
|
||||
@@ -76,8 +102,7 @@ class ProblemTDEM_b(BaseTDEMProblem):
|
||||
|
||||
Multiply G by a vector
|
||||
"""
|
||||
if u is None:
|
||||
u = self.fields(m)
|
||||
u = u or self.fields(m)
|
||||
|
||||
# Note: Fields has shape (nF/E, nTx, nT+1)
|
||||
# However, p will only really fill (:,:,1:nT+1)
|
||||
@@ -104,8 +129,7 @@ class ProblemTDEM_b(BaseTDEMProblem):
|
||||
|
||||
Multiply G.T by a vector
|
||||
"""
|
||||
if u is None:
|
||||
u = self.fields(m)
|
||||
u = u or self.fields(m)
|
||||
nTx, nE = self.survey.nTx, self.mesh.nE
|
||||
tmp = np.zeros(nE)
|
||||
# Here we can do internal multiplications of Gt*v and then multiply by MsigDeriv.T in one go.
|
||||
@@ -120,6 +144,38 @@ class ProblemTDEM_b(BaseTDEMProblem):
|
||||
return p
|
||||
|
||||
def solveAh(self, m, p):
|
||||
"""
|
||||
:param numpy.array m: Conductivity model
|
||||
:param simpegEM.TDEM.FieldsTDEM p: Fields object
|
||||
:rtype: simpegEM.TDEM.FieldsTDEM
|
||||
:return: y
|
||||
|
||||
Solve the block-matrix system \\\(\\\hat{A} \\\hat{y} = \\\hat{p}\\\):
|
||||
|
||||
.. math::
|
||||
\mathbf{\hat{A}} = \left[
|
||||
\\begin{array}{cccc}
|
||||
A & 0 & & \\\\
|
||||
B & A & & \\\\
|
||||
& \ddots & \ddots & \\\\
|
||||
& & B & A
|
||||
\end{array}
|
||||
\\right] \\\\
|
||||
\mathbf{A} =
|
||||
\left[
|
||||
\\begin{array}{cc}
|
||||
\\frac{1}{\delta t} \MfMui & \MfMui\dcurl \\\\
|
||||
\dcurl^\\top \MfMui & -\MeSig
|
||||
\end{array}
|
||||
\\right] \\\\
|
||||
\mathbf{B} =
|
||||
\left[
|
||||
\\begin{array}{cc}
|
||||
-\\frac{1}{\delta t} \MfMui & 0 \\\\
|
||||
0 & 0
|
||||
\end{array}
|
||||
\\right] \\\\
|
||||
"""
|
||||
|
||||
def AhRHS(tInd, y):
|
||||
rhs = self.MfMui*self.mesh.edgeCurl*self.MeSigmaI*p[:,'e',tInd+1] + p[:,'b',tInd+1]
|
||||
@@ -139,6 +195,38 @@ class ProblemTDEM_b(BaseTDEMProblem):
|
||||
return self.forward(m, AhRHS, AhCalcFields)
|
||||
|
||||
def solveAht(self, m, p):
|
||||
"""
|
||||
:param numpy.array m: Conductivity model
|
||||
:param simpegEM.TDEM.FieldsTDEM p: Fields object
|
||||
:rtype: simpegEM.TDEM.FieldsTDEM
|
||||
:return: y
|
||||
|
||||
Solve the block-matrix system \\\(\\\hat{A}^\\\\top \\\hat{y} = \\\hat{p}\\\):
|
||||
|
||||
.. math::
|
||||
\mathbf{\hat{A}}^\\top = \left[
|
||||
\\begin{array}{cccc}
|
||||
A & B & & \\\\
|
||||
& \ddots & \ddots & \\\\
|
||||
& & A & B \\\\
|
||||
& & 0 & A
|
||||
\end{array}
|
||||
\\right] \\\\
|
||||
\mathbf{A} =
|
||||
\left[
|
||||
\\begin{array}{cc}
|
||||
\\frac{1}{\delta t} \MfMui & \MfMui\dcurl \\\\
|
||||
\dcurl^\\top \MfMui & -\MeSig
|
||||
\end{array}
|
||||
\\right] \\\\
|
||||
\mathbf{B} =
|
||||
\left[
|
||||
\\begin{array}{cc}
|
||||
-\\frac{1}{\delta t} \MfMui & 0 \\\\
|
||||
0 & 0
|
||||
\end{array}
|
||||
\\right] \\\\
|
||||
"""
|
||||
|
||||
# Mini Example:
|
||||
#
|
||||
@@ -184,7 +272,7 @@ class ProblemTDEM_b(BaseTDEMProblem):
|
||||
# Functions for tests
|
||||
####################################################
|
||||
|
||||
def AhVec(self, m, vec):
|
||||
def _AhVec(self, m, vec):
|
||||
"""
|
||||
:param numpy.array m: Conductivity model
|
||||
:param simpegEM.TDEM.FieldsTDEM vec: Fields object
|
||||
@@ -229,7 +317,7 @@ class ProblemTDEM_b(BaseTDEMProblem):
|
||||
f[:,'e',i] = self.mesh.edgeCurl.T*self.MfMui*vec[:,'b',i] - self.MeSigma*vec[:,'e',i]
|
||||
return f
|
||||
|
||||
def AhtVec(self, m, vec):
|
||||
def _AhtVec(self, m, vec):
|
||||
"""
|
||||
:param numpy.array m: Conductivity model
|
||||
:param simpegEM.TDEM.FieldsTDEM vec: Fields object
|
||||
|
||||
@@ -48,7 +48,7 @@ class TDEM_bDerivTests(unittest.TestCase):
|
||||
sigma = self.sigma
|
||||
|
||||
u = prb.fields(sigma)
|
||||
Ahu = prb.AhVec(sigma, u)
|
||||
Ahu = prb._AhVec(sigma, u)
|
||||
|
||||
V1 = Ahu[:,'b',1]
|
||||
V2 = 1./prb.timeSteps[0]*prb.MfMui*u[:,'b',0]
|
||||
@@ -87,7 +87,7 @@ class TDEM_bDerivTests(unittest.TestCase):
|
||||
|
||||
f = prb.fields(sigma)
|
||||
u1 = A*f.tovec()
|
||||
u2 = prb.AhVec(sigma,f).tovec()
|
||||
u2 = prb._AhVec(sigma,f).tovec()
|
||||
|
||||
self.assertTrue(np.linalg.norm(u1-u2)/np.linalg.norm(u1)<1e-12)
|
||||
|
||||
@@ -130,7 +130,7 @@ class TDEM_bDerivTests(unittest.TestCase):
|
||||
for i in range(prb.nT):
|
||||
f[:,'e', i] = np.random.rand(mesh.nE, 1)
|
||||
|
||||
Ahf = prb.AhVec(sigma, f)
|
||||
Ahf = prb._AhVec(sigma, f)
|
||||
f_test = prb.solveAh(sigma, Ahf)
|
||||
|
||||
u1 = f.tovec()
|
||||
@@ -149,7 +149,7 @@ class TDEM_bDerivTests(unittest.TestCase):
|
||||
dm = 1000*np.random.rand(self.prb.mapping.nP)
|
||||
h = 0.01
|
||||
|
||||
derChk = lambda m: [self.prb.AhVec(m, f).tovec(), lambda mx: self.prb.Gvec(sigma, mx, u=f).tovec()]
|
||||
derChk = lambda m: [self.prb._AhVec(m, f).tovec(), lambda mx: self.prb.Gvec(sigma, mx, u=f).tovec()]
|
||||
print '\ntest_DerivG'
|
||||
passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20)
|
||||
self.assertTrue(passed)
|
||||
@@ -221,8 +221,8 @@ class TDEM_bDerivTests(unittest.TestCase):
|
||||
f2[:,'b',i] = np.random.rand(mesh.nF, 1)
|
||||
f2[:,'e',i] = np.random.rand(mesh.nE, 1)
|
||||
|
||||
V1 = f2.tovec().dot(prb.AhVec(sigma, f1).tovec())
|
||||
V2 = f1.tovec().dot(prb.AhtVec(sigma, f2).tovec())
|
||||
V1 = f2.tovec().dot(prb._AhVec(sigma, f1).tovec())
|
||||
V2 = f1.tovec().dot(prb._AhtVec(sigma, f2).tovec())
|
||||
self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6)
|
||||
|
||||
# def test_solveAhtVsAhtVec(self):
|
||||
@@ -236,7 +236,7 @@ class TDEM_bDerivTests(unittest.TestCase):
|
||||
# f1[:,'e',i] = np.random.rand(mesh.nE, 1)
|
||||
|
||||
# f2 = prb.solveAht(sigma, f1)
|
||||
# f3 = prb.AhtVec(sigma, f2)
|
||||
# f3 = prb._AhtVec(sigma, f2)
|
||||
|
||||
# if True:
|
||||
# import matplotlib.pyplot as plt
|
||||
|
||||
@@ -53,7 +53,7 @@ class TDEM_bDerivTests(unittest.TestCase):
|
||||
dm = 1000*np.random.rand(self.prb.mapping.nP)
|
||||
h = 0.01
|
||||
|
||||
derChk = lambda m: [self.prb.AhVec(m, f).tovec(), lambda mx: self.prb.Gvec(sigma, mx, u=f).tovec()]
|
||||
derChk = lambda m: [self.prb._AhVec(m, f).tovec(), lambda mx: self.prb.Gvec(sigma, mx, u=f).tovec()]
|
||||
print '\ntest_DerivG'
|
||||
passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20)
|
||||
self.assertTrue(passed)
|
||||
|
||||
Reference in New Issue
Block a user