From e53b1247241c857acc6f3dcd13bd9001af75f252 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 2 Jun 2015 13:54:14 -0700 Subject: [PATCH] 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()