From 2864a976c6c781e60bad67b3d2773234be8b97a6 Mon Sep 17 00:00:00 2001 From: Lindsey Date: Fri, 17 Apr 2015 11:21:07 -0700 Subject: [PATCH] derivatives for fields objects, UNTESTED. TODO: clean up how we get derivs of mass matrices wrt the physical properties --- simpegEM/Base.py | 71 +++++++++++++++++++++++++------------ simpegEM/FDEM/FDEM.py | 3 +- simpegEM/FDEM/FieldsFDEM.py | 68 +++++++++++++++++------------------ 3 files changed, 83 insertions(+), 59 deletions(-) diff --git a/simpegEM/Base.py b/simpegEM/Base.py index e88fdc5d..dc8dec51 100644 --- a/simpegEM/Base.py +++ b/simpegEM/Base.py @@ -48,59 +48,84 @@ class BaseEMProblem(Problem.BaseProblem): # Mass Matrices #################################################### - @property - def MfMui(self): - if getattr(self, '_MfMui', None) is None: - 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: self._Me = self.mesh.getEdgeInnerProduct() return self._Me + @property + def Mf(self): + if getattr(self, '_Mf', None) is None: + self._Mf = self.mesh.getFaceInnerProduct() + return self._Mf + + + # ----- Magnetic Permeability ----- # + @property + def MfMui(self): + # TODO: hardcoded to assume diagonal mu + if getattr(self, '_MfMui', None) is None: + self._MfMui = self.mesh.getFaceInnerProduct(1/self.mu) + return self._MfMui + + @property + def MeMuI(self): + if getattr(self, '_MeMuI', None) is None: + self._MeMuI = self.mesh.getEdgeInnerProduct(self.mu, invMat=True) + return self._MeMuI + + @property + def MeMu(self): + if getattr(self, '_MeMu', None) is None: + self._MeMu = self.mesh.getEdgeInnerProduct(self.mu) + return self._MeMu + + + # ----- Electrical Conductivity ----- # + #TODO: hardcoded to sigma as the model @property def MeSigma(self): - #TODO: hardcoded to sigma as the model if getattr(self, '_MeSigma', None) is None: sigma = self.curModel.transform self._MeSigma = self.mesh.getEdgeInnerProduct(sigma) return self._MeSigma + @property def MeSigmaI(self): - #TODO: hardcoded to sigma as the model if getattr(self, '_MeSigmaI', None) is None: sigma = self.curModel.transform self._MeSigmaI = self.mesh.getEdgeInnerProduct(sigma, invMat=True) return self._MeSigmaI + @property + def dMeSigmaI_dI(self): + # TODO: hardcoded that sigma is diagonal + if getattr(self, '_dMeSigmaI_dI', None) is None: + self._dMeSigmaI_dI = - self.MeSigmaI**2 + return self._dMeSigmaI_dI + @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 + @property + def dMfSigmai_dsig(self): + return self._dMfSigmai_dsig + + deleteTheseOnModelUpdate = ['_MeSigma', '_MeSigmaI','_MfSigmai'] + + #################################################### + # Fields + #################################################### + def fields(self, m): self.curModel = m F = self.forward(m, self.getRHS) diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py index 870558b6..0fe726af 100644 --- a/simpegEM/FDEM/FDEM.py +++ b/simpegEM/FDEM/FDEM.py @@ -232,8 +232,7 @@ class ProblemFDEM_b(BaseFDEMProblem): C = self.mesh.edgeCurl sig = self.curModel.transform dsig_dm = self.curModel.transformDeriv - #TODO: This only works if diagonal (no tensors)... - dMeSigmaI_dI = - self.MeSigmaI**2 + dMeSigmaI_dI = self._dMeSigmaI_dI vec = (C.T*(mui*u)) dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig)(vec) diff --git a/simpegEM/FDEM/FieldsFDEM.py b/simpegEM/FDEM/FieldsFDEM.py index bda8accf..9c27c065 100644 --- a/simpegEM/FDEM/FieldsFDEM.py +++ b/simpegEM/FDEM/FieldsFDEM.py @@ -6,7 +6,7 @@ class FieldsFDEM(Problem.Fields): """Fancy Field Storage for a FDEM survey.""" knownFields = {} dtype = complex - + class FieldsFDEM_e(FieldsFDEM): knownFields = {'e':'E'} @@ -26,7 +26,7 @@ class FieldsFDEM_e(FieldsFDEM): def _b_sec(self, e, tx): #adjoint=False return - 1./(1j*omega(tx.freq)) * (self.edgeCurl * e) - def _b_secDeriv(self, e, tx, adjoint=False): + def _b_secDeriv(self, e, tx, v, adjoint=False): return None def _b(self, e, tx): #adjoint=False @@ -34,9 +34,9 @@ class FieldsFDEM_e(FieldsFDEM): S_m,_ = self.getSource(tx.freq) return b_sec + 1./(1j*omega(tx.freq)) * S_m - def _bDeriv(self, e, tx, adjoint=False): - S_mDeriv,_ = self.getSourceDeriv(tx.freq, adjoint) - b_secDeriv = self._b_secDeriv(e,tx.freq,adjoint) + def _bDeriv(self, e, tx, v, adjoint=False): + S_mDeriv,_ = self.getSourceDeriv(tx.freq, v, adjoint) + b_secDeriv = self._b_secDeriv(e, tx.freq, v, adjoint) if S_mDeriv is None & b_secDeriv is None: return None elif b_secDeriv is None: @@ -67,7 +67,7 @@ class FieldsFDEM_b(FieldsFDEM): def _e_sec(self, b, tx): return self.MeSigmaI * ( self.edgeCurl.T * ( self.MfMui * b) ) - def _e_secDeriv(self, b, tx, adjoint=False): + def _e_secDeriv(self, b, tx, v, adjoint=False): return None def _e(self, b, tx): @@ -75,9 +75,9 @@ class FieldsFDEM_b(FieldsFDEM): _, S_e = self.getSource(tx.freq) return e_sec + S_e - def _eDeriv(self, b, tx, adjoint=False): - _,S_eDeriv = self.getSourceDeriv(tx.freq, adjoint) - e_secDeriv = self._e_secDeriv(b,tx,adjoint) + def _eDeriv(self, b, tx, v, adjoint=False): + _,S_eDeriv = self.getSourceDeriv(tx.freq, v, adjoint) + e_secDeriv = self._e_secDeriv(b, tx, v, adjoint) if S_eDeriv is None & e_secDeriv is None: return None @@ -105,33 +105,33 @@ class FieldsFDEM_j(FieldsFDEM): self.MfSigmai = self.survey.prob.MfSigmai self.getSource = self.survey.prob.getSource self.getSourceDeriv = self.survey.prob.getSourceDeriv + self.curModel = self.prob.curModel - def _h_sec(self, j, tx): #adjoint=False + def _h_sec(self, j, tx): #v, 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_secDeriv(self, j, tx, v, adjoint=False): + MeMuI = self.MeMuI + C = self.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 ) ) ) ) - def _h(self, j, tx): #adjoint=False + def _h(self, j, tx): #v, adjoint=False h_sec = self._h_sec(j,tx) S_m,_ = self.getSource(tx.freq) return h_sec + 1./(1j*omega(tx.freq)) * self.MeMuI * S_m - def _hDeriv(self, j, tx, adjoint=False): - S_mDeriv,_ = self.getSourceDeriv(tx.freq, adjoint) - h_secDeriv = self._h_secDeriv(j,tx.freq,adjoint) + def _hDeriv(self, j, tx, v, adjoint=False): + S_mDeriv,_ = self.getSourceDeriv(tx.freq, v, adjoint) + h_secDeriv = self._h_secDeriv(j,tx.freq, v, adjoint) if S_mDeriv is None & h_secDeriv is None: return None elif h_secDeriv is None: @@ -158,20 +158,20 @@ class FieldsFDEM_h(FieldsFDEM): self.getSource = self.survey.prob.getSource self.getSourceDeriv = self.survey.prob.getSourceDeriv - def _j_sec(self, h, tx): #adjoint=False + def _j_sec(self, h, tx): # adjoint=False return self.edgeCurl*h - def _j_secDeriv(self, h, tx, adjoint=False): + def _j_secDeriv(self, h, tx, v, adjoint=False): return None - def _j(self, h, tx): #adjoint=False + def _j(self, h, tx): # adjoint=False j_sec = self._j_sec(h,tx) _,S_e = self.getSource(tx.freq) return j_sec - S_e - def _jDeriv(self, h, tx, adjoint=False): - _,S_eDeriv = self.getSourceDeriv(tx.freq, adjoint) - j_secDeriv = self._j_secDeriv(j,tx.freq,adjoint) + def _jDeriv(self, h, tx, v, adjoint=False): + _,S_eDeriv = self.getSourceDeriv(tx.freq, v, adjoint) + j_secDeriv = self._j_secDeriv(j,tx.freq, v, adjoint) if S_eDeriv is None & j_secDeriv is None: return None elif j_secDeriv is None: