wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-07-05 03:04:17 +08:00
Code Issues Packages Projects Releases Wiki Activity

- moved integration of src term in to definition of src (makes it easier to trigger on / off)

Browse Source 
- added mu to testing against analytics
This commit is contained in:
Lindsey Heagy
2015-10-05 07:42:58 -07:00
parent 8a358506f1
commit b927a1daf2
4 changed files with 50 additions and 37 deletions
Download Patch File Download Diff File Expand all files Collapse all files
simpegEM/FDEM/FDEM.py
+28 -26
View File
@@ -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)
simpegEM/FDEM/FieldsFDEM.py
+3 -3
View File
@@ -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):
simpegEM/FDEM/SurveyFDEM.py
+12 -3
View File
@@ -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):
simpegEM/Tests/test_FDEM_analytics.py
+7 -5
View File
@@ -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'