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

Merge branch 'master' of https://github.com/simpeg/simpeg into dcip/dev

Browse Source
This commit is contained in:
seogi_macbook
2016-04-05 21:49:38 -07:00
parent fdc081970e aad596a8cc
commit a31319b46e
20 changed files with 1606 additions and 698 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/EM/FDEM/FDEM.py
+48 -64
View File
@@ -54,8 +54,7 @@ class BaseFDEMProblem(BaseEMProblem):
Ainv = self.Solver(A, **self.solverOpts)
sol = Ainv * rhs
Srcs = self.survey.getSrcByFreq(freq)
ftype = self._fieldType + 'Solution'
F[Srcs, ftype] = sol
F[Srcs, self._solutionType] = sol
Ainv.clean()
return F
@@ -78,30 +77,19 @@ class BaseFDEMProblem(BaseEMProblem):
Jv = self.dataPair(self.survey)
for freq in self.survey.freqs:
A = self.getA(freq) #
A = self.getA(freq)
Ainv = self.Solver(A, **self.solverOpts)
for src in self.survey.getSrcByFreq(freq):
ftype = self._fieldType + 'Solution'
u_src = u[src, ftype]
dA_dm = self.getADeriv_m(freq, u_src, v)
dRHS_dm = self.getRHSDeriv_m(freq, src, v)
du_dm = Ainv * ( - dA_dm + dRHS_dm )
u_src = u[src, self._solutionType]
dA_dm_v = self.getADeriv(freq, u_src, v)
dRHS_dm_v = self.getRHSDeriv(freq, src, v)
du_dm_v = Ainv * ( - dA_dm_v + dRHS_dm_v )
for rx in src.rxList:
df_duFun = getattr(u, '_%sDeriv_u'%rx.projField, None)
df_dudu_dm = df_duFun(src, du_dm, adjoint=False)
df_dmFun = getattr(u, '_%sDeriv_m'%rx.projField, None)
df_dm = df_dmFun(src, v, adjoint=False)
Df_Dm = np.array(df_dudu_dm + df_dm,dtype=complex)
P = lambda v: rx.evalDeriv(src, self.mesh, u, v) # wrt u, also have wrt m
Jv[src, rx] = P(Df_Dm)
df_dmFun = getattr(u, '_%sDeriv'%rx.projField, None)
df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False)
Jv[src, rx] = rx.evalDeriv(src, self.mesh, u, df_dm_v)
Ainv.clean()
return Utils.mkvc(Jv)
@@ -132,32 +120,28 @@ class BaseFDEMProblem(BaseEMProblem):
ATinv = self.Solver(AT, **self.solverOpts)
for src in self.survey.getSrcByFreq(freq):
ftype = self._fieldType + 'Solution'
u_src = u[src, ftype]
u_src = u[src, self._solutionType]
for rx in src.rxList:
PTv = rx.evalDeriv(src, self.mesh, u, v[src, rx], adjoint=True) # wrt u, need possibility wrt m
df_duTFun = getattr(u, '_%sDeriv_u'%rx.projField, None)
df_duT = df_duTFun(src, PTv, adjoint=True)
df_duTFun = getattr(u, '_%sDeriv'%rx.projField, None)
df_duT, df_dmT = df_duTFun(src, None, PTv, adjoint=True)
ATinvdf_duT = ATinv * df_duT
dA_dmT = self.getADeriv_m(freq, u_src, ATinvdf_duT, adjoint=True)
dRHS_dmT = self.getRHSDeriv_m(freq,src, ATinvdf_duT, adjoint=True)
dA_dmT = self.getADeriv(freq, u_src, ATinvdf_duT, adjoint=True)
dRHS_dmT = self.getRHSDeriv(freq, src, ATinvdf_duT, adjoint=True)
du_dmT = -dA_dmT + dRHS_dmT
df_dmFun = getattr(u, '_%sDeriv_m'%rx.projField, None)
dfT_dm = df_dmFun(src, PTv, adjoint=True)
df_dmT = df_dmT + du_dmT
du_dmT += dfT_dm
# TODO: this should be taken care of by the reciever
# TODO: this should be taken care of by the reciever?
real_or_imag = rx.projComp
if real_or_imag is 'real':
Jtv += np.array(du_dmT,dtype=complex).real
Jtv += np.array(df_dmT, dtype=complex).real
elif real_or_imag is 'imag':
Jtv += - np.array(du_dmT,dtype=complex).real
Jtv += - np.array(df_dmT, dtype=complex).real
else:
raise Exception('Must be real or imag')
@@ -174,10 +158,10 @@ class BaseFDEMProblem(BaseEMProblem):
:return: S_m, S_e (nE or nF, nSrc)
"""
Srcs = self.survey.getSrcByFreq(freq)
if self._eqLocs is 'FE':
if self._formulation is 'EB':
S_m = np.zeros((self.mesh.nF,len(Srcs)), dtype=complex)
S_e = np.zeros((self.mesh.nE,len(Srcs)), dtype=complex)
elif self._eqLocs is 'EF':
elif self._formulation is 'HJ':
S_m = np.zeros((self.mesh.nE,len(Srcs)), dtype=complex)
S_e = np.zeros((self.mesh.nF,len(Srcs)), dtype=complex)
@@ -213,9 +197,9 @@ class Problem_e(BaseFDEMProblem):
:param SimPEG.Mesh mesh: mesh
"""
_fieldType = 'e'
_eqLocs = 'FE'
fieldsPair = Fields_e
_solutionType = 'eSolution'
_formulation = 'EB'
fieldsPair = Fields_e
def __init__(self, mesh, **kwargs):
BaseFDEMProblem.__init__(self, mesh, **kwargs)
@@ -239,7 +223,7 @@ class Problem_e(BaseFDEMProblem):
return C.T*MfMui*C + 1j*omega(freq)*MeSigma
def getADeriv_m(self, freq, u, v, adjoint=False):
def getADeriv(self, freq, u, v, adjoint=False):
"""
Product of the derivative of our system matrix with respect to the model and a vector
@@ -280,7 +264,7 @@ class Problem_e(BaseFDEMProblem):
return C.T * (MfMui * S_m) -1j * omega(freq) * S_e
def getRHSDeriv_m(self, freq, src, v, adjoint=False):
def getRHSDeriv(self, freq, src, v, adjoint=False):
"""
Derivative of the right hand side with respect to the model
@@ -324,9 +308,9 @@ class Problem_b(BaseFDEMProblem):
:param SimPEG.Mesh mesh: mesh
"""
_fieldType = 'b'
_eqLocs = 'FE'
fieldsPair = Fields_b
_solutionType = 'bSolution'
_formulation = 'EB'
fieldsPair = Fields_b
def __init__(self, mesh, **kwargs):
BaseFDEMProblem.__init__(self, mesh, **kwargs)
@@ -354,7 +338,7 @@ class Problem_b(BaseFDEMProblem):
return MfMui.T*A
return A
def getADeriv_m(self, freq, u, v, adjoint=False):
def getADeriv(self, freq, u, v, adjoint=False):
"""
Product of the derivative of our system matrix with respect to the model and a vector
@@ -411,7 +395,7 @@ class Problem_b(BaseFDEMProblem):
return RHS
def getRHSDeriv_m(self, freq, src, v, adjoint=False):
def getRHSDeriv(self, freq, src, v, adjoint=False):
"""
Derivative of the right hand side with respect to the model
@@ -472,9 +456,9 @@ class Problem_j(BaseFDEMProblem):
:param SimPEG.Mesh mesh: mesh
"""
_fieldType = 'j'
_eqLocs = 'EF'
fieldsPair = Fields_j
_solutionType = 'jSolution'
_formulation = 'HJ'
fieldsPair = Fields_j
def __init__(self, mesh, **kwargs):
BaseFDEMProblem.__init__(self, mesh, **kwargs)
@@ -503,7 +487,7 @@ class Problem_j(BaseFDEMProblem):
return A
def getADeriv_m(self, freq, u, v, adjoint=False):
def getADeriv(self, freq, u, v, adjoint=False):
"""
Product of the derivative of our system matrix with respect to the model and a vector
@@ -524,16 +508,16 @@ class Problem_j(BaseFDEMProblem):
MeMuI = self.MeMuI
MfRho = self.MfRho
C = self.mesh.edgeCurl
MfRhoDeriv_m = self.MfRhoDeriv(u)
MfRhoDeriv = self.MfRhoDeriv(u)
if adjoint:
if self._makeASymmetric is True:
v = MfRho * v
return MfRhoDeriv_m.T * (C * (MeMuI.T * (C.T * v)))
return MfRhoDeriv.T * (C * (MeMuI.T * (C.T * v)))
if self._makeASymmetric is True:
return MfRho.T * (C * ( MeMuI * (C.T * (MfRhoDeriv_m * v) )))
return C * (MeMuI * (C.T * (MfRhoDeriv_m * v)))
return MfRho.T * (C * ( MeMuI * (C.T * (MfRhoDeriv * v) )))
return C * (MeMuI * (C.T * (MfRhoDeriv * v)))
def getRHS(self, freq):
@@ -560,7 +544,7 @@ class Problem_j(BaseFDEMProblem):
return RHS
def getRHSDeriv_m(self, freq, src, v, adjoint=False):
def getRHSDeriv(self, freq, src, v, adjoint=False):
"""
Derivative of the right hand side with respect to the model
@@ -610,9 +594,9 @@ class Problem_h(BaseFDEMProblem):
:param SimPEG.Mesh mesh: mesh
"""
_fieldType = 'h'
_eqLocs = 'EF'
fieldsPair = Fields_h
_solutionType = 'hSolution'
_formulation = 'HJ'
fieldsPair = Fields_h
def __init__(self, mesh, **kwargs):
BaseFDEMProblem.__init__(self, mesh, **kwargs)
@@ -635,7 +619,7 @@ class Problem_h(BaseFDEMProblem):
return C.T * (MfRho * C) + 1j*omega(freq)*MeMu
def getADeriv_m(self, freq, u, v, adjoint=False):
def getADeriv(self, freq, u, v, adjoint=False):
"""
Product of the derivative of our system matrix with respect to the model and a vector
@@ -652,11 +636,11 @@ class Problem_h(BaseFDEMProblem):
MeMu = self.MeMu
C = self.mesh.edgeCurl
MfRhoDeriv_m = self.MfRhoDeriv(C*u)
MfRhoDeriv = self.MfRhoDeriv(C*u)
if adjoint:
return MfRhoDeriv_m.T * (C * v)
return C.T * (MfRhoDeriv_m * v)
return MfRhoDeriv.T * (C * v)
return C.T * (MfRhoDeriv * v)
def getRHS(self, freq):
"""
@@ -677,7 +661,7 @@ class Problem_h(BaseFDEMProblem):
return S_m + C.T * ( MfRho * S_e )
def getRHSDeriv_m(self, freq, src, v, adjoint=False):
def getRHSDeriv(self, freq, src, v, adjoint=False):
"""
Derivative of the right hand side with respect to the model
SimPEG/EM/FDEM/FieldsFDEM.py
+698 -313
View File
File diff suppressed because it is too large Load Diff
SimPEG/EM/FDEM/SrcFDEM.py
+117 -79
View File
@@ -1,7 +1,7 @@
from SimPEG import Survey, Problem, Utils, np, sp
from scipy.constants import mu_0
from SimPEG.EM.Utils import *
from SimPEG.Utils import Zero
from SimPEG.Utils import Zero
class BaseSrc(Survey.BaseSrc):
"""
@@ -14,7 +14,7 @@ class BaseSrc(Survey.BaseSrc):
def eval(self, prob):
"""
Evaluate the source terms.
Evaluate the source terms.
- :math:`S_m` : magnetic source term
- :math:`S_e` : electric source term
@@ -36,12 +36,12 @@ class BaseSrc(Survey.BaseSrc):
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
:rtype: (numpy.ndarray, numpy.ndarray)
:return: tuple with magnetic source term and electric source term derivatives times a vector
:return: tuple with magnetic source term and electric source term derivatives times a vector
"""
if v is not None:
return self.S_mDeriv(prob,v,adjoint), self.S_eDeriv(prob,v,adjoint)
if v is not None:
return self.S_mDeriv(prob, v, adjoint), self.S_eDeriv(prob, v, adjoint)
else:
return lambda v: self.S_mDeriv(prob,v,adjoint), lambda v: self.S_eDeriv(prob,v,adjoint)
return lambda v: self.S_mDeriv(prob, v, adjoint), lambda v: self.S_eDeriv(prob, v, adjoint)
def bPrimary(self, prob):
"""
@@ -49,7 +49,7 @@ class BaseSrc(Survey.BaseSrc):
:param Problem prob: FDEM Problem
:rtype: numpy.ndarray
:return: primary magnetic flux density
:return: primary magnetic flux density
"""
return Zero()
@@ -59,7 +59,7 @@ class BaseSrc(Survey.BaseSrc):
:param Problem prob: FDEM Problem
:rtype: numpy.ndarray
:return: primary magnetic field
:return: primary magnetic field
"""
return Zero()
@@ -69,7 +69,7 @@ class BaseSrc(Survey.BaseSrc):
:param Problem prob: FDEM Problem
:rtype: numpy.ndarray
:return: primary electric field
:return: primary electric field
"""
return Zero()
@@ -79,13 +79,13 @@ class BaseSrc(Survey.BaseSrc):
:param Problem prob: FDEM Problem
:rtype: numpy.ndarray
:return: primary current density
:return: primary current density
"""
return Zero()
def S_m(self, prob):
"""
Magnetic source term
Magnetic source term
:param Problem prob: FDEM Problem
:rtype: numpy.ndarray
@@ -95,7 +95,7 @@ class BaseSrc(Survey.BaseSrc):
def S_e(self, prob):
"""
Electric source term
Electric source term
:param Problem prob: FDEM Problem
:rtype: numpy.ndarray
@@ -136,15 +136,26 @@ class RawVec_e(BaseSrc):
:param list rxList: receiver list
:param float freq: frequency
:param numpy.array S_e: electric source term
:param bool integrate: Integrate the source term (multiply by Me) [True]
"""
def __init__(self, rxList, freq, S_e): #, ePrimary=None, bPrimary=None, hPrimary=None, jPrimary=None):
self._S_e = np.array(S_e,dtype=complex)
def __init__(self, rxList, freq, S_e, integrate=True): #, ePrimary=None, bPrimary=None, hPrimary=None, jPrimary=None):
self._S_e = np.array(S_e, dtype=complex)
self.freq = float(freq)
self.integrate = integrate
BaseSrc.__init__(self, rxList)
def S_e(self, prob):
"""
Electric source term
:param Problem prob: FDEM Problem
:rtype: numpy.ndarray
:return: electric source term on mesh
"""
if prob._formulation is 'EB' and self.integrate is True:
return prob.Me * self._S_e
return self._S_e
@@ -155,10 +166,11 @@ class RawVec_m(BaseSrc):
:param float freq: frequency
:param rxList: receiver list
:param numpy.array S_m: magnetic source term
:param bool integrate: Integrate the source term (multiply by Me) [True]
"""
def __init__(self, rxList, freq, S_m, integrate = True): #ePrimary=Zero(), bPrimary=Zero(), hPrimary=Zero(), jPrimary=Zero()):
self._S_m = np.array(S_m,dtype=complex)
def __init__(self, rxList, freq, S_m, integrate=True): #ePrimary=Zero(), bPrimary=Zero(), hPrimary=Zero(), jPrimary=Zero()):
self._S_m = np.array(S_m, dtype=complex)
self.freq = float(freq)
self.integrate = integrate
@@ -166,12 +178,14 @@ class RawVec_m(BaseSrc):
def S_m(self, prob):
"""
Magnetic source term
Magnetic source term
:param Problem prob: FDEM Problem
:rtype: numpy.ndarray
:return: magnetic source term on mesh
"""
if prob._formulation is 'HJ' and self.integrate is True:
return prob.Me * self._S_m
return self._S_m
@@ -183,36 +197,51 @@ class RawVec(BaseSrc):
:param float freq: frequency
:param numpy.array S_m: magnetic source term
:param numpy.array S_e: electric source term
:param bool integrate: Integrate the source term (multiply by Me) [True]
"""
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)
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
BaseSrc.__init__(self, rxList)
def S_m(self, prob):
if prob._eqLocs is 'EF' and self.integrate is True:
"""
Magnetic source term
:param Problem prob: FDEM Problem
:rtype: numpy.ndarray
:return: magnetic source term on mesh
"""
if prob._formulation is 'HJ' 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:
"""
Electric source term
:param Problem prob: FDEM Problem
:rtype: numpy.ndarray
:return: electric source term on mesh
"""
if prob._formulation is 'EB' and self.integrate is True:
return prob.Me * self._S_e
return self._S_e
class MagDipole(BaseSrc):
"""
"""
Point magnetic dipole source calculated by taking the curl of a magnetic
vector potential. By taking the discrete curl, we ensure that the magnetic
flux density is divergence free (no magnetic monopoles!).
flux density is divergence free (no magnetic monopoles!).
This approach uses a primary-secondary in frequency. Here we show the
derivation for E-B formulation noting that similar steps are followed for
the H-J formulation.
.. math::
.. 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}}
@@ -225,17 +254,17 @@ class MagDipole(BaseSrc):
and define a zero-frequency primary problem, noting that the source is
generated by a divergence free electric current
.. math::
.. math::
\mathbf{C} \mathbf{e^P} = \mathbf{s_m^P} = 0 \\\\
{\mathbf{C}^T \mathbf{{M_{\mu^{-1}}^f}^P} \mathbf{b^P} - \mathbf{M_{\sigma}^e} \mathbf{e^P} = \mathbf{M^e} \mathbf{s_e^P}}
Since :math:`\mathbf{e^P}` is curl-free, divergence-free, we assume that there is no constant field background, the :math:`\mathbf{e^P} = 0`, so our primary problem is
Since :math:`\mathbf{e^P}` is curl-free, divergence-free, we assume that there is no constant field background, the :math:`\mathbf{e^P} = 0`, so our primary problem is
.. math::
.. math::
\mathbf{e^P} = 0 \\\\
{\mathbf{C}^T \mathbf{{M_{\mu^{-1}}^f}^P} \mathbf{b^P} = \mathbf{s_e^P}}
Our secondary problem is then
Our secondary problem is then
.. math::
\mathbf{C} \mathbf{e^S} + i \omega \mathbf{b^S} = - i \omega \mathbf{b^P} \\\\
@@ -245,15 +274,15 @@ class MagDipole(BaseSrc):
:param float freq: frequency
:param numpy.ndarray loc: source location (ie: :code:`np.r_[xloc,yloc,zloc]`)
:param string orientation: 'X', 'Y', 'Z'
:param float moment: magnetic dipole moment
:param float mu: background magnetic permeability
:param float moment: magnetic dipole moment
:param float mu: background magnetic permeability
"""
#TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that
def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu = mu_0):
def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu=mu_0):
self.freq = float(freq)
self.loc = loc
self.orientation = orientation
assert orientation in ['X','Y','Z'], "Orientation (right now) doesn't actually do anything! The methods in SrcUtils should take care of this..."
self.moment = moment
self.mu = mu
self.integrate = False
@@ -265,17 +294,17 @@ class MagDipole(BaseSrc):
:param Problem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
:return: primary magnetic field
"""
eqLocs = prob._eqLocs
formulation = prob._formulation
if eqLocs is 'FE':
if formulation is 'EB':
gridX = prob.mesh.gridEx
gridY = prob.mesh.gridEy
gridZ = prob.mesh.gridEz
C = prob.mesh.edgeCurl
elif eqLocs is 'EF':
elif formulation is 'HJ':
gridX = prob.mesh.gridFx
gridY = prob.mesh.gridFy
gridZ = prob.mesh.gridFz
@@ -303,10 +332,10 @@ class MagDipole(BaseSrc):
:param Problem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
:return: primary magnetic field
"""
b = self.bPrimary(prob)
return h_from_b(prob,b)
return 1./self.mu * b
def S_m(self, prob):
"""
@@ -314,10 +343,12 @@ class MagDipole(BaseSrc):
:param Problem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
:return: primary magnetic field
"""
b_p = self.bPrimary(prob)
if prob._formulation is 'HJ':
b_p = prob.Me * b_p
return -1j*omega(self.freq)*b_p
def S_e(self, prob):
@@ -326,21 +357,21 @@ class MagDipole(BaseSrc):
:param Problem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
:return: primary magnetic field
"""
if all(np.r_[self.mu] == np.r_[prob.curModel.mu]):
return Zero()
else:
eqLocs = prob._eqLocs
formulation = prob._formulation
if eqLocs is 'FE':
if formulation is 'EB':
mui_s = prob.curModel.mui - 1./self.mu
MMui_s = prob.mesh.getFaceInnerProduct(mui_s)
C = prob.mesh.edgeCurl
elif eqLocs is 'EF':
elif formulation is 'HJ':
mu_s = prob.curModel.mu - self.mu
MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True)
MMui_s = prob.mesh.getEdgeInnerProduct(mu_s, invMat=True)
C = prob.mesh.edgeCurl.T
return -C.T * (MMui_s * self.bPrimary(prob))
@@ -353,21 +384,20 @@ class MagDipole_Bfield(BaseSrc):
fields from a magnetic dipole. No discrete curl is taken, so the magnetic
flux density may not be strictly divergence free.
This approach uses a primary-secondary in frequency in the same fashion as the MagDipole.
This approach uses a primary-secondary in frequency in the same fashion as the MagDipole.
:param list rxList: receiver list
:param float freq: frequency
:param numpy.ndarray loc: source location (ie: :code:`np.r_[xloc,yloc,zloc]`)
:param string orientation: 'X', 'Y', 'Z'
:param float moment: magnetic dipole moment
:param float mu: background magnetic permeability
:param float moment: magnetic dipole moment
:param float mu: background magnetic permeability
"""
#TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that
#TODO: neither does moment
def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu = mu_0):
self.freq = float(freq)
self.loc = loc
assert orientation in ['X','Y','Z'], "Orientation (right now) doesn't actually do anything! The methods in SrcUtils should take care of this..."
self.orientation = orientation
self.moment = moment
self.mu = mu
@@ -379,18 +409,18 @@ class MagDipole_Bfield(BaseSrc):
:param Problem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
:return: primary magnetic field
"""
eqLocs = prob._eqLocs
formulation = prob._formulation
if eqLocs is 'FE':
if formulation is 'EB':
gridX = prob.mesh.gridFx
gridY = prob.mesh.gridFy
gridZ = prob.mesh.gridFz
C = prob.mesh.edgeCurl
elif eqLocs is 'EF':
elif formulation is 'HJ':
gridX = prob.mesh.gridEx
gridY = prob.mesh.gridEy
gridZ = prob.mesh.gridEz
@@ -418,10 +448,10 @@ class MagDipole_Bfield(BaseSrc):
:param Problem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
:return: primary magnetic field
"""
b = self.bPrimary(prob)
return h_from_b(prob, b)
return 1/self.mu * b
def S_m(self, prob):
"""
@@ -429,9 +459,11 @@ class MagDipole_Bfield(BaseSrc):
:param Problem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
:return: primary magnetic field
"""
b = self.bPrimary(prob)
if prob._formulation is 'HJ':
b = prob.Me * b
return -1j*omega(self.freq)*b
def S_e(self, prob):
@@ -440,20 +472,20 @@ class MagDipole_Bfield(BaseSrc):
:param Problem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
:return: primary magnetic field
"""
if all(np.r_[self.mu] == np.r_[prob.curModel.mu]):
return Zero()
else:
eqLocs = prob._eqLocs
formulation = prob._formulation
if eqLocs is 'FE':
if formulation is 'EB':
mui_s = prob.curModel.mui - 1./self.mu
MMui_s = prob.mesh.getFaceInnerProduct(mui_s)
C = prob.mesh.edgeCurl
elif eqLocs is 'EF':
elif formulation is 'HJ':
mu_s = prob.curModel.mu - self.mu
MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True)
MMui_s = prob.mesh.getEdgeInnerProduct(mu_s, invMat=True)
C = prob.mesh.edgeCurl.T
return -C.T * (MMui_s * self.bPrimary(prob))
@@ -463,22 +495,22 @@ class CircularLoop(BaseSrc):
"""
Circular loop magnetic source calculated by taking the curl of a magnetic
vector potential. By taking the discrete curl, we ensure that the magnetic
flux density is divergence free (no magnetic monopoles!).
flux density is divergence free (no magnetic monopoles!).
This approach uses a primary-secondary in frequency in the same fashion as the MagDipole.
This approach uses a primary-secondary in frequency in the same fashion as the MagDipole.
:param list rxList: receiver list
:param float freq: frequency
:param numpy.ndarray loc: source location (ie: :code:`np.r_[xloc,yloc,zloc]`)
:param string orientation: 'X', 'Y', 'Z'
:param float moment: magnetic dipole moment
:param float mu: background magnetic permeability
:param float moment: magnetic dipole moment
:param float mu: background magnetic permeability
"""
#TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that
def __init__(self, rxList, freq, loc, orientation='Z', radius = 1., mu=mu_0):
def __init__(self, rxList, freq, loc, orientation='Z', radius=1., mu=mu_0):
self.freq = float(freq)
self.orientation = orientation
assert orientation in ['X','Y','Z'], "Orientation (right now) doesn't actually do anything! The methods in SrcUtils should take care of this..."
self.radius = radius
self.mu = mu
self.loc = loc
@@ -491,17 +523,17 @@ class CircularLoop(BaseSrc):
:param Problem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
:return: primary magnetic field
"""
eqLocs = prob._eqLocs
formulation = prob._formulation
if eqLocs is 'FE':
if formulation is 'EB':
gridX = prob.mesh.gridEx
gridY = prob.mesh.gridEy
gridZ = prob.mesh.gridEz
C = prob.mesh.edgeCurl
elif eqLocs is 'EF':
elif formulation is 'HJ':
gridX = prob.mesh.gridFx
gridY = prob.mesh.gridFy
gridZ = prob.mesh.gridFz
@@ -528,7 +560,7 @@ class CircularLoop(BaseSrc):
:param Problem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
:return: primary magnetic field
"""
b = self.bPrimary(prob)
return 1./self.mu*b
@@ -539,9 +571,11 @@ class CircularLoop(BaseSrc):
:param Problem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
:return: primary magnetic field
"""
b = self.bPrimary(prob)
if prob._formulation is 'HJ':
b = prob.Me * b
return -1j*omega(self.freq)*b
def S_e(self, prob):
@@ -550,22 +584,26 @@ class CircularLoop(BaseSrc):
:param Problem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
:return: primary magnetic field
"""
if all(np.r_[self.mu] == np.r_[prob.curModel.mu]):
return Zero()
else:
eqLocs = prob._eqLocs
formulation = prob._formulation
if eqLocs is 'FE':
if formulation is 'EB':
mui_s = prob.curModel.mui - 1./self.mu
MMui_s = prob.mesh.getFaceInnerProduct(mui_s)
C = prob.mesh.edgeCurl
elif eqLocs is 'EF':
elif formulation is 'HJ':
mu_s = prob.curModel.mu - self.mu
MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True)
MMui_s = prob.mesh.getEdgeInnerProduct(mu_s, invMat=True)
C = prob.mesh.edgeCurl.T
return -C.T * (MMui_s * self.bPrimary(prob))
SimPEG/EM/FDEM/SurveyFDEM.py
+43 -36
View File
@@ -3,6 +3,7 @@ from SimPEG.EM.Utils import *
from scipy.constants import mu_0
from SimPEG.Utils import Zero, Identity
import SrcFDEM as Src
from SimPEG import sp
####################################################
@@ -18,33 +19,33 @@ class Rx(SimPEG.Survey.BaseRx):
"""
knownRxTypes = {
'exr':['e', 'Ex', 'real'],
'eyr':['e', 'Ey', 'real'],
'ezr':['e', 'Ez', 'real'],
'exi':['e', 'Ex', 'imag'],
'eyi':['e', 'Ey', 'imag'],
'ezi':['e', 'Ez', 'imag'],
'exr':['e', 'x', 'real'],
'eyr':['e', 'y', 'real'],
'ezr':['e', 'z', 'real'],
'exi':['e', 'x', 'imag'],
'eyi':['e', 'y', 'imag'],
'ezi':['e', 'z', 'imag'],
'bxr':['b', 'Fx', 'real'],
'byr':['b', 'Fy', 'real'],
'bzr':['b', 'Fz', 'real'],
'bxi':['b', 'Fx', 'imag'],
'byi':['b', 'Fy', 'imag'],
'bzi':['b', 'Fz', 'imag'],
'bxr':['b', 'x', 'real'],
'byr':['b', 'y', 'real'],
'bzr':['b', 'z', 'real'],
'bxi':['b', 'x', 'imag'],
'byi':['b', 'y', 'imag'],
'bzi':['b', 'z', 'imag'],
'jxr':['j', 'Fx', 'real'],
'jyr':['j', 'Fy', 'real'],
'jzr':['j', 'Fz', 'real'],
'jxi':['j', 'Fx', 'imag'],
'jyi':['j', 'Fy', 'imag'],
'jzi':['j', 'Fz', 'imag'],
'jxr':['j', 'x', 'real'],
'jyr':['j', 'y', 'real'],
'jzr':['j', 'z', 'real'],
'jxi':['j', 'x', 'imag'],
'jyi':['j', 'y', 'imag'],
'jzi':['j', 'z', 'imag'],
'hxr':['h', 'Ex', 'real'],
'hyr':['h', 'Ey', 'real'],
'hzr':['h', 'Ez', 'real'],
'hxi':['h', 'Ex', 'imag'],
'hyi':['h', 'Ey', 'imag'],
'hzi':['h', 'Ez', 'imag'],
'hxr':['h', 'x', 'real'],
'hyr':['h', 'y', 'real'],
'hzr':['h', 'z', 'real'],
'hxi':['h', 'x', 'imag'],
'hyi':['h', 'y', 'imag'],
'hzi':['h', 'z', 'imag'],
}
radius = None
@@ -56,17 +57,16 @@ class Rx(SimPEG.Survey.BaseRx):
"""Field Type projection (e.g. e b ...)"""
return self.knownRxTypes[self.rxType][0]
@property
def projGLoc(self):
"""Grid Location projection (e.g. Ex Fy ...)"""
return self.knownRxTypes[self.rxType][1]
@property
def projComp(self):
"""Component projection (real/imag)"""
return self.knownRxTypes[self.rxType][2]
def eval(self, src, mesh, f):
def projGLoc(self, u):
"""Grid Location projection (e.g. Ex Fy ...)"""
return u._GLoc(self.rxType[0]) + self.knownRxTypes[self.rxType][1]
def eval(self, src, mesh, u):
"""
Project fields to recievers to get data.
@@ -76,24 +76,30 @@ class Rx(SimPEG.Survey.BaseRx):
:rtype: numpy.ndarray
:return: fields projected to recievers
"""
P = self.getP(mesh) # get interpolation to recievers
u_part_complex = f[src, self.projField]
real_or_imag = self.projComp # get the real or imag component
# projGLoc = u._GLoc(self.knownRxTypes[self.rxType][0])
# projGLoc += self.knownRxTypes[self.rxType][1]
P = self.getP(mesh, self.projGLoc(u))
u_part_complex = u[src, self.projField]
# get the real or imag component
real_or_imag = self.projComp
u_part = getattr(u_part_complex, real_or_imag)
return P*u_part
def evalDeriv(self, src, mesh, f, v, adjoint=False):
def evalDeriv(self, src, mesh, u, v, adjoint=False):
"""
Derivative of projected fields with respect to the inversion model times a vector.
:param Source src: FDEM source
:param Mesh mesh: mesh used
:param Fields f: fields object
:param Fields u: fields object
:param numpy.ndarray v: vector to multiply
:rtype: numpy.ndarray
:return: fields projected to recievers
"""
P = self.getP(mesh)
P = self.getP(mesh, self.projGLoc(u))
if not adjoint:
Pv_complex = P * v
@@ -185,3 +191,4 @@ class Survey(SimPEG.Survey.BaseSurvey):
def evalDeriv(self, u):
raise Exception('Use Receivers to project fields deriv.')
SimPEG/EM/Utils/EMUtils.py
-33
View File
@@ -13,37 +13,4 @@ def k(freq, sigma, mu=mu_0, eps=epsilon_0):
beta = w * np.sqrt( mu*eps/2 * ( np.sqrt(1. + (sigma / (eps*w))**2 ) - 1) )
return alp - 1j*beta
# Constitutive relations
def e_from_j(prob,j):
eqLocs = prob._eqLocs
if eqLocs is 'FE':
MSigmaI = prob.MeSigmaI
elif eqLocs is 'EF':
MSigmaI = prob.MfRho
return MSigmaI*j
def j_from_e(prob,e):
eqLocs = prob._eqLocs
if eqLocs is 'FE':
MSigma = prob.MeSigma
elif eqLocs is 'EF':
MSigma = prob.MfRhoI
return MSigma*e
def b_from_h(prob,h):
eqLocs = prob._eqLocs
if eqLocs is 'FE':
MMu = prob.MfMuiI
elif eqLocs is 'EF':
MMu = prob.MeMu
return MMu*h
def h_from_b(prob,b):
eqLocs = prob._eqLocs
if eqLocs is 'FE':
MMuI = prob.MfMui
elif eqLocs is 'EF':
MMuI = prob.MeMuI
return MMuI*b
SimPEG/EM/Utils/__init__.py
+1 -4
View File
@@ -1,5 +1,2 @@
# import Sources
# import Ana
# import Solver
from EMUtils import omega, e_from_j, j_from_e, b_from_h, h_from_b
from EMUtils import omega, k
from AnalyticUtils import MagneticDipoleFields, MagneticDipoleVectorPotential, MagneticLoopVectorPotential
SimPEG/EM/Utils/testingUtils.py
+64 -13
View File
@@ -4,19 +4,28 @@ from SimPEG import EM
import sys
from scipy.constants import mu_0
def getFDEMProblem(fdemType, comp, SrcList, freq, verbose=False):
cs = 5.
ncx, ncy, ncz = 6, 6, 6
npad = 3
FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order
CONDUCTIVITY = 1e1
MU = mu_0
freq = 5e-1
def getFDEMProblem(fdemType, comp, SrcList, freq, useMu=False, verbose=False):
cs = 10.
ncx, ncy, ncz = 0, 0, 0
npad = 8
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)]
mesh = Mesh.TensorMesh([hx,hy,hz],['C','C','C'])
mapping = Maps.ExpMap(mesh)
if useMu is True:
mapping = [('sigma', Maps.ExpMap(mesh)), ('mu', Maps.IdentityMap(mesh))]
else:
mapping = Maps.ExpMap(mesh)
x = np.array([np.linspace(-30,-15,3),np.linspace(15,30,3)]) #don't sample right by the source
XYZ = Utils.ndgrid(x,x,np.r_[0.])
x = np.array([np.linspace(-5.*cs,-2.*cs,3),np.linspace(5.*cs,2.*cs,3)]) + cs/4. #don't sample right by the source, slightly off alignment from either staggered grid
XYZ = Utils.ndgrid(x,x,np.linspace(-2.*cs,2.*cs,5))
Rx0 = EM.FDEM.Rx(XYZ, comp)
Src = []
@@ -32,15 +41,15 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, verbose=False):
if fdemType is 'e' or fdemType is 'b':
S_m = np.zeros(mesh.nF)
S_e = np.zeros(mesh.nE)
S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1.
S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1.
S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1e-3
S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1e-3
Src.append(EM.FDEM.Src.RawVec([Rx0], freq, S_m, S_e))
elif fdemType is 'h' or fdemType is 'j':
S_m = np.zeros(mesh.nE)
S_e = np.zeros(mesh.nF)
S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1.
S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1.
S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1e-3
S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1e-3
Src.append(EM.FDEM.Src.RawVec([Rx0], freq, S_m, S_e))
if verbose:
@@ -70,6 +79,48 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, verbose=False):
from pymatsolver import MumpsSolver
prb.Solver = MumpsSolver
except ImportError, e:
pass
prb.Solver = SolverLU
return prb
return prb
def crossCheckTest(SrcList, fdemType1, fdemType2, comp, addrandoms = False, useMu=False, TOL=1e-5, verbose=False):
l2norm = lambda r: np.sqrt(r.dot(r))
prb1 = getFDEMProblem(fdemType1, comp, SrcList, freq, useMu, verbose)
mesh = prb1.mesh
print 'Cross Checking Forward: %s, %s formulations - %s' % (fdemType1, fdemType2, comp)
logsig = np.log(np.ones(mesh.nC)*CONDUCTIVITY)
mu = np.ones(mesh.nC)*MU
if addrandoms is True:
logsig += np.random.randn(mesh.nC)*np.log(CONDUCTIVITY)*1e-1
mu += np.random.randn(mesh.nC)*MU*1e-1
if useMu is True:
m = np.r_[logsig, mu]
else:
m = logsig
survey1 = prb1.survey
d1 = survey1.dpred(m)
if verbose:
print ' Problem 1 solved'
prb2 = getFDEMProblem(fdemType2, comp, SrcList, freq, useMu, verbose)
survey2 = prb2.survey
d2 = survey2.dpred(m)
if verbose:
print ' Problem 2 solved'
r = d2-d1
l2r = l2norm(r)
tol = np.max([TOL*(10**int(np.log10(0.5* (l2norm(d1) + l2norm(d2)) ))),FLR])
print l2norm(d1), l2norm(d2), l2r , tol, l2r < tol
return l2r < tol
SimPEG/Examples/EM_FDEM_1D_Inversion.py
+1 -2
View File
@@ -48,8 +48,7 @@ def run(plotIt=True):
freqs = np.logspace(1,3,10)
srcLoc = np.array([0., 0., 10.])
srcList = []
[srcList.append(EM.FDEM.Src.MagDipole([bzi],freq, srcLoc,orientation='Z')) for freq in freqs]
srcList = [EM.FDEM.Src.MagDipole([bzi],freq, srcLoc,orientation='Z') for freq in freqs]
survey = EM.FDEM.Survey(srcList)
prb = EM.FDEM.Problem_b(mesh, mapping=mapping)
SimPEG/Mesh/TensorMesh.py
+13
View File
@@ -234,6 +234,9 @@ class BaseTensorMesh(BaseMesh):
'Fz' -> z-component of field defined on faces
'N' -> scalar field defined on nodes
'CC' -> scalar field defined on cell centers
'CCVx' -> x-component of vector field defined on cell centers
'CCVy' -> y-component of vector field defined on cell centers
'CCVz' -> z-component of vector field defined on cell centers
"""
if self._meshType == 'CYL' and self.isSymmetric and locType in ['Ex','Ez','Fy']:
raise Exception('Symmetric CylMesh does not support %s interpolation, as this variable does not exist.' % locType)
@@ -257,6 +260,16 @@ class BaseTensorMesh(BaseMesh):
Q = sp.hstack(components)
elif locType in ['CC', 'N']:
Q = Utils.interpmat(loc, *self.getTensor(locType))
elif locType in ['CCVx', 'CCVy', 'CCVz']:
Q = Utils.interpmat(loc, *self.getTensor('CC'))
Z = Utils.spzeros(loc.shape[0],self.nC)
if locType == 'CCVx':
Q = sp.hstack([Q,Z,Z])
elif locType == 'CCVy':
Q = sp.hstack([Z,Q,Z])
elif locType == 'CCVz':
Q = sp.hstack([Z,Z,Q])
else:
raise NotImplementedError('getInterpolationMat: locType=='+locType+' and mesh.dim=='+str(self.dim))
SimPEG/Survey.py
+5 -2
View File
@@ -34,7 +34,7 @@ class BaseRx(object):
"""Number of data in the receiver."""
return self.locs.shape[0]
def getP(self, mesh):
def getP(self, mesh, projGLoc=None):
"""
Returns the projection matrices as a
list for all components collected by
@@ -47,7 +47,10 @@ class BaseRx(object):
if mesh in self._Ps:
return self._Ps[mesh]
P = mesh.getInterpolationMat(self.locs, self.projGLoc)
if projGLoc is None:
projGLoc = self.projGLoc
P = mesh.getInterpolationMat(self.locs, projGLoc)
if self.storeProjections:
self._Ps[mesh] = P
return P
SimPEG/__init__.py
+1 -1
View File
@@ -15,7 +15,7 @@ import Directives
import Inversion
import Tests
__version__ = '0.1.9'
__version__ = '0.1.10'
__author__ = 'Rowan Cockett'
__license__ = 'MIT'
__copyright__ = 'Copyright 2014 Rowan Cockett'