mirror of
https://github.com/wassname/simpeg.git
synced 2026-06-27 19:48:52 +08:00
Merge branch 'master' into docs
# Conflicts: # SimPEG/Mesh/TensorMesh.py
This commit is contained in:
+1
-1
@@ -1,4 +1,4 @@
|
||||
[bumpversion]
|
||||
current_version = 0.1.9
|
||||
current_version = 0.1.10
|
||||
files = setup.py SimPEG/__init__.py docs/conf.py
|
||||
|
||||
|
||||
+48
-64
@@ -55,8 +55,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
|
||||
|
||||
@@ -79,30 +78,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)
|
||||
|
||||
@@ -133,32 +121,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')
|
||||
|
||||
@@ -175,10 +159,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)
|
||||
|
||||
@@ -214,9 +198,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)
|
||||
@@ -240,7 +224,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
|
||||
|
||||
@@ -281,7 +265,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
|
||||
|
||||
@@ -325,9 +309,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)
|
||||
@@ -355,7 +339,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
|
||||
@@ -412,7 +396,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
|
||||
|
||||
@@ -474,9 +458,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)
|
||||
@@ -505,7 +489,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
|
||||
|
||||
@@ -526,16 +510,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):
|
||||
@@ -562,7 +546,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
|
||||
|
||||
@@ -612,9 +596,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)
|
||||
@@ -639,7 +623,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
|
||||
|
||||
@@ -656,11 +640,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):
|
||||
"""
|
||||
@@ -682,7 +666,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
|
||||
|
||||
|
||||
+698
-313
File diff suppressed because it is too large
Load Diff
+117
-79
@@ -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))
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -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.')
|
||||
|
||||
|
||||
@@ -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
|
||||
|
||||
|
||||
|
||||
@@ -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
|
||||
@@ -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
|
||||
|
||||
+572
-559
File diff suppressed because it is too large
Load Diff
+5
-2
@@ -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
|
||||
|
||||
+1
-1
@@ -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'
|
||||
|
||||
+2
-2
@@ -51,9 +51,9 @@ copyright = u'2013, SimPEG Developers'
|
||||
# built documents.
|
||||
#
|
||||
# The short X.Y version.
|
||||
version = '0.1.9'
|
||||
version = '0.1.10'
|
||||
# The full version, including alpha/beta/rc tags.
|
||||
release = '0.1.9'
|
||||
release = '0.1.10'
|
||||
|
||||
# The language for content autogenerated by Sphinx. Refer to documentation
|
||||
# for a list of supported languages.
|
||||
|
||||
@@ -77,7 +77,7 @@ with open("README.rst") as f:
|
||||
|
||||
setup(
|
||||
name = "SimPEG",
|
||||
version = "0.1.9",
|
||||
version = "0.1.10",
|
||||
packages = find_packages(),
|
||||
install_requires = ['numpy>=1.7',
|
||||
'scipy>=0.13',
|
||||
|
||||
@@ -3,125 +3,75 @@ from SimPEG import *
|
||||
from SimPEG import EM
|
||||
import sys
|
||||
from scipy.constants import mu_0
|
||||
from SimPEG.EM.Utils.testingUtils import getFDEMProblem
|
||||
from SimPEG.EM.Utils.testingUtils import getFDEMProblem, crossCheckTest
|
||||
|
||||
testEB = True
|
||||
testHJ = True
|
||||
|
||||
testEJ = True
|
||||
testBH = True
|
||||
verbose = False
|
||||
|
||||
TOL = 1e-5
|
||||
FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order
|
||||
CONDUCTIVITY = 1e1
|
||||
MU = mu_0
|
||||
freq = 1e-1
|
||||
addrandoms = True
|
||||
TOLEBHJ = 1e-5
|
||||
TOLEJHB = 1 # averaging and more sensitive to boundary condition violations (ie. the impact of violating the boundary conditions in each case is different.)
|
||||
#TODO: choose better testing parameters to lower this
|
||||
|
||||
SrcList = ['RawVec', 'MagDipole_Bfield', 'MagDipole', 'CircularLoop']
|
||||
|
||||
|
||||
def crossCheckTest(fdemType, comp):
|
||||
|
||||
l2norm = lambda r: np.sqrt(r.dot(r))
|
||||
|
||||
prb1 = getFDEMProblem(fdemType, comp, SrcList, freq, verbose)
|
||||
mesh = prb1.mesh
|
||||
print 'Cross Checking Forward: %s formulation - %s' % (fdemType, comp)
|
||||
m = np.log(np.ones(mesh.nC)*CONDUCTIVITY)
|
||||
mu = np.log(np.ones(mesh.nC)*MU)
|
||||
|
||||
if addrandoms is True:
|
||||
m = m + np.random.randn(mesh.nC)*np.log(CONDUCTIVITY)*1e-1
|
||||
mu = mu + np.random.randn(mesh.nC)*MU*1e-1
|
||||
|
||||
# prb1.PropMap.PropModel.mu = mu
|
||||
# prb1.PropMap.PropModel.mui = 1./mu
|
||||
survey1 = prb1.survey
|
||||
d1 = survey1.dpred(m)
|
||||
|
||||
if verbose:
|
||||
print ' Problem 1 solved'
|
||||
|
||||
if fdemType == 'e':
|
||||
prb2 = getFDEMProblem('b', comp, SrcList, freq, verbose)
|
||||
elif fdemType == 'b':
|
||||
prb2 = getFDEMProblem('e', comp, SrcList, freq, verbose)
|
||||
elif fdemType == 'j':
|
||||
prb2 = getFDEMProblem('h', comp, SrcList, freq, verbose)
|
||||
elif fdemType == 'h':
|
||||
prb2 = getFDEMProblem('j', comp, SrcList, freq, verbose)
|
||||
else:
|
||||
raise NotImplementedError()
|
||||
|
||||
# prb2.mu = mu
|
||||
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(l2norm(d1)))),FLR])
|
||||
print l2norm(d1), l2norm(d2), l2r , tol, l2r < tol
|
||||
return l2r < tol
|
||||
|
||||
|
||||
class FDEM_CrossCheck(unittest.TestCase):
|
||||
if testEB:
|
||||
def test_EB_CrossCheck_exr_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'exr'))
|
||||
self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'exr', verbose=verbose))
|
||||
def test_EB_CrossCheck_eyr_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'eyr'))
|
||||
self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'eyr', verbose=verbose))
|
||||
def test_EB_CrossCheck_ezr_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'ezr'))
|
||||
self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'ezr', verbose=verbose))
|
||||
def test_EB_CrossCheck_exi_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'exi'))
|
||||
self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'exi', verbose=verbose))
|
||||
def test_EB_CrossCheck_eyi_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'eyi'))
|
||||
self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'eyi', verbose=verbose))
|
||||
def test_EB_CrossCheck_ezi_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'ezi'))
|
||||
self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'ezi', verbose=verbose))
|
||||
|
||||
def test_EB_CrossCheck_bxr_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'bxr'))
|
||||
self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'bxr', verbose=verbose))
|
||||
def test_EB_CrossCheck_byr_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'byr'))
|
||||
self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'byr', verbose=verbose))
|
||||
def test_EB_CrossCheck_bzr_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'bzr'))
|
||||
self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'bzr', verbose=verbose))
|
||||
def test_EB_CrossCheck_bxi_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'bxi'))
|
||||
self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'bxi', verbose=verbose))
|
||||
def test_EB_CrossCheck_byi_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'byi'))
|
||||
self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'byi', verbose=verbose))
|
||||
def test_EB_CrossCheck_bzi_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'bzi'))
|
||||
self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'bzi', verbose=verbose))
|
||||
|
||||
if testHJ:
|
||||
def test_HJ_CrossCheck_jxr_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'jxr'))
|
||||
self.assertTrue(crossCheckTest(SrcList, 'j', 'h', 'jxr', verbose=verbose))
|
||||
def test_HJ_CrossCheck_jyr_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'jyr'))
|
||||
self.assertTrue(crossCheckTest(SrcList, 'j', 'h', 'jyr', verbose=verbose))
|
||||
def test_HJ_CrossCheck_jzr_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'jzr'))
|
||||
self.assertTrue(crossCheckTest(SrcList, 'j', 'h', 'jzr', verbose=verbose))
|
||||
def test_HJ_CrossCheck_jxi_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'jxi'))
|
||||
self.assertTrue(crossCheckTest(SrcList, 'j', 'h', 'jxi', verbose=verbose))
|
||||
def test_HJ_CrossCheck_jyi_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'jyi'))
|
||||
self.assertTrue(crossCheckTest(SrcList, 'j', 'h', 'jyi', verbose=verbose))
|
||||
def test_HJ_CrossCheck_jzi_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'jzi'))
|
||||
self.assertTrue(crossCheckTest(SrcList, 'j', 'h', 'jzi', verbose=verbose))
|
||||
|
||||
def test_HJ_CrossCheck_hxr_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'hxr'))
|
||||
self.assertTrue(crossCheckTest(SrcList, 'j', 'h', 'hxr', verbose=verbose))
|
||||
def test_HJ_CrossCheck_hyr_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'hyr'))
|
||||
self.assertTrue(crossCheckTest(SrcList, 'j', 'h', 'hyr', verbose=verbose))
|
||||
def test_HJ_CrossCheck_hzr_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'hzr'))
|
||||
self.assertTrue(crossCheckTest(SrcList, 'j', 'h', 'hzr', verbose=verbose))
|
||||
def test_HJ_CrossCheck_hxi_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'hxi'))
|
||||
self.assertTrue(crossCheckTest(SrcList, 'j', 'h', 'hxi', verbose=verbose))
|
||||
def test_HJ_CrossCheck_hyi_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'hyi'))
|
||||
self.assertTrue(crossCheckTest(SrcList, 'j', 'h', 'hyi', verbose=verbose))
|
||||
def test_HJ_CrossCheck_hzi_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'hzi'))
|
||||
self.assertTrue(crossCheckTest(SrcList, 'j', 'h', 'hzi', verbose=verbose))
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
@@ -0,0 +1,125 @@
|
||||
import unittest
|
||||
from SimPEG import *
|
||||
from SimPEG import EM
|
||||
import sys
|
||||
from scipy.constants import mu_0
|
||||
from SimPEG.EM.Utils.testingUtils import getFDEMProblem, crossCheckTest
|
||||
|
||||
testEJ = True
|
||||
testBH = True
|
||||
|
||||
TOLEJHB = 1 # averaging and more sensitive to boundary condition violations (ie. the impact of violating the boundary conditions in each case is different.)
|
||||
#TODO: choose better testing parameters to lower this
|
||||
|
||||
SrcList = ['RawVec', 'MagDipole', 'MagDipole_Bfield', 'MagDipole', 'CircularLoop']
|
||||
|
||||
|
||||
class FDEM_CrossCheck(unittest.TestCase):
|
||||
if testEJ:
|
||||
def test_EJ_CrossCheck_jxr_Jform(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'jxr', TOL=TOLEJHB))
|
||||
def test_EJ_CrossCheck_jyr_Jform(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'jyr', TOL=TOLEJHB))
|
||||
def test_EJ_CrossCheck_jzr_Jform(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'jzr', TOL=TOLEJHB))
|
||||
def test_EJ_CrossCheck_jxi_Jform(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'jxi', TOL=TOLEJHB))
|
||||
def test_EJ_CrossCheck_jyi_Jform(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'jyi', TOL=TOLEJHB))
|
||||
def test_EJ_CrossCheck_jzi_Jform(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'jzi', TOL=TOLEJHB))
|
||||
|
||||
def test_EJ_CrossCheck_exr_Jform(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'exr', TOL=TOLEJHB))
|
||||
def test_EJ_CrossCheck_eyr_Jform(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'eyr', TOL=TOLEJHB))
|
||||
def test_EJ_CrossCheck_ezr_Jform(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'ezr', TOL=TOLEJHB))
|
||||
def test_EJ_CrossCheck_exi_Jform(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'exi', TOL=TOLEJHB))
|
||||
def test_EJ_CrossCheck_eyi_Jform(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'eyi', TOL=TOLEJHB))
|
||||
def test_EJ_CrossCheck_ezi_Jform(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'ezi', TOL=TOLEJHB))
|
||||
|
||||
def test_EJ_CrossCheck_bxr_Jform(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'bxr', TOL=TOLEJHB))
|
||||
def test_EJ_CrossCheck_byr_Jform(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'byr', TOL=TOLEJHB))
|
||||
def test_EJ_CrossCheck_bzr_Jform(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'bzr', TOL=TOLEJHB))
|
||||
def test_EJ_CrossCheck_bxi_Jform(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'bxi', TOL=TOLEJHB))
|
||||
def test_EJ_CrossCheck_byi_Jform(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'byi', TOL=TOLEJHB))
|
||||
def test_EJ_CrossCheck_bzi_Jform(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'bzi', TOL=TOLEJHB))
|
||||
|
||||
def test_EJ_CrossCheck_hxr_Jform(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'hxr', TOL=TOLEJHB))
|
||||
def test_EJ_CrossCheck_hyr_Jform(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'hyr', TOL=TOLEJHB))
|
||||
def test_EJ_CrossCheck_hzr_Jform(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'hzr', TOL=TOLEJHB))
|
||||
def test_EJ_CrossCheck_hxi_Jform(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'hxi', TOL=TOLEJHB))
|
||||
def test_EJ_CrossCheck_hyi_Jform(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'hyi', TOL=TOLEJHB))
|
||||
def test_EJ_CrossCheck_hzi_Jform(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'hzi', TOL=TOLEJHB))
|
||||
|
||||
if testBH:
|
||||
def test_HB_CrossCheck_jxr_Jform(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'jxr', TOL=TOLEJHB))
|
||||
def test_HB_CrossCheck_jyr_Jform(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'jyr', TOL=TOLEJHB))
|
||||
def test_HB_CrossCheck_jzr_Jform(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'jzr', TOL=TOLEJHB))
|
||||
def test_HB_CrossCheck_jxi_Jform(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'jxi', TOL=TOLEJHB))
|
||||
def test_HB_CrossCheck_jyi_Jform(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'jyi', TOL=TOLEJHB))
|
||||
def test_HB_CrossCheck_jzi_Jform(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'jzi', TOL=TOLEJHB))
|
||||
|
||||
def test_HB_CrossCheck_exr_Jform(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'exr', TOL=TOLEJHB))
|
||||
def test_HB_CrossCheck_eyr_Jform(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'eyr', TOL=TOLEJHB))
|
||||
def test_HB_CrossCheck_ezr_Jform(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'ezr', TOL=TOLEJHB))
|
||||
def test_HB_CrossCheck_exi_Jform(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'exi', TOL=TOLEJHB))
|
||||
def test_HB_CrossCheck_eyi_Jform(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'eyi', TOL=TOLEJHB))
|
||||
def test_HB_CrossCheck_ezi_Jform(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'ezi', TOL=TOLEJHB))
|
||||
|
||||
def test_HB_CrossCheck_bxr_Jform(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'bxr', TOL=TOLEJHB))
|
||||
def test_HB_CrossCheck_byr_Jform(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'byr', TOL=TOLEJHB))
|
||||
def test_HB_CrossCheck_bzr_Jform(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'bzr', TOL=TOLEJHB))
|
||||
def test_HB_CrossCheck_bxi_Jform(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'bxi', TOL=TOLEJHB))
|
||||
def test_HB_CrossCheck_byi_Jform(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'byi', TOL=TOLEJHB))
|
||||
def test_HB_CrossCheck_bzi_Jform(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'bzi', TOL=TOLEJHB))
|
||||
|
||||
def test_HB_CrossCheck_hxr_Jform(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'hxr', TOL=TOLEJHB))
|
||||
def test_HB_CrossCheck_hyr_Jform(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'hyr', TOL=TOLEJHB))
|
||||
def test_HB_CrossCheck_hzr_Jform(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'hzr', TOL=TOLEJHB))
|
||||
def test_HB_CrossCheck_hxi_Jform(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'hxi', TOL=TOLEJHB))
|
||||
def test_HB_CrossCheck_hyi_Jform(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'hyi', TOL=TOLEJHB))
|
||||
def test_HB_CrossCheck_hzi_Jform(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'hzi', TOL=TOLEJHB))
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
@@ -0,0 +1,128 @@
|
||||
import unittest
|
||||
from SimPEG import *
|
||||
from SimPEG import EM
|
||||
import sys
|
||||
from scipy.constants import mu_0
|
||||
from SimPEG.EM.Utils.testingUtils import getFDEMProblem, crossCheckTest
|
||||
|
||||
testEB = True
|
||||
testHJ = True
|
||||
testEJ = True
|
||||
testBH = True
|
||||
verbose = False
|
||||
|
||||
TOLEJHB = 1 # averaging and more sensitive to boundary condition violations (ie. the impact of violating the boundary conditions in each case is different.)
|
||||
#TODO: choose better testing parameters to lower this
|
||||
|
||||
SrcList = ['RawVec', 'MagDipole_Bfield', 'MagDipole', 'CircularLoop']
|
||||
|
||||
|
||||
class FDEM_CrossCheck(unittest.TestCase):
|
||||
if testBH:
|
||||
def test_BH_CrossCheck_jxr(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jxr', verbose=verbose, TOL=TOLEJHB))
|
||||
def test_BH_CrossCheck_jyr(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jyr', verbose=verbose, TOL=TOLEJHB))
|
||||
def test_BH_CrossCheck_jzr(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jzr', verbose=verbose, TOL=TOLEJHB))
|
||||
def test_BH_CrossCheck_jxi(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jxi', verbose=verbose, TOL=TOLEJHB))
|
||||
def test_BH_CrossCheck_jyi(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jyi', verbose=verbose, TOL=TOLEJHB))
|
||||
def test_BH_CrossCheck_jzi(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jzi', verbose=verbose, TOL=TOLEJHB))
|
||||
|
||||
def test_BH_CrossCheck_exr(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'exr', verbose=verbose, TOL=TOLEJHB))
|
||||
def test_BH_CrossCheck_eyr(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'eyr', verbose=verbose, TOL=TOLEJHB))
|
||||
def test_BH_CrossCheck_ezr(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'ezr', verbose=verbose, TOL=TOLEJHB))
|
||||
def test_BH_CrossCheck_exi(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'exi', verbose=verbose, TOL=TOLEJHB))
|
||||
def test_BH_CrossCheck_eyi(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'eyi', verbose=verbose, TOL=TOLEJHB))
|
||||
def test_BH_CrossCheck_ezi(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'ezi', verbose=verbose, TOL=TOLEJHB))
|
||||
|
||||
def test_BH_CrossCheck_bxr(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'bxr', verbose=verbose, TOL=TOLEJHB))
|
||||
def test_BH_CrossCheck_byr(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'byr', verbose=verbose, TOL=TOLEJHB))
|
||||
def test_BH_CrossCheck_bzr(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'bzr', verbose=verbose, TOL=TOLEJHB))
|
||||
def test_BH_CrossCheck_bxi(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'bxi', verbose=verbose, TOL=TOLEJHB))
|
||||
def test_BH_CrossCheck_byi(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'byi', verbose=verbose, TOL=TOLEJHB))
|
||||
def test_BH_CrossCheck_bzi(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'bzi', verbose=verbose, TOL=TOLEJHB))
|
||||
|
||||
def test_BH_CrossCheck_hxr(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hxr', verbose=verbose, TOL=TOLEJHB))
|
||||
def test_BH_CrossCheck_hyr(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hyr', verbose=verbose, TOL=TOLEJHB))
|
||||
def test_BH_CrossCheck_hzr(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hzr', verbose=verbose, TOL=TOLEJHB))
|
||||
def test_BH_CrossCheck_hxi(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hxi', verbose=verbose, TOL=TOLEJHB))
|
||||
def test_BH_CrossCheck_hyi(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hyi', verbose=verbose, TOL=TOLEJHB))
|
||||
def test_BH_CrossCheck_hzi(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hzi', verbose=verbose, TOL=TOLEJHB))
|
||||
|
||||
if testBH:
|
||||
def test_BH_CrossCheck_jxr(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jxr', verbose=verbose, TOL=TOLEJHB))
|
||||
def test_BH_CrossCheck_jyr(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jyr', verbose=verbose, TOL=TOLEJHB))
|
||||
def test_BH_CrossCheck_jzr(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jzr', verbose=verbose, TOL=TOLEJHB))
|
||||
def test_BH_CrossCheck_jxi(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jxi', verbose=verbose, TOL=TOLEJHB))
|
||||
def test_BH_CrossCheck_jyi(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jyi', verbose=verbose, TOL=TOLEJHB))
|
||||
def test_BH_CrossCheck_jzi(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jzi', verbose=verbose, TOL=TOLEJHB))
|
||||
|
||||
def test_BH_CrossCheck_exr(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'exr', verbose=verbose, TOL=TOLEJHB))
|
||||
def test_BH_CrossCheck_eyr(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'eyr', verbose=verbose, TOL=TOLEJHB))
|
||||
def test_BH_CrossCheck_ezr(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'ezr', verbose=verbose, TOL=TOLEJHB))
|
||||
def test_BH_CrossCheck_exi(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'exi', verbose=verbose, TOL=TOLEJHB))
|
||||
def test_BH_CrossCheck_eyi(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'eyi', verbose=verbose, TOL=TOLEJHB))
|
||||
def test_BH_CrossCheck_ezi(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'ezi', verbose=verbose, TOL=TOLEJHB))
|
||||
|
||||
def test_BH_CrossCheck_bxr(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'bxr', verbose=verbose, TOL=TOLEJHB))
|
||||
def test_BH_CrossCheck_byr(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'byr', verbose=verbose, TOL=TOLEJHB))
|
||||
def test_BH_CrossCheck_bzr(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'bzr', verbose=verbose, TOL=TOLEJHB))
|
||||
def test_BH_CrossCheck_bxi(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'bxi', verbose=verbose, TOL=TOLEJHB))
|
||||
def test_BH_CrossCheck_byi(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'byi', verbose=verbose, TOL=TOLEJHB))
|
||||
def test_BH_CrossCheck_bzi(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'bzi', verbose=verbose, TOL=TOLEJHB))
|
||||
|
||||
def test_BH_CrossCheck_hxr(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hxr', verbose=verbose, TOL=TOLEJHB))
|
||||
def test_BH_CrossCheck_hyr(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hyr', verbose=verbose, TOL=TOLEJHB))
|
||||
def test_BH_CrossCheck_hzr(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hzr', verbose=verbose, TOL=TOLEJHB))
|
||||
def test_BH_CrossCheck_hxi(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hxi', verbose=verbose, TOL=TOLEJHB))
|
||||
def test_BH_CrossCheck_hyi(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hyi', verbose=verbose, TOL=TOLEJHB))
|
||||
def test_BH_CrossCheck_hzi(self):
|
||||
self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hzi', verbose=verbose, TOL=TOLEJHB))
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
+57
-57
@@ -5,8 +5,8 @@ import sys
|
||||
from scipy.constants import mu_0
|
||||
from SimPEG.EM.Utils.testingUtils import getFDEMProblem
|
||||
|
||||
testEB = True
|
||||
testHJ = True
|
||||
testE = True
|
||||
testB = True
|
||||
|
||||
verbose = False
|
||||
|
||||
@@ -17,10 +17,10 @@ MU = mu_0
|
||||
freq = 1e-1
|
||||
addrandoms = True
|
||||
|
||||
SrcType = 'RawVec' #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVec'
|
||||
SrcList = ['RawVec', 'MagDipole'] #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVec'
|
||||
|
||||
def adjointTest(fdemType, comp):
|
||||
prb = getFDEMProblem(fdemType, comp, [SrcType], freq)
|
||||
prb = getFDEMProblem(fdemType, comp, SrcList, freq)
|
||||
print 'Adjoint %s formulation - %s' % (fdemType, comp)
|
||||
|
||||
m = np.log(np.ones(prb.mapping.nP)*CONDUCTIVITY)
|
||||
@@ -45,7 +45,7 @@ def adjointTest(fdemType, comp):
|
||||
return np.abs(vJw - wJtv) < tol
|
||||
|
||||
class FDEM_AdjointTests(unittest.TestCase):
|
||||
if testEB:
|
||||
if testE:
|
||||
def test_Jtvec_adjointTest_exr_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'exr'))
|
||||
def test_Jtvec_adjointTest_eyr_Eform(self):
|
||||
@@ -72,6 +72,33 @@ class FDEM_AdjointTests(unittest.TestCase):
|
||||
def test_Jtvec_adjointTest_bzi_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'bzi'))
|
||||
|
||||
def test_Jtvec_adjointTest_jxr_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'jxr'))
|
||||
def test_Jtvec_adjointTest_jyr_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'jyr'))
|
||||
def test_Jtvec_adjointTest_jzr_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'jzr'))
|
||||
def test_Jtvec_adjointTest_jxi_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'jxi'))
|
||||
def test_Jtvec_adjointTest_jyi_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'jyi'))
|
||||
def test_Jtvec_adjointTest_jzi_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'jzi'))
|
||||
|
||||
def test_Jtvec_adjointTest_hxr_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'hxr'))
|
||||
def test_Jtvec_adjointTest_hyr_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'hyr'))
|
||||
def test_Jtvec_adjointTest_hzr_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'hzr'))
|
||||
def test_Jtvec_adjointTest_hxi_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'hxi'))
|
||||
def test_Jtvec_adjointTest_hyi_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'hyi'))
|
||||
def test_Jtvec_adjointTest_hzi_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'hzi'))
|
||||
|
||||
if testB:
|
||||
def test_Jtvec_adjointTest_exr_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'exr'))
|
||||
def test_Jtvec_adjointTest_eyr_Bform(self):
|
||||
@@ -84,6 +111,7 @@ class FDEM_AdjointTests(unittest.TestCase):
|
||||
self.assertTrue(adjointTest('b', 'eyi'))
|
||||
def test_Jtvec_adjointTest_ezi_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'ezi'))
|
||||
|
||||
def test_Jtvec_adjointTest_bxr_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'bxr'))
|
||||
def test_Jtvec_adjointTest_byr_Bform(self):
|
||||
@@ -97,59 +125,31 @@ class FDEM_AdjointTests(unittest.TestCase):
|
||||
def test_Jtvec_adjointTest_bzi_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'bzi'))
|
||||
|
||||
def test_Jtvec_adjointTest_jxr_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'jxr'))
|
||||
def test_Jtvec_adjointTest_jyr_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'jyr'))
|
||||
def test_Jtvec_adjointTest_jzr_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'jzr'))
|
||||
def test_Jtvec_adjointTest_jxi_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'jxi'))
|
||||
def test_Jtvec_adjointTest_jyi_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'jyi'))
|
||||
def test_Jtvec_adjointTest_jzi_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'jzi'))
|
||||
|
||||
if testHJ:
|
||||
def test_Jtvec_adjointTest_jxr_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'jxr'))
|
||||
def test_Jtvec_adjointTest_jyr_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'jyr'))
|
||||
def test_Jtvec_adjointTest_jzr_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'jzr'))
|
||||
def test_Jtvec_adjointTest_jxi_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'jxi'))
|
||||
def test_Jtvec_adjointTest_jyi_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'jyi'))
|
||||
def test_Jtvec_adjointTest_jzi_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'jzi'))
|
||||
|
||||
def test_Jtvec_adjointTest_hxr_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'hxr'))
|
||||
def test_Jtvec_adjointTest_hyr_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'hyr'))
|
||||
def test_Jtvec_adjointTest_hzr_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'hzr'))
|
||||
def test_Jtvec_adjointTest_hxi_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'hxi'))
|
||||
def test_Jtvec_adjointTest_hyi_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'hyi'))
|
||||
def test_Jtvec_adjointTest_hzi_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'hzi'))
|
||||
|
||||
def test_Jtvec_adjointTest_hxr_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'hxr'))
|
||||
def test_Jtvec_adjointTest_hyr_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'hyr'))
|
||||
def test_Jtvec_adjointTest_hzr_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'hzr'))
|
||||
def test_Jtvec_adjointTest_hxi_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'hxi'))
|
||||
def test_Jtvec_adjointTest_hyi_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'hyi'))
|
||||
def test_Jtvec_adjointTest_hzi_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'hzi'))
|
||||
|
||||
def test_Jtvec_adjointTest_hxr_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'jxr'))
|
||||
def test_Jtvec_adjointTest_hyr_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'jyr'))
|
||||
def test_Jtvec_adjointTest_hzr_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'jzr'))
|
||||
def test_Jtvec_adjointTest_hxi_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'jxi'))
|
||||
def test_Jtvec_adjointTest_hyi_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'jyi'))
|
||||
def test_Jtvec_adjointTest_hzi_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'jzi'))
|
||||
def test_Jtvec_adjointTest_hxr_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'hxr'))
|
||||
def test_Jtvec_adjointTest_hyr_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'hyr'))
|
||||
def test_Jtvec_adjointTest_hzr_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'hzr'))
|
||||
def test_Jtvec_adjointTest_hxi_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'hxi'))
|
||||
def test_Jtvec_adjointTest_hyi_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'hyi'))
|
||||
def test_Jtvec_adjointTest_hzi_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'hzi'))
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
@@ -0,0 +1,155 @@
|
||||
import unittest
|
||||
from SimPEG import *
|
||||
from SimPEG import EM
|
||||
import sys
|
||||
from scipy.constants import mu_0
|
||||
from SimPEG.EM.Utils.testingUtils import getFDEMProblem
|
||||
|
||||
testJ = True
|
||||
testH = True
|
||||
|
||||
verbose = False
|
||||
|
||||
TOL = 1e-5
|
||||
FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order
|
||||
CONDUCTIVITY = 1e1
|
||||
MU = mu_0
|
||||
freq = 1e-1
|
||||
addrandoms = True
|
||||
|
||||
SrcList = ['RawVec', 'MagDipole'] #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVec'
|
||||
|
||||
def adjointTest(fdemType, comp):
|
||||
prb = getFDEMProblem(fdemType, comp, SrcList, freq)
|
||||
print 'Adjoint %s formulation - %s' % (fdemType, comp)
|
||||
|
||||
m = np.log(np.ones(prb.mapping.nP)*CONDUCTIVITY)
|
||||
mu = np.ones(prb.mesh.nC)*MU
|
||||
|
||||
if addrandoms is True:
|
||||
m = m + np.random.randn(prb.mapping.nP)*np.log(CONDUCTIVITY)*1e-1
|
||||
mu = mu + np.random.randn(prb.mesh.nC)*MU*1e-1
|
||||
|
||||
survey = prb.survey
|
||||
u = prb.fields(m)
|
||||
|
||||
v = np.random.rand(survey.nD)
|
||||
w = np.random.rand(prb.mesh.nC)
|
||||
|
||||
vJw = v.dot(prb.Jvec(m, w, u))
|
||||
wJtv = w.dot(prb.Jtvec(m, v, u))
|
||||
tol = np.max([TOL*(10**int(np.log10(np.abs(vJw)))),FLR])
|
||||
print vJw, wJtv, vJw - wJtv, tol, np.abs(vJw - wJtv) < tol
|
||||
return np.abs(vJw - wJtv) < tol
|
||||
|
||||
class FDEM_AdjointTests(unittest.TestCase):
|
||||
|
||||
if testJ:
|
||||
def test_Jtvec_adjointTest_jxr_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'jxr'))
|
||||
def test_Jtvec_adjointTest_jyr_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'jyr'))
|
||||
def test_Jtvec_adjointTest_jzr_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'jzr'))
|
||||
def test_Jtvec_adjointTest_jxi_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'jxi'))
|
||||
def test_Jtvec_adjointTest_jyi_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'jyi'))
|
||||
def test_Jtvec_adjointTest_jzi_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'jzi'))
|
||||
|
||||
def test_Jtvec_adjointTest_hxr_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'hxr'))
|
||||
def test_Jtvec_adjointTest_hyr_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'hyr'))
|
||||
def test_Jtvec_adjointTest_hzr_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'hzr'))
|
||||
def test_Jtvec_adjointTest_hxi_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'hxi'))
|
||||
def test_Jtvec_adjointTest_hyi_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'hyi'))
|
||||
def test_Jtvec_adjointTest_hzi_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'hzi'))
|
||||
|
||||
def test_Jtvec_adjointTest_exr_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'exr'))
|
||||
def test_Jtvec_adjointTest_eyr_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'eyr'))
|
||||
def test_Jtvec_adjointTest_ezr_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'ezr'))
|
||||
def test_Jtvec_adjointTest_exi_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'exi'))
|
||||
def test_Jtvec_adjointTest_eyi_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'eyi'))
|
||||
def test_Jtvec_adjointTest_ezi_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'ezi'))
|
||||
|
||||
def test_Jtvec_adjointTest_bxr_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'bxr'))
|
||||
def test_Jtvec_adjointTest_byr_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'byr'))
|
||||
def test_Jtvec_adjointTest_bzr_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'bzr'))
|
||||
def test_Jtvec_adjointTest_bxi_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'bxi'))
|
||||
def test_Jtvec_adjointTest_byi_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'byi'))
|
||||
def test_Jtvec_adjointTest_bzi_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'bzi'))
|
||||
|
||||
if testH:
|
||||
def test_Jtvec_adjointTest_hxr_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'hxr'))
|
||||
def test_Jtvec_adjointTest_hyr_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'hyr'))
|
||||
def test_Jtvec_adjointTest_hzr_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'hzr'))
|
||||
def test_Jtvec_adjointTest_hxi_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'hxi'))
|
||||
def test_Jtvec_adjointTest_hyi_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'hyi'))
|
||||
def test_Jtvec_adjointTest_hzi_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'hzi'))
|
||||
|
||||
def test_Jtvec_adjointTest_jxr_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'jxr'))
|
||||
def test_Jtvec_adjointTest_jyr_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'jyr'))
|
||||
def test_Jtvec_adjointTest_jzr_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'jzr'))
|
||||
def test_Jtvec_adjointTest_jxi_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'jxi'))
|
||||
def test_Jtvec_adjointTest_jyi_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'jyi'))
|
||||
def test_Jtvec_adjointTest_jzi_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'jzi'))
|
||||
|
||||
def test_Jtvec_adjointTest_exr_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'exr'))
|
||||
def test_Jtvec_adjointTest_eyr_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'eyr'))
|
||||
def test_Jtvec_adjointTest_ezr_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'ezr'))
|
||||
def test_Jtvec_adjointTest_exi_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'exi'))
|
||||
def test_Jtvec_adjointTest_eyi_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'eyi'))
|
||||
def test_Jtvec_adjointTest_ezi_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'ezi'))
|
||||
|
||||
def test_Jtvec_adjointTest_bxr_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'bxr'))
|
||||
def test_Jtvec_adjointTest_byr_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'byr'))
|
||||
def test_Jtvec_adjointTest_bzr_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'bzr'))
|
||||
def test_Jtvec_adjointTest_bxi_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'bxi'))
|
||||
def test_Jtvec_adjointTest_byi_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'byi'))
|
||||
def test_Jtvec_adjointTest_bzi_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'bzi'))
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
@@ -5,9 +5,11 @@ import sys
|
||||
from scipy.constants import mu_0
|
||||
from SimPEG.EM.Utils.testingUtils import getFDEMProblem
|
||||
|
||||
testDerivs = True
|
||||
testEB = True
|
||||
testHJ = True
|
||||
|
||||
testE = True
|
||||
testB = True
|
||||
testH = True
|
||||
testJ = True
|
||||
|
||||
verbose = False
|
||||
|
||||
@@ -18,12 +20,12 @@ MU = mu_0
|
||||
freq = 1e-1
|
||||
addrandoms = True
|
||||
|
||||
SrcType = 'RawVec' #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVec'
|
||||
SrcType = ['MagDipole', 'RawVec'] #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVec'
|
||||
|
||||
|
||||
def derivTest(fdemType, comp):
|
||||
|
||||
prb = getFDEMProblem(fdemType, comp, [SrcType], freq)
|
||||
prb = getFDEMProblem(fdemType, comp, SrcType, freq)
|
||||
print '%s formulation - %s' % (fdemType, comp)
|
||||
x0 = np.log(np.ones(prb.mapping.nP)*CONDUCTIVITY)
|
||||
mu = np.log(np.ones(prb.mesh.nC)*MU)
|
||||
@@ -32,9 +34,6 @@ def derivTest(fdemType, comp):
|
||||
x0 = x0 + np.random.randn(prb.mapping.nP)*np.log(CONDUCTIVITY)*1e-1
|
||||
mu = mu + np.random.randn(prb.mapping.nP)*MU*1e-1
|
||||
|
||||
# prb.PropMap.PropModel.mu = mu
|
||||
# prb.PropMap.PropModel.mui = 1./mu
|
||||
|
||||
survey = prb.survey
|
||||
def fun(x):
|
||||
return survey.dpred(x), lambda x: prb.Jvec(x0, x)
|
||||
@@ -43,7 +42,7 @@ def derivTest(fdemType, comp):
|
||||
|
||||
class FDEM_DerivTests(unittest.TestCase):
|
||||
|
||||
if testEB:
|
||||
if testE:
|
||||
def test_Jvec_exr_Eform(self):
|
||||
self.assertTrue(derivTest('e', 'exr'))
|
||||
def test_Jvec_eyr_Eform(self):
|
||||
@@ -70,6 +69,33 @@ class FDEM_DerivTests(unittest.TestCase):
|
||||
def test_Jvec_bzi_Eform(self):
|
||||
self.assertTrue(derivTest('e', 'bzi'))
|
||||
|
||||
def test_Jvec_exr_Eform(self):
|
||||
self.assertTrue(derivTest('e', 'jxr'))
|
||||
def test_Jvec_eyr_Eform(self):
|
||||
self.assertTrue(derivTest('e', 'jyr'))
|
||||
def test_Jvec_ezr_Eform(self):
|
||||
self.assertTrue(derivTest('e', 'jzr'))
|
||||
def test_Jvec_exi_Eform(self):
|
||||
self.assertTrue(derivTest('e', 'jxi'))
|
||||
def test_Jvec_eyi_Eform(self):
|
||||
self.assertTrue(derivTest('e', 'jyi'))
|
||||
def test_Jvec_ezi_Eform(self):
|
||||
self.assertTrue(derivTest('e', 'jzi'))
|
||||
|
||||
def test_Jvec_bxr_Eform(self):
|
||||
self.assertTrue(derivTest('e', 'hxr'))
|
||||
def test_Jvec_byr_Eform(self):
|
||||
self.assertTrue(derivTest('e', 'hyr'))
|
||||
def test_Jvec_bzr_Eform(self):
|
||||
self.assertTrue(derivTest('e', 'hzr'))
|
||||
def test_Jvec_bxi_Eform(self):
|
||||
self.assertTrue(derivTest('e', 'hxi'))
|
||||
def test_Jvec_byi_Eform(self):
|
||||
self.assertTrue(derivTest('e', 'hyi'))
|
||||
def test_Jvec_bzi_Eform(self):
|
||||
self.assertTrue(derivTest('e', 'hzi'))
|
||||
|
||||
if testB:
|
||||
def test_Jvec_exr_Bform(self):
|
||||
self.assertTrue(derivTest('b', 'exr'))
|
||||
def test_Jvec_eyr_Bform(self):
|
||||
@@ -96,7 +122,33 @@ class FDEM_DerivTests(unittest.TestCase):
|
||||
def test_Jvec_bzi_Bform(self):
|
||||
self.assertTrue(derivTest('b', 'bzi'))
|
||||
|
||||
if testHJ:
|
||||
def test_Jvec_jxr_Bform(self):
|
||||
self.assertTrue(derivTest('b', 'jxr'))
|
||||
def test_Jvec_jyr_Bform(self):
|
||||
self.assertTrue(derivTest('b', 'jyr'))
|
||||
def test_Jvec_jzr_Bform(self):
|
||||
self.assertTrue(derivTest('b', 'jzr'))
|
||||
def test_Jvec_jxi_Bform(self):
|
||||
self.assertTrue(derivTest('b', 'jxi'))
|
||||
def test_Jvec_jyi_Bform(self):
|
||||
self.assertTrue(derivTest('b', 'jyi'))
|
||||
def test_Jvec_jzi_Bform(self):
|
||||
self.assertTrue(derivTest('b', 'jzi'))
|
||||
|
||||
def test_Jvec_hxr_Bform(self):
|
||||
self.assertTrue(derivTest('b', 'hxr'))
|
||||
def test_Jvec_hyr_Bform(self):
|
||||
self.assertTrue(derivTest('b', 'hyr'))
|
||||
def test_Jvec_hzr_Bform(self):
|
||||
self.assertTrue(derivTest('b', 'hzr'))
|
||||
def test_Jvec_hxi_Bform(self):
|
||||
self.assertTrue(derivTest('b', 'hxi'))
|
||||
def test_Jvec_hyi_Bform(self):
|
||||
self.assertTrue(derivTest('b', 'hyi'))
|
||||
def test_Jvec_hzi_Bform(self):
|
||||
self.assertTrue(derivTest('b', 'hzi'))
|
||||
|
||||
if testJ:
|
||||
def test_Jvec_jxr_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'jxr'))
|
||||
def test_Jvec_jyr_Jform(self):
|
||||
@@ -123,6 +175,34 @@ class FDEM_DerivTests(unittest.TestCase):
|
||||
def test_Jvec_hzi_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'hzi'))
|
||||
|
||||
def test_Jvec_exr_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'exr'))
|
||||
def test_Jvec_eyr_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'eyr'))
|
||||
def test_Jvec_ezr_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'ezr'))
|
||||
def test_Jvec_exi_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'exi'))
|
||||
def test_Jvec_eyi_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'eyi'))
|
||||
def test_Jvec_ezi_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'ezi'))
|
||||
|
||||
def test_Jvec_bxr_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'bxr'))
|
||||
def test_Jvec_byr_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'byr'))
|
||||
def test_Jvec_bzr_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'bzr'))
|
||||
def test_Jvec_bxi_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'bxi'))
|
||||
def test_Jvec_byi_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'byi'))
|
||||
def test_Jvec_bzi_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'bzi'))
|
||||
|
||||
|
||||
if testH:
|
||||
def test_Jvec_hxr_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'hxr'))
|
||||
def test_Jvec_hyr_Hform(self):
|
||||
@@ -149,6 +229,32 @@ class FDEM_DerivTests(unittest.TestCase):
|
||||
def test_Jvec_hzi_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'jzi'))
|
||||
|
||||
def test_Jvec_exr_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'exr'))
|
||||
def test_Jvec_eyr_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'eyr'))
|
||||
def test_Jvec_ezr_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'ezr'))
|
||||
def test_Jvec_exi_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'exi'))
|
||||
def test_Jvec_eyi_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'eyi'))
|
||||
def test_Jvec_ezi_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'ezi'))
|
||||
|
||||
def test_Jvec_bxr_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'bxr'))
|
||||
def test_Jvec_byr_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'byr'))
|
||||
def test_Jvec_bzr_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'bzr'))
|
||||
def test_Jvec_bxi_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'bxi'))
|
||||
def test_Jvec_byi_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'byi'))
|
||||
def test_Jvec_bzi_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'bzi'))
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
|
||||
Reference in New Issue
Block a user