From 7239d9cbe63b649f9d0fbd989cc5473e2d194cc2 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 16 Jun 2015 11:28:29 -0700 Subject: [PATCH] FDEM problems: e,b,h,j up and running with Jvec and Jtvec within the re-factored framework. Whenever a new element is created, we create its derive wrt u (the computed field) and wrt m (the model). Jvec and Jtvec then stitch these pieces together using chain rule --- simpegEM/Base.py | 3 +- simpegEM/FDEM/FDEM.py | 35 +++--------- simpegEM/FDEM/FieldsFDEM.py | 105 +++++++++--------------------------- simpegEM/Tests/test_FDEM.py | 10 ++-- 4 files changed, 39 insertions(+), 114 deletions(-) diff --git a/simpegEM/Base.py b/simpegEM/Base.py index acd3af6c..4205c518 100644 --- a/simpegEM/Base.py +++ b/simpegEM/Base.py @@ -130,7 +130,8 @@ class BaseEMProblem(Problem.BaseProblem): # TODO: This should take a vector def MfRhoDeriv(self,u): - return self.mesh.getFaceInnerProductDeriv(self.curModel.rho)(u) * self.curModel.rhoDeriv + return self.mesh.getFaceInnerProductDeriv(self.curModel.rho)(u) * -Utils.sdiag(self.curModel.rho**2) * self.curModel.sigmaDeriv + # self.curModel.rhoDeriv @property def MfRhoI(self): diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py index 846d8364..036a6af0 100644 --- a/simpegEM/FDEM/FDEM.py +++ b/simpegEM/FDEM/FDEM.py @@ -36,7 +36,7 @@ class BaseFDEMProblem(BaseEMProblem): def Jvec(self, m, v, f=None): if f is None: - f = self.fields(m) # rename to f? + f = self.fields(m) self.curModel = m @@ -238,8 +238,6 @@ class ProblemFDEM_e(BaseFDEMProblem): C = self.mesh.edgeCurl MfMui = self.MfMui S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint) - # # evalDeriv(MfMui.T* (C * v), adjoint) - # raise Exception('Not implemented') if adjoint: dRHS = MfMui * (C * v) @@ -303,19 +301,14 @@ class ProblemFDEM_b(BaseFDEMProblem): MeSigmaIDeriv = MeSigmaIDeriv(vec) - # dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig)(vec) - if adjoint: if self._makeASymmetric is True: v = MfMui * v return MeSigmaIDeriv.T * (C.T * v) - # dsig_dm.T * ( dMe_dsig.T * ( dMeSigmaI_dI.T * ( C.T * v ) ) ) if self._makeASymmetric is True: - # return MfMui.T * ( C * ( dMeSigmaI_dI * ( dMe_dsig * ( dsig_dm * v ) ) ) ) return MfMui.T * ( C * ( MeSigmaIDeriv * v ) ) return C * ( MeSigmaIDeriv * v ) - # C * ( dMeSigmaI_dI * ( dMe_dsig * ( dsig_dm * v ) ) ) def getRHS(self, freq): @@ -456,22 +449,13 @@ class ProblemFDEM_j(BaseFDEMProblem): C = self.mesh.edgeCurl MfRhoDeriv_m = self.MfRhoDeriv(u) - # sigi = self.curModel.rho - # dsig_dm = self.curModel.rhoDeriv - # dsigi_dsig = -Utils.sdiag(sigi)**2 - # dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(u) - # MfRhoDeriv_m = (dMf_dsigi * ( dsigi_dsig * ( dsig_dm))) - if adjoint: if self._makeASymmetric is True: v = MfRho * v - # return dsig_dm.T * ( dsigi_dsig.T *( dMf_dsigi.T * ( C * ( MeMuI.T * ( C.T * v ) ) ) ) ) - return MfRhoDeriv_m.T * (C * (MeMuI.T * v)) + return MfRhoDeriv_m.T * (C * (MeMuI.T * (C.T * v))) if self._makeASymmetric is True: - # return MfRho.T * ( C * ( MeMuI * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) ) ) ) return MfRho.T * (C * ( MeMuI * (C.T * (MfRhoDeriv_m * v) ))) - # return C * ( MeMuI * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) ) ) return C * (MeMuI * (C.T * (MfRhoDeriv_m * v))) @@ -501,6 +485,7 @@ class ProblemFDEM_j(BaseFDEMProblem): if adjoint: if self._makeASymmetric: + MfRho = self.MfRho v = MfRho*v S_mDerivv = S_mDeriv(MeMuI.T * (C.T * v)) S_eDerivv = S_eDeriv(v) @@ -525,6 +510,7 @@ class ProblemFDEM_j(BaseFDEMProblem): return None if self._makeASymmetric: + MfRho = self.MfRho return MfRho.T * RHSDeriv return RHSDeriv @@ -580,17 +566,10 @@ class ProblemFDEM_h(BaseFDEMProblem): MeMu = self.MeMu C = self.mesh.edgeCurl - # sigi = self.sigmai - # dsig_dm = self.curModel.transformDeriv - # dsigi_dsig = -Utils.sdiag(sigi)**2 - MfRhoDeriv_m = self.MfRhoDeriv(C*u) - # dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(C*u) if adjoint: - # return (dsig_dm.T * (dsigi_dsig.T * (dMf_dsigi.T * (C * v)))) return MfRhoDeriv_m.T * (C * v) - # return (C.T * (dMf_dsigi * (dsigi_dsig * (dsig_dm * v)))) return C.T * (MfRhoDeriv_m * v) def getRHS(self, freq): @@ -623,8 +602,10 @@ class ProblemFDEM_h(BaseFDEMProblem): S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint) - RHSDeriv = C.T * (MfRhoDeriv * v) - # = S_m + C.T * ( MfRho * S_e ) + if not adjoint: + RHSDeriv = C.T * (MfRhoDeriv * v) + elif adjoint: + RHSDeriv = MfRhoDeriv.T * (C * v) S_mDeriv = S_mDeriv(v) S_eDeriv = S_eDeriv(v) diff --git a/simpegEM/FDEM/FieldsFDEM.py b/simpegEM/FDEM/FieldsFDEM.py index 2b599cdd..5838ddbe 100644 --- a/simpegEM/FDEM/FieldsFDEM.py +++ b/simpegEM/FDEM/FieldsFDEM.py @@ -206,11 +206,10 @@ class FieldsFDEM_j(FieldsFDEM): self._edgeCurl = self.survey.prob.mesh.edgeCurl self._MeMuI = self.survey.prob.MeMuI self._MfRho = self.survey.prob.MfRho - self._curModel = self.survey.prob.curModel self._MfRhoDeriv = self.survey.prob.MfRhoDeriv def _jPrimary(self, jSolution, srcList): - jPrimary = np.zeros_like(jSolution) + jPrimary = np.zeros_like(jSolution,dtype = complex) for i, src in enumerate(srcList): jp = src.jPrimary(self.survey.prob) if jp is not None: @@ -257,7 +256,7 @@ class FieldsFDEM_j(FieldsFDEM): if not adjoint: return -1./(1j*omega(src.freq)) * MeMuI * (C.T * (MfRho * v) ) elif adjoint: - return -1./(1j*omega(src.freq)) * MfRho * (C * ( MeMuI .T * v)) + return -1./(1j*omega(src.freq)) * MfRho.T * (C * ( MeMuI.T * v)) def _hSecondaryDeriv_m(self, src, v, adjoint=False): jSolution = self[[src],'jSolution'] @@ -281,36 +280,18 @@ class FieldsFDEM_j(FieldsFDEM): S_mDeriv = S_mDeriv(MeMuI.T * v) if S_mDeriv is not None: hDeriv_m += 1./(1j*omega(src.freq)) * S_mDeriv - return h + return hDeriv_m def _h(self, jSolution, srcList): return self._hPrimary(jSolution, srcList) + self._hSecondary(jSolution, srcList) - # raise NotImplementedError('Fields Derivs Not Implemented Yet') - # 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._MfRho - - # S_mDeriv,_ = src.getSourceDeriv(self.survey.prob, v, adjoint) - - # if not adjoint: - # h_Deriv= -(1./(1j*omega(freq))) * MeMuI * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) ) - # else: - # h_Deriv= -(1./(1j*omega(freq))) * dsig_dm.T * ( dsigi_dsig.T * ( dMf_dsigi.T * ( C * ( MeMuI.T * v ) ) ) ) - - # if S_mDeriv is not None: - # return 1./(1j*omega(src.freq)) * S_mDeriv + h_Deriv - def _hDeriv_u(self, src, v, adjoint=False): - return _hSecondaryDeriv_u(self, src, v, adjoint) + return self._hSecondaryDeriv_u(src, v, adjoint) def _hDeriv_m(self, src, v, adjoint=False): # assuming the primary doesn't depend on the model - return _hSecondaryDeriv_u(self, src, v, adjoint) + return self._hSecondaryDeriv_m(src, v, adjoint) class FieldsFDEM_h(FieldsFDEM): @@ -333,7 +314,7 @@ class FieldsFDEM_h(FieldsFDEM): self._MfRho = self.survey.prob.MfRho def _hPrimary(self, hSolution, srcList): - hPrimary = np.zeros_like(hSolution) + hPrimary = np.zeros_like(hSolution,dtype = complex) for i, src in enumerate(srcList): hp = src.hPrimary(self.survey.prob) if hp is not None: @@ -369,63 +350,25 @@ class FieldsFDEM_h(FieldsFDEM): j[:,i] += -S_e return j + def _jSecondaryDeriv_u(self, src, v, adjoint=False): + if not adjoint: + return self._edgeCurl*v + elif adjoint: + return self._edgeCurl.T*v + + def _jSecondaryDeriv_m(self, src, v, adjoint=False): + _,S_eDeriv = src.evalDeriv(self.survey.prob, adjoint) + S_eDeriv = S_eDeriv(v) + if S_eDeriv is not None: + return -S_eDeriv + return None + def _j(self, hSolution, srcList): return self._jPrimary(hSolution, srcList) + self._jSecondary(hSolution, srcList) - def _jDeriv(self, hSolution, srcList, v, adjoint=False): - raise NotImplementedError('Fields Derivs Not Implemented Yet') - _,S_eDeriv = src.getSourceDeriv(self.survey.prob, v, adjoint) - if S_eDeriv is None: - return None - else: - return - S_eDeriv + def _jDeriv_u(self, src, v, adjoint=False): + return self._jSecondaryDeriv_u(src,v,adjoint) - - # def calcFields(self, Solution, freq, fieldType, adjoint=False): - # j = Solution - # if fieldType == 'j': - # return j - # elif fieldType == 'h': - # MeMuI = self._MeMuI - # C = self.mesh.edgeCurl - # MfRho = self._MfRho - # if not adjoint: - # h = -(1./(1j*omega(freq))) * MeMuI * ( C.T * ( MfRho * j ) ) - # else: - # h = -(1./(1j*omega(freq))) * MfRho.T * ( C * ( MeMuI.T * j ) ) - # return h - # raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) - - # def calcFieldsDeriv(self, Solution, freq, fieldType, v, adjoint=False): - # j = Solution - # 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._MfRho - # 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, Solution, freq, fieldType, adjoint=False): - # h = Solution - # 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, Solution, freq, fieldType, v, adjoint=False): - # return None \ No newline at end of file + def _jDeriv_m(self, src, v, adjoint=False): + # assuming the primary does not depend on the model + return self._jSecondaryDeriv_m(src,v,adjoint) \ No newline at end of file diff --git a/simpegEM/Tests/test_FDEM.py b/simpegEM/Tests/test_FDEM.py index bc4e9807..e56a92ba 100644 --- a/simpegEM/Tests/test_FDEM.py +++ b/simpegEM/Tests/test_FDEM.py @@ -1,7 +1,7 @@ import unittest from SimPEG import * import simpegEM as EM -import sys +import sys from scipy.constants import mu_0 import copy @@ -9,7 +9,7 @@ testDerivs = True testCrossCheck = True testAdjoint = True testEB = True -testHJ = False +testHJ = True verbose = False @@ -35,8 +35,8 @@ def getProblem(fdemType, comp): 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.]) Rx0 = EM.FDEM.RxFDEM(XYZ, comp) - Src0 = EM.FDEM.SrcFDEM_MagDipole([Rx0], loc=np.r_[0.,0.,0.], freq=freq[0]) - Src1 = EM.FDEM.SrcFDEM_MagDipole([Rx0], loc=np.r_[0.,0.,0.], freq=freq[1]) + Src0 = EM.FDEM.SrcFDEM_MagDipole([Rx0], freq=freq[0], loc=np.r_[0.,0.,0.]) + Src1 = EM.FDEM.SrcFDEM_MagDipole([Rx0], freq=freq[1], loc=np.r_[0.,0.,0.]) survey = EM.FDEM.SurveyFDEM([Src0, Src1]) @@ -69,7 +69,7 @@ def adjointTest(fdemType, comp): print 'Adjoint %s formulation - %s' % (fdemType, comp) m = np.log(np.ones(prb.mesh.nC)*CONDUCTIVITY) - mu = np.log(np.ones(prb.mesh.nC)*MU) + mu = np.log(np.ones(prb.mesh.nC))*MU if addrandoms is True: m = m + np.random.randn(prb.mesh.nC)*CONDUCTIVITY*1e-1