wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-07-17 11:32:59 +08:00
Code Issues Packages Projects Releases Wiki Activity

s_m, s_e are vectors (so they should not be capitalized)

Browse Source
This commit is contained in:
Lindsey Heagy
2016-03-20 13:04:02 -07:00
parent cc9d2e5ac7
commit c51afa4aad
3 changed files with 206 additions and 206 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/EM/FDEM/FDEM.py
+31 -31
View File
@@ -155,22 +155,22 @@ class BaseFDEMProblem(BaseEMProblem):
:param float freq: Frequency
:rtype: (numpy.ndarray, numpy.ndarray)
:return: S_m, S_e (nE or nF, nSrc)
:return: s_m, s_e (nE or nF, nSrc)
"""
Srcs = self.survey.getSrcByFreq(freq)
if self._formulation is 'EB':
S_m = np.zeros((self.mesh.nF,len(Srcs)), dtype=complex)
S_e = np.zeros((self.mesh.nE,len(Srcs)), dtype=complex)
s_m = np.zeros((self.mesh.nF,len(Srcs)), dtype=complex)
s_e = np.zeros((self.mesh.nE,len(Srcs)), dtype=complex)
elif self._formulation is 'HJ':
S_m = np.zeros((self.mesh.nE,len(Srcs)), dtype=complex)
S_e = np.zeros((self.mesh.nF,len(Srcs)), dtype=complex)
s_m = np.zeros((self.mesh.nE,len(Srcs)), dtype=complex)
s_e = np.zeros((self.mesh.nF,len(Srcs)), dtype=complex)
for i, src in enumerate(Srcs):
smi, sei = src.eval(self)
S_m[:,i] = S_m[:,i] + smi
S_e[:,i] = S_e[:,i] + sei
s_m[:,i] = s_m[:,i] + smi
s_e[:,i] = s_e[:,i] + sei
return S_m, S_e
return s_m, s_e
##########################################################################################
@@ -258,11 +258,11 @@ class Problem_e(BaseFDEMProblem):
:return: RHS (nE, nSrc)
"""
S_m, S_e = self.getSourceTerm(freq)
s_m, s_e = self.getSourceTerm(freq)
C = self.mesh.edgeCurl
MfMui = self.MfMui
return C.T * (MfMui * S_m) -1j * omega(freq) * S_e
return C.T * (MfMui * s_m) -1j * omega(freq) * s_e
def getRHSDeriv(self, freq, src, v, adjoint=False):
"""
@@ -278,14 +278,14 @@ class Problem_e(BaseFDEMProblem):
C = self.mesh.edgeCurl
MfMui = self.MfMui
S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint=adjoint)
s_mDeriv, s_eDeriv = src.evalDeriv(self, adjoint=adjoint)
if adjoint:
dRHS = MfMui * (C * v)
return S_mDeriv(dRHS) - 1j * omega(freq) * S_eDeriv(v)
return s_mDeriv(dRHS) - 1j * omega(freq) * s_eDeriv(v)
else:
return C.T * (MfMui * S_mDeriv(v)) -1j * omega(freq) * S_eDeriv(v)
return C.T * (MfMui * s_mDeriv(v)) -1j * omega(freq) * s_eDeriv(v)
class Problem_b(BaseFDEMProblem):
@@ -383,11 +383,11 @@ class Problem_b(BaseFDEMProblem):
:return: RHS (nE, nSrc)
"""
S_m, S_e = self.getSourceTerm(freq)
s_m, s_e = self.getSourceTerm(freq)
C = self.mesh.edgeCurl
MeSigmaI = self.MeSigmaI
RHS = S_m + C * ( MeSigmaI * S_e )
RHS = s_m + C * ( MeSigmaI * s_e )
if self._makeASymmetric is True:
MfMui = self.MfMui
@@ -408,21 +408,21 @@ class Problem_b(BaseFDEMProblem):
"""
C = self.mesh.edgeCurl
S_m, S_e = src.eval(self)
s_m, s_e = src.eval(self)
MfMui = self.MfMui
if self._makeASymmetric and adjoint:
v = self.MfMui * v
MeSigmaIDeriv = self.MeSigmaIDeriv(S_e)
S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint=adjoint)
MeSigmaIDeriv = self.MeSigmaIDeriv(s_e)
s_mDeriv, s_eDeriv = src.evalDeriv(self, adjoint=adjoint)
if not adjoint:
RHSderiv = C * (MeSigmaIDeriv * v)
SrcDeriv = S_mDeriv(v) + C * (self.MeSigmaI * S_eDeriv(v))
SrcDeriv = s_mDeriv(v) + C * (self.MeSigmaI * s_eDeriv(v))
elif adjoint:
RHSderiv = MeSigmaIDeriv.T * (C.T * v)
SrcDeriv = S_mDeriv(v) + self.MeSigmaI.T * (C.T * S_eDeriv(v))
SrcDeriv = s_mDeriv(v) + self.MeSigmaI.T * (C.T * s_eDeriv(v))
if self._makeASymmetric is True and not adjoint:
return MfMui.T * (SrcDeriv + RHSderiv)
@@ -533,11 +533,11 @@ class Problem_j(BaseFDEMProblem):
:return: RHS
"""
S_m, S_e = self.getSourceTerm(freq)
s_m, s_e = self.getSourceTerm(freq)
C = self.mesh.edgeCurl
MeMuI = self.MeMuI
RHS = C * (MeMuI * S_m) - 1j * omega(freq) * S_e
RHS = C * (MeMuI * s_m) - 1j * omega(freq) * s_e
if self._makeASymmetric is True:
MfRho = self.MfRho
return MfRho.T*RHS
@@ -558,16 +558,16 @@ class Problem_j(BaseFDEMProblem):
C = self.mesh.edgeCurl
MeMuI = self.MeMuI
S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint=adjoint)
s_mDeriv, s_eDeriv = src.evalDeriv(self, adjoint=adjoint)
if adjoint:
if self._makeASymmetric:
MfRho = self.MfRho
v = MfRho*v
return S_mDeriv(MeMuI.T * (C.T * v)) - 1j * omega(freq) * S_eDeriv(v)
return s_mDeriv(MeMuI.T * (C.T * v)) - 1j * omega(freq) * s_eDeriv(v)
else:
RHSDeriv = C * (MeMuI * S_mDeriv(v)) - 1j * omega(freq) * S_eDeriv(v)
RHSDeriv = C * (MeMuI * s_mDeriv(v)) - 1j * omega(freq) * s_eDeriv(v)
if self._makeASymmetric:
MfRho = self.MfRho
@@ -655,11 +655,11 @@ class Problem_h(BaseFDEMProblem):
:return: RHS (nE, nSrc)
"""
S_m, S_e = self.getSourceTerm(freq)
s_m, s_e = self.getSourceTerm(freq)
C = self.mesh.edgeCurl
MfRho = self.MfRho
return S_m + C.T * ( MfRho * S_e )
return s_m + C.T * ( MfRho * s_e )
def getRHSDeriv(self, freq, src, v, adjoint=False):
"""
@@ -673,17 +673,17 @@ class Problem_h(BaseFDEMProblem):
:return: product of rhs deriv with a vector
"""
_, S_e = src.eval(self)
_, s_e = src.eval(self)
C = self.mesh.edgeCurl
MfRho = self.MfRho
MfRhoDeriv = self.MfRhoDeriv(S_e)
MfRhoDeriv = self.MfRhoDeriv(s_e)
if not adjoint:
RHSDeriv = C.T * (MfRhoDeriv * v)
elif adjoint:
RHSDeriv = MfRhoDeriv.T * (C * v)
S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint=adjoint)
s_mDeriv, s_eDeriv = src.evalDeriv(self, adjoint=adjoint)
return RHSDeriv + S_mDeriv(v) + C.T * (MfRho * S_eDeriv(v))
return RHSDeriv + s_mDeriv(v) + C.T * (MfRho * s_eDeriv(v))
SimPEG/EM/FDEM/FieldsFDEM.py
+127 -127
View File
@@ -8,7 +8,7 @@ from SimPEG.Utils import Zero, Identity, sdiag
class Fields(SimPEG.Problem.Fields):
"""
Fancy Field Storage for a FDEM survey. Only one field type is stored for
each problem, the rest are computed. The fields obejct acts like an array and is indexed by
@@ -34,56 +34,56 @@ class Fields(SimPEG.Problem.Fields):
def _e(self, solution, srcList):
"""
Total electric field is sum of primary and secondary
Total electric field is sum of primary and secondary
:param numpy.ndarray solution: field we solved for
:param list srcList: list of sources
:rtype: numpy.ndarray
:return: total electric field
"""
if getattr(self, '_ePrimary', None) is None or getattr(self, '_eSecondary', None) is None:
if getattr(self, '_ePrimary', None) is None or getattr(self, '_eSecondary', None) is None:
raise NotImplementedError ('Getting e from %s is not implemented' %self.knownFields.keys()[0])
return self._ePrimary(solution,srcList) + self._eSecondary(solution,srcList)
def _b(self, solution, srcList):
"""
Total magnetic flux density is sum of primary and secondary
Total magnetic flux density is sum of primary and secondary
:param numpy.ndarray solution: field we solved for
:param list srcList: list of sources
:rtype: numpy.ndarray
:return: total magnetic flux density
:return: total magnetic flux density
"""
if getattr(self, '_bPrimary', None) is None or getattr(self, '_bSecondary', None) is None:
if getattr(self, '_bPrimary', None) is None or getattr(self, '_bSecondary', None) is None:
raise NotImplementedError ('Getting b from %s is not implemented' %self.knownFields.keys()[0])
return self._bPrimary(solution, srcList) + self._bSecondary(solution, srcList)
def _h(self, solution, srcList):
"""
Total magnetic field is sum of primary and secondary
Total magnetic field is sum of primary and secondary
:param numpy.ndarray solution: field we solved for
:param list srcList: list of sources
:rtype: numpy.ndarray
:return: total magnetic field
"""
if getattr(self, '_hPrimary', None) is None or getattr(self, '_hSecondary', None) is None:
if getattr(self, '_hPrimary', None) is None or getattr(self, '_hSecondary', None) is None:
raise NotImplementedError ('Getting h from %s is not implemented' %self.knownFields.keys()[0])
return self._hPrimary(solution, srcList) + self._hSecondary(solution, srcList)
def _j(self, solution, srcList):
"""
Total current density is sum of primary and secondary
Total current density is sum of primary and secondary
:param numpy.ndarray solution: field we solved for
:param list srcList: list of sources
:rtype: numpy.ndarray
:return: total current density
:return: total current density
"""
if getattr(self, '_jPrimary', None) is None or getattr(self, '_jSecondary', None) is None:
if getattr(self, '_jPrimary', None) is None or getattr(self, '_jSecondary', None) is None:
raise NotImplementedError ('Getting j from %s is not implemented' %self.knownFields.keys()[0])
return self._jPrimary(solution, srcList) + self._jSecondary(solution, srcList)
@@ -99,7 +99,7 @@ class Fields(SimPEG.Problem.Fields):
:rtype: numpy.ndarray
:return: derivative times a vector (or tuple for adjoint)
"""
if getattr(self, '_eDeriv_u', None) is None or getattr(self, '_eDeriv_m', None) is None:
if getattr(self, '_eDeriv_u', None) is None or getattr(self, '_eDeriv_m', None) is None:
raise NotImplementedError ('Getting eDerivs from %s is not implemented' %self.knownFields.keys()[0])
if adjoint:
@@ -117,12 +117,12 @@ class Fields(SimPEG.Problem.Fields):
:rtype: numpy.ndarray
:return: derivative times a vector (or tuple for adjoint)
"""
if getattr(self, '_bDeriv_u', None) is None or getattr(self, '_bDeriv_m', None) is None:
if getattr(self, '_bDeriv_u', None) is None or getattr(self, '_bDeriv_m', None) is None:
raise NotImplementedError ('Getting bDerivs from %s is not implemented' %self.knownFields.keys()[0])
if adjoint:
return self._bDeriv_u(src, v, adjoint), self._bDeriv_m(src, v, adjoint)
return np.array(self._bDeriv_u(src, du_dm_v, adjoint) + self._bDeriv_m(src, v, adjoint), dtype = complex)
return np.array(self._bDeriv_u(src, du_dm_v, adjoint) + self._bDeriv_m(src, v, adjoint), dtype = complex)
def _hDeriv(self, src, du_dm_v, v, adjoint = False):
"""
@@ -135,10 +135,10 @@ class Fields(SimPEG.Problem.Fields):
:rtype: numpy.ndarray
:return: derivative times a vector (or tuple for adjoint)
"""
if getattr(self, '_hDeriv_u', None) is None or getattr(self, '_hDeriv_m', None) is None:
if getattr(self, '_hDeriv_u', None) is None or getattr(self, '_hDeriv_m', None) is None:
raise NotImplementedError ('Getting hDerivs from %s is not implemented' %self.knownFields.keys()[0])
if adjoint:
if adjoint:
return self._hDeriv_u(src, v, adjoint), self._hDeriv_m(src, v, adjoint)
return np.array(self._hDeriv_u(src, du_dm_v, adjoint) + self._hDeriv_m(src, v, adjoint), dtype = complex)
@@ -153,7 +153,7 @@ class Fields(SimPEG.Problem.Fields):
:rtype: numpy.ndarray
:return: derivative times a vector (or tuple for adjoint)
"""
if getattr(self, '_jDeriv_u', None) is None or getattr(self, '_jDeriv_m', None) is None:
if getattr(self, '_jDeriv_u', None) is None or getattr(self, '_jDeriv_m', None) is None:
raise NotImplementedError ('Getting jDerivs from %s is not implemented' %self.knownFields.keys()[0])
if adjoint:
@@ -162,10 +162,10 @@ class Fields(SimPEG.Problem.Fields):
class Fields_e(Fields):
"""
Fields object for Problem_e.
Fields object for Problem_e.
:param Mesh mesh: mesh
:param Survey survey: survey
:param Survey survey: survey
"""
knownFields = {'eSolution':'E'}
@@ -233,9 +233,9 @@ class Fields_e(Fields):
def _eDeriv_u(self, src, v, adjoint = False):
"""
Partial derivative of the total electric field with respect to the thing we
Partial derivative of the total electric field with respect to the thing we
solved for.
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
@@ -247,8 +247,8 @@ class Fields_e(Fields):
def _eDeriv_m(self, src, v, adjoint = False):
"""
Partial derivative of the total electric field with respect to the inversion model. Here, we assume that the primary does not depend on the model. Note that this also includes derivative contributions from the sources.
Partial derivative of the total electric field with respect to the inversion model. Here, we assume that the primary does not depend on the model. Note that this also includes derivative contributions from the sources.
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
@@ -289,14 +289,14 @@ class Fields_e(Fields):
b = (C * eSolution)
for i, src in enumerate(srcList):
b[:,i] *= - 1./(1j*omega(src.freq))
S_m, _ = src.eval(self.prob)
b[:,i] = b[:,i]+ 1./(1j*omega(src.freq)) * S_m
s_m, _ = src.eval(self.prob)
b[:,i] = b[:,i]+ 1./(1j*omega(src.freq)) * s_m
return b
def _bDeriv_u(self, src, du_dm_v, adjoint = False):
"""
Derivative of the magnetic flux density with respect to the thing we solved for
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray du_dm_v: vector to take product with
:param bool adjoint: adjoint?
@@ -312,8 +312,8 @@ class Fields_e(Fields):
def _bDeriv_m(self, src, v, adjoint = False):
"""
Derivative of the magnetic flux density with respect to the inversion model.
Derivative of the magnetic flux density with respect to the inversion model.
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
@@ -321,8 +321,8 @@ class Fields_e(Fields):
:return: product of the magnetic flux density derivative with respect to the inversion model with a vector
"""
S_mDeriv, _ = src.evalDeriv(self.prob, v, adjoint)
return 1./(1j * omega(src.freq)) * S_mDeriv
s_mDeriv, _ = src.evalDeriv(self.prob, v, adjoint)
return 1./(1j * omega(src.freq)) * s_mDeriv
def _j(self, eSolution, srcList):
"""
@@ -341,7 +341,7 @@ class Fields_e(Fields):
def _jDeriv_u(self, src, du_dm_v, adjoint = False):
"""
Derivative of the current density with respect to the thing we solved for
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray du_dm_v: vector to take product with
:param bool adjoint: adjoint?
@@ -351,15 +351,15 @@ class Fields_e(Fields):
n = int(self._aveE2CCV.shape[0] / self._nC) # number of components (instead of checking if cyl or not)
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
if adjoint:
if adjoint:
return self._eDeriv_u(src, self._MeSigma.T * (self._aveE2CCV.T * (VI.T * du_dm_v) ), adjoint = adjoint)
return VI * (self._aveE2CCV * (self._MeSigma * (self._eDeriv_u(src, du_dm_v, adjoint=adjoint) ) ) )
def _jDeriv_m(self, src, v, adjoint = False):
"""
Derivative of the current density with respect to the inversion model.
Derivative of the current density with respect to the inversion model.
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
@@ -373,7 +373,7 @@ class Fields_e(Fields):
if adjoint:
return self._MeSigmaDeriv(e).T * (self._aveE2CCV.T * (VI.T * v)) + self._eDeriv_m(src, self._aveE2CCV.T * (VI.T * v), adjoint=adjoint)
return VI * (self._aveE2CCV * ( self._eDeriv_m(src, v, adjoint=adjoint) + self._MeSigmaDeriv(e) * v))
def _h(self, eSolution, srcList):
@@ -393,7 +393,7 @@ class Fields_e(Fields):
def _hDeriv_u(self, src, du_dm_v, adjoint = False):
"""
Derivative of the magnetic field with respect to the thing we solved for
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray du_dm_v: vector to take product with
:param bool adjoint: adjoint?
@@ -409,8 +409,8 @@ class Fields_e(Fields):
def _hDeriv_m(self, src, v, adjoint = False):
"""
Derivative of the magnetic field with respect to the inversion model.
Derivative of the magnetic field with respect to the inversion model.
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
@@ -428,10 +428,10 @@ class Fields_e(Fields):
class Fields_b(Fields):
"""
Fields object for Problem_b.
Fields object for Problem_b.
:param Mesh mesh: mesh
:param Survey survey: survey
:param Survey survey: survey
"""
knownFields = {'bSolution':'F'}
@@ -506,9 +506,9 @@ class Fields_b(Fields):
def _bDeriv_u(self, src, du_dm_v, adjoint=False):
"""
Partial derivative of the total magnetic flux density with respect to the thing we
Partial derivative of the total magnetic flux density with respect to the thing we
solved for.
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray du_dm_v: vector to take product with
:param bool adjoint: adjoint?
@@ -520,8 +520,8 @@ class Fields_b(Fields):
def _bDeriv_m(self, src, v, adjoint=False):
"""
Partial derivative of the total magnetic flux density with respect to the inversion model. Here, we assume that the primary does not depend on the model. Note that this also includes derivative contributions from the sources.
Partial derivative of the total magnetic flux density with respect to the inversion model. Here, we assume that the primary does not depend on the model. Note that this also includes derivative contributions from the sources.
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
@@ -560,15 +560,15 @@ class Fields_b(Fields):
e = ( self._edgeCurl.T * ( self._MfMui * bSolution))
for i,src in enumerate(srcList):
_,S_e = src.eval(self.prob)
e[:,i] = e[:,i] + - S_e
_,s_e = src.eval(self.prob)
e[:,i] = e[:,i] + - s_e
return self._MeSigmaI * e
def _eDeriv_u(self, src, du_dm_v, adjoint=False):
"""
Derivative of the electric field with respect to the thing we solved for
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
@@ -583,8 +583,8 @@ class Fields_b(Fields):
def _eDeriv_m(self, src, v, adjoint=False):
"""
Derivative of the electric field with respect to the inversion model
Derivative of the electric field with respect to the inversion model
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
@@ -593,15 +593,15 @@ class Fields_b(Fields):
"""
bSolution = Utils.mkvc(self[src, 'bSolution'])
_,S_e = src.eval(self.prob)
_,s_e = src.eval(self.prob)
w = -S_e + self._edgeCurl.T * (self._MfMui * bSolution)
_, S_eDeriv = src.evalDeriv(self.prob, v, adjoint)
w = -s_e + self._edgeCurl.T * (self._MfMui * bSolution)
_, s_eDeriv = src.evalDeriv(self.prob, v, adjoint)
if adjoint:
return self._MeSigmaIDeriv(w).T * v - self._MeSigmaI.T * S_eDeriv
return self._MeSigmaIDeriv(w) * v - self._MeSigmaI * S_eDeriv
return self._MeSigmaIDeriv(w).T * v - self._MeSigmaI.T * s_eDeriv
return self._MeSigmaIDeriv(w) * v - self._MeSigmaI * s_eDeriv
def _j(self, bSolution, srcList):
"""
@@ -617,13 +617,13 @@ class Fields_b(Fields):
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
return VI * (self._aveE2CCV * ( self._MeSigma * self._e(bSolution,srcList ) ) )
def _jDeriv_u(self, src, du_dm_v, adjoint=False):
"""
Partial derivative of the current density with respect to the thing we
Partial derivative of the current density with respect to the thing we
solved for.
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray du_dm_v: vector to take product with
:param bool adjoint: adjoint?
@@ -639,8 +639,8 @@ class Fields_b(Fields):
def _jDeriv_m(self, src, v, adjoint=False):
"""
Derivative of the current density with respect to the inversion model
Derivative of the current density with respect to the inversion model
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
@@ -664,9 +664,9 @@ class Fields_b(Fields):
def _hDeriv_u(self, src, du_dm_v, adjoint=False):
"""
Partial derivative of the magnetic field with respect to the thing we
Partial derivative of the magnetic field with respect to the thing we
solved for.
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray du_dm_v: vector to take product with
:param bool adjoint: adjoint?
@@ -682,8 +682,8 @@ class Fields_b(Fields):
def _hDeriv_m(self, src, v, adjoint=False):
"""
Derivative of the magnetic field with respect to the inversion model
Derivative of the magnetic field with respect to the inversion model
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
@@ -695,10 +695,10 @@ class Fields_b(Fields):
class Fields_j(Fields):
"""
Fields object for Problem_j.
Fields object for Problem_j.
:param Mesh mesh: mesh
:param Survey survey: survey
:param Survey survey: survey
"""
knownFields = {'jSolution':'F'}
@@ -769,12 +769,12 @@ class Fields_j(Fields):
def _j(self, jSolution, srcList):
"""
Total current density is sum of primary and secondary
Total current density is sum of primary and secondary
:param numpy.ndarray jSolution: field we solved for
:param list srcList: list of sources
:rtype: numpy.ndarray
:return: total current density
:return: total current density
"""
return self._jPrimary(jSolution, srcList) + self._jSecondary(jSolution, srcList)
@@ -782,9 +782,9 @@ class Fields_j(Fields):
def _jDeriv_u(self, src, du_dm_v, adjoint=False):
"""
Partial derivative of the total current density with respect to the thing we
Partial derivative of the total current density with respect to the thing we
solved for.
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
@@ -797,8 +797,8 @@ class Fields_j(Fields):
def _jDeriv_m(self, src, v, adjoint=False):
"""
Partial derivative of the total current density with respect to the inversion model. Here, we assume that the primary does not depend on the model. Note that this also includes derivative contributions from the sources.
Partial derivative of the total current density with respect to the inversion model. Here, we assume that the primary does not depend on the model. Note that this also includes derivative contributions from the sources.
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
@@ -837,15 +837,15 @@ class Fields_j(Fields):
h = (self._edgeCurl.T * (self._MfRho * jSolution) )
for i, src in enumerate(srcList):
h[:,i] *= -1./(1j*omega(src.freq))
S_m,_ = src.eval(self.prob)
h[:,i] = h[:,i] + 1./(1j*omega(src.freq)) * (S_m)
s_m,_ = src.eval(self.prob)
h[:,i] = h[:,i] + 1./(1j*omega(src.freq)) * (s_m)
return self._MeMuI * h
def _hDeriv_u(self, src, du_dm_v, adjoint=False):
"""
Derivative of the magnetic field with respect to the thing we solved for
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray du_dm_v: vector to take product with
:param bool adjoint: adjoint?
@@ -856,13 +856,13 @@ class Fields_j(Fields):
if adjoint:
return -1./(1j*omega(src.freq)) * self._MfRho.T * (self._edgeCurl * ( self._MeMuI.T * du_dm_v))
return -1./(1j*omega(src.freq)) * self._MeMuI * (self._edgeCurl.T * (self._MfRho * du_dm_v) )
def _hDeriv_m(self, src, v, adjoint=False):
"""
Derivative of the magnetic field with respect to the inversion model
Derivative of the magnetic field with respect to the inversion model
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
@@ -875,19 +875,19 @@ class Fields_j(Fields):
C = self._edgeCurl
MfRho = self._MfRho
MfRhoDeriv = self._MfRhoDeriv
S_mDeriv,_ = src.evalDeriv(self.prob, adjoint = adjoint)
s_mDeriv,_ = src.evalDeriv(self.prob, adjoint = adjoint)
if not adjoint:
hDeriv_m = -1./(1j*omega(src.freq)) * MeMuI * (C.T * (MfRhoDeriv(jSolution)*v ) )
S_mDeriv = S_mDeriv(v)
hDeriv_m = hDeriv_m + 1./(1j*omega(src.freq)) * MeMuI * ( S_mDeriv)
s_mDeriv = s_mDeriv(v)
hDeriv_m = hDeriv_m + 1./(1j*omega(src.freq)) * MeMuI * ( s_mDeriv)
elif adjoint:
hDeriv_m = -1./(1j*omega(src.freq)) * MfRhoDeriv(jSolution).T * ( C * (MeMuI.T * v ) )
S_mDeriv = S_mDeriv(MeMuI.T * v)
hDeriv_m = hDeriv_m + 1./(1j*omega(src.freq)) * S_mDeriv
s_mDeriv = s_mDeriv(MeMuI.T * v)
hDeriv_m = hDeriv_m + 1./(1j*omega(src.freq)) * s_mDeriv
return hDeriv_m
def _e(self, jSolution, srcList):
@@ -901,12 +901,12 @@ class Fields_j(Fields):
"""
n = int(self._aveF2CCV.shape[0] / self._nC) # number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
return VI * (self._aveF2CCV * (self._MfRho * self._j(jSolution, srcList)))
return VI * (self._aveF2CCV * (self._MfRho * self._j(jSolution, srcList)))
def _eDeriv_u(self, src, du_dm_v, adjoint=False):
"""
Derivative of the electric field with respect to the thing we solved for
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray du_dm_v: vector to take product with
:param bool adjoint: adjoint?
@@ -921,8 +921,8 @@ class Fields_j(Fields):
def _eDeriv_m(self, src, v, adjoint=False):
"""
Derivative of the electric field with respect to the inversion model
Derivative of the electric field with respect to the inversion model
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
@@ -943,17 +943,17 @@ class Fields_j(Fields):
:param numpy.ndarray hSolution: field we solved for
:param list srcList: list of sources
:rtype: numpy.ndarray
:return: secondary magnetic flux density
:return: secondary magnetic flux density
"""
n = int(self._aveE2CCV.shape[0] / self._nC) # number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
return VI * (self._aveE2CCV * ( self._MeMu * self._h(jSolution,srcList)) )
return VI * (self._aveE2CCV * ( self._MeMu * self._h(jSolution,srcList)) )
def _bDeriv_u(self, src, du_dm_v, adjoint=False):
"""
Derivative of the magnetic flux density with respect to the thing we solved for
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray du_dm_v: vector to take product with
:param bool adjoint: adjoint?
@@ -969,8 +969,8 @@ class Fields_j(Fields):
def _bDeriv_m(self, src, v, adjoint=False):
"""
Derivative of the magnetic flux density with respect to the inversion model
Derivative of the magnetic flux density with respect to the inversion model
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
@@ -980,20 +980,20 @@ class Fields_j(Fields):
jSolution = self[src,'jSolution']
n = int(self._aveE2CCV.shape[0] / self._nC) # number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
S_mDeriv,_ = src.evalDeriv(self.prob, adjoint = adjoint)
s_mDeriv,_ = src.evalDeriv(self.prob, adjoint = adjoint)
if adjoint:
v = self._aveE2CCV.T * ( VI.T * v)
return 1./(1j * omega(src.freq)) * ( S_mDeriv(v) - self._MfRhoDeriv(jSolution).T * (self._edgeCurl * v ))
return 1./(1j * omega(src.freq)) * VI * (self._aveE2CCV * ( S_mDeriv(v) - self._edgeCurl.T * ( self._MfRhoDeriv(jSolution) * v ) ) )
return 1./(1j * omega(src.freq)) * ( s_mDeriv(v) - self._MfRhoDeriv(jSolution).T * (self._edgeCurl * v ))
return 1./(1j * omega(src.freq)) * VI * (self._aveE2CCV * ( s_mDeriv(v) - self._edgeCurl.T * ( self._MfRhoDeriv(jSolution) * v ) ) )
class Fields_h(Fields):
"""
Fields object for Problem_h.
Fields object for Problem_h.
:param Mesh mesh: mesh
:param Survey survey: survey
:param Survey survey: survey
"""
knownFields = {'hSolution':'E'}
@@ -1065,9 +1065,9 @@ class Fields_h(Fields):
def _hDeriv_u(self, src, du_dm_v, adjoint=False):
"""
Partial derivative of the total magnetic field with respect to the thing we
Partial derivative of the total magnetic field with respect to the thing we
solved for.
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray du_dm_v: vector to take product with
:param bool adjoint: adjoint?
@@ -1079,8 +1079,8 @@ class Fields_h(Fields):
def _hDeriv_m(self, src, v, adjoint=False):
"""
Partial derivative of the total magnetic field with respect to the inversion model. Here, we assume that the primary does not depend on the model. Note that this also includes derivative contributions from the sources.
Partial derivative of the total magnetic field with respect to the inversion model. Here, we assume that the primary does not depend on the model. Note that this also includes derivative contributions from the sources.
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
@@ -1119,14 +1119,14 @@ class Fields_h(Fields):
j = self._edgeCurl*hSolution
for i, src in enumerate(srcList):
_,S_e = src.eval(self.prob)
j[:,i] = j[:,i]+ -S_e
_,s_e = src.eval(self.prob)
j[:,i] = j[:,i]+ -s_e
return j
def _jDeriv_u(self, src, du_dm_v, adjoint=False):
"""
Derivative of the current density with respect to the thing we solved for
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray du_dm_v: vector to take product with
:param bool adjoint: adjoint?
@@ -1142,8 +1142,8 @@ class Fields_h(Fields):
def _jDeriv_m(self, src, v, adjoint=False):
"""
Derivative of the current density with respect to the inversion model.
Derivative of the current density with respect to the inversion model.
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
@@ -1151,9 +1151,9 @@ class Fields_h(Fields):
:return: product of the current density derivative with respect to the inversion model with a vector
"""
_,S_eDeriv = src.evalDeriv(self.prob, v, adjoint)
return -S_eDeriv
_,s_eDeriv = src.evalDeriv(self.prob, v, adjoint)
return -s_eDeriv
def _e(self, hSolution, srcList):
"""
Electric field from hSolution
@@ -1165,12 +1165,12 @@ class Fields_h(Fields):
"""
n = int(self._aveF2CCV.shape[0] / self._nC) #number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
return VI * (self._aveF2CCV * (self._MfRho * self._j(hSolution, srcList)))
return VI * (self._aveF2CCV * (self._MfRho * self._j(hSolution, srcList)))
def _eDeriv_u(self, src, du_dm_v, adjoint=False):
"""
Derivative of the electric field with respect to the thing we solved for
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray du_dm_v: vector to take product with
:param bool adjoint: adjoint?
@@ -1181,12 +1181,12 @@ class Fields_h(Fields):
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
if adjoint:
return self._edgeCurl.T * ( self._MfRho.T * ( self._aveF2CCV.T * ( VI.T * du_dm_v ) ) )
return VI * (self._aveF2CCV * (self._MfRho * self._edgeCurl * du_dm_v ))
return VI * (self._aveF2CCV * (self._MfRho * self._edgeCurl * du_dm_v ))
def _eDeriv_m(self, src, v, adjoint=False):
"""
Derivative of the electric field with respect to the inversion model.
Derivative of the electric field with respect to the inversion model.
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
@@ -1196,7 +1196,7 @@ class Fields_h(Fields):
hSolution = Utils.mkvc(self[src,'hSolution'])
n = int(self._aveF2CCV.shape[0] / self._nC) #number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
if adjoint:
if adjoint:
return ( self._MfRhoDeriv(self._edgeCurl * hSolution).T * ( self._aveF2CCV.T * (VI.T * v) ) )
return VI * (self._aveF2CCV * (self._MfRhoDeriv(self._edgeCurl * hSolution) * v ))
@@ -1207,10 +1207,10 @@ class Fields_h(Fields):
:param numpy.ndarray hSolution: field we solved for
:param list srcList: list of sources
:rtype: numpy.ndarray
:return: magnetic flux density
:return: magnetic flux density
"""
h = self._h(hSolution, srcList)
n = int(self._aveE2CCV.shape[0] / self._nC) #number of components
n = int(self._aveE2CCV.shape[0] / self._nC) #number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
return VI * (self._aveE2CCV * (self._MeMu * h))
@@ -1218,14 +1218,14 @@ class Fields_h(Fields):
def _bDeriv_u(self, src, du_dm_v, adjoint=False):
"""
Derivative of the magnetic flux density with respect to the thing we solved for
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray du_dm_v: vector to take product with
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
:return: product of the derivative of the magnetic flux density with respect to the field we solved for with a vector
"""
n = int(self._aveE2CCV.shape[0] / self._nC) #number of components
n = int(self._aveE2CCV.shape[0] / self._nC) #number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
if adjoint:
return self._MeMu.T * (self._aveE2CCV.T * ( VI.T * du_dm_v ))
@@ -1233,8 +1233,8 @@ class Fields_h(Fields):
def _bDeriv_m(self, src, v, adjoint=False):
"""
Derivative of the magnetic flux density with respect to the inversion model.
Derivative of the magnetic flux density with respect to the inversion model.
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
SimPEG/EM/FDEM/SrcFDEM.py
+48 -48
View File
@@ -15,22 +15,22 @@ class BaseSrc(Survey.BaseSrc):
def eval(self, prob):
"""
Evaluate the source terms.
- :math:`S_m` : magnetic source term
- :math:`S_e` : electric source term
- :math:`s_m` : magnetic source term
- :math:`s_e` : electric source term
:param Problem prob: FDEM Problem
:rtype: (numpy.ndarray, numpy.ndarray)
:return: tuple with magnetic source term and electric source term
"""
S_m = self.S_m(prob)
S_e = self.S_e(prob)
return S_m, S_e
s_m = self.s_m(prob)
s_e = self.s_e(prob)
return s_m, s_e
def evalDeriv(self, prob, v=None, adjoint=False):
"""
Derivatives of the source terms with respect to the inversion model
- :code:`S_mDeriv` : derivative of the magnetic source term
- :code:`S_eDeriv` : derivative of the electric source term
- :code:`s_mDeriv` : derivative of the magnetic source term
- :code:`s_eDeriv` : derivative of the electric source term
:param Problem prob: FDEM Problem
:param numpy.ndarray v: vector to take product with
@@ -39,9 +39,9 @@ class BaseSrc(Survey.BaseSrc):
:return: tuple with magnetic source term and electric source term derivatives times a vector
"""
if v is not None:
return self.S_mDeriv(prob, v, adjoint), self.S_eDeriv(prob, v, adjoint)
return self.s_mDeriv(prob, v, adjoint), self.s_eDeriv(prob, v, adjoint)
else:
return lambda v: self.S_mDeriv(prob, v, adjoint), lambda v: self.S_eDeriv(prob, v, adjoint)
return lambda v: self.s_mDeriv(prob, v, adjoint), lambda v: self.s_eDeriv(prob, v, adjoint)
def bPrimary(self, prob):
"""
@@ -83,7 +83,7 @@ class BaseSrc(Survey.BaseSrc):
"""
return Zero()
def S_m(self, prob):
def s_m(self, prob):
"""
Magnetic source term
@@ -93,7 +93,7 @@ class BaseSrc(Survey.BaseSrc):
"""
return Zero()
def S_e(self, prob):
def s_e(self, prob):
"""
Electric source term
@@ -103,7 +103,7 @@ class BaseSrc(Survey.BaseSrc):
"""
return Zero()
def S_mDeriv(self, prob, v, adjoint = False):
def s_mDeriv(self, prob, v, adjoint = False):
"""
Derivative of magnetic source term with respect to the inversion model
@@ -116,7 +116,7 @@ class BaseSrc(Survey.BaseSrc):
return Zero()
def S_eDeriv(self, prob, v, adjoint = False):
def s_eDeriv(self, prob, v, adjoint = False):
"""
Derivative of electric source term with respect to the inversion model
@@ -131,22 +131,22 @@ class BaseSrc(Survey.BaseSrc):
class RawVec_e(BaseSrc):
"""
RawVec electric source. It is defined by the user provided vector S_e
RawVec electric source. It is defined by the user provided vector s_e
:param list rxList: receiver list
:param float freq: frequency
:param numpy.array S_e: electric source term
:param numpy.array s_e: electric source term
:param bool integrate: Integrate the source term (multiply by Me) [True]
"""
def __init__(self, rxList, freq, S_e, integrate=True): #, ePrimary=None, bPrimary=None, hPrimary=None, jPrimary=None):
self._S_e = np.array(S_e, dtype=complex)
def __init__(self, rxList, freq, s_e, integrate=True): #, ePrimary=None, bPrimary=None, hPrimary=None, jPrimary=None):
self._s_e = np.array(s_e, dtype=complex)
self.freq = float(freq)
self.integrate = integrate
BaseSrc.__init__(self, rxList)
def S_e(self, prob):
def s_e(self, prob):
"""
Electric source term
@@ -155,28 +155,28 @@ class RawVec_e(BaseSrc):
:return: electric source term on mesh
"""
if prob._formulation is 'EB' and self.integrate is True:
return prob.Me * self._S_e
return self._S_e
return prob.Me * self._s_e
return self._s_e
class RawVec_m(BaseSrc):
"""
RawVec magnetic source. It is defined by the user provided vector S_m
RawVec magnetic source. It is defined by the user provided vector s_m
:param float freq: frequency
:param rxList: receiver list
:param numpy.array S_m: magnetic source term
:param numpy.array s_m: magnetic source term
:param bool integrate: Integrate the source term (multiply by Me) [True]
"""
def __init__(self, rxList, freq, S_m, integrate=True): #ePrimary=Zero(), bPrimary=Zero(), hPrimary=Zero(), jPrimary=Zero()):
self._S_m = np.array(S_m, dtype=complex)
def __init__(self, rxList, freq, s_m, integrate=True): #ePrimary=Zero(), bPrimary=Zero(), hPrimary=Zero(), jPrimary=Zero()):
self._s_m = np.array(s_m, dtype=complex)
self.freq = float(freq)
self.integrate = integrate
BaseSrc.__init__(self, rxList)
def S_m(self, prob):
def s_m(self, prob):
"""
Magnetic source term
@@ -185,28 +185,28 @@ class RawVec_m(BaseSrc):
:return: magnetic source term on mesh
"""
if prob._formulation is 'HJ' and self.integrate is True:
return prob.Me * self._S_m
return self._S_m
return prob.Me * self._s_m
return self._s_m
class RawVec(BaseSrc):
"""
RawVec source. It is defined by the user provided vectors S_m, S_e
RawVec source. It is defined by the user provided vectors s_m, s_e
:param rxList: receiver list
:param float freq: frequency
:param numpy.array S_m: magnetic source term
:param numpy.array S_e: electric source term
:param numpy.array s_m: magnetic source term
:param numpy.array s_e: electric source term
:param bool integrate: Integrate the source term (multiply by Me) [True]
"""
def __init__(self, rxList, freq, S_m, S_e, integrate=True):
self._S_m = np.array(S_m, dtype=complex)
self._S_e = np.array(S_e, dtype=complex)
def __init__(self, rxList, freq, s_m, s_e, integrate=True):
self._s_m = np.array(s_m, dtype=complex)
self._s_e = np.array(s_e, dtype=complex)
self.freq = float(freq)
self.integrate = integrate
BaseSrc.__init__(self, rxList)
def S_m(self, prob):
def s_m(self, prob):
"""
Magnetic source term
@@ -215,10 +215,10 @@ class RawVec(BaseSrc):
:return: magnetic source term on mesh
"""
if prob._formulation is 'HJ' and self.integrate is True:
return prob.Me * self._S_m
return self._S_m
return prob.Me * self._s_m
return self._s_m
def S_e(self, prob):
def s_e(self, prob):
"""
Electric source term
@@ -227,8 +227,8 @@ class RawVec(BaseSrc):
:return: electric source term on mesh
"""
if prob._formulation is 'EB' and self.integrate is True:
return prob.Me * self._S_e
return self._S_e
return prob.Me * self._s_e
return self._s_e
class MagDipole(BaseSrc):
@@ -335,9 +335,9 @@ class MagDipole(BaseSrc):
:return: primary magnetic field
"""
b = self.bPrimary(prob)
return 1./self.mu * b
return 1./self.mu * b
def S_m(self, prob):
def s_m(self, prob):
"""
The magnetic source term
@@ -348,10 +348,10 @@ class MagDipole(BaseSrc):
b_p = self.bPrimary(prob)
if prob._formulation is 'HJ':
b_p = prob.Me * b_p
b_p = prob.Me * b_p
return -1j*omega(self.freq)*b_p
def S_e(self, prob):
def s_e(self, prob):
"""
The electric source term
@@ -453,7 +453,7 @@ class MagDipole_Bfield(BaseSrc):
b = self.bPrimary(prob)
return 1/self.mu * b
def S_m(self, prob):
def s_m(self, prob):
"""
The magnetic source term
@@ -466,7 +466,7 @@ class MagDipole_Bfield(BaseSrc):
b = prob.Me * b
return -1j*omega(self.freq)*b
def S_e(self, prob):
def s_e(self, prob):
"""
The electric source term
@@ -565,7 +565,7 @@ class CircularLoop(BaseSrc):
b = self.bPrimary(prob)
return 1./self.mu*b
def S_m(self, prob):
def s_m(self, prob):
"""
The magnetic source term
@@ -578,7 +578,7 @@ class CircularLoop(BaseSrc):
b = prob.Me * b
return -1j*omega(self.freq)*b
def S_e(self, prob):
def s_e(self, prob):
"""
The electric source term
@@ -604,6 +604,6 @@ class CircularLoop(BaseSrc):
return -C.T * (MMui_s * self.bPrimary(prob))