mirror of
https://github.com/wassname/simpeg.git
synced 2026-06-30 10:52:27 +08:00
Merge branch 'master' of https://github.com/simpeg/simpegem
This commit is contained in:
SEOGI KANG
@@ -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:
|
||||
|
||||
+46
-11
@@ -16,17 +16,6 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
|
||||
|
||||
surveyPair = SurveyTDEM
|
||||
|
||||
def calcFields(self, sol, tInd):
|
||||
|
||||
if self.solType == 'b':
|
||||
b = sol
|
||||
e = self.MeSigmaI*self.mesh.edgeCurl.T*self.MfMui*b
|
||||
# Todo: implement non-zero js
|
||||
else:
|
||||
raise NotImplementedError('solType "%s" is not implemented in CalcFields.' % self.solType)
|
||||
|
||||
return {'b':b, 'e':e}
|
||||
|
||||
def fields(self, m):
|
||||
self.curModel = m
|
||||
# Create a fields storage object
|
||||
@@ -73,3 +62,49 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
|
||||
F[:,:,tInd+1] = CalcFields(sol, tInd)
|
||||
return F
|
||||
|
||||
def Jvec(self, m, v, u=None):
|
||||
"""
|
||||
: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):
|
||||
"""
|
||||
: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)
|
||||
|
||||
p = self.survey.projectFieldsDeriv(u, v=v, adjoint=True)
|
||||
y = self.solveAht(m, p)
|
||||
w = self.Gtvec(m, y, u)
|
||||
return - mkvc(w)
|
||||
|
||||
|
||||
@@ -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
|
||||
|
||||
+89
-72
@@ -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
|
||||
|
||||
@@ -41,31 +41,22 @@ class ProblemTDEM_b(BaseTDEMProblem):
|
||||
RHS = (1.0/dt)*self.MfMui*B_n
|
||||
return RHS
|
||||
|
||||
def calcFields(self, sol, tInd):
|
||||
|
||||
if self.solType == 'b':
|
||||
b = sol
|
||||
e = self.MeSigmaI*self.mesh.edgeCurl.T*self.MfMui*b
|
||||
# Todo: implement non-zero js
|
||||
else:
|
||||
raise NotImplementedError('solType "%s" is not implemented in CalcFields.' % self.solType)
|
||||
|
||||
return {'b':b, 'e':e}
|
||||
|
||||
|
||||
####################################################
|
||||
# Derivatives
|
||||
####################################################
|
||||
|
||||
def Jvec(self, m, v, u=None):
|
||||
if u is None:
|
||||
u = 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)
|
||||
|
||||
if not isinstance(v, self.dataPair):
|
||||
v = self.dataPair(self.survey, v)
|
||||
|
||||
p = self.survey.projectFieldsDeriv(u, v=v, adjoint=True)
|
||||
y = self.solveAht(m, p)
|
||||
w = self.Gtvec(m, y, u)
|
||||
return - mkvc(w)
|
||||
|
||||
def Gvec(self, m, vec, u=None):
|
||||
"""
|
||||
:param numpy.array m: Conductivity model
|
||||
@@ -76,8 +67,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 +94,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 +109,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 +160,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:
|
||||
#
|
||||
@@ -148,10 +201,10 @@ class ProblemTDEM_b(BaseTDEMProblem):
|
||||
# ^
|
||||
# fLoc 0 1 2 3
|
||||
# |-----|-----|-----|
|
||||
# tInd 0 1 2 / /
|
||||
# / __/
|
||||
# tInd 0 1 2
|
||||
# / ___/
|
||||
# 2 (tInd=2 uses fields 3 and would use 4 but it doesn't exist)
|
||||
# / __/
|
||||
# / ___/
|
||||
# 1 (tInd=1 uses fields 2 and 3)
|
||||
|
||||
def AhtRHS(tInd, y):
|
||||
@@ -184,7 +237,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 +282,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
|
||||
@@ -264,46 +317,10 @@ class ProblemTDEM_b(BaseTDEMProblem):
|
||||
"""
|
||||
self.curModel = m
|
||||
f = FieldsTDEM(self.mesh, self.survey)
|
||||
for i in range(1,self.nT+1):
|
||||
b = 1.0/self.timeSteps[i-1]*self.MfMui*vec[:,'b',i] + self.MfMui*self.mesh.edgeCurl*vec[:,'e',i]
|
||||
if i < self.nT:
|
||||
b = b - 1.0/self.timeSteps[i]*self.MfMui*vec[:,'b',i+1]
|
||||
f[:,'b', i] = b
|
||||
f[:,'e', i] = self.mesh.edgeCurl.T*self.MfMui*vec[:,'b',i] - self.MeSigma*vec[:,'e',i]
|
||||
for i in range(self.nT):
|
||||
b = 1.0/self.timeSteps[i]*self.MfMui*vec[:,'b',i+1] + self.MfMui*self.mesh.edgeCurl*vec[:,'e',i+1]
|
||||
if i < self.nT-1:
|
||||
b = b - 1.0/self.timeSteps[i+1]*self.MfMui*vec[:,'b',i+2]
|
||||
f[:,'b', i+1] = b
|
||||
f[:,'e', i+1] = self.mesh.edgeCurl.T*self.MfMui*vec[:,'b',i+1] - self.MeSigma*vec[:,'e',i+1]
|
||||
return f
|
||||
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
from SimPEG import *
|
||||
import simpegEM as EM
|
||||
from simpegEM.Utils.Ana import hzAnalyticDipoleT
|
||||
from scipy.constants import mu_0
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
cs, ncx, ncz, npad = 5., 20, 6, 20
|
||||
hx = [(cs, ncx), (cs, npad, 1.3)]
|
||||
hz = [(cs, npad, -1.3), (cs, ncz), (cs, npad, 1.3)]
|
||||
mesh = Mesh.CylMesh([hx,1,hz], '00C')
|
||||
mapping = Maps.Vertical1DMap(mesh)
|
||||
|
||||
opts = {'txLoc':0.,
|
||||
'txType':'VMD_MVP',
|
||||
'rxLoc':np.r_[150., 0., 0.],
|
||||
'rxType':'bz',
|
||||
'timeCh':np.logspace(-4,-2,20),
|
||||
}
|
||||
survey = EM.TDEM.SurveyTDEM1D(**opts)
|
||||
|
||||
prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping)
|
||||
# prb.setTimes([1e-5, 5e-5, 2.5e-4], [150, 150, 150])
|
||||
# prb.setTimes([1e-5, 5e-5, 2.5e-4], [10, 10, 10])
|
||||
prb.timeSteps = [(1e-5, 10)]
|
||||
prb.pair(survey)
|
||||
m = np.random.rand(mesh.nCz)
|
||||
|
||||
print survey.dpred(m)
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -55,7 +55,7 @@ class FDEM_analyticTests(unittest.TestCase):
|
||||
|
||||
an = EM.Utils.Ana.FEM.hzAnalyticDipoleF(x, self.Tx0.freq, self.sig)
|
||||
|
||||
diff = np.log10(np.abs(P*np.imag(u[self.Tx0, 'b']) - np.abs(mu_0*np.imag(an))))
|
||||
diff = np.log10(np.abs(P*np.imag(u[self.Tx0, 'b']) - mu_0*np.imag(an)))
|
||||
|
||||
if plotIt:
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
@@ -22,9 +22,8 @@ class TDEM_bDerivTests(unittest.TestCase):
|
||||
[Maps.ExpMap, Maps.Vertical1DMap, activeMap])
|
||||
|
||||
rxOffset = 40.
|
||||
rxTypes = 'bx,bz'
|
||||
rxs = [EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20), rxType) for rxType in rxTypes.split(',')]
|
||||
tx = EM.TDEM.TxTDEM(np.array([0., 0., 0.]), 'VMD_MVP', rxs)
|
||||
rx = EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20), 'bz')
|
||||
tx = EM.TDEM.TxTDEM(np.array([0., 0., 0.]), 'VMD_MVP', [rx])
|
||||
|
||||
survey = EM.TDEM.SurveyTDEM([tx])
|
||||
|
||||
@@ -49,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]
|
||||
@@ -88,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)
|
||||
|
||||
@@ -131,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()
|
||||
@@ -150,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)
|
||||
@@ -222,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):
|
||||
@@ -237,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)
|
||||
|
||||
@@ -0,0 +1,76 @@
|
||||
import unittest
|
||||
from SimPEG import *
|
||||
import simpegEM as EM
|
||||
|
||||
plotIt = False
|
||||
|
||||
def getProb(meshType='CYL',rxTypes='bx,bz',nTx=1):
|
||||
cs = 5.
|
||||
ncx = 20
|
||||
ncy = 6
|
||||
npad = 20
|
||||
hx = [(cs,ncx), (cs,npad,1.3)]
|
||||
hy = [(cs,npad,-1.3), (cs,ncy), (cs,npad,1.3)]
|
||||
mesh = Mesh.CylMesh([hx,1,hy], '00C')
|
||||
|
||||
active = mesh.vectorCCz<0.
|
||||
activeMap = Maps.ActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
|
||||
mapping = Maps.ComboMap(mesh, [Maps.ExpMap, Maps.Vertical1DMap, activeMap])
|
||||
|
||||
rxOffset = 40.
|
||||
|
||||
txs = []
|
||||
for ii in range(nTx):
|
||||
rxs = [EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20 + ii), rxType) for rxType in rxTypes.split(',')]
|
||||
txs += [EM.TDEM.TxTDEM(np.array([0., 0., 0.]), 'VMD_MVP', rxs)]
|
||||
|
||||
survey = EM.TDEM.SurveyTDEM(txs)
|
||||
|
||||
prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping)
|
||||
# prb.timeSteps = [1e-5]
|
||||
prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)]
|
||||
# prb.timeSteps = [(1e-05, 100)]
|
||||
|
||||
sigma = np.ones(mesh.nCz)*1e-8
|
||||
sigma[mesh.vectorCCz<0] = 1e-1
|
||||
sigma = np.log(sigma[active])
|
||||
|
||||
prb.pair(survey)
|
||||
return prb, mesh, sigma
|
||||
|
||||
def dotestJvec(prb, mesh, sigma):
|
||||
prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)]
|
||||
# d_sig = 0.8*sigma #np.random.rand(mesh.nCz)
|
||||
d_sig = 10*np.random.rand(prb.mapping.nP)
|
||||
derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(sigma, mx)]
|
||||
return Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=2, eps=1e-20)
|
||||
|
||||
def dotestAdjoint(prb, mesh, sigma):
|
||||
m = np.random.rand(prb.mapping.nP)
|
||||
d = np.random.rand(prb.survey.nD)
|
||||
|
||||
V1 = d.dot(prb.Jvec(sigma, m))
|
||||
V2 = m.dot(prb.Jtvec(sigma, d))
|
||||
print 'AdjointTest', V1, V2
|
||||
return np.abs(V1-V2)/np.abs(V1), 1e-6
|
||||
|
||||
class TDEM_bDerivTests(unittest.TestCase):
|
||||
|
||||
def test_Jvec_bx(self): self.assertTrue(dotestJvec(*getProb(rxTypes='bx')))
|
||||
def test_Adjoint_bx(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='bx')))
|
||||
|
||||
def test_Jvec_bxbz(self): self.assertTrue(dotestJvec(*getProb(rxTypes='bx,bz')))
|
||||
def test_Adjoint_bxbz(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='bx,bz')))
|
||||
|
||||
def test_Jvec_bxbz_2tx(self): self.assertTrue(dotestJvec(*getProb(rxTypes='bx,bz',nTx=2)))
|
||||
def test_Adjoint_bxbz_2tx(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='bx,bz',nTx=2)))
|
||||
|
||||
def test_Jvec_bxbzbz(self): self.assertTrue(dotestJvec(*getProb(rxTypes='bx,bz,bz')))
|
||||
def test_Adjoint_bxbzbz(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='bx,bz,bz')))
|
||||
|
||||
# def test_Jvec_ey(self): self.assertTrue(dotestJvec(*getProb(rxTypes='ey')))
|
||||
# def test_Adjoint_ey(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='ey')))
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
@@ -36,7 +36,7 @@ def MagneticDipoleVectorPotential(txLoc, obsLoc, component, dipoleMoment=(0., 0.
|
||||
dR = obsLoc - txLoc[i, np.newaxis].repeat(nEdges, axis=0)
|
||||
mCr = np.cross(m, dR)
|
||||
r = np.sqrt((dR**2).sum(axis=1))
|
||||
A[:, i] = -(mu_0/(4*pi)) * mCr[:,dimInd]/(r**3)
|
||||
A[:, i] = +(mu_0/(4*pi)) * mCr[:,dimInd]/(r**3)
|
||||
if nTx == 1:
|
||||
return A.flatten()
|
||||
return A
|
||||
|
||||
Reference in New Issue
Block a user