mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-02 00:46:36 +08:00
Lindsey
+36
-1
@@ -17,6 +17,14 @@ class BaseEMProblem(Problem.BaseProblem):
|
||||
|
||||
verbose = False
|
||||
|
||||
####################################################
|
||||
# Make A Symmetric
|
||||
####################################################
|
||||
@property
|
||||
def _makeASymmetric(self):
|
||||
if getattr(self, '__makeASymmetric', None) is None:
|
||||
self.__makeASymmetric = True
|
||||
return self.__makeASymmetric
|
||||
|
||||
####################################################
|
||||
# Mu Model
|
||||
@@ -30,6 +38,10 @@ class BaseEMProblem(Problem.BaseProblem):
|
||||
def mu(self, value):
|
||||
if getattr(self, '_MfMui', None) is not None:
|
||||
del self._MfMui
|
||||
if getattr(self, '_MeMu', None) is not None:
|
||||
del delf._MeMu
|
||||
if getattr(self, '_MeMuI', None) is not None:
|
||||
del self._MeMuI
|
||||
self._mu = value
|
||||
|
||||
|
||||
@@ -44,6 +56,20 @@ class BaseEMProblem(Problem.BaseProblem):
|
||||
self._MfMui = self.mesh.getFaceInnerProduct(1/self.mu)
|
||||
return self._MfMui
|
||||
|
||||
@property
|
||||
def MeMuI(self):
|
||||
# TODO: Assuming isotropic mu
|
||||
if getattr(self, '_MeMuI', None) is None:
|
||||
self._MeMuI = self.mesh.getEdgeInnerProduct(self.mu, invMat=True)
|
||||
return self._MeMuI
|
||||
|
||||
@property
|
||||
def MeMu(self):
|
||||
#TODO: Assuming isotropic mu
|
||||
if getattr(self, '_MeMu', None) is None:
|
||||
self._MeMu = self.mesh.getEdgeInnerProduct(self.mu)
|
||||
return self._MeMu
|
||||
|
||||
@property
|
||||
def Me(self):
|
||||
if getattr(self, '_Me', None) is None:
|
||||
@@ -66,7 +92,16 @@ class BaseEMProblem(Problem.BaseProblem):
|
||||
self._MeSigmaI = self.mesh.getEdgeInnerProduct(sigma, invMat=True)
|
||||
return self._MeSigmaI
|
||||
|
||||
deleteTheseOnModelUpdate = ['_MeSigma', '_MeSigmaI']
|
||||
@property
|
||||
def MfSigmai(self):
|
||||
#TODO: hardcoded to sigma as the model
|
||||
#TODO: hardcoded to sigma diagonal
|
||||
if getattr(self, '_MfSigmai', None) is None:
|
||||
sigma = self.curModel.transform
|
||||
self._MfSigmai = self.mesh.getFaceInnerProduct(1/sigma)
|
||||
return self._MfSigmai
|
||||
|
||||
deleteTheseOnModelUpdate = ['_MeSigma', '_MeSigmaI','_MfSigmai']
|
||||
|
||||
def fields(self, m):
|
||||
self.curModel = m
|
||||
|
||||
+345
-77
@@ -8,6 +8,102 @@ def omega(freq):
|
||||
"""Change frequency to angular frequency, omega"""
|
||||
return 2.*np.pi*freq
|
||||
|
||||
def getSource(self,freq):
|
||||
"""
|
||||
:param float freq: Frequency
|
||||
:rtype: numpy.ndarray (nE, nTx)
|
||||
:return: RHS
|
||||
"""
|
||||
Txs = self.survey.getTransmitters(freq)
|
||||
rhs = range(len(Txs))
|
||||
|
||||
solType = self.solType
|
||||
|
||||
if solType == 'e' or solType == 'b':
|
||||
gridEJx = self.mesh.gridEx
|
||||
gridEJy = self.mesh.gridEy
|
||||
gridEJz = self.mesh.gridEz
|
||||
nEJ = self.mesh.nE
|
||||
|
||||
gridBHx = self.mesh.gridFx
|
||||
gridBHy = self.mesh.gridFy
|
||||
gridBHz = self.mesh.gridFz
|
||||
nBH = self.mesh.nF
|
||||
|
||||
|
||||
C = self.mesh.edgeCurl
|
||||
mui = self.MfMui
|
||||
|
||||
elif solType == 'h' or solType == 'j':
|
||||
gridEJx = self.mesh.gridFx
|
||||
gridEJy = self.mesh.gridFy
|
||||
gridEJz = self.mesh.gridFz
|
||||
nEJ = self.mesh.nF
|
||||
|
||||
gridBHx = self.mesh.gridEx
|
||||
gridBHy = self.mesh.gridEy
|
||||
gridBHz = self.mesh.gridEz
|
||||
nBH = self.mesh.nE
|
||||
|
||||
C = self.mesh.edgeCurl.T
|
||||
mui = self.MeMuI
|
||||
|
||||
else:
|
||||
NotImplementedError('Only E or F sources')
|
||||
|
||||
for i, tx in enumerate(Txs):
|
||||
if self.mesh._meshType is 'CYL':
|
||||
if self.mesh.isSymmetric:
|
||||
if tx.txType == 'VMD':
|
||||
SRC = Sources.MagneticDipoleVectorPotential(tx.loc, gridEJy, 'y')
|
||||
elif tx.txType =='CircularLoop':
|
||||
SRC = Sources.MagneticLoopVectorPotential(tx.loc, gridEJy, 'y', tx.radius)
|
||||
else:
|
||||
raise NotImplementedError('Only VMD and CircularLoop')
|
||||
else:
|
||||
raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
|
||||
|
||||
elif self.mesh._meshType is 'TENSOR':
|
||||
|
||||
if tx.txType == 'VMD':
|
||||
src = Sources.MagneticDipoleVectorPotential
|
||||
SRCx = src(tx.loc, gridEJx, 'x')
|
||||
SRCy = src(tx.loc, gridEJy, 'y')
|
||||
SRCz = src(tx.loc, gridEJz, 'z')
|
||||
|
||||
elif tx.txType == 'VMD_B':
|
||||
src = Sources.MagneticDipoleFields
|
||||
SRCx = src(tx.loc, gridBHx, 'x')
|
||||
SRCy = src(tx.loc, gridBHy, 'y')
|
||||
SRCz = src(tx.loc, gridBHz, 'z')
|
||||
|
||||
elif tx.txType == 'CircularLoop':
|
||||
src = Sources.MagneticLoopVectorPotential
|
||||
SRCx = src(tx.loc, gridEJx, 'x', tx.radius)
|
||||
SRCy = src(tx.loc, gridEJy, 'y', tx.radius)
|
||||
SRCz = src(tx.loc, gridEJz, 'z', tx.radius)
|
||||
else:
|
||||
|
||||
raise NotImplemented('%s txType is not implemented' % tx.txType)
|
||||
SRC = np.concatenate((SRCx, SRCy, SRCz))
|
||||
|
||||
else:
|
||||
raise Exception('Unknown mesh for VMD')
|
||||
|
||||
rhs[i] = SRC
|
||||
|
||||
# b-forumlation
|
||||
if tx.txType == 'VMD_B':
|
||||
b_0 = np.concatenate(rhs).reshape((nBH, len(Txs)), order='E')
|
||||
else:
|
||||
a = np.concatenate(rhs).reshape((nEJ, len(Txs)), order='F')
|
||||
b_0 = C*a
|
||||
|
||||
if solType == 'b' or solType == 'h':
|
||||
return b_0
|
||||
elif solType == 'e' or solType == 'j':
|
||||
return C.T*mui*b_0
|
||||
|
||||
class BaseFDEMProblem(BaseEMProblem):
|
||||
"""
|
||||
We start by looking at Maxwell's equations in the electric field \\(\\vec{E}\\) and the magnetic flux density \\(\\vec{B}\\):
|
||||
@@ -18,7 +114,6 @@ class BaseFDEMProblem(BaseEMProblem):
|
||||
\\nabla \\times \\mu^{-1} \\vec{B} - \\sigma \\vec{E} = \\vec{J_s}
|
||||
|
||||
"""
|
||||
|
||||
surveyPair = SurveyFDEM
|
||||
|
||||
def forward(self, m, RHS, CalcFields):
|
||||
@@ -105,6 +200,12 @@ class BaseFDEMProblem(BaseEMProblem):
|
||||
|
||||
return Jtv
|
||||
|
||||
def getSource(self,freq):
|
||||
return self.getSource(freq)
|
||||
|
||||
##########################################################################################
|
||||
################################ E-B Formulation #########################################
|
||||
##########################################################################################
|
||||
|
||||
class ProblemFDEM_e(BaseFDEMProblem):
|
||||
"""
|
||||
@@ -141,6 +242,7 @@ class ProblemFDEM_e(BaseFDEMProblem):
|
||||
|
||||
return C.T*mui*C + 1j*omega(freq)*sig
|
||||
|
||||
|
||||
def getADeriv(self, freq, u, v, adjoint=False):
|
||||
sig = self.curModel.transform
|
||||
dsig_dm = self.curModel.transformDeriv
|
||||
@@ -157,29 +259,8 @@ class ProblemFDEM_e(BaseFDEMProblem):
|
||||
:rtype: numpy.ndarray (nE, nTx)
|
||||
:return: RHS
|
||||
"""
|
||||
Txs = self.survey.getTransmitters(freq)
|
||||
rhs = range(len(Txs))
|
||||
for i, tx in enumerate(Txs):
|
||||
if tx.txType == 'VMD':
|
||||
src = Sources.MagneticDipoleVectorPotential
|
||||
SRCx = src(tx.loc, self.mesh.gridEx, 'x')
|
||||
SRCy = src(tx.loc, self.mesh.gridEy, 'y')
|
||||
SRCz = src(tx.loc, self.mesh.gridEz, 'z')
|
||||
|
||||
elif tx.txType == 'CircularLoop':
|
||||
src = Sources.MagneticLoopVectorPotential
|
||||
SRCx = src(tx.loc, self.mesh.gridEx, 'x', tx.radius)
|
||||
SRCy = src(tx.loc, self.mesh.gridEy, 'y', tx.radius)
|
||||
SRCz = src(tx.loc, self.mesh.gridEz, 'z', tx.radius)
|
||||
else:
|
||||
raise NotImplemented('%s txType is not implemented' % tx.txType)
|
||||
rhs[i] = np.concatenate((SRCx, SRCy, SRCz))
|
||||
|
||||
a = np.concatenate(rhs).reshape((self.mesh.nE, len(Txs)), order='F')
|
||||
mui = self.MfMui
|
||||
C = self.mesh.edgeCurl
|
||||
|
||||
j_s = C.T*mui*C*a
|
||||
j_s = getSource(self,freq)
|
||||
return -1j*omega(freq)*j_s
|
||||
|
||||
def calcFields(self, sol, freq, fieldType, adjoint=False):
|
||||
@@ -221,8 +302,13 @@ class ProblemFDEM_b(BaseFDEMProblem):
|
||||
mui = self.MfMui
|
||||
sigI = self.MeSigmaI
|
||||
C = self.mesh.edgeCurl
|
||||
iomega = 1j * omega(freq) * sp.eye(self.mesh.nF)
|
||||
|
||||
return mui*C*sigI*C.T*mui + 1j*omega(freq)*mui
|
||||
A = C*sigI*C.T*mui + iomega
|
||||
|
||||
if self._makeASymmetric is True:
|
||||
return mui.T*A
|
||||
return A
|
||||
|
||||
def getADeriv(self, freq, u, v, adjoint=False):
|
||||
|
||||
@@ -237,9 +323,14 @@ class ProblemFDEM_b(BaseFDEMProblem):
|
||||
dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig)(vec)
|
||||
|
||||
if adjoint:
|
||||
return dsig_dm.T * ( dMe_dsig.T * ( dMeSigmaI_dI.T * ( C.T * ( mui.T * v ) ) ) )
|
||||
if self._makeASymmetric is True:
|
||||
v = mui * v
|
||||
return dsig_dm.T * ( dMe_dsig.T * ( dMeSigmaI_dI.T * ( C.T * v ) ) )
|
||||
|
||||
if self._makeASymmetric is True:
|
||||
return mui.T * ( C * ( dMeSigmaI_dI * ( dMe_dsig * ( dsig_dm * v ) ) ) )
|
||||
return C * ( dMeSigmaI_dI * ( dMe_dsig * ( dsig_dm * v ) ) )
|
||||
|
||||
return mui * ( C * ( dMeSigmaI_dI * ( dMe_dsig * ( dsig_dm * v ) ) ) )
|
||||
|
||||
def getRHS(self, freq):
|
||||
"""
|
||||
@@ -247,59 +338,14 @@ class ProblemFDEM_b(BaseFDEMProblem):
|
||||
:rtype: numpy.ndarray (nE, nTx)
|
||||
:return: RHS
|
||||
"""
|
||||
Txs = self.survey.getTransmitters(freq)
|
||||
rhs = range(len(Txs))
|
||||
for i, tx in enumerate(Txs):
|
||||
|
||||
if self.mesh._meshType is 'CYL':
|
||||
if self.mesh.isSymmetric:
|
||||
if tx.txType == 'VMD':
|
||||
SRC = Sources.MagneticDipoleVectorPotential(tx.loc, self.mesh.gridEy, 'y')
|
||||
elif tx.txType =='CircularLoop':
|
||||
SRC = Sources.MagneticLoopVectorPotential(tx.loc, self.mesh.gridEy, 'y', tx.radius)
|
||||
else:
|
||||
raise NotImplementedError('Only VMD and CircularLoop')
|
||||
else:
|
||||
raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
|
||||
|
||||
elif self.mesh._meshType is 'TENSOR':
|
||||
b_0 = getSource(self,freq)
|
||||
|
||||
if tx.txType == 'VMD':
|
||||
src = Sources.MagneticDipoleVectorPotential
|
||||
SRCx = src(tx.loc, self.mesh.gridEx, 'x')
|
||||
SRCy = src(tx.loc, self.mesh.gridEy, 'y')
|
||||
SRCz = src(tx.loc, self.mesh.gridEz, 'z')
|
||||
|
||||
elif tx.txType == 'VMD_B':
|
||||
src = Sources.MagneticDipoleFields
|
||||
SRCx = src(tx.loc, self.mesh.gridFx, 'x')
|
||||
SRCy = src(tx.loc, self.mesh.gridFy, 'y')
|
||||
SRCz = src(tx.loc, self.mesh.gridFz, 'z')
|
||||
|
||||
elif tx.txType == 'CircularLoop':
|
||||
src = Sources.MagneticLoopVectorPotential
|
||||
SRCx = src(tx.loc, self.mesh.gridEx, 'x', tx.radius)
|
||||
SRCy = src(tx.loc, self.mesh.gridEy, 'y', tx.radius)
|
||||
SRCz = src(tx.loc, self.mesh.gridEz, 'z', tx.radius)
|
||||
else:
|
||||
|
||||
raise NotImplemented('%s txType is not implemented' % tx.txType)
|
||||
SRC = np.concatenate((SRCx, SRCy, SRCz))
|
||||
|
||||
else:
|
||||
raise Exception('Unknown mesh for VMD')
|
||||
|
||||
rhs[i] = SRC
|
||||
|
||||
mui = self.MfMui
|
||||
if tx.txType == 'VMD_B':
|
||||
b_0 = np.concatenate(rhs).reshape((self.mesh.nF, len(Txs)), order='F')
|
||||
else:
|
||||
a = np.concatenate(rhs).reshape((self.mesh.nE, len(Txs)), order='F')
|
||||
C = self.mesh.edgeCurl
|
||||
b_0 = C*a
|
||||
|
||||
return -1j*omega(freq)*mui*b_0
|
||||
rhs = -1j*omega(freq)*b_0
|
||||
if self._makeASymmetric is True:
|
||||
mui = self.MfMui
|
||||
return mui.T*rhs
|
||||
return rhs
|
||||
|
||||
def calcFields(self, sol, freq, fieldType, adjoint=False):
|
||||
b = sol
|
||||
@@ -335,3 +381,225 @@ class ProblemFDEM_b(BaseFDEMProblem):
|
||||
return None
|
||||
raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
|
||||
|
||||
|
||||
##########################################################################################
|
||||
################################ H-J Formulation #########################################
|
||||
##########################################################################################
|
||||
|
||||
|
||||
class ProblemFDEM_j(BaseFDEMProblem):
|
||||
"""
|
||||
Using the H-J formulation of Maxwell's equations
|
||||
|
||||
.. math::
|
||||
\\nabla \\times \\sigma^{-1} \\vec{J} + i\\omega\\mu\\vec{H} = 0
|
||||
\\nabla \\times \\vec{H} - \\vec{J} = \\vec{J_s}
|
||||
|
||||
Since \(\\vec{J}\) is a flux and \(\\vec{H}\) is a field, we discretize \(\\vec{J}\) on faces and \(\\vec{H}\) on edges.
|
||||
|
||||
For this implementation, we solve for J using \( \\vec{H} = - (i\\omega\\mu)^{-1} \\nabla \\times \\sigma^{-1} \\vec{J} \) :
|
||||
|
||||
.. math::
|
||||
\\nabla \\times ( \\mu^{-1} \\nabla \\times \\sigma^{-1} \\vec{J} ) + i\\omega \\vec{J} = - i\\omega\\vec{J_s}
|
||||
|
||||
We discretize this to:
|
||||
|
||||
.. math::
|
||||
(\\mathbf{C} \\mathbf{M^e_{mu^{-1}}} \\mathbf{C^T} \\mathbf{M^f_{\\sigma^{-1}}} + i\\omega ) \\mathbf{j} = - i\\omega \\mathbf{j_s}
|
||||
|
||||
.. note::
|
||||
This implementation does not yet work with full anisotropy!!
|
||||
|
||||
"""
|
||||
|
||||
solType = 'j'
|
||||
storeTheseFields = ['j','h']
|
||||
|
||||
def __init__(self, model, **kwargs):
|
||||
BaseFDEMProblem.__init__(self, model, **kwargs)
|
||||
|
||||
def getA(self, freq):
|
||||
"""
|
||||
Here, we form the operator \(\\mathbf{A}\) to solce
|
||||
.. math::
|
||||
\\mathbf{A} = \\mathbf{C} \\mathbf{M^e_{mu^{-1}}} \\mathbf{C^T} \\mathbf{M^f_{\\sigma^{-1}}} + i\\omega
|
||||
|
||||
:param float freq: Frequency
|
||||
:rtype: scipy.sparse.csr_matrix
|
||||
:return: A
|
||||
"""
|
||||
|
||||
MeMuI = self.MeMuI
|
||||
MfSigi = self.MfSigmai
|
||||
C = self.mesh.edgeCurl
|
||||
iomega = 1j * omega(freq) * sp.eye(self.mesh.nF)
|
||||
|
||||
A = C * MeMuI * C.T * MfSigi + iomega
|
||||
|
||||
if self._makeASymmetric is True:
|
||||
return MfSigi.T*A
|
||||
return A
|
||||
|
||||
|
||||
def getADeriv(self, freq, u, v, adjoint=False):
|
||||
"""
|
||||
In this case, we assume that electrical conductivity, \(\\sigma\) is the physical property of interest (i.e. \(\sigma\) = model.transform). Then we want
|
||||
.. math::
|
||||
\\frac{\mathbf{A(\\sigma)} \mathbf{v}}{d \\mathbf{m}} &= \\mathbf{C} \\mathbf{M^e_{mu^{-1}}} \\mathbf{C^T} \\frac{d \\mathbf{M^f_{\\sigma^{-1}}}}{d \\mathbf{m}}
|
||||
&= \\mathbf{C} \\mathbf{M^e_{mu}^{-1}} \\mathbf{C^T} \\frac{d \\mathbf{M^f_{\\sigma^{-1}}}}{d \\mathbf{\\sigma^{-1}}} \\frac{d \\mathbf{\\sigma^{-1}}}{d \\mathbf{\\sigma}} \\frac{d \\mathbf{\\sigma}}{d \\mathbf{m}}
|
||||
"""
|
||||
|
||||
MeMuI = self.MeMuI
|
||||
MfSigi = self.MfSigmai
|
||||
C = self.mesh.edgeCurl
|
||||
sig = self.curModel.transform
|
||||
sigi = 1/sig
|
||||
dsig_dm = self.curModel.transformDeriv
|
||||
dsigi_dsig = -Utils.sdiag(sigi)**2
|
||||
dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(u)
|
||||
|
||||
if adjoint:
|
||||
if self._makeASymmetric is True:
|
||||
v = MfSigi * v
|
||||
return dsig_dm.T * ( dsigi_dsig.T *( dMf_dsigi.T * ( C * ( MeMuI.T * ( C.T * v ) ) ) ) )
|
||||
|
||||
if self._makeASymmetric is True:
|
||||
return MfSigi.T * ( C * ( MeMuI * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) ) ) )
|
||||
return C * ( MeMuI * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) ) )
|
||||
|
||||
|
||||
def getRHS(self, freq):
|
||||
"""
|
||||
:param float freq: Frequency
|
||||
:rtype: numpy.ndarray (nE, nTx)
|
||||
:return: RHS
|
||||
"""
|
||||
j_s = getSource(self,freq)
|
||||
rhs = -1j*omega(freq)*j_s
|
||||
|
||||
if self._makeASymmetric is True:
|
||||
MfSigi = self.MfSigmai
|
||||
return MfSigi.T*rhs
|
||||
return rhs
|
||||
|
||||
def calcFields(self, sol, freq, fieldType, adjoint=False):
|
||||
j = sol
|
||||
if fieldType == 'j':
|
||||
return j
|
||||
elif fieldType == 'h':
|
||||
MeMuI = self.MeMuI
|
||||
C = self.mesh.edgeCurl
|
||||
MfSigi = self.MfSigmai
|
||||
if not adjoint:
|
||||
h = -(1./(1j*omega(freq))) * MeMuI * ( C.T * ( MfSigi * j ) )
|
||||
else:
|
||||
h = -(1./(1j*omega(freq))) * MfSigi.T * ( C * ( MeMuI.T * j ) )
|
||||
return h
|
||||
raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
|
||||
|
||||
def calcFieldsDeriv(self, sol, freq, fieldType, v, adjoint=False):
|
||||
j = sol
|
||||
if fieldType == 'j':
|
||||
return None
|
||||
elif fieldType == 'h':
|
||||
MeMuI = self.MeMuI
|
||||
C = self.mesh.edgeCurl
|
||||
sig = self.curModel.transform
|
||||
sigi = 1/sig
|
||||
dsig_dm = self.curModel.transformDeriv
|
||||
dsigi_dsig = -Utils.sdiag(sigi)**2
|
||||
dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(j)
|
||||
sigi = self.MfSigmai
|
||||
if not adjoint:
|
||||
return -(1./(1j*omega(freq))) * MeMuI * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) )
|
||||
else:
|
||||
return -(1./(1j*omega(freq))) * dsig_dm.T * ( dsigi_dsig.T * ( dMf_dsigi.T * ( C * ( MeMuI.T * v ) ) ) )
|
||||
raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
|
||||
|
||||
|
||||
|
||||
|
||||
# Solving for h! - using primary- secondary approach
|
||||
class ProblemFDEM_h(BaseFDEMProblem):
|
||||
"""
|
||||
Using the H-J formulation of Maxwell's equations
|
||||
|
||||
.. math::
|
||||
\\nabla \\times \\sigma^{-1} \\vec{J} + i\\omega\\mu\\vec{H} = 0
|
||||
\\nabla \\times \\vec{H} - \\vec{J} = \\vec{J_s}
|
||||
|
||||
Since \(\\vec{J}\) is a flux and \(\\vec{H}\) is a field, we discretize \(\\vec{J}\) on faces and \(\\vec{H}\) on edges.
|
||||
|
||||
For this implementation, we solve for J using \( \\vec{J} = \\nabla \\times \\vec{H} - \\vec{J_s} \)
|
||||
|
||||
.. math::
|
||||
\\nabla \\times \\sigma^{-1} \\nabla \\times \\vec{H} + i\\omega\\mu\\vec{H} = \\nabla \\times \\sigma^{-1} \\vec{J_s}
|
||||
|
||||
We discretize and solve
|
||||
|
||||
.. math::
|
||||
(\\mathbf{C^T} \\mathbf{M^f_{\\sigma^{-1}}} \\mathbf{C} + i\\omega \\mathbf{M_{\mu}} ) \\mathbf{h} = \\mathbf{C^T} \\mathbf{M^f_{\\sigma^{-1}}} \\vec{J_s}
|
||||
|
||||
.. note::
|
||||
This implementation does not yet work with full anisotropy!!
|
||||
|
||||
"""
|
||||
|
||||
solType = 'h'
|
||||
storeTheseFields = ['j','h']
|
||||
|
||||
def __init__(self, model, **kwargs):
|
||||
BaseFDEMProblem.__init__(self, model, **kwargs)
|
||||
|
||||
def getA(self, freq):
|
||||
"""
|
||||
:param float freq: Frequency
|
||||
:rtype: scipy.sparse.csr_matrix
|
||||
:return: A
|
||||
"""
|
||||
|
||||
MeMu = self.MeMu
|
||||
MfSigi = self.MfSigmai
|
||||
C = self.mesh.edgeCurl
|
||||
|
||||
return C.T * MfSigi * C + 1j*omega(freq)*MeMu
|
||||
|
||||
def getADeriv(self, freq, u, v, adjoint=False):
|
||||
|
||||
MeMu = self.MeMu
|
||||
C = self.mesh.edgeCurl
|
||||
sig = self.curModel.transform
|
||||
sigi = 1/sig
|
||||
dsig_dm = self.curModel.transformDeriv
|
||||
dsigi_dsig = -Utils.sdiag(sigi)**2
|
||||
|
||||
dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(C*u)
|
||||
|
||||
if adjoint:
|
||||
return (dsig_dm.T * (dsigi_dsig.T * (dMf_dsigi.T * (C * v))))
|
||||
return (C.T * (dMf_dsigi * (dsigi_dsig * (dsig_dm * v))))
|
||||
|
||||
|
||||
|
||||
def getRHS(self, freq):
|
||||
"""
|
||||
:param float freq: Frequency
|
||||
:rtype: numpy.ndarray (nE, nTx)
|
||||
:return: RHS
|
||||
"""
|
||||
b_0 = getSource(self,freq)
|
||||
return -1j*omega(freq)*b_0
|
||||
|
||||
def calcFields(self, sol, freq, fieldType, adjoint=False):
|
||||
h = sol
|
||||
if fieldType == 'j':
|
||||
C = self.mesh.edgeCurl
|
||||
if adjoint:
|
||||
return C.T*h
|
||||
return C*h
|
||||
elif fieldType == 'h':
|
||||
return h
|
||||
raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
|
||||
|
||||
def calcFieldsDeriv(self, sol, freq, fieldType, v, adjoint=False):
|
||||
return None
|
||||
@@ -16,6 +16,20 @@ class RxFDEM(Survey.BaseRx):
|
||||
'bxi':['b', 'Fx', 'imag'],
|
||||
'byi':['b', 'Fy', 'imag'],
|
||||
'bzi':['b', 'Fz', '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'],
|
||||
|
||||
'hxr':['h', 'Ex', 'real'],
|
||||
'hyr':['h', 'Ey', 'real'],
|
||||
'hzr':['h', 'Ez', 'real'],
|
||||
'hxi':['h', 'Ex', 'imag'],
|
||||
'hyi':['h', 'Ey', 'imag'],
|
||||
'hzi':['h', 'Ez', 'imag'],
|
||||
}
|
||||
radius = None
|
||||
|
||||
@@ -84,7 +98,7 @@ class TxFDEM(Survey.BaseTx):
|
||||
|
||||
class FieldsFDEM(Problem.Fields):
|
||||
"""Fancy Field Storage for a FDEM survey."""
|
||||
knownFields = {'b': 'F', 'e': 'E'}
|
||||
knownFields = {'b': 'F', 'e': 'E', 'j': 'F', 'h': 'E'}
|
||||
dtype = complex
|
||||
|
||||
|
||||
|
||||
@@ -1,2 +1,2 @@
|
||||
from SurveyFDEM import *
|
||||
from FDEM import ProblemFDEM_e, ProblemFDEM_b
|
||||
from FDEM import ProblemFDEM_e, ProblemFDEM_b, ProblemFDEM_j, ProblemFDEM_h
|
||||
|
||||
+227
-13
@@ -5,9 +5,11 @@ import sys
|
||||
from scipy.constants import mu_0
|
||||
|
||||
TOL = 1e-4
|
||||
CONDUCTIVITY = 1e3
|
||||
FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order
|
||||
CONDUCTIVITY = 1e1
|
||||
MU = mu_0
|
||||
addrandoms = True
|
||||
freq = 1e-1
|
||||
addrandoms = True
|
||||
|
||||
def getProblem(fdemType, comp):
|
||||
cs = 5.
|
||||
@@ -20,17 +22,22 @@ def getProblem(fdemType, comp):
|
||||
|
||||
mapping = Maps.ExpMap(mesh)
|
||||
|
||||
x = np.linspace(-30,30,6)
|
||||
XYZ = Utils.ndgrid(x,x,np.r_[0])
|
||||
x = np.array([np.linspace(-30,-15,3),np.linspace(15,30,3)]) #don't sample right by the transmitter
|
||||
XYZ = Utils.ndgrid(x,x,np.r_[0.])
|
||||
Rx0 = EM.FDEM.RxFDEM(XYZ, comp)
|
||||
Tx0 = EM.FDEM.TxFDEM(np.r_[4.,2.,2.], 'VMD', 1e-2, [Rx0])
|
||||
Tx0 = EM.FDEM.TxFDEM(np.r_[0.,0.,0.], 'VMD', freq, [Rx0])
|
||||
|
||||
survey = EM.FDEM.SurveyFDEM([Tx0])
|
||||
|
||||
|
||||
if fdemType == 'e':
|
||||
prb = EM.FDEM.ProblemFDEM_e(mesh, mapping=mapping)
|
||||
elif fdemType == 'b':
|
||||
prb = EM.FDEM.ProblemFDEM_b(mesh, mapping=mapping)
|
||||
elif fdemType == 'j':
|
||||
prb = EM.FDEM.ProblemFDEM_j(mesh, mapping=mapping)
|
||||
elif fdemType == 'h':
|
||||
prb = EM.FDEM.ProblemFDEM_h(mesh, mapping=mapping)
|
||||
else:
|
||||
raise NotImplementedError()
|
||||
prb.pair(survey)
|
||||
@@ -51,19 +58,20 @@ def adjointTest(fdemType, comp):
|
||||
mu = np.log(np.ones(prb.mesh.nC)*MU)
|
||||
|
||||
if addrandoms is True:
|
||||
m = m + np.random.randn(prb.mesh.nC)*CONDUCTIVITY*1e-3
|
||||
mu = mu + np.random.randn(prb.mesh.nC)*MU*1e-3
|
||||
m = m + np.random.randn(prb.mesh.nC)*CONDUCTIVITY*1e-1
|
||||
mu = mu + np.random.randn(prb.mesh.nC)*MU*1e-1
|
||||
|
||||
prb.mu = mu
|
||||
survey = prb.survey
|
||||
|
||||
u = prb.fields(m)
|
||||
|
||||
v = np.random.rand(survey.nD)
|
||||
w = np.random.rand(prb.mapping.nP)
|
||||
|
||||
u = prb.fields(m)
|
||||
vJw = v.dot(prb.Jvec(m, w, u=u))
|
||||
wJtv = w.dot(prb.Jtvec(m, v, u=u))
|
||||
tol = TOL*(10**int(np.log10(np.abs(vJw))))
|
||||
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
|
||||
|
||||
@@ -74,14 +82,64 @@ def derivTest(fdemType, comp):
|
||||
mu = np.log(np.ones(prb.mesh.nC)*MU)
|
||||
|
||||
if addrandoms is True:
|
||||
x0 = x0 + np.random.randn(prb.mesh.nC)*CONDUCTIVITY*1e-3
|
||||
mu = mu + np.random.randn(prb.mesh.nC)*MU*1e-3
|
||||
x0 = x0 + np.random.randn(prb.mesh.nC)*CONDUCTIVITY*1e-1
|
||||
mu = mu + np.random.randn(prb.mesh.nC)*MU*1e-1
|
||||
|
||||
prb.mu = mu
|
||||
prb.mu = mu
|
||||
survey = prb.survey
|
||||
def fun(x):
|
||||
return survey.dpred(x), lambda x: prb.Jvec(x0, x)
|
||||
return Tests.checkDerivative(fun, x0, num=3, plotIt=False, eps=1e-25)
|
||||
return Tests.checkDerivative(fun, x0, num=3, plotIt=False, eps=FLR)
|
||||
|
||||
|
||||
def crossCheckTest(fdemType, comp):
|
||||
|
||||
l2norm = lambda r: np.sqrt(r.dot(r))
|
||||
|
||||
prb1 = getProblem(fdemType, comp)
|
||||
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)*CONDUCTIVITY*1e-1
|
||||
mu = mu + np.random.randn(mesh.nC)*MU*1e-1
|
||||
|
||||
prb1.mu = mu
|
||||
survey1 = prb1.survey
|
||||
|
||||
u1 = prb1.fields(m)
|
||||
d1 = Utils.mkvc(survey1.projectFields(u1))
|
||||
|
||||
prb1.unpair
|
||||
|
||||
if fdemType == 'e':
|
||||
prb2 = getProblem('b', comp)
|
||||
elif fdemType == 'b':
|
||||
prb2 = getProblem('e', comp)
|
||||
elif fdemType == 'j':
|
||||
prb2 = getProblem('h', comp)
|
||||
elif fdemType == 'h':
|
||||
prb2 = getProblem('j', comp)
|
||||
else:
|
||||
raise NotImplementedError()
|
||||
|
||||
prb2.mu = mu
|
||||
survey2 = prb2.survey
|
||||
|
||||
u2 = prb2.fields(m)
|
||||
d2 = Utils.mkvc(survey2.projectFields(u2))
|
||||
|
||||
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_DerivTests(unittest.TestCase):
|
||||
@@ -189,6 +247,162 @@ class FDEM_DerivTests(unittest.TestCase):
|
||||
self.assertTrue(adjointTest('b', 'bzi'))
|
||||
|
||||
|
||||
def test_Jvec_jxr_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'jxr'))
|
||||
def test_Jvec_jyr_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'jyr'))
|
||||
def test_Jvec_jzr_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'jzr'))
|
||||
def test_Jvec_jxi_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'jxi'))
|
||||
def test_Jvec_jyi_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'jyi'))
|
||||
def test_Jvec_jzi_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'jzi'))
|
||||
|
||||
def test_Jvec_hxr_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'hxr'))
|
||||
def test_Jvec_hyr_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'hyr'))
|
||||
def test_Jvec_hzr_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'hzr'))
|
||||
def test_Jvec_hxi_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'hxi'))
|
||||
def test_Jvec_hyi_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'hyi'))
|
||||
def test_Jvec_hzi_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'hzi'))
|
||||
|
||||
def test_Jvec_hxr_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'hxr'))
|
||||
def test_Jvec_hyr_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'hyr'))
|
||||
def test_Jvec_hzr_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'hzr'))
|
||||
def test_Jvec_hxi_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'hxi'))
|
||||
def test_Jvec_hyi_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'hyi'))
|
||||
def test_Jvec_hzi_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'hzi'))
|
||||
|
||||
def test_Jvec_hxr_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'jxr'))
|
||||
def test_Jvec_hyr_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'jyr'))
|
||||
def test_Jvec_hzr_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'jzr'))
|
||||
def test_Jvec_hxi_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'jxi'))
|
||||
def test_Jvec_hyi_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'jyi'))
|
||||
def test_Jvec_hzi_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'jzi'))
|
||||
|
||||
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_EB_CrossCheck_exr_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'exr'))
|
||||
def test_EB_CrossCheck_eyr_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'eyr'))
|
||||
def test_EB_CrossCheck_ezr_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'ezr'))
|
||||
def test_EB_CrossCheck_exi_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'exi'))
|
||||
def test_EB_CrossCheck_eyi_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'eyi'))
|
||||
def test_EB_CrossCheck_ezi_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'ezi'))
|
||||
|
||||
def test_EB_CrossCheck_bxr_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'bxr'))
|
||||
def test_EB_CrossCheck_byr_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'byr'))
|
||||
# def test_EB_CrossCheck_bzr_Eform(self):
|
||||
# self.assertTrue(crossCheckTest('e', 'bzr')) # Doesn't make sense to test this for p-s approach
|
||||
def test_EB_CrossCheck_bxi_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'bxi'))
|
||||
def test_EB_CrossCheck_byi_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'byi'))
|
||||
def test_EB_CrossCheck_bzi_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'bzi'))
|
||||
|
||||
def test_HJ_CrossCheck_jxr_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'jxr'))
|
||||
def test_HJ_CrossCheck_jyr_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'jyr'))
|
||||
def test_HJ_CrossCheck_jzr_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'jzr'))
|
||||
def test_HJ_CrossCheck_jxi_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'jxi'))
|
||||
def test_HJ_CrossCheck_jyi_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'jyi'))
|
||||
def test_HJ_CrossCheck_jzi_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'jzi'))
|
||||
|
||||
def test_HJ_CrossCheck_hxr_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'hxr'))
|
||||
def test_HJ_CrossCheck_hyr_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'hyr'))
|
||||
# def test_HJ_CrossCheck_hzr_Jform(self):
|
||||
# self.assertTrue(crossCheckTest('j', 'hzr')) # Doesn't make sense to test this for p-s approach
|
||||
def test_HJ_CrossCheck_hxi_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'hxi'))
|
||||
def test_HJ_CrossCheck_hyi_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'hyi'))
|
||||
def test_HJ_CrossCheck_hzi_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'hzi'))
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
|
||||
Reference in New Issue
Block a user