From e980477031e6e2f6f9edf43d9d4c953803496bf0 Mon Sep 17 00:00:00 2001 From: Lindsey Date: Thu, 16 Apr 2015 16:52:29 -0700 Subject: [PATCH] H-J Fields objects created and consistent. Derivatives not implemented yet --- simpegEM/FDEM/FDEM.py | 127 +++--------------------- simpegEM/FDEM/FieldsFDEM.py | 189 ++++++++++++++++++++++++++++-------- simpegEM/FDEM/SurveyFDEM.py | 2 +- simpegEM/Tests/test_FDEM.py | 2 +- 4 files changed, 162 insertions(+), 158 deletions(-) diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py index acd46d8e..fe89951f 100644 --- a/simpegEM/FDEM/FDEM.py +++ b/simpegEM/FDEM/FDEM.py @@ -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 diff --git a/simpegEM/FDEM/FieldsFDEM.py b/simpegEM/FDEM/FieldsFDEM.py index ad8b7bc3..0f8c8596 100644 --- a/simpegEM/FDEM/FieldsFDEM.py +++ b/simpegEM/FDEM/FieldsFDEM.py @@ -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 \ No newline at end of file + 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 \ No newline at end of file diff --git a/simpegEM/FDEM/SurveyFDEM.py b/simpegEM/FDEM/SurveyFDEM.py index 82556277..e7553c43 100644 --- a/simpegEM/FDEM/SurveyFDEM.py +++ b/simpegEM/FDEM/SurveyFDEM.py @@ -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) diff --git a/simpegEM/Tests/test_FDEM.py b/simpegEM/Tests/test_FDEM.py index 285d9b5e..dc72f02c 100644 --- a/simpegEM/Tests/test_FDEM.py +++ b/simpegEM/Tests/test_FDEM.py @@ -9,7 +9,7 @@ testDerivs = False testCrossCheck = True testAdjoint = False testEB = True -testHJ = False +testHJ = True verbose = False