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 'dev' into dcip/dev

Browse Source
This commit is contained in:
D Fournier
2016-04-06 09:03:37 -07:00
parent e646211e7d d8aec96080
commit 09f3f7b55b
36 changed files with 2305 additions and 1044 deletions
Download Patch File Download Diff File Expand all files Collapse all files
.bumpversion.cfg
+1 -1
View File
@@ -1,4 +1,4 @@
[bumpversion]
current_version = 0.1.9
current_version = 0.1.10
files = setup.py SimPEG/__init__.py docs/conf.py
SimPEG/DCIP/BaseDC.py
+11 -13
View File
@@ -200,11 +200,11 @@ class ProblemDC_CC(Problem.BaseProblem):
return F
def Jvec(self, m, v, u=None):
def Jvec(self, m, v, f=None):
"""
:param numpy.array m: model
:param numpy.array v: vector to multiply
:param numpy.array u: fields
:param Fields f: fields
:rtype: numpy.array
:return: Jv
@@ -225,11 +225,10 @@ class ProblemDC_CC(Problem.BaseProblem):
# Set current model; clear dependent property $\mathbf{A(m)}$
self.curModel = m
sigma = self.curModel.transform # $\sigma = \mathcal{M}(\m)$
if u is None:
if f is None:
# Run forward simulation if $u$ not provided
u = self.fields(self.curModel)[self.survey.srcList, 'phi_sol']
else:
u = u[self.survey.srcList, 'phi_sol']
f = self.fields(self.curModel)
u = f[self.survey.srcList, 'phi_sol']
D = self.mesh.faceDiv
G = self.mesh.cellGrad
@@ -251,19 +250,18 @@ class ProblemDC_CC(Problem.BaseProblem):
if self.Ainv is None:
self.Ainv = self.Solver(dA_du, **self.solverOpts)
P = self.survey.getP(self.mesh)
P = self.survey.getP(self.mesh)
Jv = - P * mkvc( self.Ainv * dCdm_x_v )
return Jv
def Jtvec(self, m, v, u=None):
def Jtvec(self, m, v, f=None):
self.curModel = m
sigma = self.curModel.transform # $\sigma = \mathcal{M}(\m)$
if u is None:
# Run forward simulation if $u$ not provided
u = self.fields(self.curModel)[self.survey.srcList, 'phi_sol']
else:
u = u[self.survey.srcList, 'phi_sol']
if f is None:
# Run forward simulation if $f$ not provided
f = self.fields(self.curModel)
u = f[self.survey.srcList, 'phi_sol']
shp = u.shape
P = self.survey.getP(self.mesh)
SimPEG/DCIP/BaseIP.py
+4 -4
View File
@@ -14,12 +14,12 @@ class SurveyIP(SurveyDC):
Survey.BaseSurvey.__init__(self, **kwargs)
self._Ps = {}
def dpred(self, m, u=None):
def dpred(self, m, f=None):
"""
Predicted data.
.. math::
d_\\text{pred} = Pu(m)
d_\\text{pred} = Pf(m)
"""
return self.prob.forward(m)
@@ -143,10 +143,10 @@ class ProblemIP(Problem.BaseProblem):
J_x_v = - P * mkvc( self.Ainv * dCdm_x_v )
return -J_x_v
def Jvec(self, m, v, u=None):
def Jvec(self, m, v, f=None):
return self.forward(v)
def Jtvec(self, m, v, u=None):
def Jtvec(self, m, v, f=None):
self.curModel = m
# sigma = self.curModel.transform # $\sigma = \mathcal{M}(\m)$
SimPEG/DataMisfit.py
+22 -26
View File
@@ -22,11 +22,11 @@ class BaseDataMisfit(object):
Utils.setKwargs(self,**kwargs)
@Utils.timeIt
def eval(self, m, u=None):
"""eval(m, u=None)
def eval(self, m, f=None):
"""eval(m, f=None)
:param numpy.array m: geophysical model
:param numpy.array u: fields
:param Fields f: fields
:rtype: float
:return: data misfit
@@ -34,11 +34,11 @@ class BaseDataMisfit(object):
raise NotImplementedError('This method should be overwritten.')
@Utils.timeIt
def evalDeriv(self, m, u=None):
"""evalDeriv(m, u=None)
def evalDeriv(self, m, f=None):
"""evalDeriv(m, f=None)
:param numpy.array m: geophysical model
:param numpy.array u: fields
:param Fields f: fields
:rtype: numpy.array
:return: data misfit derivative
@@ -47,12 +47,12 @@ class BaseDataMisfit(object):
@Utils.timeIt
def eval2Deriv(self, m, v, u=None):
"""eval2Deriv(m, v, u=None)
def eval2Deriv(self, m, v, f=None):
"""eval2Deriv(m, v, f=None)
:param numpy.array m: geophysical model
:param numpy.array v: vector to multiply
:param numpy.array u: fields
:param Fields f: fields
:rtype: numpy.array
:return: data misfit derivative
@@ -89,7 +89,7 @@ class l2_DataMisfit(BaseDataMisfit):
"""
if getattr(self, '_Wd', None) is None:
survey = self.survey
if getattr(survey,'std', None) is None:
@@ -108,24 +108,20 @@ class l2_DataMisfit(BaseDataMisfit):
self._Wd = value
@Utils.timeIt
def eval(self, m, u=None):
"eval(m, u=None)"
prob = self.prob
survey = self.survey
R = self.Wd * survey.residual(m, u=u)
def eval(self, m, f=None):
"eval(m, f=None)"
if f is None: f = self.prob.fields(m)
R = self.Wd * self.survey.residual(m, f)
return 0.5*np.vdot(R, R)
@Utils.timeIt
def evalDeriv(self, m, u=None):
"evalDeriv(m, u=None)"
prob = self.prob
survey = self.survey
if u is None: u = prob.fields(m)
return prob.Jtvec(m, self.Wd * (self.Wd * survey.residual(m, u=u)), u=u)
def evalDeriv(self, m, f=None):
"evalDeriv(m, f=None)"
if f is None: f = self.prob.fields(m)
return self.prob.Jtvec(m, self.Wd * (self.Wd * self.survey.residual(m, f=f)), f=f)
@Utils.timeIt
def eval2Deriv(self, m, v, u=None):
"eval2Deriv(m, v, u=None)"
prob = self.prob
if u is None: u = prob.fields(m)
return prob.Jtvec_approx(m, self.Wd * (self.Wd * prob.Jvec_approx(m, v, u=u)), u=u)
def eval2Deriv(self, m, v, f=None):
"eval2Deriv(m, v, f=None)"
if f is None: f = self.prob.fields(m)
return self.prob.Jtvec_approx(m, self.Wd * (self.Wd * self.prob.Jvec_approx(m, v, f=f)), f=f)
SimPEG/Directives.py
+2 -2
View File
@@ -123,10 +123,10 @@ class BetaEstimate_ByEig(InversionDirective):
if self.debug: print 'Calculating the beta0 parameter.'
m = self.invProb.curModel
u = self.invProb.getFields(m, store=True, deleteWarmstart=False)
f = self.invProb.getFields(m, store=True, deleteWarmstart=False)
x0 = np.random.rand(*m.shape)
t = x0.dot(self.dmisfit.eval2Deriv(m,x0,u=u))
t = x0.dot(self.dmisfit.eval2Deriv(m,x0,f=f))
b = x0.dot(self.reg.eval2Deriv(m, v=x0))
self.beta0 = self.beta0_ratio*(t/b)
SimPEG/EM/Base.py
+37 -14
View File
@@ -2,14 +2,14 @@ from SimPEG import Survey, Problem, Utils, Models, Maps, PropMaps, np, sp, Solve
from scipy.constants import mu_0
class EMPropMap(Maps.PropMap):
"""
"""
Property Map for EM Problems. The electrical conductivity (\\(\\sigma\\)) is the default inversion property, and the default value of the magnetic permeability is that of free space (\\(\\mu = 4\\pi\\times 10^{-7} \\) H/m)
"""
sigma = Maps.Property("Electrical Conductivity", defaultInvProp = True, propertyLink=('rho',Maps.ReciprocalMap))
mu = Maps.Property("Inverse Magnetic Permeability", defaultVal = mu_0, propertyLink=('mui',Maps.ReciprocalMap))
rho = Maps.Property("Electrical Resistivity", propertyLink=('sigma', Maps.ReciprocalMap))
rho = Maps.Property("Electrical Resistivity", propertyLink=('sigma', Maps.ReciprocalMap))
mui = Maps.Property("Inverse Magnetic Permeability", defaultVal = 1./mu_0, propertyLink=('mu', Maps.ReciprocalMap))
@@ -21,7 +21,7 @@ class BaseEMProblem(Problem.BaseProblem):
surveyPair = Survey.BaseSurvey
dataPair = Survey.Data
PropMap = EMPropMap
Solver = SimpegSolver
@@ -51,7 +51,7 @@ class BaseEMProblem(Problem.BaseProblem):
if self.mapping.muMap is not None or self.mapping.muiMap is not None:
toDelete += ['_MeMu', '_MeMuI','_MfMui','_MfMuiI']
return toDelete
@property
def Me(self):
"""
@@ -71,7 +71,7 @@ class BaseEMProblem(Problem.BaseProblem):
return self._Mf
# ----- Magnetic Permeability ----- #
# ----- Magnetic Permeability ----- #
@property
def MfMui(self):
"""
@@ -109,7 +109,7 @@ class BaseEMProblem(Problem.BaseProblem):
return self._MeMuI
# ----- Electrical Conductivity ----- #
# ----- Electrical Conductivity ----- #
#TODO: hardcoded to sigma as the model
@property
def MeSigma(self):
@@ -120,18 +120,18 @@ class BaseEMProblem(Problem.BaseProblem):
self._MeSigma = self.mesh.getEdgeInnerProduct(self.curModel.sigma)
return self._MeSigma
# TODO: This should take a vector
# TODO: This should take a vector
def MeSigmaDeriv(self, u):
"""
Derivative of MeSigma with respect to the model
"""
"""
return self.mesh.getEdgeInnerProductDeriv(self.curModel.sigma)(u) * self.curModel.sigmaDeriv
@property
def MeSigmaI(self):
"""
Inverse of the edge inner product matrix for \\(\\sigma\\).
Inverse of the edge inner product matrix for \\(\\sigma\\).
"""
if getattr(self, '_MeSigmaI', None) is None:
self._MeSigmaI = self.mesh.getEdgeInnerProduct(self.curModel.sigma, invMat=True)
@@ -140,8 +140,8 @@ class BaseEMProblem(Problem.BaseProblem):
# TODO: This should take a vector
def MeSigmaIDeriv(self, u):
"""
Derivative of :code:`MeSigma` with respect to the model
"""
Derivative of :code:`MeSigma` with respect to the model
"""
# TODO: only works for diagonal tensors. getEdgeInnerProductDeriv, invMat=True should be implemented in SimPEG
dMeSigmaI_dI = -self.MeSigmaI**2
@@ -163,7 +163,7 @@ class BaseEMProblem(Problem.BaseProblem):
# TODO: This should take a vector
def MfRhoDeriv(self,u):
"""
Derivative of :code:`MfRho` with respect to the model.
Derivative of :code:`MfRho` with respect to the model.
"""
return self.mesh.getFaceInnerProductDeriv(self.curModel.rho)(u) * (-Utils.sdiag(self.curModel.rho**2) * self.curModel.sigmaDeriv)
# self.curModel.rhoDeriv
@@ -181,6 +181,29 @@ class BaseEMProblem(Problem.BaseProblem):
# TODO: This should take a vector
def MfRhoIDeriv(self,u):
"""
Derivative of :code:`MfRhoI` with respect to the model.
Derivative of :code:`MfRhoI` with respect to the model.
"""
return self.mesh.getFaceInnerProductDeriv(self.curModel.rho, invMat=True)(u) * self.curModel.rhoDeriv
class BaseEMSurvey(Survey.BaseSurvey):
def __init__(self, srcList, **kwargs):
# Sort these by frequency
self.srcList = srcList
Survey.BaseSurvey.__init__(self, **kwargs)
def eval(self, u):
"""
Project fields to receiver locations
:param Fields u: fields object
:rtype: numpy.ndarray
:return: data
"""
data = Survey.Data(self)
for src in self.srcList:
for rx in src.rxList:
data[src, rx] = rx.eval(src, self.mesh, u)
return data
def evalDeriv(self, u):
raise Exception('Use Receivers to project fields deriv.')
SimPEG/EM/FDEM/FDEM.py
+126 -142
View File
@@ -18,9 +18,9 @@ class BaseFDEMProblem(BaseEMProblem):
{\mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{b} - \mathbf{M_{\sigma}^e} \mathbf{e} = \mathbf{s_e}}
if using the E-B formulation (:code:`Problem_e`
or :code:`Problem_b`). Note that in this case, :math:`\mathbf{s_e}` is an integrated quantity.
or :code:`Problem_b`). Note that in this case, :math:`\mathbf{s_e}` is an integrated quantity.
If we write Maxwell's equations in terms of
If we write Maxwell's equations in terms of
\\\(\\\mathbf{h}\\\) and current density \\\(\\\mathbf{j}\\\)
.. math ::
@@ -28,7 +28,7 @@ class BaseFDEMProblem(BaseEMProblem):
\mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{j} + i \omega \mathbf{M_{\mu}^e} \mathbf{h} = \mathbf{s_m} \\\\
\mathbf{C} \mathbf{h} - \mathbf{j} = \mathbf{s_e}
if using the H-J formulation (:code:`Problem_j` or :code:`Problem_h`). Note that here, :math:`\mathbf{s_m}` is an integrated quantity.
if using the H-J formulation (:code:`Problem_j` or :code:`Problem_h`). Note that here, :math:`\mathbf{s_m}` is an integrated quantity.
The problem performs the elimination so that we are solving the system for \\\(\\\mathbf{e},\\\mathbf{b},\\\mathbf{j} \\\) or \\\(\\\mathbf{h}\\\)
"""
@@ -36,88 +36,76 @@ class BaseFDEMProblem(BaseEMProblem):
surveyPair = SurveyFDEM
fieldsPair = Fields
def fields(self, m=None):
def fields(self, m):
"""
Solve the forward problem for the fields.
:param numpy.array m: inversion model (nP,)
:rtype numpy.array:
:return F: forward solution
:return f: forward solution
"""
self.curModel = m
F = self.fieldsPair(self.mesh, self.survey)
f = self.fieldsPair(self.mesh, self.survey)
for freq in self.survey.freqs:
A = self.getA(freq)
rhs = self.getRHS(freq)
Ainv = self.Solver(A, **self.solverOpts)
sol = Ainv * rhs
u = Ainv * rhs
Srcs = self.survey.getSrcByFreq(freq)
ftype = self._fieldType + 'Solution'
F[Srcs, ftype] = sol
f[Srcs, self._solutionType] = u
Ainv.clean()
return F
return f
def Jvec(self, m, v, u=None):
def Jvec(self, m, v, f=None):
"""
Sensitivity times a vector.
:param numpy.array m: inversion model (nP,)
:param numpy.array v: vector which we take sensitivity product with (nP,)
:param SimPEG.EM.FDEM.Fields u: fields object
:param SimPEG.EM.FDEM.Fields u: fields object
:rtype numpy.array:
:return: Jv (ndata,)
:return: Jv (ndata,)
"""
if u is None:
u = self.fields(m)
if f is None:
f = self.fields(m)
self.curModel = m
Jv = self.dataPair(self.survey)
for freq in self.survey.freqs:
A = self.getA(freq) #
Ainv = self.Solver(A, **self.solverOpts)
A = self.getA(freq)
Ainv = self.Solver(A, **self.solverOpts) # create the concept of Ainv (actually a solve)
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 = f[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(f, '_%sDeriv'%rx.projField, None)
df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False)
Jv[src, rx] = rx.evalDeriv(src, self.mesh, f, df_dm_v)
Ainv.clean()
return Utils.mkvc(Jv)
def Jtvec(self, m, v, u=None):
def Jtvec(self, m, v, f=None):
"""
Sensitivity transpose times a vector
:param numpy.array m: inversion model (nP,)
:param numpy.array v: vector which we take adjoint product with (nP,)
:param SimPEG.EM.FDEM.Fields u: fields object
:param SimPEG.EM.FDEM.Fields u: fields object
:rtype numpy.array:
:return: Jv (ndata,)
:return: Jv (ndata,)
"""
if u is None:
u = self.fields(m)
if f is None:
f = self.fields(m)
self.curModel = m
@@ -132,35 +120,31 @@ 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 = f[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
PTv = rx.evalDeriv(src, self.mesh, f, v[src, rx], adjoint=True) # wrt f, need possibility wrt m
df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None)
df_duT, df_dmT = df_duTFun(src, None, PTv, adjoint=True)
df_duTFun = getattr(u, '_%sDeriv_u'%rx.projField, None)
df_duT = df_duTFun(src, 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')
ATinv.clean()
return Utils.mkvc(Jtv)
@@ -170,23 +154,23 @@ class BaseFDEMProblem(BaseEMProblem):
Evaluates the sources for a given frequency and puts them in matrix form
:param float freq: Frequency
:rtype: (numpy.ndarray, numpy.ndarray)
:return: S_m, S_e (nE or nF, nSrc)
:rtype: (numpy.ndarray, numpy.ndarray)
:return: s_m, s_e (nE or nF, nSrc)
"""
Srcs = self.survey.getSrcByFreq(freq)
if self._eqLocs is 'FE':
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':
S_m = np.zeros((self.mesh.nE,len(Srcs)), dtype=complex)
S_e = np.zeros((self.mesh.nF,len(Srcs)), dtype=complex)
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._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)
for i, src in enumerate(Srcs):
smi, sei = src.eval(self)
S_m[:,i] = S_m[:,i] + smi
S_e[:,i] = S_e[:,i] + sei
s_m[:,i] = s_m[:,i] + smi
s_e[:,i] = s_e[:,i] + sei
return S_m, S_e
return s_m, s_e
##########################################################################################
@@ -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)
@@ -223,7 +207,7 @@ class Problem_e(BaseFDEMProblem):
def getA(self, freq):
"""
System matrix
.. math ::
\mathbf{A} = \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{C} + i \omega \mathbf{M^e_{\sigma}}
@@ -239,19 +223,19 @@ 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
.. math ::
\\frac{\mathbf{A}(\mathbf{m}) \mathbf{v}}{d \mathbf{m}} = i \omega \\frac{d \mathbf{M^e_{\sigma}}\mathbf{v} }{d\mathbf{m}}
:param float freq: frequency
:param numpy.ndarray u: solution vector (nE,)
:param float freq: frequency
:param numpy.ndarray u: solution vector (nE,)
:param numpy.ndarray v: vector to take prodct with (nP,) or (nD,) for adjoint
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
:return: derivative of the system matrix times a vector (nP,) or adjoint (nD,)
:return: derivative of the system matrix times a vector (nP,) or adjoint (nD,)
"""
dsig_dm = self.curModel.sigmaDeriv
@@ -264,25 +248,25 @@ class Problem_e(BaseFDEMProblem):
def getRHS(self, freq):
"""
Right hand side for the system
Right hand side for the system
.. math ::
\mathbf{RHS} = \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f}\mathbf{s_m} -i\omega\mathbf{M_e}\mathbf{s_e}
:param float freq: Frequency
:rtype: numpy.ndarray
:rtype: numpy.ndarray
:return: RHS (nE, nSrc)
"""
S_m, S_e = self.getSourceTerm(freq)
s_m, s_e = self.getSourceTerm(freq)
C = self.mesh.edgeCurl
MfMui = self.MfMui
return C.T * (MfMui * S_m) -1j * omega(freq) * S_e
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
Derivative of the right hand side with respect to the model
:param float freq: frequency
:param SimPEG.EM.FDEM.Src src: FDEM source
@@ -294,14 +278,14 @@ class Problem_e(BaseFDEMProblem):
C = self.mesh.edgeCurl
MfMui = self.MfMui
S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint=adjoint)
s_mDeriv, s_eDeriv = src.evalDeriv(self, adjoint=adjoint)
if adjoint:
dRHS = MfMui * (C * v)
return S_mDeriv(dRHS) - 1j * omega(freq) * S_eDeriv(v)
return s_mDeriv(dRHS) - 1j * omega(freq) * s_eDeriv(v)
else:
return C.T * (MfMui * S_mDeriv(v)) -1j * omega(freq) * S_eDeriv(v)
return C.T * (MfMui * s_mDeriv(v)) -1j * omega(freq) * s_eDeriv(v)
class Problem_b(BaseFDEMProblem):
@@ -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
@@ -362,12 +346,12 @@ class Problem_b(BaseFDEMProblem):
.. math ::
\\frac{\mathbf{A}(\mathbf{m}) \mathbf{v}}{d \mathbf{m}} = \mathbf{C} \\frac{\mathbf{M^e_{\sigma}} \mathbf{v}}{d\mathbf{m}}
:param float freq: frequency
:param numpy.ndarray u: solution vector (nF,)
:param float freq: frequency
:param numpy.ndarray u: solution vector (nF,)
:param numpy.ndarray v: vector to take prodct with (nP,) or (nD,) for adjoint
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
:return: derivative of the system matrix times a vector (nP,) or adjoint (nD,)
:return: derivative of the system matrix times a vector (nP,) or adjoint (nD,)
"""
MfMui = self.MfMui
@@ -389,21 +373,21 @@ class Problem_b(BaseFDEMProblem):
def getRHS(self, freq):
"""
Right hand side for the system
Right hand side for the system
.. math ::
\mathbf{RHS} = \mathbf{s_m} + \mathbf{M^e_{\sigma}}^{-1}\mathbf{s_e}
:param float freq: Frequency
:rtype: numpy.ndarray
:rtype: numpy.ndarray
:return: RHS (nE, nSrc)
"""
S_m, S_e = self.getSourceTerm(freq)
s_m, s_e = self.getSourceTerm(freq)
C = self.mesh.edgeCurl
MeSigmaI = self.MeSigmaI
RHS = S_m + C * ( MeSigmaI * S_e )
RHS = s_m + C * ( MeSigmaI * s_e )
if self._makeASymmetric is True:
MfMui = self.MfMui
@@ -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
@@ -424,21 +408,21 @@ class Problem_b(BaseFDEMProblem):
"""
C = self.mesh.edgeCurl
S_m, S_e = src.eval(self)
s_m, s_e = src.eval(self)
MfMui = self.MfMui
if self._makeASymmetric and adjoint:
v = self.MfMui * v
MeSigmaIDeriv = self.MeSigmaIDeriv(S_e)
S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint=adjoint)
MeSigmaIDeriv = self.MeSigmaIDeriv(s_e)
s_mDeriv, s_eDeriv = src.evalDeriv(self, adjoint=adjoint)
if not adjoint:
RHSderiv = C * (MeSigmaIDeriv * v)
SrcDeriv = S_mDeriv(v) + C * (self.MeSigmaI * S_eDeriv(v))
SrcDeriv = s_mDeriv(v) + C * (self.MeSigmaI * s_eDeriv(v))
elif adjoint:
RHSderiv = MeSigmaIDeriv.T * (C.T * v)
SrcDeriv = S_mDeriv(v) + self.MeSigmaI.T * (C.T * S_eDeriv(v))
SrcDeriv = s_mDeriv(v) + self.MeSigmaI.T * (C.T * s_eDeriv(v))
if self._makeASymmetric is True and not adjoint:
return MfMui.T * (SrcDeriv + RHSderiv)
@@ -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
@@ -513,32 +497,32 @@ class Problem_j(BaseFDEMProblem):
\\frac{\mathbf{A(\sigma)} \mathbf{v}}{d \mathbf{m}} = \mathbf{C} \mathbf{M^e_{mu^{-1}}} \mathbf{C^{\\top}} \\frac{d \mathbf{M^f_{\sigma^{-1}}}\mathbf{v} }{d \mathbf{m}}
:param float freq: frequency
:param numpy.ndarray u: solution vector (nF,)
:param float freq: frequency
:param numpy.ndarray u: solution vector (nF,)
:param numpy.ndarray v: vector to take prodct with (nP,) or (nD,) for adjoint
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
:return: derivative of the system matrix times a vector (nP,) or adjoint (nD,)
:return: derivative of the system matrix times a vector (nP,) or adjoint (nD,)
"""
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):
"""
Right hand side for the system
Right hand side for the system
.. math ::
@@ -549,20 +533,20 @@ class Problem_j(BaseFDEMProblem):
:return: RHS
"""
S_m, S_e = self.getSourceTerm(freq)
s_m, s_e = self.getSourceTerm(freq)
C = self.mesh.edgeCurl
MeMuI = self.MeMuI
RHS = C * (MeMuI * S_m) - 1j * omega(freq) * S_e
RHS = C * (MeMuI * s_m) - 1j * omega(freq) * s_e
if self._makeASymmetric is True:
MfRho = self.MfRho
return MfRho.T*RHS
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
Derivative of the right hand side with respect to the model
:param float freq: frequency
:param SimPEG.EM.FDEM.Src src: FDEM source
@@ -574,16 +558,16 @@ class Problem_j(BaseFDEMProblem):
C = self.mesh.edgeCurl
MeMuI = self.MeMuI
S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint=adjoint)
s_mDeriv, s_eDeriv = src.evalDeriv(self, adjoint=adjoint)
if adjoint:
if self._makeASymmetric:
MfRho = self.MfRho
v = MfRho*v
return S_mDeriv(MeMuI.T * (C.T * v)) - 1j * omega(freq) * S_eDeriv(v)
return s_mDeriv(MeMuI.T * (C.T * v)) - 1j * omega(freq) * s_eDeriv(v)
else:
RHSDeriv = C * (MeMuI * S_mDeriv(v)) - 1j * omega(freq) * S_eDeriv(v)
RHSDeriv = C * (MeMuI * s_mDeriv(v)) - 1j * omega(freq) * s_eDeriv(v)
if self._makeASymmetric:
MfRho = self.MfRho
@@ -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,51 +619,51 @@ 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
.. math::
\\frac{\mathbf{A}(\mathbf{m}) \mathbf{v}}{d \mathbf{m}} = \mathbf{C}^{\\top}\\frac{d \mathbf{M^f_{\\rho}}\mathbf{v} }{d\mathbf{m}}
:param float freq: frequency
:param numpy.ndarray u: solution vector (nE,)
:param float freq: frequency
:param numpy.ndarray u: solution vector (nE,)
:param numpy.ndarray v: vector to take prodct with (nP,) or (nD,) for adjoint
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
:return: derivative of the system matrix times a vector (nP,) or adjoint (nD,)
:return: derivative of the system matrix times a vector (nP,) or adjoint (nD,)
"""
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):
"""
Right hand side for the system
Right hand side for the system
.. math ::
\mathbf{RHS} = \mathbf{M^e} \mathbf{s_m} + \mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{s_e}
:param float freq: Frequency
:rtype: numpy.ndarray
:rtype: numpy.ndarray
:return: RHS (nE, nSrc)
"""
S_m, S_e = self.getSourceTerm(freq)
s_m, s_e = self.getSourceTerm(freq)
C = self.mesh.edgeCurl
MfRho = self.MfRho
return S_m + C.T * ( MfRho * S_e )
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
Derivative of the right hand side with respect to the model
:param float freq: frequency
:param SimPEG.EM.FDEM.Src src: FDEM source
@@ -689,17 +673,17 @@ class Problem_h(BaseFDEMProblem):
:return: product of rhs deriv with a vector
"""
_, S_e = src.eval(self)
_, s_e = src.eval(self)
C = self.mesh.edgeCurl
MfRho = self.MfRho
MfRhoDeriv = self.MfRhoDeriv(S_e)
MfRhoDeriv = self.MfRhoDeriv(s_e)
if not adjoint:
RHSDeriv = C.T * (MfRhoDeriv * v)
elif adjoint:
RHSDeriv = MfRhoDeriv.T * (C * v)
S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint=adjoint)
s_mDeriv, s_eDeriv = src.evalDeriv(self, adjoint=adjoint)
return RHSDeriv + S_mDeriv(v) + C.T * (MfRho * S_eDeriv(v))
return RHSDeriv + s_mDeriv(v) + C.T * (MfRho * s_eDeriv(v))
SimPEG/EM/FDEM/FieldsFDEM.py
+730 -345
View File
File diff suppressed because it is too large Load Diff
SimPEG/EM/FDEM/SrcFDEM.py
+164 -126
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,34 +14,34 @@ class BaseSrc(Survey.BaseSrc):
def eval(self, prob):
"""
Evaluate the source terms.
- :math:`S_m` : magnetic source term
- :math:`S_e` : electric source term
Evaluate the source terms.
- :math:`s_m` : magnetic source term
- :math:`s_e` : electric source term
:param Problem prob: FDEM Problem
:rtype: (numpy.ndarray, numpy.ndarray)
:return: tuple with magnetic source term and electric source term
"""
S_m = self.S_m(prob)
S_e = self.S_e(prob)
return S_m, S_e
s_m = self.s_m(prob)
s_e = self.s_e(prob)
return s_m, s_e
def evalDeriv(self, prob, v=None, adjoint=False):
"""
Derivatives of the source terms with respect to the inversion model
- :code:`S_mDeriv` : derivative of the magnetic source term
- :code:`S_eDeriv` : derivative of the electric source term
- :code:`s_mDeriv` : derivative of the magnetic source term
- :code:`s_eDeriv` : derivative of the electric source term
:param Problem prob: FDEM Problem
: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):
def s_m(self, prob):
"""
Magnetic source term
Magnetic source term
:param Problem prob: FDEM Problem
:rtype: numpy.ndarray
@@ -93,9 +93,9 @@ class BaseSrc(Survey.BaseSrc):
"""
return Zero()
def S_e(self, prob):
def s_e(self, prob):
"""
Electric source term
Electric source term
:param Problem prob: FDEM Problem
:rtype: numpy.ndarray
@@ -103,7 +103,7 @@ class BaseSrc(Survey.BaseSrc):
"""
return Zero()
def S_mDeriv(self, prob, v, adjoint = False):
def s_mDeriv(self, prob, v, adjoint = False):
"""
Derivative of magnetic source term with respect to the inversion model
@@ -116,7 +116,7 @@ class BaseSrc(Survey.BaseSrc):
return Zero()
def S_eDeriv(self, prob, v, adjoint = False):
def s_eDeriv(self, prob, v, adjoint = False):
"""
Derivative of electric source term with respect to the inversion model
@@ -131,88 +131,117 @@ class BaseSrc(Survey.BaseSrc):
class RawVec_e(BaseSrc):
"""
RawVec electric source. It is defined by the user provided vector S_e
RawVec electric source. It is defined by the user provided vector s_e
:param list rxList: receiver list
:param float freq: frequency
:param numpy.array S_e: electric 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_e): #, ePrimary=None, bPrimary=None, hPrimary=None, jPrimary=None):
self._S_e = np.array(S_e,dtype=complex)
self.freq = float(freq)
BaseSrc.__init__(self, rxList)
def S_e(self, prob):
return self._S_e
class RawVec_m(BaseSrc):
"""
RawVec magnetic source. It is defined by the user provided vector S_m
:param float freq: frequency
:param rxList: receiver list
:param numpy.array S_m: magnetic source term
"""
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_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_m(self, prob):
def s_e(self, prob):
"""
Magnetic source term
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 RawVec_m(BaseSrc):
"""
RawVec magnetic source. It is defined by the user provided vector s_m
: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)
self.freq = float(freq)
self.integrate = integrate
BaseSrc.__init__(self, rxList)
def s_m(self, prob):
"""
Magnetic source term
:param Problem prob: FDEM Problem
:rtype: numpy.ndarray
:return: magnetic source term on mesh
"""
return self._S_m
if prob._formulation is 'HJ' and self.integrate is True:
return prob.Me * self._s_m
return self._s_m
class RawVec(BaseSrc):
"""
RawVec source. It is defined by the user provided vectors S_m, S_e
RawVec source. It is defined by the user provided vectors s_m, s_e
:param rxList: receiver list
:param float freq: frequency
:param numpy.array S_m: magnetic source term
:param numpy.array S_e: electric source term
: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:
return prob.Me * self._S_m
return self._S_m
def s_m(self, prob):
"""
Magnetic source term
def S_e(self, prob):
if prob._eqLocs is 'FE' and self.integrate is True:
return prob.Me * self._S_e
return self._S_e
: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):
"""
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,44 +332,46 @@ 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):
def s_m(self, prob):
"""
The magnetic source term
: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):
def s_e(self, prob):
"""
The electric source term
: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,42 +448,44 @@ 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):
def s_m(self, prob):
"""
The magnetic source term
: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):
def s_e(self, prob):
"""
The electric source term
: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,44 +560,50 @@ 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
def S_m(self, prob):
def s_m(self, prob):
"""
The magnetic source term
: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):
def s_e(self, prob):
"""
The electric source term
: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
+46 -54
View File
@@ -1,8 +1,10 @@
import SimPEG
from SimPEG.EM.Utils import *
from SimPEG.EM.Base import BaseEMSurvey
from scipy.constants import mu_0
from SimPEG.Utils import Zero, Identity
import SrcFDEM as Src
from SimPEG import sp
####################################################
@@ -18,33 +20,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,16 +58,15 @@ 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 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, f):
"""
Project fields to recievers to get data.
@@ -76,11 +77,16 @@ 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
u_part = getattr(u_part_complex, real_or_imag)
return P*u_part
# projGLoc = u._GLoc(self.knownRxTypes[self.rxType][0])
# projGLoc += self.knownRxTypes[self.rxType][1]
P = self.getP(mesh, self.projGLoc(f))
f_part_complex = f[src, self.projField]
# get the real or imag component
real_or_imag = self.projComp
f_part = getattr(f_part_complex, real_or_imag)
return P*f_part
def evalDeriv(self, src, mesh, f, v, adjoint=False):
"""
@@ -93,7 +99,8 @@ class Rx(SimPEG.Survey.BaseRx):
:rtype: numpy.ndarray
:return: fields projected to recievers
"""
P = self.getP(mesh)
P = self.getP(mesh, self.projGLoc(f))
if not adjoint:
Pv_complex = P * v
@@ -117,7 +124,7 @@ class Rx(SimPEG.Survey.BaseRx):
# Survey
####################################################
class Survey(SimPEG.Survey.BaseSurvey):
class Survey(BaseEMSurvey):
"""
Frequency domain electromagnetic survey
@@ -125,12 +132,12 @@ class Survey(SimPEG.Survey.BaseSurvey):
"""
srcPair = Src.BaseSrc
rxPaair = Rx
rxPair = Rx
def __init__(self, srcList, **kwargs):
# Sort these by frequency
self.srcList = srcList
SimPEG.Survey.BaseSurvey.__init__(self, **kwargs)
BaseEMSurvey.__init__(self, srcList, **kwargs)
_freqDict = {}
for src in srcList:
@@ -165,23 +172,8 @@ class Survey(SimPEG.Survey.BaseSurvey):
Returns the sources associated with a specific frequency.
:param float freq: frequency for which we look up sources
:rtype: dictionary
:return: sources at the sepcified frequency
:return: sources at the sepcified frequency
"""
assert freq in self._freqDict, "The requested frequency is not in this survey."
return self._freqDict[freq]
def eval(self, u):
"""
Project fields to receiver locations
:param Fields u: fields object
:rtype: numpy.ndarray
:return: data
"""
data = SimPEG.Survey.Data(self)
for src in self.srcList:
for rx in src.rxList:
data[src, rx] = rx.eval(src, self.mesh, u)
return data
def evalDeriv(self, u):
raise Exception('Use Receivers to project fields deriv.')
SimPEG/EM/TDEM/BaseTDEM.py
+11 -11
View File
@@ -108,11 +108,11 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
Ainv.clean()
return F
def Jvec(self, m, v, u=None):
def Jvec(self, m, v, f=None):
"""
:param numpy.array m: Conductivity model
:param numpy.ndarray v: vector (model object)
:param simpegEM.TDEM.FieldsTDEM u: Fields resulting from m
:param simpegEM.TDEM.FieldsTDEM f: Fields resulting from m
:rtype: numpy.ndarray
:return: w (data object)
@@ -125,15 +125,15 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
"""
if self.verbose: print '%s\nCalculating J(v)\n%s'%('*'*50,'*'*50)
self.curModel = m
if u is None:
u = self.fields(m)
p = self.Gvec(m, v, u)
if f is None:
f = self.fields(m)
p = self.Gvec(m, v, f)
y = self.solveAh(m, p)
Jv = self.survey.evalDeriv(u, v=y)
Jv = self.survey.evalDeriv(f, v=y)
if self.verbose: print '%s\nDone calculating J(v)\n%s'%('*'*50,'*'*50)
return - mkvc(Jv)
def Jtvec(self, m, v, u=None):
def Jtvec(self, m, v, f=None):
"""
:param numpy.array m: Conductivity model
:param numpy.ndarray,SimPEG.Survey.Data v: vector (data object)
@@ -150,15 +150,15 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
"""
if self.verbose: print '%s\nCalculating J^T(v)\n%s'%('*'*50,'*'*50)
self.curModel = m
if u is None:
u = self.fields(m)
if f is None:
f = self.fields(m)
if not isinstance(v, self.dataPair):
v = self.dataPair(self.survey, v)
p = self.survey.evalDeriv(u, v=v, adjoint=True)
p = self.survey.evalDeriv(f, v=v, adjoint=True)
y = self.solveAht(m, p)
w = self.Gtvec(m, y, u)
w = self.Gtvec(m, y, f)
if self.verbose: print '%s\nDone calculating J^T(v)\n%s'%('*'*50,'*'*50)
return - mkvc(w)
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/Examples/EM_Schenkel_Morrison_Casing.py
+275
View File
@@ -0,0 +1,275 @@
from SimPEG import *
from SimPEG.EM import FDEM, Analytics, mu_0
import time
try:
from pymatsolver import MumpsSolver
solver = MumpsSolver
except Exception:
solver = SolverLU
pass
def run(plotIt=True):
"""
EM: Schenkel and Morrison Casing Model
======================================
Here we create and run a FDEM forward simulation to calculate the vertical
current inside a steel-cased. The model is based on the Schenkel and
Morrison Casing Model, and the results are used in a 2016 SEG abstract by
Yang et al.
- Schenkel, C.J., and H.F. Morrison, 1990, Effects of well casing on potential field measurements using downhole current sources: Geophysical prospecting, 38, 663-686.
The model consists of:
- Air: Conductivity 1e-8 S/m, above z = 0
- Background: conductivity 1e-2 S/m, below z = 0
- Casing: conductivity 1e6 S/m
- 300m long
- radius of 0.1m
- thickness of 6e-3m
Inside the casing, we take the same conductivity as the background.
We are using an EM code to simulate DC, so we use frequency low enough
that the skin depth inside the casing is longer than the casing length (f
= 1e-6 Hz). The plot produced is of the current inside the casing.
These results are shown in the SEG abstract by Yang et al., 2016: 3D DC
resistivity modeling of steel casing for reservoir monitoring using
equivalent resistor network. The solver used to produce these results and
achieve the CPU time of ~30s is Mumps, which was installed using pymatsolver_
.. _pymatsolver: https://github.com/rowanc1/pymatsolver
This example is on figshare: https://dx.doi.org/10.6084/m9.figshare.3126961.v1
If you would use this example for a code comparison, or build upon it, a
citation would be much appreciated!
"""
if plotIt:
import matplotlib.pylab as plt
# ------------------ MODEL ------------------
sigmaair = 1e-8 # air
sigmaback = 1e-2 # background
sigmacasing = 1e6 # casing
sigmainside = sigmaback # inside the casing
casing_t = 0.006 # 1cm thickness
casing_l = 300 # length of the casing
casing_r = 0.1
casing_a = casing_r - casing_t/2. # inner radius
casing_b = casing_r + casing_t/2. # outer radius
casing_z = np.r_[-casing_l,0.]
# ------------------ SURVEY PARAMETERS ------------------
freqs = np.r_[1e-6] #[1e-1, 1, 5] # frequencies
dsz = -300 # down-hole z source location
src_loc = np.r_[0.,0.,dsz]
inf_loc = np.r_[0.,0.,1e4]
print 'Skin Depth: ', [(500./np.sqrt(sigmaback*_)) for _ in freqs]
# ------------------ MESH ------------------
# fine cells near well bore
csx1, csx2 = 2e-3, 60.
pfx1, pfx2 = 1.3, 1.3
ncx1 = np.ceil(casing_b/csx1+2)
# pad nicely to second cell size
npadx1 = np.floor(np.log(csx2/csx1) / np.log(pfx1))
hx1a,hx1b = Utils.meshTensor([(csx1,ncx1)]),Utils.meshTensor([(csx1,npadx1,pfx1)])
dx1 = sum(hx1a)+sum(hx1b)
dx1 = np.floor(dx1/csx2)
hx1b *= (dx1*csx2 - sum(hx1a))/sum(hx1b)
# second chunk of mesh
dx2 = 300. # uniform mesh out to here
ncx2 = np.ceil((dx2 - dx1)/csx2)
npadx2 = 45
hx2a, hx2b = Utils.meshTensor([(csx2,ncx2)]), Utils.meshTensor([(csx2,npadx2,pfx2)])
hx = np.hstack([hx1a,hx1b,hx2a,hx2b])
# z-direction
csz = 0.05
nza = 10
ncz, npadzu, npadzd = np.int(np.ceil(np.diff(casing_z)[0]/csz))+10, 68, 68 # cell size, number of core cells, number of padding cells in the x- direction
hz = Utils.meshTensor([(csz,npadzd,-1.3), (csz,ncz), (csz,npadzu,1.3)]) # vector of cell widths in the z-direction
# Mesh
mesh = Mesh.CylMesh([hx,1.,hz], [0.,0.,-np.sum(hz[:npadzu+ncz-nza])])
print 'Mesh Extent xmax: %f,: zmin: %f, zmax: %f'%(mesh.vectorCCx.max(), mesh.vectorCCz.min(), mesh.vectorCCz.max())
print 'Number of cells', mesh.nC
if plotIt is True:
fig, ax = plt.subplots(1, 1, figsize=(6, 4))
ax.set_title('Simulation Mesh')
mesh.plotGrid(ax=ax)
plt.show()
# Put the model on the mesh
sigWholespace = sigmaback*np.ones((mesh.nC))
sigBack = sigWholespace.copy()
sigBack[mesh.gridCC[:,2] > 0.] = sigmaair
sigCasing = sigBack.copy()
iCasingZ = (mesh.gridCC[:,2] <= casing_z[1]) & (mesh.gridCC[:,2] >= casing_z[0])
iCasingX = (mesh.gridCC[:,0] >= casing_a) & (mesh.gridCC[:,0] <= casing_b)
iCasing = iCasingX & iCasingZ
sigCasing[iCasing] = sigmacasing
if plotIt is True:
# plotting parameters
xlim = np.r_[0., 0.2]
zlim = np.r_[-350., 10.]
clim_sig = np.r_[-8,6]
# plot models
fig, ax = plt.subplots(1,1,figsize=(4,4))
f = plt.colorbar(mesh.plotImage(np.log10(sigCasing),ax=ax)[0], ax=ax)
ax.grid(which='both')
ax.set_title('Log_10 (Sigma)')
ax.set_xlim(xlim)
ax.set_ylim(zlim)
f.set_clim(clim_sig)
plt.show()
# -------------- Sources --------------------
# Define Custom Current Sources
# surface source
sg_x = np.zeros(mesh.vnF[0],dtype=complex)
sg_y = np.zeros(mesh.vnF[1],dtype=complex)
sg_z = np.zeros(mesh.vnF[2],dtype=complex)
nza = 2 # put the wire two cells above the surface
ncin = 2
# vertically directed wire
sgv_indx = (mesh.gridFz[:,0] > casing_a) & (mesh.gridFz[:,0] < casing_a + csx1) # hook it up to casing at the surface
sgv_indz = (mesh.gridFz[:,2] <= +csz*nza) & (mesh.gridFz[:,2] >= -csz*2)
sgv_ind = sgv_indx & sgv_indz
sg_z[sgv_ind] = -1.
# horizontally directed wire
sgh_indx = (mesh.gridFx[:,0] > casing_a) & (mesh.gridFx[:,0] <= inf_loc[2])
sgh_indz = (mesh.gridFx[:,2] > csz*(nza-0.5)) & (mesh.gridFx[:,2] < csz*(nza+0.5))
sgh_ind = sgh_indx & sgh_indz
sg_x[sgh_ind] = -1.
sgv2_indx = (mesh.gridFz[:,0] >= mesh.gridFx[sgh_ind,0].max()) & (mesh.gridFz[:,0] <= inf_loc[2]*1.2) # hook it up to casing at the surface
sgv2_indz = (mesh.gridFz[:,2] <= +csz*nza) & (mesh.gridFz[:,2] >= -csz*2)
sgv2_ind = sgv2_indx & sgv2_indz
sg_z[sgv2_ind] = 1.
# assemble the source
sg = np.hstack([sg_x,sg_y,sg_z])
sg_p = [FDEM.Src.RawVec_e([],_,sg/mesh.area) for _ in freqs]
# downhole source
dg_x = np.zeros(mesh.vnF[0],dtype=complex)
dg_y = np.zeros(mesh.vnF[1],dtype=complex)
dg_z = np.zeros(mesh.vnF[2],dtype=complex)
# vertically directed wire
dgv_indx = (mesh.gridFz[:,0] < csx1) # go through the center of the well
dgv_indz = (mesh.gridFz[:,2] <= +csz*nza) & (mesh.gridFz[:,2] > dsz + csz/2.)
dgv_ind = dgv_indx & dgv_indz
dg_z[dgv_ind] = -1.
# couple to the casing downhole
dgh_indx = mesh.gridFx[:,0] < casing_a + csx1
dgh_indz = (mesh.gridFx[:,2] < dsz + csz) & (mesh.gridFx[:,2] >= dsz)
dgh_ind = dgh_indx & dgh_indz
dg_x[dgh_ind] = 1.
# horizontal part at surface
dgh2_indx = mesh.gridFx[:,0] <= inf_loc[2]*1.2
dgh2_indz = sgh_indz.copy()
dgh2_ind = dgh2_indx & dgh2_indz
dg_x[dgh2_ind] = -1.
# vertical part at surface
dgv2_ind = sgv2_ind.copy()
dg_z[dgv2_ind] = 1.
# assemble the source
dg = np.hstack([dg_x,dg_y,dg_z])
dg_p = [FDEM.Src.RawVec_e([],_,dg/mesh.area) for _ in freqs]
# ------------ Problem and Survey ---------------
survey = FDEM.Survey(sg_p + dg_p)
mapping = [('sigma', Maps.IdentityMap(mesh))]
problem = FDEM.Problem_h(mesh, mapping=mapping)
problem.pair(survey)
# ------------- Solve ---------------------------
t0 = time.time()
fieldsCasing = problem.fields(sigCasing)
print 'Time to solve 2 sources', time.time() - t0
# Plot current
# current density
jn0 = fieldsCasing[dg_p,'j']
jn1 = fieldsCasing[sg_p,'j']
# current
in0 = [mesh.area*fieldsCasing[dg_p,'j'][:,i] for i in range(len(freqs))]
in1 = [mesh.area*fieldsCasing[sg_p,'j'][:,i] for i in range(len(freqs))]
in0 = np.vstack(in0).T
in1 = np.vstack(in1).T
# integrate to get z-current inside casing
inds_inx = (mesh.gridFz[:,0] >= casing_a) & (mesh.gridFz[:,0] <= casing_b)
inds_inz = (mesh.gridFz[:,2] >= dsz ) & (mesh.gridFz[:,2] <= 0)
inds_fz = inds_inx & inds_inz
indsx = [False]*mesh.nFx
inds = list(indsx) + list(inds_fz)
in0_in = in0[np.r_[inds]]
in1_in = in1[np.r_[inds]]
z_in = mesh.gridFz[inds_fz,2]
in0_in = in0_in.reshape([in0_in.shape[0]/3,3])
in1_in = in1_in.reshape([in1_in.shape[0]/3,3])
z_in = z_in.reshape([z_in.shape[0]/3,3])
I0 = in0_in.sum(1).real
I1 = in1_in.sum(1).real
z_in = z_in[:,0]
if plotIt is True:
fig, ax = plt.subplots(1,2,figsize=(12,4))
ax[0].plot(z_in,np.absolute(I0), z_in,np.absolute(I1))
ax[0].legend(['top casing', 'bottom casing'],loc='best')
ax[0].set_title('Magnitude of Vertical Current in Casing')
ax[1].semilogy(z_in,np.absolute(I0), z_in,np.absolute(I1))
ax[1].legend(['top casing', 'bottom casing'],loc='best')
ax[1].set_title('Magnitude of Vertical Current in Casing')
ax[1].set_ylim([1e-2, 1.])
plt.show()
if __name__ == '__main__':
run()
SimPEG/Examples/__init__.py
+2 -1
View File
@@ -5,6 +5,7 @@ import DC_Analytic_Dipole
import DC_Forward_PseudoSection
import EM_FDEM_1D_Inversion
import EM_FDEM_Analytic_MagDipoleWholespace
import EM_Schenkel_Morrison_Casing
import EM_TDEM_1D_Inversion
import FLOW_Richards_1D_Celia1990
import Forward_BasicDirectCurrent
@@ -19,7 +20,7 @@ import Mesh_Tensor_Creation
import MT_1D_ForwardAndInversion
import MT_3D_Foward
__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Forward_BasicDirectCurrent", "Inversion_Linear", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation", "MT_1D_ForwardAndInversion", "MT_3D_Foward"]
__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_Schenkel_Morrison_Casing", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Forward_BasicDirectCurrent", "Inversion_Linear", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation", "MT_1D_ForwardAndInversion", "MT_3D_Foward"]
##### AUTOIMPORTS #####
SimPEG/FLOW/Richards/RichardsProblem.py
+29 -29
View File
@@ -45,19 +45,19 @@ class RichardsSurvey(Survey.BaseSurvey):
@Utils.count
@Utils.requires('prob')
def dpred(self, m, u=None):
def dpred(self, m, f=None):
"""
Create the projected data from a model.
The field, u, (if provided) will be used for the predicted data
The field, f, (if provided) will be used for the predicted data
instead of recalculating the fields (which may be expensive!).
.. math::
d_\\text{pred} = P(u(m), m)
d_\\text{pred} = P(f(m), m)
Where P is a projection of the fields onto the data space.
"""
if u is None: u = self.prob.fields(m)
return Utils.mkvc(self.eval(u, m))
if f is None: f = self.prob.fields(m)
return Utils.mkvc(self.eval(f, m))
@Utils.requires('prob')
def eval(self, U, m):
@@ -233,16 +233,16 @@ class RichardsProblem(Problem.BaseTimeProblem):
return r, J
@Utils.timeIt
def Jfull(self, m, u=None):
if u is None:
u = self.fields(m)
def Jfull(self, m, f=None):
if f is None:
f = self.fields(m)
nn = len(u)-1
nn = len(f)-1
Asubs, Adiags, Bs = range(nn), range(nn), range(nn)
for ii in range(nn):
dt = self.timeSteps[ii]
bc = self.getBoundaryConditions(ii, u[ii])
Asubs[ii], Adiags[ii], Bs[ii] = self.diagsJacobian(m, u[ii], u[ii+1], dt, bc)
bc = self.getBoundaryConditions(ii, f[ii])
Asubs[ii], Adiags[ii], Bs[ii] = self.diagsJacobian(m, f[ii], f[ii+1], dt, bc)
Ad = sp.block_diag(Adiags)
zRight = Utils.spzeros((len(Asubs)-1)*Asubs[0].shape[0],Adiags[0].shape[1])
zTop = Utils.spzeros(Adiags[0].shape[0], len(Adiags)*Adiags[0].shape[1])
@@ -251,7 +251,7 @@ class RichardsProblem(Problem.BaseTimeProblem):
B = np.array(sp.vstack(Bs).todense())
Ainv = self.Solver(A, **self.solverOpts)
P = self.survey.evalDeriv(u, m)
P = self.survey.evalDeriv(f, m)
AinvB = Ainv * B
z = np.zeros((self.mesh.nC, B.shape[1]))
zAinvB = np.vstack((z, AinvB))
@@ -259,41 +259,41 @@ class RichardsProblem(Problem.BaseTimeProblem):
return J
@Utils.timeIt
def Jvec(self, m, v, u=None):
if u is None:
u = self.fields(m)
def Jvec(self, m, v, f=None):
if f is None:
f = self.fields(m)
JvC = range(len(u)-1) # Cell to hold each row of the long vector.
JvC = range(len(f)-1) # Cell to hold each row of the long vector.
# This is done via forward substitution.
bc = self.getBoundaryConditions(0, u[0])
temp, Adiag, B = self.diagsJacobian(m, u[0], u[1], self.timeSteps[0], bc)
bc = self.getBoundaryConditions(0, f[0])
temp, Adiag, B = self.diagsJacobian(m, f[0], f[1], self.timeSteps[0], bc)
Adiaginv = self.Solver(Adiag, **self.solverOpts)
JvC[0] = Adiaginv * (B*v)
for ii in range(1,len(u)-1):
bc = self.getBoundaryConditions(ii, u[ii])
Asub, Adiag, B = self.diagsJacobian(m, u[ii], u[ii+1], self.timeSteps[ii], bc)
for ii in range(1,len(f)-1):
bc = self.getBoundaryConditions(ii, f[ii])
Asub, Adiag, B = self.diagsJacobian(m, f[ii], f[ii+1], self.timeSteps[ii], bc)
Adiaginv = self.Solver(Adiag, **self.solverOpts)
JvC[ii] = Adiaginv * (B*v - Asub*JvC[ii-1])
P = self.survey.evalDeriv(u, m)
P = self.survey.evalDeriv(f, m)
return P * np.concatenate([np.zeros(self.mesh.nC)] + JvC)
@Utils.timeIt
def Jtvec(self, m, v, u=None):
if u is None:
u = self.field(m)
def Jtvec(self, m, v, f=None):
if f is None:
f = self.field(m)
P = self.survey.evalDeriv(u, m)
P = self.survey.evalDeriv(f, m)
PTv = P.T*v
# This is done via backward substitution.
minus = 0
BJtv = 0
for ii in range(len(u)-1,0,-1):
bc = self.getBoundaryConditions(ii-1, u[ii-1])
Asub, Adiag, B = self.diagsJacobian(m, u[ii-1], u[ii], self.timeSteps[ii-1], bc)
for ii in range(len(f)-1,0,-1):
bc = self.getBoundaryConditions(ii-1, f[ii-1])
Asub, Adiag, B = self.diagsJacobian(m, f[ii-1], f[ii], self.timeSteps[ii-1], bc)
#select the correct part of v
vpart = range((ii)*Adiag.shape[0], (ii+1)*Adiag.shape[0])
AdiaginvT = self.Solver(Adiag.T, **self.solverOpts)
SimPEG/InvProblem.py
+13 -13
View File
@@ -82,23 +82,23 @@ class BaseInvProblem(object):
self._warmstart = value
def getFields(self, m, store=False, deleteWarmstart=True):
u = None
f = None
for mtest, u_ofmtest in self.warmstart:
if m is mtest:
u = u_ofmtest
f = u_ofmtest
if self.debug: print 'InvProb is Warm Starting!'
break
if u is None:
u = self.prob.fields(m)
if f is None:
f = self.prob.fields(m)
if deleteWarmstart:
self.warmstart = []
if store:
self.warmstart += [(m,u)]
self.warmstart += [(m,f)]
return u
return f
@Utils.timeIt
def evalFunction(self, m, return_g=True, return_H=True):
@@ -109,21 +109,21 @@ class BaseInvProblem(object):
gc.collect()
# Store fields if doing a line-search
u = self.getFields(m, store=(return_g==False and return_H==False))
f = self.getFields(m, store=(return_g==False and return_H==False))
phi_d = self.dmisfit.eval(m, u=u)
phi_d = self.dmisfit.eval(m, f=f)
phi_m = self.reg.eval(m)
self.dpred = self.survey.dpred(m, u=u) # This is a cheap matrix vector calculation.
self.dpred = self.survey.dpred(m, f=f) # This is a cheap matrix vector calculation.
self.phi_d, self.phi_d_last = phi_d, self.phi_d
self.phi_m, self.phi_m_last = phi_m, self.phi_m
f = phi_d + self.beta * phi_m
phi = phi_d + self.beta * phi_m
out = (f,)
out = (phi,)
if return_g:
phi_dDeriv = self.dmisfit.evalDeriv(m, u=u)
phi_dDeriv = self.dmisfit.evalDeriv(m, f=f)
phi_mDeriv = self.reg.evalDeriv(m)
g = phi_dDeriv + self.beta * phi_mDeriv
@@ -131,7 +131,7 @@ class BaseInvProblem(object):
if return_H:
def H_fun(v):
phi_d2Deriv = self.dmisfit.eval2Deriv(m, v, u=u)
phi_d2Deriv = self.dmisfit.eval2Deriv(m, v, f=f)
phi_m2Deriv = self.reg.eval2Deriv(m, v=v)
return phi_d2Deriv + self.beta * phi_m2Deriv
SimPEG/MT/BaseMT.py
+13 -13
View File
@@ -27,7 +27,7 @@ class BaseMTProblem(BaseFDEMProblem):
# Might need to add more stuff here.
## NEED to clean up the Jvec and Jtvec to use Zero and Identities for None components.
def Jvec(self, m, v, u=None):
def Jvec(self, m, v, f=None):
"""
Function to calculate the data sensitivities dD/dm times a vector.
@@ -39,8 +39,8 @@ class BaseMTProblem(BaseFDEMProblem):
"""
# Calculate the fields
if u is None:
u = self.fields(m)
if f is None:
f= self.fields(m)
# Set current model
self.curModel = m
# Initiate the Jv object
@@ -56,9 +56,9 @@ class BaseMTProblem(BaseFDEMProblem):
# We need fDeriv_m = df/du*du/dm + df/dm
# Construct du/dm, it requires a solve
# NOTE: need to account for the 2 polarizations in the derivatives.
u_src = u[src,:]
f_src = f[src,:]
# dA_dm and dRHS_dm should be of size nE,2, so that we can multiply by dA_duI. The 2 columns are each of the polarizations.
dA_dm = self.getADeriv_m(freq, u_src, v) # Size: nE,2 (u_px,u_py) in the columns.
dA_dm = self.getADeriv_m(freq, f_src, v) # Size: nE,2 (u_px,u_py) in the columns.
dRHS_dm = self.getRHSDeriv_m(freq, v) # Size: nE,2 (u_px,u_py) in the columns.
if dRHS_dm is None:
du_dm = dA_duI * ( -dA_dm )
@@ -68,13 +68,13 @@ class BaseMTProblem(BaseFDEMProblem):
for rx in src.rxList:
# Get the projection derivative
# v should be of size 2*nE (for 2 polarizations)
PDeriv_u = lambda t: rx.evalDeriv(src, self.mesh, u, t) # wrt u, we don't have have PDeriv wrt m
PDeriv_u = lambda t: rx.evalDeriv(src, self.mesh, f, t) # wrt u, we don't have have PDeriv wrt m
Jv[src, rx] = PDeriv_u(mkvc(du_dm))
dA_duI.clean()
# Return the vectorized sensitivities
return mkvc(Jv)
def Jtvec(self, m, v, u=None):
def Jtvec(self, m, v, f=None):
"""
Function to calculate the transpose of the data sensitivities (dD/dm)^T times a vector.
@@ -85,8 +85,8 @@ class BaseMTProblem(BaseFDEMProblem):
:return: Data sensitivities wrt m
"""
if u is None:
u = self.fields(m)
if f is None:
f = self.fields(m)
self.curModel = m
@@ -103,15 +103,15 @@ class BaseMTProblem(BaseFDEMProblem):
for src in self.survey.getSrcByFreq(freq):
ftype = self._fieldType + 'Solution'
u_src = u[src, :]
f_src = f[src, :]
for rx in src.rxList:
# Get the adjoint evalDeriv
# PTv needs to be nE,
PTv = rx.evalDeriv(src, self.mesh, u, mkvc(v[src, rx],2), adjoint=True) # wrt u, need possibility wrt m
PTv = rx.evalDeriv(src, self.mesh, f, mkvc(v[src, rx],2), adjoint=True) # wrt u, need possibility wrt m
# Get the
dA_duIT = ATinv * PTv
dA_dmT = self.getADeriv_m(freq, u_src, mkvc(dA_duIT), adjoint=True)
dA_dmT = self.getADeriv_m(freq, f_src, mkvc(dA_duIT), adjoint=True)
dRHS_dmT = self.getRHSDeriv_m(freq, mkvc(dA_duIT), adjoint=True)
# Make du_dmT
if dRHS_dmT is None:
@@ -129,4 +129,4 @@ class BaseMTProblem(BaseFDEMProblem):
raise Exception('Must be real or imag')
# Clean the factorization, clear memory.
ATinv.clean()
return Jtv
return Jtv
SimPEG/MT/SurveyMT.py
+3 -3
View File
@@ -427,15 +427,15 @@ class Survey(SimPEGsurvey.BaseSurvey):
assert freq in self._freqDict, "The requested frequency is not in this survey."
return self._freqDict[freq]
def eval(self, u):
def eval(self, f):
data = Data(self)
for src in self.srcList:
sys.stdout.flush()
for rx in src.rxList:
data[src, rx] = rx.eval(src, self.mesh, u)
data[src, rx] = rx.eval(src, self.mesh, f)
return data
def evalDeriv(self, u):
def evalDeriv(self, f):
raise Exception('Use Transmitters to project fields deriv.')
#################
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/Optimization.py
+17
View File
@@ -888,6 +888,8 @@ class ProjectedGNCG(BFGS, Minimize, Remember):
maxIterCG = 5
tolCG = 1e-1
stepOffBoundsFact = 0.1 # perturbation of the inactive set off the bounds
lower = -np.inf
upper = np.inf
@@ -990,4 +992,19 @@ class ProjectedGNCG(BFGS, Minimize, Remember):
cgFlag = 1
# End CG Iterations
# Take a gradient step on the active cells if exist
if temp != self.xc.size:
rhs_a = (Active) * -self.g
dm_i = max( abs( delx ) )
dm_a = max( abs(rhs_a) )
# perturb inactive set off of bounds so that they are included in the step
delx = delx + self.stepOffBoundsFact * (rhs_a * dm_i / dm_a)
# Only keep gradients going in the right direction on the active set
indx = ((self.xc<=self.lower) & (delx < 0)) | ((self.xc>=self.upper) & (delx > 0))
delx[indx] = 0.
return delx
SimPEG/Problem.py
+16 -16
View File
@@ -88,28 +88,28 @@ class BaseProblem(object):
return self.survey is not None
@Utils.timeIt
def Jvec(self, m, v, u=None):
"""Jvec(m, v, u=None)
def Jvec(self, m, v, f=None):
"""Jvec(m, v, f=None)
Effect of J(m) on a vector v.
:param numpy.array m: model
:param numpy.array v: vector to multiply
:param numpy.array u: fields
:param Fields f: fields
:rtype: numpy.array
:return: Jv
"""
raise NotImplementedError('J is not yet implemented.')
@Utils.timeIt
def Jtvec(self, m, v, u=None):
"""Jtvec(m, v, u=None)
def Jtvec(self, m, v, f=None):
"""Jtvec(m, v, f=None)
Effect of transpose of J(m) on a vector v.
:param numpy.array m: model
:param numpy.array v: vector to multiply
:param numpy.array u: fields
:param Fields f: fields
:rtype: numpy.array
:return: JTv
"""
@@ -117,32 +117,32 @@ class BaseProblem(object):
@Utils.timeIt
def Jvec_approx(self, m, v, u=None):
"""Jvec_approx(m, v, u=None)
def Jvec_approx(self, m, v, f=None):
"""Jvec_approx(m, v, f=None)
Approximate effect of J(m) on a vector v
:param numpy.array m: model
:param numpy.array v: vector to multiply
:param numpy.array u: fields
:param Fields f: fields
:rtype: numpy.array
:return: approxJv
"""
return self.Jvec(m, v, u)
return self.Jvec(m, v, f)
@Utils.timeIt
def Jtvec_approx(self, m, v, u=None):
"""Jtvec_approx(m, v, u=None)
def Jtvec_approx(self, m, v, f=None):
"""Jtvec_approx(m, v, f=None)
Approximate effect of transpose of J(m) on a vector v.
:param numpy.array m: model
:param numpy.array v: vector to multiply
:param numpy.array u: fields
:param Fields f: fields
:rtype: numpy.array
:return: JTv
"""
return self.Jtvec(m, v, u)
return self.Jtvec(m, v, f)
def fields(self, m):
"""
@@ -224,9 +224,9 @@ class LinearProblem(BaseProblem):
def fields(self, m):
return self.G.dot(m)
def Jvec(self, m, v, u=None):
def Jvec(self, m, v, f=None):
return self.G.dot(v)
def Jtvec(self, m, v, u=None):
def Jtvec(self, m, v, f=None):
return self.G.T.dot(v)
SimPEG/Survey.py
+25 -22
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
@@ -292,38 +295,38 @@ class BaseSurvey(object):
@Utils.count
@Utils.requires('prob')
def dpred(self, m, u=None):
"""dpred(m, u=None)
def dpred(self, m, f=None):
"""dpred(m, f=None)
Create the projected data from a model.
The field, u, (if provided) will be used for the predicted data
The fields, f, (if provided) will be used for the predicted data
instead of recalculating the fields (which may be expensive!).
.. math::
d_\\text{pred} = P(u(m))
d_\\text{pred} = P(f(m))
Where P is a projection of the fields onto the data space.
"""
if u is None: u = self.prob.fields(m)
return Utils.mkvc(self.eval(u))
if f is None: f = self.prob.fields(m)
return Utils.mkvc(self.eval(f))
@Utils.count
def eval(self, u):
"""eval(u)
def eval(self, f):
"""eval(f)
This function projects the fields onto the data space.
.. math::
d_\\text{pred} = \mathbf{P} u(m)
d_\\text{pred} = \mathbf{P} f(m)
"""
raise NotImplemented('eval is not yet implemented.')
@Utils.count
def evalDeriv(self, u):
"""evalDeriv(u)
def evalDeriv(self, f):
"""evalDeriv(f)
This function s the derivative of projects the fields onto the data space.
@@ -334,11 +337,11 @@ class BaseSurvey(object):
raise NotImplemented('eval is not yet implemented.')
@Utils.count
def residual(self, m, u=None):
"""residual(m, u=None)
def residual(self, m, f=None):
"""residual(m, f=None)
:param numpy.array m: geophysical model
:param numpy.array u: fields
:param numpy.array f: fields
:rtype: numpy.array
:return: data residual
@@ -349,14 +352,14 @@ class BaseSurvey(object):
\mu_\\text{data} = \mathbf{d}_\\text{pred} - \mathbf{d}_\\text{obs}
"""
return Utils.mkvc(self.dpred(m, u=u) - self.dobs)
return Utils.mkvc(self.dpred(m, f=f) - self.dobs)
@property
def isSynthetic(self):
"Check if the data is synthetic."
return self.mtrue is not None
def makeSyntheticData(self, m, std=0.05, u=None, force=False):
def makeSyntheticData(self, m, std=0.05, f=None, force=False):
"""
Make synthetic data given a model, and a standard deviation.
@@ -369,16 +372,16 @@ class BaseSurvey(object):
if getattr(self, 'dobs', None) is not None and not force:
raise Exception('Survey already has dobs. You can use force=True to override this exception.')
self.mtrue = m
self.dtrue = self.dpred(m, u=u)
self.dtrue = self.dpred(m, f=f)
noise = std*abs(self.dtrue)*np.random.randn(*self.dtrue.shape)
self.dobs = self.dtrue+noise
self.std = self.dobs*0 + std
return self.dobs
class LinearSurvey(BaseSurvey):
def eval(self, u):
return u
def eval(self, f):
return f
@property
def nD(self):
return self.prob.G.shape[0]
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'
docs/conf.py
+2 -2
View File
@@ -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.
docs/examples/EM_Schenkel_Morrison_Casing.rst
+58
View File
@@ -0,0 +1,58 @@
.. _examples_EM_Schenkel_Morrison_Casing:
.. --------------------------------- ..
.. ..
.. THIS FILE IS AUTO GENEREATED ..
.. ..
.. SimPEG/Examples/__init__.py ..
.. ..
.. --------------------------------- ..
EM: Schenkel and Morrison Casing Model
======================================
Here we create and run a FDEM forward simulation to calculate the vertical
current inside a steel-cased. The model is based on the Schenkel and
Morrison Casing Model, and the results are used in a 2016 SEG abstract by
Yang et al.
- Schenkel, C.J., and H.F. Morrison, 1990, Effects of well casing on potential field measurements using downhole current sources: Geophysical prospecting, 38, 663-686.
The model consists of:
- Air: Conductivity 1e-8 S/m, above z = 0
- Background: conductivity 1e-2 S/m, below z = 0
- Casing: conductivity 1e6 S/m
- 300m long
- radius of 0.1m
- thickness of 6e-3m
Inside the casing, we take the same conductivity as the background.
We are using an EM code to simulate DC, so we use frequency low enough
that the skin depth inside the casing is longer than the casing length (f
= 1e-6 Hz). The plot produced is of the current inside the casing.
These results are shown in the SEG abstract by Yang et al., 2016: 3D DC
resistivity modeling of steel casing for reservoir monitoring using
equivalent resistor network. The solver used to produce these results and
achieve the CPU time of ~30s is Mumps, which was installed using pymatsolver_
.. _pymatsolver: https://github.com/rowanc1/pymatsolver
This example is on figshare: https://dx.doi.org/10.6084/m9.figshare.3126961.v1
If you would use this example for a code comparison, or build upon it, a
citation would be much appreciated!
.. plot::
from SimPEG import Examples
Examples.EM_Schenkel_Morrison_Casing.run()
.. literalinclude:: ../../SimPEG/Examples/EM_Schenkel_Morrison_Casing.py
:language: python
:linenos:
setup.py
+1 -1
View File
@@ -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',
tests/em/fdem/forward/test_FDEM_forward.py
+30 -80
View File
@@ -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()
tests/em/fdem/forward/test_FDEM_forwardEJHB.py
+125
View File
@@ -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()
tests/em/fdem/forward/test_FDEM_forwardHB.py
+128
View File
@@ -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()
tests/em/fdem/inverse/adjoint/test_FDEM_adjoint.py → tests/em/fdem/inverse/adjoint/test_FDEM_adjointEB.py
+57 -57
View File
@@ -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__':
tests/em/fdem/inverse/adjoint/test_FDEM_adjointHJ.py
+155
View File
@@ -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()
tests/em/fdem/inverse/derivs/test_FDEM_derivs.py
+116 -10
View File
@@ -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()
tests/flow/test_Richards.py
+6 -6
View File
@@ -116,8 +116,8 @@ class RichardsTests1D(unittest.TestCase):
v = np.random.rand(self.survey.nD)
z = np.random.rand(self.M.nC)
Hs = self.prob.fields(self.Ks)
vJz = v.dot(self.prob.Jvec(self.Ks,z,u=Hs))
zJv = z.dot(self.prob.Jtvec(self.Ks,v,u=Hs))
vJz = v.dot(self.prob.Jvec(self.Ks,z,f=Hs))
zJv = z.dot(self.prob.Jtvec(self.Ks,v,f=Hs))
tol = TOL*(10**int(np.log10(np.abs(zJv))))
passed = np.abs(vJz - zJv) < tol
print 'Richards Adjoint Test - PressureHead'
@@ -188,8 +188,8 @@ class RichardsTests2D(unittest.TestCase):
v = np.random.rand(self.survey.nD)
z = np.random.rand(self.M.nC)
Hs = self.prob.fields(self.Ks)
vJz = v.dot(self.prob.Jvec(self.Ks,z,u=Hs))
zJv = z.dot(self.prob.Jtvec(self.Ks,v,u=Hs))
vJz = v.dot(self.prob.Jvec(self.Ks,z,f=Hs))
zJv = z.dot(self.prob.Jtvec(self.Ks,v,f=Hs))
tol = TOL*(10**int(np.log10(np.abs(zJv))))
passed = np.abs(vJz - zJv) < tol
print '2D: Richards Adjoint Test - PressureHead'
@@ -260,8 +260,8 @@ class RichardsTests3D(unittest.TestCase):
v = np.random.rand(self.survey.nD)
z = np.random.rand(self.M.nC)
Hs = self.prob.fields(self.Ks)
vJz = v.dot(self.prob.Jvec(self.Ks,z,u=Hs))
zJv = z.dot(self.prob.Jtvec(self.Ks,v,u=Hs))
vJz = v.dot(self.prob.Jvec(self.Ks,z,f=Hs))
zJv = z.dot(self.prob.Jtvec(self.Ks,v,f=Hs))
tol = TOL*(10**int(np.log10(np.abs(zJv))))
passed = np.abs(vJz - zJv) < tol
print '3D: Richards Adjoint Test - PressureHead'