mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-19 11:28:00 +08:00
H-J Fields objects created and consistent. Derivatives not implemented yet
This commit is contained in:
Lindsey
+11
-116
@@ -1,7 +1,7 @@
|
||||
from SimPEG import Survey, Problem, Utils, np, sp, Solver as SimpegSolver
|
||||
from scipy.constants import mu_0
|
||||
from SurveyFDEM import SurveyFDEM
|
||||
from FieldsFDEM import FieldsFDEM, FieldsFDEM_e, FieldsFDEM_b
|
||||
from FieldsFDEM import FieldsFDEM, FieldsFDEM_e, FieldsFDEM_b, FieldsFDEM_h, FieldsFDEM_j
|
||||
from simpegEM.Base import BaseEMProblem
|
||||
from simpegEM.Utils.EMUtils import omega
|
||||
|
||||
@@ -12,8 +12,8 @@ class BaseFDEMProblem(BaseEMProblem):
|
||||
|
||||
.. math::
|
||||
|
||||
\\nabla \\times \\vec{E} + i \\omega \\vec{B} = 0 \\\\
|
||||
\\nabla \\times \\mu^{-1} \\vec{B} - \\sigma \\vec{E} = \\vec{J_s}
|
||||
\\nabla \\times \\vec{E} + i \\omega \\vec{B} = \\vec{S_m} \\\\
|
||||
\\nabla \\times \\mu^{-1} \\vec{B} - \\sigma \\vec{E} = \\vec{S_e}
|
||||
|
||||
"""
|
||||
surveyPair = SurveyFDEM
|
||||
@@ -28,10 +28,6 @@ class BaseFDEMProblem(BaseEMProblem):
|
||||
rhs = RHS(freq)
|
||||
Ainv = self.Solver(A, **self.solverOpts)
|
||||
sol = Ainv * rhs
|
||||
# for fieldType in self.storeTheseFields:
|
||||
# Txs = self.survey.getTransmitters(freq)
|
||||
# F[Txs, fieldType] = CalcFields(sol, freq, fieldType)
|
||||
|
||||
Txs = self.survey.getTransmitters(freq)
|
||||
F[Txs, self._fieldType] = sol
|
||||
|
||||
@@ -149,6 +145,7 @@ class ProblemFDEM_e(BaseFDEMProblem):
|
||||
"""
|
||||
|
||||
_fieldType = 'e'
|
||||
_eqLocs = 'FE'
|
||||
fieldsPair = FieldsFDEM_e
|
||||
|
||||
def __init__(self, model, **kwargs):
|
||||
@@ -202,33 +199,12 @@ class ProblemFDEM_e(BaseFDEMProblem):
|
||||
return RHS
|
||||
|
||||
|
||||
# def calcFields(self, sol, freq, fieldType, adjoint=False):
|
||||
# e = sol
|
||||
# if fieldType == 'e':
|
||||
# return e
|
||||
# elif fieldType == 'b':
|
||||
# if not adjoint:
|
||||
# b = - self.mesh.edgeCurl * e
|
||||
# b = 1./(1j*omega(freq)) * b
|
||||
# else:
|
||||
# b = -(1./(1j*omega(freq))) * ( self.mesh.edgeCurl.T * e )
|
||||
# return b
|
||||
# raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
|
||||
|
||||
# def calcFieldsDeriv(self, sol, freq, fieldType, v, adjoint=False):
|
||||
# e = sol
|
||||
# if fieldType == 'e':
|
||||
# return None
|
||||
# elif fieldType == 'b':
|
||||
# return None
|
||||
# raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
|
||||
|
||||
|
||||
class ProblemFDEM_b(BaseFDEMProblem):
|
||||
"""
|
||||
Solving for b!
|
||||
"""
|
||||
_fieldType = 'b'
|
||||
_eqLocs = 'FE'
|
||||
fieldsPair = FieldsFDEM_b
|
||||
|
||||
def __init__(self, model, **kwargs):
|
||||
@@ -302,41 +278,6 @@ class ProblemFDEM_b(BaseFDEMProblem):
|
||||
return RHS
|
||||
|
||||
|
||||
# def calcFields(self, sol, freq, fieldType, adjoint=False):
|
||||
# b = sol
|
||||
# if fieldType == 'e':
|
||||
# if not adjoint:
|
||||
# e = self.MeSigmaI * ( self.mesh.edgeCurl.T * ( self.MfMui * b ) )
|
||||
# else:
|
||||
# e = self.MfMui.T * ( self.mesh.edgeCurl * ( self.MeSigmaI.T * b ) )
|
||||
# return e
|
||||
# elif fieldType == 'b':
|
||||
# return b
|
||||
# raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
|
||||
|
||||
|
||||
# def calcFieldsDeriv(self, sol, freq, fieldType, v, adjoint=False):
|
||||
# b = sol
|
||||
# if fieldType == 'e':
|
||||
# sig = self.curModel.transform
|
||||
# dsig_dm = self.curModel.transformDeriv
|
||||
|
||||
# C = self.mesh.edgeCurl
|
||||
# mui = self.MfMui
|
||||
|
||||
# #TODO: This only works if diagonal (no tensors)...
|
||||
# dMeSigmaI_dI = - self.MeSigmaI**2
|
||||
|
||||
# vec = C.T * ( mui * b )
|
||||
# dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig)(vec)
|
||||
# if not adjoint:
|
||||
# return dMeSigmaI_dI * ( dMe_dsig * ( dsig_dm * v ) )
|
||||
# else:
|
||||
# return dsig_dm.T * ( dMe_dsig.T * ( dMeSigmaI_dI.T * v ) )
|
||||
# elif fieldType == 'b':
|
||||
# return None
|
||||
# raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
|
||||
|
||||
|
||||
##########################################################################################
|
||||
################################ H-J Formulation #########################################
|
||||
@@ -368,8 +309,9 @@ class ProblemFDEM_j(BaseFDEMProblem):
|
||||
|
||||
"""
|
||||
|
||||
solType = 'j'
|
||||
storeTheseFields = ['j','h']
|
||||
_fieldType = 'j'
|
||||
_eqLocs = 'EF'
|
||||
fieldsPair = FieldsFDEM_j
|
||||
|
||||
def __init__(self, model, **kwargs):
|
||||
BaseFDEMProblem.__init__(self, model, **kwargs)
|
||||
@@ -453,42 +395,6 @@ class ProblemFDEM_j(BaseFDEMProblem):
|
||||
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):
|
||||
@@ -516,8 +422,9 @@ class ProblemFDEM_h(BaseFDEMProblem):
|
||||
|
||||
"""
|
||||
|
||||
solType = 'h'
|
||||
storeTheseFields = ['j','h']
|
||||
_fieldType = 'h'
|
||||
_eqLocs = 'EF'
|
||||
fieldsPair = FieldsFDEM_h
|
||||
|
||||
def __init__(self, model, **kwargs):
|
||||
BaseFDEMProblem.__init__(self, model, **kwargs)
|
||||
@@ -575,16 +482,4 @@ class ProblemFDEM_h(BaseFDEMProblem):
|
||||
return RHS
|
||||
|
||||
|
||||
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
|
||||
|
||||
+149
-40
@@ -7,31 +7,6 @@ class FieldsFDEM(Problem.Fields):
|
||||
knownFields = None
|
||||
dtype = complex
|
||||
|
||||
# def calcFields(self,sol,tx,fieldType):
|
||||
# if fieldType == 'e':
|
||||
# return self._e(sol,tx)
|
||||
# elif fieldType == 'e_sec':
|
||||
# return self._e_sec(sol,tx)
|
||||
# elif fieldType == 'b':
|
||||
# return self._b(sol,tx)
|
||||
# elif fieldType == 'b_sec':
|
||||
# return self._b_sec(sol,tx)
|
||||
# else:
|
||||
# raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
|
||||
|
||||
# def calcFieldsDeriv(self,sol,tx,fieldType,adjoint=False):
|
||||
# if fieldType == 'e':
|
||||
# return self._eDeriv(sol,tx,adjoint)
|
||||
# elif fieldType == 'e_sec':
|
||||
# return self._e_secDeriv(sol,tx,adjoint)
|
||||
# elif fieldType == 'b':
|
||||
# return self._bDeriv(sol,tx,adjoint)
|
||||
# elif fieldType == 'b_sec':
|
||||
# return self._b_secDeriv(sol,tx,adjoint)
|
||||
# else:
|
||||
# raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
|
||||
|
||||
|
||||
|
||||
class FieldsFDEM_e(FieldsFDEM):
|
||||
knownFields = {'e':'E'}
|
||||
@@ -48,12 +23,6 @@ class FieldsFDEM_e(FieldsFDEM):
|
||||
self.getSource = self.survey.prob.getSource
|
||||
self.getSourceDeriv = self.survey.prob.getSourceDeriv
|
||||
|
||||
# def _e(self, e, tx):
|
||||
# return e
|
||||
|
||||
# def _eDeriv(self, e, tx, adjoint=False):
|
||||
# return None
|
||||
|
||||
def _b_sec(self, e, tx): #adjoint=False
|
||||
return - 1./(1j*omega(tx.freq)) * (self.edgeCurl * e)
|
||||
|
||||
@@ -62,13 +31,12 @@ class FieldsFDEM_e(FieldsFDEM):
|
||||
|
||||
def _b(self, e, tx): #adjoint=False
|
||||
b = self._b_sec(e,tx)
|
||||
print b.shape
|
||||
j_m,_ = self.getSource(tx.freq)
|
||||
if j_m[0] is not None:
|
||||
b += 1./(1j*omega(tx.freq)) * np.array([j_m[0]]).T
|
||||
return b
|
||||
|
||||
def _bDreiv(self, e, tx, adjoint=False):
|
||||
def _bDeriv(self, e, tx, adjoint=False):
|
||||
j_mDeriv,_ = self.getSourceDeriv(tx.freq, adjoint)
|
||||
b_secDeriv = self._b_secDeriv(e,tx.freq,adjoint)
|
||||
if j_mDeriv is None & b_secDeriv is None:
|
||||
@@ -98,12 +66,6 @@ class FieldsFDEM_b(FieldsFDEM):
|
||||
self.getSource = self.survey.prob.getSource
|
||||
self.getSourceDeriv = self.survey.prob.getSourceDeriv
|
||||
|
||||
# def _b(self, b, tx):
|
||||
# return b
|
||||
|
||||
# def _bDeriv(self, b, tx, adjoint=False):
|
||||
# return None
|
||||
|
||||
def _e_sec(self, b, tx):
|
||||
return self.MeSigmaI * ( self.edgeCurl.T * ( self.MfMui * b) )
|
||||
|
||||
@@ -128,4 +90,151 @@ class FieldsFDEM_b(FieldsFDEM):
|
||||
elif j_gDeriv is None:
|
||||
return e_secDeriv
|
||||
else:
|
||||
return e_secDeriv - j_gDeriv
|
||||
return e_secDeriv - j_gDeriv
|
||||
|
||||
|
||||
class FieldsFDEM_j(FieldsFDEM):
|
||||
knownFields = {'j':'F'}
|
||||
aliasFields = {
|
||||
'h_sec' : ['j','E','_h_sec'],
|
||||
'h' : ['j','E','_h']
|
||||
}
|
||||
|
||||
def __init__(self,mesh,survey,**kwargs):
|
||||
FieldsFDEM.__init__(self,mesh,survey,**kwargs)
|
||||
|
||||
def startup(self):
|
||||
self.edgeCurl = self.survey.prob.mesh.edgeCurl
|
||||
self.MeMuI = self.survey.prob.MeMuI
|
||||
self.MfSigmai = self.survey.prob.MfSigmai
|
||||
self.getSource = self.survey.prob.getSource
|
||||
self.getSourceDeriv = self.survey.prob.getSourceDeriv
|
||||
|
||||
def _h_sec(self, j, tx): #adjoint=False
|
||||
return - 1./(1j*omega(tx.freq)) * self.MeMuI * (self.edgeCurl.T * (self.MfSigmai * j) )
|
||||
|
||||
def _h_secDeriv(self, j, tx, adjoint=False):
|
||||
# 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)
|
||||
|
||||
def _h(self, j, tx): #adjoint=False
|
||||
h = self._h_sec(j,tx)
|
||||
j_m,_ = self.getSource(tx.freq)
|
||||
if j_m[0] is not None:
|
||||
h += 1./(1j*omega(tx.freq)) * self.MeMuI * np.array([j_m[0]]).T
|
||||
return h
|
||||
|
||||
def _hDeriv(self, j, tx, adjoint=False):
|
||||
j_mDeriv,_ = self.getSourceDeriv(tx.freq, adjoint)
|
||||
h_secDeriv = self._h_secDeriv(j,tx.freq,adjoint)
|
||||
if j_mDeriv is None & h_secDeriv is None:
|
||||
return None
|
||||
elif h_secDeriv is None:
|
||||
return 1./(1j*omega(tx.freq)) * j_mDeriv
|
||||
elif j_mDeriv is None:
|
||||
return h_secDeriv
|
||||
else:
|
||||
return 1./(1j*omega(tx.freq)) * j_mDeriv + h_secDeriv
|
||||
|
||||
class FieldsFDEM_h(FieldsFDEM):
|
||||
knownFields = {'h':'E'}
|
||||
aliasFields = {
|
||||
'j_sec' : ['h','F','_j_sec'],
|
||||
'j' : ['h','F','_j']
|
||||
}
|
||||
|
||||
def __init__(self,mesh,survey,**kwargs):
|
||||
FieldsFDEM.__init__(self,mesh,survey,**kwargs)
|
||||
|
||||
def startup(self):
|
||||
self.edgeCurl = self.survey.prob.mesh.edgeCurl
|
||||
self.MeMuI = self.survey.prob.MeMuI
|
||||
self.MfSigmai = self.survey.prob.MfSigmai
|
||||
self.getSource = self.survey.prob.getSource
|
||||
self.getSourceDeriv = self.survey.prob.getSourceDeriv
|
||||
|
||||
def _j_sec(self, h, tx): #adjoint=False
|
||||
return self.edgeCurl*h
|
||||
|
||||
def _j_secDeriv(self, h, tx, adjoint=False):
|
||||
return None
|
||||
|
||||
def _j(self, h, tx): #adjoint=False
|
||||
j = self._j_sec(h,tx)
|
||||
_,j_g = self.getSource(tx.freq)
|
||||
if j_g[0] is not None:
|
||||
j += -np.array([j_g[0]]).T
|
||||
return j
|
||||
|
||||
def _jDeriv(self, h, tx, adjoint=False):
|
||||
_,j_gDeriv = self.getSourceDeriv(tx.freq, adjoint)
|
||||
j_secDeriv = self._j_secDeriv(j,tx.freq,adjoint)
|
||||
if j_gDeriv is None & j_secDeriv is None:
|
||||
return None
|
||||
elif j_secDeriv is None:
|
||||
return - j_gDeriv
|
||||
elif j_gDeriv is None:
|
||||
return j_secDeriv
|
||||
else:
|
||||
return - j_gDeriv + j_secDeriv
|
||||
|
||||
|
||||
# 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
|
||||
# MfSigmai = self.MfSigmai
|
||||
# if not adjoint:
|
||||
# h = -(1./(1j*omega(freq))) * MeMuI * ( C.T * ( MfSigmai * j ) )
|
||||
# else:
|
||||
# h = -(1./(1j*omega(freq))) * MfSigmai.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)
|
||||
|
||||
|
||||
# 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
|
||||
@@ -83,7 +83,7 @@ class RxFDEM(Survey.BaseRx):
|
||||
return Pv
|
||||
|
||||
# SrcFDEM
|
||||
class SrcFDEM(Survey.BaseTx):
|
||||
class TxFDEM(Survey.BaseTx):
|
||||
#TODO: Break these out into Classes of Sources.
|
||||
|
||||
freq = None #: Frequency (float)
|
||||
|
||||
@@ -9,7 +9,7 @@ testDerivs = False
|
||||
testCrossCheck = True
|
||||
testAdjoint = False
|
||||
testEB = True
|
||||
testHJ = False
|
||||
testHJ = True
|
||||
|
||||
verbose = False
|
||||
|
||||
|
||||
Reference in New Issue
Block a user