diff --git a/simpegEM/Base.py b/simpegEM/Base.py index 00b628ca..27de846c 100644 --- a/simpegEM/Base.py +++ b/simpegEM/Base.py @@ -25,8 +25,10 @@ class BaseEMProblem(Problem.BaseProblem): return self.__makeASymmetric #################################################### - # Mu Model + # Phys Props #################################################### + + # Mu @property def mu(self): if getattr(self, '_mu', None) is None: @@ -36,12 +38,72 @@ class BaseEMProblem(Problem.BaseProblem): def mu(self, value): if getattr(self, '_MfMui', None) is not None: del self._MfMui + if getattr(self, '_MfMuiI', None) is not None: + del self._MfMuiI 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 - + + # TODO: hardcoded to assume diagonal mu + @property + def mui(self): + if getattr(self, '_mui', None) is None: + self._mui = 1./mu_0 + return self._mui + @mui.setter + def mui(self, value): + if getattr(self, '_MfMui', None) is not None: + del self._MfMui + if getattr(self, '_MfMuiI', None) is not None: + del self._MfMuiI + if getattr(self, '_MeMu', None) is not None: + del delf._MeMu + if getattr(self, '_MeMuI', None) is not None: + del self._MeMuI + self._mui = value + + # Sigma + # deleteTheseOnModelUpdate = ['_MeSigma', '_MeSigmaI','_MfSigmai','_MfSigmaiI'] + @property + def sigma(self): + if getattr(self, '_sigma', None) is None: + self._sigma = self.curModel.transform + return self._sigma + @sigma.setter + def sigma(self, value): + if getattr(self, '_MeSigma', None) is not None: + del self._MeSigma + if getattr(self, '_MeSigmaI', None) is not None: + del self._MeSigmaI + if getattr(self, '_MfSigmai', None) is not None: + del delf._MfSigmai + if getattr(self, '_MfSigmaiI', None) is not None: + del self._MfSigmaiI + self._sigma = value + + # def dsigma_dm(self): + # return self.curModel.transformDeriv + + + # TODO: hardcoded to assume diagonal sigma + @property + def sigmai(self): + if getattr(self, '_sigmai', None) is None: + self._sigmai = 1./self.curModel.transform + return self._sigmai + @sigmai.setter + def sigmai(self, value): + if getattr(self, '_MeSigma', None) is not None: + del self._MeSigma + if getattr(self, '_MeSigmaI', None) is not None: + del self._MeSigmaI + if getattr(self, '_MfSigmai', None) is not None: + del delf._MfSigmai + if getattr(self, '_MfSigmaiI', None) is not None: + del self._MfSigmaiI + self._sigma = value #################################################### @@ -64,16 +126,15 @@ class BaseEMProblem(Problem.BaseProblem): # ----- 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) + self._MfMui = self.mesh.getFaceInnerProduct(self.mui) 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 + def MfMuiI(self): + if getattr(self, '_MfMuiI', None) is None: + self._MfMuiI = self.mesh.getFaceInnerProduct(self.mui, invMat=True) + return self._MfMuiI @property def MeMu(self): @@ -81,43 +142,48 @@ class BaseEMProblem(Problem.BaseProblem): self._MeMu = self.mesh.getEdgeInnerProduct(self.mu) return self._MeMu + @property + def MeMuI(self): + if getattr(self, '_MeMuI', None) is None: + self._MeMuI = self.mesh.getEdgeInnerProduct(self.mu, invMat=True) + return self._MeMuI # ----- Electrical Conductivity ----- # #TODO: hardcoded to sigma as the model @property def MeSigma(self): if getattr(self, '_MeSigma', None) is None: - sigma = self.curModel.transform - self._MeSigma = self.mesh.getEdgeInnerProduct(sigma) + self._MeSigma = self.mesh.getEdgeInnerProduct(self.sigma) return self._MeSigma + # def dMeSigma_dsigma(self, u): + # return self.mesh.getEdgeInnerProductDeriv(self.sigma)(u) @property def MeSigmaI(self): if getattr(self, '_MeSigmaI', None) is None: - sigma = self.curModel.transform - self._MeSigmaI = self.mesh.getEdgeInnerProduct(sigma, invMat=True) + self._MeSigmaI = self.mesh.getEdgeInnerProduct(self.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 + # def dMeSigmaI_dsigma(self,u) @property def MfSigmai(self): - #TODO: hardcoded to sigma diagonal if getattr(self, '_MfSigmai', None) is None: - sigma = self.curModel.transform - self._MfSigmai = self.mesh.getFaceInnerProduct(1/sigma) + self._MfSigmai = self.mesh.getFaceInnerProduct(self.sigmai) return self._MfSigmai + + # def dMfSigmai_dsigmai(self,u) + + @property + def MfSigmaiI(self): + if getattr(self, '_MfSigmaiI', None) is None: + self._MfSigmaiI = self.mesh.getFaceInnerProduct(self.sigmai, invMat=True) + return self._MfSigmaiI + + # def dMfSigmaiI(self,u) - deleteTheseOnModelUpdate = ['_MeSigma', '_MeSigmaI','_MfSigmai'] - - #################################################### # Fields #################################################### diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py index c9020c80..e83fe0e0 100644 --- a/simpegEM/FDEM/FDEM.py +++ b/simpegEM/FDEM/FDEM.py @@ -167,15 +167,15 @@ class ProblemFDEM_e(BaseFDEMProblem): :rtype: scipy.sparse.csr_matrix :return: A """ - mui = self.MfMui - sig = self.MeSigma + MfMui = self.MfMui + MeSigma = self.MeSigma C = self.mesh.edgeCurl - return C.T*mui*C + 1j*omega(freq)*sig + return C.T*MfMui*C + 1j*omega(freq)*MeSigma def getADeriv(self, freq, u, v, adjoint=False): - sig = self.curModel.transform + sig = self.sigma dsig_dm = self.curModel.transformDeriv dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig)(u) @@ -221,35 +221,35 @@ class ProblemFDEM_b(BaseFDEMProblem): :rtype: scipy.sparse.csr_matrix :return: A """ - mui = self.MfMui - sigI = self.MeSigmaI + MfMui = self.MfMui + MeSigmaI = self.MeSigmaI C = self.mesh.edgeCurl iomega = 1j * omega(freq) * sp.eye(self.mesh.nF) - A = C*sigI*C.T*mui + iomega + A = C*MeSigmaI*C.T*MfMui + iomega if self._makeASymmetric is True: - return mui.T*A + return MfMui.T*A return A def getADeriv(self, freq, u, v, adjoint=False): - mui = self.MfMui + MfMui = self.MfMui C = self.mesh.edgeCurl - sig = self.curModel.transform + sig = self.sigma dsig_dm = self.curModel.transformDeriv dMeSigmaI_dI = self._dMeSigmaI_dI - vec = (C.T*(mui*u)) + vec = (C.T*(MfMui*u)) dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig)(vec) if adjoint: if self._makeASymmetric is True: - v = mui * v + v = MfMui * 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 MfMui.T * ( C * ( dMeSigmaI_dI * ( dMe_dsig * ( dsig_dm * v ) ) ) ) return C * ( dMeSigmaI_dI * ( dMe_dsig * ( dsig_dm * v ) ) ) @@ -267,8 +267,8 @@ class ProblemFDEM_b(BaseFDEMProblem): RHS = S_m + C * ( MeSigmaI * S_e ) if self._makeASymmetric is True: - mui = self.MfMui - return mui.T*RHS + MfMui = self.MfMui + return MfMui.T*RHS return RHS @@ -327,14 +327,14 @@ class ProblemFDEM_j(BaseFDEMProblem): """ MeMuI = self.MeMuI - MfSigi = self.MfSigmai + MfSigmai = self.MfSigmai C = self.mesh.edgeCurl iomega = 1j * omega(freq) * sp.eye(self.mesh.nF) - A = C * MeMuI * C.T * MfSigi + iomega + A = C * MeMuI * C.T * MfSigmai + iomega if self._makeASymmetric is True: - return MfSigi.T*A + return MfSigmai.T*A return A @@ -347,21 +347,20 @@ class ProblemFDEM_j(BaseFDEMProblem): """ MeMuI = self.MeMuI - MfSigi = self.MfSigmai + MfSigmai = self.MfSigmai C = self.mesh.edgeCurl - sig = self.curModel.transform - sigi = 1/sig + sigi = self.sigmai 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 + v = MfSigmai * 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 MfSigmai.T * ( C * ( MeMuI * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) ) ) ) return C * ( MeMuI * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) ) ) @@ -379,8 +378,8 @@ class ProblemFDEM_j(BaseFDEMProblem): RHS = C * (MeMuI * S_m) - 1j * omega(freq) * S_e if self._makeASymmetric is True: - MfSigi = self.MfSigmai - return MfSigi.T*RHS + MfSigmai = self.MfSigmai + return MfSigmai.T*RHS return RHS @@ -430,17 +429,16 @@ class ProblemFDEM_h(BaseFDEMProblem): """ MeMu = self.MeMu - MfSigi = self.MfSigmai + MfSigmai = self.MfSigmai C = self.mesh.edgeCurl - return C.T * MfSigi * C + 1j*omega(freq)*MeMu + return C.T * MfSigmai * 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 + sigi = self.sigmai dsig_dm = self.curModel.transformDeriv dsigi_dsig = -Utils.sdiag(sigi)**2