From 27c8da341df6ad8364302ec9a280a06889479e66 Mon Sep 17 00:00:00 2001 From: Lindsey Date: Mon, 1 Jun 2015 14:52:39 -0700 Subject: [PATCH 1/7] start of using PropMaps for EB formulation --- simpegEM/Base.py | 165 ++++++++++++++++++++++++++--------------------- 1 file changed, 93 insertions(+), 72 deletions(-) diff --git a/simpegEM/Base.py b/simpegEM/Base.py index 7f938bf6..9e213825 100644 --- a/simpegEM/Base.py +++ b/simpegEM/Base.py @@ -1,6 +1,16 @@ -from SimPEG import Survey, Problem, Utils, Models, np, sp, Solver as SimpegSolver +from SimPEG import Survey, Problem, Utils, Models, Maps, PropMaps, np, sp, Solver as SimpegSolver from scipy.constants import mu_0 +class EMPropMap(Maps.PropMap): + sigma = Maps.Property("Electrical Conductivity", defaultInvProp = True) + mui = Maps.Property("Inverse Magnetic Permeability", defaultVal = 1./mu_0) + + # rho = Maps.Property("Electrical Resistivity") + # mu = Maps.Property("Inverse Magnetic Permeability", defaultVal = 1./mu_0) + + # Do some error checking: only 1 of sigma, rho can be InvProp similar story with mu and mui + # Also ensure that sigma and rho are reciprocals of one another "" + class BaseEMProblem(Problem.BaseProblem): def __init__(self, mesh, **kwargs): @@ -9,6 +19,8 @@ class BaseEMProblem(Problem.BaseProblem): surveyPair = Survey.BaseSurvey dataPair = Survey.Data + + PropMap = EMPropMap Solver = SimpegSolver solverOpts = {} @@ -29,89 +41,98 @@ class BaseEMProblem(Problem.BaseProblem): #################################################### # Mu - @property - def mu(self): - if getattr(self, '_mu', None) is None: - # if getattr(self, '_mui', None) is not None: - # self._mu = sel - self._mu = mu_0 - return self._mu - @mu.setter - 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 + # @property + # def mu(self): + # if getattr(self, '_mu', None) is None: + # # if getattr(self, '_mui', None) is not None: + # # self._mu = sel + # self._mu = mu_0 + # return self._mu + # @mu.setter + # 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 + # @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 + # @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 + # @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 #################################################### # Mass Matrices #################################################### + # TODO: Link to EMPropMap + # if Prop + # deleteTheseOnModelUpdate = ['_MeSigma', '_MeSigmaI','_MfSigmai','_MfSigmaiI'] + @property + def deleteTheseOnModelUpdate(self): + toDelete = [] + if self.mapping.sigmaMap is not None: + toDelete += ['_MeSigma', '_MeSigmaI','_MfSigmai','_MfSigmaiI'] + return toDelete + @property def Me(self): if getattr(self, '_Me', None) is None: @@ -129,13 +150,13 @@ class BaseEMProblem(Problem.BaseProblem): @property def MfMui(self): if getattr(self, '_MfMui', None) is None: - self._MfMui = self.mesh.getFaceInnerProduct(self.mui) + self._MfMui = self.mesh.getFaceInnerProduct(self.curModel.mui) return self._MfMui @property def MfMuiI(self): if getattr(self, '_MfMuiI', None) is None: - self._MfMuiI = self.mesh.getFaceInnerProduct(self.mui, invMat=True) + self._MfMuiI = self.mesh.getFaceInnerProduct(self.curModel.mui, invMat=True) return self._MfMuiI @property @@ -155,7 +176,7 @@ class BaseEMProblem(Problem.BaseProblem): @property def MeSigma(self): if getattr(self, '_MeSigma', None) is None: - self._MeSigma = self.mesh.getEdgeInnerProduct(self.sigma) + self._MeSigma = self.mesh.getEdgeInnerProduct(self.curModel.sigma) return self._MeSigma # def dMeSigma_dsigma(self, u): @@ -164,7 +185,7 @@ class BaseEMProblem(Problem.BaseProblem): @property def MeSigmaI(self): if getattr(self, '_MeSigmaI', None) is None: - self._MeSigmaI = self.mesh.getEdgeInnerProduct(self.sigma, invMat=True) + self._MeSigmaI = self.mesh.getEdgeInnerProduct(self.curModel.sigma, invMat=True) return self._MeSigmaI # def dMeSigmaI_dsigma(self,u) From b7a3b4e5e6555b70e8c47d691dd49daacb256fa7 Mon Sep 17 00:00:00 2001 From: Lindsey Date: Mon, 1 Jun 2015 16:21:15 -0700 Subject: [PATCH 2/7] Running EB, HJ with PropMap --- simpegEM/Base.py | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/simpegEM/Base.py b/simpegEM/Base.py index 9e213825..d9ea7984 100644 --- a/simpegEM/Base.py +++ b/simpegEM/Base.py @@ -2,11 +2,14 @@ from SimPEG import Survey, Problem, Utils, Models, Maps, PropMaps, np, sp, Solve from scipy.constants import mu_0 class EMPropMap(Maps.PropMap): - sigma = Maps.Property("Electrical Conductivity", defaultInvProp = True) - mui = Maps.Property("Inverse Magnetic Permeability", defaultVal = 1./mu_0) + sigma = Maps.Property("Electrical Conductivity", defaultInvProp = True, propertyLink=('rho',Maps.ReciprocalMap)) + mu = Maps.Property("Inverse Magnetic Permeability", defaultVal = mu_0, propertyLink=('mui',Maps.ReciprocalMap)) - # rho = Maps.Property("Electrical Resistivity") - # mu = Maps.Property("Inverse Magnetic Permeability", defaultVal = 1./mu_0) + rho = Maps.Property("Electrical Resistivity", propertyLink=('sigma', Maps.ReciprocalMap)) + mui = Maps.Property("Inverse Magnetic Permeability", defaultVal = 1./mu_0, propertyLink=('mu', Maps.ReciprocalMap)) + + + # Do some error checking: only 1 of sigma, rho can be InvProp similar story with mu and mui # Also ensure that sigma and rho are reciprocals of one another "" @@ -162,13 +165,13 @@ class BaseEMProblem(Problem.BaseProblem): @property def MeMu(self): if getattr(self, '_MeMu', None) is None: - self._MeMu = self.mesh.getEdgeInnerProduct(self.mu) + self._MeMu = self.mesh.getEdgeInnerProduct(self.curModel.mu) return self._MeMu @property def MeMuI(self): if getattr(self, '_MeMuI', None) is None: - self._MeMuI = self.mesh.getEdgeInnerProduct(self.mu, invMat=True) + self._MeMuI = self.mesh.getEdgeInnerProduct(self.curModel.mu, invMat=True) return self._MeMuI # ----- Electrical Conductivity ----- # @@ -193,7 +196,7 @@ class BaseEMProblem(Problem.BaseProblem): @property def MfSigmai(self): if getattr(self, '_MfSigmai', None) is None: - self._MfSigmai = self.mesh.getFaceInnerProduct(self.sigmai) + self._MfSigmai = self.mesh.getFaceInnerProduct(self.curModel.rho) return self._MfSigmai # def dMfSigmai_dsigmai(self,u) @@ -201,7 +204,7 @@ class BaseEMProblem(Problem.BaseProblem): @property def MfSigmaiI(self): if getattr(self, '_MfSigmaiI', None) is None: - self._MfSigmaiI = self.mesh.getFaceInnerProduct(self.sigmai, invMat=True) + self._MfSigmaiI = self.mesh.getFaceInnerProduct(self.curModel.rho, invMat=True) return self._MfSigmaiI # def dMfSigmaiI(self,u) From e0387978f29b907b4cbef636bfdec01153ddc518 Mon Sep 17 00:00:00 2001 From: Lindsey Date: Mon, 1 Jun 2015 16:32:23 -0700 Subject: [PATCH 3/7] Simgai --> Rho --- simpegEM/Base.py | 111 +++++------------------------------- simpegEM/FDEM/FDEM.py | 24 ++++---- simpegEM/FDEM/FieldsFDEM.py | 18 +++--- simpegEM/Utils/EMUtils.py | 6 +- 4 files changed, 38 insertions(+), 121 deletions(-) diff --git a/simpegEM/Base.py b/simpegEM/Base.py index d9ea7984..3bf0dec7 100644 --- a/simpegEM/Base.py +++ b/simpegEM/Base.py @@ -39,101 +39,18 @@ class BaseEMProblem(Problem.BaseProblem): self.__makeASymmetric = True return self.__makeASymmetric - #################################################### - # Phys Props - #################################################### - - # Mu - # @property - # def mu(self): - # if getattr(self, '_mu', None) is None: - # # if getattr(self, '_mui', None) is not None: - # # self._mu = sel - # self._mu = mu_0 - # return self._mu - # @mu.setter - # 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 - # @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 - #################################################### # Mass Matrices #################################################### - # TODO: Link to EMPropMap - # if Prop - # deleteTheseOnModelUpdate = ['_MeSigma', '_MeSigmaI','_MfSigmai','_MfSigmaiI'] @property def deleteTheseOnModelUpdate(self): toDelete = [] - if self.mapping.sigmaMap is not None: - toDelete += ['_MeSigma', '_MeSigmaI','_MfSigmai','_MfSigmaiI'] + if self.mapping.sigmaMap is not None or self.mapping.rhoMap is not None: + toDelete += ['_MeSigma', '_MeSigmaI','_MfRho','_MfRhoI'] + if self.mapping.muMap is not None or self.mapping.muiMap is not None: + toDelete += ['_MeMu', '_MeMuI','_MfMui','_MfMuiI'] return toDelete @property @@ -194,20 +111,20 @@ class BaseEMProblem(Problem.BaseProblem): # def dMeSigmaI_dsigma(self,u) @property - def MfSigmai(self): - if getattr(self, '_MfSigmai', None) is None: - self._MfSigmai = self.mesh.getFaceInnerProduct(self.curModel.rho) - return self._MfSigmai + def MfRho(self): + if getattr(self, '_MfRho', None) is None: + self._MfRho = self.mesh.getFaceInnerProduct(self.curModel.rho) + return self._MfRho - # def dMfSigmai_dsigmai(self,u) + # def dMfRho_dsigmai(self,u) @property - def MfSigmaiI(self): - if getattr(self, '_MfSigmaiI', None) is None: - self._MfSigmaiI = self.mesh.getFaceInnerProduct(self.curModel.rho, invMat=True) - return self._MfSigmaiI + def MfRhoI(self): + if getattr(self, '_MfRhoI', None) is None: + self._MfRhoI = self.mesh.getFaceInnerProduct(self.curModel.rho, invMat=True) + return self._MfRhoI - # def dMfSigmaiI(self,u) + # def dMfRhoI(self,u) #################################################### diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py index e83fe0e0..7ba70da2 100644 --- a/simpegEM/FDEM/FDEM.py +++ b/simpegEM/FDEM/FDEM.py @@ -327,14 +327,14 @@ class ProblemFDEM_j(BaseFDEMProblem): """ MeMuI = self.MeMuI - MfSigmai = self.MfSigmai + MfRho = self.MfRho C = self.mesh.edgeCurl iomega = 1j * omega(freq) * sp.eye(self.mesh.nF) - A = C * MeMuI * C.T * MfSigmai + iomega + A = C * MeMuI * C.T * MfRho + iomega if self._makeASymmetric is True: - return MfSigmai.T*A + return MfRho.T*A return A @@ -347,7 +347,7 @@ class ProblemFDEM_j(BaseFDEMProblem): """ MeMuI = self.MeMuI - MfSigmai = self.MfSigmai + MfRho = self.MfRho C = self.mesh.edgeCurl sigi = self.sigmai dsig_dm = self.curModel.transformDeriv @@ -356,11 +356,11 @@ class ProblemFDEM_j(BaseFDEMProblem): if adjoint: if self._makeASymmetric is True: - v = MfSigmai * v + v = MfRho * v return dsig_dm.T * ( dsigi_dsig.T *( dMf_dsigi.T * ( C * ( MeMuI.T * ( C.T * v ) ) ) ) ) if self._makeASymmetric is True: - return MfSigmai.T * ( C * ( MeMuI * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) ) ) ) + return MfRho.T * ( C * ( MeMuI * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) ) ) ) return C * ( MeMuI * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) ) ) @@ -378,8 +378,8 @@ class ProblemFDEM_j(BaseFDEMProblem): RHS = C * (MeMuI * S_m) - 1j * omega(freq) * S_e if self._makeASymmetric is True: - MfSigmai = self.MfSigmai - return MfSigmai.T*RHS + MfRho = self.MfRho + return MfRho.T*RHS return RHS @@ -429,10 +429,10 @@ class ProblemFDEM_h(BaseFDEMProblem): """ MeMu = self.MeMu - MfSigmai = self.MfSigmai + MfRho = self.MfRho C = self.mesh.edgeCurl - return C.T * MfSigmai * C + 1j*omega(freq)*MeMu + return C.T * MfRho * C + 1j*omega(freq)*MeMu def getADeriv(self, freq, u, v, adjoint=False): @@ -458,9 +458,9 @@ class ProblemFDEM_h(BaseFDEMProblem): S_m, S_e = self.getSourceTerm(freq) C = self.mesh.edgeCurl - MfSigmai = self.MfSigmai + MfRho = self.MfRho - RHS = S_m + C.T * ( MfSigmai * S_e ) + RHS = S_m + C.T * ( MfRho * S_e ) return RHS diff --git a/simpegEM/FDEM/FieldsFDEM.py b/simpegEM/FDEM/FieldsFDEM.py index e12752b6..9e93d3d5 100644 --- a/simpegEM/FDEM/FieldsFDEM.py +++ b/simpegEM/FDEM/FieldsFDEM.py @@ -114,7 +114,7 @@ class FieldsFDEM_j(FieldsFDEM): def startup(self): self._edgeCurl = self.survey.prob.mesh.edgeCurl self._MeMuI = self.survey.prob.MeMuI - self._MfSigmai = self.survey.prob.MfSigmai + self._MfRho = self.survey.prob.MfRho self._curModel = self.survey.prob.curModel def _j(self, j_sol, srcList): @@ -128,9 +128,9 @@ class FieldsFDEM_j(FieldsFDEM): def _h(self, j_sol, srcList): MeMuI = self._MeMuI C = self._edgeCurl - MfSigmai = self._MfSigmai + MfRho = self._MfRho - h = MeMuI * (C.T * (MfSigmai * j_sol) ) + h = MeMuI * (C.T * (MfRho * j_sol) ) for i, src in enumerate(srcList): h[:,i] *= -1./(1j*omega(src.freq)) @@ -151,7 +151,7 @@ class FieldsFDEM_j(FieldsFDEM): dsig_dm = self._curModel.transformDeriv dsigi_dsig = -Utils.sdiag(sigi)**2 dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(j) - sigi = self._MfSigmai + sigi = self._MfRho S_mDeriv,_ = src.getSourceDeriv(self.survey.prob, v, adjoint) @@ -177,7 +177,7 @@ class FieldsFDEM_h(FieldsFDEM): def startup(self): self._edgeCurl = self.survey.prob.mesh.edgeCurl self._MeMuI = self.survey.prob.MeMuI - self._MfSigmai = self.survey.prob.MfSigmai + self._MfRho = self.survey.prob.MfRho def _h(self, h_sol, srcList): h = h_sol @@ -215,11 +215,11 @@ class FieldsFDEM_h(FieldsFDEM): # elif fieldType == 'h': # MeMuI = self._MeMuI # C = self.mesh.edgeCurl - # MfSigmai = self._MfSigmai + # MfRho = self._MfRho # if not adjoint: - # h = -(1./(1j*omega(freq))) * MeMuI * ( C.T * ( MfSigmai * j ) ) + # h = -(1./(1j*omega(freq))) * MeMuI * ( C.T * ( MfRho * j ) ) # else: - # h = -(1./(1j*omega(freq))) * MfSigmai.T * ( C * ( MeMuI.T * j ) ) + # h = -(1./(1j*omega(freq))) * MfRho.T * ( C * ( MeMuI.T * j ) ) # return h # raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) @@ -235,7 +235,7 @@ class FieldsFDEM_h(FieldsFDEM): # dsig_dm = self._curModel.transformDeriv # dsigi_dsig = -Utils.sdiag(sigi)**2 # dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(j) - # sigi = self._MfSigmai + # sigi = self._MfRho # if not adjoint: # return -(1./(1j*omega(freq))) * MeMuI * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) ) # else: diff --git a/simpegEM/Utils/EMUtils.py b/simpegEM/Utils/EMUtils.py index 5cdaf150..916d6ae2 100644 --- a/simpegEM/Utils/EMUtils.py +++ b/simpegEM/Utils/EMUtils.py @@ -16,7 +16,7 @@ def e_from_j(prob,j): if eqLocs is 'FE': MSigmaI = prob.MeSigmaI elif eqLocs is 'EF': - MSigmaI = prob.MfSigmai + MSigmaI = prob.MfRho return MSigmaI*j def j_from_e(prob,e): @@ -24,13 +24,13 @@ def j_from_e(prob,e): if eqLocs is 'FE': MSigma = prob.MeSigma elif eqLocs is 'EF': - MSigma = prob.MfSigma + MSigma = prob.MfRhoI return MSigma*e def b_from_h(prob,h): eqLocs = prob._eqLocs if eqLocs is 'FE': - MMu = prob.MfMu + MMu = prob.MfMuiI elif eqLocs is 'EF': MMu = prob.MeMu return MMu*h From 80a7dfb51bee5747f8040bcc95a50a5d8095ba10 Mon Sep 17 00:00:00 2001 From: Lindsey Date: Mon, 1 Jun 2015 16:56:11 -0700 Subject: [PATCH 4/7] start of moving mass matrix derivs onto base.py --- simpegEM/Base.py | 10 ++++-- simpegEM/Tests/test_Base.py | 65 +++++++++++++++++++++++++++++++++++++ 2 files changed, 73 insertions(+), 2 deletions(-) create mode 100644 simpegEM/Tests/test_Base.py diff --git a/simpegEM/Base.py b/simpegEM/Base.py index 3bf0dec7..2c89fc59 100644 --- a/simpegEM/Base.py +++ b/simpegEM/Base.py @@ -99,8 +99,14 @@ class BaseEMProblem(Problem.BaseProblem): self._MeSigma = self.mesh.getEdgeInnerProduct(self.curModel.sigma) return self._MeSigma - # def dMeSigma_dsigma(self, u): - # return self.mesh.getEdgeInnerProductDeriv(self.sigma)(u) + @property + def MeSigmaDeriv(self, u): + """ + Deriv of MeSigma wrt sigma + """ + if getattr(self, 'MeSigmaDeriv', None) is None: + self._MeSigmaDeriv = self.mesh.getEdgeInnerProductDeriv(self.curModel.sigma)(u) + return self._MeSigmaDeriv @property def MeSigmaI(self): diff --git a/simpegEM/Tests/test_Base.py b/simpegEM/Tests/test_Base.py new file mode 100644 index 00000000..7d168148 --- /dev/null +++ b/simpegEM/Tests/test_Base.py @@ -0,0 +1,65 @@ +import unittest +from SimPEG import * +import simpegEM as EM +import sys +from scipy.constants import mu_0 +import copy + +testDerivs = False +testCrossCheck = True +testAdjoint = False +testEB = True +testHJ = True + +verbose = False + +TOL = 1e-4 +FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order +CONDUCTIVITY = 1e1 +MU = mu_0 +freq = [1e-1, 2e-1] +addrandoms = True + +def getProblem(fdemType): + cs = 5. + ncx, ncy, ncz = 6, 6, 6 + npad = 3 + hx = [(cs,npad,-1.3), (cs,ncx), (cs,npad,1.3)] + hy = [(cs,npad,-1.3), (cs,ncy), (cs,npad,1.3)] + hz = [(cs,npad,-1.3), (cs,ncz), (cs,npad,1.3)] + mesh = Mesh.TensorMesh([hx,hy,hz],['C','C','C']) + + mapping = Maps.ExpMap(mesh) + + 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.]) + Src0 = EM.FDEM.SrcFDEM_MagDipole([], loc=np.r_[0.,0.,0.], freq=freq[0]) + Src1 = EM.FDEM.SrcFDEM_MagDipole([], loc=np.r_[0.,0.,0.], freq=freq[1]) + + survey = EM.FDEM.SurveyFDEM([Src0]) + + + if verbose: + print ' Fetching %s problem' % (fdemType) + + if fdemType == 'e': + prb = EM.FDEM.ProblemFDEM_e(mesh, mapping=mapping) + elif fdemType == 'b': + prb = EM.FDEM.ProblemFDEM_b(mesh, mapping=mapping) + elif fdemType == 'j': + prb = EM.FDEM.ProblemFDEM_j(mesh, mapping=mapping) + elif fdemType == 'h': + prb = EM.FDEM.ProblemFDEM_h(mesh, mapping=mapping) + else: + raise NotImplementedError() + prb.pair(survey) + + try: + from pymatsolver import MumpsSolver + prb.Solver = MumpsSolver + except ImportError, e: + pass + + return prb + +def test_MassMatDeriv() From e53b1247241c857acc6f3dcd13bd9001af75f252 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 2 Jun 2015 13:54:14 -0700 Subject: [PATCH 5/7] Mass matrix derivs, and replaced call of forward with fields --- simpegEM/Base.py | 36 ++++++++------------ simpegEM/FDEM/FDEM.py | 12 ++++--- simpegEM/Tests/test_Base.py | 65 ------------------------------------- 3 files changed, 20 insertions(+), 93 deletions(-) delete mode 100644 simpegEM/Tests/test_Base.py diff --git a/simpegEM/Base.py b/simpegEM/Base.py index 2c89fc59..7f6c9766 100644 --- a/simpegEM/Base.py +++ b/simpegEM/Base.py @@ -8,11 +8,6 @@ class EMPropMap(Maps.PropMap): rho = Maps.Property("Electrical Resistivity", propertyLink=('sigma', Maps.ReciprocalMap)) mui = Maps.Property("Inverse Magnetic Permeability", defaultVal = 1./mu_0, propertyLink=('mu', Maps.ReciprocalMap)) - - - - # Do some error checking: only 1 of sigma, rho can be InvProp similar story with mu and mui - # Also ensure that sigma and rho are reciprocals of one another "" class BaseEMProblem(Problem.BaseProblem): @@ -99,22 +94,25 @@ class BaseEMProblem(Problem.BaseProblem): self._MeSigma = self.mesh.getEdgeInnerProduct(self.curModel.sigma) return self._MeSigma - @property def MeSigmaDeriv(self, u): """ Deriv of MeSigma wrt sigma - """ - if getattr(self, 'MeSigmaDeriv', None) is None: - self._MeSigmaDeriv = self.mesh.getEdgeInnerProductDeriv(self.curModel.sigma)(u) - return self._MeSigmaDeriv + """ + return self.mesh.getEdgeInnerProductDeriv(self.curModel.sigma)(u) + @property def MeSigmaI(self): if getattr(self, '_MeSigmaI', None) is None: self._MeSigmaI = self.mesh.getEdgeInnerProduct(self.curModel.sigma, invMat=True) return self._MeSigmaI - # def dMeSigmaI_dsigma(self,u) + def MeSigmaIDeriv(self, u): + """ + Deriv of MeSigma wrt sigma + """ + return self.mesh.getEdgeInnerProductDeriv(self.curModel.sigma, invMat=True)(u) + @property def MfRho(self): @@ -122,7 +120,8 @@ class BaseEMProblem(Problem.BaseProblem): self._MfRho = self.mesh.getFaceInnerProduct(self.curModel.rho) return self._MfRho - # def dMfRho_dsigmai(self,u) + def MfRhoDeriv(self,u): + return self.mesh.getFaceInnerProductDeriv(self.curModel.rho)(u) @property def MfRhoI(self): @@ -130,14 +129,5 @@ class BaseEMProblem(Problem.BaseProblem): self._MfRhoI = self.mesh.getFaceInnerProduct(self.curModel.rho, invMat=True) return self._MfRhoI - # def dMfRhoI(self,u) - - - #################################################### - # Fields - #################################################### - - def fields(self, m): - self.curModel = m - F = self.forward(m) - return F + def dMfRhoIDeriv(self,u): + return self.mesh.getFaceInnerProductDeriv(self.curModel.rho, invMat=True) diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py index 7ba70da2..6564183d 100644 --- a/simpegEM/FDEM/FDEM.py +++ b/simpegEM/FDEM/FDEM.py @@ -19,8 +19,8 @@ class BaseFDEMProblem(BaseEMProblem): surveyPair = SurveyFDEM fieldsPair = FieldsFDEM - def forward(self, m): - + def fields(self, m): + self.curModel = m F = self.fieldsPair(self.mesh, self.survey) for freq in self.survey.freqs: @@ -175,9 +175,11 @@ class ProblemFDEM_e(BaseFDEMProblem): def getADeriv(self, freq, u, v, adjoint=False): - sig = self.sigma - dsig_dm = self.curModel.transformDeriv - dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig)(u) + # sig = self.sigma + # dsig_dm = self.curModel.transformDeriv + # dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig)(u) + dsig_dm = self.curModel.sigmaDeriv + dMe_dsig = self.MeSigmaDeriv(u) if adjoint: return 1j * omega(freq) * ( dsig_dm.T * ( dMe_dsig.T * v ) ) diff --git a/simpegEM/Tests/test_Base.py b/simpegEM/Tests/test_Base.py deleted file mode 100644 index 7d168148..00000000 --- a/simpegEM/Tests/test_Base.py +++ /dev/null @@ -1,65 +0,0 @@ -import unittest -from SimPEG import * -import simpegEM as EM -import sys -from scipy.constants import mu_0 -import copy - -testDerivs = False -testCrossCheck = True -testAdjoint = False -testEB = True -testHJ = True - -verbose = False - -TOL = 1e-4 -FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order -CONDUCTIVITY = 1e1 -MU = mu_0 -freq = [1e-1, 2e-1] -addrandoms = True - -def getProblem(fdemType): - cs = 5. - ncx, ncy, ncz = 6, 6, 6 - npad = 3 - hx = [(cs,npad,-1.3), (cs,ncx), (cs,npad,1.3)] - hy = [(cs,npad,-1.3), (cs,ncy), (cs,npad,1.3)] - hz = [(cs,npad,-1.3), (cs,ncz), (cs,npad,1.3)] - mesh = Mesh.TensorMesh([hx,hy,hz],['C','C','C']) - - mapping = Maps.ExpMap(mesh) - - 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.]) - Src0 = EM.FDEM.SrcFDEM_MagDipole([], loc=np.r_[0.,0.,0.], freq=freq[0]) - Src1 = EM.FDEM.SrcFDEM_MagDipole([], loc=np.r_[0.,0.,0.], freq=freq[1]) - - survey = EM.FDEM.SurveyFDEM([Src0]) - - - if verbose: - print ' Fetching %s problem' % (fdemType) - - if fdemType == 'e': - prb = EM.FDEM.ProblemFDEM_e(mesh, mapping=mapping) - elif fdemType == 'b': - prb = EM.FDEM.ProblemFDEM_b(mesh, mapping=mapping) - elif fdemType == 'j': - prb = EM.FDEM.ProblemFDEM_j(mesh, mapping=mapping) - elif fdemType == 'h': - prb = EM.FDEM.ProblemFDEM_h(mesh, mapping=mapping) - else: - raise NotImplementedError() - prb.pair(survey) - - try: - from pymatsolver import MumpsSolver - prb.Solver = MumpsSolver - except ImportError, e: - pass - - return prb - -def test_MassMatDeriv() From 1424d071d4baf9f5a697307e66be87f12d58c6e4 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 2 Jun 2015 14:07:41 -0700 Subject: [PATCH 6/7] _p --> Primary --- simpegEM/FDEM/FieldsFDEM.py | 48 +++++++++++++++++------------------ simpegEM/FDEM/SurveyFDEM.py | 50 +++++++++++++------------------------ 2 files changed, 41 insertions(+), 57 deletions(-) diff --git a/simpegEM/FDEM/FieldsFDEM.py b/simpegEM/FDEM/FieldsFDEM.py index 9e93d3d5..ba4dbdd3 100644 --- a/simpegEM/FDEM/FieldsFDEM.py +++ b/simpegEM/FDEM/FieldsFDEM.py @@ -24,9 +24,9 @@ class FieldsFDEM_e(FieldsFDEM): def _e(self, e_sol, srcList): e = e_sol for i, src in enumerate(srcList): - e_p = src.e_p(self.survey.prob) - if e_p is not None: - e[:,i] += e_p + ePrimary = src.ePrimary(self.survey.prob) + if ePrimary is not None: + e[:,i] += ePrimary return e def _b(self, e_sol, srcList): @@ -38,9 +38,9 @@ class FieldsFDEM_e(FieldsFDEM): if S_m is not None: b[:,i] += 1./(1j*omega(src.freq)) * S_m - b_p = src.b_p(self.survey.prob) - if b_p is not None: - b[:,i] += b_p + bPrimary = src.bPrimary(self.survey.prob) + if bPrimary is not None: + b[:,i] += bPrimary return b @@ -72,9 +72,9 @@ class FieldsFDEM_b(FieldsFDEM): b = b_sol for i, src in enumerate(srcList): - b_p = src.b_p(self.survey.prob) - if b_p is not None: - b[:,i] += b_p + bPrimary = src.bPrimary(self.survey.prob) + if bPrimary is not None: + b[:,i] += bPrimary return b def _e(self, b_sol, srcList): @@ -85,9 +85,9 @@ class FieldsFDEM_b(FieldsFDEM): if S_e is not None: e += -self._MeSigmaI*S_e - e_p = src.e_p(self.survey.prob) - if e_p is not None: - e[:,i] += e_p + ePrimary = src.ePrimary(self.survey.prob) + if ePrimary is not None: + e[:,i] += ePrimary return e @@ -120,9 +120,9 @@ class FieldsFDEM_j(FieldsFDEM): def _j(self, j_sol, srcList): j = j_sol for i, src in enumerate(srcList): - j_p = src.j_p(self.survey.prob) - if j_p is not None: - j[:,i] += j_p + jPrimary = src.jPrimary(self.survey.prob) + if jPrimary is not None: + j[:,i] += jPrimary return j def _h(self, j_sol, srcList): @@ -138,9 +138,9 @@ class FieldsFDEM_j(FieldsFDEM): if S_m is not None: h[:,i] += 1./(1j*omega(src.freq)) * MeMuI * S_m - h_p = src.h_p(self.survey.prob) - if h_p is not None: - h[:,i] += h_p + hPrimary = src.hPrimary(self.survey.prob) + if hPrimary is not None: + h[:,i] += hPrimary return h @@ -182,9 +182,9 @@ class FieldsFDEM_h(FieldsFDEM): def _h(self, h_sol, srcList): h = h_sol for i, src in enumerate(srcList): - h_p = src.h_p(self.survey.prob) - if h_p is not None: - h[:,i] += h_p + hPrimary = src.hPrimary(self.survey.prob) + if hPrimary is not None: + h[:,i] += hPrimary return h def _j(self, h_sol, srcList): @@ -194,9 +194,9 @@ class FieldsFDEM_h(FieldsFDEM): if S_e is not None: j[:,i] += -S_e - j_p = src.j_p(self.survey.prob) - if j_p is not None: - j[:,i] += j_p + jPrimary = src.jPrimary(self.survey.prob) + if jPrimary is not None: + j[:,i] += jPrimary return j def _jDeriv(self, h_sol, srcList, v, adjoint=False): diff --git a/simpegEM/FDEM/SurveyFDEM.py b/simpegEM/FDEM/SurveyFDEM.py index ff1f50fe..56105903 100644 --- a/simpegEM/FDEM/SurveyFDEM.py +++ b/simpegEM/FDEM/SurveyFDEM.py @@ -103,32 +103,16 @@ class SrcFDEM(Survey.BaseSrc): def evalDeriv(self, prob, v, adjoint=None): return self._getS_mDeriv(prob,v,adjoint), self._getS_eDeriv(prob,v,adjoint) - def b_p(self,prob): - b_p = self._getb_p(prob) - return b_p + def bPrimary(self,prob): + return None - def h_p(self,prob): - h_p = self._geth_p(prob) - return h_p - - def e_p(self,prob): - e_p = self._gete_p(prob) - return e_p - - def j_p(self,prob): - j_p = self._getj_p(prob) - return j_p - - def _getb_p(self,prob): + def hPrimary(self,prob): return None - def _geth_p(self,prob): + def ePrimary(self,prob): return None - def _gete_p(self,prob): - return None - - def _getj_p(self,prob): + def jPrimary(self,prob): return None def _getS_m(self,prob): @@ -212,7 +196,7 @@ class SrcFDEM_MagDipole(SrcFDEM): self.moment = moment SrcFDEM.__init__(self, rxList) - def _getb_p(self,prob): + def bPrimary(self,prob): eqLocs = prob._eqLocs if eqLocs is 'FE': @@ -243,12 +227,12 @@ class SrcFDEM_MagDipole(SrcFDEM): return C*a - def _geth_p(self,prob): - b = self._getb_p(prob) + def hPrimary(self,prob): + b = self.bPrimary(prob) return h_from_b(prob,b) def _getS_m(self,prob): - b_p = self._getb_p(prob) + b_p = self.bPrimary(prob) return -1j*omega(self.freq)*b_p @@ -264,7 +248,7 @@ class SrcFDEM_MagDipole_Bfield(SrcFDEM): self.moment = moment SrcFDEM.__init__(self, rxList) - def _getb_p(self,prob): + def bPrimary(self,prob): eqLocs = prob._eqLocs if eqLocs is 'FE': @@ -295,12 +279,12 @@ class SrcFDEM_MagDipole_Bfield(SrcFDEM): return b - def _geth_p(self,prob): - b = self._getb_p(prob) + def hPrimary(self,prob): + b = self.bPrimary(prob) return h_from_b(prob, b) def _getS_m(self,prob): - b = self._getb_p(prob) + b = self.bPrimary(prob) return -1j*omega(self.freq)*b @@ -313,7 +297,7 @@ class SrcFDEM_CircularLoop(SrcFDEM): self.radius = radius SrcFDEM.__init__(self, rxList) - def _getb_p(self,prob): + def bPrimary(self,prob): eqLocs = prob._eqLocs if eqLocs is 'FE': @@ -343,12 +327,12 @@ class SrcFDEM_CircularLoop(SrcFDEM): return C*a - def _geth_p(self,prob): - b = self._getb_p(prob) + def hPrimary(self,prob): + b = self.bPrimary(prob) return h_from_b def _getS_m(self, prob): - b = self._getb_p(prob) + b = self.bPrimary(prob) return -1j*omega(self.freq)*b From 66b26a1b5e1aac493233f7c757d35d5691e761cb Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 2 Jun 2015 14:14:34 -0700 Subject: [PATCH 7/7] _getS_m now S_m --- simpegEM/FDEM/SurveyFDEM.py | 44 ++++++++++++++++++------------------- 1 file changed, 22 insertions(+), 22 deletions(-) diff --git a/simpegEM/FDEM/SurveyFDEM.py b/simpegEM/FDEM/SurveyFDEM.py index 56105903..e2df08c9 100644 --- a/simpegEM/FDEM/SurveyFDEM.py +++ b/simpegEM/FDEM/SurveyFDEM.py @@ -96,12 +96,12 @@ class SrcFDEM(Survey.BaseSrc): rxPair = RxFDEM def eval(self, prob): - S_m = self._getS_m(prob) - S_e = self._getS_e(prob) + S_m = self.S_m(prob) + S_e = self.S_e(prob) return S_m, S_e def evalDeriv(self, prob, v, adjoint=None): - return self._getS_mDeriv(prob,v,adjoint), self._getS_eDeriv(prob,v,adjoint) + return self.S_mDeriv(prob,v,adjoint), self.S_eDeriv(prob,v,adjoint) def bPrimary(self,prob): return None @@ -115,16 +115,16 @@ class SrcFDEM(Survey.BaseSrc): def jPrimary(self,prob): return None - def _getS_m(self,prob): + def S_m(self,prob): return None - def _getS_e(self,prob): + def S_e(self,prob): return None - def _getS_mDeriv(self, prob, v, adjoint = False): + def S_mDeriv(self, prob, v, adjoint = False): return None - def _getS_eDeriv(self, prob, v, adjoint = False): + def S_eDeriv(self, prob, v, adjoint = False): return None @@ -138,12 +138,12 @@ class SrcFDEM_RawVec_e(SrcFDEM): """ def __init__(self, rxList, freq, S_e): - self.S_e = np.array(S_e,dtype=float) + self._S_e = np.array(S_e,dtype=float) self.freq = float(freq) SrcFDEM.__init__(self, rxList) - def _getS_e(self, prob): - return self.S_e + def S_e(self, prob): + return self._S_e class SrcFDEM_RawVec_m(SrcFDEM): @@ -156,12 +156,12 @@ class SrcFDEM_RawVec_m(SrcFDEM): """ def __init__(self, rxList, freq, S_m): - self.S_m = np.array(S_m,dtype=float) + self._S_m = np.array(S_m,dtype=float) self.freq = float(freq) SrcFDEM.__init__(self, rxList) - def _getS_m(self, prob): - return self.S_m + def S_m(self, prob): + return self._S_m class SrcFDEM_RawVec(SrcFDEM): @@ -174,16 +174,16 @@ class SrcFDEM_RawVec(SrcFDEM): :param rxList: receiver list """ def __init__(self, rxList, freq, S_m, S_e): - self.S_m = np.array(S_m,dtype=float) - self.S_e = np.array(S_e,dtype=float) + self._S_m = np.array(S_m,dtype=float) + self._S_e = np.array(S_e,dtype=float) self.freq = float(freq) SrcFDEM.__init__(self, rxList) - def _getS_m(self,prob): - return self.S_m + def S_m(self,prob): + return self._S_m - def _getS_e(self,prob): - return self.S_e + def S_e(self,prob): + return self._S_e class SrcFDEM_MagDipole(SrcFDEM): @@ -231,7 +231,7 @@ class SrcFDEM_MagDipole(SrcFDEM): b = self.bPrimary(prob) return h_from_b(prob,b) - def _getS_m(self,prob): + def S_m(self,prob): b_p = self.bPrimary(prob) return -1j*omega(self.freq)*b_p @@ -283,7 +283,7 @@ class SrcFDEM_MagDipole_Bfield(SrcFDEM): b = self.bPrimary(prob) return h_from_b(prob, b) - def _getS_m(self,prob): + def S_m(self,prob): b = self.bPrimary(prob) return -1j*omega(self.freq)*b @@ -331,7 +331,7 @@ class SrcFDEM_CircularLoop(SrcFDEM): b = self.bPrimary(prob) return h_from_b - def _getS_m(self, prob): + def S_m(self, prob): b = self.bPrimary(prob) return -1j*omega(self.freq)*b