mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-02 07:56:09 +08:00
Merge pull request #244 from simpeg/em/FDEMfieldsDerivs
Em/fdem fields derivs
This commit is contained in:
Lindsey
+31
-42
@@ -84,23 +84,15 @@ class BaseFDEMProblem(BaseEMProblem):
|
||||
for src in self.survey.getSrcByFreq(freq):
|
||||
ftype = self._fieldType + 'Solution'
|
||||
u_src = u[src, ftype]
|
||||
dA_dm = self.getADeriv_m(freq, u_src, v)
|
||||
dRHS_dm = self.getRHSDeriv_m(freq, src, v)
|
||||
du_dm = Ainv * ( - dA_dm + dRHS_dm )
|
||||
dA_dm_v = self.getADeriv(freq, u_src, v)
|
||||
dRHS_dm_v = self.getRHSDeriv(freq, src, v)
|
||||
du_dm_v = Ainv * ( - dA_dm_v + dRHS_dm_v )
|
||||
|
||||
for rx in src.rxList:
|
||||
df_duFun = getattr(u, '_%sDeriv_u'%rx.projField, None)
|
||||
df_dudu_dm = df_duFun(src, du_dm, adjoint=False)
|
||||
|
||||
df_dmFun = getattr(u, '_%sDeriv_m'%rx.projField, None)
|
||||
df_dm = df_dmFun(src, v, adjoint=False)
|
||||
|
||||
|
||||
Df_Dm = np.array(df_dudu_dm + df_dm,dtype=complex)
|
||||
|
||||
P = lambda v: rx.projectFieldsDeriv(src, self.mesh, u, v) # wrt u, also have wrt m
|
||||
|
||||
Jv[src, rx] = P(Df_Dm)
|
||||
df_dmFun = getattr(u, '_%sDeriv'%rx.projField, None)
|
||||
df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False)
|
||||
df_dm_v = np.array(df_dm_v,dtype=complex)
|
||||
Jv[src, rx] = rx.projectFieldsDeriv(src, self.mesh, u, df_dm_v)
|
||||
|
||||
Ainv.clean()
|
||||
return Utils.mkvc(Jv)
|
||||
@@ -138,26 +130,23 @@ class BaseFDEMProblem(BaseEMProblem):
|
||||
for rx in src.rxList:
|
||||
PTv = rx.projectFieldsDeriv(src, self.mesh, u, v[src, rx], adjoint=True) # wrt u, need possibility wrt m
|
||||
|
||||
df_duTFun = getattr(u, '_%sDeriv_u'%rx.projField, None)
|
||||
df_duT = df_duTFun(src, PTv, adjoint=True)
|
||||
|
||||
df_duTFun = getattr(u, '_%sDeriv'%rx.projField, None)
|
||||
df_duT, df_dmT = df_duTFun(src, None, PTv, adjoint=True)
|
||||
|
||||
ATinvdf_duT = ATinv * df_duT
|
||||
|
||||
dA_dmT = self.getADeriv_m(freq, u_src, ATinvdf_duT, adjoint=True)
|
||||
dRHS_dmT = self.getRHSDeriv_m(freq,src, ATinvdf_duT, adjoint=True)
|
||||
dA_dmT = self.getADeriv(freq, u_src, ATinvdf_duT, adjoint=True)
|
||||
dRHS_dmT = self.getRHSDeriv(freq, src, ATinvdf_duT, adjoint=True)
|
||||
du_dmT = -dA_dmT + dRHS_dmT
|
||||
|
||||
df_dmFun = getattr(u, '_%sDeriv_m'%rx.projField, None)
|
||||
dfT_dm = df_dmFun(src, PTv, adjoint=True)
|
||||
df_dmT += du_dmT
|
||||
|
||||
du_dmT += dfT_dm
|
||||
|
||||
# TODO: this should be taken care of by the reciever
|
||||
# TODO: this should be taken care of by the reciever?
|
||||
real_or_imag = rx.projComp
|
||||
if real_or_imag is 'real':
|
||||
Jtv += np.array(du_dmT,dtype=complex).real
|
||||
Jtv += np.array(df_dmT,dtype=complex).real
|
||||
elif real_or_imag is 'imag':
|
||||
Jtv += - np.array(du_dmT,dtype=complex).real
|
||||
Jtv += - np.array(df_dmT,dtype=complex).real
|
||||
else:
|
||||
raise Exception('Must be real or imag')
|
||||
|
||||
@@ -239,7 +228,7 @@ class Problem_e(BaseFDEMProblem):
|
||||
return C.T*MfMui*C + 1j*omega(freq)*MeSigma
|
||||
|
||||
|
||||
def getADeriv_m(self, freq, u, v, adjoint=False):
|
||||
def getADeriv(self, freq, u, v, adjoint=False):
|
||||
"""
|
||||
Product of the derivative of our system matrix with respect to the model and a vector
|
||||
|
||||
@@ -280,7 +269,7 @@ class Problem_e(BaseFDEMProblem):
|
||||
|
||||
return C.T * (MfMui * S_m) -1j * omega(freq) * S_e
|
||||
|
||||
def getRHSDeriv_m(self, freq, src, v, adjoint=False):
|
||||
def getRHSDeriv(self, freq, src, v, adjoint=False):
|
||||
"""
|
||||
Derivative of the right hand side with respect to the model
|
||||
|
||||
@@ -354,7 +343,7 @@ class Problem_b(BaseFDEMProblem):
|
||||
return MfMui.T*A
|
||||
return A
|
||||
|
||||
def getADeriv_m(self, freq, u, v, adjoint=False):
|
||||
def getADeriv(self, freq, u, v, adjoint=False):
|
||||
|
||||
"""
|
||||
Product of the derivative of our system matrix with respect to the model and a vector
|
||||
@@ -411,7 +400,7 @@ class Problem_b(BaseFDEMProblem):
|
||||
|
||||
return RHS
|
||||
|
||||
def getRHSDeriv_m(self, freq, src, v, adjoint=False):
|
||||
def getRHSDeriv(self, freq, src, v, adjoint=False):
|
||||
"""
|
||||
Derivative of the right hand side with respect to the model
|
||||
|
||||
@@ -503,7 +492,7 @@ class Problem_j(BaseFDEMProblem):
|
||||
return A
|
||||
|
||||
|
||||
def getADeriv_m(self, freq, u, v, adjoint=False):
|
||||
def getADeriv(self, freq, u, v, adjoint=False):
|
||||
"""
|
||||
Product of the derivative of our system matrix with respect to the model and a vector
|
||||
|
||||
@@ -524,16 +513,16 @@ class Problem_j(BaseFDEMProblem):
|
||||
MeMuI = self.MeMuI
|
||||
MfRho = self.MfRho
|
||||
C = self.mesh.edgeCurl
|
||||
MfRhoDeriv_m = self.MfRhoDeriv(u)
|
||||
MfRhoDeriv = self.MfRhoDeriv(u)
|
||||
|
||||
if adjoint:
|
||||
if self._makeASymmetric is True:
|
||||
v = MfRho * v
|
||||
return MfRhoDeriv_m.T * (C * (MeMuI.T * (C.T * v)))
|
||||
return MfRhoDeriv.T * (C * (MeMuI.T * (C.T * v)))
|
||||
|
||||
if self._makeASymmetric is True:
|
||||
return MfRho.T * (C * ( MeMuI * (C.T * (MfRhoDeriv_m * v) )))
|
||||
return C * (MeMuI * (C.T * (MfRhoDeriv_m * v)))
|
||||
return MfRho.T * (C * ( MeMuI * (C.T * (MfRhoDeriv * v) )))
|
||||
return C * (MeMuI * (C.T * (MfRhoDeriv * v)))
|
||||
|
||||
|
||||
def getRHS(self, freq):
|
||||
@@ -560,7 +549,7 @@ class Problem_j(BaseFDEMProblem):
|
||||
|
||||
return RHS
|
||||
|
||||
def getRHSDeriv_m(self, freq, src, v, adjoint=False):
|
||||
def getRHSDeriv(self, freq, src, v, adjoint=False):
|
||||
"""
|
||||
Derivative of the right hand side with respect to the model
|
||||
|
||||
@@ -635,7 +624,7 @@ class Problem_h(BaseFDEMProblem):
|
||||
|
||||
return C.T * (MfRho * C) + 1j*omega(freq)*MeMu
|
||||
|
||||
def getADeriv_m(self, freq, u, v, adjoint=False):
|
||||
def getADeriv(self, freq, u, v, adjoint=False):
|
||||
"""
|
||||
Product of the derivative of our system matrix with respect to the model and a vector
|
||||
|
||||
@@ -652,11 +641,11 @@ class Problem_h(BaseFDEMProblem):
|
||||
|
||||
MeMu = self.MeMu
|
||||
C = self.mesh.edgeCurl
|
||||
MfRhoDeriv_m = self.MfRhoDeriv(C*u)
|
||||
MfRhoDeriv = self.MfRhoDeriv(C*u)
|
||||
|
||||
if adjoint:
|
||||
return MfRhoDeriv_m.T * (C * v)
|
||||
return C.T * (MfRhoDeriv_m * v)
|
||||
return MfRhoDeriv.T * (C * v)
|
||||
return C.T * (MfRhoDeriv * v)
|
||||
|
||||
def getRHS(self, freq):
|
||||
"""
|
||||
@@ -677,7 +666,7 @@ class Problem_h(BaseFDEMProblem):
|
||||
|
||||
return S_m + C.T * ( MfRho * S_e )
|
||||
|
||||
def getRHSDeriv_m(self, freq, src, v, adjoint=False):
|
||||
def getRHSDeriv(self, freq, src, v, adjoint=False):
|
||||
"""
|
||||
Derivative of the right hand side with respect to the model
|
||||
|
||||
|
||||
+189
-136
@@ -32,6 +32,135 @@ class Fields(SimPEG.Problem.Fields):
|
||||
knownFields = {}
|
||||
dtype = complex
|
||||
|
||||
def _e(self, solution, srcList):
|
||||
"""
|
||||
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:
|
||||
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
|
||||
|
||||
:param numpy.ndarray solution: field we solved for
|
||||
:param list srcList: list of sources
|
||||
:rtype: numpy.ndarray
|
||||
:return: total magnetic flux density
|
||||
"""
|
||||
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
|
||||
|
||||
: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:
|
||||
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
|
||||
|
||||
:param numpy.ndarray solution: field we solved for
|
||||
:param list srcList: list of sources
|
||||
:rtype: numpy.ndarray
|
||||
:return: total current density
|
||||
"""
|
||||
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)
|
||||
|
||||
def _eDeriv(self, src, du_dm_v, v, adjoint = False):
|
||||
"""
|
||||
Total derivative of e with respect to the inversion model. Returns :math:`d\mathbf{e}/d\mathbf{m}` for forward and (:math:`d\mathbf{e}/d\mathbf{u}`, :math:`d\mathb{u}/d\mathbf{m}`) for the adjoint
|
||||
|
||||
:param Src src: sorce
|
||||
:param numpy.ndarray du_dm_v: derivative of the solution vector with respect to the model times a vector (is None for adjoint)
|
||||
:param numpy.ndarray v: vector to take sensitivity product with
|
||||
:param bool adjoint: adjoint?
|
||||
: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:
|
||||
raise NotImplementedError ('Getting eDerivs from %s is not implemented' %self.knownFields.keys()[0])
|
||||
|
||||
if adjoint:
|
||||
return self._eDeriv_u(src, v, adjoint), self._eDeriv_m(src, v, adjoint)
|
||||
return self._eDeriv_u(src, du_dm_v, adjoint) + self._eDeriv_m(src, v, adjoint)
|
||||
|
||||
def _bDeriv(self, src, du_dm_v, v, adjoint = False):
|
||||
"""
|
||||
Total derivative of b with respect to the inversion model. Returns :math:`d\mathbf{b}/d\mathbf{m}` for forward and (:math:`d\mathbf{b}/d\mathbf{u}`, :math:`d\mathb{u}/d\mathbf{m}`) for the adjoint
|
||||
|
||||
:param Src src: sorce
|
||||
:param numpy.ndarray du_dm_v: derivative of the solution vector with respect to the model times a vector (is None for adjoint)
|
||||
:param numpy.ndarray v: vector to take sensitivity product with
|
||||
:param bool adjoint: adjoint?
|
||||
: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:
|
||||
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 self._bDeriv_u(src, du_dm_v, adjoint) + self._bDeriv_m(src, v, adjoint)
|
||||
|
||||
def _hDeriv(self, src, du_dm_v, v, adjoint = False):
|
||||
"""
|
||||
Total derivative of h with respect to the inversion model. Returns :math:`d\mathbf{h}/d\mathbf{m}` for forward and (:math:`d\mathbf{h}/d\mathbf{u}`, :math:`d\mathb{u}/d\mathbf{m}`) for the adjoint
|
||||
|
||||
:param Src src: sorce
|
||||
:param numpy.ndarray du_dm_v: derivative of the solution vector with respect to the model times a vector (is None for adjoint)
|
||||
:param numpy.ndarray v: vector to take sensitivity product with
|
||||
:param bool adjoint: adjoint?
|
||||
: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:
|
||||
raise NotImplementedError ('Getting hDerivs from %s is not implemented' %self.knownFields.keys()[0])
|
||||
|
||||
if adjoint:
|
||||
return self._hDeriv_u(src, v, adjoint), self._hDeriv_m(src, v, adjoint)
|
||||
return self._hDeriv_u(src, du_dm_v, adjoint) + self._hDeriv_m(src, v, adjoint)
|
||||
|
||||
def _jDeriv(self, src, du_dm_v, v, adjoint = False):
|
||||
"""
|
||||
Total derivative of j with respect to the inversion model. Returns :math:`d\mathbf{j}/d\mathbf{m}` for forward and (:math:`d\mathbf{j}/d\mathbf{u}`, :math:`d\mathb{u}/d\mathbf{m}`) for the adjoint
|
||||
|
||||
:param Src src: sorce
|
||||
:param numpy.ndarray du_dm_v: derivative of the solution vector with respect to the model times a vector (is None for adjoint)
|
||||
:param numpy.ndarray v: vector to take sensitivity product with
|
||||
:param bool adjoint: adjoint?
|
||||
: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:
|
||||
raise NotImplementedError ('Getting jDerivs from %s is not implemented' %self.knownFields.keys()[0])
|
||||
|
||||
if adjoint:
|
||||
return self._jDeriv_u(src, v, adjoint), self._jDeriv_m(src, v, adjoint)
|
||||
return self._jDeriv_u(src, du_dm_v, adjoint) + self._jDeriv_m(src, v, adjoint)
|
||||
|
||||
|
||||
class Fields_e(Fields):
|
||||
"""
|
||||
Fields object for Problem_e.
|
||||
@@ -85,22 +214,11 @@ class Fields_e(Fields):
|
||||
|
||||
return eSolution
|
||||
|
||||
def _e(self, eSolution, srcList):
|
||||
"""
|
||||
Total electric field is sum of primary and secondary
|
||||
|
||||
:param numpy.ndarray eSolution: field we solved for
|
||||
:param list srcList: list of sources
|
||||
:rtype: numpy.ndarray
|
||||
:return: total electric field
|
||||
"""
|
||||
|
||||
return self._ePrimary(eSolution,srcList) + self._eSecondary(eSolution,srcList)
|
||||
|
||||
def _eDeriv_u(self, src, v, adjoint = False):
|
||||
"""
|
||||
Derivative of the total electric field with respect to the thing we
|
||||
solved for
|
||||
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
|
||||
@@ -113,7 +231,7 @@ class Fields_e(Fields):
|
||||
|
||||
def _eDeriv_m(self, src, v, adjoint = False):
|
||||
"""
|
||||
Derivative of the total electric field with respect to the inversion model. Here, we assume that the primary does not depend on the model.
|
||||
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
|
||||
@@ -159,12 +277,12 @@ class Fields_e(Fields):
|
||||
b[:,i] = b[:,i]+ 1./(1j*omega(src.freq)) * S_m
|
||||
return b
|
||||
|
||||
def _bSecondaryDeriv_u(self, src, v, adjoint = False):
|
||||
def _bSecondaryDeriv_u(self, src, du_dm_v, adjoint = False):
|
||||
"""
|
||||
Derivative of the secondary magnetic flux 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 numpy.ndarray du_dm_v: vector to take product with
|
||||
:param bool adjoint: adjoint?
|
||||
:rtype: numpy.ndarray
|
||||
:return: product of the derivative of the secondary magnetic flux density with respect to the field we solved for with a vector
|
||||
@@ -172,8 +290,8 @@ class Fields_e(Fields):
|
||||
|
||||
C = self._edgeCurl
|
||||
if adjoint:
|
||||
return - 1./(1j*omega(src.freq)) * (C.T * v)
|
||||
return - 1./(1j*omega(src.freq)) * (C * v)
|
||||
return - 1./(1j*omega(src.freq)) * (C.T * du_dm_v)
|
||||
return - 1./(1j*omega(src.freq)) * (C * du_dm_v)
|
||||
|
||||
def _bSecondaryDeriv_m(self, src, v, adjoint = False):
|
||||
"""
|
||||
@@ -189,35 +307,23 @@ class Fields_e(Fields):
|
||||
S_mDeriv, _ = src.evalDeriv(self.prob, v, adjoint)
|
||||
return 1./(1j * omega(src.freq)) * S_mDeriv
|
||||
|
||||
def _b(self, eSolution, srcList):
|
||||
"""
|
||||
Total magnetic flux density is sum of primary and secondary
|
||||
|
||||
:param numpy.ndarray eSolution: field we solved for
|
||||
:param list srcList: list of sources
|
||||
:rtype: numpy.ndarray
|
||||
:return: total magnetic flux density
|
||||
"""
|
||||
|
||||
return self._bPrimary(eSolution, srcList) + self._bSecondary(eSolution, srcList)
|
||||
|
||||
def _bDeriv_u(self, src, v, adjoint=False):
|
||||
def _bDeriv_u(self, src, du_dm_v, adjoint=False):
|
||||
"""
|
||||
Derivative of the total magnetic flux density with respect to the thing we solved for
|
||||
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 v: vector to take product with
|
||||
: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
|
||||
"""
|
||||
|
||||
# Primary does not depend on u
|
||||
return self._bSecondaryDeriv_u(src, v, adjoint)
|
||||
return self._bSecondaryDeriv_u(src, du_dm_v, adjoint)
|
||||
|
||||
def _bDeriv_m(self, src, v, adjoint=False):
|
||||
"""
|
||||
Derivative of the total magnetic flux density with respect to the inversion model.
|
||||
Partial derivative of the total 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
|
||||
@@ -287,34 +393,23 @@ class Fields_b(Fields):
|
||||
|
||||
return bSolution
|
||||
|
||||
def _b(self, bSolution, srcList):
|
||||
def _bDeriv_u(self, src, du_dm_v, adjoint=False):
|
||||
"""
|
||||
Total magnetic flux density is sum of primary and secondary
|
||||
|
||||
:param numpy.ndarray bSolution: field we solved for
|
||||
:param list srcList: list of sources
|
||||
:rtype: numpy.ndarray
|
||||
:return: total magnetic flux density
|
||||
"""
|
||||
|
||||
return self._bPrimary(bSolution, srcList) + self._bSecondary(bSolution, srcList)
|
||||
|
||||
def _bDeriv_u(self, src, v, adjoint=False):
|
||||
"""
|
||||
Derivative of the total magnetic flux density with respect to the thing we
|
||||
solved for
|
||||
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 v: vector to take product with
|
||||
: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
|
||||
"""
|
||||
return Identity()*v
|
||||
|
||||
return Identity()*du_dm_v
|
||||
|
||||
def _bDeriv_m(self, src, v, adjoint=False):
|
||||
"""
|
||||
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.
|
||||
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
|
||||
@@ -358,7 +453,7 @@ class Fields_b(Fields):
|
||||
e[:,i] = e[:,i]+ -self._MeSigmaI * S_e
|
||||
return e
|
||||
|
||||
def _eSecondaryDeriv_u(self, src, v, adjoint=False):
|
||||
def _eSecondaryDeriv_u(self, src, du_dm_v, adjoint=False):
|
||||
"""
|
||||
Derivative of the secondary electric field with respect to the thing we solved for
|
||||
|
||||
@@ -370,9 +465,9 @@ class Fields_b(Fields):
|
||||
"""
|
||||
|
||||
if not adjoint:
|
||||
return self._MeSigmaI * ( self._edgeCurl.T * ( self._MfMui * v) )
|
||||
return self._MeSigmaI * ( self._edgeCurl.T * ( self._MfMui * du_dm_v) )
|
||||
else:
|
||||
return self._MfMui.T * (self._edgeCurl * (self._MeSigmaI.T * v))
|
||||
return self._MfMui.T * (self._edgeCurl * (self._MeSigmaI.T * du_dm_v))
|
||||
|
||||
def _eSecondaryDeriv_m(self, src, v, adjoint=False):
|
||||
"""
|
||||
@@ -406,34 +501,22 @@ class Fields_b(Fields):
|
||||
|
||||
return de_dm
|
||||
|
||||
def _e(self, bSolution, srcList):
|
||||
def _eDeriv_u(self, src, du_dm_v, adjoint=False):
|
||||
"""
|
||||
Total electric field is sum of primary and secondary
|
||||
|
||||
:param numpy.ndarray eSolution: field we solved for
|
||||
:param list srcList: list of sources
|
||||
:rtype: numpy.ndarray
|
||||
:return: total electric field
|
||||
"""
|
||||
|
||||
return self._ePrimary(bSolution, srcList) + self._eSecondary(bSolution, srcList)
|
||||
|
||||
def _eDeriv_u(self, src, v, adjoint=False):
|
||||
"""
|
||||
Derivative of the total electric field with respect to the thing we solved for
|
||||
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 numpy.ndarray du_dm_v: vector to take product with
|
||||
:param bool adjoint: adjoint?
|
||||
:rtype: numpy.ndarray
|
||||
:return: product of the derivative of the electric field with respect to the field we solved for with a vector
|
||||
"""
|
||||
|
||||
return self._eSecondaryDeriv_u(src, v, adjoint)
|
||||
return self._eSecondaryDeriv_u(src, du_dm_v, adjoint)
|
||||
|
||||
def _eDeriv_m(self, src, v, adjoint=False):
|
||||
"""
|
||||
Derivative of the total electric field density with respect to the inversion model.
|
||||
Partial derivative of the total electric field density with respect to the inversion model.
|
||||
|
||||
:param SimPEG.EM.FDEM.Src src: source
|
||||
:param numpy.ndarray v: vector to take product with
|
||||
@@ -515,10 +598,11 @@ class Fields_j(Fields):
|
||||
|
||||
return self._jPrimary(jSolution, srcList) + self._jSecondary(jSolution, srcList)
|
||||
|
||||
def _jDeriv_u(self, src, v, adjoint=False):
|
||||
|
||||
def _jDeriv_u(self, src, du_dm_v, adjoint=False):
|
||||
"""
|
||||
Derivative of the total current density with respect to the thing we
|
||||
solved for
|
||||
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
|
||||
@@ -527,11 +611,12 @@ class Fields_j(Fields):
|
||||
:return: product of the derivative of the current density with respect to the field we solved for with a vector
|
||||
"""
|
||||
|
||||
return Identity()*v
|
||||
return Identity()*du_dm_v
|
||||
|
||||
|
||||
def _jDeriv_m(self, src, v, adjoint=False):
|
||||
"""
|
||||
Derivative of the total current density with respect to the inversion model. Here, we assume that the primary does not depend on the model.
|
||||
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
|
||||
@@ -575,21 +660,22 @@ class Fields_j(Fields):
|
||||
h[:,i] = h[:,i]+ 1./(1j*omega(src.freq)) * self._MeMuI * (S_m)
|
||||
return h
|
||||
|
||||
def _hSecondaryDeriv_u(self, src, v, adjoint=False):
|
||||
|
||||
def _hSecondaryDeriv_u(self, src, du_dm_v, adjoint=False):
|
||||
"""
|
||||
Derivative of the secondary magnetic 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 numpy.ndarray du_dm_v: vector to take product with
|
||||
:param bool adjoint: adjoint?
|
||||
:rtype: numpy.ndarray
|
||||
:return: product of the derivative of the secondary magnetic field with respect to the field we solved for with a vector
|
||||
"""
|
||||
|
||||
if not adjoint:
|
||||
return -1./(1j*omega(src.freq)) * self._MeMuI * (self._edgeCurl.T * (self._MfRho * v) )
|
||||
return -1./(1j*omega(src.freq)) * self._MeMuI * (self._edgeCurl.T * (self._MfRho * du_dm_v) )
|
||||
elif adjoint:
|
||||
return -1./(1j*omega(src.freq)) * self._MfRho.T * (self._edgeCurl * ( self._MeMuI.T * v))
|
||||
return -1./(1j*omega(src.freq)) * self._MfRho.T * (self._edgeCurl * ( self._MeMuI.T * du_dm_v))
|
||||
|
||||
def _hSecondaryDeriv_m(self, src, v, adjoint=False):
|
||||
"""
|
||||
@@ -625,34 +711,22 @@ class Fields_j(Fields):
|
||||
return hDeriv_m
|
||||
|
||||
|
||||
def _h(self, jSolution, srcList):
|
||||
def _hDeriv_u(self, src, du_dm_v, adjoint=False):
|
||||
"""
|
||||
Total magnetic field is sum of primary and secondary
|
||||
|
||||
:param numpy.ndarray eSolution: field we solved for
|
||||
:param list srcList: list of sources
|
||||
:rtype: numpy.ndarray
|
||||
:return: total magnetic field
|
||||
"""
|
||||
|
||||
return self._hPrimary(jSolution, srcList) + self._hSecondary(jSolution, srcList)
|
||||
|
||||
def _hDeriv_u(self, src, v, adjoint=False):
|
||||
"""
|
||||
Derivative of the total magnetic field with respect to the thing we solved for
|
||||
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 v: vector to take product with
|
||||
: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 field with respect to the field we solved for with a vector
|
||||
"""
|
||||
|
||||
return self._hSecondaryDeriv_u(src, v, adjoint)
|
||||
return self._hSecondaryDeriv_u(src, du_dm_v, adjoint)
|
||||
|
||||
def _hDeriv_m(self, src, v, adjoint=False):
|
||||
"""
|
||||
Derivative of the total magnetic field density with respect to the inversion model.
|
||||
Partial derivative of the total magnetic field density with respect to the inversion model.
|
||||
|
||||
:param SimPEG.EM.FDEM.Src src: source
|
||||
:param numpy.ndarray v: vector to take product with
|
||||
@@ -720,35 +794,24 @@ class Fields_h(Fields):
|
||||
|
||||
return hSolution
|
||||
|
||||
def _h(self, hSolution, srcList):
|
||||
"""
|
||||
Total magnetic field is sum of primary and secondary
|
||||
|
||||
:param numpy.ndarray hSolution: field we solved for
|
||||
:param list srcList: list of sources
|
||||
:rtype: numpy.ndarray
|
||||
:return: total magnetic field
|
||||
"""
|
||||
|
||||
return self._hPrimary(hSolution, srcList) + self._hSecondary(hSolution, srcList)
|
||||
|
||||
def _hDeriv_u(self, src, v, adjoint=False):
|
||||
def _hDeriv_u(self, src, du_dm_v, adjoint=False):
|
||||
"""
|
||||
Derivative of the total magnetic field with respect to the thing we
|
||||
solved for
|
||||
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 v: vector to take product with
|
||||
: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 field with respect to the field we solved for with a vector
|
||||
"""
|
||||
|
||||
return Identity()*v
|
||||
return Identity()*du_dm_v
|
||||
|
||||
def _hDeriv_m(self, src, v, adjoint=False):
|
||||
"""
|
||||
Derivative of the total magnetic field with respect to the inversion model. Here, we assume that the primary does not depend on the model.
|
||||
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
|
||||
@@ -792,21 +855,21 @@ class Fields_h(Fields):
|
||||
j[:,i] = j[:,i]+ -S_e
|
||||
return j
|
||||
|
||||
def _jSecondaryDeriv_u(self, src, v, adjoint=False):
|
||||
def _jSecondaryDeriv_u(self, src, du_dm_v, adjoint=False):
|
||||
"""
|
||||
Derivative of the secondary 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 numpy.ndarray du_dm_v: vector to take product with
|
||||
:param bool adjoint: adjoint?
|
||||
:rtype: numpy.ndarray
|
||||
:return: product of the derivative of the secondary current density with respect to the field we solved for with a vector
|
||||
"""
|
||||
|
||||
if not adjoint:
|
||||
return self._edgeCurl*v
|
||||
return self._edgeCurl*du_dm_v
|
||||
elif adjoint:
|
||||
return self._edgeCurl.T*v
|
||||
return self._edgeCurl.T*du_dm_v
|
||||
|
||||
def _jSecondaryDeriv_m(self, src, v, adjoint=False):
|
||||
"""
|
||||
@@ -822,33 +885,23 @@ class Fields_h(Fields):
|
||||
_,S_eDeriv = src.evalDeriv(self.prob, v, adjoint)
|
||||
return -S_eDeriv
|
||||
|
||||
def _j(self, hSolution, srcList):
|
||||
"""
|
||||
Total current density is sum of primary and secondary
|
||||
|
||||
:param numpy.ndarray eSolution: field we solved for
|
||||
:param list srcList: list of sources
|
||||
:rtype: numpy.ndarray
|
||||
:return: total current density
|
||||
"""
|
||||
|
||||
return self._jPrimary(hSolution, srcList) + self._jSecondary(hSolution, srcList)
|
||||
|
||||
def _jDeriv_u(self, src, v, adjoint=False):
|
||||
def _jDeriv_u(self, src, du_dm_v, adjoint=False):
|
||||
"""
|
||||
Derivative of the total current density with respect to the thing we solved for
|
||||
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 numpy.ndarray du_dm_v: vector to take product with
|
||||
:param bool adjoint: adjoint?
|
||||
:rtype: numpy.ndarray
|
||||
:return: product of the derivative of the current density with respect to the field we solved for with a vector
|
||||
"""
|
||||
return self._jSecondaryDeriv_u(src,v,adjoint)
|
||||
|
||||
return self._jSecondaryDeriv_u(src,du_dm_v,adjoint)
|
||||
|
||||
def _jDeriv_m(self, src, v, adjoint=False):
|
||||
"""
|
||||
Derivative of the total current density with respect to the inversion model.
|
||||
Partial derivative of the total current density with respect to the inversion model.
|
||||
|
||||
:param SimPEG.EM.FDEM.Src src: source
|
||||
:param numpy.ndarray v: vector to take product with
|
||||
|
||||
+85
-57
@@ -1,7 +1,7 @@
|
||||
from SimPEG import Survey, Problem, Utils, np, sp
|
||||
from scipy.constants import mu_0
|
||||
from SimPEG.EM.Utils import *
|
||||
from SimPEG.Utils import Zero
|
||||
from SimPEG.Utils import Zero
|
||||
|
||||
class BaseSrc(Survey.BaseSrc):
|
||||
"""
|
||||
@@ -14,7 +14,7 @@ class BaseSrc(Survey.BaseSrc):
|
||||
|
||||
def eval(self, prob):
|
||||
"""
|
||||
Evaluate the source terms.
|
||||
Evaluate the source terms.
|
||||
- :math:`S_m` : magnetic source term
|
||||
- :math:`S_e` : electric source term
|
||||
|
||||
@@ -36,12 +36,12 @@ class BaseSrc(Survey.BaseSrc):
|
||||
:param numpy.ndarray v: vector to take product with
|
||||
:param bool adjoint: adjoint?
|
||||
:rtype: (numpy.ndarray, numpy.ndarray)
|
||||
:return: tuple with magnetic source term and electric source term derivatives times a vector
|
||||
: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)
|
||||
if v is not None:
|
||||
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):
|
||||
"""
|
||||
@@ -49,7 +49,7 @@ class BaseSrc(Survey.BaseSrc):
|
||||
|
||||
:param Problem prob: FDEM Problem
|
||||
:rtype: numpy.ndarray
|
||||
:return: primary magnetic flux density
|
||||
:return: primary magnetic flux density
|
||||
"""
|
||||
return Zero()
|
||||
|
||||
@@ -59,7 +59,7 @@ class BaseSrc(Survey.BaseSrc):
|
||||
|
||||
:param Problem prob: FDEM Problem
|
||||
:rtype: numpy.ndarray
|
||||
:return: primary magnetic field
|
||||
:return: primary magnetic field
|
||||
"""
|
||||
return Zero()
|
||||
|
||||
@@ -69,7 +69,7 @@ class BaseSrc(Survey.BaseSrc):
|
||||
|
||||
:param Problem prob: FDEM Problem
|
||||
:rtype: numpy.ndarray
|
||||
:return: primary electric field
|
||||
:return: primary electric field
|
||||
"""
|
||||
return Zero()
|
||||
|
||||
@@ -79,13 +79,13 @@ class BaseSrc(Survey.BaseSrc):
|
||||
|
||||
:param Problem prob: FDEM Problem
|
||||
:rtype: numpy.ndarray
|
||||
:return: primary current density
|
||||
:return: primary current density
|
||||
"""
|
||||
return Zero()
|
||||
|
||||
def S_m(self, prob):
|
||||
"""
|
||||
Magnetic source term
|
||||
Magnetic source term
|
||||
|
||||
:param Problem prob: FDEM Problem
|
||||
:rtype: numpy.ndarray
|
||||
@@ -95,7 +95,7 @@ class BaseSrc(Survey.BaseSrc):
|
||||
|
||||
def S_e(self, prob):
|
||||
"""
|
||||
Electric source term
|
||||
Electric source term
|
||||
|
||||
:param Problem prob: FDEM Problem
|
||||
:rtype: numpy.ndarray
|
||||
@@ -136,15 +136,26 @@ class RawVec_e(BaseSrc):
|
||||
:param list rxList: receiver list
|
||||
:param float freq: frequency
|
||||
: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): #, 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):
|
||||
"""
|
||||
Electric source term
|
||||
|
||||
:param Problem prob: FDEM Problem
|
||||
:rtype: numpy.ndarray
|
||||
:return: electric source term on mesh
|
||||
"""
|
||||
if prob._eqLocs is 'FE' and self.integrate is True:
|
||||
return prob.Me * self._S_e
|
||||
return self._S_e
|
||||
|
||||
|
||||
@@ -155,10 +166,11 @@ class RawVec_m(BaseSrc):
|
||||
:param float freq: frequency
|
||||
:param rxList: receiver list
|
||||
: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
|
||||
|
||||
@@ -166,12 +178,14 @@ class RawVec_m(BaseSrc):
|
||||
|
||||
def S_m(self, prob):
|
||||
"""
|
||||
Magnetic source term
|
||||
Magnetic source term
|
||||
|
||||
:param Problem prob: FDEM Problem
|
||||
:rtype: numpy.ndarray
|
||||
:return: magnetic source term on mesh
|
||||
"""
|
||||
if prob._eqLocs is 'EF' and self.integrate is True:
|
||||
return prob.Me * self._S_m
|
||||
return self._S_m
|
||||
|
||||
|
||||
@@ -183,36 +197,51 @@ class RawVec(BaseSrc):
|
||||
:param float freq: frequency
|
||||
: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):
|
||||
"""
|
||||
Magnetic source term
|
||||
|
||||
:param Problem prob: FDEM Problem
|
||||
:rtype: numpy.ndarray
|
||||
:return: magnetic source term on mesh
|
||||
"""
|
||||
if prob._eqLocs is 'EF' and self.integrate is True:
|
||||
return prob.Me * self._S_m
|
||||
return self._S_m
|
||||
|
||||
def S_e(self, prob):
|
||||
"""
|
||||
Electric source term
|
||||
|
||||
:param Problem prob: FDEM Problem
|
||||
:rtype: numpy.ndarray
|
||||
:return: electric source term on mesh
|
||||
"""
|
||||
if prob._eqLocs is 'FE' and self.integrate is True:
|
||||
return prob.Me * self._S_e
|
||||
return self._S_e
|
||||
|
||||
|
||||
class MagDipole(BaseSrc):
|
||||
"""
|
||||
"""
|
||||
Point magnetic dipole source calculated by taking the curl of a magnetic
|
||||
vector potential. By taking the discrete curl, we ensure that the magnetic
|
||||
flux density is divergence free (no magnetic monopoles!).
|
||||
flux density is divergence free (no magnetic monopoles!).
|
||||
|
||||
This approach uses a primary-secondary in frequency. Here we show the
|
||||
derivation for E-B formulation noting that similar steps are followed for
|
||||
the H-J formulation.
|
||||
|
||||
.. math::
|
||||
.. math::
|
||||
\mathbf{C} \mathbf{e} + i \omega \mathbf{b} = \mathbf{s_m} \\\\
|
||||
{\mathbf{C}^T \mathbf{M_{\mu^{-1}}^f} \mathbf{b} - \mathbf{M_{\sigma}^e} \mathbf{e} = \mathbf{s_e}}
|
||||
|
||||
@@ -225,17 +254,17 @@ class MagDipole(BaseSrc):
|
||||
and define a zero-frequency primary problem, noting that the source is
|
||||
generated by a divergence free electric current
|
||||
|
||||
.. math::
|
||||
.. math::
|
||||
\mathbf{C} \mathbf{e^P} = \mathbf{s_m^P} = 0 \\\\
|
||||
{\mathbf{C}^T \mathbf{{M_{\mu^{-1}}^f}^P} \mathbf{b^P} - \mathbf{M_{\sigma}^e} \mathbf{e^P} = \mathbf{M^e} \mathbf{s_e^P}}
|
||||
|
||||
Since :math:`\mathbf{e^P}` is curl-free, divergence-free, we assume that there is no constant field background, the :math:`\mathbf{e^P} = 0`, so our primary problem is
|
||||
Since :math:`\mathbf{e^P}` is curl-free, divergence-free, we assume that there is no constant field background, the :math:`\mathbf{e^P} = 0`, so our primary problem is
|
||||
|
||||
.. math::
|
||||
.. math::
|
||||
\mathbf{e^P} = 0 \\\\
|
||||
{\mathbf{C}^T \mathbf{{M_{\mu^{-1}}^f}^P} \mathbf{b^P} = \mathbf{s_e^P}}
|
||||
|
||||
Our secondary problem is then
|
||||
Our secondary problem is then
|
||||
|
||||
.. math::
|
||||
\mathbf{C} \mathbf{e^S} + i \omega \mathbf{b^S} = - i \omega \mathbf{b^P} \\\\
|
||||
@@ -245,15 +274,15 @@ class MagDipole(BaseSrc):
|
||||
:param float freq: frequency
|
||||
:param numpy.ndarray loc: source location (ie: :code:`np.r_[xloc,yloc,zloc]`)
|
||||
:param string orientation: 'X', 'Y', 'Z'
|
||||
:param float moment: magnetic dipole moment
|
||||
:param float mu: background magnetic permeability
|
||||
:param float moment: magnetic dipole moment
|
||||
:param float mu: background magnetic permeability
|
||||
"""
|
||||
|
||||
#TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that
|
||||
def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu = mu_0):
|
||||
def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu=mu_0):
|
||||
self.freq = float(freq)
|
||||
self.loc = loc
|
||||
self.orientation = orientation
|
||||
assert orientation in ['X','Y','Z'], "Orientation (right now) doesn't actually do anything! The methods in SrcUtils should take care of this..."
|
||||
self.moment = moment
|
||||
self.mu = mu
|
||||
self.integrate = False
|
||||
@@ -265,7 +294,7 @@ class MagDipole(BaseSrc):
|
||||
|
||||
:param Problem prob: FDEM problem
|
||||
:rtype: numpy.ndarray
|
||||
:return: primary magnetic field
|
||||
:return: primary magnetic field
|
||||
"""
|
||||
eqLocs = prob._eqLocs
|
||||
|
||||
@@ -303,7 +332,7 @@ class MagDipole(BaseSrc):
|
||||
|
||||
:param Problem prob: FDEM problem
|
||||
:rtype: numpy.ndarray
|
||||
:return: primary magnetic field
|
||||
:return: primary magnetic field
|
||||
"""
|
||||
b = self.bPrimary(prob)
|
||||
return h_from_b(prob,b)
|
||||
@@ -314,7 +343,7 @@ class MagDipole(BaseSrc):
|
||||
|
||||
:param Problem prob: FDEM problem
|
||||
:rtype: numpy.ndarray
|
||||
:return: primary magnetic field
|
||||
:return: primary magnetic field
|
||||
"""
|
||||
|
||||
b_p = self.bPrimary(prob)
|
||||
@@ -326,7 +355,7 @@ class MagDipole(BaseSrc):
|
||||
|
||||
:param Problem prob: FDEM problem
|
||||
:rtype: numpy.ndarray
|
||||
:return: primary magnetic field
|
||||
:return: primary magnetic field
|
||||
"""
|
||||
|
||||
if all(np.r_[self.mu] == np.r_[prob.curModel.mu]):
|
||||
@@ -340,7 +369,7 @@ class MagDipole(BaseSrc):
|
||||
C = prob.mesh.edgeCurl
|
||||
elif eqLocs is 'EF':
|
||||
mu_s = prob.curModel.mu - self.mu
|
||||
MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True)
|
||||
MMui_s = prob.mesh.getEdgeInnerProduct(mu_s, invMat=True)
|
||||
C = prob.mesh.edgeCurl.T
|
||||
|
||||
return -C.T * (MMui_s * self.bPrimary(prob))
|
||||
@@ -353,21 +382,20 @@ class MagDipole_Bfield(BaseSrc):
|
||||
fields from a magnetic dipole. No discrete curl is taken, so the magnetic
|
||||
flux density may not be strictly divergence free.
|
||||
|
||||
This approach uses a primary-secondary in frequency in the same fashion as the MagDipole.
|
||||
This approach uses a primary-secondary in frequency in the same fashion as the MagDipole.
|
||||
|
||||
:param list rxList: receiver list
|
||||
:param float freq: frequency
|
||||
:param numpy.ndarray loc: source location (ie: :code:`np.r_[xloc,yloc,zloc]`)
|
||||
:param string orientation: 'X', 'Y', 'Z'
|
||||
:param float moment: magnetic dipole moment
|
||||
:param float mu: background magnetic permeability
|
||||
:param float moment: magnetic dipole moment
|
||||
:param float mu: background magnetic permeability
|
||||
"""
|
||||
|
||||
#TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that
|
||||
#TODO: neither does moment
|
||||
def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu = mu_0):
|
||||
self.freq = float(freq)
|
||||
self.loc = loc
|
||||
assert orientation in ['X','Y','Z'], "Orientation (right now) doesn't actually do anything! The methods in SrcUtils should take care of this..."
|
||||
self.orientation = orientation
|
||||
self.moment = moment
|
||||
self.mu = mu
|
||||
@@ -379,7 +407,7 @@ class MagDipole_Bfield(BaseSrc):
|
||||
|
||||
:param Problem prob: FDEM problem
|
||||
:rtype: numpy.ndarray
|
||||
:return: primary magnetic field
|
||||
:return: primary magnetic field
|
||||
"""
|
||||
|
||||
eqLocs = prob._eqLocs
|
||||
@@ -418,7 +446,7 @@ class MagDipole_Bfield(BaseSrc):
|
||||
|
||||
:param Problem prob: FDEM problem
|
||||
:rtype: numpy.ndarray
|
||||
:return: primary magnetic field
|
||||
:return: primary magnetic field
|
||||
"""
|
||||
b = self.bPrimary(prob)
|
||||
return h_from_b(prob, b)
|
||||
@@ -429,7 +457,7 @@ class MagDipole_Bfield(BaseSrc):
|
||||
|
||||
:param Problem prob: FDEM problem
|
||||
:rtype: numpy.ndarray
|
||||
:return: primary magnetic field
|
||||
:return: primary magnetic field
|
||||
"""
|
||||
b = self.bPrimary(prob)
|
||||
return -1j*omega(self.freq)*b
|
||||
@@ -440,7 +468,7 @@ class MagDipole_Bfield(BaseSrc):
|
||||
|
||||
:param Problem prob: FDEM problem
|
||||
:rtype: numpy.ndarray
|
||||
:return: primary magnetic field
|
||||
:return: primary magnetic field
|
||||
"""
|
||||
if all(np.r_[self.mu] == np.r_[prob.curModel.mu]):
|
||||
return Zero()
|
||||
@@ -453,7 +481,7 @@ class MagDipole_Bfield(BaseSrc):
|
||||
C = prob.mesh.edgeCurl
|
||||
elif eqLocs is 'EF':
|
||||
mu_s = prob.curModel.mu - self.mu
|
||||
MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True)
|
||||
MMui_s = prob.mesh.getEdgeInnerProduct(mu_s, invMat=True)
|
||||
C = prob.mesh.edgeCurl.T
|
||||
|
||||
return -C.T * (MMui_s * self.bPrimary(prob))
|
||||
@@ -463,22 +491,22 @@ class CircularLoop(BaseSrc):
|
||||
"""
|
||||
Circular loop magnetic source calculated by taking the curl of a magnetic
|
||||
vector potential. By taking the discrete curl, we ensure that the magnetic
|
||||
flux density is divergence free (no magnetic monopoles!).
|
||||
flux density is divergence free (no magnetic monopoles!).
|
||||
|
||||
This approach uses a primary-secondary in frequency in the same fashion as the MagDipole.
|
||||
This approach uses a primary-secondary in frequency in the same fashion as the MagDipole.
|
||||
|
||||
:param list rxList: receiver list
|
||||
:param float freq: frequency
|
||||
:param numpy.ndarray loc: source location (ie: :code:`np.r_[xloc,yloc,zloc]`)
|
||||
:param string orientation: 'X', 'Y', 'Z'
|
||||
:param float moment: magnetic dipole moment
|
||||
:param float mu: background magnetic permeability
|
||||
:param float moment: magnetic dipole moment
|
||||
:param float mu: background magnetic permeability
|
||||
"""
|
||||
|
||||
#TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that
|
||||
def __init__(self, rxList, freq, loc, orientation='Z', radius = 1., mu=mu_0):
|
||||
def __init__(self, rxList, freq, loc, orientation='Z', radius=1., mu=mu_0):
|
||||
self.freq = float(freq)
|
||||
self.orientation = orientation
|
||||
assert orientation in ['X','Y','Z'], "Orientation (right now) doesn't actually do anything! The methods in SrcUtils should take care of this..."
|
||||
self.radius = radius
|
||||
self.mu = mu
|
||||
self.loc = loc
|
||||
@@ -491,7 +519,7 @@ class CircularLoop(BaseSrc):
|
||||
|
||||
:param Problem prob: FDEM problem
|
||||
:rtype: numpy.ndarray
|
||||
:return: primary magnetic field
|
||||
:return: primary magnetic field
|
||||
"""
|
||||
eqLocs = prob._eqLocs
|
||||
|
||||
@@ -528,7 +556,7 @@ class CircularLoop(BaseSrc):
|
||||
|
||||
:param Problem prob: FDEM problem
|
||||
:rtype: numpy.ndarray
|
||||
:return: primary magnetic field
|
||||
:return: primary magnetic field
|
||||
"""
|
||||
b = self.bPrimary(prob)
|
||||
return 1./self.mu*b
|
||||
@@ -539,7 +567,7 @@ class CircularLoop(BaseSrc):
|
||||
|
||||
:param Problem prob: FDEM problem
|
||||
:rtype: numpy.ndarray
|
||||
:return: primary magnetic field
|
||||
:return: primary magnetic field
|
||||
"""
|
||||
b = self.bPrimary(prob)
|
||||
return -1j*omega(self.freq)*b
|
||||
@@ -550,7 +578,7 @@ class CircularLoop(BaseSrc):
|
||||
|
||||
:param Problem prob: FDEM problem
|
||||
:rtype: numpy.ndarray
|
||||
:return: primary magnetic field
|
||||
:return: primary magnetic field
|
||||
"""
|
||||
if all(np.r_[self.mu] == np.r_[prob.curModel.mu]):
|
||||
return Zero()
|
||||
@@ -563,7 +591,7 @@ class CircularLoop(BaseSrc):
|
||||
C = prob.mesh.edgeCurl
|
||||
elif eqLocs is 'EF':
|
||||
mu_s = prob.curModel.mu - self.mu
|
||||
MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True)
|
||||
MMui_s = prob.mesh.getEdgeInnerProduct(mu_s, invMat=True)
|
||||
C = prob.mesh.edgeCurl.T
|
||||
|
||||
return -C.T * (MMui_s * self.bPrimary(prob))
|
||||
|
||||
Reference in New Issue
Block a user