diff --git a/SimPEG/EM/FDEM/FieldsFDEM.py b/SimPEG/EM/FDEM/FieldsFDEM.py index 8f6fafe9..d3d7f3d5 100644 --- a/SimPEG/EM/FDEM/FieldsFDEM.py +++ b/SimPEG/EM/FDEM/FieldsFDEM.py @@ -3,7 +3,7 @@ import scipy.sparse as sp import SimPEG from SimPEG import Utils from SimPEG.EM.Utils import omega -from SimPEG.Utils import Zero, Identity +from SimPEG.Utils import Zero, Identity, sdiag class Fields(SimPEG.Problem.Fields): @@ -19,7 +19,9 @@ class Fields_e(Fields): 'eSecondary' : ['eSolution','E','_eSecondary'], 'b' : ['eSolution','F','_b'], 'bPrimary' : ['eSolution','F','_bPrimary'], - 'bSecondary' : ['eSolution','F','_bSecondary'] + 'bSecondary' : ['eSolution','F','_bSecondary'], + 'j' : ['eSolution','CC','_j'], + 'h' : ['eSolution','CC','_h'], } def __init__(self,mesh,survey,**kwargs): @@ -28,6 +30,21 @@ class Fields_e(Fields): def startup(self): self.prob = self.survey.prob self._edgeCurl = self.survey.prob.mesh.edgeCurl + self._aveE2CCV = self.survey.prob.mesh.aveE2CCV + self._aveF2CCV = self.survey.prob.mesh.aveF2CCV + self._sigma = self.survey.prob.curModel.sigma + self._sigmaDeriv = self.survey.prob.curModel.sigmaDeriv + self._nC = self.survey.prob.mesh.nC + + def _GLoc(self,fieldType): + if fieldType == 'e': + return 'E' + elif fieldType == 'b': + return 'F' + elif (fieldType == 'h') or (fieldType == 'j'): + return 'CC' + else: + raise Exception('Field type must be e, b, h, j') def _ePrimary(self, eSolution, srcList): ePrimary = np.zeros_like(eSolution) @@ -87,6 +104,41 @@ class Fields_e(Fields): # Assuming the primary does not depend on the model return self._bSecondaryDeriv_m(src, v, adjoint) + def _j(self, eSolution, srcList): + sigma = self._sigma + aveE2CCV = self._aveE2CCV + n = int(aveE2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy + VI = sdiag(1./np.kron(np.ones(n), self.prob.mesh.vol)) + Sigma = sdiag(np.kron(np.ones(n), sigma)) + + e = self._e(eSolution, srcList) + + return Sigma * (aveE2CCV * e) + + def _jDeriv_u(self, src, v, adjoint=False): + raise NotImplementedError + sigma = self._sigma + aveE2CCV = self._aveE2CCV + n = int(aveE2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy + Sigma = sdiag(sp.kron(np.ones(n), sigma)) + + if not adjoint: + return Sigma * (aveE2CCV * (v + self._eDeriv_u(src, v, adjoint))) + return aveE2CCV.T * Sigma.T * v + + def _jDeriv_m(self, src, v, adjoint=False): + raise NotImplementedError + sigma = self._sigma + aveE2CCV = self._aveE2CCV + n = int(aveE2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy + Sigma = sdiag(sp.kron(np.ones(n), sigma)) + + if not adjoint: + dsigma_dm = self._sigmaDeriv(v) + dSigma_dm = sdiag(sp.kron(np.ones(n), dsigma_dm)) + + + class Fields_b(Fields): knownFields = {'bSolution':'F'} @@ -110,6 +162,16 @@ class Fields_b(Fields): self._MeSigmaIDeriv = self.survey.prob.MeSigmaIDeriv self._Me = self.survey.prob.Me + def _GLoc(self,fieldType): + if fieldType == 'e': + return 'E' + elif fieldType == 'b': + return 'F' + elif (fieldType == 'h') or (fieldType == 'j'): + return'CC' + else: + raise Exception('Field type must be e, b, h, j') + def _bPrimary(self, bSolution, srcList): bPrimary = np.zeros_like(bSolution) for i, src in enumerate(srcList): @@ -193,6 +255,7 @@ class Fields_j(Fields): 'h' : ['jSolution','E','_h'], 'hPrimary' : ['jSolution','E','_hPrimary'], 'hSecondary' : ['jSolution','E','_hSecondary'], + 'e' : ['jSolution','C','_e'], } def __init__(self,mesh,survey,**kwargs): @@ -205,6 +268,19 @@ class Fields_j(Fields): self._MfRho = self.survey.prob.MfRho self._MfRhoDeriv = self.survey.prob.MfRhoDeriv self._Me = self.survey.prob.Me + self._rho = self.survey.prob.curModel.rho + self._aveF2CCV = self.survey.prob.mesh.aveF2CCV + self._nC = self.survey.prob.mesh.nC + + def _GLoc(self,fieldType): + if fieldType == 'h': + return 'E' + elif fieldType == 'j': + return 'F' + elif (fieldType == 'e') or (fieldType == 'b'): + return 'CC' + else: + raise Exception('Field type must be e, b, h, j') def _jPrimary(self, jSolution, srcList): jPrimary = np.zeros_like(jSolution,dtype = complex) @@ -281,6 +357,22 @@ class Fields_j(Fields): # assuming the primary doesn't depend on the model return self._hSecondaryDeriv_m(src, v, adjoint) + def _e(self, jSolution, srcList): + rho = self._rho + aveF2CCV = self._aveF2CCV + n = int(aveF2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy + Rho = sdiag(np.kron(np.ones(n), rho)) + + j = self._j(jSolution, srcList) + + return Rho * (aveF2CCV * j) + + def _eDeriv_u(self, src, v, adjoint=False): + raise NotImplementedError + + def _eDeriv_m(self, src, v, adjoint=False): + raise NotImplementedError + class Fields_h(Fields): knownFields = {'hSolution':'E'} diff --git a/SimPEG/EM/FDEM/SurveyFDEM.py b/SimPEG/EM/FDEM/SurveyFDEM.py index f60cbfdf..36be9ba5 100644 --- a/SimPEG/EM/FDEM/SurveyFDEM.py +++ b/SimPEG/EM/FDEM/SurveyFDEM.py @@ -3,6 +3,7 @@ from SimPEG.EM.Utils import * from scipy.constants import mu_0 from SimPEG.Utils import Zero, Identity import SrcFDEM as Src +from SimPEG import sp #################################################### @@ -12,33 +13,33 @@ import SrcFDEM as Src class Rx(SimPEG.Survey.BaseRx): knownRxTypes = { - 'exr':['e', 'Ex', 'real'], - 'eyr':['e', 'Ey', 'real'], - 'ezr':['e', 'Ez', 'real'], - 'exi':['e', 'Ex', 'imag'], - 'eyi':['e', 'Ey', 'imag'], - 'ezi':['e', 'Ez', 'imag'], + 'exr':['e', 'x', 'real'], + 'eyr':['e', 'y', 'real'], + 'ezr':['e', 'z', 'real'], + 'exi':['e', 'x', 'imag'], + 'eyi':['e', 'y', 'imag'], + 'ezi':['e', 'z', 'imag'], - 'bxr':['b', 'Fx', 'real'], - 'byr':['b', 'Fy', 'real'], - 'bzr':['b', 'Fz', 'real'], - 'bxi':['b', 'Fx', 'imag'], - 'byi':['b', 'Fy', 'imag'], - 'bzi':['b', 'Fz', 'imag'], + 'bxr':['b', 'x', 'real'], + 'byr':['b', 'y', 'real'], + 'bzr':['b', 'z', 'real'], + 'bxi':['b', 'x', 'imag'], + 'byi':['b', 'y', 'imag'], + 'bzi':['b', 'z', 'imag'], - 'jxr':['j', 'Fx', 'real'], - 'jyr':['j', 'Fy', 'real'], - 'jzr':['j', 'Fz', 'real'], - 'jxi':['j', 'Fx', 'imag'], - 'jyi':['j', 'Fy', 'imag'], - 'jzi':['j', 'Fz', 'imag'], + 'jxr':['j', 'x', 'real'], + 'jyr':['j', 'y', 'real'], + 'jzr':['j', 'z', 'real'], + 'jxi':['j', 'x', 'imag'], + 'jyi':['j', 'y', 'imag'], + 'jzi':['j', 'z', 'imag'], - 'hxr':['h', 'Ex', 'real'], - 'hyr':['h', 'Ey', 'real'], - 'hzr':['h', 'Ez', 'real'], - 'hxi':['h', 'Ex', 'imag'], - 'hyi':['h', 'Ey', 'imag'], - 'hzi':['h', 'Ez', 'imag'], + 'hxr':['h', 'x', 'real'], + 'hyr':['h', 'y', 'real'], + 'hzr':['h', 'z', 'real'], + 'hxi':['h', 'x', 'imag'], + 'hyi':['h', 'y', 'imag'], + 'hzi':['h', 'z', 'imag'], } radius = None @@ -50,10 +51,11 @@ class Rx(SimPEG.Survey.BaseRx): """Field Type projection (e.g. e b ...)""" return self.knownRxTypes[self.rxType][0] - @property - def projGLoc(self): - """Grid Location projection (e.g. Ex Fy ...)""" - return self.knownRxTypes[self.rxType][1] + # @property + # def projGLoc(self, u): + # """Grid Location projection (e.g. Ex Fy ...)""" + # return u._GLoc(self.rxType[0]) + # return self.knownRxTypes[self.rxType][1] @property def projComp(self): @@ -61,15 +63,46 @@ class Rx(SimPEG.Survey.BaseRx): return self.knownRxTypes[self.rxType][2] def projectFields(self, src, mesh, u): - P = self.getP(mesh) + u_part_complex = u[src, self.projField] # get the real or imag component real_or_imag = self.projComp u_part = getattr(u_part_complex, real_or_imag) + + projGLoc = u._GLoc(self.knownRxTypes[self.rxType][0]) + + if projGLoc == 'CC': + P = self.getP(mesh, projGLoc) + Z = 0.*P + if mesh.dim == 3: + if mesh._meshType == 'CYL' and mesh.isSymmetric and u_part.size > mesh.nC: # TODO: there must be a better way to do this! + if self.knownRxTypes[self.rxType][1] == 'x': + P = sp.hstack([P,Z]) + elif self.knownRxTypes[self.rxType][1] == 'z': + P = sp.hstack([Z,P]) + elif self.knownRxTypes[self.rxType][1] == 'y': + raise Exception('Symmetric CylMesh does not support y interpolation, as this variable does not exist.') + else: + if self.knownRxTypes[self.rxType][1] == 'x': + P = sp.hstack([P,Z,Z]) + elif self.knownRxTypes[self.rxType][1] == 'y': + P = sp.hstack([Z,P,Z]) + elif self.knownRxTypes[self.rxType][1] == 'z': + P = sp.hstack([Z,Z,P]) + else: + projGLoc += self.knownRxTypes[self.rxType][1] + P = self.getP(mesh, projGLoc) + return P*u_part def projectFieldsDeriv(self, src, mesh, u, v, adjoint=False): - P = self.getP(mesh) + + projGLoc = u._GLoc(self.knownRxTypes[self.rxType][0]) + + if projGLoc != 'CC': + projGLoc += self.knownRxTypes[self.rxType][1] + + P = self.getP(mesh, projGLoc) if not adjoint: Pv_complex = P * v diff --git a/SimPEG/EM/Utils/testingUtils.py b/SimPEG/EM/Utils/testingUtils.py index 8c703083..cc630d94 100644 --- a/SimPEG/EM/Utils/testingUtils.py +++ b/SimPEG/EM/Utils/testingUtils.py @@ -5,9 +5,9 @@ import sys from scipy.constants import mu_0 def getFDEMProblem(fdemType, comp, SrcList, freq, verbose=False): - cs = 5. + cs = 10. ncx, ncy, ncz = 6, 6, 6 - npad = 3 + npad = 5 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)] diff --git a/SimPEG/Mesh/TensorMesh.py b/SimPEG/Mesh/TensorMesh.py index c76306fe..874a1f9f 100644 --- a/SimPEG/Mesh/TensorMesh.py +++ b/SimPEG/Mesh/TensorMesh.py @@ -233,6 +233,10 @@ class BaseTensorMesh(BaseMesh): 'Fz' -> z-component of field defined on faces 'N' -> scalar field defined on nodes 'CC' -> scalar field defined on cell centers + # 'CCVx' -> x-component of a field defined on cell centers + # 'CCVy' -> y-component of a field defined on cell centers + # 'CCVz' -> z-component of a field defined on cell centers + """ if self._meshType == 'CYL' and self.isSymmetric and locType in ['Ex','Ez','Fy']: raise Exception('Symmetric CylMesh does not support %s interpolation, as this variable does not exist.' % locType) @@ -256,6 +260,12 @@ class BaseTensorMesh(BaseMesh): Q = sp.hstack(components) elif locType in ['CC', 'N']: Q = Utils.interpmat(loc, *self.getTensor(locType)) + # elif locType in ['CCVx', 'CCVy', 'CCVz']: + # Q = Utils.interpmat(loc, 'CC') + # Zero = 0.*Q + + + else: raise NotImplementedError('getInterpolationMat: locType=='+locType+' and mesh.dim=='+str(self.dim)) diff --git a/SimPEG/Survey.py b/SimPEG/Survey.py index 88355df1..957b2c73 100644 --- a/SimPEG/Survey.py +++ b/SimPEG/Survey.py @@ -35,7 +35,7 @@ class BaseRx(object): """Number of data in the receiver.""" return self.locs.shape[0] - def getP(self, mesh): + def getP(self, mesh, projGLoc=None): """ Returns the projection matrices as a list for all components collected by @@ -48,7 +48,10 @@ class BaseRx(object): if mesh in self._Ps: return self._Ps[mesh] - P = mesh.getInterpolationMat(self.locs, self.projGLoc) + if projGLoc is None: + projGLoc = self.projGLoc + + P = mesh.getInterpolationMat(self.locs, projGLoc) if self.storeProjections: self._Ps[mesh] = P return P diff --git a/tests/em/fdem/forward/test_FDEM_forward.py b/tests/em/fdem/forward/test_FDEM_forward.py index 437f3708..cf1ad608 100644 --- a/tests/em/fdem/forward/test_FDEM_forward.py +++ b/tests/em/fdem/forward/test_FDEM_forward.py @@ -7,26 +7,28 @@ from SimPEG.EM.Utils.testingUtils import getFDEMProblem testEB = True testHJ = True - +testEJ = False verbose = False -TOL = 1e-5 +TOLEBHJ = 1e-5 +TOLEJHB = 1e-1 + FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order CONDUCTIVITY = 1e1 MU = mu_0 freq = 1e-1 -addrandoms = True +addrandoms = False SrcList = ['RawVec', 'MagDipole_Bfield', 'MagDipole', 'CircularLoop'] -def crossCheckTest(fdemType, comp): +def crossCheckTest(fdemType1, fdemType2, comp, TOL=TOLEBHJ): l2norm = lambda r: np.sqrt(r.dot(r)) - prb1 = getFDEMProblem(fdemType, comp, SrcList, freq, verbose) + prb1 = getFDEMProblem(fdemType1, comp, SrcList, freq, verbose) mesh = prb1.mesh - print 'Cross Checking Forward: %s formulation - %s' % (fdemType, comp) + print 'Cross Checking Forward: %s formulation - %s' % (fdemType1, comp) m = np.log(np.ones(mesh.nC)*CONDUCTIVITY) mu = np.log(np.ones(mesh.nC)*MU) @@ -42,16 +44,8 @@ def crossCheckTest(fdemType, comp): if verbose: print ' Problem 1 solved' - if fdemType == 'e': - prb2 = getFDEMProblem('b', comp, SrcList, freq, verbose) - elif fdemType == 'b': - prb2 = getFDEMProblem('e', comp, SrcList, freq, verbose) - elif fdemType == 'j': - prb2 = getFDEMProblem('h', comp, SrcList, freq, verbose) - elif fdemType == 'h': - prb2 = getFDEMProblem('j', comp, SrcList, freq, verbose) - else: - raise NotImplementedError() + + prb2 = getFDEMProblem(fdemType2, comp, SrcList, freq, verbose) # prb2.mu = mu survey2 = prb2.survey @@ -71,57 +65,84 @@ def crossCheckTest(fdemType, comp): class FDEM_CrossCheck(unittest.TestCase): if testEB: def test_EB_CrossCheck_exr_Eform(self): - self.assertTrue(crossCheckTest('e', 'exr')) + self.assertTrue(crossCheckTest('e', 'b', 'exr')) def test_EB_CrossCheck_eyr_Eform(self): - self.assertTrue(crossCheckTest('e', 'eyr')) + self.assertTrue(crossCheckTest('e', 'b', 'eyr')) def test_EB_CrossCheck_ezr_Eform(self): - self.assertTrue(crossCheckTest('e', 'ezr')) + self.assertTrue(crossCheckTest('e', 'b', 'ezr')) def test_EB_CrossCheck_exi_Eform(self): - self.assertTrue(crossCheckTest('e', 'exi')) + self.assertTrue(crossCheckTest('e', 'b', 'exi')) def test_EB_CrossCheck_eyi_Eform(self): - self.assertTrue(crossCheckTest('e', 'eyi')) + self.assertTrue(crossCheckTest('e', 'b', 'eyi')) def test_EB_CrossCheck_ezi_Eform(self): - self.assertTrue(crossCheckTest('e', 'ezi')) + self.assertTrue(crossCheckTest('e', 'b', 'ezi')) def test_EB_CrossCheck_bxr_Eform(self): - self.assertTrue(crossCheckTest('e', 'bxr')) + self.assertTrue(crossCheckTest('e', 'b', 'bxr')) def test_EB_CrossCheck_byr_Eform(self): - self.assertTrue(crossCheckTest('e', 'byr')) + self.assertTrue(crossCheckTest('e', 'b', 'byr')) def test_EB_CrossCheck_bzr_Eform(self): - self.assertTrue(crossCheckTest('e', 'bzr')) + self.assertTrue(crossCheckTest('e', 'b', 'bzr')) def test_EB_CrossCheck_bxi_Eform(self): - self.assertTrue(crossCheckTest('e', 'bxi')) + self.assertTrue(crossCheckTest('e', 'b', 'bxi')) def test_EB_CrossCheck_byi_Eform(self): - self.assertTrue(crossCheckTest('e', 'byi')) + self.assertTrue(crossCheckTest('e', 'b', 'byi')) def test_EB_CrossCheck_bzi_Eform(self): - self.assertTrue(crossCheckTest('e', 'bzi')) + self.assertTrue(crossCheckTest('e', 'b', 'bzi')) if testHJ: def test_HJ_CrossCheck_jxr_Jform(self): - self.assertTrue(crossCheckTest('j', 'jxr')) + self.assertTrue(crossCheckTest('j', 'h', 'jxr')) def test_HJ_CrossCheck_jyr_Jform(self): - self.assertTrue(crossCheckTest('j', 'jyr')) + self.assertTrue(crossCheckTest('j', 'h', 'jyr')) def test_HJ_CrossCheck_jzr_Jform(self): - self.assertTrue(crossCheckTest('j', 'jzr')) + self.assertTrue(crossCheckTest('j', 'h', 'jzr')) def test_HJ_CrossCheck_jxi_Jform(self): - self.assertTrue(crossCheckTest('j', 'jxi')) + self.assertTrue(crossCheckTest('j', 'h', 'jxi')) def test_HJ_CrossCheck_jyi_Jform(self): - self.assertTrue(crossCheckTest('j', 'jyi')) + self.assertTrue(crossCheckTest('j', 'h', 'jyi')) def test_HJ_CrossCheck_jzi_Jform(self): - self.assertTrue(crossCheckTest('j', 'jzi')) + self.assertTrue(crossCheckTest('j', 'h', 'jzi')) def test_HJ_CrossCheck_hxr_Jform(self): - self.assertTrue(crossCheckTest('j', 'hxr')) + self.assertTrue(crossCheckTest('j', 'h', 'hxr')) def test_HJ_CrossCheck_hyr_Jform(self): - self.assertTrue(crossCheckTest('j', 'hyr')) + self.assertTrue(crossCheckTest('j', 'h', 'hyr')) def test_HJ_CrossCheck_hzr_Jform(self): - self.assertTrue(crossCheckTest('j', 'hzr')) + self.assertTrue(crossCheckTest('j', 'h', 'hzr')) def test_HJ_CrossCheck_hxi_Jform(self): - self.assertTrue(crossCheckTest('j', 'hxi')) + self.assertTrue(crossCheckTest('j', 'h', 'hxi')) def test_HJ_CrossCheck_hyi_Jform(self): - self.assertTrue(crossCheckTest('j', 'hyi')) + self.assertTrue(crossCheckTest('j', 'h', 'hyi')) def test_HJ_CrossCheck_hzi_Jform(self): - self.assertTrue(crossCheckTest('j', 'hzi')) + self.assertTrue(crossCheckTest('j', 'h', 'hzi')) + + if testEJ: + # def test_EJ_CrossCheck_jxr_Jform(self): + # self.assertTrue(crossCheckTest('e', 'j', 'jxr')) + # def test_EJ_CrossCheck_jyr_Jform(self): + # self.assertTrue(crossCheckTest('e', 'j', 'jyr')) + # def test_EJ_CrossCheck_jzr_Jform(self): + # self.assertTrue(crossCheckTest('e', 'j', 'jzr')) + # def test_EJ_CrossCheck_jxi_Jform(self): + # self.assertTrue(crossCheckTest('e', 'j', 'jxi')) + # def test_EJ_CrossCheck_jyi_Jform(self): + # self.assertTrue(crossCheckTest('e', 'j', 'jyi')) + # def test_EJ_CrossCheck_jzi_Jform(self): + # self.assertTrue(crossCheckTest('e', 'j', 'jzi')) + + def test_EJ_CrossCheck_jxr_Jform(self): + self.assertTrue(crossCheckTest('e', 'j', 'exr', TOL=TOLEJHB)) + # def test_EJ_CrossCheck_jyr_Jform(self): + # self.assertTrue(crossCheckTest('e', 'j', 'eyr')) + # def test_EJ_CrossCheck_jzr_Jform(self): + # self.assertTrue(crossCheckTest('e', 'j', 'ezr')) + # def test_EJ_CrossCheck_jxi_Jform(self): + # self.assertTrue(crossCheckTest('e', 'j', 'exi')) + # def test_EJ_CrossCheck_jyi_Jform(self): + # self.assertTrue(crossCheckTest('e', 'j', 'eyi')) + # def test_EJ_CrossCheck_jzi_Jform(self): + # self.assertTrue(crossCheckTest('e', 'j', 'ezi')) if __name__ == '__main__': unittest.main() \ No newline at end of file