From a478b976bcb4e3a382514b4a7bbb4b8db7b4457e Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Thu, 31 Mar 2016 13:27:29 -0700 Subject: [PATCH] prim sec for different meshes with EB formulation --- SimPEG/EM/FDEM/SrcFDEM.py | 89 +++++++++++++------ SimPEG/EM/Utils/testingUtils.py | 19 ++-- .../fdem/inverse/derivs/test_FDEM_derivs.py | 2 +- tests/mesh/test_InterpolationMesh2Mesh.py | 5 +- 4 files changed, 77 insertions(+), 38 deletions(-) diff --git a/SimPEG/EM/FDEM/SrcFDEM.py b/SimPEG/EM/FDEM/SrcFDEM.py index d361bd18..e8622eb2 100644 --- a/SimPEG/EM/FDEM/SrcFDEM.py +++ b/SimPEG/EM/FDEM/SrcFDEM.py @@ -608,54 +608,85 @@ class CircularLoop(BaseSrc): class PrimSecSigma(BaseSrc): # TODO: This will only work for E-B formulation - def __init__(self, rxList, freq, mPrimary, primaryProblem, primarySurvey): + def __init__(self, rxList, freq, mPrimary, primaryProblem=None, primarySurvey=None, primaryFields=None): self.freq = float(freq) self.mPrimary = mPrimary - self.primaryProblem = primaryProblem - self.primarySurvey = primarySurvey - self.primaryMesh = self.primaryProblem.mesh + + if primaryFields is None: + assert primaryProblem is not None and primarySurvey is not None, 'If no primary fields are provided, a primaryProblem and primarySurvey must be provided' + else: + self._fields = primaryFields + + if primaryProblem is not None and primarySurvey is not None: + self.prob = primaryProblem + self.survey = primarySurvey + if self.survey.ispaired: + if self.survey.prob is not self.prob: + raise Exception('The survey object is already paired to a problem. Use survey.unpair()') + else: + self.prob.pair(self.survey) + + + self.mesh = self.prob.mesh self.integrate = False BaseSrc.__init__(self, rxList) - @property + # @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 = + sigmaprimary = self.prob.mapping.sigmaMap * self.mPrimary + if self.mesh != prob.mesh: + P = self.mesh.getInterpolationMatMesh2Mesh(prob.mesh, locType='CC') + sigmaprimary = P * sigmaprimary self._MeSigmaPrimary = prob.mesh.getEdgeInnerProduct(sigmaprimary) return self._MeSigmaPrimary + # @property + def MfRhoIPrimary(self, prob): + if getattr(self, '_MfRhoIPrimary', None) is None: + rhoprimary = self.prob.mapping.rhoMap * self.mPrimary + if self.mesh != prob.mesh: + P = self.mesh.getInterpolationMatMesh2Mesh(prob.mesh, locType='CC') + rhoprimary = P * rhoprimary + self._MfRhoIPrimary = prob.mesh.getFaceInnerProduct(rhoprimary, invMat=True) + return self._MfRhoIPrimary + @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 fields(self): + if getattr(self, '_fields', None) is None: + # check if I have fields if not, solve the primary to get fields object + self._fields = self.prob.fields(self.mPrimary) + return self._fields 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.') + ePrimary = self.fields[:,'e'] - return Utils.mkvc(self._primaryFields[:,'e']) + if self.mesh != prob.mesh: + P = self.mesh.getInterpolationMatMesh2Mesh(prob.mesh, locType='E') + ePrimary = Utils.mkvc(P * ePrimary) + self._ePrimary = Utils.mkvc(ePrimary) - def S_e(self,prob): + return self._ePrimary + + def jPrimary(self,prob): + # check if a primary problem is defined + # if getattr(self, '_jPrimary', None) is None: + # if self.prob is None or self.survey is None: + # raise Exception('if Not specifying a jPrimary, a primarySurvey and primaryProblem must be provided.') + jPrimary = self.prob.fields(self.mPrimary)[:,'j'] + if self.mesh != prob.mesh: + P = self.mesh.getInterpolationMatMesh2Mesh(prob.mesh, locType='F') + jPrimary = P * jPrimary + return jPrimary + + def s_e(self,prob): # todo --> see if the primary and secondary are on the same mesh or not - MeSigma = prob.MeSigma - print type(prob.mesh) - return Utils.mkvc((MeSigma - self.MeSigmaPrimary(prob)) * self.ePrimary(prob)) + if prob._formulation == 'EB': + return Utils.mkvc(prob.MeSigma * self.ePrimary(prob) - self.MeSigmaPrimary(prob) * self.ePrimary(prob)) - def S_eDeriv(self, prob, v, adjoint = False): + def s_eDeriv(self, prob, v, adjoint = False): MeSigmaDeriv = prob.MeSigmaDeriv if adjoint is not True: return MeSigmaDeriv(self.ePrimary(prob)) * v diff --git a/SimPEG/EM/Utils/testingUtils.py b/SimPEG/EM/Utils/testingUtils.py index 7b4ee2fd..b913c558 100644 --- a/SimPEG/EM/Utils/testingUtils.py +++ b/SimPEG/EM/Utils/testingUtils.py @@ -53,15 +53,24 @@ 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. 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)) + Src.append(EM.FDEM.Src.PrimSecSigma([Rx0], freq, mPrimary, primaryProblem=primaryProblem, primarySurvey=primarySurvey)) + + elif SrcType is 'PrimSecCyl': + if fdemType is 'e' or fdemType is 'b': + + hx = [(cs,ncx + 2), (cs,npad + 2,1.3)] + hz = [(cs,npad + 2 ,-1.3), (cs,ncz+2), (cs,npad+2,1.3)] + primmesh = Mesh.CylMesh([hx,1,hz], '00C') + + primSrc = EM.FDEM.Src.MagDipole([], freq, np.r_[0.,0.,0.]) + primarySurvey = EM.FDEM.Survey([primSrc]) + primaryProblem = EM.FDEM.Problem_e(primmesh) + mPrimary = np.ones(primmesh.nC)*CONDUCTIVITY + Src.append(EM.FDEM.Src.PrimSecSigma([Rx0], freq, mPrimary, primaryProblem=primaryProblem, primarySurvey=primarySurvey)) if verbose: print ' Fetching %s problem' % (fdemType) diff --git a/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py b/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py index fdaabe86..0a455515 100644 --- a/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py +++ b/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py @@ -21,7 +21,7 @@ freq = 1e-1 addrandoms = True # SrcType = ['MagDipole', 'RawVec'] #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVec' -SrcType = ['PrimSecSigma'] +SrcType = ['PrimSecCyl'] def derivTest(fdemType, comp): diff --git a/tests/mesh/test_InterpolationMesh2Mesh.py b/tests/mesh/test_InterpolationMesh2Mesh.py index 79d3e665..7b4ddf45 100644 --- a/tests/mesh/test_InterpolationMesh2Mesh.py +++ b/tests/mesh/test_InterpolationMesh2Mesh.py @@ -4,9 +4,9 @@ from SimPEG.Utils import mkvc from SimPEG import Mesh, Tests import unittest -test1D = False +test1D = True test2D = True -test3D = False +test3D = True call1 = lambda fun, xyz: fun(xyz) call2 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, -1]) @@ -203,7 +203,6 @@ if test2D: v = call2(funX, self.M.gridN) P = self.M.getInterpolationMatMesh2Mesh(mesh2, locType=self.type) - print P.shape, v.shape num = P*v return np.linalg.norm((num - ana), np.inf)