wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-08-13 12:40:08 +08:00
Code Issues Packages Projects Releases Wiki Activity

explicit integration of sources (s_e for E-B formulation, s_m for H-J formulation) in the calculation of the RHS and associated updates to the documentation

Browse Source
This commit is contained in:
Lindsey Heagy
2015-10-01 21:45:06 -07:00
parent c7c1126e9b
commit 0d88ec3652
5 changed files with 67 additions and 54 deletions
Download Patch File Download Diff File Expand all files Collapse all files
docs/api_FDEM.rst
+2 -2
View File
@@ -124,13 +124,13 @@ E-B Formulation:
.. math ::
\mathbf{C} \mathbf{e} + i \omega \mathbf{b} = \mathbf{s_m} \\
\mathbf{C^T} \mathbf{M^f_{\mu^{-1}}} \mathbf{b} - \mathbf{M^e_\sigma} \mathbf{e} = \mathbf{s_e}
\mathbf{C^T} \mathbf{M^f_{\mu^{-1}}} \mathbf{b} - \mathbf{M^e_\sigma} \mathbf{e} = \mathbf{M^e} \mathbf{s_e}
H-J Formulation:
****************
.. math ::
\mathbf{C^T} \mathbf{M^f_\rho} \mathbf{j} + i \omega \mathbf{M^e_\mu} \mathbf{h} = \mathbf{s_m} \\
\mathbf{C^T} \mathbf{M^f_\rho} \mathbf{j} + i \omega \mathbf{M^e_\mu} \mathbf{h} = \mathbf{M^e} \mathbf{s_m} \\
\mathbf{C} \mathbf{h} - \mathbf{j} = \mathbf{s_e}
simpegEM/FDEM/FDEM.py
+46 -34
View File
@@ -15,7 +15,7 @@ class BaseFDEMProblem(BaseEMProblem):
.. 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}}
{\mathbf{C}^T \mathbf{M_{\mu^{-1}}^f} \mathbf{b} - \mathbf{M_{\sigma}^e} \mathbf{e} = \mathbf{M^e} \mathbf{s_e}}
if using the E-B formulation (:code:`ProblemFDEM_e`
or :code:`ProblemFDEM_b`) or the magnetic field
@@ -23,7 +23,7 @@ class BaseFDEMProblem(BaseEMProblem):
.. math ::
\mathbf{C}^T \mathbf{M_{\\rho}^f} \mathbf{j} + i \omega \mathbf{M_{\mu}^e} \mathbf{h} = \mathbf{s_m} \\\\
\mathbf{C}^T \mathbf{M_{\\rho}^f} \mathbf{j} + i \omega \mathbf{M_{\mu}^e} \mathbf{h} = \mathbf{M^e} \mathbf{s_m} \\\\
\mathbf{C} \mathbf{h} - \mathbf{j} = \mathbf{s_e}
if using the H-J formulation (:code:`ProblemFDEM_j` or :code:`ProblemFDEM_h`).
@@ -198,7 +198,7 @@ class ProblemFDEM_e(BaseFDEMProblem):
.. math ::
\\left(\mathbf{C}^T \mathbf{M_{\mu^{-1}}^f} \mathbf{C}+ i \omega \mathbf{M^e_{\sigma}} \\right)\mathbf{e} = \mathbf{C}^T \mathbf{M_{\mu^{-1}}^f}\mathbf{s_m} -i\omega\mathbf{s_e}
\\left(\mathbf{C}^T \mathbf{M_{\mu^{-1}}^f} \mathbf{C}+ i \omega \mathbf{M^e_{\sigma}} \\right)\mathbf{e} = \mathbf{C}^T \mathbf{M_{\mu^{-1}}^f}\mathbf{s_m} -i\omega\mathbf{M^e}\mathbf{s_e}
which we solve for \\\(\\\mathbf{e}\\\).
"""
@@ -238,7 +238,7 @@ class ProblemFDEM_e(BaseFDEMProblem):
def getRHS(self, freq):
"""
.. math ::
\mathbf{RHS} = \mathbf{C}^T \mathbf{M_{\mu^{-1}}^f}\mathbf{s_m} -i\omega\mathbf{s_e}
\mathbf{RHS} = \mathbf{C}^T \mathbf{M_{\mu^{-1}}^f}\mathbf{s_m} -i\omega\mathbf{M_e}\mathbf{s_e}
:param float freq: Frequency
:rtype: numpy.ndarray (nE, nSrc)
@@ -246,16 +246,18 @@ class ProblemFDEM_e(BaseFDEMProblem):
"""
S_m, S_e = self.getSourceTerm(freq)
Me = self.Me
C = self.mesh.edgeCurl
MfMui = self.MfMui
RHS = C.T * (MfMui * S_m) -1j*omega(freq)*S_e
RHS = C.T * (MfMui * S_m) -1j * omega(freq) * Me * S_e
return RHS
def getRHSDeriv_m(self, src, v, adjoint=False):
C = self.mesh.edgeCurl
MfMui = self.MfMui
Me = self.Me
S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint)
if adjoint:
@@ -263,22 +265,22 @@ class ProblemFDEM_e(BaseFDEMProblem):
S_mDerivv = S_mDeriv(dRHS)
S_eDerivv = S_eDeriv(v)
if S_mDerivv is not None and S_eDerivv is not None:
return S_mDerivv - 1j*omega(freq)*S_eDerivv
return S_mDerivv - 1j * omega(freq) * Me.T * S_eDerivv
elif S_mDerivv is not None:
return S_mDerivv
elif S_eDerivv is not None:
return - 1j*omega(freq)*S_eDerivv
return - 1j * omega(freq) * Me.T * S_eDerivv
else:
return None
else:
S_mDerivv, S_eDerivv = S_mDeriv(v), S_eDeriv(v)
if S_mDerivv is not None and S_eDerivv is not None:
return C.T * (MfMui * S_mDerivv) -1j*omega(freq)*S_eDerivv
return C.T * (MfMui * S_mDerivv) -1j * omega(freq) * Me * S_eDerivv
elif S_mDerivv is not None:
return C.T * (MfMui * S_mDerivv)
elif S_eDerivv is not None:
return -1j*omega(freq)*S_eDerivv
return -1j * omega(freq) * Me * S_eDerivv
else:
return None
@@ -295,7 +297,7 @@ class ProblemFDEM_b(BaseFDEMProblem):
.. math ::
\\left(\mathbf{C} \mathbf{M^e_{\sigma}}^{-1} \mathbf{C}^T \mathbf{M_{\mu^{-1}}^f} + i \omega \\right)\mathbf{b} = \mathbf{s_m} + \mathbf{M^e_{\sigma}}^{-1}\mathbf{s_e}
\\left(\mathbf{C} \mathbf{M^e_{\sigma}}^{-1} \mathbf{C}^T \mathbf{M_{\mu^{-1}}^f} + i \omega \\right)\mathbf{b} = \mathbf{s_m} + \mathbf{M^e_{\sigma}}^{-1}\mathbf{M^e}\mathbf{s_e}
.. note ::
The inverse problem will not work with full anisotropy
@@ -323,7 +325,7 @@ class ProblemFDEM_b(BaseFDEMProblem):
C = self.mesh.edgeCurl
iomega = 1j * omega(freq) * sp.eye(self.mesh.nF)
A = C*MeSigmaI*C.T*MfMui + iomega
A = C * (MeSigmaI * (C.T * MfMui)) + iomega
if self._makeASymmetric is True:
return MfMui.T*A
@@ -334,7 +336,7 @@ class ProblemFDEM_b(BaseFDEMProblem):
MfMui = self.MfMui
C = self.mesh.edgeCurl
MeSigmaIDeriv = self.MeSigmaIDeriv
vec = C.T*(MfMui*u)
vec = C.T * (MfMui * u)
MeSigmaIDeriv = MeSigmaIDeriv(vec)
@@ -361,12 +363,13 @@ class ProblemFDEM_b(BaseFDEMProblem):
S_m, S_e = self.getSourceTerm(freq)
C = self.mesh.edgeCurl
MeSigmaI = self.MeSigmaI
Me = self.Me
RHS = S_m + C * ( MeSigmaI * S_e )
RHS = S_m + C * ( MeSigmaI * Me * S_e )
if self._makeASymmetric is True:
MfMui = self.MfMui
return MfMui.T*RHS
return MfMui.T * RHS
return RHS
@@ -374,12 +377,13 @@ class ProblemFDEM_b(BaseFDEMProblem):
C = self.mesh.edgeCurl
S_m, S_e = src.eval(self)
MfMui = self.MfMui
Me = self.Me
if self._makeASymmetric and adjoint:
v = self.MfMui * v
if S_e is not None:
MeSigmaIDeriv = self.MeSigmaIDeriv(Utils.mkvc(S_e))
MeSigmaIDeriv = self.MeSigmaIDeriv(Utils.mkvc( Me * S_e))
if not adjoint:
RHSderiv = C * (MeSigmaIDeriv * v)
elif adjoint:
@@ -391,16 +395,16 @@ class ProblemFDEM_b(BaseFDEMProblem):
S_mDeriv, S_eDeriv = S_mDeriv(v), S_eDeriv(v)
if S_mDeriv is not None and S_eDeriv is not None:
if not adjoint:
SrcDeriv = S_mDeriv + C * (self.MeSigmaI * S_eDeriv)
SrcDeriv = S_mDeriv + C * (self.MeSigmaI * (Me * S_eDeriv))
elif adjoint:
SrcDeriv = S_mDeriv + Self.MeSigmaI.T * ( C.T * S_eDeriv)
SrcDeriv = S_mDeriv + Me.T * (Self.MeSigmaI.T * ( C.T * S_eDeriv))
elif S_mDeriv is not None:
SrcDeriv = S_mDeriv
elif S_eDeriv is not None:
if not adjoint:
SrcDeriv = C * (self.MeSigmaI * S_eDeriv)
SrcDeriv = C * (self.MeSigmaI * (Me * S_eDeriv))
elif adjoint:
SrcDeriv = self.MeSigmaI.T * ( C.T * S_eDeriv)
SrcDeriv = Me.T * (self.MeSigmaI.T * ( C.T * S_eDeriv))
else:
SrcDeriv = None
@@ -428,13 +432,13 @@ class ProblemFDEM_j(BaseFDEMProblem):
.. math ::
\mathbf{h} = \\frac{1}{i \omega} \mathbf{M_{\mu}^e}^{-1} \\left(-\mathbf{C}^T \mathbf{M_{\\rho}^f} \mathbf{j} + \mathbf{s_m} \\right)
\mathbf{h} = \\frac{1}{i \omega} \mathbf{M_{\mu}^e}^{-1} \\left(-\mathbf{C}^T \mathbf{M_{\\rho}^f} \mathbf{j} + \mathbf{M^e} \mathbf{s_m} \\right)
and solve for \\\(\\\mathbf{j}\\\) using
.. math ::
\\left(\mathbf{C} \mathbf{M_{\mu}^e}^{-1} \mathbf{C}^T \mathbf{M_{\\rho}^f} + i \omega\\right)\mathbf{j} = \mathbf{C} \mathbf{M_{\mu}^e}^{-1}\mathbf{s_m} -i\omega\mathbf{s_e}
\\left(\mathbf{C} \mathbf{M_{\mu}^e}^{-1} \mathbf{C}^T \mathbf{M_{\\rho}^f} + i \omega\\right)\mathbf{j} = \mathbf{C} \mathbf{M_{\mu}^e}^{-1} \mathbf{M^e} \mathbf{s_m} -i\omega\mathbf{s_e}
.. note::
This implementation does not yet work with full anisotropy!!
@@ -507,10 +511,11 @@ class ProblemFDEM_j(BaseFDEMProblem):
S_m, S_e = self.getSourceTerm(freq)
C = self.mesh.edgeCurl
MeMuI = self.MeMuI
MeMuI = self.MeMuI
Me = self.Me
RHS = C * (MeMuI * S_m) - 1j * omega(freq) * S_e
RHS = C * (MeMuI * (Me * S_m)) - 1j * omega(freq) * S_e
if self._makeASymmetric is True:
MfRho = self.MfRho
return MfRho.T*RHS
@@ -520,6 +525,7 @@ class ProblemFDEM_j(BaseFDEMProblem):
def getRHSDeriv_m(self, src, v, adjoint=False):
C = self.mesh.edgeCurl
MeMuI = self.MeMuI
Me = self.Me
S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint)
if adjoint:
@@ -529,20 +535,20 @@ class ProblemFDEM_j(BaseFDEMProblem):
S_mDerivv = S_mDeriv(MeMuI.T * (C.T * v))
S_eDerivv = S_eDeriv(v)
if S_mDerivv is not None and S_eDerivv is not None:
return S_mDerivv - 1j*omega(freq)*S_eDerivv
return Me.T * S_mDerivv - 1j * omega(freq) * S_eDerivv
elif S_mDerivv is not None:
return S_mDerivv
return Me.T * S_mDerivv
elif S_eDerivv is not None:
return - 1j*omega(freq)*S_eDerivv
return - 1j * omega(freq) * S_eDerivv
else:
return None
else:
S_mDerivv, S_eDerivv = S_mDeriv(v), S_eDeriv(v)
if S_mDerivv is not None and S_eDerivv is not None:
RHSDeriv = C * (MeMuI * S_mDerivv) - 1j * omega(freq) * S_eDerivv
RHSDeriv = C * (MeMuI * (Me * S_mDerivv)) - 1j * omega(freq) * S_eDerivv
elif S_mDerivv is not None:
RHSDeriv = C * (MeMuI * S_mDerivv)
RHSDeriv = C * (MeMuI * (Me * S_mDerivv))
elif S_eDerivv is not None:
RHSDeriv = - 1j * omega(freq) * S_eDerivv
else:
@@ -568,7 +574,7 @@ class ProblemFDEM_h(BaseFDEMProblem):
.. math ::
\\left(\mathbf{C}^T \mathbf{M_{\\rho}^f} \mathbf{C} + i \omega \mathbf{M_{\mu}^e}\\right) \mathbf{h} = \mathbf{s_m} + \mathbf{C}^T \mathbf{M_{\\rho}^f} \mathbf{s_e}
\\left(\mathbf{C}^T \mathbf{M_{\\rho}^f} \mathbf{C} + i \omega \mathbf{M_{\mu}^e}\\right) \mathbf{h} = \mathbf{M^e} \mathbf{s_m} + \mathbf{C}^T \mathbf{M_{\\rho}^f} \mathbf{s_e}
"""
@@ -594,7 +600,7 @@ class ProblemFDEM_h(BaseFDEMProblem):
MfRho = self.MfRho
C = self.mesh.edgeCurl
return C.T * MfRho * C + 1j*omega(freq)*MeMu
return C.T * (MfRho * C) + 1j*omega(freq)*MeMu
def getADeriv_m(self, freq, u, v, adjoint=False):
@@ -610,7 +616,7 @@ class ProblemFDEM_h(BaseFDEMProblem):
"""
.. math ::
\mathbf{RHS} = \mathbf{s_m} + \mathbf{C}^T \mathbf{M_{\\rho}^f} \mathbf{s_e}
\mathbf{RHS} = \mathbf{M^e} \mathbf{s_m} + \mathbf{C}^T \mathbf{M_{\\rho}^f} \mathbf{s_e}
:param float freq: Frequency
:rtype: numpy.ndarray (nE, nSrc)
@@ -620,8 +626,9 @@ class ProblemFDEM_h(BaseFDEMProblem):
S_m, S_e = self.getSourceTerm(freq)
C = self.mesh.edgeCurl
MfRho = self.MfRho
Me = self.Me
RHS = S_m + C.T * ( MfRho * S_e )
RHS = Me * S_m + C.T * ( MfRho * S_e )
return RHS
@@ -629,6 +636,7 @@ class ProblemFDEM_h(BaseFDEMProblem):
_, S_e = src.eval(self)
C = self.mesh.edgeCurl
MfRho = self.MfRho
Me = self.Me
RHSDeriv = None
@@ -643,11 +651,15 @@ class ProblemFDEM_h(BaseFDEMProblem):
S_mDeriv = S_mDeriv(v)
S_eDeriv = S_eDeriv(v)
if adjoint:
Me = Me.T
if S_mDeriv is not None:
if RHSDeriv is not None:
RHSDeriv += S_mDeriv(v)
RHSDeriv += Me * S_mDeriv(v)
else:
RHSDeriv = S_mDeriv(v)
RHSDeriv = Me * S_mDeriv(v)
if S_eDeriv is not None:
if RHSDeriv is not None:
RHSDeriv += C.T * (MfRho * S_e)
simpegEM/FDEM/FieldsFDEM.py
+16 -15
View File
@@ -109,6 +109,7 @@ class FieldsFDEM_b(FieldsFDEM):
self._MeSigmaI = self.survey.prob.MeSigmaI
self._MfMui = self.survey.prob.MfMui
self._MeSigmaIDeriv = self.survey.prob.MeSigmaIDeriv
self._Me = self.survey.prob.Me
def _bPrimary(self, bSolution, srcList):
bPrimary = np.zeros_like(bSolution)
@@ -144,7 +145,7 @@ class FieldsFDEM_b(FieldsFDEM):
for i,src in enumerate(srcList):
_,S_e = src.eval(self.prob)
if S_e is not None:
e[:,i] += -self._MeSigmaI*S_e
e[:,i] += -self._MeSigmaI * (self._Me * S_e)
return e
def _eSecondaryDeriv_u(self, src, v, adjoint=False):
@@ -156,10 +157,14 @@ class FieldsFDEM_b(FieldsFDEM):
def _eSecondaryDeriv_m(self, src, v, adjoint=False):
bSolution = self[[src],'bSolution']
_,S_e = src.eval(self.prob)
Me = self._Me
if adjoint:
Me = Me.T
w = self._edgeCurl.T * (self._MfMui * bSolution)
if S_e is not None:
w += -Utils.mkvc(S_e,2)
w += -Utils.mkvc(Me * S_e,2)
if not adjoint:
de_dm = self._MeSigmaIDeriv(w) * v
@@ -170,7 +175,7 @@ class FieldsFDEM_b(FieldsFDEM):
Se_Deriv = S_eDeriv(v)
if Se_Deriv is not None:
de_dm += -self._MeSigmaI * Se_Deriv
de_dm += -self._MeSigmaI * (self._Me * Se_Deriv)
return de_dm
@@ -205,6 +210,7 @@ class FieldsFDEM_j(FieldsFDEM):
self._MeMuI = self.survey.prob.MeMuI
self._MfRho = self.survey.prob.MfRho
self._MfRhoDeriv = self.survey.prob.MfRhoDeriv
self._Me = self.survey.prob.Me
def _jPrimary(self, jSolution, srcList):
jPrimary = np.zeros_like(jSolution,dtype = complex)
@@ -236,25 +242,19 @@ class FieldsFDEM_j(FieldsFDEM):
return hPrimary
def _hSecondary(self, jSolution, srcList):
MeMuI = self._MeMuI
C = self._edgeCurl
MfRho = self._MfRho
h = MeMuI * (C.T * (MfRho * jSolution) )
h = self._MeMuI * (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)
if S_m is not None:
h[:,i] += 1./(1j*omega(src.freq)) * MeMuI * S_m
h[:,i] += 1./(1j*omega(src.freq)) * self._MeMuI * (self._Me * S_m)
return h
def _hSecondaryDeriv_u(self, src, v, adjoint=False):
MeMuI = self._MeMuI
C = self._edgeCurl
MfRho = self._MfRho
if not adjoint:
return -1./(1j*omega(src.freq)) * MeMuI * (C.T * (MfRho * v) )
return -1./(1j*omega(src.freq)) * self._MeMuI * (self._edgeCurl.T * (self._MfRho * v) )
elif adjoint:
return -1./(1j*omega(src.freq)) * MfRho.T * (C * ( MeMuI.T * v))
return -1./(1j*omega(src.freq)) * self._MfRho.T * (self._edgeCurl * ( self._MeMuI.T * v))
def _hSecondaryDeriv_m(self, src, v, adjoint=False):
jSolution = self[[src],'jSolution']
@@ -262,6 +262,7 @@ class FieldsFDEM_j(FieldsFDEM):
C = self._edgeCurl
MfRho = self._MfRho
MfRhoDeriv = self._MfRhoDeriv
Me = self._Me
if not adjoint:
hDeriv_m = -1./(1j*omega(src.freq)) * MeMuI * (C.T * (MfRhoDeriv(jSolution)*v ) )
@@ -273,9 +274,9 @@ class FieldsFDEM_j(FieldsFDEM):
if not adjoint:
S_mDeriv = S_mDeriv(v)
if S_mDeriv is not None:
hDeriv_m += 1./(1j*omega(src.freq)) * MeMuI * S_mDeriv
hDeriv_m += 1./(1j*omega(src.freq)) * MeMuI * (Me * S_mDeriv)
elif adjoint:
S_mDeriv = S_mDeriv(MeMuI.T * v)
S_mDeriv = S_mDeriv(Me.T * (MeMuI.T * v))
if S_mDeriv is not None:
hDeriv_m += 1./(1j*omega(src.freq)) * S_mDeriv
return hDeriv_m
simpegEM/TDEM/SurveyTDEM.py
+1 -1
View File
@@ -102,7 +102,7 @@ class SrcTDEM_CircularLoop_MVP(SrcTDEM):
def __init__(self,rxList,loc,radius):
self.loc = loc
self.radius =radius
self.radius = radius
SrcTDEM.__init__(self,rxList)
def getInitialFields(self, mesh):
simpegEM/Tests/test_FDEM.py
+2 -2
View File
@@ -12,14 +12,14 @@ testHJ = True
verbose = False
TOL = 1e-4
TOL = 1e-6
FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order
CONDUCTIVITY = 1e1
MU = mu_0
freq = 1e-1
addrandoms = True
SrcType = 'MagDipole' #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVec'
SrcType = 'RawVec' #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVec'
def getProblem(fdemType, comp):