From b927a1daf2f23c284813fde154fb0844bf5f32f3 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Mon, 5 Oct 2015 07:42:58 -0700 Subject: [PATCH] - 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'