From 884d27b5418693b28020ccea32b2e8865fe8f8a9 Mon Sep 17 00:00:00 2001 From: Lindsey Date: Tue, 8 Dec 2015 18:11:01 -0800 Subject: [PATCH 01/17] start of getting any field from any formulation --- SimPEG/EM/FDEM/FieldsFDEM.py | 96 +++++++++++++++++++- SimPEG/EM/FDEM/SurveyFDEM.py | 93 +++++++++++++------ SimPEG/EM/Utils/testingUtils.py | 4 +- SimPEG/Mesh/TensorMesh.py | 10 ++ SimPEG/Survey.py | 7 +- tests/em/fdem/forward/test_FDEM_forward.py | 101 +++++++++++++-------- 6 files changed, 235 insertions(+), 76 deletions(-) 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 From adda7a43dc5aa1c8f381709f1c34d9ba294b373f Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 8 Dec 2015 19:14:28 -0800 Subject: [PATCH 02/17] E from J and J from E, and a test (need to choose better parameters) --- SimPEG/EM/FDEM/FieldsFDEM.py | 12 ++++- SimPEG/EM/Utils/testingUtils.py | 14 +++--- tests/em/fdem/forward/test_FDEM_forward.py | 55 +++++++++++----------- 3 files changed, 46 insertions(+), 35 deletions(-) diff --git a/SimPEG/EM/FDEM/FieldsFDEM.py b/SimPEG/EM/FDEM/FieldsFDEM.py index d3d7f3d5..22838e17 100644 --- a/SimPEG/EM/FDEM/FieldsFDEM.py +++ b/SimPEG/EM/FDEM/FieldsFDEM.py @@ -365,7 +365,7 @@ class Fields_j(Fields): j = self._j(jSolution, srcList) - return Rho * (aveF2CCV * j) + return 1./self.survey.prob.mesh.dim * Rho * (aveF2CCV * j) def _eDeriv_u(self, src, v, adjoint=False): raise NotImplementedError @@ -394,6 +394,16 @@ class Fields_h(Fields): self._MeMuI = self.survey.prob.MeMuI self._MfRho = self.survey.prob.MfRho + 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 _hPrimary(self, hSolution, srcList): hPrimary = np.zeros_like(hSolution,dtype = complex) for i, src in enumerate(srcList): diff --git a/SimPEG/EM/Utils/testingUtils.py b/SimPEG/EM/Utils/testingUtils.py index cc630d94..bd6e0459 100644 --- a/SimPEG/EM/Utils/testingUtils.py +++ b/SimPEG/EM/Utils/testingUtils.py @@ -6,17 +6,17 @@ from scipy.constants import mu_0 def getFDEMProblem(fdemType, comp, SrcList, freq, verbose=False): cs = 10. - ncx, ncy, ncz = 6, 6, 6 - 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)] + ncx, ncy, ncz = 2, 2, 2 + npad = 8 + hx = [(cs,npad,-1.5), (cs,ncx), (cs,npad,1.5)] + hy = [(cs,npad,-1.5), (cs,ncy), (cs,npad,1.5)] + hz = [(cs,npad,-1.5), (cs,ncz), (cs,npad,1.5)] 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.]) + x = np.array([np.linspace(-5.5*cs,-2.5*cs,4),np.linspace(3.5*cs,2.5*cs,4)]) #don't sample right by the source + XYZ = Utils.ndgrid(x,x,np.linspace(-2.25*cs,2.25*cs,5)) Rx0 = EM.FDEM.Rx(XYZ, comp) Src = [] diff --git a/tests/em/fdem/forward/test_FDEM_forward.py b/tests/em/fdem/forward/test_FDEM_forward.py index cf1ad608..21c0d114 100644 --- a/tests/em/fdem/forward/test_FDEM_forward.py +++ b/tests/em/fdem/forward/test_FDEM_forward.py @@ -7,16 +7,17 @@ from SimPEG.EM.Utils.testingUtils import getFDEMProblem testEB = True testHJ = True -testEJ = False +testEJ = True verbose = False TOLEBHJ = 1e-5 -TOLEJHB = 1e-1 +TOLEJHB = 1 # averaging and more sensitive to boundary condition violations (ie. the impact of violating the boundary conditions in each case is different.) +#TODO: choose better testing parameters to lower this FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order CONDUCTIVITY = 1e1 MU = mu_0 -freq = 1e-1 +freq = 5e-1 addrandoms = False SrcList = ['RawVec', 'MagDipole_Bfield', 'MagDipole', 'CircularLoop'] @@ -28,7 +29,7 @@ def crossCheckTest(fdemType1, fdemType2, comp, TOL=TOLEBHJ): prb1 = getFDEMProblem(fdemType1, comp, SrcList, freq, verbose) mesh = prb1.mesh - print 'Cross Checking Forward: %s formulation - %s' % (fdemType1, comp) + print 'Cross Checking Forward: %s, %s formulations - %s' % (fdemType1, fdemType2, comp) m = np.log(np.ones(mesh.nC)*CONDUCTIVITY) mu = np.log(np.ones(mesh.nC)*MU) @@ -118,31 +119,31 @@ class FDEM_CrossCheck(unittest.TestCase): 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', 'jxr', TOL=TOLEJHB)) + def test_EJ_CrossCheck_jyr_Jform(self): + self.assertTrue(crossCheckTest('e', 'j', 'jyr', TOL=TOLEJHB)) + def test_EJ_CrossCheck_jzr_Jform(self): + self.assertTrue(crossCheckTest('e', 'j', 'jzr', TOL=TOLEJHB)) + def test_EJ_CrossCheck_jxi_Jform(self): + self.assertTrue(crossCheckTest('e', 'j', 'jxi', TOL=TOLEJHB)) + def test_EJ_CrossCheck_jyi_Jform(self): + self.assertTrue(crossCheckTest('e', 'j', 'jyi', TOL=TOLEJHB)) + def test_EJ_CrossCheck_jzi_Jform(self): + self.assertTrue(crossCheckTest('e', 'j', 'jzi', TOL=TOLEJHB)) + + def test_EJ_CrossCheck_exr_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')) + def test_EJ_CrossCheck_eyr_Jform(self): + self.assertTrue(crossCheckTest('e', 'j', 'eyr', TOL=TOLEJHB)) + def test_EJ_CrossCheck_ezr_Jform(self): + self.assertTrue(crossCheckTest('e', 'j', 'ezr', TOL=TOLEJHB)) + def test_EJ_CrossCheck_exi_Jform(self): + self.assertTrue(crossCheckTest('e', 'j', 'exi', TOL=TOLEJHB)) + def test_EJ_CrossCheck_eyi_Jform(self): + self.assertTrue(crossCheckTest('e', 'j', 'eyi', TOL=TOLEJHB)) + def test_EJ_CrossCheck_ezi_Jform(self): + self.assertTrue(crossCheckTest('e', 'j', 'ezi', TOL=TOLEJHB)) if __name__ == '__main__': unittest.main() \ No newline at end of file From b620c4286b624721ae73a88ca961491d501c29c2 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 8 Dec 2015 21:48:01 -0800 Subject: [PATCH 03/17] smaller test mesh --- SimPEG/EM/Utils/testingUtils.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/SimPEG/EM/Utils/testingUtils.py b/SimPEG/EM/Utils/testingUtils.py index bd6e0459..59971192 100644 --- a/SimPEG/EM/Utils/testingUtils.py +++ b/SimPEG/EM/Utils/testingUtils.py @@ -6,17 +6,17 @@ from scipy.constants import mu_0 def getFDEMProblem(fdemType, comp, SrcList, freq, verbose=False): cs = 10. - ncx, ncy, ncz = 2, 2, 2 + ncx, ncy, ncz = 0, 0, 0 npad = 8 - hx = [(cs,npad,-1.5), (cs,ncx), (cs,npad,1.5)] - hy = [(cs,npad,-1.5), (cs,ncy), (cs,npad,1.5)] - hz = [(cs,npad,-1.5), (cs,ncz), (cs,npad,1.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)] mesh = Mesh.TensorMesh([hx,hy,hz],['C','C','C']) mapping = Maps.ExpMap(mesh) - x = np.array([np.linspace(-5.5*cs,-2.5*cs,4),np.linspace(3.5*cs,2.5*cs,4)]) #don't sample right by the source - XYZ = Utils.ndgrid(x,x,np.linspace(-2.25*cs,2.25*cs,5)) + x = np.array([np.linspace(-5.*cs,-2.*cs,3),np.linspace(5.*cs,2.*cs,3)]) #don't sample right by the source + XYZ = Utils.ndgrid(x,x,np.linspace(-2.*cs,2.*cs,5)) Rx0 = EM.FDEM.Rx(XYZ, comp) Src = [] @@ -70,6 +70,6 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, verbose=False): from pymatsolver import MumpsSolver prb.Solver = MumpsSolver except ImportError, e: - pass + prb.Solver = SolverLU return prb \ No newline at end of file From 9e69455d6e1116e67fe3ad889e1693effb8b02ac Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 8 Dec 2015 23:35:06 -0800 Subject: [PATCH 04/17] b from h and h from b --- SimPEG/EM/FDEM/FieldsFDEM.py | 108 ++++++++++++++++++++- SimPEG/EM/Utils/testingUtils.py | 8 +- tests/em/fdem/forward/test_FDEM_forward.py | 80 +++++++++++++++ 3 files changed, 187 insertions(+), 9 deletions(-) diff --git a/SimPEG/EM/FDEM/FieldsFDEM.py b/SimPEG/EM/FDEM/FieldsFDEM.py index 22838e17..fb05b8c6 100644 --- a/SimPEG/EM/FDEM/FieldsFDEM.py +++ b/SimPEG/EM/FDEM/FieldsFDEM.py @@ -34,6 +34,7 @@ class Fields_e(Fields): self._aveF2CCV = self.survey.prob.mesh.aveF2CCV self._sigma = self.survey.prob.curModel.sigma self._sigmaDeriv = self.survey.prob.curModel.sigmaDeriv + self._mui = self.survey.prob.curModel.mui self._nC = self.survey.prob.mesh.nC def _GLoc(self,fieldType): @@ -108,12 +109,25 @@ class Fields_e(Fields): sigma = self._sigma aveE2CCV = self._aveE2CCV n = int(aveE2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy + # Sigma = sdiag(np.kron(np.ones(n), sigma)) + Sigma = self.prob.MeSigma 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) + return VI * (aveE2CCV * (Sigma *e) ) + + def _h(self, eolution, srcList): + b = self._b(eSolution, srcList) + Mui = self.survey.prob.MfMui + aveF2CCV = self._aveF2CCV + n = int(aveF2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy + # Mui = sdiag(sp.kron(np.ones(n), mui)) + VI = sdiag(1./np.kron(np.ones(n), self.prob.mesh.vol)) + + + return VI * (aveF2CCV * (Mui * b)) + def _jDeriv_u(self, src, v, adjoint=False): raise NotImplementedError @@ -149,6 +163,8 @@ class Fields_b(Fields): 'e' : ['bSolution','E','_e'], 'ePrimary' : ['bSolution','E','_ePrimary'], 'eSecondary' : ['bSolution','E','_eSecondary'], + 'j' : ['bSolution','C','_j'], + 'h' : ['bSolution','C','_h'], } def __init__(self,mesh,survey,**kwargs): @@ -161,6 +177,13 @@ class Fields_b(Fields): self._MfMui = self.survey.prob.MfMui self._MeSigmaIDeriv = self.survey.prob.MeSigmaIDeriv self._Me = self.survey.prob.Me + self._aveF2CCV = self.survey.prob.mesh.aveF2CCV + self._aveE2CCV = self.survey.prob.mesh.aveE2CCV + self._sigma = self.survey.prob.curModel.sigma + self._mui = self.survey.prob.curModel.mui + self._nC = self.survey.prob.mesh.nC + + def _GLoc(self,fieldType): if fieldType == 'e': @@ -245,6 +268,29 @@ class Fields_b(Fields): # assuming primary doesn't depend on model return self._eSecondaryDeriv_m(src, v, adjoint) + def _j(self, bSolution, srcList): + sigma = self._sigma + aveE2CCV = self._aveE2CCV + n = int(aveE2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy + # Sigma = sdiag(np.kron(np.ones(n), sigma)) + Sigma = self.prob.MeSigma + VI = sdiag(1./np.kron(np.ones(n), self.prob.mesh.vol)) + + e = self._e(bSolution, srcList) + + return VI * (aveE2CCV * (Sigma *e) ) + + def _h(self, bSolution, srcList): + b = self._b(bSolution, srcList) + Mui = self.survey.prob.MfMui + aveF2CCV = self._aveF2CCV + n = int(aveF2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy + # Mui = sdiag(sp.kron(np.ones(n), mui)) + VI = sdiag(1./np.kron(np.ones(n), self.prob.mesh.vol)) + + + return VI * (aveF2CCV * (Mui * b)) + class Fields_j(Fields): knownFields = {'jSolution':'F'} @@ -256,6 +302,7 @@ class Fields_j(Fields): 'hPrimary' : ['jSolution','E','_hPrimary'], 'hSecondary' : ['jSolution','E','_hSecondary'], 'e' : ['jSolution','C','_e'], + 'b' : ['jSolution','C','_b'], } def __init__(self,mesh,survey,**kwargs): @@ -269,7 +316,9 @@ class Fields_j(Fields): self._MfRhoDeriv = self.survey.prob.MfRhoDeriv self._Me = self.survey.prob.Me self._rho = self.survey.prob.curModel.rho + self._mu = self.survey.prob.curModel.mui self._aveF2CCV = self.survey.prob.mesh.aveF2CCV + self._aveE2CCV = self.survey.prob.mesh.aveE2CCV self._nC = self.survey.prob.mesh.nC def _GLoc(self,fieldType): @@ -361,11 +410,14 @@ class Fields_j(Fields): 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)) + + Rho = self.prob.MfRho + VI = sdiag(1./np.kron(np.ones(n), self.prob.mesh.vol)) + j = self._j(jSolution, srcList) - return 1./self.survey.prob.mesh.dim * Rho * (aveF2CCV * j) + return VI * (aveF2CCV * (Rho * j)) def _eDeriv_u(self, src, v, adjoint=False): raise NotImplementedError @@ -373,6 +425,16 @@ class Fields_j(Fields): def _eDeriv_m(self, src, v, adjoint=False): raise NotImplementedError + def _b(self, jSolution, srcList): + h = self._h(jSolution, srcList) + Mu = self.prob.MeMu + aveE2CCV = self._aveE2CCV + n = int(aveE2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy + # Mu = sdiag(sp.kron(np.ones(n), mu)) + VI = sdiag(1./np.kron(np.ones(n), self.prob.mesh.vol)) + + return VI * (aveE2CCV * (Mu * h)) + class Fields_h(Fields): knownFields = {'hSolution':'E'} @@ -382,7 +444,9 @@ class Fields_h(Fields): 'hSecondary' : ['hSolution','E','_hSecondary'], 'j' : ['hSolution','F','_j'], 'jPrimary' : ['hSolution','F','_jPrimary'], - 'jSecondary' : ['hSolution','F','_jSecondary'] + 'jSecondary' : ['hSolution','F','_jSecondary'], + 'e' : ['hSolution','C','_e'], + 'b' : ['hSolution','C','_b'], } def __init__(self,mesh,survey,**kwargs): @@ -393,6 +457,11 @@ class Fields_h(Fields): self._edgeCurl = self.survey.prob.mesh.edgeCurl self._MeMuI = self.survey.prob.MeMuI self._MfRho = self.survey.prob.MfRho + self._rho = self.survey.prob.curModel.rho + self._mu = self.survey.prob.curModel.mui + self._aveF2CCV = self.survey.prob.mesh.aveF2CCV + self._aveE2CCV = self.survey.prob.mesh.aveE2CCV + self._nC = self.survey.prob.mesh.nC def _GLoc(self,fieldType): if fieldType == 'h': @@ -458,3 +527,32 @@ class Fields_h(Fields): def _jDeriv_m(self, src, v, adjoint=False): # assuming the primary does not depend on the model return self._jSecondaryDeriv_m(src,v,adjoint) + + def _e(self, hSolution, srcList): + rho = self._rho + aveF2CCV = self._aveF2CCV + n = int(aveF2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy + + Rho = self.prob.MfRho + VI = sdiag(1./np.kron(np.ones(n), self.prob.mesh.vol)) + + + j = self._j(hSolution, srcList) + + return VI * (aveF2CCV * (Rho * j)) + + def _eDeriv_u(self, src, v, adjoint=False): + raise NotImplementedError + + def _eDeriv_m(self, src, v, adjoint=False): + raise NotImplementedError + + def _b(self, hSolution, srcList): + h = self._h(hSolution, srcList) + Mu = self.prob.MeMu + aveE2CCV = self._aveE2CCV + n = int(aveE2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy + # Mu = sdiag(sp.kron(np.ones(n), mu)) + VI = sdiag(1./np.kron(np.ones(n), self.prob.mesh.vol)) + + return VI * (aveE2CCV * (Mu * h)) diff --git a/SimPEG/EM/Utils/testingUtils.py b/SimPEG/EM/Utils/testingUtils.py index 59971192..9979417d 100644 --- a/SimPEG/EM/Utils/testingUtils.py +++ b/SimPEG/EM/Utils/testingUtils.py @@ -32,15 +32,15 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, verbose=False): 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. + S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1e-3 + S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1e-3 Src.append(EM.FDEM.Src.RawVec([Rx0], freq, S_m, S_e)) elif fdemType is 'h' or fdemType is 'j': S_m = np.zeros(mesh.nE) S_e = np.zeros(mesh.nF) - S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1. - S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1. + S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1e-3 + S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1e-3 Src.append(EM.FDEM.Src.RawVec([Rx0], freq, S_m, S_e)) if verbose: diff --git a/tests/em/fdem/forward/test_FDEM_forward.py b/tests/em/fdem/forward/test_FDEM_forward.py index 21c0d114..da22c950 100644 --- a/tests/em/fdem/forward/test_FDEM_forward.py +++ b/tests/em/fdem/forward/test_FDEM_forward.py @@ -8,6 +8,7 @@ from SimPEG.EM.Utils.testingUtils import getFDEMProblem testEB = True testHJ = True testEJ = True +testBH = True verbose = False TOLEBHJ = 1e-5 @@ -145,5 +146,84 @@ class FDEM_CrossCheck(unittest.TestCase): def test_EJ_CrossCheck_ezi_Jform(self): self.assertTrue(crossCheckTest('e', 'j', 'ezi', TOL=TOLEJHB)) + def test_EJ_CrossCheck_bxr_Jform(self): + self.assertTrue(crossCheckTest('e', 'j', 'bxr', TOL=TOLEJHB)) + def test_EJ_CrossCheck_byr_Jform(self): + self.assertTrue(crossCheckTest('e', 'j', 'byr', TOL=TOLEJHB)) + def test_EJ_CrossCheck_bzr_Jform(self): + self.assertTrue(crossCheckTest('e', 'j', 'bzr', TOL=TOLEJHB)) + def test_EJ_CrossCheck_bxi_Jform(self): + self.assertTrue(crossCheckTest('e', 'j', 'bxi', TOL=TOLEJHB)) + def test_EJ_CrossCheck_byi_Jform(self): + self.assertTrue(crossCheckTest('e', 'j', 'byi', TOL=TOLEJHB)) + def test_EJ_CrossCheck_bzi_Jform(self): + self.assertTrue(crossCheckTest('e', 'j', 'bzi', TOL=TOLEJHB)) + + def test_EJ_CrossCheck_hxr_Jform(self): + self.assertTrue(crossCheckTest('e', 'j', 'hxr', TOL=TOLEJHB)) + def test_EJ_CrossCheck_hyr_Jform(self): + self.assertTrue(crossCheckTest('e', 'j', 'hyr', TOL=TOLEJHB)) + def test_EJ_CrossCheck_hzr_Jform(self): + self.assertTrue(crossCheckTest('e', 'j', 'hzr', TOL=TOLEJHB)) + def test_EJ_CrossCheck_hxi_Jform(self): + self.assertTrue(crossCheckTest('e', 'j', 'hxi', TOL=TOLEJHB)) + def test_EJ_CrossCheck_hyi_Jform(self): + self.assertTrue(crossCheckTest('e', 'j', 'hyi', TOL=TOLEJHB)) + def test_EJ_CrossCheck_hzi_Jform(self): + self.assertTrue(crossCheckTest('e', 'j', 'hzi', TOL=TOLEJHB)) + + if testBH: + def test_BH_CrossCheck_jxr_Jform(self): + self.assertTrue(crossCheckTest('b', 'h', 'jxr', TOL=TOLEJHB)) + def test_BH_CrossCheck_jyr_Jform(self): + self.assertTrue(crossCheckTest('b', 'h', 'jyr', TOL=TOLEJHB)) + def test_BH_CrossCheck_jzr_Jform(self): + self.assertTrue(crossCheckTest('b', 'h', 'jzr', TOL=TOLEJHB)) + def test_BH_CrossCheck_jxi_Jform(self): + self.assertTrue(crossCheckTest('b', 'h', 'jxi', TOL=TOLEJHB)) + def test_BH_CrossCheck_jyi_Jform(self): + self.assertTrue(crossCheckTest('b', 'h', 'jyi', TOL=TOLEJHB)) + def test_BH_CrossCheck_jzi_Jform(self): + self.assertTrue(crossCheckTest('b', 'h', 'jzi', TOL=TOLEJHB)) + + def test_BH_CrossCheck_exr_Jform(self): + self.assertTrue(crossCheckTest('b', 'h', 'exr', TOL=TOLEJHB)) + def test_BH_CrossCheck_eyr_Jform(self): + self.assertTrue(crossCheckTest('b', 'h', 'eyr', TOL=TOLEJHB)) + def test_BH_CrossCheck_ezr_Jform(self): + self.assertTrue(crossCheckTest('b', 'h', 'ezr', TOL=TOLEJHB)) + def test_BH_CrossCheck_exi_Jform(self): + self.assertTrue(crossCheckTest('b', 'h', 'exi', TOL=TOLEJHB)) + def test_BH_CrossCheck_eyi_Jform(self): + self.assertTrue(crossCheckTest('b', 'h', 'eyi', TOL=TOLEJHB)) + def test_BH_CrossCheck_ezi_Jform(self): + self.assertTrue(crossCheckTest('b', 'h', 'ezi', TOL=TOLEJHB)) + + def test_BH_CrossCheck_bxr_Jform(self): + self.assertTrue(crossCheckTest('b', 'h', 'bxr', TOL=TOLEJHB)) + def test_BH_CrossCheck_byr_Jform(self): + self.assertTrue(crossCheckTest('b', 'h', 'byr', TOL=TOLEJHB)) + def test_BH_CrossCheck_bzr_Jform(self): + self.assertTrue(crossCheckTest('b', 'h', 'bzr', TOL=TOLEJHB)) + def test_BH_CrossCheck_bxi_Jform(self): + self.assertTrue(crossCheckTest('b', 'h', 'bxi', TOL=TOLEJHB)) + def test_BH_CrossCheck_byi_Jform(self): + self.assertTrue(crossCheckTest('b', 'h', 'byi', TOL=TOLEJHB)) + def test_BH_CrossCheck_bzi_Jform(self): + self.assertTrue(crossCheckTest('b', 'h', 'bzi', TOL=TOLEJHB)) + + def test_BH_CrossCheck_hxr_Jform(self): + self.assertTrue(crossCheckTest('b', 'h', 'hxr', TOL=TOLEJHB)) + def test_BH_CrossCheck_hyr_Jform(self): + self.assertTrue(crossCheckTest('b', 'h', 'hyr', TOL=TOLEJHB)) + def test_BH_CrossCheck_hzr_Jform(self): + self.assertTrue(crossCheckTest('b', 'h', 'hzr', TOL=TOLEJHB)) + def test_BH_CrossCheck_hxi_Jform(self): + self.assertTrue(crossCheckTest('b', 'h', 'hxi', TOL=TOLEJHB)) + def test_BH_CrossCheck_hyi_Jform(self): + self.assertTrue(crossCheckTest('b', 'h', 'hyi', TOL=TOLEJHB)) + def test_BH_CrossCheck_hzi_Jform(self): + self.assertTrue(crossCheckTest('b', 'h', 'hzi', TOL=TOLEJHB)) + if __name__ == '__main__': unittest.main() \ No newline at end of file From 22a6310a594df388dc89e90dbb9deaa4de1b6ef8 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Thu, 10 Dec 2015 08:51:06 -0800 Subject: [PATCH 05/17] breaking up tests, also pass in u to df_dm --- SimPEG/EM/FDEM/FDEM.py | 4 +- SimPEG/EM/FDEM/FieldsFDEM.py | 74 ++++--- SimPEG/EM/Utils/testingUtils.py | 48 ++++- tests/em/fdem/forward/test_FDEM_forward.py | 202 +++--------------- .../em/fdem/forward/test_FDEM_forwardEJHB.py | 76 +++++++ tests/em/fdem/forward/test_FDEM_forwardHB.py | 76 +++++++ .../inverse/adjoint/test_FDEM_adjointEB.py | 130 +++++++++++ ...FDEM_adjoint.py => test_FDEM_adjointHJ.py} | 52 ----- .../fdem/inverse/derivs/test_FDEM_derivs.py | 80 ++++--- 9 files changed, 446 insertions(+), 296 deletions(-) create mode 100644 tests/em/fdem/forward/test_FDEM_forwardEJHB.py create mode 100644 tests/em/fdem/forward/test_FDEM_forwardHB.py create mode 100644 tests/em/fdem/inverse/adjoint/test_FDEM_adjointEB.py rename tests/em/fdem/inverse/adjoint/{test_FDEM_adjoint.py => test_FDEM_adjointHJ.py} (60%) diff --git a/SimPEG/EM/FDEM/FDEM.py b/SimPEG/EM/FDEM/FDEM.py index f2167fd8..3d6c88a7 100644 --- a/SimPEG/EM/FDEM/FDEM.py +++ b/SimPEG/EM/FDEM/FDEM.py @@ -78,10 +78,10 @@ class BaseFDEMProblem(BaseEMProblem): for rx in src.rxList: df_duFun = getattr(f, '_%sDeriv_u'%rx.projField, None) - df_dudu_dm = df_duFun(src, du_dm, adjoint=False) + df_dudu_dm = df_duFun(src, du_dm, u_src, adjoint=False) df_dmFun = getattr(f, '_%sDeriv_m'%rx.projField, None) - df_dm = df_dmFun(src, v, adjoint=False) + df_dm = df_dmFun(src, v, u_src, adjoint=False) Df_Dm = np.array(df_dudu_dm + df_dm,dtype=complex) diff --git a/SimPEG/EM/FDEM/FieldsFDEM.py b/SimPEG/EM/FDEM/FieldsFDEM.py index fb05b8c6..e7d96dbe 100644 --- a/SimPEG/EM/FDEM/FieldsFDEM.py +++ b/SimPEG/EM/FDEM/FieldsFDEM.py @@ -32,10 +32,9 @@ class Fields_e(Fields): 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._mui = self.survey.prob.curModel.mui self._nC = self.survey.prob.mesh.nC + self._MeSigma = self.survey.prob.MeSigma + self._MeSigmaDeriv = self.survey.prob.MeSigmaDeriv def _GLoc(self,fieldType): if fieldType == 'e': @@ -60,10 +59,10 @@ class Fields_e(Fields): def _e(self, eSolution, srcList): return self._ePrimary(eSolution,srcList) + self._eSecondary(eSolution,srcList) - def _eDeriv_u(self, src, v, adjoint = False): + def _eDeriv_u(self, src, v, eSolution, adjoint = False): return Identity()*v - def _eDeriv_m(self, src, v, adjoint = False): + def _eDeriv_m(self, src, v, eSolution, adjoint = False): # assuming primary does not depend on the model return Zero() @@ -83,13 +82,13 @@ class Fields_e(Fields): b[:,i] = b[:,i]+ 1./(1j*omega(src.freq)) * S_m return b - def _bSecondaryDeriv_u(self, src, v, adjoint = False): + def _bSecondaryDeriv_u(self, src, v, eSolution, adjoint = False): C = self._edgeCurl if adjoint: return - 1./(1j*omega(src.freq)) * (C.T * v) return - 1./(1j*omega(src.freq)) * (C * v) - def _bSecondaryDeriv_m(self, src, v, adjoint = False): + def _bSecondaryDeriv_m(self, src, v, eSolution, adjoint = False): S_mDeriv, _ = src.evalDeriv(self.prob, adjoint) S_mDeriv = S_mDeriv(v) return 1./(1j * omega(src.freq)) * S_mDeriv @@ -97,27 +96,50 @@ class Fields_e(Fields): def _b(self, eSolution, srcList): return self._bPrimary(eSolution, srcList) + self._bSecondary(eSolution, srcList) - def _bDeriv_u(self, src, v, adjoint=False): + def _bDeriv_u(self, src, v, eSolution, adjoint = False): # Primary does not depend on u return self._bSecondaryDeriv_u(src, v, adjoint) - def _bDeriv_m(self, src, v, adjoint=False): + def _bDeriv_m(self, src, v, eSolution, adjoint = False): # 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 - # Sigma = sdiag(np.kron(np.ones(n), sigma)) - Sigma = self.prob.MeSigma + Sigma = self._MeSigma VI = sdiag(1./np.kron(np.ones(n), self.prob.mesh.vol)) e = self._e(eSolution, srcList) - return VI * (aveE2CCV * (Sigma *e) ) - def _h(self, eolution, srcList): + def _jDeriv_u(self, src, eSolution, v, adjoint = False): + aveE2CCV = self._aveE2CCV + n = int(aveE2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy + Sigma = self._MeSigma + + VI = sdiag(1./np.kron(np.ones(n), self.prob.mesh.vol)) + + if not adjoint: + return VI * (aveE2CCV * (Sigma * (self._eDeriv_u(src, v, adjoint) ) ) ) + return self._eDeriv_u(src, Sigma.T * (aveE2CCV.T * (VI.T * v) ), adjoint) + + def _jDeriv_m(self, src, v, eSolution, adjoint = False): + aveE2CCV = self._aveE2CCV + Sigma = self._MeSigma + SigmaDeriv = self._MeSigmaDeriv + e = self._e(eSolution, [src]) + + 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)) + + if not adjoint: + return VI * (aveE2CCV * ( SigmaDeriv(e) * v + self._eDeriv_m(src, v, adjoint) )) + return SigmaDeriv(aveE2CCV.T * (VI.T * e), adjoint) * v + self._eDeriv_m(src, aveE2CCV.T * (VI.T * v), adjoint) + + + def _h(self, eSolution, srcList): b = self._b(eSolution, srcList) Mui = self.survey.prob.MfMui aveF2CCV = self._aveF2CCV @@ -129,30 +151,6 @@ class Fields_e(Fields): return VI * (aveF2CCV * (Mui * b)) - 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'} diff --git a/SimPEG/EM/Utils/testingUtils.py b/SimPEG/EM/Utils/testingUtils.py index 9979417d..ae2b7321 100644 --- a/SimPEG/EM/Utils/testingUtils.py +++ b/SimPEG/EM/Utils/testingUtils.py @@ -4,6 +4,13 @@ from SimPEG import EM import sys from scipy.constants import mu_0 +FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order +CONDUCTIVITY = 1e1 +MU = mu_0 +freq = 5e-1 +addrandoms = False + + def getFDEMProblem(fdemType, comp, SrcList, freq, verbose=False): cs = 10. ncx, ncy, ncz = 0, 0, 0 @@ -72,4 +79,43 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, verbose=False): except ImportError, e: prb.Solver = SolverLU - return prb \ No newline at end of file + return prb + +def crossCheckTest(SrcList, fdemType1, fdemType2, comp, TOL=1e-5, verbose=False): + + l2norm = lambda r: np.sqrt(r.dot(r)) + + prb1 = getFDEMProblem(fdemType1, comp, SrcList, freq, verbose) + mesh = prb1.mesh + print 'Cross Checking Forward: %s, %s formulations - %s' % (fdemType1, fdemType2, comp) + m = np.log(np.ones(mesh.nC)*CONDUCTIVITY) + mu = np.log(np.ones(mesh.nC)*MU) + + if addrandoms is True: + m = m + np.random.randn(mesh.nC)*np.log(CONDUCTIVITY)*1e-1 + mu = mu + np.random.randn(mesh.nC)*MU*1e-1 + + # prb1.PropMap.PropModel.mu = mu + # prb1.PropMap.PropModel.mui = 1./mu + survey1 = prb1.survey + d1 = survey1.dpred(m) + + if verbose: + print ' Problem 1 solved' + + + prb2 = getFDEMProblem(fdemType2, comp, SrcList, freq, verbose) + + # prb2.mu = mu + survey2 = prb2.survey + d2 = survey2.dpred(m) + + if verbose: + print ' Problem 2 solved' + + r = d2-d1 + l2r = l2norm(r) + + tol = np.max([TOL*(10**int(np.log10(l2norm(d1)))),FLR]) + print l2norm(d1), l2norm(d2), l2r , tol, l2r < tol + return l2r < tol diff --git a/tests/em/fdem/forward/test_FDEM_forward.py b/tests/em/fdem/forward/test_FDEM_forward.py index da22c950..da446675 100644 --- a/tests/em/fdem/forward/test_FDEM_forward.py +++ b/tests/em/fdem/forward/test_FDEM_forward.py @@ -3,7 +3,7 @@ from SimPEG import * from SimPEG import EM import sys from scipy.constants import mu_0 -from SimPEG.EM.Utils.testingUtils import getFDEMProblem +from SimPEG.EM.Utils.testingUtils import getFDEMProblem, crossCheckTest testEB = True testHJ = True @@ -15,215 +15,63 @@ TOLEBHJ = 1e-5 TOLEJHB = 1 # averaging and more sensitive to boundary condition violations (ie. the impact of violating the boundary conditions in each case is different.) #TODO: choose better testing parameters to lower this -FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order -CONDUCTIVITY = 1e1 -MU = mu_0 -freq = 5e-1 -addrandoms = False - SrcList = ['RawVec', 'MagDipole_Bfield', 'MagDipole', 'CircularLoop'] -def crossCheckTest(fdemType1, fdemType2, comp, TOL=TOLEBHJ): - - l2norm = lambda r: np.sqrt(r.dot(r)) - - prb1 = getFDEMProblem(fdemType1, comp, SrcList, freq, verbose) - mesh = prb1.mesh - print 'Cross Checking Forward: %s, %s formulations - %s' % (fdemType1, fdemType2, comp) - m = np.log(np.ones(mesh.nC)*CONDUCTIVITY) - mu = np.log(np.ones(mesh.nC)*MU) - - if addrandoms is True: - m = m + np.random.randn(mesh.nC)*np.log(CONDUCTIVITY)*1e-1 - mu = mu + np.random.randn(mesh.nC)*MU*1e-1 - - # prb1.PropMap.PropModel.mu = mu - # prb1.PropMap.PropModel.mui = 1./mu - survey1 = prb1.survey - d1 = survey1.dpred(m) - - if verbose: - print ' Problem 1 solved' - - - prb2 = getFDEMProblem(fdemType2, comp, SrcList, freq, verbose) - - # prb2.mu = mu - survey2 = prb2.survey - d2 = survey2.dpred(m) - - if verbose: - print ' Problem 2 solved' - - r = d2-d1 - l2r = l2norm(r) - - tol = np.max([TOL*(10**int(np.log10(l2norm(d1)))),FLR]) - print l2norm(d1), l2norm(d2), l2r , tol, l2r < tol - return l2r < tol - - class FDEM_CrossCheck(unittest.TestCase): if testEB: def test_EB_CrossCheck_exr_Eform(self): - self.assertTrue(crossCheckTest('e', 'b', 'exr')) + self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'exr', verbose=verbose)) def test_EB_CrossCheck_eyr_Eform(self): - self.assertTrue(crossCheckTest('e', 'b', 'eyr')) + self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'eyr', verbose=verbose)) def test_EB_CrossCheck_ezr_Eform(self): - self.assertTrue(crossCheckTest('e', 'b', 'ezr')) + self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'ezr', verbose=verbose)) def test_EB_CrossCheck_exi_Eform(self): - self.assertTrue(crossCheckTest('e', 'b', 'exi')) + self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'exi', verbose=verbose)) def test_EB_CrossCheck_eyi_Eform(self): - self.assertTrue(crossCheckTest('e', 'b', 'eyi')) + self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'eyi', verbose=verbose)) def test_EB_CrossCheck_ezi_Eform(self): - self.assertTrue(crossCheckTest('e', 'b', 'ezi')) + self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'ezi', verbose=verbose)) def test_EB_CrossCheck_bxr_Eform(self): - self.assertTrue(crossCheckTest('e', 'b', 'bxr')) + self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'bxr', verbose=verbose)) def test_EB_CrossCheck_byr_Eform(self): - self.assertTrue(crossCheckTest('e', 'b', 'byr')) + self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'byr', verbose=verbose)) def test_EB_CrossCheck_bzr_Eform(self): - self.assertTrue(crossCheckTest('e', 'b', 'bzr')) + self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'bzr', verbose=verbose)) def test_EB_CrossCheck_bxi_Eform(self): - self.assertTrue(crossCheckTest('e', 'b', 'bxi')) + self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'bxi', verbose=verbose)) def test_EB_CrossCheck_byi_Eform(self): - self.assertTrue(crossCheckTest('e', 'b', 'byi')) + self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'byi', verbose=verbose)) def test_EB_CrossCheck_bzi_Eform(self): - self.assertTrue(crossCheckTest('e', 'b', 'bzi')) + self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'bzi', verbose=verbose)) if testHJ: def test_HJ_CrossCheck_jxr_Jform(self): - self.assertTrue(crossCheckTest('j', 'h', 'jxr')) + self.assertTrue(crossCheckTest(SrcList, 'j', 'h', 'jxr', verbose=verbose)) def test_HJ_CrossCheck_jyr_Jform(self): - self.assertTrue(crossCheckTest('j', 'h', 'jyr')) + self.assertTrue(crossCheckTest(SrcList, 'j', 'h', 'jyr', verbose=verbose)) def test_HJ_CrossCheck_jzr_Jform(self): - self.assertTrue(crossCheckTest('j', 'h', 'jzr')) + self.assertTrue(crossCheckTest(SrcList, 'j', 'h', 'jzr', verbose=verbose)) def test_HJ_CrossCheck_jxi_Jform(self): - self.assertTrue(crossCheckTest('j', 'h', 'jxi')) + self.assertTrue(crossCheckTest(SrcList, 'j', 'h', 'jxi', verbose=verbose)) def test_HJ_CrossCheck_jyi_Jform(self): - self.assertTrue(crossCheckTest('j', 'h', 'jyi')) + self.assertTrue(crossCheckTest(SrcList, 'j', 'h', 'jyi', verbose=verbose)) def test_HJ_CrossCheck_jzi_Jform(self): - self.assertTrue(crossCheckTest('j', 'h', 'jzi')) + self.assertTrue(crossCheckTest(SrcList, 'j', 'h', 'jzi', verbose=verbose)) def test_HJ_CrossCheck_hxr_Jform(self): - self.assertTrue(crossCheckTest('j', 'h', 'hxr')) + self.assertTrue(crossCheckTest(SrcList, 'j', 'h', 'hxr', verbose=verbose)) def test_HJ_CrossCheck_hyr_Jform(self): - self.assertTrue(crossCheckTest('j', 'h', 'hyr')) + self.assertTrue(crossCheckTest(SrcList, 'j', 'h', 'hyr', verbose=verbose)) def test_HJ_CrossCheck_hzr_Jform(self): - self.assertTrue(crossCheckTest('j', 'h', 'hzr')) + self.assertTrue(crossCheckTest(SrcList, 'j', 'h', 'hzr', verbose=verbose)) def test_HJ_CrossCheck_hxi_Jform(self): - self.assertTrue(crossCheckTest('j', 'h', 'hxi')) + self.assertTrue(crossCheckTest(SrcList, 'j', 'h', 'hxi', verbose=verbose)) def test_HJ_CrossCheck_hyi_Jform(self): - self.assertTrue(crossCheckTest('j', 'h', 'hyi')) + self.assertTrue(crossCheckTest(SrcList, 'j', 'h', 'hyi', verbose=verbose)) def test_HJ_CrossCheck_hzi_Jform(self): - self.assertTrue(crossCheckTest('j', 'h', 'hzi')) - - if testEJ: - def test_EJ_CrossCheck_jxr_Jform(self): - self.assertTrue(crossCheckTest('e', 'j', 'jxr', TOL=TOLEJHB)) - def test_EJ_CrossCheck_jyr_Jform(self): - self.assertTrue(crossCheckTest('e', 'j', 'jyr', TOL=TOLEJHB)) - def test_EJ_CrossCheck_jzr_Jform(self): - self.assertTrue(crossCheckTest('e', 'j', 'jzr', TOL=TOLEJHB)) - def test_EJ_CrossCheck_jxi_Jform(self): - self.assertTrue(crossCheckTest('e', 'j', 'jxi', TOL=TOLEJHB)) - def test_EJ_CrossCheck_jyi_Jform(self): - self.assertTrue(crossCheckTest('e', 'j', 'jyi', TOL=TOLEJHB)) - def test_EJ_CrossCheck_jzi_Jform(self): - self.assertTrue(crossCheckTest('e', 'j', 'jzi', TOL=TOLEJHB)) - - def test_EJ_CrossCheck_exr_Jform(self): - self.assertTrue(crossCheckTest('e', 'j', 'exr', TOL=TOLEJHB)) - def test_EJ_CrossCheck_eyr_Jform(self): - self.assertTrue(crossCheckTest('e', 'j', 'eyr', TOL=TOLEJHB)) - def test_EJ_CrossCheck_ezr_Jform(self): - self.assertTrue(crossCheckTest('e', 'j', 'ezr', TOL=TOLEJHB)) - def test_EJ_CrossCheck_exi_Jform(self): - self.assertTrue(crossCheckTest('e', 'j', 'exi', TOL=TOLEJHB)) - def test_EJ_CrossCheck_eyi_Jform(self): - self.assertTrue(crossCheckTest('e', 'j', 'eyi', TOL=TOLEJHB)) - def test_EJ_CrossCheck_ezi_Jform(self): - self.assertTrue(crossCheckTest('e', 'j', 'ezi', TOL=TOLEJHB)) - - def test_EJ_CrossCheck_bxr_Jform(self): - self.assertTrue(crossCheckTest('e', 'j', 'bxr', TOL=TOLEJHB)) - def test_EJ_CrossCheck_byr_Jform(self): - self.assertTrue(crossCheckTest('e', 'j', 'byr', TOL=TOLEJHB)) - def test_EJ_CrossCheck_bzr_Jform(self): - self.assertTrue(crossCheckTest('e', 'j', 'bzr', TOL=TOLEJHB)) - def test_EJ_CrossCheck_bxi_Jform(self): - self.assertTrue(crossCheckTest('e', 'j', 'bxi', TOL=TOLEJHB)) - def test_EJ_CrossCheck_byi_Jform(self): - self.assertTrue(crossCheckTest('e', 'j', 'byi', TOL=TOLEJHB)) - def test_EJ_CrossCheck_bzi_Jform(self): - self.assertTrue(crossCheckTest('e', 'j', 'bzi', TOL=TOLEJHB)) - - def test_EJ_CrossCheck_hxr_Jform(self): - self.assertTrue(crossCheckTest('e', 'j', 'hxr', TOL=TOLEJHB)) - def test_EJ_CrossCheck_hyr_Jform(self): - self.assertTrue(crossCheckTest('e', 'j', 'hyr', TOL=TOLEJHB)) - def test_EJ_CrossCheck_hzr_Jform(self): - self.assertTrue(crossCheckTest('e', 'j', 'hzr', TOL=TOLEJHB)) - def test_EJ_CrossCheck_hxi_Jform(self): - self.assertTrue(crossCheckTest('e', 'j', 'hxi', TOL=TOLEJHB)) - def test_EJ_CrossCheck_hyi_Jform(self): - self.assertTrue(crossCheckTest('e', 'j', 'hyi', TOL=TOLEJHB)) - def test_EJ_CrossCheck_hzi_Jform(self): - self.assertTrue(crossCheckTest('e', 'j', 'hzi', TOL=TOLEJHB)) - - if testBH: - def test_BH_CrossCheck_jxr_Jform(self): - self.assertTrue(crossCheckTest('b', 'h', 'jxr', TOL=TOLEJHB)) - def test_BH_CrossCheck_jyr_Jform(self): - self.assertTrue(crossCheckTest('b', 'h', 'jyr', TOL=TOLEJHB)) - def test_BH_CrossCheck_jzr_Jform(self): - self.assertTrue(crossCheckTest('b', 'h', 'jzr', TOL=TOLEJHB)) - def test_BH_CrossCheck_jxi_Jform(self): - self.assertTrue(crossCheckTest('b', 'h', 'jxi', TOL=TOLEJHB)) - def test_BH_CrossCheck_jyi_Jform(self): - self.assertTrue(crossCheckTest('b', 'h', 'jyi', TOL=TOLEJHB)) - def test_BH_CrossCheck_jzi_Jform(self): - self.assertTrue(crossCheckTest('b', 'h', 'jzi', TOL=TOLEJHB)) - - def test_BH_CrossCheck_exr_Jform(self): - self.assertTrue(crossCheckTest('b', 'h', 'exr', TOL=TOLEJHB)) - def test_BH_CrossCheck_eyr_Jform(self): - self.assertTrue(crossCheckTest('b', 'h', 'eyr', TOL=TOLEJHB)) - def test_BH_CrossCheck_ezr_Jform(self): - self.assertTrue(crossCheckTest('b', 'h', 'ezr', TOL=TOLEJHB)) - def test_BH_CrossCheck_exi_Jform(self): - self.assertTrue(crossCheckTest('b', 'h', 'exi', TOL=TOLEJHB)) - def test_BH_CrossCheck_eyi_Jform(self): - self.assertTrue(crossCheckTest('b', 'h', 'eyi', TOL=TOLEJHB)) - def test_BH_CrossCheck_ezi_Jform(self): - self.assertTrue(crossCheckTest('b', 'h', 'ezi', TOL=TOLEJHB)) - - def test_BH_CrossCheck_bxr_Jform(self): - self.assertTrue(crossCheckTest('b', 'h', 'bxr', TOL=TOLEJHB)) - def test_BH_CrossCheck_byr_Jform(self): - self.assertTrue(crossCheckTest('b', 'h', 'byr', TOL=TOLEJHB)) - def test_BH_CrossCheck_bzr_Jform(self): - self.assertTrue(crossCheckTest('b', 'h', 'bzr', TOL=TOLEJHB)) - def test_BH_CrossCheck_bxi_Jform(self): - self.assertTrue(crossCheckTest('b', 'h', 'bxi', TOL=TOLEJHB)) - def test_BH_CrossCheck_byi_Jform(self): - self.assertTrue(crossCheckTest('b', 'h', 'byi', TOL=TOLEJHB)) - def test_BH_CrossCheck_bzi_Jform(self): - self.assertTrue(crossCheckTest('b', 'h', 'bzi', TOL=TOLEJHB)) - - def test_BH_CrossCheck_hxr_Jform(self): - self.assertTrue(crossCheckTest('b', 'h', 'hxr', TOL=TOLEJHB)) - def test_BH_CrossCheck_hyr_Jform(self): - self.assertTrue(crossCheckTest('b', 'h', 'hyr', TOL=TOLEJHB)) - def test_BH_CrossCheck_hzr_Jform(self): - self.assertTrue(crossCheckTest('b', 'h', 'hzr', TOL=TOLEJHB)) - def test_BH_CrossCheck_hxi_Jform(self): - self.assertTrue(crossCheckTest('b', 'h', 'hxi', TOL=TOLEJHB)) - def test_BH_CrossCheck_hyi_Jform(self): - self.assertTrue(crossCheckTest('b', 'h', 'hyi', TOL=TOLEJHB)) - def test_BH_CrossCheck_hzi_Jform(self): - self.assertTrue(crossCheckTest('b', 'h', 'hzi', TOL=TOLEJHB)) + self.assertTrue(crossCheckTest(SrcList, 'j', 'h', 'hzi', verbose=verbose)) if __name__ == '__main__': unittest.main() \ No newline at end of file diff --git a/tests/em/fdem/forward/test_FDEM_forwardEJHB.py b/tests/em/fdem/forward/test_FDEM_forwardEJHB.py new file mode 100644 index 00000000..be389282 --- /dev/null +++ b/tests/em/fdem/forward/test_FDEM_forwardEJHB.py @@ -0,0 +1,76 @@ +import unittest +from SimPEG import * +from SimPEG import EM +import sys +from scipy.constants import mu_0 +from SimPEG.EM.Utils.testingUtils import getFDEMProblem, crossCheckTest + +testEB = True +testHJ = True +testEJ = True +testBH = True +verbose = False + +TOLEBHJ = 1e-5 +TOLEJHB = 1 # averaging and more sensitive to boundary condition violations (ie. the impact of violating the boundary conditions in each case is different.) +#TODO: choose better testing parameters to lower this + +SrcList = ['RawVec', 'MagDipole_Bfield', 'MagDipole', 'CircularLoop'] + + +class FDEM_CrossCheck(unittest.TestCase): + if testEJ: + def test_EJ_CrossCheck_jxr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'jxr', verbose=verbose, TOL=TOLEJHB)) + def test_EJ_CrossCheck_jyr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'jyr', verbose=verbose, TOL=TOLEJHB)) + def test_EJ_CrossCheck_jzr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'jzr', verbose=verbose, TOL=TOLEJHB)) + def test_EJ_CrossCheck_jxi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'jxi', verbose=verbose, TOL=TOLEJHB)) + def test_EJ_CrossCheck_jyi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'jyi', verbose=verbose, TOL=TOLEJHB)) + def test_EJ_CrossCheck_jzi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'jzi', verbose=verbose, TOL=TOLEJHB)) + + def test_EJ_CrossCheck_exr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'exr', verbose=verbose, TOL=TOLEJHB)) + def test_EJ_CrossCheck_eyr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'eyr', verbose=verbose, TOL=TOLEJHB)) + def test_EJ_CrossCheck_ezr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'ezr', verbose=verbose, TOL=TOLEJHB)) + def test_EJ_CrossCheck_exi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'exi', verbose=verbose, TOL=TOLEJHB)) + def test_EJ_CrossCheck_eyi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'eyi', verbose=verbose, TOL=TOLEJHB)) + def test_EJ_CrossCheck_ezi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'ezi', verbose=verbose, TOL=TOLEJHB)) + + def test_EJ_CrossCheck_bxr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'bxr', verbose=verbose, TOL=TOLEJHB)) + def test_EJ_CrossCheck_byr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'byr', verbose=verbose, TOL=TOLEJHB)) + def test_EJ_CrossCheck_bzr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'bzr', verbose=verbose, TOL=TOLEJHB)) + def test_EJ_CrossCheck_bxi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'bxi', verbose=verbose, TOL=TOLEJHB)) + def test_EJ_CrossCheck_byi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'byi', verbose=verbose, TOL=TOLEJHB)) + def test_EJ_CrossCheck_bzi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'bzi', verbose=verbose, TOL=TOLEJHB)) + + def test_EJ_CrossCheck_hxr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'hxr', verbose=verbose, TOL=TOLEJHB)) + def test_EJ_CrossCheck_hyr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'hyr', verbose=verbose, TOL=TOLEJHB)) + def test_EJ_CrossCheck_hzr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'hzr', verbose=verbose, TOL=TOLEJHB)) + def test_EJ_CrossCheck_hxi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'hxi', verbose=verbose, TOL=TOLEJHB)) + def test_EJ_CrossCheck_hyi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'hyi', verbose=verbose, TOL=TOLEJHB)) + def test_EJ_CrossCheck_hzi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'hzi', verbose=verbose, TOL=TOLEJHB)) + +if __name__ == '__main__': + unittest.main() \ No newline at end of file diff --git a/tests/em/fdem/forward/test_FDEM_forwardHB.py b/tests/em/fdem/forward/test_FDEM_forwardHB.py new file mode 100644 index 00000000..78b3e376 --- /dev/null +++ b/tests/em/fdem/forward/test_FDEM_forwardHB.py @@ -0,0 +1,76 @@ +import unittest +from SimPEG import * +from SimPEG import EM +import sys +from scipy.constants import mu_0 +from SimPEG.EM.Utils.testingUtils import getFDEMProblem, crossCheckTest + +testEB = True +testHJ = True +testEJ = True +testBH = True +verbose = False + +TOLEBHJ = 1e-5 +TOLEJHB = 1 # averaging and more sensitive to boundary condition violations (ie. the impact of violating the boundary conditions in each case is different.) +#TODO: choose better testing parameters to lower this + +SrcList = ['RawVec', 'MagDipole_Bfield', 'MagDipole', 'CircularLoop'] + + +class FDEM_CrossCheck(unittest.TestCase): + if testBH: + def test_BH_CrossCheck_jxr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jxr', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_jyr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jyr', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_jzr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jzr', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_jxi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jxi', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_jyi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jyi', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_jzi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jzi', verbose=verbose, TOL=TOLEJHB)) + + def test_BH_CrossCheck_exr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'exr', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_eyr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'eyr', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_ezr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'ezr', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_exi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'exi', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_eyi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'eyi', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_ezi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'ezi', verbose=verbose, TOL=TOLEJHB)) + + def test_BH_CrossCheck_bxr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'bxr', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_byr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'byr', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_bzr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'bzr', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_bxi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'bxi', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_byi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'byi', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_bzi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'bzi', verbose=verbose, TOL=TOLEJHB)) + + def test_BH_CrossCheck_hxr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hxr', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_hyr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hyr', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_hzr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hzr', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_hxi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hxi', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_hyi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hyi', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_hzi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hzi', verbose=verbose, TOL=TOLEJHB)) + +if __name__ == '__main__': + unittest.main() \ No newline at end of file diff --git a/tests/em/fdem/inverse/adjoint/test_FDEM_adjointEB.py b/tests/em/fdem/inverse/adjoint/test_FDEM_adjointEB.py new file mode 100644 index 00000000..f917845a --- /dev/null +++ b/tests/em/fdem/inverse/adjoint/test_FDEM_adjointEB.py @@ -0,0 +1,130 @@ +import unittest +from SimPEG import * +from SimPEG import EM +import sys +from scipy.constants import mu_0 +from SimPEG.EM.Utils.testingUtils import getFDEMProblem + +testE = True +testB = False + +verbose = False + +TOL = 1e-5 +FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order +CONDUCTIVITY = 1e1 +MU = mu_0 +freq = 1e-1 +addrandoms = True + +SrcType = 'RawVec' #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVec' + +def adjointTest(fdemType, comp): + prb = getFDEMProblem(fdemType, comp, [SrcType], freq) + print 'Adjoint %s formulation - %s' % (fdemType, comp) + + m = np.log(np.ones(prb.mapping.nP)*CONDUCTIVITY) + mu = np.ones(prb.mesh.nC)*MU + + if addrandoms is True: + m = m + np.random.randn(prb.mapping.nP)*np.log(CONDUCTIVITY)*1e-1 + mu = mu + np.random.randn(prb.mesh.nC)*MU*1e-1 + + survey = prb.survey + # prb.PropMap.PropModel.mu = mu + # prb.PropMap.PropModel.mui = 1./mu + u = prb.fields(m) + + v = np.random.rand(survey.nD) + w = np.random.rand(prb.mesh.nC) + + vJw = v.dot(prb.Jvec(m, w, u)) + wJtv = w.dot(prb.Jtvec(m, v, u)) + tol = np.max([TOL*(10**int(np.log10(np.abs(vJw)))),FLR]) + print vJw, wJtv, vJw - wJtv, tol, np.abs(vJw - wJtv) < tol + return np.abs(vJw - wJtv) < tol + +class FDEM_AdjointTests(unittest.TestCase): + if testE: + # def test_Jtvec_adjointTest_exr_Eform(self): + # self.assertTrue(adjointTest('e', 'exr')) + # def test_Jtvec_adjointTest_eyr_Eform(self): + # self.assertTrue(adjointTest('e', 'eyr')) + # def test_Jtvec_adjointTest_ezr_Eform(self): + # self.assertTrue(adjointTest('e', 'ezr')) + # def test_Jtvec_adjointTest_exi_Eform(self): + # self.assertTrue(adjointTest('e', 'exi')) + # def test_Jtvec_adjointTest_eyi_Eform(self): + # self.assertTrue(adjointTest('e', 'eyi')) + # def test_Jtvec_adjointTest_ezi_Eform(self): + # self.assertTrue(adjointTest('e', 'ezi')) + + # def test_Jtvec_adjointTest_bxr_Eform(self): + # self.assertTrue(adjointTest('e', 'bxr')) + # def test_Jtvec_adjointTest_byr_Eform(self): + # self.assertTrue(adjointTest('e', 'byr')) + # def test_Jtvec_adjointTest_bzr_Eform(self): + # self.assertTrue(adjointTest('e', 'bzr')) + # def test_Jtvec_adjointTest_bxi_Eform(self): + # self.assertTrue(adjointTest('e', 'bxi')) + # def test_Jtvec_adjointTest_byi_Eform(self): + # self.assertTrue(adjointTest('e', 'byi')) + # def test_Jtvec_adjointTest_bzi_Eform(self): + # self.assertTrue(adjointTest('e', 'bzi')) + + def test_Jtvec_adjointTest_exr_Eform(self): + self.assertTrue(adjointTest('e', 'jxr')) + # def test_Jtvec_adjointTest_eyr_Eform(self): + # self.assertTrue(adjointTest('e', 'jyr')) + # def test_Jtvec_adjointTest_ezr_Eform(self): + # self.assertTrue(adjointTest('e', 'jzr')) + # def test_Jtvec_adjointTest_exi_Eform(self): + # self.assertTrue(adjointTest('e', 'jxi')) + # def test_Jtvec_adjointTest_eyi_Eform(self): + # self.assertTrue(adjointTest('e', 'jyi')) + # def test_Jtvec_adjointTest_ezi_Eform(self): + # self.assertTrue(adjointTest('e', 'jzi')) + + # def test_Jtvec_adjointTest_bxr_Eform(self): + # self.assertTrue(adjointTest('e', 'hxr')) + # def test_Jtvec_adjointTest_byr_Eform(self): + # self.assertTrue(adjointTest('e', 'hyr')) + # def test_Jtvec_adjointTest_bzr_Eform(self): + # self.assertTrue(adjointTest('e', 'hzr')) + # def test_Jtvec_adjointTest_bxi_Eform(self): + # self.assertTrue(adjointTest('e', 'hxi')) + # def test_Jtvec_adjointTest_byi_Eform(self): + # self.assertTrue(adjointTest('e', 'hyi')) + # def test_Jtvec_adjointTest_bzi_Eform(self): + # self.assertTrue(adjointTest('e', 'hzi')) + + if testB: + def test_Jtvec_adjointTest_exr_Bform(self): + self.assertTrue(adjointTest('b', 'exr')) + def test_Jtvec_adjointTest_eyr_Bform(self): + self.assertTrue(adjointTest('b', 'eyr')) + def test_Jtvec_adjointTest_ezr_Bform(self): + self.assertTrue(adjointTest('b', 'ezr')) + def test_Jtvec_adjointTest_exi_Bform(self): + self.assertTrue(adjointTest('b', 'exi')) + def test_Jtvec_adjointTest_eyi_Bform(self): + self.assertTrue(adjointTest('b', 'eyi')) + def test_Jtvec_adjointTest_ezi_Bform(self): + self.assertTrue(adjointTest('b', 'ezi')) + + def test_Jtvec_adjointTest_bxr_Bform(self): + self.assertTrue(adjointTest('b', 'bxr')) + def test_Jtvec_adjointTest_byr_Bform(self): + self.assertTrue(adjointTest('b', 'byr')) + def test_Jtvec_adjointTest_bzr_Bform(self): + self.assertTrue(adjointTest('b', 'bzr')) + def test_Jtvec_adjointTest_bxi_Bform(self): + self.assertTrue(adjointTest('b', 'bxi')) + def test_Jtvec_adjointTest_byi_Bform(self): + self.assertTrue(adjointTest('b', 'byi')) + def test_Jtvec_adjointTest_bzi_Bform(self): + self.assertTrue(adjointTest('b', 'bzi')) + + +if __name__ == '__main__': + unittest.main() diff --git a/tests/em/fdem/inverse/adjoint/test_FDEM_adjoint.py b/tests/em/fdem/inverse/adjoint/test_FDEM_adjointHJ.py similarity index 60% rename from tests/em/fdem/inverse/adjoint/test_FDEM_adjoint.py rename to tests/em/fdem/inverse/adjoint/test_FDEM_adjointHJ.py index f77f2131..68224e34 100644 --- a/tests/em/fdem/inverse/adjoint/test_FDEM_adjoint.py +++ b/tests/em/fdem/inverse/adjoint/test_FDEM_adjointHJ.py @@ -45,58 +45,6 @@ def adjointTest(fdemType, comp): return np.abs(vJw - wJtv) < tol class FDEM_AdjointTests(unittest.TestCase): - if testEB: - def test_Jtvec_adjointTest_exr_Eform(self): - self.assertTrue(adjointTest('e', 'exr')) - def test_Jtvec_adjointTest_eyr_Eform(self): - self.assertTrue(adjointTest('e', 'eyr')) - def test_Jtvec_adjointTest_ezr_Eform(self): - self.assertTrue(adjointTest('e', 'ezr')) - def test_Jtvec_adjointTest_exi_Eform(self): - self.assertTrue(adjointTest('e', 'exi')) - def test_Jtvec_adjointTest_eyi_Eform(self): - self.assertTrue(adjointTest('e', 'eyi')) - def test_Jtvec_adjointTest_ezi_Eform(self): - self.assertTrue(adjointTest('e', 'ezi')) - - def test_Jtvec_adjointTest_bxr_Eform(self): - self.assertTrue(adjointTest('e', 'bxr')) - def test_Jtvec_adjointTest_byr_Eform(self): - self.assertTrue(adjointTest('e', 'byr')) - def test_Jtvec_adjointTest_bzr_Eform(self): - self.assertTrue(adjointTest('e', 'bzr')) - def test_Jtvec_adjointTest_bxi_Eform(self): - self.assertTrue(adjointTest('e', 'bxi')) - def test_Jtvec_adjointTest_byi_Eform(self): - self.assertTrue(adjointTest('e', 'byi')) - def test_Jtvec_adjointTest_bzi_Eform(self): - self.assertTrue(adjointTest('e', 'bzi')) - - def test_Jtvec_adjointTest_exr_Bform(self): - self.assertTrue(adjointTest('b', 'exr')) - def test_Jtvec_adjointTest_eyr_Bform(self): - self.assertTrue(adjointTest('b', 'eyr')) - def test_Jtvec_adjointTest_ezr_Bform(self): - self.assertTrue(adjointTest('b', 'ezr')) - def test_Jtvec_adjointTest_exi_Bform(self): - self.assertTrue(adjointTest('b', 'exi')) - def test_Jtvec_adjointTest_eyi_Bform(self): - self.assertTrue(adjointTest('b', 'eyi')) - def test_Jtvec_adjointTest_ezi_Bform(self): - self.assertTrue(adjointTest('b', 'ezi')) - def test_Jtvec_adjointTest_bxr_Bform(self): - self.assertTrue(adjointTest('b', 'bxr')) - def test_Jtvec_adjointTest_byr_Bform(self): - self.assertTrue(adjointTest('b', 'byr')) - def test_Jtvec_adjointTest_bzr_Bform(self): - self.assertTrue(adjointTest('b', 'bzr')) - def test_Jtvec_adjointTest_bxi_Bform(self): - self.assertTrue(adjointTest('b', 'bxi')) - def test_Jtvec_adjointTest_byi_Bform(self): - self.assertTrue(adjointTest('b', 'byi')) - def test_Jtvec_adjointTest_bzi_Bform(self): - self.assertTrue(adjointTest('b', 'bzi')) - if testHJ: def test_Jtvec_adjointTest_jxr_Jform(self): diff --git a/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py b/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py index 52108c4e..9831d279 100644 --- a/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py +++ b/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py @@ -6,8 +6,9 @@ from scipy.constants import mu_0 from SimPEG.EM.Utils.testingUtils import getFDEMProblem testDerivs = True -testEB = True -testHJ = True +testE = True +testB = False +testHJ = False verbose = False @@ -43,33 +44,60 @@ def derivTest(fdemType, comp): class FDEM_DerivTests(unittest.TestCase): - if testEB: + if testE: + # def test_Jvec_exr_Eform(self): + # self.assertTrue(derivTest('e', 'exr')) + # def test_Jvec_eyr_Eform(self): + # self.assertTrue(derivTest('e', 'eyr')) + # def test_Jvec_ezr_Eform(self): + # self.assertTrue(derivTest('e', 'ezr')) + # def test_Jvec_exi_Eform(self): + # self.assertTrue(derivTest('e', 'exi')) + # def test_Jvec_eyi_Eform(self): + # self.assertTrue(derivTest('e', 'eyi')) + # def test_Jvec_ezi_Eform(self): + # self.assertTrue(derivTest('e', 'ezi')) + + # def test_Jvec_bxr_Eform(self): + # self.assertTrue(derivTest('e', 'bxr')) + # def test_Jvec_byr_Eform(self): + # self.assertTrue(derivTest('e', 'byr')) + # def test_Jvec_bzr_Eform(self): + # self.assertTrue(derivTest('e', 'bzr')) + # def test_Jvec_bxi_Eform(self): + # self.assertTrue(derivTest('e', 'bxi')) + # def test_Jvec_byi_Eform(self): + # self.assertTrue(derivTest('e', 'byi')) + # def test_Jvec_bzi_Eform(self): + # self.assertTrue(derivTest('e', 'bzi')) + def test_Jvec_exr_Eform(self): - self.assertTrue(derivTest('e', 'exr')) - def test_Jvec_eyr_Eform(self): - self.assertTrue(derivTest('e', 'eyr')) - def test_Jvec_ezr_Eform(self): - self.assertTrue(derivTest('e', 'ezr')) - def test_Jvec_exi_Eform(self): - self.assertTrue(derivTest('e', 'exi')) - def test_Jvec_eyi_Eform(self): - self.assertTrue(derivTest('e', 'eyi')) - def test_Jvec_ezi_Eform(self): - self.assertTrue(derivTest('e', 'ezi')) + self.assertTrue(derivTest('e', 'jxr')) + # def test_Jvec_eyr_Eform(self): + # self.assertTrue(derivTest('e', 'jyr')) + # def test_Jvec_ezr_Eform(self): + # self.assertTrue(derivTest('e', 'jzr')) + # def test_Jvec_exi_Eform(self): + # self.assertTrue(derivTest('e', 'jxi')) + # def test_Jvec_eyi_Eform(self): + # self.assertTrue(derivTest('e', 'jyi')) + # def test_Jvec_ezi_Eform(self): + # self.assertTrue(derivTest('e', 'jzi')) - def test_Jvec_bxr_Eform(self): - self.assertTrue(derivTest('e', 'bxr')) - def test_Jvec_byr_Eform(self): - self.assertTrue(derivTest('e', 'byr')) - def test_Jvec_bzr_Eform(self): - self.assertTrue(derivTest('e', 'bzr')) - def test_Jvec_bxi_Eform(self): - self.assertTrue(derivTest('e', 'bxi')) - def test_Jvec_byi_Eform(self): - self.assertTrue(derivTest('e', 'byi')) - def test_Jvec_bzi_Eform(self): - self.assertTrue(derivTest('e', 'bzi')) + # def test_Jvec_bxr_Eform(self): + # self.assertTrue(derivTest('e', 'hxr')) + # def test_Jvec_byr_Eform(self): + # self.assertTrue(derivTest('e', 'hyr')) + # def test_Jvec_bzr_Eform(self): + # self.assertTrue(derivTest('e', 'hzr')) + # def test_Jvec_bxi_Eform(self): + # self.assertTrue(derivTest('e', 'hxi')) + # def test_Jvec_byi_Eform(self): + # self.assertTrue(derivTest('e', 'hyi')) + # def test_Jvec_bzi_Eform(self): + # self.assertTrue(derivTest('e', 'hzi')) + if testB: def test_Jvec_exr_Bform(self): self.assertTrue(derivTest('b', 'exr')) def test_Jvec_eyr_Bform(self): From bb3f9a6a872a0a79793333c3bb5f936924337318 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sun, 17 Jan 2016 14:54:02 -0800 Subject: [PATCH 06/17] make projGLoc a method not a property --- SimPEG/EM/FDEM/FDEM.py | 18 +++-- SimPEG/EM/FDEM/FieldsFDEM.py | 78 +++++++++---------- SimPEG/EM/FDEM/SurveyFDEM.py | 78 +++++++++++-------- SimPEG/Mesh/TensorMesh.py | 12 ++- .../fdem/inverse/derivs/test_FDEM_derivs.py | 44 +++++------ 5 files changed, 128 insertions(+), 102 deletions(-) diff --git a/SimPEG/EM/FDEM/FDEM.py b/SimPEG/EM/FDEM/FDEM.py index 96fae3c5..7f530814 100644 --- a/SimPEG/EM/FDEM/FDEM.py +++ b/SimPEG/EM/FDEM/FDEM.py @@ -64,30 +64,36 @@ class BaseFDEMProblem(BaseEMProblem): self.curModel = m Jv = self.dataPair(self.survey) + utype = self._fieldType + 'Solution' for freq in self.survey.freqs: A = self.getA(freq) # Ainv = self.Solver(A, **self.solverOpts) for src in self.survey.getSrcByFreq(freq): - ftype = self._fieldType + 'Solution' - u_src = f[src, ftype] + u_src = f[src, utype] dA_dm = self.getADeriv_m(freq, u_src, v) dRHS_dm = self.getRHSDeriv_m(freq, src, v) du_dm = Ainv * ( - dA_dm + dRHS_dm ) for rx in src.rxList: df_duFun = getattr(f, '_%sDeriv_u'%rx.projField, None) - df_dudu_dm = df_duFun(src, du_dm, u_src, adjoint=False) + df_dudu_dm = df_duFun(src, u_src, du_dm, adjoint=False) df_dmFun = getattr(f, '_%sDeriv_m'%rx.projField, None) - df_dm = df_dmFun(src, v, u_src, adjoint=False) + df_dm = df_dmFun(src, u_src, v, adjoint=False) + + # print df_dudu_dm.shape, df_dm.shape, du_dm.shape Df_Dm = np.array(df_dudu_dm + df_dm,dtype=complex) - P = lambda v: rx.projectFieldsDeriv(src, self.mesh, f, v) # wrt u, also have wrt m + print 'getting to P', u_src.shape + # P = lambda v: rx.projectFieldsDeriv(src, self.mesh, f, v) # wrt u, also have wrt m - Jv[src, rx] = P(Df_Dm) + # Jv[src, rx] = P(Df_Dm) + print Df_Dm.shape + print 'should be', m.shape + Jv[src,rx] = rx.projectFieldsDeriv(src, self.mesh, f, Df_Dm) Ainv.clean() return Utils.mkvc(Jv) diff --git a/SimPEG/EM/FDEM/FieldsFDEM.py b/SimPEG/EM/FDEM/FieldsFDEM.py index e7d96dbe..8cbc21c0 100644 --- a/SimPEG/EM/FDEM/FieldsFDEM.py +++ b/SimPEG/EM/FDEM/FieldsFDEM.py @@ -36,13 +36,13 @@ class Fields_e(Fields): self._MeSigma = self.survey.prob.MeSigma self._MeSigmaDeriv = self.survey.prob.MeSigmaDeriv - def _GLoc(self,fieldType): + def _GLoc(self, fieldType): if fieldType == 'e': return 'E' elif fieldType == 'b': return 'F' elif (fieldType == 'h') or (fieldType == 'j'): - return 'CC' + return 'CCV' else: raise Exception('Field type must be e, b, h, j') @@ -59,10 +59,10 @@ class Fields_e(Fields): def _e(self, eSolution, srcList): return self._ePrimary(eSolution,srcList) + self._eSecondary(eSolution,srcList) - def _eDeriv_u(self, src, v, eSolution, adjoint = False): + def _eDeriv_u(self, src, eSolution, v, adjoint = False): return Identity()*v - def _eDeriv_m(self, src, v, eSolution, adjoint = False): + def _eDeriv_m(self, src, eSolution, v, adjoint = False): # assuming primary does not depend on the model return Zero() @@ -82,13 +82,13 @@ class Fields_e(Fields): b[:,i] = b[:,i]+ 1./(1j*omega(src.freq)) * S_m return b - def _bSecondaryDeriv_u(self, src, v, eSolution, adjoint = False): + def _bSecondaryDeriv_u(self, src, eSolution, v, adjoint = False): C = self._edgeCurl if adjoint: return - 1./(1j*omega(src.freq)) * (C.T * v) return - 1./(1j*omega(src.freq)) * (C * v) - def _bSecondaryDeriv_m(self, src, v, eSolution, adjoint = False): + def _bSecondaryDeriv_m(self, src, eSolution, v, adjoint = False): S_mDeriv, _ = src.evalDeriv(self.prob, adjoint) S_mDeriv = S_mDeriv(v) return 1./(1j * omega(src.freq)) * S_mDeriv @@ -96,13 +96,13 @@ class Fields_e(Fields): def _b(self, eSolution, srcList): return self._bPrimary(eSolution, srcList) + self._bSecondary(eSolution, srcList) - def _bDeriv_u(self, src, v, eSolution, adjoint = False): + def _bDeriv_u(self, src, eSolution, v, adjoint = False): # Primary does not depend on u - return self._bSecondaryDeriv_u(src, v, adjoint) + return self._bSecondaryDeriv_u(src, eSolution, v, adjoint) - def _bDeriv_m(self, src, v, eSolution, adjoint = False): + def _bDeriv_m(self, src, eSolution, v, adjoint = False): # Assuming the primary does not depend on the model - return self._bSecondaryDeriv_m(src, v, adjoint) + return self._bSecondaryDeriv_m(src, eSolution, v, adjoint) def _j(self, eSolution, srcList): aveE2CCV = self._aveE2CCV @@ -121,10 +121,10 @@ class Fields_e(Fields): VI = sdiag(1./np.kron(np.ones(n), self.prob.mesh.vol)) if not adjoint: - return VI * (aveE2CCV * (Sigma * (self._eDeriv_u(src, v, adjoint) ) ) ) - return self._eDeriv_u(src, Sigma.T * (aveE2CCV.T * (VI.T * v) ), adjoint) + return VI * (aveE2CCV * (Sigma * (self._eDeriv_u(src, eSolution, v, adjoint) ) ) ) + return self._eDeriv_u(src, eSolution, Sigma.T * (aveE2CCV.T * (VI.T * v) ), adjoint) - def _jDeriv_m(self, src, v, eSolution, adjoint = False): + def _jDeriv_m(self, src, eSolution, v, adjoint = False): aveE2CCV = self._aveE2CCV Sigma = self._MeSigma SigmaDeriv = self._MeSigmaDeriv @@ -135,7 +135,7 @@ class Fields_e(Fields): VI = sdiag(1./np.kron(np.ones(n), self.prob.mesh.vol)) if not adjoint: - return VI * (aveE2CCV * ( SigmaDeriv(e) * v + self._eDeriv_m(src, v, adjoint) )) + return VI * (aveE2CCV * ( self._eDeriv_m(src, v, adjoint) + SigmaDeriv(e) * v)) return SigmaDeriv(aveE2CCV.T * (VI.T * e), adjoint) * v + self._eDeriv_m(src, aveE2CCV.T * (VI.T * v), adjoint) @@ -206,10 +206,10 @@ class Fields_b(Fields): def _b(self, bSolution, srcList): return self._bPrimary(bSolution, srcList) + self._bSecondary(bSolution, srcList) - def _bDeriv_u(self, src, v, adjoint=False): + def _bDeriv_u(self, src, bSolution, v, adjoint=False): return Identity()*v - def _bDeriv_m(self, src, v, adjoint=False): + def _bDeriv_m(self, src, bSolution, v, adjoint=False): # assuming primary does not depend on the model return Zero() @@ -227,13 +227,13 @@ class Fields_b(Fields): e[:,i] = e[:,i]+ -self._MeSigmaI * S_e return e - def _eSecondaryDeriv_u(self, src, v, adjoint=False): + def _eSecondaryDeriv_u(self, src, bSolution, v, adjoint=False): if not adjoint: return self._MeSigmaI * ( self._edgeCurl.T * ( self._MfMui * v) ) else: return self._MfMui.T * (self._edgeCurl * (self._MeSigmaI.T * v)) - def _eSecondaryDeriv_m(self, src, v, adjoint=False): + def _eSecondaryDeriv_m(self, src, bSolution, v, adjoint=False): bSolution = self[[src],'bSolution'] _,S_e = src.eval(self.prob) Me = self._Me @@ -259,12 +259,12 @@ class Fields_b(Fields): def _e(self, bSolution, srcList): return self._ePrimary(bSolution, srcList) + self._eSecondary(bSolution, srcList) - def _eDeriv_u(self, src, v, adjoint=False): - return self._eSecondaryDeriv_u(src, v, adjoint) + def _eDeriv_u(self, src, bSolution, v, adjoint=False): + return self._eSecondaryDeriv_u(src, bSolution, v, adjoint) - def _eDeriv_m(self, src, v, adjoint=False): + def _eDeriv_m(self, src, bSolution, v, adjoint=False): # assuming primary doesn't depend on model - return self._eSecondaryDeriv_m(src, v, adjoint) + return self._eSecondaryDeriv_m(src, bSolution, v, adjoint) def _j(self, bSolution, srcList): sigma = self._sigma @@ -342,10 +342,10 @@ class Fields_j(Fields): def _j(self, jSolution, srcList): return self._jPrimary(jSolution, srcList) + self._jSecondary(jSolution, srcList) - def _jDeriv_u(self, src, v, adjoint=False): + def _jDeriv_u(self, src, u, v, adjoint=False): return Identity()*v - def _jDeriv_m(self, src, v, adjoint=False): + def _jDeriv_m(self, src, u, v, adjoint=False): # assuming primary does not depend on the model return Zero() @@ -364,13 +364,13 @@ class Fields_j(Fields): h[:,i] = h[:,i]+ 1./(1j*omega(src.freq)) * self._MeMuI * (S_m) return h - def _hSecondaryDeriv_u(self, src, v, adjoint=False): + def _hSecondaryDeriv_u(self, src, u, v, adjoint=False): if not adjoint: return -1./(1j*omega(src.freq)) * self._MeMuI * (self._edgeCurl.T * (self._MfRho * v) ) elif adjoint: return -1./(1j*omega(src.freq)) * self._MfRho.T * (self._edgeCurl * ( self._MeMuI.T * v)) - def _hSecondaryDeriv_m(self, src, v, adjoint=False): + def _hSecondaryDeriv_m(self, src, u, v, adjoint=False): jSolution = self[[src],'jSolution'] MeMuI = self._MeMuI C = self._edgeCurl @@ -397,12 +397,12 @@ class Fields_j(Fields): def _h(self, jSolution, srcList): return self._hPrimary(jSolution, srcList) + self._hSecondary(jSolution, srcList) - def _hDeriv_u(self, src, v, adjoint=False): + def _hDeriv_u(self, src, u, v, adjoint=False): return self._hSecondaryDeriv_u(src, v, adjoint) - def _hDeriv_m(self, src, v, adjoint=False): + def _hDeriv_m(self, src, u, v, adjoint=False): # assuming the primary doesn't depend on the model - return self._hSecondaryDeriv_m(src, v, adjoint) + return self._hSecondaryDeriv_m(src, u, v, adjoint) def _e(self, jSolution, srcList): rho = self._rho @@ -417,10 +417,10 @@ class Fields_j(Fields): return VI * (aveF2CCV * (Rho * j)) - def _eDeriv_u(self, src, v, adjoint=False): + def _eDeriv_u(self, src, u, v, adjoint=False): raise NotImplementedError - def _eDeriv_m(self, src, v, adjoint=False): + def _eDeriv_m(self, src, u, v, adjoint=False): raise NotImplementedError def _b(self, jSolution, srcList): @@ -484,10 +484,10 @@ class Fields_h(Fields): def _h(self, hSolution, srcList): return self._hPrimary(hSolution, srcList) + self._hSecondary(hSolution, srcList) - def _hDeriv_u(self, src, v, adjoint=False): + def _hDeriv_u(self, src, u, v, adjoint=False): return Identity()*v - def _hDeriv_m(self, src, v, adjoint=False): + def _hDeriv_m(self, src, u, v, adjoint=False): # assuming primary does not depend on the model return Zero() @@ -505,13 +505,13 @@ class Fields_h(Fields): j[:,i] = j[:,i]+ -S_e return j - def _jSecondaryDeriv_u(self, src, v, adjoint=False): + def _jSecondaryDeriv_u(self, src, u, v, adjoint=False): if not adjoint: return self._edgeCurl*v elif adjoint: return self._edgeCurl.T*v - def _jSecondaryDeriv_m(self, src, v, adjoint=False): + def _jSecondaryDeriv_m(self, src, u, v, adjoint=False): _,S_eDeriv = src.evalDeriv(self.prob, adjoint) S_eDeriv = S_eDeriv(v) return -S_eDeriv @@ -519,10 +519,10 @@ class Fields_h(Fields): def _j(self, hSolution, srcList): return self._jPrimary(hSolution, srcList) + self._jSecondary(hSolution, srcList) - def _jDeriv_u(self, src, v, adjoint=False): + def _jDeriv_u(self, src, u, v, adjoint=False): return self._jSecondaryDeriv_u(src,v,adjoint) - def _jDeriv_m(self, src, v, adjoint=False): + def _jDeriv_m(self, src, u, v, adjoint=False): # assuming the primary does not depend on the model return self._jSecondaryDeriv_m(src,v,adjoint) @@ -539,10 +539,10 @@ class Fields_h(Fields): return VI * (aveF2CCV * (Rho * j)) - def _eDeriv_u(self, src, v, adjoint=False): + def _eDeriv_u(self, src, u, v, adjoint=False): raise NotImplementedError - def _eDeriv_m(self, src, v, adjoint=False): + def _eDeriv_m(self, src, u, v, adjoint=False): raise NotImplementedError def _b(self, hSolution, srcList): diff --git a/SimPEG/EM/FDEM/SurveyFDEM.py b/SimPEG/EM/FDEM/SurveyFDEM.py index 36be9ba5..7e4c77d0 100644 --- a/SimPEG/EM/FDEM/SurveyFDEM.py +++ b/SimPEG/EM/FDEM/SurveyFDEM.py @@ -51,47 +51,51 @@ class Rx(SimPEG.Survey.BaseRx): """Field Type projection (e.g. e b ...)""" return self.knownRxTypes[self.rxType][0] - # @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): """Component projection (real/imag)""" return self.knownRxTypes[self.rxType][2] + def projGLoc(self, u): + """Grid Location projection (e.g. Ex Fy ...)""" + print 'here', u._GLoc(self.rxType[0]) + self.knownRxTypes[self.rxType][1] + return u._GLoc(self.rxType[0]) + self.knownRxTypes[self.rxType][1] + def projectFields(self, src, mesh, u): + # projGLoc = u._GLoc(self.knownRxTypes[self.rxType][0]) + # projGLoc += self.knownRxTypes[self.rxType][1] + + P = self.getP(mesh, self.projGLoc(u)) + 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) + # 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 @@ -99,12 +103,19 @@ class Rx(SimPEG.Survey.BaseRx): projGLoc = u._GLoc(self.knownRxTypes[self.rxType][0]) - if projGLoc != 'CC': - projGLoc += self.knownRxTypes[self.rxType][1] + print self.knownRxTypes[self.rxType][:2], 'Deriv', projGLoc + projGLoc += self.knownRxTypes[self.rxType][1] - P = self.getP(mesh, projGLoc) + # if projGLoc = 'CC': + # P = self.getP(mesh) + # if sel + + # else projGLoc != 'CC': + # projGLoc += self.knownRxTypes[self.rxType][1] + P = self.getP(mesh) if not adjoint: + print P.shape, v.shape Pv_complex = P * v real_or_imag = self.projComp Pv = getattr(Pv_complex, real_or_imag) @@ -174,8 +185,11 @@ class Survey(SimPEG.Survey.BaseSurvey): data = SimPEG.Survey.Data(self) for src in self.srcList: for rx in src.rxList: + print rx.nD + dat = rx.projectFields(src, self.mesh, u) + print dat.shape data[src, rx] = rx.projectFields(src, self.mesh, u) return data def projectFieldsDeriv(self, u): - raise Exception('Use Sources to project fields deriv.') + raise Exception('Use Source to project fields deriv.') diff --git a/SimPEG/Mesh/TensorMesh.py b/SimPEG/Mesh/TensorMesh.py index 874a1f9f..eb2a6965 100644 --- a/SimPEG/Mesh/TensorMesh.py +++ b/SimPEG/Mesh/TensorMesh.py @@ -260,9 +260,15 @@ 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 + elif locType in ['CCVx', 'CCVy', 'CCVz']: + Q = Utils.interpmat(loc, *self.getTensor('CC')) + Zero = 0.*Q + if locType == 'CCVx': + Q = np.r_[Q,Zero,Zero] + elif locType == 'CCVy': + Q = np.r_[Zero,Q,Zero] + elif locType == 'CCVz': + Q = np.r_[Zero,Zero,Q] diff --git a/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py b/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py index f6d0da90..8530d843 100644 --- a/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py +++ b/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py @@ -100,29 +100,29 @@ class FDEM_DerivTests(unittest.TestCase): if testB: def test_Jvec_exr_Bform(self): self.assertTrue(derivTest('b', 'exr')) - def test_Jvec_eyr_Bform(self): - self.assertTrue(derivTest('b', 'eyr')) - def test_Jvec_ezr_Bform(self): - self.assertTrue(derivTest('b', 'ezr')) - def test_Jvec_exi_Bform(self): - self.assertTrue(derivTest('b', 'exi')) - def test_Jvec_eyi_Bform(self): - self.assertTrue(derivTest('b', 'eyi')) - def test_Jvec_ezi_Bform(self): - self.assertTrue(derivTest('b', 'ezi')) + # def test_Jvec_eyr_Bform(self): + # self.assertTrue(derivTest('b', 'eyr')) + # def test_Jvec_ezr_Bform(self): + # self.assertTrue(derivTest('b', 'ezr')) + # def test_Jvec_exi_Bform(self): + # self.assertTrue(derivTest('b', 'exi')) + # def test_Jvec_eyi_Bform(self): + # self.assertTrue(derivTest('b', 'eyi')) + # def test_Jvec_ezi_Bform(self): + # self.assertTrue(derivTest('b', 'ezi')) - def test_Jvec_bxr_Bform(self): - self.assertTrue(derivTest('b', 'bxr')) - def test_Jvec_byr_Bform(self): - self.assertTrue(derivTest('b', 'byr')) - def test_Jvec_bzr_Bform(self): - self.assertTrue(derivTest('b', 'bzr')) - def test_Jvec_bxi_Bform(self): - self.assertTrue(derivTest('b', 'bxi')) - def test_Jvec_byi_Bform(self): - self.assertTrue(derivTest('b', 'byi')) - def test_Jvec_bzi_Bform(self): - self.assertTrue(derivTest('b', 'bzi')) + # def test_Jvec_bxr_Bform(self): + # self.assertTrue(derivTest('b', 'bxr')) + # def test_Jvec_byr_Bform(self): + # self.assertTrue(derivTest('b', 'byr')) + # def test_Jvec_bzr_Bform(self): + # self.assertTrue(derivTest('b', 'bzr')) + # def test_Jvec_bxi_Bform(self): + # self.assertTrue(derivTest('b', 'bxi')) + # def test_Jvec_byi_Bform(self): + # self.assertTrue(derivTest('b', 'byi')) + # def test_Jvec_bzi_Bform(self): + # self.assertTrue(derivTest('b', 'bzi')) if testHJ: def test_Jvec_jxr_Jform(self): From 73edef5eb7a929d564576a1117acd771801d8169 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Thu, 4 Feb 2016 19:34:16 -0800 Subject: [PATCH 07/17] the forward for all fields and fluxes can be computed from any formulation (todo: derive) --- SimPEG/EM/FDEM/FieldsFDEM.py | 10 ++-- SimPEG/EM/FDEM/SurveyFDEM.py | 33 ------------- SimPEG/Mesh/TensorMesh.py | 13 +++++ .../em/fdem/forward/test_FDEM_forwardEJHB.py | 49 +++++++++---------- 4 files changed, 42 insertions(+), 63 deletions(-) diff --git a/SimPEG/EM/FDEM/FieldsFDEM.py b/SimPEG/EM/FDEM/FieldsFDEM.py index 8cbc21c0..914e86f3 100644 --- a/SimPEG/EM/FDEM/FieldsFDEM.py +++ b/SimPEG/EM/FDEM/FieldsFDEM.py @@ -20,8 +20,8 @@ class Fields_e(Fields): 'b' : ['eSolution','F','_b'], 'bPrimary' : ['eSolution','F','_bPrimary'], 'bSecondary' : ['eSolution','F','_bSecondary'], - 'j' : ['eSolution','CC','_j'], - 'h' : ['eSolution','CC','_h'], + 'j' : ['eSolution','CCV','_j'], + 'h' : ['eSolution','CCV','_h'], } def __init__(self,mesh,survey,**kwargs): @@ -189,7 +189,7 @@ class Fields_b(Fields): elif fieldType == 'b': return 'F' elif (fieldType == 'h') or (fieldType == 'j'): - return'CC' + return'CCV' else: raise Exception('Field type must be e, b, h, j') @@ -325,7 +325,7 @@ class Fields_j(Fields): elif fieldType == 'j': return 'F' elif (fieldType == 'e') or (fieldType == 'b'): - return 'CC' + return 'CCV' else: raise Exception('Field type must be e, b, h, j') @@ -467,7 +467,7 @@ class Fields_h(Fields): elif fieldType == 'j': return 'F' elif (fieldType == 'e') or (fieldType == 'b'): - return 'CC' + return 'CCV' else: raise Exception('Field type must be e, b, h, j') diff --git a/SimPEG/EM/FDEM/SurveyFDEM.py b/SimPEG/EM/FDEM/SurveyFDEM.py index 7e4c77d0..dde72497 100644 --- a/SimPEG/EM/FDEM/SurveyFDEM.py +++ b/SimPEG/EM/FDEM/SurveyFDEM.py @@ -58,7 +58,6 @@ class Rx(SimPEG.Survey.BaseRx): def projGLoc(self, u): """Grid Location projection (e.g. Ex Fy ...)""" - print 'here', u._GLoc(self.rxType[0]) + self.knownRxTypes[self.rxType][1] return u._GLoc(self.rxType[0]) + self.knownRxTypes[self.rxType][1] def projectFields(self, src, mesh, u): @@ -72,30 +71,6 @@ class Rx(SimPEG.Survey.BaseRx): # get the real or imag component real_or_imag = self.projComp u_part = getattr(u_part_complex, real_or_imag) - - - - # 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 @@ -106,12 +81,6 @@ class Rx(SimPEG.Survey.BaseRx): print self.knownRxTypes[self.rxType][:2], 'Deriv', projGLoc projGLoc += self.knownRxTypes[self.rxType][1] - # if projGLoc = 'CC': - # P = self.getP(mesh) - # if sel - - # else projGLoc != 'CC': - # projGLoc += self.knownRxTypes[self.rxType][1] P = self.getP(mesh) if not adjoint: @@ -185,9 +154,7 @@ class Survey(SimPEG.Survey.BaseSurvey): data = SimPEG.Survey.Data(self) for src in self.srcList: for rx in src.rxList: - print rx.nD dat = rx.projectFields(src, self.mesh, u) - print dat.shape data[src, rx] = rx.projectFields(src, self.mesh, u) return data diff --git a/SimPEG/Mesh/TensorMesh.py b/SimPEG/Mesh/TensorMesh.py index 508f015c..1650b5cb 100644 --- a/SimPEG/Mesh/TensorMesh.py +++ b/SimPEG/Mesh/TensorMesh.py @@ -234,6 +234,9 @@ 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 vector field defined on cell centers + 'CCVy' -> y-component of vector field defined on cell centers + 'CCVz' -> z-component of vector 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) @@ -257,6 +260,16 @@ 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, *self.getTensor('CC')) + Z = Utils.spzeros(loc.shape[0],self.nC) + if locType == 'CCVx': + Q = sp.hstack([Q,Z,Z]) + elif locType == 'CCVy': + Q = sp.hstack([Z,Q,Z]) + elif locType == 'CCVz': + Q = sp.hstack([Z,Z,Q]) + else: raise NotImplementedError('getInterpolationMat: locType=='+locType+' and mesh.dim=='+str(self.dim)) diff --git a/tests/em/fdem/forward/test_FDEM_forwardEJHB.py b/tests/em/fdem/forward/test_FDEM_forwardEJHB.py index be389282..1157720a 100644 --- a/tests/em/fdem/forward/test_FDEM_forwardEJHB.py +++ b/tests/em/fdem/forward/test_FDEM_forwardEJHB.py @@ -9,7 +9,6 @@ testEB = True testHJ = True testEJ = True testBH = True -verbose = False TOLEBHJ = 1e-5 TOLEJHB = 1 # averaging and more sensitive to boundary condition violations (ie. the impact of violating the boundary conditions in each case is different.) @@ -21,56 +20,56 @@ SrcList = ['RawVec', 'MagDipole_Bfield', 'MagDipole', 'CircularLoop'] class FDEM_CrossCheck(unittest.TestCase): if testEJ: def test_EJ_CrossCheck_jxr_Jform(self): - self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'jxr', verbose=verbose, TOL=TOLEJHB)) + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'jxr', TOL=TOLEJHB)) def test_EJ_CrossCheck_jyr_Jform(self): - self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'jyr', verbose=verbose, TOL=TOLEJHB)) + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'jyr', TOL=TOLEJHB)) def test_EJ_CrossCheck_jzr_Jform(self): - self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'jzr', verbose=verbose, TOL=TOLEJHB)) + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'jzr', TOL=TOLEJHB)) def test_EJ_CrossCheck_jxi_Jform(self): - self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'jxi', verbose=verbose, TOL=TOLEJHB)) + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'jxi', TOL=TOLEJHB)) def test_EJ_CrossCheck_jyi_Jform(self): - self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'jyi', verbose=verbose, TOL=TOLEJHB)) + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'jyi', TOL=TOLEJHB)) def test_EJ_CrossCheck_jzi_Jform(self): - self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'jzi', verbose=verbose, TOL=TOLEJHB)) + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'jzi', TOL=TOLEJHB)) def test_EJ_CrossCheck_exr_Jform(self): - self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'exr', verbose=verbose, TOL=TOLEJHB)) + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'exr', TOL=TOLEJHB)) def test_EJ_CrossCheck_eyr_Jform(self): - self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'eyr', verbose=verbose, TOL=TOLEJHB)) + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'eyr', TOL=TOLEJHB)) def test_EJ_CrossCheck_ezr_Jform(self): - self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'ezr', verbose=verbose, TOL=TOLEJHB)) + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'ezr', TOL=TOLEJHB)) def test_EJ_CrossCheck_exi_Jform(self): - self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'exi', verbose=verbose, TOL=TOLEJHB)) + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'exi', TOL=TOLEJHB)) def test_EJ_CrossCheck_eyi_Jform(self): - self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'eyi', verbose=verbose, TOL=TOLEJHB)) + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'eyi', TOL=TOLEJHB)) def test_EJ_CrossCheck_ezi_Jform(self): - self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'ezi', verbose=verbose, TOL=TOLEJHB)) + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'ezi', TOL=TOLEJHB)) def test_EJ_CrossCheck_bxr_Jform(self): - self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'bxr', verbose=verbose, TOL=TOLEJHB)) + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'bxr', TOL=TOLEJHB)) def test_EJ_CrossCheck_byr_Jform(self): - self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'byr', verbose=verbose, TOL=TOLEJHB)) + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'byr', TOL=TOLEJHB)) def test_EJ_CrossCheck_bzr_Jform(self): - self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'bzr', verbose=verbose, TOL=TOLEJHB)) + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'bzr', TOL=TOLEJHB)) def test_EJ_CrossCheck_bxi_Jform(self): - self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'bxi', verbose=verbose, TOL=TOLEJHB)) + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'bxi', TOL=TOLEJHB)) def test_EJ_CrossCheck_byi_Jform(self): - self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'byi', verbose=verbose, TOL=TOLEJHB)) + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'byi', TOL=TOLEJHB)) def test_EJ_CrossCheck_bzi_Jform(self): - self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'bzi', verbose=verbose, TOL=TOLEJHB)) + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'bzi', TOL=TOLEJHB)) def test_EJ_CrossCheck_hxr_Jform(self): - self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'hxr', verbose=verbose, TOL=TOLEJHB)) + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'hxr', TOL=TOLEJHB)) def test_EJ_CrossCheck_hyr_Jform(self): - self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'hyr', verbose=verbose, TOL=TOLEJHB)) + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'hyr', TOL=TOLEJHB)) def test_EJ_CrossCheck_hzr_Jform(self): - self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'hzr', verbose=verbose, TOL=TOLEJHB)) + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'hzr', TOL=TOLEJHB)) def test_EJ_CrossCheck_hxi_Jform(self): - self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'hxi', verbose=verbose, TOL=TOLEJHB)) + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'hxi', TOL=TOLEJHB)) def test_EJ_CrossCheck_hyi_Jform(self): - self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'hyi', verbose=verbose, TOL=TOLEJHB)) + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'hyi', TOL=TOLEJHB)) def test_EJ_CrossCheck_hzi_Jform(self): - self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'hzi', verbose=verbose, TOL=TOLEJHB)) + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'hzi', TOL=TOLEJHB)) if __name__ == '__main__': unittest.main() \ No newline at end of file From ce49249664a8e83be123f41e8a7edf6eb11c3ab6 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sat, 20 Feb 2016 11:05:25 -0800 Subject: [PATCH 08/17] e,b,h,j with deriv and adjoint from e formulation --- SimPEG/EM/FDEM/FieldsFDEM.py | 134 ++++++++---------- SimPEG/EM/FDEM/SurveyFDEM.py | 2 +- .../inverse/adjoint/test_FDEM_adjointEB.py | 92 ++++++------ .../fdem/inverse/derivs/test_FDEM_derivs.py | 68 ++++----- 4 files changed, 144 insertions(+), 152 deletions(-) diff --git a/SimPEG/EM/FDEM/FieldsFDEM.py b/SimPEG/EM/FDEM/FieldsFDEM.py index fc170d46..65873188 100644 --- a/SimPEG/EM/FDEM/FieldsFDEM.py +++ b/SimPEG/EM/FDEM/FieldsFDEM.py @@ -178,7 +178,11 @@ class Fields_e(Fields): 'bPrimary' : ['eSolution','F','_bPrimary'], 'bSecondary' : ['eSolution','F','_bSecondary'], 'j' : ['eSolution','CCV','_j'], + 'jPrimary' : ['eSolution','CCV','_jPrimary'], + 'jSecondary' : ['eSolution','CCV','_jSecondary'], 'h' : ['eSolution','CCV','_h'], + 'hPrimary' : ['eSolution','CCV','_hPrimary'], + 'hSecondary' : ['eSolution','CCV','_hSecondary'], } def __init__(self, mesh, survey, **kwargs): @@ -192,6 +196,7 @@ class Fields_e(Fields): self._nC = self.survey.prob.mesh.nC self._MeSigma = self.survey.prob.MeSigma self._MeSigmaDeriv = self.survey.prob.MeSigmaDeriv + self._MfMui = self.survey.prob.MfMui def _GLoc(self, fieldType): if fieldType == 'e': @@ -203,6 +208,7 @@ class Fields_e(Fields): else: raise Exception('Field type must be e, b, h, j') + def _ePrimary(self, eSolution, srcList): """ Primary electric field from source @@ -293,9 +299,9 @@ class Fields_e(Fields): b[:,i] = b[:,i]+ 1./(1j*omega(src.freq)) * S_m return b - def _bSecondaryDeriv_u(self, src, du_dm_v, adjoint = False): + def _bDeriv_u(self, src, du_dm_v, adjoint = False): """ - Derivative of the secondary magnetic flux density with respect to the thing we solved for + Derivative of the magnetic flux density with respect to the thing we solved for :param SimPEG.EM.FDEM.Src src: source :param numpy.ndarray du_dm_v: vector to take product with @@ -310,94 +316,80 @@ class Fields_e(Fields): return - 1./(1j*omega(src.freq)) * (C * du_dm_v) - def _bSecondaryDeriv_m(self, src, v, adjoint = False): + def _bDeriv_m(self, src, v, adjoint = False): """ - Derivative of the secondary magnetic flux density with respect to the inversion model. + Derivative of the magnetic flux density with respect to the inversion model. :param SimPEG.EM.FDEM.Src src: source :param numpy.ndarray v: vector to take product with :param bool adjoint: adjoint? :rtype: numpy.ndarray - :return: product of the secondary magnetic flux density derivative with respect to the inversion model with a vector + :return: product of the magnetic flux density derivative with respect to the inversion model with a vector """ S_mDeriv, _ = src.evalDeriv(self.prob, v, adjoint) return 1./(1j * omega(src.freq)) * S_mDeriv - def _bDeriv_u(self, src, du_dm_v, adjoint=False): - """ - Partial derivative of the total magnetic flux density with respect to the thing we solved for - - :param SimPEG.EM.FDEM.Src src: source - :param numpy.ndarray du_dm_v: vector to take product with - :param bool adjoint: adjoint? - :rtype: numpy.ndarray - :return: product of the derivative of the magnetic flux density with respect to the field we solved for with a vector - """ - - return self._bSecondaryDeriv_u(src, du_dm_v, adjoint) - - def _bDeriv_m(self, src, v, adjoint=False): - """ - Partial derivative of the total magnetic flux density with respect to the inversion model. - - :param SimPEG.EM.FDEM.Src src: source - :param numpy.ndarray v: vector to take product with - :param bool adjoint: adjoint? - :rtype: SimPEG.Utils.Zero - :return: product of the magnetic flux density derivative with respect to the inversion model with a vector - """ - - # Assuming the primary does not depend on the model - return self._bSecondaryDeriv_m(src, eSolution, v, adjoint) - - def _j(self, eSolution, srcList): + def _jSecondary(self, eSolution, srcList): aveE2CCV = self._aveE2CCV - n = int(aveE2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy - Sigma = self._MeSigma - VI = sdiag(1./np.kron(np.ones(n), self.prob.mesh.vol)) - return VI * (aveE2CCV * (Sigma *eSolution) ) + n = int(aveE2CCV.shape[0] / self._nC) # number of components (instead of checking if cyl or not) + VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) + return VI * (aveE2CCV * (self._MeSigma * eSolution) ) - def _jDeriv_u(self, src, v, adjoint = False): - eSolution = self._eSecondary + def _jPrimary(self, eSolution, srcList): aveE2CCV = self._aveE2CCV - n = int(aveE2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy - Sigma = self._MeSigma + n = int(aveE2CCV.shape[0] / self._nC) # number of components (instead of checking if cyl or not) + VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) + return VI * (aveE2CCV * (self._MeSigma * self._ePrimary(eSolution, srcList)) ) - VI = sdiag(1./np.kron(np.ones(n), self.prob.mesh.vol)) + def _jDeriv_u(self, src, du_dm_v, adjoint = False): + n = int(self._aveE2CCV.shape[0] / self._nC) # number of components (instead of checking if cyl or not) + VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) - if not adjoint: - return VI * (aveE2CCV * (Sigma * (self._eDeriv_u(src, v, adjoint) ) ) ) - return self._eDeriv_u(src, Sigma.T * (aveE2CCV.T * (VI.T * v) ), adjoint) + if adjoint: + return self._eDeriv_u(src, self._MeSigma.T * (self._aveE2CCV.T * (VI.T * du_dm_v) ), adjoint) + return VI * (self._aveE2CCV * (self._MeSigma * (self._eDeriv_u(src, du_dm_v, adjoint) ) ) ) + def _jDeriv_m(self, src, v, adjoint = False): - eSolution = self._fields['eSolution'] + e = self._e(self._fields['eSolution'], [src]) + n = int(self._aveE2CCV.shape[0] / self._nC) #number of components + VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) - aveE2CCV = self._aveE2CCV - Sigma = self._MeSigma - SigmaDeriv = self._MeSigmaDeriv - e = self._e(eSolution, [src]) - - n = int(aveE2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy + if adjoint: + return self._MeSigmaDeriv(e).T * (self._aveE2CCV.T * (VI.T * v)) + self._eDeriv_m(src, self._aveE2CCV.T * (VI.T * v), adjoint=adjoint) + return VI * (self._aveE2CCV * ( self._eDeriv_m(src, v, adjoint=adjoint) + self._MeSigmaDeriv(e) * v)) - VI = sdiag(1./np.kron(np.ones(n), self.prob.mesh.vol)) - - if not adjoint: - return VI * (aveE2CCV * ( self._eDeriv_m(src, v, adjoint=adjoint) + SigmaDeriv(e) * v)) - return SigmaDeriv(aveE2CCV.T * (VI.T * e), adjoint=adjoint) * v + self._eDeriv_m(src, aveE2CCV.T * (VI.T * v), adjoint=adjoint) - def _h(self, eSolution, srcList): - b = self._b(eSolution, srcList) - Mui = self.survey.prob.MfMui - aveF2CCV = self._aveF2CCV - n = int(aveF2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy - # Mui = sdiag(sp.kron(np.ones(n), mui)) - VI = sdiag(1./np.kron(np.ones(n), self.prob.mesh.vol)) + def _hPrimary(self, eSolution, srcList): + n = int(self._aveF2CCV.shape[0] / self._nC) # Number of Components + VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) + return VI * (self._aveF2CCV * (self._MfMui * self._bPrimary(eSolution, srcList))) - return VI * (aveF2CCV * (Mui * b)) + def _hSecondary(self, eSolution, srcList): + n = int(self._aveF2CCV.shape[0] / self._nC) # Number of Components + VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) + + return VI * (self._aveF2CCV * (self._MfMui * self._bSecondary(eSolution, srcList))) + + def _hDeriv_u(self, src, du_dm_v, adjoint = False): + n = int(self._aveF2CCV.shape[0] / self._nC) # Number of Components + VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) + if adjoint: + v = self._MfMui.T * (self._aveF2CCV.T * (VI.T * du_dm_v)) + return self._bDeriv_u(src, v, adjoint=adjoint) + return VI * (self._aveF2CCV * (self._MfMui * self._bDeriv_u(src, du_dm_v, adjoint = adjoint))) + + def _hDeriv_m(self, src, v, adjoint = False): + n = int(self._aveF2CCV.shape[0] / self._nC) # Number of Components + VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) + if adjoint: + v = self._MfMui.T * (self._aveF2CCV.T * (VI.T * v)) + return self._bDeriv_m(src, v, adjoint=adjoint) + return VI * (self._aveF2CCV * (self._MfMui * self._bDeriv_m(src, v, adjoint = adjoint))) @@ -617,7 +609,7 @@ class Fields_b(Fields): n = int(aveE2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy # Sigma = sdiag(np.kron(np.ones(n), sigma)) Sigma = self.prob.MeSigma - VI = sdiag(1./np.kron(np.ones(n), self.prob.mesh.vol)) + VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) e = self._e(bSolution, srcList) @@ -629,7 +621,7 @@ class Fields_b(Fields): aveF2CCV = self._aveF2CCV n = int(aveF2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy # Mui = sdiag(sp.kron(np.ones(n), mui)) - VI = sdiag(1./np.kron(np.ones(n), self.prob.mesh.vol)) + VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) return VI * (aveF2CCV * (Mui * b)) @@ -868,7 +860,7 @@ class Fields_j(Fields): n = int(aveF2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy Rho = self.prob.MfRho - VI = sdiag(1./np.kron(np.ones(n), self.prob.mesh.vol)) + VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) j = self._j(jSolution, srcList) @@ -887,7 +879,7 @@ class Fields_j(Fields): aveE2CCV = self._aveE2CCV n = int(aveE2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy # Mu = sdiag(sp.kron(np.ones(n), mu)) - VI = sdiag(1./np.kron(np.ones(n), self.prob.mesh.vol)) + VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) return VI * (aveE2CCV * (Mu * h)) @@ -1088,7 +1080,7 @@ class Fields_h(Fields): n = int(aveF2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy Rho = self.prob.MfRho - VI = sdiag(1./np.kron(np.ones(n), self.prob.mesh.vol)) + VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) j = self._j(hSolution, srcList) @@ -1107,6 +1099,6 @@ class Fields_h(Fields): aveE2CCV = self._aveE2CCV n = int(aveE2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy # Mu = sdiag(sp.kron(np.ones(n), mu)) - VI = sdiag(1./np.kron(np.ones(n), self.prob.mesh.vol)) + VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) return VI * (aveE2CCV * (Mu * h)) diff --git a/SimPEG/EM/FDEM/SurveyFDEM.py b/SimPEG/EM/FDEM/SurveyFDEM.py index 1ba63988..111961ef 100644 --- a/SimPEG/EM/FDEM/SurveyFDEM.py +++ b/SimPEG/EM/FDEM/SurveyFDEM.py @@ -104,7 +104,7 @@ class Rx(SimPEG.Survey.BaseRx): # print self.knownRxTypes[self.rxType][:2], 'Deriv', projGLoc # projGLoc += self.knownRxTypes[self.rxType][1] - P = self.getP(mesh) + P = self.getP(mesh, self.projGLoc(u)) if not adjoint: Pv_complex = P * v diff --git a/tests/em/fdem/inverse/adjoint/test_FDEM_adjointEB.py b/tests/em/fdem/inverse/adjoint/test_FDEM_adjointEB.py index f917845a..96b6957d 100644 --- a/tests/em/fdem/inverse/adjoint/test_FDEM_adjointEB.py +++ b/tests/em/fdem/inverse/adjoint/test_FDEM_adjointEB.py @@ -46,57 +46,57 @@ def adjointTest(fdemType, comp): class FDEM_AdjointTests(unittest.TestCase): if testE: - # def test_Jtvec_adjointTest_exr_Eform(self): - # self.assertTrue(adjointTest('e', 'exr')) - # def test_Jtvec_adjointTest_eyr_Eform(self): - # self.assertTrue(adjointTest('e', 'eyr')) - # def test_Jtvec_adjointTest_ezr_Eform(self): - # self.assertTrue(adjointTest('e', 'ezr')) - # def test_Jtvec_adjointTest_exi_Eform(self): - # self.assertTrue(adjointTest('e', 'exi')) - # def test_Jtvec_adjointTest_eyi_Eform(self): - # self.assertTrue(adjointTest('e', 'eyi')) - # def test_Jtvec_adjointTest_ezi_Eform(self): - # self.assertTrue(adjointTest('e', 'ezi')) + def test_Jtvec_adjointTest_exr_Eform(self): + self.assertTrue(adjointTest('e', 'exr')) + def test_Jtvec_adjointTest_eyr_Eform(self): + self.assertTrue(adjointTest('e', 'eyr')) + def test_Jtvec_adjointTest_ezr_Eform(self): + self.assertTrue(adjointTest('e', 'ezr')) + def test_Jtvec_adjointTest_exi_Eform(self): + self.assertTrue(adjointTest('e', 'exi')) + def test_Jtvec_adjointTest_eyi_Eform(self): + self.assertTrue(adjointTest('e', 'eyi')) + def test_Jtvec_adjointTest_ezi_Eform(self): + self.assertTrue(adjointTest('e', 'ezi')) - # def test_Jtvec_adjointTest_bxr_Eform(self): - # self.assertTrue(adjointTest('e', 'bxr')) - # def test_Jtvec_adjointTest_byr_Eform(self): - # self.assertTrue(adjointTest('e', 'byr')) - # def test_Jtvec_adjointTest_bzr_Eform(self): - # self.assertTrue(adjointTest('e', 'bzr')) - # def test_Jtvec_adjointTest_bxi_Eform(self): - # self.assertTrue(adjointTest('e', 'bxi')) - # def test_Jtvec_adjointTest_byi_Eform(self): - # self.assertTrue(adjointTest('e', 'byi')) - # def test_Jtvec_adjointTest_bzi_Eform(self): - # self.assertTrue(adjointTest('e', 'bzi')) + def test_Jtvec_adjointTest_bxr_Eform(self): + self.assertTrue(adjointTest('e', 'bxr')) + def test_Jtvec_adjointTest_byr_Eform(self): + self.assertTrue(adjointTest('e', 'byr')) + def test_Jtvec_adjointTest_bzr_Eform(self): + self.assertTrue(adjointTest('e', 'bzr')) + def test_Jtvec_adjointTest_bxi_Eform(self): + self.assertTrue(adjointTest('e', 'bxi')) + def test_Jtvec_adjointTest_byi_Eform(self): + self.assertTrue(adjointTest('e', 'byi')) + def test_Jtvec_adjointTest_bzi_Eform(self): + self.assertTrue(adjointTest('e', 'bzi')) def test_Jtvec_adjointTest_exr_Eform(self): self.assertTrue(adjointTest('e', 'jxr')) - # def test_Jtvec_adjointTest_eyr_Eform(self): - # self.assertTrue(adjointTest('e', 'jyr')) - # def test_Jtvec_adjointTest_ezr_Eform(self): - # self.assertTrue(adjointTest('e', 'jzr')) - # def test_Jtvec_adjointTest_exi_Eform(self): - # self.assertTrue(adjointTest('e', 'jxi')) - # def test_Jtvec_adjointTest_eyi_Eform(self): - # self.assertTrue(adjointTest('e', 'jyi')) - # def test_Jtvec_adjointTest_ezi_Eform(self): - # self.assertTrue(adjointTest('e', 'jzi')) + def test_Jtvec_adjointTest_eyr_Eform(self): + self.assertTrue(adjointTest('e', 'jyr')) + def test_Jtvec_adjointTest_ezr_Eform(self): + self.assertTrue(adjointTest('e', 'jzr')) + def test_Jtvec_adjointTest_exi_Eform(self): + self.assertTrue(adjointTest('e', 'jxi')) + def test_Jtvec_adjointTest_eyi_Eform(self): + self.assertTrue(adjointTest('e', 'jyi')) + def test_Jtvec_adjointTest_ezi_Eform(self): + self.assertTrue(adjointTest('e', 'jzi')) - # def test_Jtvec_adjointTest_bxr_Eform(self): - # self.assertTrue(adjointTest('e', 'hxr')) - # def test_Jtvec_adjointTest_byr_Eform(self): - # self.assertTrue(adjointTest('e', 'hyr')) - # def test_Jtvec_adjointTest_bzr_Eform(self): - # self.assertTrue(adjointTest('e', 'hzr')) - # def test_Jtvec_adjointTest_bxi_Eform(self): - # self.assertTrue(adjointTest('e', 'hxi')) - # def test_Jtvec_adjointTest_byi_Eform(self): - # self.assertTrue(adjointTest('e', 'hyi')) - # def test_Jtvec_adjointTest_bzi_Eform(self): - # self.assertTrue(adjointTest('e', 'hzi')) + def test_Jtvec_adjointTest_bxr_Eform(self): + self.assertTrue(adjointTest('e', 'hxr')) + def test_Jtvec_adjointTest_byr_Eform(self): + self.assertTrue(adjointTest('e', 'hyr')) + def test_Jtvec_adjointTest_bzr_Eform(self): + self.assertTrue(adjointTest('e', 'hzr')) + def test_Jtvec_adjointTest_bxi_Eform(self): + self.assertTrue(adjointTest('e', 'hxi')) + def test_Jtvec_adjointTest_byi_Eform(self): + self.assertTrue(adjointTest('e', 'hyi')) + def test_Jtvec_adjointTest_bzi_Eform(self): + self.assertTrue(adjointTest('e', 'hzi')) if testB: def test_Jtvec_adjointTest_exr_Bform(self): diff --git a/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py b/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py index 721c2747..08fc8ae7 100644 --- a/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py +++ b/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py @@ -84,45 +84,45 @@ class FDEM_DerivTests(unittest.TestCase): def test_Jvec_ezi_Eform(self): self.assertTrue(derivTest('e', 'jzi')) - # def test_Jvec_bxr_Eform(self): - # self.assertTrue(derivTest('e', 'hxr')) - # def test_Jvec_byr_Eform(self): - # self.assertTrue(derivTest('e', 'hyr')) - # def test_Jvec_bzr_Eform(self): - # self.assertTrue(derivTest('e', 'hzr')) - # def test_Jvec_bxi_Eform(self): - # self.assertTrue(derivTest('e', 'hxi')) - # def test_Jvec_byi_Eform(self): - # self.assertTrue(derivTest('e', 'hyi')) - # def test_Jvec_bzi_Eform(self): - # self.assertTrue(derivTest('e', 'hzi')) + def test_Jvec_bxr_Eform(self): + self.assertTrue(derivTest('e', 'hxr')) + def test_Jvec_byr_Eform(self): + self.assertTrue(derivTest('e', 'hyr')) + def test_Jvec_bzr_Eform(self): + self.assertTrue(derivTest('e', 'hzr')) + def test_Jvec_bxi_Eform(self): + self.assertTrue(derivTest('e', 'hxi')) + def test_Jvec_byi_Eform(self): + self.assertTrue(derivTest('e', 'hyi')) + def test_Jvec_bzi_Eform(self): + self.assertTrue(derivTest('e', 'hzi')) if testB: def test_Jvec_exr_Bform(self): self.assertTrue(derivTest('b', 'exr')) - # def test_Jvec_eyr_Bform(self): - # self.assertTrue(derivTest('b', 'eyr')) - # def test_Jvec_ezr_Bform(self): - # self.assertTrue(derivTest('b', 'ezr')) - # def test_Jvec_exi_Bform(self): - # self.assertTrue(derivTest('b', 'exi')) - # def test_Jvec_eyi_Bform(self): - # self.assertTrue(derivTest('b', 'eyi')) - # def test_Jvec_ezi_Bform(self): - # self.assertTrue(derivTest('b', 'ezi')) + def test_Jvec_eyr_Bform(self): + self.assertTrue(derivTest('b', 'eyr')) + def test_Jvec_ezr_Bform(self): + self.assertTrue(derivTest('b', 'ezr')) + def test_Jvec_exi_Bform(self): + self.assertTrue(derivTest('b', 'exi')) + def test_Jvec_eyi_Bform(self): + self.assertTrue(derivTest('b', 'eyi')) + def test_Jvec_ezi_Bform(self): + self.assertTrue(derivTest('b', 'ezi')) - # def test_Jvec_bxr_Bform(self): - # self.assertTrue(derivTest('b', 'bxr')) - # def test_Jvec_byr_Bform(self): - # self.assertTrue(derivTest('b', 'byr')) - # def test_Jvec_bzr_Bform(self): - # self.assertTrue(derivTest('b', 'bzr')) - # def test_Jvec_bxi_Bform(self): - # self.assertTrue(derivTest('b', 'bxi')) - # def test_Jvec_byi_Bform(self): - # self.assertTrue(derivTest('b', 'byi')) - # def test_Jvec_bzi_Bform(self): - # self.assertTrue(derivTest('b', 'bzi')) + def test_Jvec_bxr_Bform(self): + self.assertTrue(derivTest('b', 'bxr')) + def test_Jvec_byr_Bform(self): + self.assertTrue(derivTest('b', 'byr')) + def test_Jvec_bzr_Bform(self): + self.assertTrue(derivTest('b', 'bzr')) + def test_Jvec_bxi_Bform(self): + self.assertTrue(derivTest('b', 'bxi')) + def test_Jvec_byi_Bform(self): + self.assertTrue(derivTest('b', 'byi')) + def test_Jvec_bzi_Bform(self): + self.assertTrue(derivTest('b', 'bzi')) if testHJ: def test_Jvec_jxr_Jform(self): From 0edbc9f6ca408856b4e1bab125751cbe50d2ffef Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sat, 20 Feb 2016 11:13:09 -0800 Subject: [PATCH 09/17] docs for fields_e --- SimPEG/EM/FDEM/FieldsFDEM.py | 82 +++++++++++++++++++++++++++++++++--- 1 file changed, 75 insertions(+), 7 deletions(-) diff --git a/SimPEG/EM/FDEM/FieldsFDEM.py b/SimPEG/EM/FDEM/FieldsFDEM.py index 65873188..267fbb9f 100644 --- a/SimPEG/EM/FDEM/FieldsFDEM.py +++ b/SimPEG/EM/FDEM/FieldsFDEM.py @@ -307,7 +307,7 @@ class Fields_e(Fields): :param numpy.ndarray du_dm_v: vector to take product with :param bool adjoint: adjoint? :rtype: numpy.ndarray - :return: product of the derivative of the secondary magnetic flux density with respect to the field we solved for with a vector + :return: product of the derivative of the magnetic flux density with respect to the field we solved for with a vector """ C = self._edgeCurl @@ -331,19 +331,44 @@ class Fields_e(Fields): return 1./(1j * omega(src.freq)) * S_mDeriv - def _jSecondary(self, eSolution, srcList): - aveE2CCV = self._aveE2CCV - n = int(aveE2CCV.shape[0] / self._nC) # number of components (instead of checking if cyl or not) - VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) - return VI * (aveE2CCV * (self._MeSigma * eSolution) ) - def _jPrimary(self, eSolution, srcList): + """ + Primary current density + + :param numpy.ndarray eSolution: field we solved for + :param list srcList: list of sources + :rtype: numpy.ndarray + :return: primary current density + """ aveE2CCV = self._aveE2CCV n = int(aveE2CCV.shape[0] / self._nC) # number of components (instead of checking if cyl or not) VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) return VI * (aveE2CCV * (self._MeSigma * self._ePrimary(eSolution, srcList)) ) + def _jSecondary(self, eSolution, srcList): + """ + Secondary current density from eSolution + + :param numpy.ndarray eSolution: field we solved for + :param list srcList: list of sources + :rtype: numpy.ndarray + :return: secondary current density + """ + aveE2CCV = self._aveE2CCV + n = int(aveE2CCV.shape[0] / self._nC) # number of components (instead of checking if cyl or not) + VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) + return VI * (aveE2CCV * (self._MeSigma * eSolution) ) + def _jDeriv_u(self, src, du_dm_v, adjoint = False): + """ + Derivative of the current density with respect to the thing we solved for + + :param SimPEG.EM.FDEM.Src src: source + :param numpy.ndarray du_dm_v: vector to take product with + :param bool adjoint: adjoint? + :rtype: numpy.ndarray + :return: product of the derivative of the current density with respect to the field we solved for with a vector + """ n = int(self._aveE2CCV.shape[0] / self._nC) # number of components (instead of checking if cyl or not) VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) @@ -353,6 +378,15 @@ class Fields_e(Fields): def _jDeriv_m(self, src, v, adjoint = False): + """ + Derivative of the current density with respect to the inversion model. + + :param SimPEG.EM.FDEM.Src src: source + :param numpy.ndarray v: vector to take product with + :param bool adjoint: adjoint? + :rtype: numpy.ndarray + :return: product of the current density derivative with respect to the inversion model with a vector + """ e = self._e(self._fields['eSolution'], [src]) n = int(self._aveE2CCV.shape[0] / self._nC) #number of components VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) @@ -364,18 +398,43 @@ class Fields_e(Fields): def _hPrimary(self, eSolution, srcList): + """ + Primary magnetic field + + :param numpy.ndarray eSolution: field we solved for + :param list srcList: list of sources + :rtype: numpy.ndarray + :return: primary magnetic field + """ n = int(self._aveF2CCV.shape[0] / self._nC) # Number of Components VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) return VI * (self._aveF2CCV * (self._MfMui * self._bPrimary(eSolution, srcList))) def _hSecondary(self, eSolution, srcList): + """ + Secondary magnetic field from eSolution + + :param numpy.ndarray eSolution: field we solved for + :param list srcList: list of sources + :rtype: numpy.ndarray + :return: secondary magnetic field + """ n = int(self._aveF2CCV.shape[0] / self._nC) # Number of Components VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) return VI * (self._aveF2CCV * (self._MfMui * self._bSecondary(eSolution, srcList))) def _hDeriv_u(self, src, du_dm_v, adjoint = False): + """ + Derivative of the magnetic field with respect to the thing we solved for + + :param SimPEG.EM.FDEM.Src src: source + :param numpy.ndarray du_dm_v: vector to take product with + :param bool adjoint: adjoint? + :rtype: numpy.ndarray + :return: product of the derivative of the magnetic field with respect to the field we solved for with a vector + """ n = int(self._aveF2CCV.shape[0] / self._nC) # Number of Components VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) if adjoint: @@ -384,6 +443,15 @@ class Fields_e(Fields): return VI * (self._aveF2CCV * (self._MfMui * self._bDeriv_u(src, du_dm_v, adjoint = adjoint))) def _hDeriv_m(self, src, v, adjoint = False): + """ + Derivative of the magnetic field with respect to the inversion model. + + :param SimPEG.EM.FDEM.Src src: source + :param numpy.ndarray v: vector to take product with + :param bool adjoint: adjoint? + :rtype: numpy.ndarray + :return: product of the magnetic field derivative with respect to the inversion model with a vector + """ n = int(self._aveF2CCV.shape[0] / self._nC) # Number of Components VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) if adjoint: From 0646a930ab2ab642b80162ac24ecd9c0f59c938f Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sat, 20 Feb 2016 13:27:51 -0800 Subject: [PATCH 10/17] e,b,h,j from b formulation --- SimPEG/EM/FDEM/FieldsFDEM.py | 179 ++++++++++++------ tests/em/fdem/forward/test_FDEM_forward.py | 46 ++--- .../em/fdem/forward/test_FDEM_forwardEJHB.py | 8 +- tests/em/fdem/forward/test_FDEM_forwardHB.py | 106 ++++++++--- .../inverse/adjoint/test_FDEM_adjointEB.py | 52 +++-- .../fdem/inverse/derivs/test_FDEM_derivs.py | 76 +++++--- 6 files changed, 317 insertions(+), 150 deletions(-) diff --git a/SimPEG/EM/FDEM/FieldsFDEM.py b/SimPEG/EM/FDEM/FieldsFDEM.py index 267fbb9f..b8971611 100644 --- a/SimPEG/EM/FDEM/FieldsFDEM.py +++ b/SimPEG/EM/FDEM/FieldsFDEM.py @@ -234,8 +234,7 @@ class Fields_e(Fields): :rtype: numpy.ndarray :return: secondary electric field """ - ind = self.prob.survey.getSourceIndex(srcList) - return eSolution[:,ind] + return eSolution def _eDeriv_u(self, src, v, adjoint = False): """ @@ -373,8 +372,8 @@ class Fields_e(Fields): VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) if adjoint: - return self._eDeriv_u(src, self._MeSigma.T * (self._aveE2CCV.T * (VI.T * du_dm_v) ), adjoint) - return VI * (self._aveE2CCV * (self._MeSigma * (self._eDeriv_u(src, du_dm_v, adjoint) ) ) ) + return self._eDeriv_u(src, self._MeSigma.T * (self._aveE2CCV.T * (VI.T * du_dm_v) ), adjoint = adjoint) + return VI * (self._aveE2CCV * (self._MeSigma * (self._eDeriv_u(src, du_dm_v, adjoint=adjoint) ) ) ) def _jDeriv_m(self, src, v, adjoint = False): @@ -387,7 +386,7 @@ class Fields_e(Fields): :rtype: numpy.ndarray :return: product of the current density derivative with respect to the inversion model with a vector """ - e = self._e(self._fields['eSolution'], [src]) + e = self[src, 'e'] n = int(self._aveE2CCV.shape[0] / self._nC) #number of components VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) @@ -477,8 +476,12 @@ class Fields_b(Fields): 'e' : ['bSolution','E','_e'], 'ePrimary' : ['bSolution','E','_ePrimary'], 'eSecondary' : ['bSolution','E','_eSecondary'], - 'j' : ['bSolution','C','_j'], - 'h' : ['bSolution','C','_h'], + 'j' : ['bSolution','CCV','_j'], + 'jPrimary' : ['bSolution','CCV','_jPrimary'], + 'jSecondary' : ['bSolution','CCV','_jSecondary'], + 'h' : ['bSolution','CCV','_h'], + 'hPrimary' : ['bSolution','CCV','_hPrimary'], + 'hSecondary' : ['bSolution','CCV','_hSecondary'], } def __init__(self,mesh,survey,**kwargs): @@ -519,7 +522,7 @@ class Fields_b(Fields): :return: primary electric field as defined by the sources """ - bPrimary = np.zeros_like(bSolution) + bPrimary = np.zeros([self.prob.mesh.nF,len(srcList)]) for i, src in enumerate(srcList): bp = src.bPrimary(self.prob) bPrimary[:,i] = bPrimary[:,i] + bp @@ -594,105 +597,165 @@ class Fields_b(Fields): e = self._MeSigmaI * ( self._edgeCurl.T * ( self._MfMui * bSolution)) for i,src in enumerate(srcList): _,S_e = src.eval(self.prob) - e[:,i] = e[:,i]+ -self._MeSigmaI * S_e + e[:,i] = e[:,i] + -self._MeSigmaI * S_e return e - def _eSecondaryDeriv_u(self, src, du_dm_v, adjoint=False): + def _eDeriv_u(self, src, du_dm_v, adjoint=False): """ - Derivative of the secondary electric field with respect to the thing we solved for + Derivative of the electric field with respect to the thing we solved for :param SimPEG.EM.FDEM.Src src: source :param numpy.ndarray v: vector to take product with :param bool adjoint: adjoint? :rtype: numpy.ndarray - :return: product of the derivative of the secondary electric field with respect to the field we solved for with a vector + :return: product of the derivative of the electric field with respect to the field we solved for with a vector """ if not adjoint: return self._MeSigmaI * ( self._edgeCurl.T * ( self._MfMui * du_dm_v) ) - else: - return self._MfMui.T * (self._edgeCurl * (self._MeSigmaI.T * du_dm_v)) + return self._MfMui.T * (self._edgeCurl * (self._MeSigmaI.T * du_dm_v)) - def _eSecondaryDeriv_m(self, src, v, adjoint=False): + def _eDeriv_m(self, src, v, adjoint=False): """ - Derivative of the secondary electric field with respect to the inversion model + Derivative of the electric field with respect to the inversion model :param SimPEG.EM.FDEM.Src src: source :param numpy.ndarray v: vector to take product with :param bool adjoint: adjoint? :rtype: numpy.ndarray - :return: product of the derivative of the secondary electric field with respect to the model with a vector + :return: product of the derivative of the electric field with respect to the model with a vector """ - bSolution = self[[src],'bSolution'] + bSolution = Utils.mkvc(self[src, 'bSolution']) _,S_e = src.eval(self.prob) - Me = self._Me - - if adjoint: - Me = Me.T - - w = self._edgeCurl.T * (self._MfMui * bSolution) - w = w - Utils.mkvc(Me * S_e,2) - - if not adjoint: - de_dm = self._MeSigmaIDeriv(w) * v - elif adjoint: - de_dm = self._MeSigmaIDeriv(w).T * v + w = -S_e + self._edgeCurl.T * (self._MfMui * bSolution) _, S_eDeriv = src.evalDeriv(self.prob, v, adjoint) - de_dm = de_dm - self._MeSigmaI * S_eDeriv - return de_dm + if adjoint: + return self._MeSigmaIDeriv(w).T * v - self._MeSigmaI.T * S_eDeriv + return self._MeSigmaIDeriv(w) * v - self._MeSigmaI *S_eDeriv - def _eDeriv_u(self, src, du_dm_v, adjoint=False): + + def _jPrimary(self, bSolution, srcList): """ - Partial derivative of the total electric field with respect to the thing we solved for + Primary current density + + :param numpy.ndarray bSolution: field we solved for + :param list srcList: list of sources + :rtype: numpy.ndarray + :return: primary current density + """ + n = int(self._aveE2CCV.shape[0] / self._nC) # number of components + VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) + + return VI * (self._aveE2CCV * (self.prob.MeSigma * self._ePrimary(bSolution, srcList))) + + def _jSecondary(self, bSolution, srcList): + """ + Secondary current density from bSolution + + :param numpy.ndarray bSolution: field we solved for + :param list srcList: list of sources + :rtype: numpy.ndarray + :return: primary current density + """ + + n = int(self._aveE2CCV.shape[0] / self._nC) # number of components + VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) + + return VI * (self._aveE2CCV * ( self._edgeCurl.T * ( self._MfMui * bSolution) ) ) + + def _jDeriv_u(self, src, du_dm_v, adjoint=False): + """ + Partial derivative of the current density with respect to the thing we + solved for. :param SimPEG.EM.FDEM.Src src: source :param numpy.ndarray du_dm_v: vector to take product with :param bool adjoint: adjoint? :rtype: numpy.ndarray - :return: product of the derivative of the electric field with respect to the field we solved for with a vector + :return: product of the derivative of the current density with respect to the field we solved for with a vector """ - return self._eSecondaryDeriv_u(src, du_dm_v, adjoint) - def _eDeriv_m(self, src, v, adjoint=False): + n = int(self._aveE2CCV.shape[0] / self._nC) # number of components + VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) + + if adjoint: + return self._MfMui.T * ( self._edgeCurl * ( self._aveE2CCV.T * ( VI.T * du_dm_v ) ) ) + return VI * ( self._aveE2CCV * ( self._edgeCurl.T * ( self._MfMui * du_dm_v ) ) ) + + def _jDeriv_m(self, src, v, adjoint=False): """ - Partial derivative of the total electric field density with respect to the inversion model. + Derivative of the current density with respect to the inversion model :param SimPEG.EM.FDEM.Src src: source :param numpy.ndarray v: vector to take product with :param bool adjoint: adjoint? :rtype: numpy.ndarray - :return: product of the electric field derivative with respect to the inversion model with a vector + :return: product of the derivative of the current density with respect to the model with a vector """ - # assuming primary doesn't depend on model - return self._eSecondaryDeriv_m(src, bSolution, v, adjoint) + return Zero() - def _j(self, bSolution, srcList): - sigma = self._sigma - aveE2CCV = self._aveE2CCV - n = int(aveE2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy - # Sigma = sdiag(np.kron(np.ones(n), sigma)) - Sigma = self.prob.MeSigma + def _hPrimary(self, bSolution, srcList): + """ + Primary magnetic field + + :param numpy.ndarray bSolution: field we solved for + :param list srcList: list of sources + :rtype: numpy.ndarray + :return: primary magnetic field + """ + n = int(self._aveF2CCV.shape[0] / self._nC) #number of components + VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) + return VI * (self._aveF2CCV * (self._MfMui * self._bPrimary(bSolution, srcList))) + + def _hSecondary(self, bSolution, srcList): + """ + Secondary magnetic field from bSolution + + :param numpy.ndarray bSolution: field we solved for + :param list srcList: list of sources + :rtype: numpy.ndarray + :return: secondary magnetic field + """ + n = int(self._aveF2CCV.shape[0] / self._nC) #number of components + VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) + return VI * (self._aveF2CCV * (self._MfMui * self._bSecondary(bSolution, srcList))) + + def _hDeriv_u(self, src, du_dm_v, adjoint=False): + """ + Partial derivative of the magnetic field with respect to the thing we + solved for. + + :param SimPEG.EM.FDEM.Src src: source + :param numpy.ndarray du_dm_v: vector to take product with + :param bool adjoint: adjoint? + :rtype: numpy.ndarray + :return: product of the derivative of the magnetic field with respect to the field we solved for with a vector + """ + n = int(self._aveF2CCV.shape[0] / self._nC) #number of components VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) - e = self._e(bSolution, srcList) + if adjoint: + return self._MfMui.T * ( self._aveF2CCV.T * ( VI.T * du_dm_v) ) + return VI * (self._aveF2CCV * (self._MfMui * du_dm_v)) - return VI * (aveE2CCV * (Sigma *e) ) - - def _h(self, bSolution, srcList): - b = self._b(bSolution, srcList) - Mui = self.survey.prob.MfMui - aveF2CCV = self._aveF2CCV - n = int(aveF2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy - # Mui = sdiag(sp.kron(np.ones(n), mui)) - VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) + def _hDeriv_m(self, src, v, adjoint=False): + """ + Derivative of the magnetic field with respect to the inversion model + + :param SimPEG.EM.FDEM.Src src: source + :param numpy.ndarray v: vector to take product with + :param bool adjoint: adjoint? + :rtype: numpy.ndarray + :return: product of the derivative of the magnetic field with respect to the model with a vector + """ + return Zero() - return VI * (aveF2CCV * (Mui * b)) class Fields_j(Fields): diff --git a/tests/em/fdem/forward/test_FDEM_forward.py b/tests/em/fdem/forward/test_FDEM_forward.py index da446675..d89612d8 100644 --- a/tests/em/fdem/forward/test_FDEM_forward.py +++ b/tests/em/fdem/forward/test_FDEM_forward.py @@ -6,7 +6,7 @@ from scipy.constants import mu_0 from SimPEG.EM.Utils.testingUtils import getFDEMProblem, crossCheckTest testEB = True -testHJ = True +testHJ = False testEJ = True testBH = True verbose = False @@ -22,29 +22,29 @@ class FDEM_CrossCheck(unittest.TestCase): if testEB: def test_EB_CrossCheck_exr_Eform(self): self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'exr', verbose=verbose)) - def test_EB_CrossCheck_eyr_Eform(self): - self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'eyr', verbose=verbose)) - def test_EB_CrossCheck_ezr_Eform(self): - self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'ezr', verbose=verbose)) - def test_EB_CrossCheck_exi_Eform(self): - self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'exi', verbose=verbose)) - def test_EB_CrossCheck_eyi_Eform(self): - self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'eyi', verbose=verbose)) - def test_EB_CrossCheck_ezi_Eform(self): - self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'ezi', verbose=verbose)) + # def test_EB_CrossCheck_eyr_Eform(self): + # self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'eyr', verbose=verbose)) + # def test_EB_CrossCheck_ezr_Eform(self): + # self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'ezr', verbose=verbose)) + # def test_EB_CrossCheck_exi_Eform(self): + # self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'exi', verbose=verbose)) + # def test_EB_CrossCheck_eyi_Eform(self): + # self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'eyi', verbose=verbose)) + # def test_EB_CrossCheck_ezi_Eform(self): + # self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'ezi', verbose=verbose)) - def test_EB_CrossCheck_bxr_Eform(self): - self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'bxr', verbose=verbose)) - def test_EB_CrossCheck_byr_Eform(self): - self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'byr', verbose=verbose)) - def test_EB_CrossCheck_bzr_Eform(self): - self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'bzr', verbose=verbose)) - def test_EB_CrossCheck_bxi_Eform(self): - self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'bxi', verbose=verbose)) - def test_EB_CrossCheck_byi_Eform(self): - self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'byi', verbose=verbose)) - def test_EB_CrossCheck_bzi_Eform(self): - self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'bzi', verbose=verbose)) + # def test_EB_CrossCheck_bxr_Eform(self): + # self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'bxr', verbose=verbose)) + # def test_EB_CrossCheck_byr_Eform(self): + # self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'byr', verbose=verbose)) + # def test_EB_CrossCheck_bzr_Eform(self): + # self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'bzr', verbose=verbose)) + # def test_EB_CrossCheck_bxi_Eform(self): + # self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'bxi', verbose=verbose)) + # def test_EB_CrossCheck_byi_Eform(self): + # self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'byi', verbose=verbose)) + # def test_EB_CrossCheck_bzi_Eform(self): + # self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'bzi', verbose=verbose)) if testHJ: def test_HJ_CrossCheck_jxr_Jform(self): diff --git a/tests/em/fdem/forward/test_FDEM_forwardEJHB.py b/tests/em/fdem/forward/test_FDEM_forwardEJHB.py index 1157720a..c01e70db 100644 --- a/tests/em/fdem/forward/test_FDEM_forwardEJHB.py +++ b/tests/em/fdem/forward/test_FDEM_forwardEJHB.py @@ -5,10 +5,10 @@ import sys from scipy.constants import mu_0 from SimPEG.EM.Utils.testingUtils import getFDEMProblem, crossCheckTest -testEB = True -testHJ = True -testEJ = True -testBH = True +testEB = False +testHJ = False +testEJ = False +testBH = False TOLEBHJ = 1e-5 TOLEJHB = 1 # averaging and more sensitive to boundary condition violations (ie. the impact of violating the boundary conditions in each case is different.) diff --git a/tests/em/fdem/forward/test_FDEM_forwardHB.py b/tests/em/fdem/forward/test_FDEM_forwardHB.py index 78b3e376..1f8be91c 100644 --- a/tests/em/fdem/forward/test_FDEM_forwardHB.py +++ b/tests/em/fdem/forward/test_FDEM_forwardHB.py @@ -5,13 +5,12 @@ import sys from scipy.constants import mu_0 from SimPEG.EM.Utils.testingUtils import getFDEMProblem, crossCheckTest -testEB = True +testEB = False testHJ = True -testEJ = True +testEJ = False testBH = True verbose = False -TOLEBHJ = 1e-5 TOLEJHB = 1 # averaging and more sensitive to boundary condition violations (ie. the impact of violating the boundary conditions in each case is different.) #TODO: choose better testing parameters to lower this @@ -20,56 +19,109 @@ SrcList = ['RawVec', 'MagDipole_Bfield', 'MagDipole', 'CircularLoop'] class FDEM_CrossCheck(unittest.TestCase): if testBH: - def test_BH_CrossCheck_jxr_Jform(self): + def test_BH_CrossCheck_jxr(self): self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jxr', verbose=verbose, TOL=TOLEJHB)) - def test_BH_CrossCheck_jyr_Jform(self): + def test_BH_CrossCheck_jyr(self): self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jyr', verbose=verbose, TOL=TOLEJHB)) - def test_BH_CrossCheck_jzr_Jform(self): + def test_BH_CrossCheck_jzr(self): self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jzr', verbose=verbose, TOL=TOLEJHB)) - def test_BH_CrossCheck_jxi_Jform(self): + def test_BH_CrossCheck_jxi(self): self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jxi', verbose=verbose, TOL=TOLEJHB)) - def test_BH_CrossCheck_jyi_Jform(self): + def test_BH_CrossCheck_jyi(self): self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jyi', verbose=verbose, TOL=TOLEJHB)) - def test_BH_CrossCheck_jzi_Jform(self): + def test_BH_CrossCheck_jzi(self): self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jzi', verbose=verbose, TOL=TOLEJHB)) - def test_BH_CrossCheck_exr_Jform(self): + def test_BH_CrossCheck_exr(self): self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'exr', verbose=verbose, TOL=TOLEJHB)) - def test_BH_CrossCheck_eyr_Jform(self): + def test_BH_CrossCheck_eyr(self): self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'eyr', verbose=verbose, TOL=TOLEJHB)) - def test_BH_CrossCheck_ezr_Jform(self): + def test_BH_CrossCheck_ezr(self): self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'ezr', verbose=verbose, TOL=TOLEJHB)) - def test_BH_CrossCheck_exi_Jform(self): + def test_BH_CrossCheck_exi(self): self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'exi', verbose=verbose, TOL=TOLEJHB)) - def test_BH_CrossCheck_eyi_Jform(self): + def test_BH_CrossCheck_eyi(self): self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'eyi', verbose=verbose, TOL=TOLEJHB)) - def test_BH_CrossCheck_ezi_Jform(self): + def test_BH_CrossCheck_ezi(self): self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'ezi', verbose=verbose, TOL=TOLEJHB)) - def test_BH_CrossCheck_bxr_Jform(self): + def test_BH_CrossCheck_bxr(self): self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'bxr', verbose=verbose, TOL=TOLEJHB)) - def test_BH_CrossCheck_byr_Jform(self): + def test_BH_CrossCheck_byr(self): self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'byr', verbose=verbose, TOL=TOLEJHB)) - def test_BH_CrossCheck_bzr_Jform(self): + def test_BH_CrossCheck_bzr(self): self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'bzr', verbose=verbose, TOL=TOLEJHB)) - def test_BH_CrossCheck_bxi_Jform(self): + def test_BH_CrossCheck_bxi(self): self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'bxi', verbose=verbose, TOL=TOLEJHB)) - def test_BH_CrossCheck_byi_Jform(self): + def test_BH_CrossCheck_byi(self): self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'byi', verbose=verbose, TOL=TOLEJHB)) - def test_BH_CrossCheck_bzi_Jform(self): + def test_BH_CrossCheck_bzi(self): self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'bzi', verbose=verbose, TOL=TOLEJHB)) - def test_BH_CrossCheck_hxr_Jform(self): + def test_BH_CrossCheck_hxr(self): self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hxr', verbose=verbose, TOL=TOLEJHB)) - def test_BH_CrossCheck_hyr_Jform(self): + def test_BH_CrossCheck_hyr(self): self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hyr', verbose=verbose, TOL=TOLEJHB)) - def test_BH_CrossCheck_hzr_Jform(self): + def test_BH_CrossCheck_hzr(self): self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hzr', verbose=verbose, TOL=TOLEJHB)) - def test_BH_CrossCheck_hxi_Jform(self): + def test_BH_CrossCheck_hxi(self): self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hxi', verbose=verbose, TOL=TOLEJHB)) - def test_BH_CrossCheck_hyi_Jform(self): + def test_BH_CrossCheck_hyi(self): self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hyi', verbose=verbose, TOL=TOLEJHB)) - def test_BH_CrossCheck_hzi_Jform(self): + def test_BH_CrossCheck_hzi(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hzi', verbose=verbose, TOL=TOLEJHB)) + + if testBH: + def test_BH_CrossCheck_jxr(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jxr', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_jyr(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jyr', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_jzr(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jzr', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_jxi(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jxi', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_jyi(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jyi', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_jzi(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jzi', verbose=verbose, TOL=TOLEJHB)) + + def test_BH_CrossCheck_exr(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'exr', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_eyr(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'eyr', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_ezr(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'ezr', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_exi(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'exi', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_eyi(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'eyi', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_ezi(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'ezi', verbose=verbose, TOL=TOLEJHB)) + + def test_BH_CrossCheck_bxr(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'bxr', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_byr(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'byr', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_bzr(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'bzr', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_bxi(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'bxi', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_byi(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'byi', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_bzi(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'bzi', verbose=verbose, TOL=TOLEJHB)) + + def test_BH_CrossCheck_hxr(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hxr', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_hyr(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hyr', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_hzr(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hzr', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_hxi(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hxi', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_hyi(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hyi', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_hzi(self): self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hzi', verbose=verbose, TOL=TOLEJHB)) if __name__ == '__main__': diff --git a/tests/em/fdem/inverse/adjoint/test_FDEM_adjointEB.py b/tests/em/fdem/inverse/adjoint/test_FDEM_adjointEB.py index 96b6957d..f8d2dc9d 100644 --- a/tests/em/fdem/inverse/adjoint/test_FDEM_adjointEB.py +++ b/tests/em/fdem/inverse/adjoint/test_FDEM_adjointEB.py @@ -6,7 +6,7 @@ from scipy.constants import mu_0 from SimPEG.EM.Utils.testingUtils import getFDEMProblem testE = True -testB = False +testB = True verbose = False @@ -72,30 +72,30 @@ class FDEM_AdjointTests(unittest.TestCase): def test_Jtvec_adjointTest_bzi_Eform(self): self.assertTrue(adjointTest('e', 'bzi')) - def test_Jtvec_adjointTest_exr_Eform(self): + def test_Jtvec_adjointTest_jxr_Eform(self): self.assertTrue(adjointTest('e', 'jxr')) - def test_Jtvec_adjointTest_eyr_Eform(self): + def test_Jtvec_adjointTest_jyr_Eform(self): self.assertTrue(adjointTest('e', 'jyr')) - def test_Jtvec_adjointTest_ezr_Eform(self): + def test_Jtvec_adjointTest_jzr_Eform(self): self.assertTrue(adjointTest('e', 'jzr')) - def test_Jtvec_adjointTest_exi_Eform(self): + def test_Jtvec_adjointTest_jxi_Eform(self): self.assertTrue(adjointTest('e', 'jxi')) - def test_Jtvec_adjointTest_eyi_Eform(self): + def test_Jtvec_adjointTest_jyi_Eform(self): self.assertTrue(adjointTest('e', 'jyi')) - def test_Jtvec_adjointTest_ezi_Eform(self): + def test_Jtvec_adjointTest_jzi_Eform(self): self.assertTrue(adjointTest('e', 'jzi')) - def test_Jtvec_adjointTest_bxr_Eform(self): + def test_Jtvec_adjointTest_hxr_Eform(self): self.assertTrue(adjointTest('e', 'hxr')) - def test_Jtvec_adjointTest_byr_Eform(self): + def test_Jtvec_adjointTest_hyr_Eform(self): self.assertTrue(adjointTest('e', 'hyr')) - def test_Jtvec_adjointTest_bzr_Eform(self): + def test_Jtvec_adjointTest_hzr_Eform(self): self.assertTrue(adjointTest('e', 'hzr')) - def test_Jtvec_adjointTest_bxi_Eform(self): + def test_Jtvec_adjointTest_hxi_Eform(self): self.assertTrue(adjointTest('e', 'hxi')) - def test_Jtvec_adjointTest_byi_Eform(self): + def test_Jtvec_adjointTest_hyi_Eform(self): self.assertTrue(adjointTest('e', 'hyi')) - def test_Jtvec_adjointTest_bzi_Eform(self): + def test_Jtvec_adjointTest_hzi_Eform(self): self.assertTrue(adjointTest('e', 'hzi')) if testB: @@ -125,6 +125,32 @@ class FDEM_AdjointTests(unittest.TestCase): def test_Jtvec_adjointTest_bzi_Bform(self): self.assertTrue(adjointTest('b', 'bzi')) + def test_Jtvec_adjointTest_jxr_Bform(self): + self.assertTrue(adjointTest('b', 'jxr')) + def test_Jtvec_adjointTest_jyr_Bform(self): + self.assertTrue(adjointTest('b', 'jyr')) + def test_Jtvec_adjointTest_jzr_Bform(self): + self.assertTrue(adjointTest('b', 'jzr')) + def test_Jtvec_adjointTest_jxi_Bform(self): + self.assertTrue(adjointTest('b', 'jxi')) + def test_Jtvec_adjointTest_jyi_Bform(self): + self.assertTrue(adjointTest('b', 'jyi')) + def test_Jtvec_adjointTest_jzi_Bform(self): + self.assertTrue(adjointTest('b', 'jzi')) + + def test_Jtvec_adjointTest_hxr_Bform(self): + self.assertTrue(adjointTest('b', 'hxr')) + def test_Jtvec_adjointTest_hyr_Bform(self): + self.assertTrue(adjointTest('b', 'hyr')) + def test_Jtvec_adjointTest_hzr_Bform(self): + self.assertTrue(adjointTest('b', 'hzr')) + def test_Jtvec_adjointTest_hxi_Bform(self): + self.assertTrue(adjointTest('b', 'hxi')) + def test_Jtvec_adjointTest_hyi_Bform(self): + self.assertTrue(adjointTest('b', 'hyi')) + def test_Jtvec_adjointTest_hzi_Bform(self): + self.assertTrue(adjointTest('b', 'hzi')) + if __name__ == '__main__': unittest.main() diff --git a/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py b/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py index 08fc8ae7..19f3e9e6 100644 --- a/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py +++ b/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py @@ -7,7 +7,7 @@ from SimPEG.EM.Utils.testingUtils import getFDEMProblem testDerivs = True testE = True -testB = False +testB = True testHJ = False verbose = False @@ -98,31 +98,57 @@ class FDEM_DerivTests(unittest.TestCase): self.assertTrue(derivTest('e', 'hzi')) if testB: - def test_Jvec_exr_Bform(self): - self.assertTrue(derivTest('b', 'exr')) - def test_Jvec_eyr_Bform(self): - self.assertTrue(derivTest('b', 'eyr')) - def test_Jvec_ezr_Bform(self): - self.assertTrue(derivTest('b', 'ezr')) - def test_Jvec_exi_Bform(self): - self.assertTrue(derivTest('b', 'exi')) - def test_Jvec_eyi_Bform(self): - self.assertTrue(derivTest('b', 'eyi')) - def test_Jvec_ezi_Bform(self): - self.assertTrue(derivTest('b', 'ezi')) + # def test_Jvec_exr_Bform(self): + # self.assertTrue(derivTest('b', 'exr')) + # def test_Jvec_eyr_Bform(self): + # self.assertTrue(derivTest('b', 'eyr')) + # def test_Jvec_ezr_Bform(self): + # self.assertTrue(derivTest('b', 'ezr')) + # def test_Jvec_exi_Bform(self): + # self.assertTrue(derivTest('b', 'exi')) + # def test_Jvec_eyi_Bform(self): + # self.assertTrue(derivTest('b', 'eyi')) + # def test_Jvec_ezi_Bform(self): + # self.assertTrue(derivTest('b', 'ezi')) - def test_Jvec_bxr_Bform(self): - self.assertTrue(derivTest('b', 'bxr')) - def test_Jvec_byr_Bform(self): - self.assertTrue(derivTest('b', 'byr')) - def test_Jvec_bzr_Bform(self): - self.assertTrue(derivTest('b', 'bzr')) - def test_Jvec_bxi_Bform(self): - self.assertTrue(derivTest('b', 'bxi')) - def test_Jvec_byi_Bform(self): - self.assertTrue(derivTest('b', 'byi')) - def test_Jvec_bzi_Bform(self): - self.assertTrue(derivTest('b', 'bzi')) + # def test_Jvec_bxr_Bform(self): + # self.assertTrue(derivTest('b', 'bxr')) + # def test_Jvec_byr_Bform(self): + # self.assertTrue(derivTest('b', 'byr')) + # def test_Jvec_bzr_Bform(self): + # self.assertTrue(derivTest('b', 'bzr')) + # def test_Jvec_bxi_Bform(self): + # self.assertTrue(derivTest('b', 'bxi')) + # def test_Jvec_byi_Bform(self): + # self.assertTrue(derivTest('b', 'byi')) + # def test_Jvec_bzi_Bform(self): + # self.assertTrue(derivTest('b', 'bzi')) + + def test_Jvec_jxr_Bform(self): + self.assertTrue(derivTest('b', 'jxr')) + def test_Jvec_jyr_Bform(self): + self.assertTrue(derivTest('b', 'jyr')) + def test_Jvec_jzr_Bform(self): + self.assertTrue(derivTest('b', 'jzr')) + def test_Jvec_jxi_Bform(self): + self.assertTrue(derivTest('b', 'jxi')) + def test_Jvec_jyi_Bform(self): + self.assertTrue(derivTest('b', 'jyi')) + def test_Jvec_jzi_Bform(self): + self.assertTrue(derivTest('b', 'jzi')) + + def test_Jvec_hxr_Bform(self): + self.assertTrue(derivTest('b', 'hxr')) + def test_Jvec_hyr_Bform(self): + self.assertTrue(derivTest('b', 'hyr')) + def test_Jvec_hzr_Bform(self): + self.assertTrue(derivTest('b', 'hzr')) + def test_Jvec_hxi_Bform(self): + self.assertTrue(derivTest('b', 'hxi')) + def test_Jvec_hyi_Bform(self): + self.assertTrue(derivTest('b', 'hyi')) + def test_Jvec_hzi_Bform(self): + self.assertTrue(derivTest('b', 'hzi')) if testHJ: def test_Jvec_jxr_Jform(self): From f1527f994bd6cc547edb62ea8e037e327f73c3fb Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sat, 20 Feb 2016 14:30:55 -0800 Subject: [PATCH 11/17] e,b,h,j from j formulation --- SimPEG/EM/FDEM/FieldsFDEM.py | 232 +++++++++++------- .../em/fdem/forward/test_FDEM_forwardEJHB.py | 60 ++++- .../inverse/adjoint/test_FDEM_adjointHJ.py | 76 ++++-- .../fdem/inverse/derivs/test_FDEM_derivs.py | 109 ++++++-- 4 files changed, 338 insertions(+), 139 deletions(-) diff --git a/SimPEG/EM/FDEM/FieldsFDEM.py b/SimPEG/EM/FDEM/FieldsFDEM.py index b8971611..dbf4c77c 100644 --- a/SimPEG/EM/FDEM/FieldsFDEM.py +++ b/SimPEG/EM/FDEM/FieldsFDEM.py @@ -774,8 +774,12 @@ class Fields_j(Fields): 'h' : ['jSolution','E','_h'], 'hPrimary' : ['jSolution','E','_hPrimary'], 'hSecondary' : ['jSolution','E','_hSecondary'], - 'e' : ['jSolution','C','_e'], - 'b' : ['jSolution','C','_b'], + 'e' : ['jSolution','CCV','_e'], + 'ePrimary' : ['jSolution','CCV','_ePrimary'], + 'eSecondary' : ['jSolution','CCV','_eSecondary'], + 'b' : ['jSolution','CCV','_b'], + 'bPrimary' : ['jSolution','CCV','_bPrimary'], + 'bSecondary' : ['jSolution','CCV','_bSecondary'], } def __init__(self,mesh,survey,**kwargs): @@ -784,10 +788,10 @@ class Fields_j(Fields): def startup(self): self.prob = self.survey.prob self._edgeCurl = self.survey.prob.mesh.edgeCurl + self._MeMu = self.survey.prob.MeMu self._MeMuI = self.survey.prob.MeMuI 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._mu = self.survey.prob.curModel.mui self._aveF2CCV = self.survey.prob.mesh.aveF2CCV @@ -907,60 +911,9 @@ class Fields_j(Fields): return h - def _hSecondaryDeriv_u(self, src, du_dm_v, adjoint=False): - """ - Derivative of the secondary magnetic field with respect to the thing we solved for - - :param SimPEG.EM.FDEM.Src src: source - :param numpy.ndarray du_dm_v: vector to take product with - :param bool adjoint: adjoint? - :rtype: numpy.ndarray - :return: product of the derivative of the secondary magnetic field with respect to the field we solved for with a vector - """ - - if not adjoint: - return -1./(1j*omega(src.freq)) * self._MeMuI * (self._edgeCurl.T * (self._MfRho * du_dm_v) ) - elif adjoint: - return -1./(1j*omega(src.freq)) * self._MfRho.T * (self._edgeCurl * ( self._MeMuI.T * du_dm_v)) - - - def _hSecondaryDeriv_m(self, src, v, adjoint=False): - """ - Derivative of the secondary magnetic field with respect to the inversion model - - :param SimPEG.EM.FDEM.Src src: source - :param numpy.ndarray v: vector to take product with - :param bool adjoint: adjoint? - :rtype: numpy.ndarray - :return: product of the derivative of the secondary magnetic field with respect to the model with a vector - """ - - jSolution = self[[src],'jSolution'] - MeMuI = self._MeMuI - C = self._edgeCurl - MfRho = self._MfRho - MfRhoDeriv = self._MfRhoDeriv - Me = self._Me - - if not adjoint: - hDeriv_m = -1./(1j*omega(src.freq)) * MeMuI * (C.T * (MfRhoDeriv(jSolution)*v ) ) - elif adjoint: - hDeriv_m = -1./(1j*omega(src.freq)) * MfRhoDeriv(jSolution).T * ( C * (MeMuI.T * v ) ) - - S_mDeriv,_ = src.evalDeriv(self.prob, adjoint = adjoint) - - if not adjoint: - S_mDeriv = S_mDeriv(v) - hDeriv_m = hDeriv_m + 1./(1j*omega(src.freq)) * MeMuI * (Me * S_mDeriv) - elif adjoint: - S_mDeriv = S_mDeriv(Me.T * (MeMuI.T * v)) - hDeriv_m = hDeriv_m + 1./(1j*omega(src.freq)) * S_mDeriv - return hDeriv_m - - def _hDeriv_u(self, src, du_dm_v, adjoint=False): """ - Partial derivative of the total magnetic field with respect to the thing we solved for + Derivative of the magnetic field with respect to the thing we solved for :param SimPEG.EM.FDEM.Src src: source :param numpy.ndarray du_dm_v: vector to take product with @@ -969,50 +922,165 @@ class Fields_j(Fields): :return: product of the derivative of the magnetic field with respect to the field we solved for with a vector """ - return self._hSecondaryDeriv_u(src, du_dm_v, adjoint) + if not adjoint: + return -1./(1j*omega(src.freq)) * self._MeMuI * (self._edgeCurl.T * (self._MfRho * du_dm_v) ) + return -1./(1j*omega(src.freq)) * self._MfRho.T * (self._edgeCurl * ( self._MeMuI.T * du_dm_v)) + def _hDeriv_m(self, src, v, adjoint=False): """ - Partial derivative of the total magnetic field density with respect to the inversion model. + Derivative of the magnetic field with respect to the inversion model :param SimPEG.EM.FDEM.Src src: source :param numpy.ndarray v: vector to take product with :param bool adjoint: adjoint? :rtype: numpy.ndarray - :return: product of the magnetic field derivative with respect to the inversion model with a vector + :return: product of the derivative of the magnetic field with respect to the model with a vector """ - # assuming the primary doesn't depend on the model - return self._hSecondaryDeriv_m(src, u, v, adjoint) + jSolution = Utils.mkvc(self[[src],'jSolution']) + MeMuI = self._MeMuI + C = self._edgeCurl + MfRho = self._MfRho + MfRhoDeriv = self._MfRhoDeriv + S_mDeriv,_ = src.evalDeriv(self.prob, adjoint = 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 + if not adjoint: + hDeriv_m = -1./(1j*omega(src.freq)) * MeMuI * (C.T * (MfRhoDeriv(jSolution)*v ) ) + S_mDeriv = S_mDeriv(v) + hDeriv_m = hDeriv_m + 1./(1j*omega(src.freq)) * MeMuI * ( S_mDeriv) + + elif adjoint: + hDeriv_m = -1./(1j*omega(src.freq)) * MfRhoDeriv(jSolution).T * ( C * (MeMuI.T * v ) ) + + S_mDeriv = S_mDeriv(MeMuI.T * v) + hDeriv_m = hDeriv_m + 1./(1j*omega(src.freq)) * S_mDeriv - Rho = self.prob.MfRho + return hDeriv_m + + + def _ePrimary(self, jSolution, srcList): + """ + Primary electric field + + :param numpy.ndarray hSolution: field we solved for + :param list srcList: list of sources + :rtype: numpy.ndarray + :return: primary electric field as defined by the sources + """ + n = int(self._aveF2CCV.shape[0] / self._nC) # number of components + VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) + return VI * (self._aveF2CCV * (self._MfRho * self._jPrimary(jSolution, srcList))) + + def _eSecondary(self, jSolution, srcList): + """ + Secondary electric field from jSolution + + :param numpy.ndarray hSolution: field we solved for + :param list srcList: list of sources + :rtype: numpy.ndarray + :return: secondary electric field + """ + n = int(self._aveF2CCV.shape[0] / self._nC) # number of components + VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) + return VI * (self._aveF2CCV * (self._MfRho * self._jSecondary(jSolution, srcList))) + + def _eDeriv_u(self, src, du_dm_v, adjoint=False): + """ + Derivative of the electric field with respect to the thing we solved for + + :param SimPEG.EM.FDEM.Src src: source + :param numpy.ndarray du_dm_v: vector to take product with + :param bool adjoint: adjoint? + :rtype: numpy.ndarray + :return: product of the derivative of the electric field with respect to the field we solved for with a vector + """ + n = int(self._aveF2CCV.shape[0] / self._nC) # number of components + VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) + if adjoint: + return self._MfRho.T * ( self._aveF2CCV.T * ( VI.T * du_dm_v ) ) + return VI * (self._aveF2CCV * (self._MfRho * du_dm_v)) + + def _eDeriv_m(self, src, v, adjoint=False): + """ + Derivative of the electric field with respect to the inversion model + + :param SimPEG.EM.FDEM.Src src: source + :param numpy.ndarray v: vector to take product with + :param bool adjoint: adjoint? + :rtype: numpy.ndarray + :return: product of the derivative of the electric field with respect to the model with a vector + """ + jSolution = Utils.mkvc(self[src,'jSolution']) + n = int(self._aveF2CCV.shape[0] / self._nC) # number of components + VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) + if adjoint: + return self._MfRhoDeriv(jSolution).T * ( self._aveF2CCV.T * ( VI.T * v ) ) + return VI * (self._aveF2CCV * (self._MfRhoDeriv(jSolution) * v)) + + def _bPrimary(self, jSolution, srcList): + """ + Primary magnetic flux density + + :param numpy.ndarray hSolution: field we solved for + :param list srcList: list of sources + :rtype: numpy.ndarray + :return: primary magnetic flux density + """ + hPrimary = self._hPrimary(jSolution, srcList) + n = int(self._aveE2CCV.shape[0] / self._nC) # number of components VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) + return VI * (self._aveE2CCV * (self._MeMu * hPrimary)) - j = self._j(jSolution, srcList) + def _bSecondary(self, jSolution, srcList): + """ + Secondary magnetic flux density from jSolution - return VI * (aveF2CCV * (Rho * j)) - - def _eDeriv_u(self, src, u, v, adjoint=False): - raise NotImplementedError - - def _eDeriv_m(self, src, u, v, adjoint=False): - raise NotImplementedError - - def _b(self, jSolution, srcList): - h = self._h(jSolution, srcList) - Mu = self.prob.MeMu - aveE2CCV = self._aveE2CCV - n = int(aveE2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy - # Mu = sdiag(sp.kron(np.ones(n), mu)) + :param numpy.ndarray hSolution: field we solved for + :param list srcList: list of sources + :rtype: numpy.ndarray + :return: secondary magnetic flux density + """ + n = int(self._aveE2CCV.shape[0] / self._nC) # number of components VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) - return VI * (aveE2CCV * (Mu * h)) + return VI * (self._aveE2CCV * (self._edgeCurl.T * (self._MfRho * jSolution))) + + def _bDeriv_u(self, src, du_dm_v, adjoint=False): + """ + Derivative of the magnetic flux density with respect to the thing we solved for + + :param SimPEG.EM.FDEM.Src src: source + :param numpy.ndarray du_dm_v: vector to take product with + :param bool adjoint: adjoint? + :rtype: numpy.ndarray + :return: product of the derivative of the magnetic flux density with respect to the field we solved for with a vector + """ + n = int(self._aveF2CCV.shape[0] / self._nC) # number of components + VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) + + if adjoint: + return self._MfRho.T * ( self._edgeCurl * ( self._aveE2CCV.T * (VI.T * du_dm_v) ) ) + return VI * (self._aveE2CCV * (self._edgeCurl.T * (self._MfRho * du_dm_v))) + + def _bDeriv_m(self, src, v, adjoint=False): + """ + Derivative of the magnetic flux density with respect to the inversion model + + :param SimPEG.EM.FDEM.Src src: source + :param numpy.ndarray v: vector to take product with + :param bool adjoint: adjoint? + :rtype: numpy.ndarray + :return: product of the derivative of the magnetic flux density with respect to the model with a vector + """ + jSolution = self[src,'jSolution'] + n = int(self._aveE2CCV.shape[0] / self._nC) # number of components + VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) + + if adjoint: + return self._MfRhoDeriv(jSolution).T * ( self._edgeCurl * (self._aveE2CCV.T * ( VI.T * v)) ) + return VI * (self._aveE2CCV * (self._edgeCurl.T * (self._MfRhoDeriv(jSolution) * v ))) class Fields_h(Fields): @@ -1208,7 +1276,7 @@ class Fields_h(Fields): def _e(self, hSolution, srcList): rho = self._rho aveF2CCV = self._aveF2CCV - n = int(aveF2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy + n = int(self._aveF2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy Rho = self.prob.MfRho VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) diff --git a/tests/em/fdem/forward/test_FDEM_forwardEJHB.py b/tests/em/fdem/forward/test_FDEM_forwardEJHB.py index c01e70db..d4218b86 100644 --- a/tests/em/fdem/forward/test_FDEM_forwardEJHB.py +++ b/tests/em/fdem/forward/test_FDEM_forwardEJHB.py @@ -5,12 +5,9 @@ import sys from scipy.constants import mu_0 from SimPEG.EM.Utils.testingUtils import getFDEMProblem, crossCheckTest -testEB = False -testHJ = False -testEJ = False -testBH = False +testEJ = True +testBH = True -TOLEBHJ = 1e-5 TOLEJHB = 1 # averaging and more sensitive to boundary condition violations (ie. the impact of violating the boundary conditions in each case is different.) #TODO: choose better testing parameters to lower this @@ -71,5 +68,58 @@ class FDEM_CrossCheck(unittest.TestCase): def test_EJ_CrossCheck_hzi_Jform(self): self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'hzi', TOL=TOLEJHB)) + if testBH: + def test_HB_CrossCheck_jxr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'jxr', TOL=TOLEJHB)) + def test_HB_CrossCheck_jyr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'jyr', TOL=TOLEJHB)) + def test_HB_CrossCheck_jzr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'jzr', TOL=TOLEJHB)) + def test_HB_CrossCheck_jxi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'jxi', TOL=TOLEJHB)) + def test_HB_CrossCheck_jyi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'jyi', TOL=TOLEJHB)) + def test_HB_CrossCheck_jzi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'jzi', TOL=TOLEJHB)) + + def test_HB_CrossCheck_exr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'exr', TOL=TOLEJHB)) + def test_HB_CrossCheck_eyr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'eyr', TOL=TOLEJHB)) + def test_HB_CrossCheck_ezr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'ezr', TOL=TOLEJHB)) + def test_HB_CrossCheck_exi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'exi', TOL=TOLEJHB)) + def test_HB_CrossCheck_eyi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'eyi', TOL=TOLEJHB)) + def test_HB_CrossCheck_ezi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'ezi', TOL=TOLEJHB)) + + def test_HB_CrossCheck_bxr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'bxr', TOL=TOLEJHB)) + def test_HB_CrossCheck_byr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'byr', TOL=TOLEJHB)) + def test_HB_CrossCheck_bzr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'bzr', TOL=TOLEJHB)) + def test_HB_CrossCheck_bxi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'bxi', TOL=TOLEJHB)) + def test_HB_CrossCheck_byi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'byi', TOL=TOLEJHB)) + def test_HB_CrossCheck_bzi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'bzi', TOL=TOLEJHB)) + + def test_HB_CrossCheck_hxr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'hxr', TOL=TOLEJHB)) + def test_HB_CrossCheck_hyr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'hyr', TOL=TOLEJHB)) + def test_HB_CrossCheck_hzr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'hzr', TOL=TOLEJHB)) + def test_HB_CrossCheck_hxi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'hxi', TOL=TOLEJHB)) + def test_HB_CrossCheck_hyi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'hyi', TOL=TOLEJHB)) + def test_HB_CrossCheck_hzi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'hzi', TOL=TOLEJHB)) + if __name__ == '__main__': unittest.main() \ No newline at end of file diff --git a/tests/em/fdem/inverse/adjoint/test_FDEM_adjointHJ.py b/tests/em/fdem/inverse/adjoint/test_FDEM_adjointHJ.py index 68224e34..242b33f2 100644 --- a/tests/em/fdem/inverse/adjoint/test_FDEM_adjointHJ.py +++ b/tests/em/fdem/inverse/adjoint/test_FDEM_adjointHJ.py @@ -71,33 +71,59 @@ class FDEM_AdjointTests(unittest.TestCase): def test_Jtvec_adjointTest_hyi_Jform(self): self.assertTrue(adjointTest('j', 'hyi')) def test_Jtvec_adjointTest_hzi_Jform(self): - self.assertTrue(adjointTest('j', 'hzi')) + self.assertTrue(adjointTest('j', 'hzi')) - def test_Jtvec_adjointTest_hxr_Hform(self): - self.assertTrue(adjointTest('h', 'hxr')) - def test_Jtvec_adjointTest_hyr_Hform(self): - self.assertTrue(adjointTest('h', 'hyr')) - def test_Jtvec_adjointTest_hzr_Hform(self): - self.assertTrue(adjointTest('h', 'hzr')) - def test_Jtvec_adjointTest_hxi_Hform(self): - self.assertTrue(adjointTest('h', 'hxi')) - def test_Jtvec_adjointTest_hyi_Hform(self): - self.assertTrue(adjointTest('h', 'hyi')) - def test_Jtvec_adjointTest_hzi_Hform(self): - self.assertTrue(adjointTest('h', 'hzi')) + def test_Jtvec_adjointTest_exr_Jform(self): + self.assertTrue(adjointTest('j', 'exr')) + def test_Jtvec_adjointTest_eyr_Jform(self): + self.assertTrue(adjointTest('j', 'eyr')) + def test_Jtvec_adjointTest_ezr_Jform(self): + self.assertTrue(adjointTest('j', 'ezr')) + def test_Jtvec_adjointTest_exi_Jform(self): + self.assertTrue(adjointTest('j', 'exi')) + def test_Jtvec_adjointTest_eyi_Jform(self): + self.assertTrue(adjointTest('j', 'eyi')) + def test_Jtvec_adjointTest_ezi_Jform(self): + self.assertTrue(adjointTest('j', 'ezi')) - def test_Jtvec_adjointTest_hxr_Hform(self): - self.assertTrue(adjointTest('h', 'jxr')) - def test_Jtvec_adjointTest_hyr_Hform(self): - self.assertTrue(adjointTest('h', 'jyr')) - def test_Jtvec_adjointTest_hzr_Hform(self): - self.assertTrue(adjointTest('h', 'jzr')) - def test_Jtvec_adjointTest_hxi_Hform(self): - self.assertTrue(adjointTest('h', 'jxi')) - def test_Jtvec_adjointTest_hyi_Hform(self): - self.assertTrue(adjointTest('h', 'jyi')) - def test_Jtvec_adjointTest_hzi_Hform(self): - self.assertTrue(adjointTest('h', 'jzi')) + def test_Jtvec_adjointTest_bxr_Jform(self): + self.assertTrue(adjointTest('j', 'bxr')) + def test_Jtvec_adjointTest_byr_Jform(self): + self.assertTrue(adjointTest('j', 'byr')) + def test_Jtvec_adjointTest_bzr_Jform(self): + self.assertTrue(adjointTest('j', 'bzr')) + def test_Jtvec_adjointTest_bxi_Jform(self): + self.assertTrue(adjointTest('j', 'bxi')) + def test_Jtvec_adjointTest_byi_Jform(self): + self.assertTrue(adjointTest('j', 'byi')) + def test_Jtvec_adjointTest_bzi_Jform(self): + self.assertTrue(adjointTest('j', 'bzi')) + + # def test_Jtvec_adjointTest_hxr_Hform(self): + # self.assertTrue(adjointTest('h', 'hxr')) + # def test_Jtvec_adjointTest_hyr_Hform(self): + # self.assertTrue(adjointTest('h', 'hyr')) + # def test_Jtvec_adjointTest_hzr_Hform(self): + # self.assertTrue(adjointTest('h', 'hzr')) + # def test_Jtvec_adjointTest_hxi_Hform(self): + # self.assertTrue(adjointTest('h', 'hxi')) + # def test_Jtvec_adjointTest_hyi_Hform(self): + # self.assertTrue(adjointTest('h', 'hyi')) + # def test_Jtvec_adjointTest_hzi_Hform(self): + # self.assertTrue(adjointTest('h', 'hzi')) + + # def test_Jtvec_adjointTest_hxr_Hform(self): + # self.assertTrue(adjointTest('h', 'jxr')) + # def test_Jtvec_adjointTest_hyr_Hform(self): + # self.assertTrue(adjointTest('h', 'jyr')) + # def test_Jtvec_adjointTest_hzr_Hform(self): + # self.assertTrue(adjointTest('h', 'jzr')) + # def test_Jtvec_adjointTest_hxi_Hform(self): + # self.assertTrue(adjointTest('h', 'jxi')) + # def test_Jtvec_adjointTest_hyi_Hform(self): + # self.assertTrue(adjointTest('h', 'jyi')) + # def test_Jtvec_adjointTest_hzi_Hform(self): + # self.assertTrue(adjointTest('h', 'jzi')) if __name__ == '__main__': diff --git a/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py b/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py index 19f3e9e6..8c233dfa 100644 --- a/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py +++ b/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py @@ -5,10 +5,11 @@ import sys from scipy.constants import mu_0 from SimPEG.EM.Utils.testingUtils import getFDEMProblem -testDerivs = True + testE = True testB = True -testHJ = False +testH = True +testJ = True verbose = False @@ -98,31 +99,31 @@ class FDEM_DerivTests(unittest.TestCase): self.assertTrue(derivTest('e', 'hzi')) if testB: - # def test_Jvec_exr_Bform(self): - # self.assertTrue(derivTest('b', 'exr')) - # def test_Jvec_eyr_Bform(self): - # self.assertTrue(derivTest('b', 'eyr')) - # def test_Jvec_ezr_Bform(self): - # self.assertTrue(derivTest('b', 'ezr')) - # def test_Jvec_exi_Bform(self): - # self.assertTrue(derivTest('b', 'exi')) - # def test_Jvec_eyi_Bform(self): - # self.assertTrue(derivTest('b', 'eyi')) - # def test_Jvec_ezi_Bform(self): - # self.assertTrue(derivTest('b', 'ezi')) + def test_Jvec_exr_Bform(self): + self.assertTrue(derivTest('b', 'exr')) + def test_Jvec_eyr_Bform(self): + self.assertTrue(derivTest('b', 'eyr')) + def test_Jvec_ezr_Bform(self): + self.assertTrue(derivTest('b', 'ezr')) + def test_Jvec_exi_Bform(self): + self.assertTrue(derivTest('b', 'exi')) + def test_Jvec_eyi_Bform(self): + self.assertTrue(derivTest('b', 'eyi')) + def test_Jvec_ezi_Bform(self): + self.assertTrue(derivTest('b', 'ezi')) - # def test_Jvec_bxr_Bform(self): - # self.assertTrue(derivTest('b', 'bxr')) - # def test_Jvec_byr_Bform(self): - # self.assertTrue(derivTest('b', 'byr')) - # def test_Jvec_bzr_Bform(self): - # self.assertTrue(derivTest('b', 'bzr')) - # def test_Jvec_bxi_Bform(self): - # self.assertTrue(derivTest('b', 'bxi')) - # def test_Jvec_byi_Bform(self): - # self.assertTrue(derivTest('b', 'byi')) - # def test_Jvec_bzi_Bform(self): - # self.assertTrue(derivTest('b', 'bzi')) + def test_Jvec_bxr_Bform(self): + self.assertTrue(derivTest('b', 'bxr')) + def test_Jvec_byr_Bform(self): + self.assertTrue(derivTest('b', 'byr')) + def test_Jvec_bzr_Bform(self): + self.assertTrue(derivTest('b', 'bzr')) + def test_Jvec_bxi_Bform(self): + self.assertTrue(derivTest('b', 'bxi')) + def test_Jvec_byi_Bform(self): + self.assertTrue(derivTest('b', 'byi')) + def test_Jvec_bzi_Bform(self): + self.assertTrue(derivTest('b', 'bzi')) def test_Jvec_jxr_Bform(self): self.assertTrue(derivTest('b', 'jxr')) @@ -150,7 +151,7 @@ class FDEM_DerivTests(unittest.TestCase): def test_Jvec_hzi_Bform(self): self.assertTrue(derivTest('b', 'hzi')) - if testHJ: + if testJ: def test_Jvec_jxr_Jform(self): self.assertTrue(derivTest('j', 'jxr')) def test_Jvec_jyr_Jform(self): @@ -177,6 +178,34 @@ class FDEM_DerivTests(unittest.TestCase): def test_Jvec_hzi_Jform(self): self.assertTrue(derivTest('j', 'hzi')) + def test_Jvec_exr_Jform(self): + self.assertTrue(derivTest('j', 'exr')) + def test_Jvec_eyr_Jform(self): + self.assertTrue(derivTest('j', 'eyr')) + def test_Jvec_ezr_Jform(self): + self.assertTrue(derivTest('j', 'ezr')) + def test_Jvec_exi_Jform(self): + self.assertTrue(derivTest('j', 'exi')) + def test_Jvec_eyi_Jform(self): + self.assertTrue(derivTest('j', 'eyi')) + def test_Jvec_ezi_Jform(self): + self.assertTrue(derivTest('j', 'ezi')) + + def test_Jvec_bxr_Jform(self): + self.assertTrue(derivTest('j', 'bxr')) + def test_Jvec_byr_Jform(self): + self.assertTrue(derivTest('j', 'byr')) + def test_Jvec_bzr_Jform(self): + self.assertTrue(derivTest('j', 'bzr')) + def test_Jvec_bxi_Jform(self): + self.assertTrue(derivTest('j', 'bxi')) + def test_Jvec_byi_Jform(self): + self.assertTrue(derivTest('j', 'byi')) + def test_Jvec_bzi_Jform(self): + self.assertTrue(derivTest('j', 'bzi')) + + + if testH: def test_Jvec_hxr_Hform(self): self.assertTrue(derivTest('h', 'hxr')) def test_Jvec_hyr_Hform(self): @@ -203,6 +232,32 @@ class FDEM_DerivTests(unittest.TestCase): def test_Jvec_hzi_Hform(self): self.assertTrue(derivTest('h', 'jzi')) + def test_Jvec_exr_Hform(self): + self.assertTrue(derivTest('h', 'exr')) + def test_Jvec_eyr_Hform(self): + self.assertTrue(derivTest('h', 'eyr')) + def test_Jvec_ezr_Hform(self): + self.assertTrue(derivTest('h', 'ezr')) + def test_Jvec_exi_Hform(self): + self.assertTrue(derivTest('h', 'exi')) + def test_Jvec_eyi_Hform(self): + self.assertTrue(derivTest('h', 'eyi')) + def test_Jvec_ezi_Hform(self): + self.assertTrue(derivTest('h', 'ezi')) + + def test_Jvec_bxr_Hform(self): + self.assertTrue(derivTest('h', 'bxr')) + def test_Jvec_byr_Hform(self): + self.assertTrue(derivTest('h', 'byr')) + def test_Jvec_bzr_Hform(self): + self.assertTrue(derivTest('h', 'bzr')) + def test_Jvec_bxi_Hform(self): + self.assertTrue(derivTest('h', 'bxi')) + def test_Jvec_byi_Hform(self): + self.assertTrue(derivTest('h', 'byi')) + def test_Jvec_bzi_Hform(self): + self.assertTrue(derivTest('h', 'bzi')) + if __name__ == '__main__': unittest.main() From 3e673d34f180206eb393cc12ad078a67e77fb4db Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sat, 20 Feb 2016 17:05:43 -0800 Subject: [PATCH 12/17] - each of e,b,h,j from every formulation. Currently, b from j is first order - removed CCV primary and secondary (dangerous the way it was previously done) - NOTE: Source derive may not be properly taken care of yet --- SimPEG/EM/FDEM/FieldsFDEM.py | 336 +++++++----------- tests/em/fdem/forward/test_FDEM_forward.py | 46 +-- tests/em/fdem/forward/test_FDEM_forwardHB.py | 4 +- .../inverse/adjoint/test_FDEM_adjointHJ.py | 81 +++-- 4 files changed, 217 insertions(+), 250 deletions(-) diff --git a/SimPEG/EM/FDEM/FieldsFDEM.py b/SimPEG/EM/FDEM/FieldsFDEM.py index dbf4c77c..68843d06 100644 --- a/SimPEG/EM/FDEM/FieldsFDEM.py +++ b/SimPEG/EM/FDEM/FieldsFDEM.py @@ -178,11 +178,7 @@ class Fields_e(Fields): 'bPrimary' : ['eSolution','F','_bPrimary'], 'bSecondary' : ['eSolution','F','_bSecondary'], 'j' : ['eSolution','CCV','_j'], - 'jPrimary' : ['eSolution','CCV','_jPrimary'], - 'jSecondary' : ['eSolution','CCV','_jSecondary'], 'h' : ['eSolution','CCV','_h'], - 'hPrimary' : ['eSolution','CCV','_hPrimary'], - 'hSecondary' : ['eSolution','CCV','_hSecondary'], } def __init__(self, mesh, survey, **kwargs): @@ -219,7 +215,7 @@ class Fields_e(Fields): :return: primary electric field as defined by the sources """ - ePrimary = np.zeros([self.prob.mesh.nE,len(srcList)]) + ePrimary = np.zeros([self.prob.mesh.nE,len(srcList)], dtype = complex) for i, src in enumerate(srcList): ep = src.ePrimary(self.prob) ePrimary[:,i] = ePrimary[:,i] + ep @@ -274,7 +270,7 @@ class Fields_e(Fields): :return: primary magnetic flux density as defined by the sources """ - bPrimary = np.zeros([self._edgeCurl.shape[0],eSolution.shape[1]],dtype = complex) + bPrimary = np.zeros([self._edgeCurl.shape[0],eSolution.shape[1]], dtype = complex) for i, src in enumerate(srcList): bp = src.bPrimary(self.prob) bPrimary[:,i] = bPrimary[:,i] + bp @@ -329,34 +325,19 @@ class Fields_e(Fields): S_mDeriv, _ = src.evalDeriv(self.prob, v, adjoint) return 1./(1j * omega(src.freq)) * S_mDeriv - - def _jPrimary(self, eSolution, srcList): + def _j(self, eSolution, srcList): """ - Primary current density + Current density from eSolution :param numpy.ndarray eSolution: field we solved for :param list srcList: list of sources :rtype: numpy.ndarray - :return: primary current density + :return: current density """ aveE2CCV = self._aveE2CCV n = int(aveE2CCV.shape[0] / self._nC) # number of components (instead of checking if cyl or not) VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) - return VI * (aveE2CCV * (self._MeSigma * self._ePrimary(eSolution, srcList)) ) - - def _jSecondary(self, eSolution, srcList): - """ - Secondary current density from eSolution - - :param numpy.ndarray eSolution: field we solved for - :param list srcList: list of sources - :rtype: numpy.ndarray - :return: secondary current density - """ - aveE2CCV = self._aveE2CCV - n = int(aveE2CCV.shape[0] / self._nC) # number of components (instead of checking if cyl or not) - VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) - return VI * (aveE2CCV * (self._MeSigma * eSolution) ) + return VI * (aveE2CCV * (self._MeSigma * self._e(eSolution,srcList) ) ) def _jDeriv_u(self, src, du_dm_v, adjoint = False): """ @@ -396,33 +377,19 @@ class Fields_e(Fields): - def _hPrimary(self, eSolution, srcList): + def _h(self, eSolution, srcList): """ - Primary magnetic field + Magnetic field from eSolution :param numpy.ndarray eSolution: field we solved for :param list srcList: list of sources :rtype: numpy.ndarray - :return: primary magnetic field + :return: magnetic field """ n = int(self._aveF2CCV.shape[0] / self._nC) # Number of Components VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) - return VI * (self._aveF2CCV * (self._MfMui * self._bPrimary(eSolution, srcList))) - - def _hSecondary(self, eSolution, srcList): - """ - Secondary magnetic field from eSolution - - :param numpy.ndarray eSolution: field we solved for - :param list srcList: list of sources - :rtype: numpy.ndarray - :return: secondary magnetic field - """ - n = int(self._aveF2CCV.shape[0] / self._nC) # Number of Components - VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) - - return VI * (self._aveF2CCV * (self._MfMui * self._bSecondary(eSolution, srcList))) + return VI * (self._aveF2CCV * (self._MfMui * self._b(eSolution, srcList))) def _hDeriv_u(self, src, du_dm_v, adjoint = False): """ @@ -477,11 +444,7 @@ class Fields_b(Fields): 'ePrimary' : ['bSolution','E','_ePrimary'], 'eSecondary' : ['bSolution','E','_eSecondary'], 'j' : ['bSolution','CCV','_j'], - 'jPrimary' : ['bSolution','CCV','_jPrimary'], - 'jSecondary' : ['bSolution','CCV','_jSecondary'], 'h' : ['bSolution','CCV','_h'], - 'hPrimary' : ['bSolution','CCV','_hPrimary'], - 'hSecondary' : ['bSolution','CCV','_hSecondary'], } def __init__(self,mesh,survey,**kwargs): @@ -490,8 +453,10 @@ class Fields_b(Fields): def startup(self): self.prob = self.survey.prob self._edgeCurl = self.survey.prob.mesh.edgeCurl + self._MeSigma = self.survey.prob.MeSigma self._MeSigmaI = self.survey.prob.MeSigmaI self._MfMui = self.survey.prob.MfMui + self._MeSigmaDeriv = self.survey.prob.MeSigmaDeriv self._MeSigmaIDeriv = self.survey.prob.MeSigmaIDeriv self._Me = self.survey.prob.Me self._aveF2CCV = self.survey.prob.mesh.aveF2CCV @@ -522,7 +487,7 @@ class Fields_b(Fields): :return: primary electric field as defined by the sources """ - bPrimary = np.zeros([self.prob.mesh.nF,len(srcList)]) + bPrimary = np.zeros([self.prob.mesh.nF,len(srcList)], dtype = complex) for i, src in enumerate(srcList): bp = src.bPrimary(self.prob) bPrimary[:,i] = bPrimary[:,i] + bp @@ -578,7 +543,7 @@ class Fields_b(Fields): :return: primary electric field as defined by the sources """ - ePrimary = np.zeros([self._edgeCurl.shape[1],bSolution.shape[1]],dtype = complex) + ePrimary = np.zeros([self._edgeCurl.shape[1],bSolution.shape[1]], dtype = complex) for i,src in enumerate(srcList): ep = src.ePrimary(self.prob) ePrimary[:,i] = ePrimary[:,i] + ep @@ -594,11 +559,12 @@ class Fields_b(Fields): :return: secondary electric field """ - e = self._MeSigmaI * ( self._edgeCurl.T * ( self._MfMui * bSolution)) + v = ( self._edgeCurl.T * ( self._MfMui * bSolution)) for i,src in enumerate(srcList): _,S_e = src.eval(self.prob) - e[:,i] = e[:,i] + -self._MeSigmaI * S_e - return e + v[:,i] = v[:,i] + - S_e + + return self._MeSigmaI * v def _eDeriv_u(self, src, du_dm_v, adjoint=False): """ @@ -636,24 +602,9 @@ class Fields_b(Fields): if adjoint: return self._MeSigmaIDeriv(w).T * v - self._MeSigmaI.T * S_eDeriv - return self._MeSigmaIDeriv(w) * v - self._MeSigmaI *S_eDeriv + return self._MeSigmaIDeriv(w) * v - self._MeSigmaI * S_eDeriv - - def _jPrimary(self, bSolution, srcList): - """ - Primary current density - - :param numpy.ndarray bSolution: field we solved for - :param list srcList: list of sources - :rtype: numpy.ndarray - :return: primary current density - """ - n = int(self._aveE2CCV.shape[0] / self._nC) # number of components - VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) - - return VI * (self._aveE2CCV * (self.prob.MeSigma * self._ePrimary(bSolution, srcList))) - - def _jSecondary(self, bSolution, srcList): + def _j(self, bSolution, srcList): """ Secondary current density from bSolution @@ -666,7 +617,8 @@ class Fields_b(Fields): n = int(self._aveE2CCV.shape[0] / self._nC) # number of components VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) - return VI * (self._aveE2CCV * ( self._edgeCurl.T * ( self._MfMui * bSolution) ) ) + return VI * (self._aveE2CCV * ( self._MeSigma * self._e(bSolution,srcList ) ) ) + def _jDeriv_u(self, src, du_dm_v, adjoint=False): """ @@ -679,13 +631,10 @@ class Fields_b(Fields): :rtype: numpy.ndarray :return: product of the derivative of the current density with respect to the field we solved for with a vector """ - n = int(self._aveE2CCV.shape[0] / self._nC) # number of components VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) + return VI * (self._aveE2CCV * (self._edgeCurl.T * ( self._MfMui * du_dm_v ) ) ) - if adjoint: - return self._MfMui.T * ( self._edgeCurl * ( self._aveE2CCV.T * ( VI.T * du_dm_v ) ) ) - return VI * ( self._aveE2CCV * ( self._edgeCurl.T * ( self._MfMui * du_dm_v ) ) ) def _jDeriv_m(self, src, v, adjoint=False): """ @@ -699,31 +648,18 @@ class Fields_b(Fields): """ return Zero() - def _hPrimary(self, bSolution, srcList): + def _h(self, bSolution, srcList): """ - Primary magnetic field + Magnetic field from bSolution :param numpy.ndarray bSolution: field we solved for :param list srcList: list of sources :rtype: numpy.ndarray - :return: primary magnetic field + :return: magnetic field """ n = int(self._aveF2CCV.shape[0] / self._nC) #number of components VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) - return VI * (self._aveF2CCV * (self._MfMui * self._bPrimary(bSolution, srcList))) - - def _hSecondary(self, bSolution, srcList): - """ - Secondary magnetic field from bSolution - - :param numpy.ndarray bSolution: field we solved for - :param list srcList: list of sources - :rtype: numpy.ndarray - :return: secondary magnetic field - """ - n = int(self._aveF2CCV.shape[0] / self._nC) #number of components - VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) - return VI * (self._aveF2CCV * (self._MfMui * self._bSecondary(bSolution, srcList))) + return VI * (self._aveF2CCV * (self._MfMui * self._b(bSolution, srcList))) def _hDeriv_u(self, src, du_dm_v, adjoint=False): """ @@ -756,8 +692,6 @@ class Fields_b(Fields): return Zero() - - class Fields_j(Fields): """ Fields object for Problem_j. @@ -775,11 +709,7 @@ class Fields_j(Fields): 'hPrimary' : ['jSolution','E','_hPrimary'], 'hSecondary' : ['jSolution','E','_hSecondary'], 'e' : ['jSolution','CCV','_e'], - 'ePrimary' : ['jSolution','CCV','_ePrimary'], - 'eSecondary' : ['jSolution','CCV','_eSecondary'], 'b' : ['jSolution','CCV','_b'], - 'bPrimary' : ['jSolution','CCV','_bPrimary'], - 'bSecondary' : ['jSolution','CCV','_bSecondary'], } def __init__(self,mesh,survey,**kwargs): @@ -818,7 +748,7 @@ class Fields_j(Fields): :return: primary current density as defined by the sources """ - jPrimary = np.zeros_like(jSolution,dtype = complex) + jPrimary = np.zeros_like(jSolution, dtype = complex) for i, src in enumerate(srcList): jp = src.jPrimary(self.prob) jPrimary[:,i] = jPrimary[:,i] + jp @@ -903,12 +833,12 @@ class Fields_j(Fields): :return: secondary magnetic field """ - h = self._MeMuI * (self._edgeCurl.T * (self._MfRho * jSolution) ) + v = (self._edgeCurl.T * (self._MfRho * jSolution) ) for i, src in enumerate(srcList): h[:,i] *= -1./(1j*omega(src.freq)) S_m,_ = src.eval(self.prob) - h[:,i] = h[:,i]+ 1./(1j*omega(src.freq)) * self._MeMuI * (S_m) - return h + h[:,i] = h[:,i]+ 1./(1j*omega(src.freq)) * (S_m) + return self._MeMuI * v def _hDeriv_u(self, src, du_dm_v, adjoint=False): @@ -922,9 +852,10 @@ class Fields_j(Fields): :return: product of the derivative of the magnetic field with respect to the field we solved for with a vector """ - if not adjoint: - return -1./(1j*omega(src.freq)) * self._MeMuI * (self._edgeCurl.T * (self._MfRho * du_dm_v) ) - return -1./(1j*omega(src.freq)) * self._MfRho.T * (self._edgeCurl * ( self._MeMuI.T * du_dm_v)) + if adjoint: + return -1./(1j*omega(src.freq)) * self._MfRho.T * (self._edgeCurl * ( self._MeMuI.T * du_dm_v)) + return -1./(1j*omega(src.freq)) * self._MeMuI * (self._edgeCurl.T * (self._MfRho * du_dm_v) ) + def _hDeriv_m(self, src, v, adjoint=False): @@ -958,32 +889,18 @@ class Fields_j(Fields): return hDeriv_m - - def _ePrimary(self, jSolution, srcList): + def _e(self, jSolution, srcList): """ - Primary electric field + Electric field from jSolution :param numpy.ndarray hSolution: field we solved for :param list srcList: list of sources :rtype: numpy.ndarray - :return: primary electric field as defined by the sources - """ - n = int(self._aveF2CCV.shape[0] / self._nC) # number of components - VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) - return VI * (self._aveF2CCV * (self._MfRho * self._jPrimary(jSolution, srcList))) - - def _eSecondary(self, jSolution, srcList): - """ - Secondary electric field from jSolution - - :param numpy.ndarray hSolution: field we solved for - :param list srcList: list of sources - :rtype: numpy.ndarray - :return: secondary electric field + :return: electric field """ n = int(self._aveF2CCV.shape[0] / self._nC) # number of components VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) - return VI * (self._aveF2CCV * (self._MfRho * self._jSecondary(jSolution, srcList))) + return VI * (self._aveF2CCV * (self._MfRho * self._j(jSolution, srcList))) def _eDeriv_u(self, src, du_dm_v, adjoint=False): """ @@ -1018,22 +935,7 @@ class Fields_j(Fields): return self._MfRhoDeriv(jSolution).T * ( self._aveF2CCV.T * ( VI.T * v ) ) return VI * (self._aveF2CCV * (self._MfRhoDeriv(jSolution) * v)) - def _bPrimary(self, jSolution, srcList): - """ - Primary magnetic flux density - - :param numpy.ndarray hSolution: field we solved for - :param list srcList: list of sources - :rtype: numpy.ndarray - :return: primary magnetic flux density - """ - hPrimary = self._hPrimary(jSolution, srcList) - n = int(self._aveE2CCV.shape[0] / self._nC) # number of components - VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) - - return VI * (self._aveE2CCV * (self._MeMu * hPrimary)) - - def _bSecondary(self, jSolution, srcList): + def _b(self, jSolution, srcList): """ Secondary magnetic flux density from jSolution @@ -1045,7 +947,7 @@ class Fields_j(Fields): n = int(self._aveE2CCV.shape[0] / self._nC) # number of components VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) - return VI * (self._aveE2CCV * (self._edgeCurl.T * (self._MfRho * jSolution))) + return VI * (self._aveE2CCV * ( self._MeMu * self._h(jSolution,srcList)) ) def _bDeriv_u(self, src, du_dm_v, adjoint=False): """ @@ -1099,8 +1001,8 @@ class Fields_h(Fields): 'j' : ['hSolution','F','_j'], 'jPrimary' : ['hSolution','F','_jPrimary'], 'jSecondary' : ['hSolution','F','_jSecondary'], - 'e' : ['hSolution','C','_e'], - 'b' : ['hSolution','C','_b'], + 'e' : ['hSolution','CCV','_e'], + 'b' : ['hSolution','CCV','_b'], } def __init__(self,mesh,survey,**kwargs): @@ -1109,8 +1011,10 @@ class Fields_h(Fields): def startup(self): self.prob = self.survey.prob self._edgeCurl = self.survey.prob.mesh.edgeCurl + self._MeMu = self.survey.prob.MeMu self._MeMuI = self.survey.prob.MeMuI self._MfRho = self.survey.prob.MfRho + self._MfRhoDeriv = self.survey.prob.MfRhoDeriv self._rho = self.survey.prob.curModel.rho self._mu = self.survey.prob.curModel.mui self._aveF2CCV = self.survey.prob.mesh.aveF2CCV @@ -1202,7 +1106,7 @@ class Fields_h(Fields): def _jSecondary(self, hSolution, srcList): """ - Secondary current density from eSolution + Secondary current density from hSolution :param numpy.ndarray hSolution: field we solved for :param list srcList: list of sources @@ -1216,39 +1120,9 @@ class Fields_h(Fields): j[:,i] = j[:,i]+ -S_e return j - def _jSecondaryDeriv_u(self, src, du_dm_v, adjoint=False): - """ - Derivative of the secondary current density with respect to the thing we solved for - - :param SimPEG.EM.FDEM.Src src: source - :param numpy.ndarray du_dm_v: vector to take product with - :param bool adjoint: adjoint? - :rtype: numpy.ndarray - :return: product of the derivative of the secondary current density with respect to the field we solved for with a vector - """ - - if not adjoint: - return self._edgeCurl*du_dm_v - elif adjoint: - return self._edgeCurl.T*du_dm_v - - - def _jSecondaryDeriv_m(self, src, v, adjoint=False): - """ - Derivative of the secondary current density with respect to the inversion model. - - :param SimPEG.EM.FDEM.Src src: source - :param numpy.ndarray v: vector to take product with - :param bool adjoint: adjoint? - :rtype: numpy.ndarray - :return: product of the secondary current density derivative with respect to the inversion model with a vector - """ - - _,S_eDeriv = src.evalDeriv(self.prob, v, adjoint) - def _jDeriv_u(self, src, du_dm_v, adjoint=False): """ - Partial derivative of the total current density with respect to the thing we solved for + Derivative of the current density with respect to the thing we solved for :param SimPEG.EM.FDEM.Src src: source :param numpy.ndarray du_dm_v: vector to take product with @@ -1257,47 +1131,113 @@ class Fields_h(Fields): :return: product of the derivative of the current density with respect to the field we solved for with a vector """ - return self._jSecondaryDeriv_u(src,du_dm_v,adjoint) + if not adjoint: + return self._edgeCurl*du_dm_v + elif adjoint: + return self._edgeCurl.T*du_dm_v + def _jDeriv_m(self, src, v, adjoint=False): """ - Partial derivative of the total current density with respect to the inversion model. + Derivative of the current density with respect to the inversion model. :param SimPEG.EM.FDEM.Src src: source :param numpy.ndarray v: vector to take product with :param bool adjoint: adjoint? - :rtype: SimPEG.Utils.Zero - :return: product of the current density with respect to the inversion model with a vector + :rtype: numpy.ndarray + :return: product of the current density derivative with respect to the inversion model with a vector """ - # assuming the primary does not depend on the model - return self._jSecondaryDeriv_m(src,v,adjoint) + _,S_eDeriv = src.evalDeriv(self.prob, v, adjoint) + return -S_eDeriv def _e(self, hSolution, srcList): - rho = self._rho - aveF2CCV = self._aveF2CCV - n = int(self._aveF2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy - - Rho = self.prob.MfRho + """ + Electric field from hSolution + + :param numpy.ndarray hSolution: field we solved for + :param list srcList: list of sources + :rtype: numpy.ndarray + :return: electric field + """ + n = int(self._aveF2CCV.shape[0] / self._nC) #number of components VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) + return VI * (self._aveF2CCV * (self._MfRho * self._j(hSolution, srcList))) + def _eDeriv_u(self, src, du_dm_v, adjoint=False): + """ + Derivative of the electric field with respect to the thing we solved for + + :param SimPEG.EM.FDEM.Src src: source + :param numpy.ndarray du_dm_v: vector to take product with + :param bool adjoint: adjoint? + :rtype: numpy.ndarray + :return: product of the derivative of the electric field with respect to the field we solved for with a vector + """ + n = int(self._aveF2CCV.shape[0] / self._nC) #number of components + VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) + if adjoint: + return self._edgeCurl.T * ( self._MfRho.T * ( self._aveF2CCV.T * ( VI.T * du_dm_v ) ) ) + return VI * (self._aveF2CCV * (self._MfRho * self._edgeCurl * du_dm_v )) - j = self._j(hSolution, srcList) - - return VI * (aveF2CCV * (Rho * j)) - - def _eDeriv_u(self, src, u, v, adjoint=False): - raise NotImplementedError - - def _eDeriv_m(self, src, u, v, adjoint=False): - raise NotImplementedError + def _eDeriv_m(self, src, v, adjoint=False): + """ + Derivative of the electric field with respect to the inversion model. + + :param SimPEG.EM.FDEM.Src src: source + :param numpy.ndarray v: vector to take product with + :param bool adjoint: adjoint? + :rtype: numpy.ndarray + :return: product of the electric field derivative with respect to the inversion model with a vector + """ + hSolution = Utils.mkvc(self[src,'hSolution']) + n = int(self._aveF2CCV.shape[0] / self._nC) #number of components + VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) + if adjoint: + return ( self._MfRhoDeriv(self._edgeCurl * hSolution).T * ( self._aveF2CCV.T * (VI.T * v) ) ) + return VI * (self._aveF2CCV * (self._MfRhoDeriv(self._edgeCurl * hSolution) * v )) def _b(self, hSolution, srcList): + """ + Magnetic flux density from hSolution + + :param numpy.ndarray hSolution: field we solved for + :param list srcList: list of sources + :rtype: numpy.ndarray + :return: magnetic flux density + """ h = self._h(hSolution, srcList) - Mu = self.prob.MeMu - aveE2CCV = self._aveE2CCV - n = int(aveE2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy - # Mu = sdiag(sp.kron(np.ones(n), mu)) + n = int(self._aveE2CCV.shape[0] / self._nC) #number of components VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) - return VI * (aveE2CCV * (Mu * h)) + return VI * (self._aveE2CCV * (self._MeMu * h)) + + def _bDeriv_u(self, src, du_dm_v, adjoint=False): + """ + Derivative of the magnetic flux density with respect to the thing we solved for + + :param SimPEG.EM.FDEM.Src src: source + :param numpy.ndarray du_dm_v: vector to take product with + :param bool adjoint: adjoint? + :rtype: numpy.ndarray + :return: product of the derivative of the magnetic flux density with respect to the field we solved for with a vector + """ + n = int(self._aveE2CCV.shape[0] / self._nC) #number of components + VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) + if adjoint: + return self._MeMu.T * (self._aveE2CCV.T * ( VI.T * du_dm_v )) + return VI * (self._aveE2CCV * (self._MeMu * du_dm_v)) + + def _bDeriv_m(self, src, v, adjoint=False): + """ + Derivative of the magnetic flux density with respect to the inversion model. + + :param SimPEG.EM.FDEM.Src src: source + :param numpy.ndarray v: vector to take product with + :param bool adjoint: adjoint? + :rtype: numpy.ndarray + :return: product of the magnetic flux density derivative with respect to the inversion model with a vector + """ + return Zero() + + diff --git a/tests/em/fdem/forward/test_FDEM_forward.py b/tests/em/fdem/forward/test_FDEM_forward.py index d89612d8..da446675 100644 --- a/tests/em/fdem/forward/test_FDEM_forward.py +++ b/tests/em/fdem/forward/test_FDEM_forward.py @@ -6,7 +6,7 @@ from scipy.constants import mu_0 from SimPEG.EM.Utils.testingUtils import getFDEMProblem, crossCheckTest testEB = True -testHJ = False +testHJ = True testEJ = True testBH = True verbose = False @@ -22,29 +22,29 @@ class FDEM_CrossCheck(unittest.TestCase): if testEB: def test_EB_CrossCheck_exr_Eform(self): self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'exr', verbose=verbose)) - # def test_EB_CrossCheck_eyr_Eform(self): - # self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'eyr', verbose=verbose)) - # def test_EB_CrossCheck_ezr_Eform(self): - # self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'ezr', verbose=verbose)) - # def test_EB_CrossCheck_exi_Eform(self): - # self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'exi', verbose=verbose)) - # def test_EB_CrossCheck_eyi_Eform(self): - # self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'eyi', verbose=verbose)) - # def test_EB_CrossCheck_ezi_Eform(self): - # self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'ezi', verbose=verbose)) + def test_EB_CrossCheck_eyr_Eform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'eyr', verbose=verbose)) + def test_EB_CrossCheck_ezr_Eform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'ezr', verbose=verbose)) + def test_EB_CrossCheck_exi_Eform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'exi', verbose=verbose)) + def test_EB_CrossCheck_eyi_Eform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'eyi', verbose=verbose)) + def test_EB_CrossCheck_ezi_Eform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'ezi', verbose=verbose)) - # def test_EB_CrossCheck_bxr_Eform(self): - # self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'bxr', verbose=verbose)) - # def test_EB_CrossCheck_byr_Eform(self): - # self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'byr', verbose=verbose)) - # def test_EB_CrossCheck_bzr_Eform(self): - # self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'bzr', verbose=verbose)) - # def test_EB_CrossCheck_bxi_Eform(self): - # self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'bxi', verbose=verbose)) - # def test_EB_CrossCheck_byi_Eform(self): - # self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'byi', verbose=verbose)) - # def test_EB_CrossCheck_bzi_Eform(self): - # self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'bzi', verbose=verbose)) + def test_EB_CrossCheck_bxr_Eform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'bxr', verbose=verbose)) + def test_EB_CrossCheck_byr_Eform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'byr', verbose=verbose)) + def test_EB_CrossCheck_bzr_Eform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'bzr', verbose=verbose)) + def test_EB_CrossCheck_bxi_Eform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'bxi', verbose=verbose)) + def test_EB_CrossCheck_byi_Eform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'byi', verbose=verbose)) + def test_EB_CrossCheck_bzi_Eform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'bzi', verbose=verbose)) if testHJ: def test_HJ_CrossCheck_jxr_Jform(self): diff --git a/tests/em/fdem/forward/test_FDEM_forwardHB.py b/tests/em/fdem/forward/test_FDEM_forwardHB.py index 1f8be91c..545a5014 100644 --- a/tests/em/fdem/forward/test_FDEM_forwardHB.py +++ b/tests/em/fdem/forward/test_FDEM_forwardHB.py @@ -5,9 +5,9 @@ import sys from scipy.constants import mu_0 from SimPEG.EM.Utils.testingUtils import getFDEMProblem, crossCheckTest -testEB = False +testEB = True testHJ = True -testEJ = False +testEJ = True testBH = True verbose = False diff --git a/tests/em/fdem/inverse/adjoint/test_FDEM_adjointHJ.py b/tests/em/fdem/inverse/adjoint/test_FDEM_adjointHJ.py index 242b33f2..9853c2ae 100644 --- a/tests/em/fdem/inverse/adjoint/test_FDEM_adjointHJ.py +++ b/tests/em/fdem/inverse/adjoint/test_FDEM_adjointHJ.py @@ -5,8 +5,8 @@ import sys from scipy.constants import mu_0 from SimPEG.EM.Utils.testingUtils import getFDEMProblem -testEB = True -testHJ = True +testJ = True +testH = True verbose = False @@ -46,7 +46,7 @@ def adjointTest(fdemType, comp): class FDEM_AdjointTests(unittest.TestCase): - if testHJ: + if testJ: def test_Jtvec_adjointTest_jxr_Jform(self): self.assertTrue(adjointTest('j', 'jxr')) def test_Jtvec_adjointTest_jyr_Jform(self): @@ -99,31 +99,58 @@ class FDEM_AdjointTests(unittest.TestCase): def test_Jtvec_adjointTest_bzi_Jform(self): self.assertTrue(adjointTest('j', 'bzi')) - # def test_Jtvec_adjointTest_hxr_Hform(self): - # self.assertTrue(adjointTest('h', 'hxr')) - # def test_Jtvec_adjointTest_hyr_Hform(self): - # self.assertTrue(adjointTest('h', 'hyr')) - # def test_Jtvec_adjointTest_hzr_Hform(self): - # self.assertTrue(adjointTest('h', 'hzr')) - # def test_Jtvec_adjointTest_hxi_Hform(self): - # self.assertTrue(adjointTest('h', 'hxi')) - # def test_Jtvec_adjointTest_hyi_Hform(self): - # self.assertTrue(adjointTest('h', 'hyi')) - # def test_Jtvec_adjointTest_hzi_Hform(self): - # self.assertTrue(adjointTest('h', 'hzi')) + if testH: + def test_Jtvec_adjointTest_hxr_Hform(self): + self.assertTrue(adjointTest('h', 'hxr')) + def test_Jtvec_adjointTest_hyr_Hform(self): + self.assertTrue(adjointTest('h', 'hyr')) + def test_Jtvec_adjointTest_hzr_Hform(self): + self.assertTrue(adjointTest('h', 'hzr')) + def test_Jtvec_adjointTest_hxi_Hform(self): + self.assertTrue(adjointTest('h', 'hxi')) + def test_Jtvec_adjointTest_hyi_Hform(self): + self.assertTrue(adjointTest('h', 'hyi')) + def test_Jtvec_adjointTest_hzi_Hform(self): + self.assertTrue(adjointTest('h', 'hzi')) - # def test_Jtvec_adjointTest_hxr_Hform(self): - # self.assertTrue(adjointTest('h', 'jxr')) - # def test_Jtvec_adjointTest_hyr_Hform(self): - # self.assertTrue(adjointTest('h', 'jyr')) - # def test_Jtvec_adjointTest_hzr_Hform(self): - # self.assertTrue(adjointTest('h', 'jzr')) - # def test_Jtvec_adjointTest_hxi_Hform(self): - # self.assertTrue(adjointTest('h', 'jxi')) - # def test_Jtvec_adjointTest_hyi_Hform(self): - # self.assertTrue(adjointTest('h', 'jyi')) - # def test_Jtvec_adjointTest_hzi_Hform(self): - # self.assertTrue(adjointTest('h', 'jzi')) + def test_Jtvec_adjointTest_jxr_Hform(self): + self.assertTrue(adjointTest('h', 'jxr')) + def test_Jtvec_adjointTest_jyr_Hform(self): + self.assertTrue(adjointTest('h', 'jyr')) + def test_Jtvec_adjointTest_jzr_Hform(self): + self.assertTrue(adjointTest('h', 'jzr')) + def test_Jtvec_adjointTest_jxi_Hform(self): + self.assertTrue(adjointTest('h', 'jxi')) + def test_Jtvec_adjointTest_jyi_Hform(self): + self.assertTrue(adjointTest('h', 'jyi')) + def test_Jtvec_adjointTest_jzi_Hform(self): + self.assertTrue(adjointTest('h', 'jzi')) + + def test_Jtvec_adjointTest_exr_Hform(self): + self.assertTrue(adjointTest('h', 'exr')) + def test_Jtvec_adjointTest_eyr_Hform(self): + self.assertTrue(adjointTest('h', 'eyr')) + def test_Jtvec_adjointTest_ezr_Hform(self): + self.assertTrue(adjointTest('h', 'ezr')) + def test_Jtvec_adjointTest_exi_Hform(self): + self.assertTrue(adjointTest('h', 'exi')) + def test_Jtvec_adjointTest_eyi_Hform(self): + self.assertTrue(adjointTest('h', 'eyi')) + def test_Jtvec_adjointTest_ezi_Hform(self): + self.assertTrue(adjointTest('h', 'ezi')) + + def test_Jtvec_adjointTest_bxr_Hform(self): + self.assertTrue(adjointTest('h', 'bxr')) + def test_Jtvec_adjointTest_byr_Hform(self): + self.assertTrue(adjointTest('h', 'byr')) + def test_Jtvec_adjointTest_bzr_Hform(self): + self.assertTrue(adjointTest('h', 'bzr')) + def test_Jtvec_adjointTest_bxi_Hform(self): + self.assertTrue(adjointTest('h', 'bxi')) + def test_Jtvec_adjointTest_byi_Hform(self): + self.assertTrue(adjointTest('h', 'byi')) + def test_Jtvec_adjointTest_bzi_Hform(self): + self.assertTrue(adjointTest('h', 'bzi')) if __name__ == '__main__': From 4df932ccfc5d26e97918b8568f4bb89382360ce1 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sun, 21 Feb 2016 10:35:27 -0800 Subject: [PATCH 13/17] notation cleanup --- SimPEG/EM/FDEM/FieldsFDEM.py | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/SimPEG/EM/FDEM/FieldsFDEM.py b/SimPEG/EM/FDEM/FieldsFDEM.py index 68843d06..74cd1a24 100644 --- a/SimPEG/EM/FDEM/FieldsFDEM.py +++ b/SimPEG/EM/FDEM/FieldsFDEM.py @@ -559,12 +559,12 @@ class Fields_b(Fields): :return: secondary electric field """ - v = ( self._edgeCurl.T * ( self._MfMui * bSolution)) + e = ( self._edgeCurl.T * ( self._MfMui * bSolution)) for i,src in enumerate(srcList): _,S_e = src.eval(self.prob) - v[:,i] = v[:,i] + - S_e + e[:,i] = e[:,i] + - S_e - return self._MeSigmaI * v + return self._MeSigmaI * e def _eDeriv_u(self, src, du_dm_v, adjoint=False): """ @@ -833,12 +833,12 @@ class Fields_j(Fields): :return: secondary magnetic field """ - v = (self._edgeCurl.T * (self._MfRho * jSolution) ) + h = (self._edgeCurl.T * (self._MfRho * jSolution) ) for i, src in enumerate(srcList): h[:,i] *= -1./(1j*omega(src.freq)) S_m,_ = src.eval(self.prob) h[:,i] = h[:,i]+ 1./(1j*omega(src.freq)) * (S_m) - return self._MeMuI * v + return self._MeMuI * h def _hDeriv_u(self, src, du_dm_v, adjoint=False): @@ -962,9 +962,12 @@ class Fields_j(Fields): n = int(self._aveF2CCV.shape[0] / self._nC) # number of components VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) + # if adjoint: + # return self._MfRho.T * ( self._edgeCurl * ( self._aveE2CCV.T * (VI.T * du_dm_v) ) ) + # return VI * (self._aveE2CCV * (self._edgeCurl.T * (self._MfRho * du_dm_v))) if adjoint: - return self._MfRho.T * ( self._edgeCurl * ( self._aveE2CCV.T * (VI.T * du_dm_v) ) ) - return VI * (self._aveE2CCV * (self._edgeCurl.T * (self._MfRho * du_dm_v))) + return self._hDeriv_u(src, self._MeMu.T*(self._aveE2CCV.T * (VI.T * du_dm_v)), adjoint=adjoint) + return VI * (self._aveE2CCV * ( self._MeMu * self._hDeriv_u(src, du_dm_v)) ) def _bDeriv_m(self, src, v, adjoint=False): """ From 6aa9b533ba2240c5e38af5af6dbb48bfba3d93ae Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sun, 21 Feb 2016 10:36:46 -0800 Subject: [PATCH 14/17] include mu in testing, cleanup and debugging in dipole sources --- SimPEG/EM/FDEM/SrcFDEM.py | 14 ++++- SimPEG/EM/Utils/EMUtils.py | 58 +++++++++---------- SimPEG/EM/Utils/__init__.py | 2 +- SimPEG/EM/Utils/testingUtils.py | 35 ++++++----- .../em/fdem/forward/test_FDEM_forwardEJHB.py | 2 +- 5 files changed, 63 insertions(+), 48 deletions(-) diff --git a/SimPEG/EM/FDEM/SrcFDEM.py b/SimPEG/EM/FDEM/SrcFDEM.py index 6193b2f4..41614b0e 100644 --- a/SimPEG/EM/FDEM/SrcFDEM.py +++ b/SimPEG/EM/FDEM/SrcFDEM.py @@ -335,7 +335,7 @@ class MagDipole(BaseSrc): :return: primary magnetic field """ b = self.bPrimary(prob) - return h_from_b(prob,b) + return 1./self.mu * b def S_m(self, prob): """ @@ -347,6 +347,8 @@ class MagDipole(BaseSrc): """ b_p = self.bPrimary(prob) + if prob._eqLocs is 'EF': + b_p = prob.Me * b_p return -1j*omega(self.freq)*b_p def S_e(self, prob): @@ -449,7 +451,7 @@ class MagDipole_Bfield(BaseSrc): :return: primary magnetic field """ b = self.bPrimary(prob) - return h_from_b(prob, b) + return 1/self.mu * b def S_m(self, prob): """ @@ -460,6 +462,8 @@ class MagDipole_Bfield(BaseSrc): :return: primary magnetic field """ b = self.bPrimary(prob) + if prob._eqLocs is 'EF': + b = prob.Me * b return -1j*omega(self.freq)*b def S_e(self, prob): @@ -570,6 +574,8 @@ class CircularLoop(BaseSrc): :return: primary magnetic field """ b = self.bPrimary(prob) + if prob._eqLocs is 'EF': + b = prob.Me * b return -1j*omega(self.freq)*b def S_e(self, prob): @@ -589,6 +595,8 @@ class CircularLoop(BaseSrc): mui_s = prob.curModel.mui - 1./self.mu MMui_s = prob.mesh.getFaceInnerProduct(mui_s) C = prob.mesh.edgeCurl + + elif eqLocs is 'EF': mu_s = prob.curModel.mu - self.mu MMui_s = prob.mesh.getEdgeInnerProduct(mu_s, invMat=True) @@ -596,4 +604,6 @@ class CircularLoop(BaseSrc): return -C.T * (MMui_s * self.bPrimary(prob)) + + diff --git a/SimPEG/EM/Utils/EMUtils.py b/SimPEG/EM/Utils/EMUtils.py index 4a342acb..660ef117 100644 --- a/SimPEG/EM/Utils/EMUtils.py +++ b/SimPEG/EM/Utils/EMUtils.py @@ -13,37 +13,37 @@ def k(freq, sigma, mu=mu_0, eps=epsilon_0): beta = w * np.sqrt( mu*eps/2 * ( np.sqrt(1. + (sigma / (eps*w))**2 ) - 1) ) return alp - 1j*beta -# Constitutive relations -def e_from_j(prob,j): - eqLocs = prob._eqLocs - if eqLocs is 'FE': - MSigmaI = prob.MeSigmaI - elif eqLocs is 'EF': - MSigmaI = prob.MfRho - return MSigmaI*j +# # Constitutive relations +# def e_from_j(prob,j): +# eqLocs = prob._eqLocs +# if eqLocs is 'FE': +# MSigmaI = prob.MeSigmaI +# elif eqLocs is 'EF': +# MSigmaI = prob.MfRho +# return MSigmaI*j -def j_from_e(prob,e): - eqLocs = prob._eqLocs - if eqLocs is 'FE': - MSigma = prob.MeSigma - elif eqLocs is 'EF': - MSigma = prob.MfRhoI - return MSigma*e +# def j_from_e(prob,e): +# eqLocs = prob._eqLocs +# if eqLocs is 'FE': +# MSigma = prob.MeSigma +# elif eqLocs is 'EF': +# MSigma = prob.MfRhoI +# return MSigma*e -def b_from_h(prob,h): - eqLocs = prob._eqLocs - if eqLocs is 'FE': - MMu = prob.MfMuiI - elif eqLocs is 'EF': - MMu = prob.MeMu - return MMu*h +# def b_from_h(prob,h): +# eqLocs = prob._eqLocs +# if eqLocs is 'FE': +# MMu = prob.MfMuiI +# elif eqLocs is 'EF': +# MMu = prob.MeMu +# return MMu*h -def h_from_b(prob,b): - eqLocs = prob._eqLocs - if eqLocs is 'FE': - MMuI = prob.MfMui - elif eqLocs is 'EF': - MMuI = prob.MeMuI - return MMuI*b +# def h_from_b(prob,b): +# eqLocs = prob._eqLocs +# if eqLocs is 'FE': +# MMuI = prob.MfMui +# elif eqLocs is 'EF': +# MMuI = prob.MeMuI +# return MMuI*b diff --git a/SimPEG/EM/Utils/__init__.py b/SimPEG/EM/Utils/__init__.py index 18dddde9..a60d9701 100644 --- a/SimPEG/EM/Utils/__init__.py +++ b/SimPEG/EM/Utils/__init__.py @@ -1,5 +1,5 @@ # import Sources # import Ana # import Solver -from EMUtils import omega, e_from_j, j_from_e, b_from_h, h_from_b +from EMUtils import omega, k from AnalyticUtils import MagneticDipoleFields, MagneticDipoleVectorPotential, MagneticLoopVectorPotential \ No newline at end of file diff --git a/SimPEG/EM/Utils/testingUtils.py b/SimPEG/EM/Utils/testingUtils.py index ae2b7321..569f8e6d 100644 --- a/SimPEG/EM/Utils/testingUtils.py +++ b/SimPEG/EM/Utils/testingUtils.py @@ -8,10 +8,9 @@ FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order CONDUCTIVITY = 1e1 MU = mu_0 freq = 5e-1 -addrandoms = False -def getFDEMProblem(fdemType, comp, SrcList, freq, verbose=False): +def getFDEMProblem(fdemType, comp, SrcList, freq, useMu=False, verbose=False): cs = 10. ncx, ncy, ncz = 0, 0, 0 npad = 8 @@ -20,9 +19,12 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, verbose=False): hz = [(cs,npad,-1.3), (cs,ncz), (cs,npad,1.3)] mesh = Mesh.TensorMesh([hx,hy,hz],['C','C','C']) - mapping = Maps.ExpMap(mesh) + if useMu is True: + mapping = [('sigma', Maps.ExpMap(mesh)), ('mu', Maps.IdentityMap(mesh))] + else: + mapping = Maps.ExpMap(mesh) - x = np.array([np.linspace(-5.*cs,-2.*cs,3),np.linspace(5.*cs,2.*cs,3)]) #don't sample right by the source + x = np.array([np.linspace(-5.*cs,-2.*cs,3),np.linspace(5.*cs,2.*cs,3)]) + cs/4. #don't sample right by the source, slightly off alignment from either staggered grid XYZ = Utils.ndgrid(x,x,np.linspace(-2.*cs,2.*cs,5)) Rx0 = EM.FDEM.Rx(XYZ, comp) @@ -81,22 +83,26 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, verbose=False): return prb -def crossCheckTest(SrcList, fdemType1, fdemType2, comp, TOL=1e-5, verbose=False): +def crossCheckTest(SrcList, fdemType1, fdemType2, comp, addrandoms = False, useMu=False, TOL=1e-5, verbose=False): l2norm = lambda r: np.sqrt(r.dot(r)) - prb1 = getFDEMProblem(fdemType1, comp, SrcList, freq, verbose) + prb1 = getFDEMProblem(fdemType1, comp, SrcList, freq, useMu, verbose) mesh = prb1.mesh print 'Cross Checking Forward: %s, %s formulations - %s' % (fdemType1, fdemType2, comp) - m = np.log(np.ones(mesh.nC)*CONDUCTIVITY) - mu = np.log(np.ones(mesh.nC)*MU) + + logsig = np.log(np.ones(mesh.nC)*CONDUCTIVITY) + mu = np.ones(mesh.nC)*MU if addrandoms is True: - m = m + np.random.randn(mesh.nC)*np.log(CONDUCTIVITY)*1e-1 - mu = mu + np.random.randn(mesh.nC)*MU*1e-1 + logsig += np.random.randn(mesh.nC)*np.log(CONDUCTIVITY)*1e-1 + mu += np.random.randn(mesh.nC)*MU*1e-1 + + if useMu is True: + m = np.r_[logsig, mu] + else: + m = logsig - # prb1.PropMap.PropModel.mu = mu - # prb1.PropMap.PropModel.mui = 1./mu survey1 = prb1.survey d1 = survey1.dpred(m) @@ -104,9 +110,8 @@ def crossCheckTest(SrcList, fdemType1, fdemType2, comp, TOL=1e-5, verbose=False) print ' Problem 1 solved' - prb2 = getFDEMProblem(fdemType2, comp, SrcList, freq, verbose) + prb2 = getFDEMProblem(fdemType2, comp, SrcList, freq, useMu, verbose) - # prb2.mu = mu survey2 = prb2.survey d2 = survey2.dpred(m) @@ -116,6 +121,6 @@ def crossCheckTest(SrcList, fdemType1, fdemType2, comp, TOL=1e-5, verbose=False) r = d2-d1 l2r = l2norm(r) - tol = np.max([TOL*(10**int(np.log10(l2norm(d1)))),FLR]) + tol = np.max([TOL*(10**int(np.log10(0.5* (l2norm(d1) + l2norm(d2)) ))),FLR]) print l2norm(d1), l2norm(d2), l2r , tol, l2r < tol return l2r < tol diff --git a/tests/em/fdem/forward/test_FDEM_forwardEJHB.py b/tests/em/fdem/forward/test_FDEM_forwardEJHB.py index d4218b86..e6319dfc 100644 --- a/tests/em/fdem/forward/test_FDEM_forwardEJHB.py +++ b/tests/em/fdem/forward/test_FDEM_forwardEJHB.py @@ -11,7 +11,7 @@ testBH = True TOLEJHB = 1 # averaging and more sensitive to boundary condition violations (ie. the impact of violating the boundary conditions in each case is different.) #TODO: choose better testing parameters to lower this -SrcList = ['RawVec', 'MagDipole_Bfield', 'MagDipole', 'CircularLoop'] +SrcList = ['RawVec', 'MagDipole', 'MagDipole_Bfield', 'MagDipole', 'CircularLoop'] class FDEM_CrossCheck(unittest.TestCase): From feff936c928b36e3921a07b955821898993a27fc Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sun, 21 Feb 2016 11:02:29 -0800 Subject: [PATCH 15/17] cleanup of tests and debugging derivs --- SimPEG/EM/FDEM/FieldsFDEM.py | 15 +++++++-------- .../fdem/inverse/adjoint/test_FDEM_adjointEB.py | 4 ++-- .../fdem/inverse/adjoint/test_FDEM_adjointHJ.py | 6 ++---- tests/em/fdem/inverse/derivs/test_FDEM_derivs.py | 5 +---- 4 files changed, 12 insertions(+), 18 deletions(-) diff --git a/SimPEG/EM/FDEM/FieldsFDEM.py b/SimPEG/EM/FDEM/FieldsFDEM.py index 74cd1a24..25415126 100644 --- a/SimPEG/EM/FDEM/FieldsFDEM.py +++ b/SimPEG/EM/FDEM/FieldsFDEM.py @@ -837,7 +837,7 @@ class Fields_j(Fields): for i, src in enumerate(srcList): h[:,i] *= -1./(1j*omega(src.freq)) S_m,_ = src.eval(self.prob) - h[:,i] = h[:,i]+ 1./(1j*omega(src.freq)) * (S_m) + h[:,i] = h[:,i] + 1./(1j*omega(src.freq)) * (S_m) return self._MeMuI * h @@ -962,12 +962,9 @@ class Fields_j(Fields): n = int(self._aveF2CCV.shape[0] / self._nC) # number of components VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) - # if adjoint: - # return self._MfRho.T * ( self._edgeCurl * ( self._aveE2CCV.T * (VI.T * du_dm_v) ) ) - # return VI * (self._aveE2CCV * (self._edgeCurl.T * (self._MfRho * du_dm_v))) if adjoint: - return self._hDeriv_u(src, self._MeMu.T*(self._aveE2CCV.T * (VI.T * du_dm_v)), adjoint=adjoint) - return VI * (self._aveE2CCV * ( self._MeMu * self._hDeriv_u(src, du_dm_v)) ) + return -1./(1j*omega(src.freq)) * self._MfRho.T * ( self._edgeCurl * ( self._aveE2CCV.T * (VI.T * du_dm_v) ) ) + return -1./(1j*omega(src.freq)) * VI * (self._aveE2CCV * (self._edgeCurl.T * (self._MfRho * du_dm_v))) def _bDeriv_m(self, src, v, adjoint=False): """ @@ -982,10 +979,12 @@ class Fields_j(Fields): jSolution = self[src,'jSolution'] n = int(self._aveE2CCV.shape[0] / self._nC) # number of components VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) + S_mDeriv,_ = src.evalDeriv(self.prob, adjoint = adjoint) if adjoint: - return self._MfRhoDeriv(jSolution).T * ( self._edgeCurl * (self._aveE2CCV.T * ( VI.T * v)) ) - return VI * (self._aveE2CCV * (self._edgeCurl.T * (self._MfRhoDeriv(jSolution) * v ))) + v = self._aveE2CCV.T * ( VI.T * v) + return 1./(1j * omega(src.freq)) * ( S_mDeriv(v) - self._MfRhoDeriv(jSolution).T * (self._edgeCurl * v )) + return 1./(1j * omega(src.freq)) * VI * (self._aveE2CCV * ( S_mDeriv(v) - self._edgeCurl.T * ( self._MfRhoDeriv(jSolution) * v ) ) ) class Fields_h(Fields): diff --git a/tests/em/fdem/inverse/adjoint/test_FDEM_adjointEB.py b/tests/em/fdem/inverse/adjoint/test_FDEM_adjointEB.py index f8d2dc9d..25762368 100644 --- a/tests/em/fdem/inverse/adjoint/test_FDEM_adjointEB.py +++ b/tests/em/fdem/inverse/adjoint/test_FDEM_adjointEB.py @@ -17,10 +17,10 @@ MU = mu_0 freq = 1e-1 addrandoms = True -SrcType = 'RawVec' #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVec' +SrcList = ['RawVec', 'MagDipole'] #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVec' def adjointTest(fdemType, comp): - prb = getFDEMProblem(fdemType, comp, [SrcType], freq) + prb = getFDEMProblem(fdemType, comp, SrcList, freq) print 'Adjoint %s formulation - %s' % (fdemType, comp) m = np.log(np.ones(prb.mapping.nP)*CONDUCTIVITY) diff --git a/tests/em/fdem/inverse/adjoint/test_FDEM_adjointHJ.py b/tests/em/fdem/inverse/adjoint/test_FDEM_adjointHJ.py index 9853c2ae..c3fb3d37 100644 --- a/tests/em/fdem/inverse/adjoint/test_FDEM_adjointHJ.py +++ b/tests/em/fdem/inverse/adjoint/test_FDEM_adjointHJ.py @@ -17,10 +17,10 @@ MU = mu_0 freq = 1e-1 addrandoms = True -SrcType = 'RawVec' #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVec' +SrcList = ['RawVec', 'MagDipole'] #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVec' def adjointTest(fdemType, comp): - prb = getFDEMProblem(fdemType, comp, [SrcType], freq) + prb = getFDEMProblem(fdemType, comp, SrcList, freq) print 'Adjoint %s formulation - %s' % (fdemType, comp) m = np.log(np.ones(prb.mapping.nP)*CONDUCTIVITY) @@ -31,8 +31,6 @@ def adjointTest(fdemType, comp): mu = mu + np.random.randn(prb.mesh.nC)*MU*1e-1 survey = prb.survey - # prb.PropMap.PropModel.mu = mu - # prb.PropMap.PropModel.mui = 1./mu u = prb.fields(m) v = np.random.rand(survey.nD) diff --git a/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py b/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py index 8c233dfa..0a2e8b82 100644 --- a/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py +++ b/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py @@ -20,7 +20,7 @@ MU = mu_0 freq = 1e-1 addrandoms = True -SrcType = ['RawVec', 'MagDipole'] #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVec' +SrcType = ['MagDipole', 'RawVec'] #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVec' def derivTest(fdemType, comp): @@ -34,9 +34,6 @@ def derivTest(fdemType, comp): x0 = x0 + np.random.randn(prb.mapping.nP)*np.log(CONDUCTIVITY)*1e-1 mu = mu + np.random.randn(prb.mapping.nP)*MU*1e-1 - # prb.PropMap.PropModel.mu = mu - # prb.PropMap.PropModel.mui = 1./mu - survey = prb.survey def fun(x): return survey.dpred(x), lambda x: prb.Jvec(x0, x) From f37235973bf67ea139f0e6bde56303fe600e913e Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sun, 21 Feb 2016 11:33:15 -0800 Subject: [PATCH 16/17] adjoint debugging --- SimPEG/EM/FDEM/FDEM.py | 6 +++--- SimPEG/EM/FDEM/FieldsFDEM.py | 2 ++ 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/SimPEG/EM/FDEM/FDEM.py b/SimPEG/EM/FDEM/FDEM.py index d19f2ad1..f591e9c2 100644 --- a/SimPEG/EM/FDEM/FDEM.py +++ b/SimPEG/EM/FDEM/FDEM.py @@ -139,14 +139,14 @@ class BaseFDEMProblem(BaseEMProblem): dRHS_dmT = self.getRHSDeriv(freq, src, ATinvdf_duT, adjoint=True) du_dmT = -dA_dmT + dRHS_dmT - df_dmT += du_dmT + Df_DmT = df_dmT + du_dmT # TODO: this should be taken care of by the reciever? real_or_imag = rx.projComp if real_or_imag is 'real': - Jtv += np.array(df_dmT,dtype=complex).real + Jtv += np.array(Df_DmT,dtype=complex).real elif real_or_imag is 'imag': - Jtv += - np.array(df_dmT,dtype=complex).real + Jtv += - np.array(Df_DmT,dtype=complex).real else: raise Exception('Must be real or imag') diff --git a/SimPEG/EM/FDEM/FieldsFDEM.py b/SimPEG/EM/FDEM/FieldsFDEM.py index 25415126..37a13f24 100644 --- a/SimPEG/EM/FDEM/FieldsFDEM.py +++ b/SimPEG/EM/FDEM/FieldsFDEM.py @@ -633,6 +633,8 @@ class Fields_b(Fields): """ n = int(self._aveE2CCV.shape[0] / self._nC) # number of components VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) + if adjoint: + return self._MfMui.T * ( self._edgeCurl * ( self._aveE2CCV.T * (VI.T * du_dm_v) ) ) return VI * (self._aveE2CCV * (self._edgeCurl.T * ( self._MfMui * du_dm_v ) ) ) From 5c8fba424251b066d4893a827ea15089e9f0e068 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sat, 27 Feb 2016 15:44:07 -0800 Subject: [PATCH 17/17] changing btwn e,j and h,b depends on which formulation we are going between and is now in the fields object --- SimPEG/EM/Utils/EMUtils.py | 33 --------------------------------- SimPEG/EM/Utils/__init__.py | 3 --- 2 files changed, 36 deletions(-) diff --git a/SimPEG/EM/Utils/EMUtils.py b/SimPEG/EM/Utils/EMUtils.py index 660ef117..e7dbf441 100644 --- a/SimPEG/EM/Utils/EMUtils.py +++ b/SimPEG/EM/Utils/EMUtils.py @@ -13,37 +13,4 @@ def k(freq, sigma, mu=mu_0, eps=epsilon_0): beta = w * np.sqrt( mu*eps/2 * ( np.sqrt(1. + (sigma / (eps*w))**2 ) - 1) ) return alp - 1j*beta -# # Constitutive relations -# def e_from_j(prob,j): -# eqLocs = prob._eqLocs -# if eqLocs is 'FE': -# MSigmaI = prob.MeSigmaI -# elif eqLocs is 'EF': -# MSigmaI = prob.MfRho -# return MSigmaI*j - -# def j_from_e(prob,e): -# eqLocs = prob._eqLocs -# if eqLocs is 'FE': -# MSigma = prob.MeSigma -# elif eqLocs is 'EF': -# MSigma = prob.MfRhoI -# return MSigma*e - -# def b_from_h(prob,h): -# eqLocs = prob._eqLocs -# if eqLocs is 'FE': -# MMu = prob.MfMuiI -# elif eqLocs is 'EF': -# MMu = prob.MeMu -# return MMu*h - -# def h_from_b(prob,b): -# eqLocs = prob._eqLocs -# if eqLocs is 'FE': -# MMuI = prob.MfMui -# elif eqLocs is 'EF': -# MMuI = prob.MeMuI -# return MMuI*b - diff --git a/SimPEG/EM/Utils/__init__.py b/SimPEG/EM/Utils/__init__.py index a60d9701..ef779042 100644 --- a/SimPEG/EM/Utils/__init__.py +++ b/SimPEG/EM/Utils/__init__.py @@ -1,5 +1,2 @@ -# import Sources -# import Ana -# import Solver from EMUtils import omega, k from AnalyticUtils import MagneticDipoleFields, MagneticDipoleVectorPotential, MagneticLoopVectorPotential \ No newline at end of file