From d11f5c736ba8766439cb30db33818068e33ff148 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 29 Mar 2016 09:59:13 -0700 Subject: [PATCH] prim-sec for EB formulation on same meshes --- SimPEG/EM/Base.py | 30 +++++----- SimPEG/EM/FDEM/SrcFDEM.py | 58 +++++++++++-------- SimPEG/EM/Utils/testingUtils.py | 18 +++--- .../fdem/inverse/derivs/test_FDEM_derivs.py | 8 +-- 4 files changed, 65 insertions(+), 49 deletions(-) diff --git a/SimPEG/EM/Base.py b/SimPEG/EM/Base.py index 32018f7e..3d4acfa0 100644 --- a/SimPEG/EM/Base.py +++ b/SimPEG/EM/Base.py @@ -2,14 +2,14 @@ from SimPEG import Survey, Problem, Utils, Models, Maps, PropMaps, np, sp, Solve from scipy.constants import mu_0 class EMPropMap(Maps.PropMap): - """ + """ Property Map for EM Problems. The electrical conductivity (\\(\\sigma\\)) is the default inversion property, and the default value of the magnetic permeability is that of free space (\\(\\mu = 4\\pi\\times 10^{-7} \\) H/m) """ sigma = Maps.Property("Electrical Conductivity", defaultInvProp = True, propertyLink=('rho',Maps.ReciprocalMap)) - mu = Maps.Property("Inverse Magnetic Permeability", defaultVal = mu_0, propertyLink=('mui',Maps.ReciprocalMap)) + mu = Maps.Property("Magnetic Permeability", defaultVal = mu_0, propertyLink=('mui',Maps.ReciprocalMap)) - rho = Maps.Property("Electrical Resistivity", propertyLink=('sigma', Maps.ReciprocalMap)) + rho = Maps.Property("Electrical Resistivity", propertyLink=('sigma', Maps.ReciprocalMap)) mui = Maps.Property("Inverse Magnetic Permeability", defaultVal = 1./mu_0, propertyLink=('mu', Maps.ReciprocalMap)) @@ -21,7 +21,7 @@ class BaseEMProblem(Problem.BaseProblem): surveyPair = Survey.BaseSurvey dataPair = Survey.Data - + PropMap = EMPropMap Solver = SimpegSolver @@ -51,7 +51,7 @@ class BaseEMProblem(Problem.BaseProblem): if self.mapping.muMap is not None or self.mapping.muiMap is not None: toDelete += ['_MeMu', '_MeMuI','_MfMui','_MfMuiI'] return toDelete - + @property def Me(self): """ @@ -71,7 +71,7 @@ class BaseEMProblem(Problem.BaseProblem): return self._Mf - # ----- Magnetic Permeability ----- # + # ----- Magnetic Permeability ----- # @property def MfMui(self): """ @@ -109,7 +109,7 @@ class BaseEMProblem(Problem.BaseProblem): return self._MeMuI - # ----- Electrical Conductivity ----- # + # ----- Electrical Conductivity ----- # #TODO: hardcoded to sigma as the model @property def MeSigma(self): @@ -120,18 +120,18 @@ class BaseEMProblem(Problem.BaseProblem): self._MeSigma = self.mesh.getEdgeInnerProduct(self.curModel.sigma) return self._MeSigma - # TODO: This should take a vector + # TODO: This should take a vector def MeSigmaDeriv(self, u): """ Derivative of MeSigma with respect to the model - """ + """ return self.mesh.getEdgeInnerProductDeriv(self.curModel.sigma)(u) * self.curModel.sigmaDeriv - + @property def MeSigmaI(self): """ - Inverse of the edge inner product matrix for \\(\\sigma\\). + Inverse of the edge inner product matrix for \\(\\sigma\\). """ if getattr(self, '_MeSigmaI', None) is None: self._MeSigmaI = self.mesh.getEdgeInnerProduct(self.curModel.sigma, invMat=True) @@ -140,8 +140,8 @@ class BaseEMProblem(Problem.BaseProblem): # TODO: This should take a vector def MeSigmaIDeriv(self, u): """ - Derivative of :code:`MeSigma` with respect to the model - """ + Derivative of :code:`MeSigma` with respect to the model + """ # TODO: only works for diagonal tensors. getEdgeInnerProductDeriv, invMat=True should be implemented in SimPEG dMeSigmaI_dI = -self.MeSigmaI**2 @@ -163,7 +163,7 @@ class BaseEMProblem(Problem.BaseProblem): # TODO: This should take a vector def MfRhoDeriv(self,u): """ - Derivative of :code:`MfRho` with respect to the model. + Derivative of :code:`MfRho` with respect to the model. """ return self.mesh.getFaceInnerProductDeriv(self.curModel.rho)(u) * (-Utils.sdiag(self.curModel.rho**2) * self.curModel.sigmaDeriv) # self.curModel.rhoDeriv @@ -181,6 +181,6 @@ class BaseEMProblem(Problem.BaseProblem): # TODO: This should take a vector def MfRhoIDeriv(self,u): """ - Derivative of :code:`MfRhoI` with respect to the model. + Derivative of :code:`MfRhoI` with respect to the model. """ return self.mesh.getFaceInnerProductDeriv(self.curModel.rho, invMat=True)(u) * self.curModel.rhoDeriv diff --git a/SimPEG/EM/FDEM/SrcFDEM.py b/SimPEG/EM/FDEM/SrcFDEM.py index 8ec9c318..5ed531d1 100644 --- a/SimPEG/EM/FDEM/SrcFDEM.py +++ b/SimPEG/EM/FDEM/SrcFDEM.py @@ -608,46 +608,58 @@ class CircularLoop(BaseSrc): class PrimSecSigma(BaseSrc): # TODO: This will only work for E-B formulation - def __init__(self, rxList, freq, sigmaPrimary, ePrimary=None, primaryProblem=None, primarySurvey=None): + def __init__(self, rxList, freq, mPrimary, primaryProblem, primarySurvey): self.freq = float(freq) - self._ePrimary = ePrimary - self.sigmaPrimary = sigmaPrimary - self._primaryProblem = primaryProblem - self._primarySurvey = primarySurvey + self.mPrimary = mPrimary + self.primaryProblem = primaryProblem + self.primarySurvey = primarySurvey + self.primaryMesh = self.primaryProblem.mesh self.integrate = False BaseSrc.__init__(self, rxList) - def _solvePrimary(self): - # check if I have fields - # if not, solve the primary to get fields object - if self._primarySurvey.ispaired: - if self._primarySurvey.prob is not self._primaryProblem: - raise Exception "The survey object is already paired to a problem. Use survey.unpair()" - else: - self._primaryProblem.pair(self._primarySurvey) + @property + def MeSigmaPrimary(self, prob): + if getattr(self, '_MeSigmaPrimary', None) is None: + sigmaprimary = self.primaryProblem.mapping.sigmaMap * self.mPrimary + # if self.primaryMesh != prob.mesh: + # if self.mesh._meshType == 'CYL' and prob.mesh._meshType != 'CYL': + # P = self.mesh.getInterpolationMatCartMesh(prob.mesh,'CC') + # sigmaprimary = + self._MeSigmaPrimary = prob.mesh.getEdgeInnerProduct(sigmaprimary) + return self._MeSigmaPrimary - return + @property + def _primaryFields(self): + if getattr(self, '__primaryFields', None) is None: + # check if I have fields + # if not, solve the primary to get fields object + if self.primarySurvey.ispaired: + if self.primarySurvey.prob is not self.primaryProblem: + raise Exception('The survey object is already paired to a problem. Use survey.unpair()') + else: + self.primaryProblem.pair(self.primarySurvey) + self.__primaryFields = self.primaryProblem.fields(self.mPrimary) + return self.__primaryFields def ePrimary(self,prob): # check if a primary problem is defined if getattr(self, '_ePrimary', None) is None: - if self._primaryProblem is None or self._primarySurvey is None: - raise Exception "if Not specifying a ePrimary, a primarySurvey and primaryProblem must be provided." + if self.primaryProblem is None or self.primarySurvey is None: + raise Exception('if Not specifying a ePrimary, a primarySurvey and primaryProblem must be provided.') - self._ePrimary = self._solvePrimary[:,'e'] - - return self._ePrimary + return Utils.mkvc(self._primaryFields[:,'e']) def S_e(self,prob): + # todo --> see if the primary and secondary are on the same mesh or not MeSigma = prob.MeSigma - MeSigmaPrimary = prob.mesh.getEdgeInnerProduct(self.sigmaPrimary) - return (MeSigma - MeSigmaPrimary) * self._ePrimary + print type(prob.mesh) + return Utils.mkvc((MeSigma - self.MeSigmaPrimary(prob)) * self.ePrimary(prob)) def S_eDeriv(self, prob, v, adjoint = False): MeSigmaDeriv = prob.MeSigmaDeriv if adjoint is not True: - return MeSigmaDeriv(Utils.mkvc(self._ePrimary)) * v - return MeSigmaDeriv(Utils.mkvc(self._ePrimary)) * v + return MeSigmaDeriv(self.ePrimary(prob)) * v + return MeSigmaDeriv(self.ePrimary(prob)) * v diff --git a/SimPEG/EM/Utils/testingUtils.py b/SimPEG/EM/Utils/testingUtils.py index dd711c29..7b4ee2fd 100644 --- a/SimPEG/EM/Utils/testingUtils.py +++ b/SimPEG/EM/Utils/testingUtils.py @@ -20,7 +20,7 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, useMu=False, verbose=False): mesh = Mesh.TensorMesh([hx,hy,hz],['C','C','C']) if useMu is True: - mapping = [('sigma', Maps.ExpMap(mesh)), ('mu', Maps.IdentityMap(mesh))] + mapping = [('sigma', Maps.ExpMap(mesh)), ('mu', Maps.IdentityMap(mesh))] else: mapping = Maps.ExpMap(mesh) @@ -53,11 +53,15 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, useMu=False, verbose=False): Src.append(EM.FDEM.Src.RawVec([Rx0], freq, S_m, S_e)) elif SrcType is 'PrimSecSigma': if fdemType is 'e' or fdemType is 'b': - S_m = np.zeros(mesh.nF) - S_e = np.zeros(mesh.nE) - S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1. - S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1. - Src.append(EM.FDEM.Src.RawVec([Rx0], freq, S_m, S_e)) + # S_m = np.zeros(mesh.nF) + # S_e = np.zeros(mesh.nE) + # S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1. + # S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1. + primSrc = EM.FDEM.Src.MagDipole([], freq, np.r_[0.,0.,0.]) + primarySurvey = EM.FDEM.Survey([primSrc]) + primaryProblem = EM.FDEM.Problem_e(mesh,mapping=mapping) + mPrimary = np.ones(mapping.nP)*np.log(CONDUCTIVITY) + Src.append(EM.FDEM.Src.PrimSecSigma([Rx0], freq, mPrimary, primaryProblem, primarySurvey)) if verbose: print ' Fetching %s problem' % (fdemType) @@ -97,7 +101,7 @@ def crossCheckTest(SrcList, fdemType1, fdemType2, comp, addrandoms = False, useM prb1 = getFDEMProblem(fdemType1, comp, SrcList, freq, useMu, verbose) mesh = prb1.mesh print 'Cross Checking Forward: %s, %s formulations - %s' % (fdemType1, fdemType2, comp) - + logsig = np.log(np.ones(mesh.nC)*CONDUCTIVITY) mu = np.ones(mesh.nC)*MU diff --git a/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py b/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py index 0a2e8b82..fdaabe86 100644 --- a/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py +++ b/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py @@ -8,8 +8,8 @@ from SimPEG.EM.Utils.testingUtils import getFDEMProblem testE = True testB = True -testH = True -testJ = True +testH = False +testJ = False verbose = False @@ -20,8 +20,8 @@ MU = mu_0 freq = 1e-1 addrandoms = True -SrcType = ['MagDipole', 'RawVec'] #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVec' - +# SrcType = ['MagDipole', 'RawVec'] #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVec' +SrcType = ['PrimSecSigma'] def derivTest(fdemType, comp):