mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-25 13:30:06 +08:00
Merge pull request #34 from simpeg/SrcIntegration
Src integration: Inside of the source definitions, there is now the option of whether to integrate the source term by multiplying with an edge inner product matrix
This commit is contained in:
Lindsey
+4
-2
@@ -2,6 +2,8 @@ language: python
|
||||
python:
|
||||
- 2.7
|
||||
|
||||
sudo: false
|
||||
|
||||
# Setup anaconda
|
||||
before_install:
|
||||
- if [ ${TRAVIS_PYTHON_VERSION:0:1} == "2" ]; then wget http://repo.continuum.io/miniconda/Miniconda-3.8.3-Linux-x86_64.sh -O miniconda.sh; else wget http://repo.continuum.io/miniconda/Miniconda3-3.8.3-Linux-x86_64.sh -O miniconda.sh; fi
|
||||
@@ -10,8 +12,8 @@ before_install:
|
||||
- export PATH=/home/travis/anaconda/bin:/home/travis/miniconda/bin:$PATH
|
||||
- conda update --yes conda
|
||||
# The next couple lines fix a crash with multiprocessing on Travis and are not specific to using Miniconda
|
||||
- sudo rm -rf /dev/shm
|
||||
- sudo ln -s /run/shm /dev/shm
|
||||
# - sudo rm -rf /dev/shm
|
||||
# - sudo ln -s /run/shm /dev/shm
|
||||
|
||||
# Install packages
|
||||
install:
|
||||
|
||||
+2
-2
@@ -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}
|
||||
|
||||
|
||||
|
||||
@@ -42,7 +42,7 @@ def hzAnalyticDipoleF(r, freq, sigma, secondary=True, mu=mu_0):
|
||||
|
||||
return hz
|
||||
|
||||
def AnalyticMagDipoleWholeSpace(XYZ, srcLoc, sig, f, moment=1., orientation='X', mu = mu_0):
|
||||
def MagneticDipoleWholeSpace(XYZ, srcLoc, sig, f, moment=1., orientation='X', mu = mu_0):
|
||||
"""
|
||||
Analytical solution for a dipole in a whole-space.
|
||||
|
||||
|
||||
+37
-23
@@ -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,20 @@ class ProblemFDEM_e(BaseFDEMProblem):
|
||||
"""
|
||||
|
||||
S_m, S_e = self.getSourceTerm(freq)
|
||||
# if self.survey.
|
||||
# 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
|
||||
RHS = C.T * (MfMui * S_m) -1j * omega(freq) * 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 +267,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) * 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) * 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) * 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) * S_eDerivv
|
||||
else:
|
||||
return None
|
||||
|
||||
@@ -295,7 +299,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 +327,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 +338,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 +365,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 )
|
||||
|
||||
if self._makeASymmetric is True:
|
||||
MfMui = self.MfMui
|
||||
return MfMui.T*RHS
|
||||
return MfMui.T * RHS
|
||||
|
||||
return RHS
|
||||
|
||||
@@ -374,12 +379,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(S_e)
|
||||
if not adjoint:
|
||||
RHSderiv = C * (MeSigmaIDeriv * v)
|
||||
elif adjoint:
|
||||
@@ -428,13 +434,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,7 +513,8 @@ 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
|
||||
@@ -520,6 +527,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,11 +537,11 @@ 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 S_mDerivv - 1j * omega(freq) * 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) * S_eDerivv
|
||||
else:
|
||||
return None
|
||||
else:
|
||||
@@ -568,7 +576,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 +602,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 +618,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,6 +628,7 @@ 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 )
|
||||
|
||||
@@ -629,6 +638,7 @@ class ProblemFDEM_h(BaseFDEMProblem):
|
||||
_, S_e = src.eval(self)
|
||||
C = self.mesh.edgeCurl
|
||||
MfRho = self.MfRho
|
||||
Me = self.Me
|
||||
|
||||
RHSDeriv = None
|
||||
|
||||
@@ -643,11 +653,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)
|
||||
else:
|
||||
RHSDeriv = S_mDeriv(v)
|
||||
RHSDeriv = S_mDeriv(v)
|
||||
if S_eDeriv is not None:
|
||||
if RHSDeriv is not None:
|
||||
RHSDeriv += C.T * (MfRho * S_e)
|
||||
|
||||
+15
-14
@@ -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 * 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
|
||||
@@ -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 * (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
|
||||
|
||||
@@ -94,6 +94,7 @@ class RxFDEM(Survey.BaseRx):
|
||||
class SrcFDEM(Survey.BaseSrc):
|
||||
freq = None
|
||||
rxPair = RxFDEM
|
||||
integrate = True
|
||||
|
||||
def eval(self, prob):
|
||||
S_m = self.S_m(prob)
|
||||
@@ -155,9 +156,10 @@ class SrcFDEM_RawVec_m(SrcFDEM):
|
||||
:param rxList: receiver list
|
||||
"""
|
||||
|
||||
def __init__(self, rxList, freq, S_m):
|
||||
def __init__(self, rxList, freq, S_m, integrate = True):
|
||||
self._S_m = np.array(S_m,dtype=complex)
|
||||
self.freq = float(freq)
|
||||
self.integrate = integrate
|
||||
SrcFDEM.__init__(self, rxList)
|
||||
|
||||
def S_m(self, prob):
|
||||
@@ -173,16 +175,21 @@ class SrcFDEM_RawVec(SrcFDEM):
|
||||
:param float freq: frequency
|
||||
:param rxList: receiver list
|
||||
"""
|
||||
def __init__(self, rxList, freq, S_m, S_e):
|
||||
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
|
||||
SrcFDEM.__init__(self, rxList)
|
||||
|
||||
def S_m(self, prob):
|
||||
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):
|
||||
if prob._eqLocs is 'FE' and self.integrate is True:
|
||||
return prob.Me * self._S_e
|
||||
return self._S_e
|
||||
|
||||
|
||||
@@ -194,7 +201,8 @@ class SrcFDEM_MagDipole(SrcFDEM):
|
||||
self.loc = loc
|
||||
self.orientation = orientation
|
||||
self.moment = moment
|
||||
self.mu = mu
|
||||
self.mu = mu
|
||||
self.integrate = False
|
||||
SrcFDEM.__init__(self, rxList)
|
||||
|
||||
def bPrimary(self, prob):
|
||||
@@ -332,6 +340,7 @@ class SrcFDEM_CircularLoop(SrcFDEM):
|
||||
self.radius = radius
|
||||
self.mu = mu
|
||||
self.loc = loc
|
||||
self.integrate = False
|
||||
SrcFDEM.__init__(self, rxList)
|
||||
|
||||
def bPrimary(self, prob):
|
||||
|
||||
@@ -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):
|
||||
|
||||
@@ -12,14 +12,14 @@ testHJ = True
|
||||
|
||||
verbose = False
|
||||
|
||||
TOL = 1e-4
|
||||
TOL = 1e-5
|
||||
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):
|
||||
|
||||
@@ -4,7 +4,11 @@ import simpegEM as EM
|
||||
from scipy.constants import mu_0
|
||||
|
||||
plotIt = False
|
||||
freq = 1e2
|
||||
tol_EBdipole = 1e-2
|
||||
|
||||
if plotIt:
|
||||
import matplotlib.pylab
|
||||
|
||||
|
||||
class FDEM_analyticTests(unittest.TestCase):
|
||||
|
||||
@@ -13,6 +17,8 @@ class FDEM_analyticTests(unittest.TestCase):
|
||||
cs = 10.
|
||||
ncx, ncy, ncz = 10, 10, 10
|
||||
npad = 4
|
||||
freq = 1e2
|
||||
|
||||
hx = [(cs,npad,-1.3), (cs,ncx), (cs,npad,1.3)]
|
||||
hy = [(cs,npad,-1.3), (cs,ncy), (cs,npad,1.3)]
|
||||
hz = [(cs,npad,-1.3), (cs,ncz), (cs,npad,1.3)]
|
||||
@@ -78,5 +84,160 @@ class FDEM_analyticTests(unittest.TestCase):
|
||||
self.assertTrue(passed)
|
||||
|
||||
|
||||
def test_CylMeshEBDipoles(self):
|
||||
print 'Testing CylMesh Electric and Magnetic Dipoles in a wholespace- Analytic: J-formulation'
|
||||
sigmaback = 1.
|
||||
mur = 2.
|
||||
freq = 1.
|
||||
skdpth = 500./np.sqrt(sigmaback*freq)
|
||||
|
||||
csx, ncx, npadx = 5, 50, 25
|
||||
csz, ncz, npadz = 5, 50, 25
|
||||
hx = Utils.meshTensor([(csx,ncx), (csx,npadx,1.3)])
|
||||
hz = Utils.meshTensor([(csz,npadz,-1.3), (csz,ncz), (csz,npadz,1.3)])
|
||||
mesh = Mesh.CylMesh([hx,1,hz], [0.,0.,-hz.sum()/2]) # define the cylindrical mesh
|
||||
|
||||
if plotIt:
|
||||
mesh.plotGrid()
|
||||
|
||||
# make sure mesh is big enough
|
||||
self.assertTrue(mesh.hz.sum() > skdpth*2.)
|
||||
self.assertTrue(mesh.hx.sum() > skdpth*2.)
|
||||
|
||||
SigmaBack = sigmaback*np.ones((mesh.nC))
|
||||
MuBack = mur*mu_0*np.ones((mesh.nC))
|
||||
|
||||
# set up source
|
||||
# test electric dipole
|
||||
src_loc = np.r_[0.,0.,0.]
|
||||
s_ind = Utils.closestPoints(mesh,src_loc,'Fz') + mesh.nFx
|
||||
|
||||
de = np.zeros(mesh.nF,dtype=complex)
|
||||
de[s_ind] = 1./csz
|
||||
de_p = [EM.FDEM.SrcFDEM_RawVec_e([],freq,de/mesh.area)]
|
||||
|
||||
dm_p = [EM.FDEM.SrcFDEM_MagDipole([],freq,src_loc)]
|
||||
|
||||
|
||||
# Pair the problem and survey
|
||||
surveye = EM.FDEM.SurveyFDEM(de_p)
|
||||
surveym = EM.FDEM.SurveyFDEM(dm_p)
|
||||
|
||||
mapping = [('sigma', Maps.IdentityMap(mesh)),('mu', Maps.IdentityMap(mesh))]
|
||||
|
||||
prbe = EM.FDEM.ProblemFDEM_h(mesh, mapping=mapping)
|
||||
prbm = EM.FDEM.ProblemFDEM_e(mesh, mapping=mapping)
|
||||
|
||||
prbe.pair(surveye) # pair problem and survey
|
||||
prbm.pair(surveym)
|
||||
|
||||
# solve
|
||||
fieldsBackE = prbe.fields(np.r_[SigmaBack, MuBack]) # Done
|
||||
fieldsBackM = prbm.fields(np.r_[SigmaBack, MuBack]) # Done
|
||||
|
||||
|
||||
rlim = [20.,500.]
|
||||
lookAtTx = de_p
|
||||
r = mesh.vectorCCx[np.argmin(np.abs(mesh.vectorCCx-rlim[0])):np.argmin(np.abs(mesh.vectorCCx-rlim[1]))]
|
||||
z = 100.
|
||||
|
||||
# where we choose to measure
|
||||
XYZ = Utils.ndgrid(r, np.r_[0.], np.r_[z])
|
||||
|
||||
Pf = mesh.getInterpolationMat(XYZ, 'CC')
|
||||
Zero = sp.csr_matrix(Pf.shape)
|
||||
Pfx,Pfz = sp.hstack([Pf,Zero]),sp.hstack([Zero,Pf])
|
||||
|
||||
jn = fieldsBackE[de_p,'j']
|
||||
bn = fieldsBackM[dm_p,'b']
|
||||
|
||||
Rho = Utils.sdiag(1./SigmaBack)
|
||||
Rho = sp.block_diag([Rho,Rho])
|
||||
|
||||
en = Rho*mesh.aveF2CCV*jn
|
||||
bn = mesh.aveF2CCV*bn
|
||||
|
||||
ex,ez = Pfx*en, Pfz*en
|
||||
bx,bz = Pfx*bn, Pfz*bn
|
||||
|
||||
# get analytic solution
|
||||
exa, eya, eza = EM.Analytics.FDEM.ElectricDipoleWholeSpace(XYZ, src_loc, sigmaback, freq,orientation='Z',mu= mur*mu_0)
|
||||
exa, eya, eza = Utils.mkvc(exa,2), Utils.mkvc(eya,2), Utils.mkvc(eza,2)
|
||||
|
||||
bxa, bya, bza = EM.Analytics.FDEM.MagneticDipoleWholeSpace(XYZ, src_loc, sigmaback, freq,orientation='Z',mu= mur*mu_0)
|
||||
bxa, bya, bza = Utils.mkvc(bxa,2), Utils.mkvc(bya,2), Utils.mkvc(bza,2)
|
||||
|
||||
print ' comp, anayltic, numeric, num - ana, (num - ana)/ana'
|
||||
print ' ex:', np.linalg.norm(exa), np.linalg.norm(ex), np.linalg.norm(exa-ex), np.linalg.norm(exa-ex)/np.linalg.norm(exa)
|
||||
print ' ez:', np.linalg.norm(eza), np.linalg.norm(ez), np.linalg.norm(eza-ez), np.linalg.norm(eza-ez)/np.linalg.norm(eza)
|
||||
|
||||
print ' bx:', np.linalg.norm(bxa), np.linalg.norm(bx), np.linalg.norm(bxa-bx), np.linalg.norm(bxa-bx)/np.linalg.norm(bxa)
|
||||
print ' bz:', np.linalg.norm(bza), np.linalg.norm(bz), np.linalg.norm(bza-bz), np.linalg.norm(bza-bz)/np.linalg.norm(bza)
|
||||
|
||||
if plotIt:
|
||||
# Edipole
|
||||
plt.subplot(221)
|
||||
plt.plot(r,ex.real,'o',r,exa.real,linewidth=2)
|
||||
plt.grid(which='both')
|
||||
plt.title('Ex Real')
|
||||
plt.xlabel('r (m)')
|
||||
|
||||
plt.subplot(222)
|
||||
plt.plot(r,ex.imag,'o',r,exa.imag,linewidth=2)
|
||||
plt.grid(which='both')
|
||||
plt.title('Ex Imag')
|
||||
plt.legend(['Num','Ana'],bbox_to_anchor=(1.5,0.5))
|
||||
plt.xlabel('r (m)')
|
||||
|
||||
plt.subplot(223)
|
||||
plt.plot(r,ez.real,'o',r,eza.real,linewidth=2)
|
||||
plt.grid(which='both')
|
||||
plt.title('Ez Real')
|
||||
plt.xlabel('r (m)')
|
||||
|
||||
plt.subplot(224)
|
||||
plt.plot(r,ez.imag,'o',r,eza.imag,linewidth=2)
|
||||
plt.grid(which='both')
|
||||
plt.title('Ez Imag')
|
||||
plt.xlabel('r (m)')
|
||||
|
||||
plt.tight_layout()
|
||||
|
||||
# Bdipole
|
||||
plt.subplot(221)
|
||||
plt.plot(r,bx.real,'o',r,bxa.real,linewidth=2)
|
||||
plt.grid(which='both')
|
||||
plt.title('Bx Real')
|
||||
plt.xlabel('r (m)')
|
||||
|
||||
plt.subplot(222)
|
||||
plt.plot(r,bx.imag,'o',r,bxa.imag,linewidth=2)
|
||||
plt.grid(which='both')
|
||||
plt.title('Bx Imag')
|
||||
plt.legend(['Num','Ana'],bbox_to_anchor=(1.5,0.5))
|
||||
plt.xlabel('r (m)')
|
||||
|
||||
plt.subplot(223)
|
||||
plt.plot(r,bz.real,'o',r,bza.real,linewidth=2)
|
||||
plt.grid(which='both')
|
||||
plt.title('Bz Real')
|
||||
plt.xlabel('r (m)')
|
||||
|
||||
plt.subplot(224)
|
||||
plt.plot(r,bz.imag,'o',r,bza.imag,linewidth=2)
|
||||
plt.grid(which='both')
|
||||
plt.title('Bz Imag')
|
||||
plt.xlabel('r (m)')
|
||||
|
||||
plt.tight_layout()
|
||||
|
||||
self.assertTrue(np.linalg.norm(exa-ex)/np.linalg.norm(exa) < tol_EBdipole)
|
||||
self.assertTrue(np.linalg.norm(eza-ez)/np.linalg.norm(eza) < tol_EBdipole)
|
||||
|
||||
self.assertTrue(np.linalg.norm(bxa-bx)/np.linalg.norm(bxa) < tol_EBdipole)
|
||||
self.assertTrue(np.linalg.norm(bza-bz)/np.linalg.norm(bza) < tol_EBdipole)
|
||||
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
|
||||
Reference in New Issue
Block a user