From 0d88ec3652d6658b230a6c14e0443126c91929d3 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Thu, 1 Oct 2015 21:45:06 -0700 Subject: [PATCH 1/6] 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 --- docs/api_FDEM.rst | 4 +- simpegEM/FDEM/FDEM.py | 80 +++++++++++++++++++++---------------- simpegEM/FDEM/FieldsFDEM.py | 31 +++++++------- simpegEM/TDEM/SurveyTDEM.py | 2 +- simpegEM/Tests/test_FDEM.py | 4 +- 5 files changed, 67 insertions(+), 54 deletions(-) diff --git a/docs/api_FDEM.rst b/docs/api_FDEM.rst index 4c173568..28317098 100644 --- a/docs/api_FDEM.rst +++ b/docs/api_FDEM.rst @@ -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} diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py index 546c2ac2..38de9a66 100644 --- a/simpegEM/FDEM/FDEM.py +++ b/simpegEM/FDEM/FDEM.py @@ -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) diff --git a/simpegEM/FDEM/FieldsFDEM.py b/simpegEM/FDEM/FieldsFDEM.py index 99866bc3..1c57c1ca 100644 --- a/simpegEM/FDEM/FieldsFDEM.py +++ b/simpegEM/FDEM/FieldsFDEM.py @@ -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 diff --git a/simpegEM/TDEM/SurveyTDEM.py b/simpegEM/TDEM/SurveyTDEM.py index 63cc008a..bcd83962 100644 --- a/simpegEM/TDEM/SurveyTDEM.py +++ b/simpegEM/TDEM/SurveyTDEM.py @@ -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): diff --git a/simpegEM/Tests/test_FDEM.py b/simpegEM/Tests/test_FDEM.py index 24c9205f..7f91236b 100644 --- a/simpegEM/Tests/test_FDEM.py +++ b/simpegEM/Tests/test_FDEM.py @@ -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): From bac28ba133c192a8dd88608c8201ff7c31177b5d Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Thu, 1 Oct 2015 22:45:57 -0700 Subject: [PATCH 2/6] updated tolerance for travis --- simpegEM/Tests/test_FDEM.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/simpegEM/Tests/test_FDEM.py b/simpegEM/Tests/test_FDEM.py index 7f91236b..938af8b9 100644 --- a/simpegEM/Tests/test_FDEM.py +++ b/simpegEM/Tests/test_FDEM.py @@ -12,7 +12,7 @@ testHJ = True verbose = False -TOL = 1e-6 +TOL = 1e-5 FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order CONDUCTIVITY = 1e1 MU = mu_0 From 97f27c76d1b9beb67158c8e5ab2c5bdec5d93558 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Thu, 1 Oct 2015 23:41:10 -0700 Subject: [PATCH 3/6] added test for H-J formulation against analytic soon for an electric dipole in a whole-space on cyl mesh --- simpegEM/Tests/test_FDEM_analytics.py | 108 +++++++++++++++++++++++++- 1 file changed, 107 insertions(+), 1 deletion(-) diff --git a/simpegEM/Tests/test_FDEM_analytics.py b/simpegEM/Tests/test_FDEM_analytics.py index c0fc4210..21aa9fda 100644 --- a/simpegEM/Tests/test_FDEM_analytics.py +++ b/simpegEM/Tests/test_FDEM_analytics.py @@ -4,7 +4,11 @@ import simpegEM as EM from scipy.constants import mu_0 plotIt = False -freq = 1e2 +tol_Edipole = 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,105 @@ class FDEM_analyticTests(unittest.TestCase): self.assertTrue(passed) + def test_CylMeshEDipole(self): + print 'Testing CylMesh E Dipole in a wholespace- Analytic: J-formulation' + sigmaback = 1. + 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)) + + # 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)] + + # Pair the problem and survey + surveye = EM.FDEM.SurveyFDEM(de_p) #+deg_p) # set survey + mapping = [('sigma', Maps.IdentityMap(mesh))] + prbe = EM.FDEM.ProblemFDEM_j(mesh, mapping=mapping) + prbe.pair(surveye) # pair problem and survey + + # solve + fieldsBack = prbe.fields(np.r_[SigmaBack]) # 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 = fieldsBack[lookAtTx,'j'] + Rho = Utils.sdiag(1./SigmaBack) + Rho = sp.block_diag([Rho,Rho]) + + en = Rho*mesh.aveF2CCV*jn + + ex,ez = Pfx*en, Pfz*en + + # get analytic solution + exa, eya, eza = EM.Analytics.FDEM.ElectricDipoleWholeSpace(XYZ, src_loc, sigmaback, freq,orientation='Z') + exa, eya, eza = Utils.mkvc(exa,2), Utils.mkvc(eya,2), Utils.mkvc(eza,2) + + print ' ex:', np.linalg.norm(exa), np.linalg.norm(ex), np.linalg.norm(exa-ex) + print ' ez:', np.linalg.norm(eza), np.linalg.norm(ez), np.linalg.norm(eza-ez) + + if plotIt: + if plotit: + 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() + + self.assertTrue(np.linalg.norm(exa-ex)/np.linalg.norm(exa) < tol_Edipole) + self.assertTrue(np.linalg.norm(eza-ez)/np.linalg.norm(eza) < tol_Edipole) + + + if __name__ == '__main__': unittest.main() From 075cea488fc6ecac38b9d15e3e441912165f9381 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Thu, 1 Oct 2015 23:52:49 -0700 Subject: [PATCH 4/6] test and see if we can run travis on updated infrastructure --- .travis.yml | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/.travis.yml b/.travis.yml index c64ba7f2..ee0b416f 100644 --- a/.travis.yml +++ b/.travis.yml @@ -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: From 8a358506f122c19be3066ce2ac57ca4627b8597d Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sat, 3 Oct 2015 14:06:06 -0700 Subject: [PATCH 5/6] add test for mag dipole in wholespace --- simpegEM/Analytics/FDEM.py | 2 +- simpegEM/Tests/test_FDEM_analytics.py | 121 ++++++++++++++++++-------- 2 files changed, 88 insertions(+), 35 deletions(-) diff --git a/simpegEM/Analytics/FDEM.py b/simpegEM/Analytics/FDEM.py index 9abb0a15..f1a0233e 100644 --- a/simpegEM/Analytics/FDEM.py +++ b/simpegEM/Analytics/FDEM.py @@ -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. diff --git a/simpegEM/Tests/test_FDEM_analytics.py b/simpegEM/Tests/test_FDEM_analytics.py index 21aa9fda..7d1128d5 100644 --- a/simpegEM/Tests/test_FDEM_analytics.py +++ b/simpegEM/Tests/test_FDEM_analytics.py @@ -4,7 +4,7 @@ import simpegEM as EM from scipy.constants import mu_0 plotIt = False -tol_Edipole = 1e-2 +tol_EBdipole = 1e-2 if plotIt: import matplotlib.pylab @@ -84,8 +84,8 @@ class FDEM_analyticTests(unittest.TestCase): self.assertTrue(passed) - def test_CylMeshEDipole(self): - print 'Testing CylMesh E Dipole in a wholespace- Analytic: J-formulation' + def test_CylMeshEBDipoles(self): + print 'Testing CylMesh Electric and Magnetic Dipoles in a wholespace- Analytic: J-formulation' sigmaback = 1. freq = 1. skdpth = 500./np.sqrt(sigmaback*freq) @@ -114,14 +114,25 @@ class FDEM_analyticTests(unittest.TestCase): 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) #+deg_p) # set survey + surveye = EM.FDEM.SurveyFDEM(de_p) + surveym = EM.FDEM.SurveyFDEM(dm_p) + mapping = [('sigma', Maps.IdentityMap(mesh))] - prbe = EM.FDEM.ProblemFDEM_j(mesh, mapping=mapping) + + 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 - fieldsBack = prbe.fields(np.r_[SigmaBack]) # Done + fieldsBackE = prbe.fields(np.r_[SigmaBack]) # Done + fieldsBackM = prbm.fields(np.r_[SigmaBack]) # Done + rlim = [20.,500.] lookAtTx = de_p @@ -135,52 +146,94 @@ class FDEM_analyticTests(unittest.TestCase): Zero = sp.csr_matrix(Pf.shape) Pfx,Pfz = sp.hstack([Pf,Zero]),sp.hstack([Zero,Pf]) - jn = fieldsBack[lookAtTx,'j'] + 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') exa, eya, eza = Utils.mkvc(exa,2), Utils.mkvc(eya,2), Utils.mkvc(eza,2) - print ' ex:', np.linalg.norm(exa), np.linalg.norm(ex), np.linalg.norm(exa-ex) - print ' ez:', np.linalg.norm(eza), np.linalg.norm(ez), np.linalg.norm(eza-ez) + bxa, bya, bza = EM.Analytics.FDEM.MagneticDipoleWholeSpace(XYZ, src_loc, sigmaback, freq,orientation='Z') + 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: - if plotit: - 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)') + # 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(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(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.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() + plt.tight_layout() - self.assertTrue(np.linalg.norm(exa-ex)/np.linalg.norm(exa) < tol_Edipole) - self.assertTrue(np.linalg.norm(eza-ez)/np.linalg.norm(eza) < tol_Edipole) + # 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) From b927a1daf2f23c284813fde154fb0844bf5f32f3 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Mon, 5 Oct 2015 07:42:58 -0700 Subject: [PATCH 6/6] - moved integration of src term in to definition of src (makes it easier to trigger on / off) - added mu to testing against analytics --- simpegEM/FDEM/FDEM.py | 54 ++++++++++++++------------- simpegEM/FDEM/FieldsFDEM.py | 6 +-- simpegEM/FDEM/SurveyFDEM.py | 15 ++++++-- simpegEM/Tests/test_FDEM_analytics.py | 12 +++--- 4 files changed, 50 insertions(+), 37 deletions(-) diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py index 38de9a66..4060e771 100644 --- a/simpegEM/FDEM/FDEM.py +++ b/simpegEM/FDEM/FDEM.py @@ -246,18 +246,20 @@ class ProblemFDEM_e(BaseFDEMProblem): """ S_m, S_e = self.getSourceTerm(freq) - Me = self.Me + # if self.survey. + # Me = self.Me C = self.mesh.edgeCurl MfMui = self.MfMui - RHS = C.T * (MfMui * S_m) -1j * omega(freq) * Me * 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 + # Me = self.Me S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint) if adjoint: @@ -265,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) * Me.T * 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) * Me.T * 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) * Me * 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) * Me * S_eDerivv + return -1j * omega(freq) * S_eDerivv else: return None @@ -363,9 +365,9 @@ class ProblemFDEM_b(BaseFDEMProblem): S_m, S_e = self.getSourceTerm(freq) C = self.mesh.edgeCurl MeSigmaI = self.MeSigmaI - Me = self.Me + # Me = self.Me - RHS = S_m + C * ( MeSigmaI * Me * S_e ) + RHS = S_m + C * ( MeSigmaI * S_e ) if self._makeASymmetric is True: MfMui = self.MfMui @@ -377,13 +379,13 @@ class ProblemFDEM_b(BaseFDEMProblem): C = self.mesh.edgeCurl S_m, S_e = src.eval(self) MfMui = self.MfMui - Me = self.Me + # Me = self.Me if self._makeASymmetric and adjoint: v = self.MfMui * v if S_e is not None: - MeSigmaIDeriv = self.MeSigmaIDeriv(Utils.mkvc( Me * S_e)) + MeSigmaIDeriv = self.MeSigmaIDeriv(S_e) if not adjoint: RHSderiv = C * (MeSigmaIDeriv * v) elif adjoint: @@ -395,16 +397,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 * (Me * S_eDeriv)) + SrcDeriv = S_mDeriv + C * (self.MeSigmaI * S_eDeriv) elif adjoint: - SrcDeriv = S_mDeriv + Me.T * (Self.MeSigmaI.T * ( C.T * S_eDeriv)) + SrcDeriv = S_mDeriv + 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 * (Me * S_eDeriv)) + SrcDeriv = C * (self.MeSigmaI * S_eDeriv) elif adjoint: - SrcDeriv = Me.T * (self.MeSigmaI.T * ( C.T * S_eDeriv)) + SrcDeriv = self.MeSigmaI.T * ( C.T * S_eDeriv) else: SrcDeriv = None @@ -512,10 +514,10 @@ class ProblemFDEM_j(BaseFDEMProblem): S_m, S_e = self.getSourceTerm(freq) C = self.mesh.edgeCurl MeMuI = self.MeMuI - Me = self.Me + # Me = self.Me - RHS = C * (MeMuI * (Me * 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 @@ -525,7 +527,7 @@ class ProblemFDEM_j(BaseFDEMProblem): def getRHSDeriv_m(self, src, v, adjoint=False): C = self.mesh.edgeCurl MeMuI = self.MeMuI - Me = self.Me + # Me = self.Me S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint) if adjoint: @@ -535,9 +537,9 @@ 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 Me.T * S_mDerivv - 1j * omega(freq) * S_eDerivv + return S_mDerivv - 1j * omega(freq) * S_eDerivv elif S_mDerivv is not None: - return Me.T * S_mDerivv + return S_mDerivv elif S_eDerivv is not None: return - 1j * omega(freq) * S_eDerivv else: @@ -546,9 +548,9 @@ class ProblemFDEM_j(BaseFDEMProblem): 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 * (Me * S_mDerivv)) - 1j * omega(freq) * S_eDerivv + RHSDeriv = C * (MeMuI * S_mDerivv) - 1j * omega(freq) * S_eDerivv elif S_mDerivv is not None: - RHSDeriv = C * (MeMuI * (Me * S_mDerivv)) + RHSDeriv = C * (MeMuI * S_mDerivv) elif S_eDerivv is not None: RHSDeriv = - 1j * omega(freq) * S_eDerivv else: @@ -626,9 +628,9 @@ class ProblemFDEM_h(BaseFDEMProblem): S_m, S_e = self.getSourceTerm(freq) C = self.mesh.edgeCurl MfRho = self.MfRho - Me = self.Me + # Me = self.Me - RHS = Me * S_m + C.T * ( MfRho * S_e ) + RHS = S_m + C.T * ( MfRho * S_e ) return RHS @@ -657,9 +659,9 @@ class ProblemFDEM_h(BaseFDEMProblem): if S_mDeriv is not None: if RHSDeriv is not None: - RHSDeriv += Me * S_mDeriv(v) + RHSDeriv += S_mDeriv(v) else: - RHSDeriv = Me * S_mDeriv(v) + RHSDeriv = S_mDeriv(v) if S_eDeriv is not None: if RHSDeriv is not None: RHSDeriv += C.T * (MfRho * S_e) diff --git a/simpegEM/FDEM/FieldsFDEM.py b/simpegEM/FDEM/FieldsFDEM.py index 1c57c1ca..21083264 100644 --- a/simpegEM/FDEM/FieldsFDEM.py +++ b/simpegEM/FDEM/FieldsFDEM.py @@ -145,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 * (self._Me * S_e) + e[:,i] += -self._MeSigmaI * S_e return e def _eSecondaryDeriv_u(self, src, v, adjoint=False): @@ -175,7 +175,7 @@ class FieldsFDEM_b(FieldsFDEM): Se_Deriv = S_eDeriv(v) if Se_Deriv is not None: - de_dm += -self._MeSigmaI * (self._Me * Se_Deriv) + de_dm += -self._MeSigmaI * Se_Deriv return de_dm @@ -247,7 +247,7 @@ class FieldsFDEM_j(FieldsFDEM): 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)) * self._MeMuI * (self._Me * S_m) + h[:,i] += 1./(1j*omega(src.freq)) * self._MeMuI * (S_m) return h def _hSecondaryDeriv_u(self, src, v, adjoint=False): diff --git a/simpegEM/FDEM/SurveyFDEM.py b/simpegEM/FDEM/SurveyFDEM.py index b47266cf..f65c5302 100644 --- a/simpegEM/FDEM/SurveyFDEM.py +++ b/simpegEM/FDEM/SurveyFDEM.py @@ -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): diff --git a/simpegEM/Tests/test_FDEM_analytics.py b/simpegEM/Tests/test_FDEM_analytics.py index 7d1128d5..ad643f50 100644 --- a/simpegEM/Tests/test_FDEM_analytics.py +++ b/simpegEM/Tests/test_FDEM_analytics.py @@ -87,6 +87,7 @@ class FDEM_analyticTests(unittest.TestCase): 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) @@ -104,6 +105,7 @@ class FDEM_analyticTests(unittest.TestCase): 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 @@ -121,7 +123,7 @@ class FDEM_analyticTests(unittest.TestCase): surveye = EM.FDEM.SurveyFDEM(de_p) surveym = EM.FDEM.SurveyFDEM(dm_p) - mapping = [('sigma', Maps.IdentityMap(mesh))] + 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) @@ -130,8 +132,8 @@ class FDEM_analyticTests(unittest.TestCase): prbm.pair(surveym) # solve - fieldsBackE = prbe.fields(np.r_[SigmaBack]) # Done - fieldsBackM = prbm.fields(np.r_[SigmaBack]) # Done + fieldsBackE = prbe.fields(np.r_[SigmaBack, MuBack]) # Done + fieldsBackM = prbm.fields(np.r_[SigmaBack, MuBack]) # Done rlim = [20.,500.] @@ -159,10 +161,10 @@ class FDEM_analyticTests(unittest.TestCase): bx,bz = Pfx*bn, Pfz*bn # get analytic solution - exa, eya, eza = EM.Analytics.FDEM.ElectricDipoleWholeSpace(XYZ, src_loc, sigmaback, freq,orientation='Z') + 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') + 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'