From ed1457df86a4d22e30c312bb48525686ce1e7452 Mon Sep 17 00:00:00 2001 From: Lindsey Date: Mon, 23 Feb 2015 17:42:25 -0800 Subject: [PATCH 01/25] start of j implementation, UNTESTED --- simpegEM/Base.py | 17 +++++++++- simpegEM/FDEM/FDEM.py | 69 +++++++++++++++++++++++++++++++++++++++ simpegEM/FDEM/__init__.py | 2 +- 3 files changed, 86 insertions(+), 2 deletions(-) diff --git a/simpegEM/Base.py b/simpegEM/Base.py index 373cee75..da7e4a86 100644 --- a/simpegEM/Base.py +++ b/simpegEM/Base.py @@ -28,6 +28,13 @@ class BaseEMProblem(Problem.BaseProblem): self._MfMui = self.mesh.getFaceInnerProduct(1/mu_0) return self._MfMui + @property + def MeMuI(self): + #TODO: assuming constant mu + if getattr(self, '_MeMuI', None) is None: + self._MeMuI = self.mesh.getEdgeInnerProduct(1/mu_0) + return self._MeMuI + @property def Me(self): if getattr(self, '_Me', None) is None: @@ -50,7 +57,15 @@ class BaseEMProblem(Problem.BaseProblem): self._MeSigmaI = self.mesh.getEdgeInnerProduct(sigma, invMat=True) return self._MeSigmaI - deleteTheseOnModelUpdate = ['_MeSigma', '_MeSigmaI'] + @property + def MfSigmaI(self): + #TODO: hardcoded to sigma as the model + if getattr(self, '_MfSigmaI', None) is None: + sigma = self.curModel.transform + self._MfSigmaI = self.mesh.getFaceInnerProduct(sigma, invMat=True) + return self._MfSigmaI + + deleteTheseOnModelUpdate = ['_MeSigma', '_MeSigmaI','_MfSigmaI'] def fields(self, m): self.curModel = m diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py index c4137486..3f0af1fb 100644 --- a/simpegEM/FDEM/FDEM.py +++ b/simpegEM/FDEM/FDEM.py @@ -335,3 +335,72 @@ class ProblemFDEM_b(BaseFDEMProblem): return None raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) + +class ProblemFDEM_j(BaseFDEMProblem): + """ + Using the H-J formulation of Maxwell's equations + + .. math:: + \\nabla \\times \\sigma^{-1} \\vec{J} + i\\omega\\mu\\vec{H} = 0 + \\nabla \\times \\vec{H} - \\vec{J} = \\vec{J_s} + + Since \(\\vec{J}\) is a flux and \(\\vec{H}\) is a field, we discretize \(\\vec{J}\) on faces and \(\\vec{H}\) on edges. + + For this implementation, we solve for J, using \( \\vec{H} = - (i\\omega\\mu)^{-1} \\nabla \\times \\sigma^{-1} \\vec{J} \) : + + .. math:: + \\nabla \\times ( \\mu^{-1} \\nabla \\times \\sigma^{-1} \\vec{J} ) + i\\omega \\vec{J} = - i\\omega\\vec{J_s} + """ + + solType = 'j' + + def __init__(self, model, **kwargs): + BaseFDEMProblem.__init__(self, model, **kwargs) + + def getA(self, freq): + """ + :param float freq: Frequency + :rtype: scipy.sparse.csr_matrix + :return: A + """ + + muI = self.MeMuI + sigI = self.MfSigmaI + C = self.mesh.edgeCurl + + return C * muI * C.T * sigI + 1j * omega(freq) + + def getADeriv(self, freq, u, v, adjoint=False): + + muI = self.MeMuI + C = self.mesh.edgeCurl + sig = self.curModel.transform + dsig_dm = self.curModel.transformDeriv + #TODO: This only works if diagonal (no tensors)... + dMfSigmaI_dI = - self.MfSigmaI**2 + + dMf_dsig = self.mesh.getFaceInnerProductDeriv(sig)(u) + + if adjoint: + raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) + + return C * ( muI * ( C.T * ( dMfSigmaI_dI * ( dMf_dsig * ( dsig_dm * v ) ) ) ) ) + + + raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) + + + + def getRHS(self, freq): + """ + :param float freq: Frequency + :rtype: numpy.ndarray (nE, nTx) + :return: RHS + """ + raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) + + def calcFields(self, sol, freq, fieldType, adjoint=False): + raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) + + def calcFieldsDeriv(self, sol, freq, fieldType, v, adjoint=False): + raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) diff --git a/simpegEM/FDEM/__init__.py b/simpegEM/FDEM/__init__.py index 70124686..4f7fcf4c 100644 --- a/simpegEM/FDEM/__init__.py +++ b/simpegEM/FDEM/__init__.py @@ -1,2 +1,2 @@ from SurveyFDEM import * -from FDEM import ProblemFDEM_e, ProblemFDEM_b +from FDEM import ProblemFDEM_e, ProblemFDEM_b, ProblemFDEM_j From 2a6f0ab4404dd8ade107ad5dd7ad42d5af1f8907 Mon Sep 17 00:00:00 2001 From: Lindsey Date: Tue, 24 Feb 2015 18:00:34 -0800 Subject: [PATCH 02/25] Start of HJ formulation: FAILING TESTS right now --- simpegEM/Base.py | 23 ++-- simpegEM/FDEM/FDEM.py | 66 ++++++++-- simpegEM/FDEM/SurveyFDEM.py | 16 ++- simpegEM/Tests/test_FDEM.py | 250 ++++++++++++++++++++++-------------- 4 files changed, 234 insertions(+), 121 deletions(-) diff --git a/simpegEM/Base.py b/simpegEM/Base.py index e2298ed0..e0b30203 100644 --- a/simpegEM/Base.py +++ b/simpegEM/Base.py @@ -7,7 +7,7 @@ class BaseEMProblem(Problem.BaseProblem): Problem.BaseProblem.__init__(self, mesh, **kwargs) solType = None - storeTheseFields = ['e', 'b'] + storeTheseFields = ['e', 'b', 'j', 'h'] surveyPair = Survey.BaseSurvey dataPair = Survey.Data @@ -30,6 +30,8 @@ class BaseEMProblem(Problem.BaseProblem): def mu(self, value): if getattr(self, '_MfMui', None) is not None: del self._MfMui + if getattr(self, '_MeMui', None) is not None: + del self._MeMui self._mu = value @@ -45,11 +47,11 @@ class BaseEMProblem(Problem.BaseProblem): return self._MfMui @property - def MeMuI(self): + def MeMui(self): #TODO: assuming constant mu - if getattr(self, '_MeMuI', None) is None: - self._MeMuI = self.mesh.getEdgeInnerProduct(1/mu_0) - return self._MeMuI + if getattr(self, '_MeMui', None) is None: + self._MeMui = self.mesh.getEdgeInnerProduct(1/mu_0) + return self._MeMui @property def Me(self): @@ -74,14 +76,15 @@ class BaseEMProblem(Problem.BaseProblem): return self._MeSigmaI @property - def MfSigmaI(self): + def MfSigmai(self): #TODO: hardcoded to sigma as the model - if getattr(self, '_MfSigmaI', None) is None: + #TODO: hardcoded to sigma diagonal + if getattr(self, '_MfSigmai', None) is None: sigma = self.curModel.transform - self._MfSigmaI = self.mesh.getFaceInnerProduct(sigma, invMat=True) - return self._MfSigmaI + self._MfSigmai = self.mesh.getFaceInnerProduct(1/sigma) + return self._MfSigmai - deleteTheseOnModelUpdate = ['_MeSigma', '_MeSigmaI','_MfSigmaI'] + deleteTheseOnModelUpdate = ['_MeSigma', '_MeSigmaI','_MfSigmai'] def fields(self, m): self.curModel = m diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py index 3f0af1fb..a8181081 100644 --- a/simpegEM/FDEM/FDEM.py +++ b/simpegEM/FDEM/FDEM.py @@ -350,6 +350,10 @@ class ProblemFDEM_j(BaseFDEMProblem): .. math:: \\nabla \\times ( \\mu^{-1} \\nabla \\times \\sigma^{-1} \\vec{J} ) + i\\omega \\vec{J} = - i\\omega\\vec{J_s} + + .. note:: + This implementation does not yet work with full anisotropy!! + """ solType = 'j' @@ -364,43 +368,77 @@ class ProblemFDEM_j(BaseFDEMProblem): :return: A """ - muI = self.MeMuI - sigI = self.MfSigmaI + mui = self.MeMui + sigi = self.MfSigmai C = self.mesh.edgeCurl + iomega = 1j * omega(freq) * sp.eye(self.mesh.nF) - return C * muI * C.T * sigI + 1j * omega(freq) + return C * mui * C.T * sigi + iomega def getADeriv(self, freq, u, v, adjoint=False): - muI = self.MeMuI + mui = self.MeMui C = self.mesh.edgeCurl sig = self.curModel.transform dsig_dm = self.curModel.transformDeriv - #TODO: This only works if diagonal (no tensors)... - dMfSigmaI_dI = - self.MfSigmaI**2 - dMf_dsig = self.mesh.getFaceInnerProductDeriv(sig)(u) + dMfSigi_di = - self.MfSigmai**2 if adjoint: - raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) + return dsig_dm.T * ( dMf_dsig.T * ( dMfSigi_di.T * ( C * ( mui.T * ( C.T * v ) ) ) ) ) - return C * ( muI * ( C.T * ( dMfSigmaI_dI * ( dMf_dsig * ( dsig_dm * v ) ) ) ) ) + return C * ( mui * ( C.T * ( dMfSigi_di * ( dMf_dsig * ( dsig_dm * v ) ) ) ) ) - raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) - - - def getRHS(self, freq): """ :param float freq: Frequency :rtype: numpy.ndarray (nE, nTx) :return: RHS """ - raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) + Txs = self.survey.getTransmitters(freq) + rhs = range(len(Txs)) + for i, tx in enumerate(Txs): + if tx.txType == 'VMD': + src = Sources.MagneticDipoleVectorPotential + SRCx = src(tx.loc, self.mesh.gridFx, 'x') + SRCy = src(tx.loc, self.mesh.gridFy, 'y') + SRCz = src(tx.loc, self.mesh.gridFz, 'z') + + elif tx.txType == 'CircularLoop': + src = Sources.MagneticLoopVectorPotential + SRCx = src(tx.loc, self.mesh.gridFx, 'x', tx.radius) + SRCy = src(tx.loc, self.mesh.gridFy, 'y', tx.radius) + SRCz = src(tx.loc, self.mesh.gridFz, 'z', tx.radius) + else: + raise NotImplemented('%s txType is not implemented' % tx.txType) + rhs[i] = np.concatenate((SRCx, SRCy, SRCz)) + + a = np.concatenate(rhs).reshape((self.mesh.nF, len(Txs)), order='F') + mui = self.MeMui + C = self.mesh.edgeCurl + + j_s = C*mui*C.T*a + return -1j*omega(freq)*j_s def calcFields(self, sol, freq, fieldType, adjoint=False): + j = sol + if fieldType == 'j': + return j + elif fieldType == 'h': + mui = self.MeMui + C = self.mesh.edgeCurl + if not adjoint: + h = -(1./(1j*omega(freq))) * mui * ( C.T * j ) + else: + h = -(1./(1j*omega(freq))) * C * ( mui.T * j ) + return h raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) def calcFieldsDeriv(self, sol, freq, fieldType, v, adjoint=False): + j = sol + if fieldType == 'j': + return None + elif fieldType == 'h': + return None raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) diff --git a/simpegEM/FDEM/SurveyFDEM.py b/simpegEM/FDEM/SurveyFDEM.py index 9e87465f..e9000b7e 100644 --- a/simpegEM/FDEM/SurveyFDEM.py +++ b/simpegEM/FDEM/SurveyFDEM.py @@ -16,6 +16,20 @@ class RxFDEM(Survey.BaseRx): 'bxi':['b', 'Fx', 'imag'], 'byi':['b', 'Fy', 'imag'], 'bzi':['b', 'Fz', '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'], + + 'hxr':['h', 'Ex', 'real'], + 'hyr':['h', 'Ey', 'real'], + 'hzr':['h', 'Ez', 'real'], + 'hxi':['h', 'Ex', 'imag'], + 'hyi':['h', 'Ey', 'imag'], + 'hzi':['h', 'Ez', 'imag'], } radius = None @@ -84,7 +98,7 @@ class TxFDEM(Survey.BaseTx): class FieldsFDEM(Problem.Fields): """Fancy Field Storage for a FDEM survey.""" - knownFields = {'b': 'F', 'e': 'E'} + knownFields = {'b': 'F', 'e': 'E', 'j': 'F', 'h': 'E'} dtype = complex diff --git a/simpegEM/Tests/test_FDEM.py b/simpegEM/Tests/test_FDEM.py index e1d8f1c5..b514f9b2 100644 --- a/simpegEM/Tests/test_FDEM.py +++ b/simpegEM/Tests/test_FDEM.py @@ -27,10 +27,14 @@ def getProblem(fdemType, comp): survey = EM.FDEM.SurveyFDEM([Tx0]) + print fdemType + if fdemType == 'e': prb = EM.FDEM.ProblemFDEM_e(mesh, mapping=mapping) elif fdemType == 'b': prb = EM.FDEM.ProblemFDEM_b(mesh, mapping=mapping) + elif fdemType == 'j': + prb = EM.FDEM.ProblemFDEM_j(mesh, mapping=mapping) else: raise NotImplementedError() prb.pair(survey) @@ -86,109 +90,163 @@ def derivTest(fdemType, comp): class FDEM_DerivTests(unittest.TestCase): - def test_Jvec_exr_Eform(self): - self.assertTrue(derivTest('e', 'exr')) - def test_Jvec_exr_Bform(self): - self.assertTrue(derivTest('b', 'exr')) - def test_Jvec_eyr_Eform(self): - self.assertTrue(derivTest('e', 'eyr')) - def test_Jvec_eyr_Bform(self): - self.assertTrue(derivTest('b', 'eyr')) - def test_Jvec_ezr_Eform(self): - self.assertTrue(derivTest('e', 'ezr')) - def test_Jvec_ezr_Bform(self): - self.assertTrue(derivTest('b', 'ezr')) - def test_Jvec_exi_Eform(self): - self.assertTrue(derivTest('e', 'exi')) - def test_Jvec_exi_Bform(self): - self.assertTrue(derivTest('b', 'exi')) - def test_Jvec_eyi_Eform(self): - self.assertTrue(derivTest('e', 'eyi')) - def test_Jvec_eyi_Bform(self): - self.assertTrue(derivTest('b', 'eyi')) - def test_Jvec_ezi_Eform(self): - self.assertTrue(derivTest('e', 'ezi')) - def test_Jvec_ezi_Bform(self): - self.assertTrue(derivTest('b', 'ezi')) + # def test_Jvec_exr_Eform(self): + # self.assertTrue(derivTest('e', 'exr')) + # def test_Jvec_exr_Bform(self): + # self.assertTrue(derivTest('b', 'exr')) + # def test_Jvec_eyr_Eform(self): + # self.assertTrue(derivTest('e', 'eyr')) + # def test_Jvec_eyr_Bform(self): + # self.assertTrue(derivTest('b', 'eyr')) + # def test_Jvec_ezr_Eform(self): + # self.assertTrue(derivTest('e', 'ezr')) + # def test_Jvec_ezr_Bform(self): + # self.assertTrue(derivTest('b', 'ezr')) + # def test_Jvec_exi_Eform(self): + # self.assertTrue(derivTest('e', 'exi')) + # def test_Jvec_exi_Bform(self): + # self.assertTrue(derivTest('b', 'exi')) + # def test_Jvec_eyi_Eform(self): + # self.assertTrue(derivTest('e', 'eyi')) + # def test_Jvec_eyi_Bform(self): + # self.assertTrue(derivTest('b', 'eyi')) + # def test_Jvec_ezi_Eform(self): + # self.assertTrue(derivTest('e', 'ezi')) + # def test_Jvec_ezi_Bform(self): + # self.assertTrue(derivTest('b', 'ezi')) - def test_Jvec_bxr_Eform(self): - self.assertTrue(derivTest('e', 'bxr')) - def test_Jvec_bxr_Bform(self): - self.assertTrue(derivTest('b', 'bxr')) - def test_Jvec_byr_Eform(self): - self.assertTrue(derivTest('e', 'byr')) - def test_Jvec_byr_Bform(self): - self.assertTrue(derivTest('b', 'byr')) - def test_Jvec_bzr_Eform(self): - self.assertTrue(derivTest('e', 'bzr')) - def test_Jvec_bzr_Bform(self): - self.assertTrue(derivTest('b', 'bzr')) - def test_Jvec_bxi_Eform(self): - self.assertTrue(derivTest('e', 'bxi')) - def test_Jvec_bxi_Bform(self): - self.assertTrue(derivTest('b', 'bxi')) - def test_Jvec_byi_Eform(self): - self.assertTrue(derivTest('e', 'byi')) - def test_Jvec_byi_Bform(self): - self.assertTrue(derivTest('b', 'byi')) - def test_Jvec_bzi_Eform(self): - self.assertTrue(derivTest('e', 'bzi')) - def test_Jvec_bzi_Bform(self): - self.assertTrue(derivTest('b', 'bzi')) + # def test_Jvec_bxr_Eform(self): + # self.assertTrue(derivTest('e', 'bxr')) + # def test_Jvec_bxr_Bform(self): + # self.assertTrue(derivTest('b', 'bxr')) + # def test_Jvec_byr_Eform(self): + # self.assertTrue(derivTest('e', 'byr')) + # def test_Jvec_byr_Bform(self): + # self.assertTrue(derivTest('b', 'byr')) + # def test_Jvec_bzr_Eform(self): + # self.assertTrue(derivTest('e', 'bzr')) + # def test_Jvec_bzr_Bform(self): + # self.assertTrue(derivTest('b', 'bzr')) + # def test_Jvec_bxi_Eform(self): + # self.assertTrue(derivTest('e', 'bxi')) + # def test_Jvec_bxi_Bform(self): + # self.assertTrue(derivTest('b', 'bxi')) + # def test_Jvec_byi_Eform(self): + # self.assertTrue(derivTest('e', 'byi')) + # def test_Jvec_byi_Bform(self): + # self.assertTrue(derivTest('b', 'byi')) + # def test_Jvec_bzi_Eform(self): + # self.assertTrue(derivTest('e', 'bzi')) + # def test_Jvec_bzi_Bform(self): + # self.assertTrue(derivTest('b', 'bzi')) - def test_Jtvec_adjointTest_exr_Eform(self): - self.assertTrue(adjointTest('e', 'exr')) - def test_Jtvec_adjointTest_exr_Bform(self): - self.assertTrue(adjointTest('b', 'exr')) - def test_Jtvec_adjointTest_eyr_Eform(self): - self.assertTrue(adjointTest('e', 'eyr')) - def test_Jtvec_adjointTest_eyr_Bform(self): - self.assertTrue(adjointTest('b', 'eyr')) - def test_Jtvec_adjointTest_ezr_Eform(self): - self.assertTrue(adjointTest('e', 'ezr')) - def test_Jtvec_adjointTest_ezr_Bform(self): - self.assertTrue(adjointTest('b', 'ezr')) - def test_Jtvec_adjointTest_exi_Eform(self): - self.assertTrue(adjointTest('e', 'exi')) - def test_Jtvec_adjointTest_exi_Bform(self): - self.assertTrue(adjointTest('b', 'exi')) - def test_Jtvec_adjointTest_eyi_Eform(self): - self.assertTrue(adjointTest('e', 'eyi')) - def test_Jtvec_adjointTest_eyi_Bform(self): - self.assertTrue(adjointTest('b', 'eyi')) - def test_Jtvec_adjointTest_ezi_Eform(self): - self.assertTrue(adjointTest('e', 'ezi')) - def test_Jtvec_adjointTest_ezi_Bform(self): - self.assertTrue(adjointTest('b', 'ezi')) + # def test_Jtvec_adjointTest_exr_Eform(self): + # self.assertTrue(adjointTest('e', 'exr')) + # def test_Jtvec_adjointTest_exr_Bform(self): + # self.assertTrue(adjointTest('b', 'exr')) + # def test_Jtvec_adjointTest_eyr_Eform(self): + # self.assertTrue(adjointTest('e', 'eyr')) + # def test_Jtvec_adjointTest_eyr_Bform(self): + # self.assertTrue(adjointTest('b', 'eyr')) + # def test_Jtvec_adjointTest_ezr_Eform(self): + # self.assertTrue(adjointTest('e', 'ezr')) + # def test_Jtvec_adjointTest_ezr_Bform(self): + # self.assertTrue(adjointTest('b', 'ezr')) + # def test_Jtvec_adjointTest_exi_Eform(self): + # self.assertTrue(adjointTest('e', 'exi')) + # def test_Jtvec_adjointTest_exi_Bform(self): + # self.assertTrue(adjointTest('b', 'exi')) + # def test_Jtvec_adjointTest_eyi_Eform(self): + # self.assertTrue(adjointTest('e', 'eyi')) + # def test_Jtvec_adjointTest_eyi_Bform(self): + # self.assertTrue(adjointTest('b', 'eyi')) + # def test_Jtvec_adjointTest_ezi_Eform(self): + # self.assertTrue(adjointTest('e', 'ezi')) + # def test_Jtvec_adjointTest_ezi_Bform(self): + # self.assertTrue(adjointTest('b', 'ezi')) - def test_Jtvec_adjointTest_bxr_Eform(self): - self.assertTrue(adjointTest('e', 'bxr')) - def test_Jtvec_adjointTest_bxr_Bform(self): - self.assertTrue(adjointTest('b', 'bxr')) - def test_Jtvec_adjointTest_byr_Eform(self): - self.assertTrue(adjointTest('e', 'byr')) - def test_Jtvec_adjointTest_byr_Bform(self): - self.assertTrue(adjointTest('b', 'byr')) - def test_Jtvec_adjointTest_bzr_Eform(self): - self.assertTrue(adjointTest('e', 'bzr')) - def test_Jtvec_adjointTest_bzr_Bform(self): - self.assertTrue(adjointTest('b', 'bzr')) - def test_Jtvec_adjointTest_bxi_Eform(self): - self.assertTrue(adjointTest('e', 'bxi')) - def test_Jtvec_adjointTest_bxi_Bform(self): - self.assertTrue(adjointTest('b', 'bxi')) - def test_Jtvec_adjointTest_byi_Eform(self): - self.assertTrue(adjointTest('e', 'byi')) - def test_Jtvec_adjointTest_byi_Bform(self): - self.assertTrue(adjointTest('b', 'byi')) - def test_Jtvec_adjointTest_bzi_Eform(self): - self.assertTrue(adjointTest('e', 'bzi')) - def test_Jtvec_adjointTest_bzi_Bform(self): - self.assertTrue(adjointTest('b', 'bzi')) + # def test_Jtvec_adjointTest_bxr_Eform(self): + # self.assertTrue(adjointTest('e', 'bxr')) + # def test_Jtvec_adjointTest_bxr_Bform(self): + # self.assertTrue(adjointTest('b', 'bxr')) + # def test_Jtvec_adjointTest_byr_Eform(self): + # self.assertTrue(adjointTest('e', 'byr')) + # def test_Jtvec_adjointTest_byr_Bform(self): + # self.assertTrue(adjointTest('b', 'byr')) + # def test_Jtvec_adjointTest_bzr_Eform(self): + # self.assertTrue(adjointTest('e', 'bzr')) + # def test_Jtvec_adjointTest_bzr_Bform(self): + # self.assertTrue(adjointTest('b', 'bzr')) + # def test_Jtvec_adjointTest_bxi_Eform(self): + # self.assertTrue(adjointTest('e', 'bxi')) + # def test_Jtvec_adjointTest_bxi_Bform(self): + # self.assertTrue(adjointTest('b', 'bxi')) + # def test_Jtvec_adjointTest_byi_Eform(self): + # self.assertTrue(adjointTest('e', 'byi')) + # def test_Jtvec_adjointTest_byi_Bform(self): + # self.assertTrue(adjointTest('b', 'byi')) + # def test_Jtvec_adjointTest_bzi_Eform(self): + # self.assertTrue(adjointTest('e', 'bzi')) + # def test_Jtvec_adjointTest_bzi_Bform(self): + # self.assertTrue(adjointTest('b', 'bzi')) + # def test_Jvec_jxr_Jform(self): + # self.assertTrue(derivTest('j', 'jxr')) + # def test_Jvec_jyr_Jform(self): + # self.assertTrue(derivTest('j', 'jyr')) + # def test_Jvec_jzr_Jform(self): + # self.assertTrue(derivTest('j', 'jzr')) + # def test_Jvec_jxi_Jform(self): + # self.assertTrue(derivTest('j', 'jxi')) + # def test_Jvec_jyi_Jform(self): + # self.assertTrue(derivTest('j', 'jyi')) + # def test_Jvec_jzi_Jform(self): + # self.assertTrue(derivTest('j', 'jzi')) + + # def test_Jvec_hxr_Jform(self): + # self.assertTrue(derivTest('j', 'hxr')) + # def test_Jvec_hyr_Jform(self): + # self.assertTrue(derivTest('j', 'hyr')) + # def test_Jvec_hzr_Jform(self): + # self.assertTrue(derivTest('j', 'hzr')) + # def test_Jvec_hxi_Jform(self): + # self.assertTrue(derivTest('j', 'hxi')) + # def test_Jvec_hyi_Jform(self): + # self.assertTrue(derivTest('j', 'hyi')) + # def test_Jvec_hzi_Jform(self): + # self.assertTrue(derivTest('j', 'hzi')) + + + + # def test_Jtvec_adjointTest_jxr_Jform(self): + # self.assertTrue(adjointTest('j', 'jxr')) + # def test_Jtvec_adjointTest_jyr_Jform(self): + # self.assertTrue(adjointTest('j', 'jyr')) + # def test_Jtvec_adjointTest_jzr_Jform(self): + # self.assertTrue(adjointTest('j', 'jzr')) + # def test_Jtvec_adjointTest_jxi_Jform(self): + # self.assertTrue(adjointTest('j', 'jxi')) + # def test_Jtvec_adjointTest_jyi_Jform(self): + # self.assertTrue(adjointTest('j', 'jyi')) + # def test_Jtvec_adjointTest_jzi_Jform(self): + # self.assertTrue(adjointTest('j', 'jzi'))R + + def test_Jtvec_adjointTest_hxr_Jform(self): + self.assertTrue(adjointTest('j', 'hxr')) + # def test_Jtvec_adjointTest_hyr_Jform(self): + # self.assertTrue(adjointTest('j', 'hyr')) + # def test_Jtvec_adjointTest_hzr_Jform(self): + # self.assertTrue(adjointTest('j', 'hzr')) + # def test_Jtvec_adjointTest_hxi_Jform(self): + # self.assertTrue(adjointTest('j', 'hxi')) + # def test_Jtvec_adjointTest_hyi_Jform(self): + # self.assertTrue(adjointTest('j', 'hyi')) + # def test_Jtvec_adjointTest_hzi_Jform(self): + # self.assertTrue(adjointTest('j', 'hzi')) + if __name__ == '__main__': unittest.main() From 57380ea2d8ed53ef9929844adf66b79a447e3523 Mon Sep 17 00:00:00 2001 From: Lindsey Date: Wed, 25 Feb 2015 16:28:54 -0800 Subject: [PATCH 03/25] Tested and passing HJ formulation, solving for J. Only works with diagonal anisotropy --- simpegEM/Base.py | 2 +- simpegEM/FDEM/FDEM.py | 37 ++++- simpegEM/Tests/test_FDEM.py | 318 ++++++++++++++++++------------------ 3 files changed, 189 insertions(+), 168 deletions(-) diff --git a/simpegEM/Base.py b/simpegEM/Base.py index e0b30203..e0373762 100644 --- a/simpegEM/Base.py +++ b/simpegEM/Base.py @@ -7,7 +7,7 @@ class BaseEMProblem(Problem.BaseProblem): Problem.BaseProblem.__init__(self, mesh, **kwargs) solType = None - storeTheseFields = ['e', 'b', 'j', 'h'] + storeTheseFields = ['e', 'b'] surveyPair = Survey.BaseSurvey dataPair = Survey.Data diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py index a8181081..c6304203 100644 --- a/simpegEM/FDEM/FDEM.py +++ b/simpegEM/FDEM/FDEM.py @@ -105,6 +105,9 @@ class BaseFDEMProblem(BaseEMProblem): return Jtv +########################################################################################## +################################ E-B Formulation ######################################### +########################################################################################## class ProblemFDEM_e(BaseFDEMProblem): """ @@ -335,6 +338,10 @@ class ProblemFDEM_b(BaseFDEMProblem): return None raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) +########################################################################################## +################################ H-J Formulation ######################################### +########################################################################################## + class ProblemFDEM_j(BaseFDEMProblem): """ @@ -360,6 +367,7 @@ class ProblemFDEM_j(BaseFDEMProblem): def __init__(self, model, **kwargs): BaseFDEMProblem.__init__(self, model, **kwargs) + BaseFDEMProblem.storeTheseFields = ['j','h'] def getA(self, freq): """ @@ -377,17 +385,18 @@ class ProblemFDEM_j(BaseFDEMProblem): def getADeriv(self, freq, u, v, adjoint=False): - mui = self.MeMui + MeMui = self.MeMui C = self.mesh.edgeCurl sig = self.curModel.transform + sigi = 1/sig dsig_dm = self.curModel.transformDeriv - dMf_dsig = self.mesh.getFaceInnerProductDeriv(sig)(u) - dMfSigi_di = - self.MfSigmai**2 + dsigi_dsig = -Utils.sdiag(sigi)**2 + dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(u) if adjoint: - return dsig_dm.T * ( dMf_dsig.T * ( dMfSigi_di.T * ( C * ( mui.T * ( C.T * v ) ) ) ) ) + return dsig_dm.T * ( dsigi_dsig.T *( dMf_dsigi.T * ( C * ( MeMui.T * ( C.T * v ) ) ) ) ) - return C * ( mui * ( C.T * ( dMfSigi_di * ( dMf_dsig * ( dsig_dm * v ) ) ) ) ) + return C * ( MeMui * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) ) ) def getRHS(self, freq): @@ -428,10 +437,11 @@ class ProblemFDEM_j(BaseFDEMProblem): elif fieldType == 'h': mui = self.MeMui C = self.mesh.edgeCurl + MfSigi = self.MfSigmai if not adjoint: - h = -(1./(1j*omega(freq))) * mui * ( C.T * j ) + h = -(1./(1j*omega(freq))) * mui * ( C.T * ( MfSigi * j ) ) else: - h = -(1./(1j*omega(freq))) * C * ( mui.T * j ) + h = -(1./(1j*omega(freq))) * MfSigi * ( C * ( mui.T * j ) ) return h raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) @@ -440,5 +450,16 @@ class ProblemFDEM_j(BaseFDEMProblem): if fieldType == 'j': return None elif fieldType == 'h': - return None + MeMui = self.MeMui + C = self.mesh.edgeCurl + sig = self.curModel.transform + sigi = 1/sig + dsig_dm = self.curModel.transformDeriv + dsigi_dsig = -Utils.sdiag(sigi)**2 + dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(j) + sigi = self.MfSigmai + if not adjoint: + return -(1./(1j*omega(freq))) * MeMui * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) ) + else: + return -(1./(1j*omega(freq))) * dsig_dm.T * ( dsigi_dsig.T * ( dMf_dsigi.T * ( C * ( MeMui.T * v ) ) ) ) raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) diff --git a/simpegEM/Tests/test_FDEM.py b/simpegEM/Tests/test_FDEM.py index b514f9b2..2091a43d 100644 --- a/simpegEM/Tests/test_FDEM.py +++ b/simpegEM/Tests/test_FDEM.py @@ -5,9 +5,10 @@ import sys from scipy.constants import mu_0 TOL = 1e-4 -CONDUCTIVITY = 1e3 +FLR = 1e-15 # "zero", so if residual below this --> pass regardless of order +CONDUCTIVITY = 1e1 MU = mu_0 -addrandoms = True +addrandoms = True # important to addrandoms if testing HJ formulation with VMD source! (or else jz ~ 0) def getProblem(fdemType, comp): cs = 5. @@ -23,7 +24,7 @@ def getProblem(fdemType, comp): x = np.linspace(-30,30,6) XYZ = Utils.ndgrid(x,x,np.r_[0]) Rx0 = EM.FDEM.RxFDEM(XYZ, comp) - Tx0 = EM.FDEM.TxFDEM(np.r_[4.,2.,2.], 'VMD', 1e-2, [Rx0]) + Tx0 = EM.FDEM.TxFDEM(np.r_[0.,0.,0.], 'VMD', 1, [Rx0]) survey = EM.FDEM.SurveyFDEM([Tx0]) @@ -39,11 +40,11 @@ def getProblem(fdemType, comp): raise NotImplementedError() prb.pair(survey) - try: - from pymatsolver import MumpsSolver - prb.Solver = MumpsSolver - except ImportError, e: - pass + # try: + # from pymatsolver import MumpsSolver + # prb.Solver = MumpsSolver + # except ImportError, e: + # pass return prb @@ -55,19 +56,20 @@ def adjointTest(fdemType, comp): mu = np.log(np.ones(prb.mesh.nC)*MU) if addrandoms is True: - m = m + np.random.randn(prb.mesh.nC)*CONDUCTIVITY*1e-3 - mu = mu + np.random.randn(prb.mesh.nC)*MU*1e-3 + m = m + np.random.randn(prb.mesh.nC)*CONDUCTIVITY*1e-1 + mu = mu + np.random.randn(prb.mesh.nC)*MU*1e-1 prb.mu = mu survey = prb.survey + u = prb.fields(m) + v = np.random.rand(survey.nD) w = np.random.rand(prb.mapping.nP) - u = prb.fields(m) vJw = v.dot(prb.Jvec(m, w, u=u)) wJtv = w.dot(prb.Jtvec(m, v, u=u)) - tol = TOL*(10**int(np.log10(np.abs(vJw)))) + 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 @@ -78,174 +80,172 @@ def derivTest(fdemType, comp): mu = np.log(np.ones(prb.mesh.nC)*MU) if addrandoms is True: - x0 = x0 + np.random.randn(prb.mesh.nC)*CONDUCTIVITY*1e-3 - mu = mu + np.random.randn(prb.mesh.nC)*MU*1e-3 + x0 = x0 + np.random.randn(prb.mesh.nC)*CONDUCTIVITY*1e-1 + mu = mu + np.random.randn(prb.mesh.nC)*MU*1e-1 prb.mu = mu survey = prb.survey def fun(x): return survey.dpred(x), lambda x: prb.Jvec(x0, x) - return Tests.checkDerivative(fun, x0, num=3, plotIt=False, eps=1e-25) + return Tests.checkDerivative(fun, x0, num=3, plotIt=False, eps=FLR) class FDEM_DerivTests(unittest.TestCase): - # def test_Jvec_exr_Eform(self): - # self.assertTrue(derivTest('e', 'exr')) - # def test_Jvec_exr_Bform(self): - # self.assertTrue(derivTest('b', 'exr')) - # def test_Jvec_eyr_Eform(self): - # self.assertTrue(derivTest('e', 'eyr')) - # def test_Jvec_eyr_Bform(self): - # self.assertTrue(derivTest('b', 'eyr')) - # def test_Jvec_ezr_Eform(self): - # self.assertTrue(derivTest('e', 'ezr')) - # def test_Jvec_ezr_Bform(self): - # self.assertTrue(derivTest('b', 'ezr')) - # def test_Jvec_exi_Eform(self): - # self.assertTrue(derivTest('e', 'exi')) - # def test_Jvec_exi_Bform(self): - # self.assertTrue(derivTest('b', 'exi')) - # def test_Jvec_eyi_Eform(self): - # self.assertTrue(derivTest('e', 'eyi')) - # def test_Jvec_eyi_Bform(self): - # self.assertTrue(derivTest('b', 'eyi')) - # def test_Jvec_ezi_Eform(self): - # self.assertTrue(derivTest('e', 'ezi')) - # def test_Jvec_ezi_Bform(self): - # self.assertTrue(derivTest('b', 'ezi')) + def test_Jvec_exr_Eform(self): + self.assertTrue(derivTest('e', 'exr')) + def test_Jvec_exr_Bform(self): + self.assertTrue(derivTest('b', 'exr')) + def test_Jvec_eyr_Eform(self): + self.assertTrue(derivTest('e', 'eyr')) + def test_Jvec_eyr_Bform(self): + self.assertTrue(derivTest('b', 'eyr')) + def test_Jvec_ezr_Eform(self): + self.assertTrue(derivTest('e', 'ezr')) + def test_Jvec_ezr_Bform(self): + self.assertTrue(derivTest('b', 'ezr')) + def test_Jvec_exi_Eform(self): + self.assertTrue(derivTest('e', 'exi')) + def test_Jvec_exi_Bform(self): + self.assertTrue(derivTest('b', 'exi')) + def test_Jvec_eyi_Eform(self): + self.assertTrue(derivTest('e', 'eyi')) + def test_Jvec_eyi_Bform(self): + self.assertTrue(derivTest('b', 'eyi')) + def test_Jvec_ezi_Eform(self): + self.assertTrue(derivTest('e', 'ezi')) + def test_Jvec_ezi_Bform(self): + self.assertTrue(derivTest('b', 'ezi')) - # def test_Jvec_bxr_Eform(self): - # self.assertTrue(derivTest('e', 'bxr')) - # def test_Jvec_bxr_Bform(self): - # self.assertTrue(derivTest('b', 'bxr')) - # def test_Jvec_byr_Eform(self): - # self.assertTrue(derivTest('e', 'byr')) - # def test_Jvec_byr_Bform(self): - # self.assertTrue(derivTest('b', 'byr')) - # def test_Jvec_bzr_Eform(self): - # self.assertTrue(derivTest('e', 'bzr')) - # def test_Jvec_bzr_Bform(self): - # self.assertTrue(derivTest('b', 'bzr')) - # def test_Jvec_bxi_Eform(self): - # self.assertTrue(derivTest('e', 'bxi')) - # def test_Jvec_bxi_Bform(self): - # self.assertTrue(derivTest('b', 'bxi')) - # def test_Jvec_byi_Eform(self): - # self.assertTrue(derivTest('e', 'byi')) - # def test_Jvec_byi_Bform(self): - # self.assertTrue(derivTest('b', 'byi')) - # def test_Jvec_bzi_Eform(self): - # self.assertTrue(derivTest('e', 'bzi')) - # def test_Jvec_bzi_Bform(self): - # self.assertTrue(derivTest('b', 'bzi')) + def test_Jvec_bxr_Eform(self): + self.assertTrue(derivTest('e', 'bxr')) + def test_Jvec_bxr_Bform(self): + self.assertTrue(derivTest('b', 'bxr')) + def test_Jvec_byr_Eform(self): + self.assertTrue(derivTest('e', 'byr')) + def test_Jvec_byr_Bform(self): + self.assertTrue(derivTest('b', 'byr')) + def test_Jvec_bzr_Eform(self): + self.assertTrue(derivTest('e', 'bzr')) + def test_Jvec_bzr_Bform(self): + self.assertTrue(derivTest('b', 'bzr')) + def test_Jvec_bxi_Eform(self): + self.assertTrue(derivTest('e', 'bxi')) + def test_Jvec_bxi_Bform(self): + self.assertTrue(derivTest('b', 'bxi')) + def test_Jvec_byi_Eform(self): + self.assertTrue(derivTest('e', 'byi')) + def test_Jvec_byi_Bform(self): + self.assertTrue(derivTest('b', 'byi')) + def test_Jvec_bzi_Eform(self): + self.assertTrue(derivTest('e', 'bzi')) + def test_Jvec_bzi_Bform(self): + self.assertTrue(derivTest('b', 'bzi')) - # def test_Jtvec_adjointTest_exr_Eform(self): - # self.assertTrue(adjointTest('e', 'exr')) - # def test_Jtvec_adjointTest_exr_Bform(self): - # self.assertTrue(adjointTest('b', 'exr')) - # def test_Jtvec_adjointTest_eyr_Eform(self): - # self.assertTrue(adjointTest('e', 'eyr')) - # def test_Jtvec_adjointTest_eyr_Bform(self): - # self.assertTrue(adjointTest('b', 'eyr')) - # def test_Jtvec_adjointTest_ezr_Eform(self): - # self.assertTrue(adjointTest('e', 'ezr')) - # def test_Jtvec_adjointTest_ezr_Bform(self): - # self.assertTrue(adjointTest('b', 'ezr')) - # def test_Jtvec_adjointTest_exi_Eform(self): - # self.assertTrue(adjointTest('e', 'exi')) - # def test_Jtvec_adjointTest_exi_Bform(self): - # self.assertTrue(adjointTest('b', 'exi')) - # def test_Jtvec_adjointTest_eyi_Eform(self): - # self.assertTrue(adjointTest('e', 'eyi')) - # def test_Jtvec_adjointTest_eyi_Bform(self): - # self.assertTrue(adjointTest('b', 'eyi')) - # def test_Jtvec_adjointTest_ezi_Eform(self): - # self.assertTrue(adjointTest('e', 'ezi')) - # def test_Jtvec_adjointTest_ezi_Bform(self): - # self.assertTrue(adjointTest('b', 'ezi')) + def test_Jtvec_adjointTest_exr_Eform(self): + self.assertTrue(adjointTest('e', 'exr')) + def test_Jtvec_adjointTest_exr_Bform(self): + self.assertTrue(adjointTest('b', 'exr')) + def test_Jtvec_adjointTest_eyr_Eform(self): + self.assertTrue(adjointTest('e', 'eyr')) + def test_Jtvec_adjointTest_eyr_Bform(self): + self.assertTrue(adjointTest('b', 'eyr')) + def test_Jtvec_adjointTest_ezr_Eform(self): + self.assertTrue(adjointTest('e', 'ezr')) + def test_Jtvec_adjointTest_ezr_Bform(self): + self.assertTrue(adjointTest('b', 'ezr')) + def test_Jtvec_adjointTest_exi_Eform(self): + self.assertTrue(adjointTest('e', 'exi')) + def test_Jtvec_adjointTest_exi_Bform(self): + self.assertTrue(adjointTest('b', 'exi')) + def test_Jtvec_adjointTest_eyi_Eform(self): + self.assertTrue(adjointTest('e', 'eyi')) + def test_Jtvec_adjointTest_eyi_Bform(self): + self.assertTrue(adjointTest('b', 'eyi')) + def test_Jtvec_adjointTest_ezi_Eform(self): + self.assertTrue(adjointTest('e', 'ezi')) + def test_Jtvec_adjointTest_ezi_Bform(self): + self.assertTrue(adjointTest('b', 'ezi')) - # def test_Jtvec_adjointTest_bxr_Eform(self): - # self.assertTrue(adjointTest('e', 'bxr')) - # def test_Jtvec_adjointTest_bxr_Bform(self): - # self.assertTrue(adjointTest('b', 'bxr')) - # def test_Jtvec_adjointTest_byr_Eform(self): - # self.assertTrue(adjointTest('e', 'byr')) - # def test_Jtvec_adjointTest_byr_Bform(self): - # self.assertTrue(adjointTest('b', 'byr')) - # def test_Jtvec_adjointTest_bzr_Eform(self): - # self.assertTrue(adjointTest('e', 'bzr')) - # def test_Jtvec_adjointTest_bzr_Bform(self): - # self.assertTrue(adjointTest('b', 'bzr')) - # def test_Jtvec_adjointTest_bxi_Eform(self): - # self.assertTrue(adjointTest('e', 'bxi')) - # def test_Jtvec_adjointTest_bxi_Bform(self): - # self.assertTrue(adjointTest('b', 'bxi')) - # def test_Jtvec_adjointTest_byi_Eform(self): - # self.assertTrue(adjointTest('e', 'byi')) - # def test_Jtvec_adjointTest_byi_Bform(self): - # self.assertTrue(adjointTest('b', 'byi')) - # def test_Jtvec_adjointTest_bzi_Eform(self): - # self.assertTrue(adjointTest('e', 'bzi')) - # def test_Jtvec_adjointTest_bzi_Bform(self): - # self.assertTrue(adjointTest('b', 'bzi')) + def test_Jtvec_adjointTest_bxr_Eform(self): + self.assertTrue(adjointTest('e', 'bxr')) + def test_Jtvec_adjointTest_bxr_Bform(self): + self.assertTrue(adjointTest('b', 'bxr')) + def test_Jtvec_adjointTest_byr_Eform(self): + self.assertTrue(adjointTest('e', 'byr')) + def test_Jtvec_adjointTest_byr_Bform(self): + self.assertTrue(adjointTest('b', 'byr')) + def test_Jtvec_adjointTest_bzr_Eform(self): + self.assertTrue(adjointTest('e', 'bzr')) + def test_Jtvec_adjointTest_bzr_Bform(self): + self.assertTrue(adjointTest('b', 'bzr')) + def test_Jtvec_adjointTest_bxi_Eform(self): + self.assertTrue(adjointTest('e', 'bxi')) + def test_Jtvec_adjointTest_bxi_Bform(self): + self.assertTrue(adjointTest('b', 'bxi')) + def test_Jtvec_adjointTest_byi_Eform(self): + self.assertTrue(adjointTest('e', 'byi')) + def test_Jtvec_adjointTest_byi_Bform(self): + self.assertTrue(adjointTest('b', 'byi')) + def test_Jtvec_adjointTest_bzi_Eform(self): + self.assertTrue(adjointTest('e', 'bzi')) + def test_Jtvec_adjointTest_bzi_Bform(self): + self.assertTrue(adjointTest('b', 'bzi')) - # def test_Jvec_jxr_Jform(self): - # self.assertTrue(derivTest('j', 'jxr')) - # def test_Jvec_jyr_Jform(self): - # self.assertTrue(derivTest('j', 'jyr')) - # def test_Jvec_jzr_Jform(self): - # self.assertTrue(derivTest('j', 'jzr')) - # def test_Jvec_jxi_Jform(self): - # self.assertTrue(derivTest('j', 'jxi')) - # def test_Jvec_jyi_Jform(self): - # self.assertTrue(derivTest('j', 'jyi')) - # def test_Jvec_jzi_Jform(self): - # self.assertTrue(derivTest('j', 'jzi')) + def test_Jvec_jxr_Jform(self): + self.assertTrue(derivTest('j', 'jxr')) + def test_Jvec_jyr_Jform(self): + self.assertTrue(derivTest('j', 'jyr')) + def test_Jvec_jzr_Jform(self): + self.assertTrue(derivTest('j', 'jzr')) + def test_Jvec_jxi_Jform(self): + self.assertTrue(derivTest('j', 'jxi')) + def test_Jvec_jyi_Jform(self): + self.assertTrue(derivTest('j', 'jyi')) + def test_Jvec_jzi_Jform(self): + self.assertTrue(derivTest('j', 'jzi')) - # def test_Jvec_hxr_Jform(self): - # self.assertTrue(derivTest('j', 'hxr')) - # def test_Jvec_hyr_Jform(self): - # self.assertTrue(derivTest('j', 'hyr')) - # def test_Jvec_hzr_Jform(self): - # self.assertTrue(derivTest('j', 'hzr')) - # def test_Jvec_hxi_Jform(self): - # self.assertTrue(derivTest('j', 'hxi')) - # def test_Jvec_hyi_Jform(self): - # self.assertTrue(derivTest('j', 'hyi')) - # def test_Jvec_hzi_Jform(self): - # self.assertTrue(derivTest('j', 'hzi')) + def test_Jvec_hxr_Jform(self): + self.assertTrue(derivTest('j', 'hxr')) + def test_Jvec_hyr_Jform(self): + self.assertTrue(derivTest('j', 'hyr')) + def test_Jvec_hzr_Jform(self): + self.assertTrue(derivTest('j', 'hzr')) + def test_Jvec_hxi_Jform(self): + self.assertTrue(derivTest('j', 'hxi')) + def test_Jvec_hyi_Jform(self): + self.assertTrue(derivTest('j', 'hyi')) + def test_Jvec_hzi_Jform(self): + self.assertTrue(derivTest('j', 'hzi')) - - - # def test_Jtvec_adjointTest_jxr_Jform(self): - # self.assertTrue(adjointTest('j', 'jxr')) - # def test_Jtvec_adjointTest_jyr_Jform(self): - # self.assertTrue(adjointTest('j', 'jyr')) - # def test_Jtvec_adjointTest_jzr_Jform(self): - # self.assertTrue(adjointTest('j', 'jzr')) - # def test_Jtvec_adjointTest_jxi_Jform(self): - # self.assertTrue(adjointTest('j', 'jxi')) - # def test_Jtvec_adjointTest_jyi_Jform(self): - # self.assertTrue(adjointTest('j', 'jyi')) - # def test_Jtvec_adjointTest_jzi_Jform(self): - # self.assertTrue(adjointTest('j', 'jzi'))R + def test_Jtvec_adjointTest_jxr_Jform(self): + self.assertTrue(adjointTest('j', 'jxr')) + def test_Jtvec_adjointTest_jyr_Jform(self): + self.assertTrue(adjointTest('j', 'jyr')) + def test_Jtvec_adjointTest_jzr_Jform(self): + self.assertTrue(adjointTest('j', 'jzr')) + def test_Jtvec_adjointTest_jxi_Jform(self): + self.assertTrue(adjointTest('j', 'jxi')) + def test_Jtvec_adjointTest_jyi_Jform(self): + self.assertTrue(adjointTest('j', 'jyi')) + def test_Jtvec_adjointTest_jzi_Jform(self): + self.assertTrue(adjointTest('j', 'jzi')) def test_Jtvec_adjointTest_hxr_Jform(self): self.assertTrue(adjointTest('j', 'hxr')) - # def test_Jtvec_adjointTest_hyr_Jform(self): - # self.assertTrue(adjointTest('j', 'hyr')) - # def test_Jtvec_adjointTest_hzr_Jform(self): - # self.assertTrue(adjointTest('j', 'hzr')) - # def test_Jtvec_adjointTest_hxi_Jform(self): - # self.assertTrue(adjointTest('j', 'hxi')) - # def test_Jtvec_adjointTest_hyi_Jform(self): - # self.assertTrue(adjointTest('j', 'hyi')) - # def test_Jtvec_adjointTest_hzi_Jform(self): - # self.assertTrue(adjointTest('j', 'hzi')) + def test_Jtvec_adjointTest_hyr_Jform(self): + self.assertTrue(adjointTest('j', 'hyr')) + def test_Jtvec_adjointTest_hzr_Jform(self): + self.assertTrue(adjointTest('j', 'hzr')) + def test_Jtvec_adjointTest_hxi_Jform(self): + self.assertTrue(adjointTest('j', 'hxi')) + def test_Jtvec_adjointTest_hyi_Jform(self): + self.assertTrue(adjointTest('j', 'hyi')) + def test_Jtvec_adjointTest_hzi_Jform(self): + self.assertTrue(adjointTest('j', 'hzi')) if __name__ == '__main__': From 31c697eafa091858b7a41a6209239615cf267e75 Mon Sep 17 00:00:00 2001 From: Lindsey Date: Wed, 25 Feb 2015 17:06:38 -0800 Subject: [PATCH 04/25] added MeMu, and am using self.mu everywhere instead of mu_0 --- simpegEM/Base.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/simpegEM/Base.py b/simpegEM/Base.py index e0373762..c53447d2 100644 --- a/simpegEM/Base.py +++ b/simpegEM/Base.py @@ -48,11 +48,17 @@ class BaseEMProblem(Problem.BaseProblem): @property def MeMui(self): - #TODO: assuming constant mu if getattr(self, '_MeMui', None) is None: - self._MeMui = self.mesh.getEdgeInnerProduct(1/mu_0) + self._MeMui = self.mesh.getEdgeInnerProduct(1/self.mu) return self._MeMui + @property + def MeMu(self): + #TODO: assuming constant mu + if getattr(self, '_MeMu', None) is None: + self._MeMu = self.mesh.getEdgeInnerProduct(self.mu) + return self._MeMu + @property def Me(self): if getattr(self, '_Me', None) is None: From 3201d1fc59a24e73164f5dd2c34063cbc4d05edc Mon Sep 17 00:00:00 2001 From: Lindsey Date: Wed, 25 Feb 2015 17:58:54 -0800 Subject: [PATCH 05/25] Moved call of store these fields before the __init__ so that we don't mess with the base problem, updated names of mass matrices so it is clear they are mass matrices not phys props --- simpegEM/FDEM/FDEM.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py index c6304203..dee3796f 100644 --- a/simpegEM/FDEM/FDEM.py +++ b/simpegEM/FDEM/FDEM.py @@ -18,7 +18,6 @@ class BaseFDEMProblem(BaseEMProblem): \\nabla \\times \\mu^{-1} \\vec{B} - \\sigma \\vec{E} = \\vec{J_s} """ - surveyPair = SurveyFDEM def forward(self, m, RHS, CalcFields): @@ -338,6 +337,7 @@ class ProblemFDEM_b(BaseFDEMProblem): return None raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) + ########################################################################################## ################################ H-J Formulation ######################################### ########################################################################################## @@ -353,7 +353,7 @@ class ProblemFDEM_j(BaseFDEMProblem): Since \(\\vec{J}\) is a flux and \(\\vec{H}\) is a field, we discretize \(\\vec{J}\) on faces and \(\\vec{H}\) on edges. - For this implementation, we solve for J, using \( \\vec{H} = - (i\\omega\\mu)^{-1} \\nabla \\times \\sigma^{-1} \\vec{J} \) : + For this implementation, we solve for J using \( \\vec{H} = - (i\\omega\\mu)^{-1} \\nabla \\times \\sigma^{-1} \\vec{J} \) : .. math:: \\nabla \\times ( \\mu^{-1} \\nabla \\times \\sigma^{-1} \\vec{J} ) + i\\omega \\vec{J} = - i\\omega\\vec{J_s} @@ -364,10 +364,10 @@ class ProblemFDEM_j(BaseFDEMProblem): """ solType = 'j' + storeTheseFields = ['j','h'] def __init__(self, model, **kwargs): BaseFDEMProblem.__init__(self, model, **kwargs) - BaseFDEMProblem.storeTheseFields = ['j','h'] def getA(self, freq): """ @@ -376,12 +376,12 @@ class ProblemFDEM_j(BaseFDEMProblem): :return: A """ - mui = self.MeMui - sigi = self.MfSigmai + MeMui = self.MeMui + MfSigi = self.MfSigmai C = self.mesh.edgeCurl iomega = 1j * omega(freq) * sp.eye(self.mesh.nF) - return C * mui * C.T * sigi + iomega + return C * MeMui * C.T * MfSigi + iomega def getADeriv(self, freq, u, v, adjoint=False): @@ -424,10 +424,10 @@ class ProblemFDEM_j(BaseFDEMProblem): rhs[i] = np.concatenate((SRCx, SRCy, SRCz)) a = np.concatenate(rhs).reshape((self.mesh.nF, len(Txs)), order='F') - mui = self.MeMui + MeMui = self.MeMui C = self.mesh.edgeCurl - j_s = C*mui*C.T*a + j_s = C*MeMui*C.T*a return -1j*omega(freq)*j_s def calcFields(self, sol, freq, fieldType, adjoint=False): From 5bb5e8f6b3f6c3bac2d44bd9dfd751318f3ef4a1 Mon Sep 17 00:00:00 2001 From: Lindsey Date: Wed, 25 Feb 2015 20:21:44 -0800 Subject: [PATCH 06/25] Start of JH implementation, solving for h. Derivatives for ProblemFDEM_h NOT working yet --- simpegEM/FDEM/FDEM.py | 152 ++++++++++++++++++++++++++++++++++++ simpegEM/FDEM/__init__.py | 2 +- simpegEM/Tests/test_FDEM.py | 60 +++++++++++++- 3 files changed, 211 insertions(+), 3 deletions(-) diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py index dee3796f..e4b94cf2 100644 --- a/simpegEM/FDEM/FDEM.py +++ b/simpegEM/FDEM/FDEM.py @@ -338,6 +338,7 @@ class ProblemFDEM_b(BaseFDEMProblem): raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) + ########################################################################################## ################################ H-J Formulation ######################################### ########################################################################################## @@ -424,6 +425,7 @@ class ProblemFDEM_j(BaseFDEMProblem): rhs[i] = np.concatenate((SRCx, SRCy, SRCz)) a = np.concatenate(rhs).reshape((self.mesh.nF, len(Txs)), order='F') + a = Utils.mkvc(a) MeMui = self.MeMui C = self.mesh.edgeCurl @@ -463,3 +465,153 @@ class ProblemFDEM_j(BaseFDEMProblem): else: return -(1./(1j*omega(freq))) * dsig_dm.T * ( dsigi_dsig.T * ( dMf_dsigi.T * ( C * ( MeMui.T * v ) ) ) ) raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) + + + + +# Solving for h! - NOTE: WE ARE GOING TO NEED dRHS / dm ! +class ProblemFDEM_h(BaseFDEMProblem): + """ + Using the H-J formulation of Maxwell's equations + + .. math:: + \\nabla \\times \\sigma^{-1} \\vec{J} + i\\omega\\mu\\vec{H} = 0 + \\nabla \\times \\vec{H} - \\vec{J} = \\vec{J_s} + + Since \(\\vec{J}\) is a flux and \(\\vec{H}\) is a field, we discretize \(\\vec{J}\) on faces and \(\\vec{H}\) on edges. + + For this implementation, we solve for J using \( \\vec{J} = \\nabla \\times \\vec{H} - \\vec{J_s} \) + + .. math:: + \\nabla \\times \\sigma^{-1} \\nabla \\times \\vec{H} + i\\omega\\mu\\vec{H} = \\nabla \\times \\sigma^{-1} \\vec{J_s} + + .. note:: + This implementation does not yet work with full anisotropy!! + + """ + + solType = 'h' + storeTheseFields = ['j','h'] + + def __init__(self, model, **kwargs): + BaseFDEMProblem.__init__(self, model, **kwargs) + + def getA(self, freq): + """ + :param float freq: Frequency + :rtype: scipy.sparse.csr_matrix + :return: A + """ + + MeMu = self.MeMu + MfSigi = self.MfSigmai + C = self.mesh.edgeCurl + + return C.T * MfSigi * C + 1j*omega(freq)*MeMu + + def getADeriv(self, freq, u, v, adjoint=False): + + MeMu = self.MeMu + C = self.mesh.edgeCurl + sig = self.curModel.transform + sigi = 1/sig + dsig_dm = self.curModel.transformDeriv + dsigi_dsig = -Utils.sdiag(sigi)**2 + + dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(C*u) + + if adjoint: + return (dsig_dm.T * (dsigi_dsig.T * (dMf_dsigi.T * (C * v)))) + return (C.T * (dMf_dsigi * (dsigi_dsig * (dsig_dm * v)))) + + + def getjs(self,freq): + """ + :param float freq: Frequency + :rtype: numpy.ndarray (nE, nTx) + :return: j_s + """ + Txs = self.survey.getTransmitters(freq) + rhs = range(len(Txs)) + for i, tx in enumerate(Txs): + if tx.txType == 'VMD': + src = Sources.MagneticDipoleVectorPotential + SRCx = src(tx.loc, self.mesh.gridFx, 'x') + SRCy = src(tx.loc, self.mesh.gridFy, 'y') + SRCz = src(tx.loc, self.mesh.gridFz, 'z') + + elif tx.txType == 'CircularLoop': + src = Sources.MagneticLoopVectorPotential + SRCx = src(tx.loc, self.mesh.gridFx, 'x', tx.radius) + SRCy = src(tx.loc, self.mesh.gridFy, 'y', tx.radius) + SRCz = src(tx.loc, self.mesh.gridFz, 'z', tx.radius) + else: + raise NotImplemented('%s txType is not implemented' % tx.txType) + rhs[i] = np.concatenate((SRCx, SRCy, SRCz)) + + a = np.concatenate(rhs).reshape((self.mesh.nF, len(Txs)), order='F') + a = Utils.mkvc(a) + + MeMui = self.MeMui + C = self.mesh.edgeCurl + + return C*MeMui*C.T*a + + def getRHS(self, freq): + """ + :param float freq: Frequency + :rtype: numpy.ndarray (nE, nTx) + :return: RHS + """ + MfSigi = self.MfSigmai + C = self.mesh.edgeCurl + j_s = self.getjs(freq) + return C.T*MfSigi*j_s + + def getRHSDeriv(self, freq, v, adjoint=False): + """ + :param float freq: Frequency + :rtype: numpy.ndarray (nE, nTx) + :return: RHSDeriv + """ + C = self.mesh.edgeCurl + sig = self.curModel.transform + sigi = 1/sig + j_s = self.getjs(freq) + dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(j_s) + dsig_dm = self.curModel.transformDeriv + dsigi_dsig = -Utils.sdiag(sigi)**2 # only works for diagonal matrices + + if adjoint: + return dsig_dm.T * dsigi_dsig.T * dMf_dsigi.T * C * v + return C.T * dMf_dsigi * dsigi_dsig * dsig_dm * v + + def calcFields(self, sol, freq, fieldType, adjoint=False): + h = sol + if fieldType == 'j': + NotImplementedError('fieldType "%s" is not implemented.' % fieldType) + C = self.mesh.edgeCurl + j_s = self.getjs(freq) + if adjoint: + NotImplementedError('fieldType "%s" is not implemented.' % fieldType) + # return C.T*h # TODO: THIS IS WRONG + NotImplementedError('fieldType "%s" is not implemented.' % fieldType) + return C*h # - j_s + elif fieldType == 'h': + return h + raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) + + def calcFieldsDeriv(self, sol, freq, fieldType, v, adjoint=False): + h = sol + if fieldType == 'j': + NotImplementedError('fieldType "%s" is not implemented.' % fieldType) + C = self.mesh.edgeCurl + drhs = self.getRHSDeriv(freq,v,adjoint) + if adjoint: + NotImplementedError('fieldType "%s" is not implemented.' % fieldType) + # return -C.T * drhs + NotImplementedError('fieldType "%s" is not implemented.' % fieldType) + # return -C * drhs + elif fieldType == 'h': + return -self.getRHSDeriv(freq,v,adjoint) + raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) \ No newline at end of file diff --git a/simpegEM/FDEM/__init__.py b/simpegEM/FDEM/__init__.py index 4f7fcf4c..562b9218 100644 --- a/simpegEM/FDEM/__init__.py +++ b/simpegEM/FDEM/__init__.py @@ -1,2 +1,2 @@ from SurveyFDEM import * -from FDEM import ProblemFDEM_e, ProblemFDEM_b, ProblemFDEM_j +from FDEM import ProblemFDEM_e, ProblemFDEM_b, ProblemFDEM_j, ProblemFDEM_h diff --git a/simpegEM/Tests/test_FDEM.py b/simpegEM/Tests/test_FDEM.py index 2091a43d..7894dbab 100644 --- a/simpegEM/Tests/test_FDEM.py +++ b/simpegEM/Tests/test_FDEM.py @@ -8,6 +8,7 @@ TOL = 1e-4 FLR = 1e-15 # "zero", so if residual below this --> pass regardless of order CONDUCTIVITY = 1e1 MU = mu_0 +freq = 1 addrandoms = True # important to addrandoms if testing HJ formulation with VMD source! (or else jz ~ 0) def getProblem(fdemType, comp): @@ -24,7 +25,7 @@ def getProblem(fdemType, comp): x = np.linspace(-30,30,6) XYZ = Utils.ndgrid(x,x,np.r_[0]) Rx0 = EM.FDEM.RxFDEM(XYZ, comp) - Tx0 = EM.FDEM.TxFDEM(np.r_[0.,0.,0.], 'VMD', 1, [Rx0]) + Tx0 = EM.FDEM.TxFDEM(np.r_[0.,0.,0.], 'VMD', freq, [Rx0]) survey = EM.FDEM.SurveyFDEM([Tx0]) @@ -36,6 +37,8 @@ def getProblem(fdemType, comp): prb = EM.FDEM.ProblemFDEM_b(mesh, mapping=mapping) elif fdemType == 'j': prb = EM.FDEM.ProblemFDEM_j(mesh, mapping=mapping) + elif fdemType == 'h': + prb = EM.FDEM.ProblemFDEM_h(mesh, mapping=mapping) else: raise NotImplementedError() prb.pair(survey) @@ -83,7 +86,7 @@ def derivTest(fdemType, comp): x0 = x0 + np.random.randn(prb.mesh.nC)*CONDUCTIVITY*1e-1 mu = mu + np.random.randn(prb.mesh.nC)*MU*1e-1 - prb.mu = mu + prb.mu = mu survey = prb.survey def fun(x): return survey.dpred(x), lambda x: prb.Jvec(x0, x) @@ -221,6 +224,32 @@ class FDEM_DerivTests(unittest.TestCase): def test_Jvec_hzi_Jform(self): self.assertTrue(derivTest('j', 'hzi')) + # def test_Jvec_hxr_Hform(self): + # self.assertTrue(derivTest('h', 'hxr')) + # def test_Jvec_hyr_Hform(self): + # self.assertTrue(derivTest('h', 'hyr')) + # def test_Jvec_hzr_Hform(self): + # self.assertTrue(derivTest('h', 'hzr')) + # def test_Jvec_hxi_Hform(self): + # self.assertTrue(derivTest('h', 'hxi')) + # def test_Jvec_hyi_Hform(self): + # self.assertTrue(derivTest('h', 'hyi')) + # def test_Jvec_hzi_Hform(self): + # self.assertTrue(derivTest('h', 'hzi')) + + # def test_Jvec_hxr_Hform(self): + # self.assertTrue(derivTest('h', 'jxr')) + # def test_Jvec_hyr_Hform(self): + # self.assertTrue(derivTest('h', 'jyr')) + # def test_Jvec_hzr_Hform(self): + # self.assertTrue(derivTest('h', 'jzr')) + # def test_Jvec_hxi_Hform(self): + # self.assertTrue(derivTest('h', 'jxi')) + # def test_Jvec_hyi_Hform(self): + # self.assertTrue(derivTest('h', 'jyi')) + # def test_Jvec_hzi_Hform(self): + # self.assertTrue(derivTest('h', 'jzi')) + def test_Jtvec_adjointTest_jxr_Jform(self): self.assertTrue(adjointTest('j', 'jxr')) def test_Jtvec_adjointTest_jyr_Jform(self): @@ -247,6 +276,33 @@ class FDEM_DerivTests(unittest.TestCase): def test_Jtvec_adjointTest_hzi_Jform(self): 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_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__': unittest.main() From 5e323591b93730a231fff53f0e107afa171d3280 Mon Sep 17 00:00:00 2001 From: Lindsey Date: Thu, 26 Feb 2015 10:39:04 -0800 Subject: [PATCH 07/25] expanded documentation for solving for J --- simpegEM/FDEM/FDEM.py | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py index e4b94cf2..5343c45d 100644 --- a/simpegEM/FDEM/FDEM.py +++ b/simpegEM/FDEM/FDEM.py @@ -359,6 +359,11 @@ class ProblemFDEM_j(BaseFDEMProblem): .. math:: \\nabla \\times ( \\mu^{-1} \\nabla \\times \\sigma^{-1} \\vec{J} ) + i\\omega \\vec{J} = - i\\omega\\vec{J_s} + We discretize this to: + + .. math:: + (\\mathbf{C} \\mathbf{M^e_{mu^{-1}}} \\mathbf{C^T} \\mathbf{M^f_{\\sigma^{-1}}} + i\\omega ) \\mathbf{j} = - i\\omega \\mathbf{j_s} + .. note:: This implementation does not yet work with full anisotropy!! @@ -372,6 +377,10 @@ class ProblemFDEM_j(BaseFDEMProblem): def getA(self, freq): """ + Here, we form the operator \(\\mathbf{A}\) to solce + .. math:: + \\mathbf{A} = \\mathbf{C} \\mathbf{M^e_{mu^{-1}}} \\mathbf{C^T} \\mathbf{M^f_{\\sigma^{-1}}} + i\\omega + :param float freq: Frequency :rtype: scipy.sparse.csr_matrix :return: A @@ -384,7 +393,14 @@ class ProblemFDEM_j(BaseFDEMProblem): return C * MeMui * C.T * MfSigi + iomega + def getADeriv(self, freq, u, v, adjoint=False): + """ + In this case, we assume that electrical conductivity, \(\\sigma\) is the physical property of interest (i.e. \(\sigma\) = model.transform). Then we want + .. math:: + \\frac{\mathbf{A(\\sigma)} \mathbf{v}}{d \\mathbf{m}} &= \\mathbf{C} \\mathbf{M^e_{mu^{-1}}} \\mathbf{C^T} \\frac{d \\mathbf{M^f_{\\sigma^{-1}}}}{d \\mathbf{m}} + &= \\mathbf{C} \\mathbf{M^e_{mu^{-1}}} \\mathbf{C^T} \\frac{d \\mathbf{M^f_{\\sigma^{-1}}}}{d \\mathbf{\\sigma^{-1}}} \\frac{d \\mathbf{\\sigma^{-1}}}{d \\mathbf{\\sigma}} \\frac{d \\mathbf{\\sigma}}{d \\mathbf{m}} + """ MeMui = self.MeMui C = self.mesh.edgeCurl From 8d4e001301ca7a97db6d779a518c9e4ccc87f14d Mon Sep 17 00:00:00 2001 From: Lindsey Date: Thu, 26 Feb 2015 13:51:42 -0800 Subject: [PATCH 08/25] HJ formulation solving for h, with h data --- simpegEM/FDEM/FDEM.py | 66 +++++++++++++++++++++++++++---------- simpegEM/Tests/test_FDEM.py | 48 +++++++++++++-------------- 2 files changed, 73 insertions(+), 41 deletions(-) diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py index 5343c45d..56edd731 100644 --- a/simpegEM/FDEM/FDEM.py +++ b/simpegEM/FDEM/FDEM.py @@ -538,7 +538,7 @@ class ProblemFDEM_h(BaseFDEMProblem): if adjoint: return (dsig_dm.T * (dsigi_dsig.T * (dMf_dsigi.T * (C * v)))) - return (C.T * (dMf_dsigi * (dsigi_dsig * (dsig_dm * v)))) + return (C.T * (dMf_dsigi * (dsigi_dsig * (dsig_dm * v)))) def getjs(self,freq): @@ -606,28 +606,60 @@ class ProblemFDEM_h(BaseFDEMProblem): h = sol if fieldType == 'j': NotImplementedError('fieldType "%s" is not implemented.' % fieldType) - C = self.mesh.edgeCurl - j_s = self.getjs(freq) - if adjoint: - NotImplementedError('fieldType "%s" is not implemented.' % fieldType) - # return C.T*h # TODO: THIS IS WRONG - NotImplementedError('fieldType "%s" is not implemented.' % fieldType) - return C*h # - j_s + # C = self.mesh.edgeCurl + # j_s = self.getjs(freq) + # if adjoint: + # MeMu = self.MeMu + # MfSigi = self.MfSigmai + # MfSigmaiinv = self.Solver(MfSigi.T, **self.solverOpts) + # Cinv = self.Solver(C, **self.solverOpts) + # return -1j * omega(freq) * MeMu.T * (MfSigmaiinv * (CTinv * h)) + # return C*h - j_s # -iomega(freq) inv(MfSigmai) inv(C.T) MeMu elif fieldType == 'h': return h raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) def calcFieldsDeriv(self, sol, freq, fieldType, v, adjoint=False): h = sol + A = self.getA(freq) + if fieldType == 'j': - NotImplementedError('fieldType "%s" is not implemented.' % fieldType) - C = self.mesh.edgeCurl - drhs = self.getRHSDeriv(freq,v,adjoint) - if adjoint: - NotImplementedError('fieldType "%s" is not implemented.' % fieldType) - # return -C.T * drhs - NotImplementedError('fieldType "%s" is not implemented.' % fieldType) - # return -C * drhs + # C = self.mesh.edgeCurl + # MeMu = self.MeMu + # MfSigi = self.MfSigmai + # MfSigmaiinv = self.Solver(MfSigi, **self.solverOpts) + # CTinv = self.Solver(C.T, **self.solverOpts) + # hDeriv = self.calcFieldsDeriv(h,freq,'h',v,adjoint) + + # pt1 = -1j * omega(freq) * (MfSigmaiinv * (CTinv * (MeMu * hDeriv))) + + # sig = self.curModel.transform + # sigi = 1/sig + # dsig_dm = self.curModel.transformDeriv + # dsigi_dsig = -Utils.sdiag(sigi)**2 + + # v1 = MfSigmaiinv *(CTinv * (MeMu * h)) + + # dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(v1) + + # pt2 = 1j * omega(freq) * (MfSigmaiinv * (dMf_dsigi * v1)) + + # return pt1+pt2 + + + # if adjoint: + # MfSigmaiTinv = self.Solver(MfSigi.T, **self.solverOpts) + # Cinv = self.Solver(C, **self.solverOpts) + # MeMu.T * (Cinv * (MfSigmaiTinv * h)) + # pt1 = -1j * omega(freq) * hDeriv.T * MeMu.T * Cinv *(MfSigmaiTinv * v) + raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) + elif fieldType == 'h': - return -self.getRHSDeriv(freq,v,adjoint) + if adjoint: + ATinv = self.Solver(A.T, **self.solverOpts) + ATinvv = ATinv*v + return self.getRHSDeriv(freq,ATinvv,adjoint=True) + dRHSh = self.getRHSDeriv(freq,v,adjoint) + Ainv = self.Solver(A, **self.solverOpts) + return Ainv*dRHSh raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) \ No newline at end of file diff --git a/simpegEM/Tests/test_FDEM.py b/simpegEM/Tests/test_FDEM.py index 7894dbab..453cd756 100644 --- a/simpegEM/Tests/test_FDEM.py +++ b/simpegEM/Tests/test_FDEM.py @@ -224,18 +224,18 @@ class FDEM_DerivTests(unittest.TestCase): def test_Jvec_hzi_Jform(self): self.assertTrue(derivTest('j', 'hzi')) - # def test_Jvec_hxr_Hform(self): - # self.assertTrue(derivTest('h', 'hxr')) - # def test_Jvec_hyr_Hform(self): - # self.assertTrue(derivTest('h', 'hyr')) - # def test_Jvec_hzr_Hform(self): - # self.assertTrue(derivTest('h', 'hzr')) - # def test_Jvec_hxi_Hform(self): - # self.assertTrue(derivTest('h', 'hxi')) - # def test_Jvec_hyi_Hform(self): - # self.assertTrue(derivTest('h', 'hyi')) - # def test_Jvec_hzi_Hform(self): - # self.assertTrue(derivTest('h', 'hzi')) + def test_Jvec_hxr_Hform(self): + self.assertTrue(derivTest('h', 'hxr')) + def test_Jvec_hyr_Hform(self): + self.assertTrue(derivTest('h', 'hyr')) + def test_Jvec_hzr_Hform(self): + self.assertTrue(derivTest('h', 'hzr')) + def test_Jvec_hxi_Hform(self): + self.assertTrue(derivTest('h', 'hxi')) + def test_Jvec_hyi_Hform(self): + self.assertTrue(derivTest('h', 'hyi')) + def test_Jvec_hzi_Hform(self): + self.assertTrue(derivTest('h', 'hzi')) # def test_Jvec_hxr_Hform(self): # self.assertTrue(derivTest('h', 'jxr')) @@ -276,18 +276,18 @@ class FDEM_DerivTests(unittest.TestCase): def test_Jtvec_adjointTest_hzi_Jform(self): 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_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')) From 5d74fc2d1bd5d0a592469ea710a53ff2ba78749b Mon Sep 17 00:00:00 2001 From: Lindsey Date: Thu, 26 Feb 2015 14:27:49 -0800 Subject: [PATCH 09/25] Changed MeMui to MeMuI, it should be the full inverse --- simpegEM/Base.py | 3 ++- simpegEM/FDEM/FDEM.py | 42 +++++++++++++++++++++++++----------------- 2 files changed, 27 insertions(+), 18 deletions(-) diff --git a/simpegEM/Base.py b/simpegEM/Base.py index c53447d2..1f2b9e45 100644 --- a/simpegEM/Base.py +++ b/simpegEM/Base.py @@ -47,7 +47,8 @@ class BaseEMProblem(Problem.BaseProblem): return self._MfMui @property - def MeMui(self): + def MeMuI(self): + # Assuming isotropic mu! if getattr(self, '_MeMui', None) is None: self._MeMui = self.mesh.getEdgeInnerProduct(1/self.mu) return self._MeMui diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py index 56edd731..275b82d2 100644 --- a/simpegEM/FDEM/FDEM.py +++ b/simpegEM/FDEM/FDEM.py @@ -386,12 +386,12 @@ class ProblemFDEM_j(BaseFDEMProblem): :return: A """ - MeMui = self.MeMui + MeMuI = self.MeMuI MfSigi = self.MfSigmai C = self.mesh.edgeCurl iomega = 1j * omega(freq) * sp.eye(self.mesh.nF) - return C * MeMui * C.T * MfSigi + iomega + return C * MeMuI * C.T * MfSigi + iomega def getADeriv(self, freq, u, v, adjoint=False): @@ -399,10 +399,10 @@ class ProblemFDEM_j(BaseFDEMProblem): In this case, we assume that electrical conductivity, \(\\sigma\) is the physical property of interest (i.e. \(\sigma\) = model.transform). Then we want .. math:: \\frac{\mathbf{A(\\sigma)} \mathbf{v}}{d \\mathbf{m}} &= \\mathbf{C} \\mathbf{M^e_{mu^{-1}}} \\mathbf{C^T} \\frac{d \\mathbf{M^f_{\\sigma^{-1}}}}{d \\mathbf{m}} - &= \\mathbf{C} \\mathbf{M^e_{mu^{-1}}} \\mathbf{C^T} \\frac{d \\mathbf{M^f_{\\sigma^{-1}}}}{d \\mathbf{\\sigma^{-1}}} \\frac{d \\mathbf{\\sigma^{-1}}}{d \\mathbf{\\sigma}} \\frac{d \\mathbf{\\sigma}}{d \\mathbf{m}} + &= \\mathbf{C} \\mathbf{M^e_{mu}^{-1}} \\mathbf{C^T} \\frac{d \\mathbf{M^f_{\\sigma^{-1}}}}{d \\mathbf{\\sigma^{-1}}} \\frac{d \\mathbf{\\sigma^{-1}}}{d \\mathbf{\\sigma}} \\frac{d \\mathbf{\\sigma}}{d \\mathbf{m}} """ - MeMui = self.MeMui + MeMuI = self.MeMuI C = self.mesh.edgeCurl sig = self.curModel.transform sigi = 1/sig @@ -411,9 +411,9 @@ class ProblemFDEM_j(BaseFDEMProblem): dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(u) if adjoint: - return dsig_dm.T * ( dsigi_dsig.T *( dMf_dsigi.T * ( C * ( MeMui.T * ( C.T * v ) ) ) ) ) + return dsig_dm.T * ( dsigi_dsig.T *( dMf_dsigi.T * ( C * ( MeMuI.T * ( C.T * v ) ) ) ) ) - return C * ( MeMui * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) ) ) + return C * ( MeMuI * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) ) ) def getRHS(self, freq): @@ -442,10 +442,10 @@ class ProblemFDEM_j(BaseFDEMProblem): a = np.concatenate(rhs).reshape((self.mesh.nF, len(Txs)), order='F') a = Utils.mkvc(a) - MeMui = self.MeMui + MeMuI = self.MeMuI C = self.mesh.edgeCurl - j_s = C*MeMui*C.T*a + j_s = C*MeMuI*C.T*a return -1j*omega(freq)*j_s def calcFields(self, sol, freq, fieldType, adjoint=False): @@ -453,7 +453,7 @@ class ProblemFDEM_j(BaseFDEMProblem): if fieldType == 'j': return j elif fieldType == 'h': - mui = self.MeMui + mui = self.MeMuI C = self.mesh.edgeCurl MfSigi = self.MfSigmai if not adjoint: @@ -468,7 +468,7 @@ class ProblemFDEM_j(BaseFDEMProblem): if fieldType == 'j': return None elif fieldType == 'h': - MeMui = self.MeMui + MeMuI = self.MeMuI C = self.mesh.edgeCurl sig = self.curModel.transform sigi = 1/sig @@ -477,9 +477,9 @@ class ProblemFDEM_j(BaseFDEMProblem): dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(j) sigi = self.MfSigmai if not adjoint: - return -(1./(1j*omega(freq))) * MeMui * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) ) + return -(1./(1j*omega(freq))) * MeMuI * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) ) else: - return -(1./(1j*omega(freq))) * dsig_dm.T * ( dsigi_dsig.T * ( dMf_dsigi.T * ( C * ( MeMui.T * v ) ) ) ) + return -(1./(1j*omega(freq))) * dsig_dm.T * ( dsigi_dsig.T * ( dMf_dsigi.T * ( C * ( MeMuI.T * v ) ) ) ) raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) @@ -501,6 +501,11 @@ class ProblemFDEM_h(BaseFDEMProblem): .. math:: \\nabla \\times \\sigma^{-1} \\nabla \\times \\vec{H} + i\\omega\\mu\\vec{H} = \\nabla \\times \\sigma^{-1} \\vec{J_s} + We discretize and solve + + .. math:: + (\\mathbf{C^T} \\mathbf{M^f_{\\sigma^{-1}}} \\mathbf{C} + i\\omega \\mathbf{M_{\mu}} ) \\mathbf{h} = \\mathbf{C^T} \\mathbf{M^f_{\\sigma^{-1}}} \\vec{J_s} + .. note:: This implementation does not yet work with full anisotropy!! @@ -568,10 +573,10 @@ class ProblemFDEM_h(BaseFDEMProblem): a = np.concatenate(rhs).reshape((self.mesh.nF, len(Txs)), order='F') a = Utils.mkvc(a) - MeMui = self.MeMui + MeMuI = self.MeMuI C = self.mesh.edgeCurl - return C*MeMui*C.T*a + return C*MeMuI*C.T*a def getRHS(self, freq): """ @@ -629,7 +634,7 @@ class ProblemFDEM_h(BaseFDEMProblem): # MfSigi = self.MfSigmai # MfSigmaiinv = self.Solver(MfSigi, **self.solverOpts) # CTinv = self.Solver(C.T, **self.solverOpts) - # hDeriv = self.calcFieldsDeriv(h,freq,'h',v,adjoint) + # hDeriv = self.calcFieldsDeriv(h,freq,'h',v,adjoint=False) # pt1 = -1j * omega(freq) * (MfSigmaiinv * (CTinv * (MeMu * hDeriv))) @@ -642,7 +647,7 @@ class ProblemFDEM_h(BaseFDEMProblem): # dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(v1) - # pt2 = 1j * omega(freq) * (MfSigmaiinv * (dMf_dsigi * v1)) + # pt2 = 1j * omega(freq) * (MfSigmaiinv * (dMf_dsigi * (dsigi_dsig * (dsig_dm * v)))) # return pt1+pt2 @@ -651,7 +656,10 @@ class ProblemFDEM_h(BaseFDEMProblem): # MfSigmaiTinv = self.Solver(MfSigi.T, **self.solverOpts) # Cinv = self.Solver(C, **self.solverOpts) # MeMu.T * (Cinv * (MfSigmaiTinv * h)) - # pt1 = -1j * omega(freq) * hDeriv.T * MeMu.T * Cinv *(MfSigmaiTinv * v) + # v1 = MeMu.T *( Cinv *(MfSigmaiTinv * v)) + # pt1 = -1j * omega(freq) * self.calcFieldsDeriv(h,freq,'h',v,adjoint=True) + + # pt2 = raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) elif fieldType == 'h': From 5284b91f33a745642e9fa7327afc8ef02c7d46b1 Mon Sep 17 00:00:00 2001 From: Lindsey Date: Thu, 26 Feb 2015 14:57:08 -0800 Subject: [PATCH 10/25] changed _MeMui to _MeMuI --- simpegEM/Base.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/simpegEM/Base.py b/simpegEM/Base.py index 1f2b9e45..6d7550d3 100644 --- a/simpegEM/Base.py +++ b/simpegEM/Base.py @@ -49,9 +49,9 @@ class BaseEMProblem(Problem.BaseProblem): @property def MeMuI(self): # Assuming isotropic mu! - if getattr(self, '_MeMui', None) is None: + if getattr(self, '_MeMuI', None) is None: self._MeMui = self.mesh.getEdgeInnerProduct(1/self.mu) - return self._MeMui + return self._MeMuI @property def MeMu(self): From 5ac746f31fde277d837e64159395f3395a722c29 Mon Sep 17 00:00:00 2001 From: Lindsey Date: Thu, 26 Feb 2015 15:03:04 -0800 Subject: [PATCH 11/25] added MeMu and changed MeMui to MeMuI in mu setter --- simpegEM/Base.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/simpegEM/Base.py b/simpegEM/Base.py index 6d7550d3..419a25ff 100644 --- a/simpegEM/Base.py +++ b/simpegEM/Base.py @@ -30,8 +30,10 @@ class BaseEMProblem(Problem.BaseProblem): def mu(self, value): if getattr(self, '_MfMui', None) is not None: del self._MfMui - if getattr(self, '_MeMui', None) is not None: - del self._MeMui + if getattr(self, '_MeMu', None) is not None: + del delf._MeMu + if getattr(self, '_MeMuI', None) is not None: + del self._MeMuI self._mu = value From d5eef78b57825d6a8013bab3a229cc71c599a6c8 Mon Sep 17 00:00:00 2001 From: Lindsey Date: Thu, 26 Feb 2015 15:04:14 -0800 Subject: [PATCH 12/25] fixed typo on MeMuI property --- simpegEM/Base.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/simpegEM/Base.py b/simpegEM/Base.py index 419a25ff..38fca561 100644 --- a/simpegEM/Base.py +++ b/simpegEM/Base.py @@ -52,7 +52,7 @@ class BaseEMProblem(Problem.BaseProblem): def MeMuI(self): # Assuming isotropic mu! if getattr(self, '_MeMuI', None) is None: - self._MeMui = self.mesh.getEdgeInnerProduct(1/self.mu) + self._MeMuI = self.mesh.getEdgeInnerProduct(1/self.mu) return self._MeMuI @property From 0e0033eb0aa39f5ff12f5d029d8db287f3edb4b5 Mon Sep 17 00:00:00 2001 From: Lindsey Date: Thu, 26 Feb 2015 15:40:03 -0800 Subject: [PATCH 13/25] HJ formulation, Problem_h solving for j implemented. I don't yet trust the solution --- simpegEM/FDEM/FDEM.py | 68 +++++++++++++++++++++++-------------------- 1 file changed, 36 insertions(+), 32 deletions(-) diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py index 275b82d2..4638936b 100644 --- a/simpegEM/FDEM/FDEM.py +++ b/simpegEM/FDEM/FDEM.py @@ -610,16 +610,16 @@ class ProblemFDEM_h(BaseFDEMProblem): def calcFields(self, sol, freq, fieldType, adjoint=False): h = sol if fieldType == 'j': - NotImplementedError('fieldType "%s" is not implemented.' % fieldType) - # C = self.mesh.edgeCurl - # j_s = self.getjs(freq) - # if adjoint: - # MeMu = self.MeMu - # MfSigi = self.MfSigmai - # MfSigmaiinv = self.Solver(MfSigi.T, **self.solverOpts) - # Cinv = self.Solver(C, **self.solverOpts) - # return -1j * omega(freq) * MeMu.T * (MfSigmaiinv * (CTinv * h)) - # return C*h - j_s # -iomega(freq) inv(MfSigmai) inv(C.T) MeMu + # NotImplementedError('fieldType "%s" is not implemented.' % fieldType) + C = self.mesh.edgeCurl + j_s = self.getjs(freq) + if adjoint: + # MeMuI = self.MeMuI + # MfSigi = self.MfSigmai + + return C.T*h + # return -1j * omega(freq) * MeMu.T * (MfSigmaiinv * (CTinv * h)) + return C*h #- j_s # -iomega(freq) inv(MfSigmai) inv(C.T) MeMu elif fieldType == 'h': return h raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) @@ -629,37 +629,41 @@ class ProblemFDEM_h(BaseFDEMProblem): A = self.getA(freq) if fieldType == 'j': - # C = self.mesh.edgeCurl - # MeMu = self.MeMu - # MfSigi = self.MfSigmai - # MfSigmaiinv = self.Solver(MfSigi, **self.solverOpts) - # CTinv = self.Solver(C.T, **self.solverOpts) - # hDeriv = self.calcFieldsDeriv(h,freq,'h',v,adjoint=False) + C = self.mesh.edgeCurl + # MeMu = self.MeMu + # MfSigi = self.MfSigmai - # pt1 = -1j * omega(freq) * (MfSigmaiinv * (CTinv * (MeMu * hDeriv))) + # if adjoint: + # MfSigiTCinv = self.Solver(MfSigi.T*C, **self.solverOpts) + # MeMu.T * (Cinv * (MfSigmaiTinv * h)) + # v1 = MeMu.T *( Cinv *(MfSigmaiTinv * v)) + # pt1 = -1j * omega(freq) * self.calcFieldsDeriv(h,freq,'h',v,adjoint=True) - # sig = self.curModel.transform - # sigi = 1/sig - # dsig_dm = self.curModel.transformDeriv - # dsigi_dsig = -Utils.sdiag(sigi)**2 + # pt2 = 1j * omega(freq) * dsig_dm.T * ( dsigi_dsig.T * ( dMf_dsigi.T * ( MfSigiTCinv * v) ) ) + # return pt1 + pt2 - # v1 = MfSigmaiinv *(CTinv * (MeMu * h)) + # CTMfSigiinv = self.Solver(C.T*MfSigi, **self.solverOpts) + # hDeriv = self.calcFieldsDeriv(h,freq,'h',v,adjoint=False) - # dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(v1) + # pt1 = -1j * omega(freq) * (CTMfSigiinv * (MeMu * hDeriv)) - # pt2 = 1j * omega(freq) * (MfSigmaiinv * (dMf_dsigi * (dsigi_dsig * (dsig_dm * v)))) + # sig = self.curModel.transform + # sigi = 1/sig + # dsig_dm = self.curModel.transformDeriv + # dsigi_dsig = -Utils.sdiag(sigi)**2 - # return pt1+pt2 + # v1 = CTMfSigiinv * (MeMu * h) + # dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(v1) - # if adjoint: - # MfSigmaiTinv = self.Solver(MfSigi.T, **self.solverOpts) - # Cinv = self.Solver(C, **self.solverOpts) - # MeMu.T * (Cinv * (MfSigmaiTinv * h)) - # v1 = MeMu.T *( Cinv *(MfSigmaiTinv * v)) - # pt1 = -1j * omega(freq) * self.calcFieldsDeriv(h,freq,'h',v,adjoint=True) + # pt2 = 1j * omega(freq) * (CTMfSigiinv * (dMf_dsigi * (dsigi_dsig * (dsig_dm * v)))) - # pt2 = + # return pt1+pt2 + if adjoint: + dh = self.calcFieldsDeriv(h,freq,'h',C.T*v,adjoint=True) + return dh + dh = self.calcFieldsDeriv(h,freq,'h',v) + return C*dh raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) elif fieldType == 'h': From 94d12b257af8d8da46f109193171a5613dc9d9bc Mon Sep 17 00:00:00 2001 From: Lindsey Date: Thu, 26 Feb 2015 15:40:51 -0800 Subject: [PATCH 14/25] Test everything --- simpegEM/Tests/test_FDEM.py | 48 ++++++++++++++++++------------------- 1 file changed, 24 insertions(+), 24 deletions(-) diff --git a/simpegEM/Tests/test_FDEM.py b/simpegEM/Tests/test_FDEM.py index 453cd756..eb0dccda 100644 --- a/simpegEM/Tests/test_FDEM.py +++ b/simpegEM/Tests/test_FDEM.py @@ -237,18 +237,18 @@ class FDEM_DerivTests(unittest.TestCase): def test_Jvec_hzi_Hform(self): self.assertTrue(derivTest('h', 'hzi')) - # def test_Jvec_hxr_Hform(self): - # self.assertTrue(derivTest('h', 'jxr')) - # def test_Jvec_hyr_Hform(self): - # self.assertTrue(derivTest('h', 'jyr')) - # def test_Jvec_hzr_Hform(self): - # self.assertTrue(derivTest('h', 'jzr')) - # def test_Jvec_hxi_Hform(self): - # self.assertTrue(derivTest('h', 'jxi')) - # def test_Jvec_hyi_Hform(self): - # self.assertTrue(derivTest('h', 'jyi')) - # def test_Jvec_hzi_Hform(self): - # self.assertTrue(derivTest('h', 'jzi')) + def test_Jvec_hxr_Hform(self): + self.assertTrue(derivTest('h', 'jxr')) + def test_Jvec_hyr_Hform(self): + self.assertTrue(derivTest('h', 'jyr')) + def test_Jvec_hzr_Hform(self): + self.assertTrue(derivTest('h', 'jzr')) + def test_Jvec_hxi_Hform(self): + self.assertTrue(derivTest('h', 'jxi')) + def test_Jvec_hyi_Hform(self): + self.assertTrue(derivTest('h', 'jyi')) + def test_Jvec_hzi_Hform(self): + self.assertTrue(derivTest('h', 'jzi')) def test_Jtvec_adjointTest_jxr_Jform(self): self.assertTrue(adjointTest('j', 'jxr')) @@ -289,18 +289,18 @@ class FDEM_DerivTests(unittest.TestCase): 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_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')) From 564fa14826353d7ca2f575835d7356c60eb407ac Mon Sep 17 00:00:00 2001 From: Lindsey Date: Thu, 26 Feb 2015 16:30:11 -0800 Subject: [PATCH 15/25] Added CrossCheck test to make sure that both formulations give the same results. HJ is going to make travis fail... --- simpegEM/Tests/test_FDEM.py | 495 ++++++++++++++++++++++-------------- 1 file changed, 298 insertions(+), 197 deletions(-) diff --git a/simpegEM/Tests/test_FDEM.py b/simpegEM/Tests/test_FDEM.py index eb0dccda..ffbed6a6 100644 --- a/simpegEM/Tests/test_FDEM.py +++ b/simpegEM/Tests/test_FDEM.py @@ -43,11 +43,11 @@ def getProblem(fdemType, comp): raise NotImplementedError() prb.pair(survey) - # try: - # from pymatsolver import MumpsSolver - # prb.Solver = MumpsSolver - # except ImportError, e: - # pass + try: + from pymatsolver import MumpsSolver + prb.Solver = MumpsSolver + except ImportError, e: + pass return prb @@ -93,216 +93,317 @@ def derivTest(fdemType, comp): return Tests.checkDerivative(fun, x0, num=3, plotIt=False, eps=FLR) +def crossCheckTest(fdemType, comp): + + l2norm = lambda r: np.sqrt(r.dot(r)) + + prb1 = getProblem(fdemType, comp) + mesh = prb1.mesh + print '%s formulation - %s' % (fdemType, 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)*CONDUCTIVITY*1e-1 + mu = mu + np.random.randn(mesh.nC)*MU*1e-1 + + prb1.mu = mu + survey1 = prb1.survey + + u1 = prb1.fields(m) + d1 = Utils.mkvc(survey1.projectFields(u1)) + + prb1.unpair + + if fdemType == 'e': + prb2 = getProblem('b', comp) + elif fdemType == 'b': + prb2 = getProblem('e', comp) + elif fdemType == 'j': + prb2 = getProblem('h', comp) + elif fdemType == 'h': + prb2 = getProblem('j', comp) + else: + raise NotImplementedError() + + prb2.mu = mu + survey2 = prb2.survey + + u2 = prb2.fields(m) + d2 = Utils.mkvc(survey2.projectFields(u2)) + + 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_DerivTests(unittest.TestCase): - def test_Jvec_exr_Eform(self): - self.assertTrue(derivTest('e', 'exr')) - def test_Jvec_exr_Bform(self): - self.assertTrue(derivTest('b', 'exr')) - def test_Jvec_eyr_Eform(self): - self.assertTrue(derivTest('e', 'eyr')) - def test_Jvec_eyr_Bform(self): - self.assertTrue(derivTest('b', 'eyr')) - def test_Jvec_ezr_Eform(self): - self.assertTrue(derivTest('e', 'ezr')) - def test_Jvec_ezr_Bform(self): - self.assertTrue(derivTest('b', 'ezr')) - def test_Jvec_exi_Eform(self): - self.assertTrue(derivTest('e', 'exi')) - def test_Jvec_exi_Bform(self): - self.assertTrue(derivTest('b', 'exi')) - def test_Jvec_eyi_Eform(self): - self.assertTrue(derivTest('e', 'eyi')) - def test_Jvec_eyi_Bform(self): - self.assertTrue(derivTest('b', 'eyi')) - def test_Jvec_ezi_Eform(self): - self.assertTrue(derivTest('e', 'ezi')) - def test_Jvec_ezi_Bform(self): - self.assertTrue(derivTest('b', 'ezi')) + # def test_Jvec_exr_Eform(self): + # self.assertTrue(derivTest('e', 'exr')) + # def test_Jvec_exr_Bform(self): + # self.assertTrue(derivTest('b', 'exr')) + # def test_Jvec_eyr_Eform(self): + # self.assertTrue(derivTest('e', 'eyr')) + # def test_Jvec_eyr_Bform(self): + # self.assertTrue(derivTest('b', 'eyr')) + # def test_Jvec_ezr_Eform(self): + # self.assertTrue(derivTest('e', 'ezr')) + # def test_Jvec_ezr_Bform(self): + # self.assertTrue(derivTest('b', 'ezr')) + # def test_Jvec_exi_Eform(self): + # self.assertTrue(derivTest('e', 'exi')) + # def test_Jvec_exi_Bform(self): + # self.assertTrue(derivTest('b', 'exi')) + # def test_Jvec_eyi_Eform(self): + # self.assertTrue(derivTest('e', 'eyi')) + # def test_Jvec_eyi_Bform(self): + # self.assertTrue(derivTest('b', 'eyi')) + # def test_Jvec_ezi_Eform(self): + # self.assertTrue(derivTest('e', 'ezi')) + # def test_Jvec_ezi_Bform(self): + # self.assertTrue(derivTest('b', 'ezi')) - def test_Jvec_bxr_Eform(self): - self.assertTrue(derivTest('e', 'bxr')) - def test_Jvec_bxr_Bform(self): - self.assertTrue(derivTest('b', 'bxr')) - def test_Jvec_byr_Eform(self): - self.assertTrue(derivTest('e', 'byr')) - def test_Jvec_byr_Bform(self): - self.assertTrue(derivTest('b', 'byr')) - def test_Jvec_bzr_Eform(self): - self.assertTrue(derivTest('e', 'bzr')) - def test_Jvec_bzr_Bform(self): - self.assertTrue(derivTest('b', 'bzr')) - def test_Jvec_bxi_Eform(self): - self.assertTrue(derivTest('e', 'bxi')) - def test_Jvec_bxi_Bform(self): - self.assertTrue(derivTest('b', 'bxi')) - def test_Jvec_byi_Eform(self): - self.assertTrue(derivTest('e', 'byi')) - def test_Jvec_byi_Bform(self): - self.assertTrue(derivTest('b', 'byi')) - def test_Jvec_bzi_Eform(self): - self.assertTrue(derivTest('e', 'bzi')) - def test_Jvec_bzi_Bform(self): - self.assertTrue(derivTest('b', 'bzi')) + # def test_Jvec_bxr_Eform(self): + # self.assertTrue(derivTest('e', 'bxr')) + # def test_Jvec_bxr_Bform(self): + # self.assertTrue(derivTest('b', 'bxr')) + # def test_Jvec_byr_Eform(self): + # self.assertTrue(derivTest('e', 'byr')) + # def test_Jvec_byr_Bform(self): + # self.assertTrue(derivTest('b', 'byr')) + # def test_Jvec_bzr_Eform(self): + # self.assertTrue(derivTest('e', 'bzr')) + # def test_Jvec_bzr_Bform(self): + # self.assertTrue(derivTest('b', 'bzr')) + # def test_Jvec_bxi_Eform(self): + # self.assertTrue(derivTest('e', 'bxi')) + # def test_Jvec_bxi_Bform(self): + # self.assertTrue(derivTest('b', 'bxi')) + # def test_Jvec_byi_Eform(self): + # self.assertTrue(derivTest('e', 'byi')) + # def test_Jvec_byi_Bform(self): + # self.assertTrue(derivTest('b', 'byi')) + # def test_Jvec_bzi_Eform(self): + # self.assertTrue(derivTest('e', 'bzi')) + # def test_Jvec_bzi_Bform(self): + # self.assertTrue(derivTest('b', 'bzi')) - def test_Jtvec_adjointTest_exr_Eform(self): - self.assertTrue(adjointTest('e', 'exr')) - def test_Jtvec_adjointTest_exr_Bform(self): - self.assertTrue(adjointTest('b', 'exr')) - def test_Jtvec_adjointTest_eyr_Eform(self): - self.assertTrue(adjointTest('e', 'eyr')) - def test_Jtvec_adjointTest_eyr_Bform(self): - self.assertTrue(adjointTest('b', 'eyr')) - def test_Jtvec_adjointTest_ezr_Eform(self): - self.assertTrue(adjointTest('e', 'ezr')) - def test_Jtvec_adjointTest_ezr_Bform(self): - self.assertTrue(adjointTest('b', 'ezr')) - def test_Jtvec_adjointTest_exi_Eform(self): - self.assertTrue(adjointTest('e', 'exi')) - def test_Jtvec_adjointTest_exi_Bform(self): - self.assertTrue(adjointTest('b', 'exi')) - def test_Jtvec_adjointTest_eyi_Eform(self): - self.assertTrue(adjointTest('e', 'eyi')) - def test_Jtvec_adjointTest_eyi_Bform(self): - self.assertTrue(adjointTest('b', 'eyi')) - def test_Jtvec_adjointTest_ezi_Eform(self): - self.assertTrue(adjointTest('e', 'ezi')) - def test_Jtvec_adjointTest_ezi_Bform(self): - self.assertTrue(adjointTest('b', 'ezi')) + # def test_Jtvec_adjointTest_exr_Eform(self): + # self.assertTrue(adjointTest('e', 'exr')) + # def test_Jtvec_adjointTest_exr_Bform(self): + # self.assertTrue(adjointTest('b', 'exr')) + # def test_Jtvec_adjointTest_eyr_Eform(self): + # self.assertTrue(adjointTest('e', 'eyr')) + # def test_Jtvec_adjointTest_eyr_Bform(self): + # self.assertTrue(adjointTest('b', 'eyr')) + # def test_Jtvec_adjointTest_ezr_Eform(self): + # self.assertTrue(adjointTest('e', 'ezr')) + # def test_Jtvec_adjointTest_ezr_Bform(self): + # self.assertTrue(adjointTest('b', 'ezr')) + # def test_Jtvec_adjointTest_exi_Eform(self): + # self.assertTrue(adjointTest('e', 'exi')) + # def test_Jtvec_adjointTest_exi_Bform(self): + # self.assertTrue(adjointTest('b', 'exi')) + # def test_Jtvec_adjointTest_eyi_Eform(self): + # self.assertTrue(adjointTest('e', 'eyi')) + # def test_Jtvec_adjointTest_eyi_Bform(self): + # self.assertTrue(adjointTest('b', 'eyi')) + # def test_Jtvec_adjointTest_ezi_Eform(self): + # self.assertTrue(adjointTest('e', 'ezi')) + # def test_Jtvec_adjointTest_ezi_Bform(self): + # self.assertTrue(adjointTest('b', 'ezi')) - def test_Jtvec_adjointTest_bxr_Eform(self): - self.assertTrue(adjointTest('e', 'bxr')) - def test_Jtvec_adjointTest_bxr_Bform(self): - self.assertTrue(adjointTest('b', 'bxr')) - def test_Jtvec_adjointTest_byr_Eform(self): - self.assertTrue(adjointTest('e', 'byr')) - def test_Jtvec_adjointTest_byr_Bform(self): - self.assertTrue(adjointTest('b', 'byr')) - def test_Jtvec_adjointTest_bzr_Eform(self): - self.assertTrue(adjointTest('e', 'bzr')) - def test_Jtvec_adjointTest_bzr_Bform(self): - self.assertTrue(adjointTest('b', 'bzr')) - def test_Jtvec_adjointTest_bxi_Eform(self): - self.assertTrue(adjointTest('e', 'bxi')) - def test_Jtvec_adjointTest_bxi_Bform(self): - self.assertTrue(adjointTest('b', 'bxi')) - def test_Jtvec_adjointTest_byi_Eform(self): - self.assertTrue(adjointTest('e', 'byi')) - def test_Jtvec_adjointTest_byi_Bform(self): - self.assertTrue(adjointTest('b', 'byi')) - def test_Jtvec_adjointTest_bzi_Eform(self): - self.assertTrue(adjointTest('e', 'bzi')) - def test_Jtvec_adjointTest_bzi_Bform(self): - self.assertTrue(adjointTest('b', 'bzi')) + # def test_Jtvec_adjointTest_bxr_Eform(self): + # self.assertTrue(adjointTest('e', 'bxr')) + # def test_Jtvec_adjointTest_bxr_Bform(self): + # self.assertTrue(adjointTest('b', 'bxr')) + # def test_Jtvec_adjointTest_byr_Eform(self): + # self.assertTrue(adjointTest('e', 'byr')) + # def test_Jtvec_adjointTest_byr_Bform(self): + # self.assertTrue(adjointTest('b', 'byr')) + # def test_Jtvec_adjointTest_bzr_Eform(self): + # self.assertTrue(adjointTest('e', 'bzr')) + # def test_Jtvec_adjointTest_bzr_Bform(self): + # self.assertTrue(adjointTest('b', 'bzr')) + # def test_Jtvec_adjointTest_bxi_Eform(self): + # self.assertTrue(adjointTest('e', 'bxi')) + # def test_Jtvec_adjointTest_bxi_Bform(self): + # self.assertTrue(adjointTest('b', 'bxi')) + # def test_Jtvec_adjointTest_byi_Eform(self): + # self.assertTrue(adjointTest('e', 'byi')) + # def test_Jtvec_adjointTest_byi_Bform(self): + # self.assertTrue(adjointTest('b', 'byi')) + # def test_Jtvec_adjointTest_bzi_Eform(self): + # self.assertTrue(adjointTest('e', 'bzi')) + # def test_Jtvec_adjointTest_bzi_Bform(self): + # self.assertTrue(adjointTest('b', 'bzi')) - def test_Jvec_jxr_Jform(self): - self.assertTrue(derivTest('j', 'jxr')) - def test_Jvec_jyr_Jform(self): - self.assertTrue(derivTest('j', 'jyr')) - def test_Jvec_jzr_Jform(self): - self.assertTrue(derivTest('j', 'jzr')) - def test_Jvec_jxi_Jform(self): - self.assertTrue(derivTest('j', 'jxi')) - def test_Jvec_jyi_Jform(self): - self.assertTrue(derivTest('j', 'jyi')) - def test_Jvec_jzi_Jform(self): - self.assertTrue(derivTest('j', 'jzi')) + # def test_Jvec_jxr_Jform(self): + # self.assertTrue(derivTest('j', 'jxr')) + # def test_Jvec_jyr_Jform(self): + # self.assertTrue(derivTest('j', 'jyr')) + # def test_Jvec_jzr_Jform(self): + # self.assertTrue(derivTest('j', 'jzr')) + # def test_Jvec_jxi_Jform(self): + # self.assertTrue(derivTest('j', 'jxi')) + # def test_Jvec_jyi_Jform(self): + # self.assertTrue(derivTest('j', 'jyi')) + # def test_Jvec_jzi_Jform(self): + # self.assertTrue(derivTest('j', 'jzi')) - def test_Jvec_hxr_Jform(self): - self.assertTrue(derivTest('j', 'hxr')) - def test_Jvec_hyr_Jform(self): - self.assertTrue(derivTest('j', 'hyr')) - def test_Jvec_hzr_Jform(self): - self.assertTrue(derivTest('j', 'hzr')) - def test_Jvec_hxi_Jform(self): - self.assertTrue(derivTest('j', 'hxi')) - def test_Jvec_hyi_Jform(self): - self.assertTrue(derivTest('j', 'hyi')) - def test_Jvec_hzi_Jform(self): - self.assertTrue(derivTest('j', 'hzi')) + # def test_Jvec_hxr_Jform(self): + # self.assertTrue(derivTest('j', 'hxr')) + # def test_Jvec_hyr_Jform(self): + # self.assertTrue(derivTest('j', 'hyr')) + # def test_Jvec_hzr_Jform(self): + # self.assertTrue(derivTest('j', 'hzr')) + # def test_Jvec_hxi_Jform(self): + # self.assertTrue(derivTest('j', 'hxi')) + # def test_Jvec_hyi_Jform(self): + # self.assertTrue(derivTest('j', 'hyi')) + # def test_Jvec_hzi_Jform(self): + # self.assertTrue(derivTest('j', 'hzi')) - def test_Jvec_hxr_Hform(self): - self.assertTrue(derivTest('h', 'hxr')) - def test_Jvec_hyr_Hform(self): - self.assertTrue(derivTest('h', 'hyr')) - def test_Jvec_hzr_Hform(self): - self.assertTrue(derivTest('h', 'hzr')) - def test_Jvec_hxi_Hform(self): - self.assertTrue(derivTest('h', 'hxi')) - def test_Jvec_hyi_Hform(self): - self.assertTrue(derivTest('h', 'hyi')) - def test_Jvec_hzi_Hform(self): - self.assertTrue(derivTest('h', 'hzi')) + # def test_Jvec_hxr_Hform(self): + # self.assertTrue(derivTest('h', 'hxr')) + # def test_Jvec_hyr_Hform(self): + # self.assertTrue(derivTest('h', 'hyr')) + # def test_Jvec_hzr_Hform(self): + # self.assertTrue(derivTest('h', 'hzr')) + # def test_Jvec_hxi_Hform(self): + # self.assertTrue(derivTest('h', 'hxi')) + # def test_Jvec_hyi_Hform(self): + # self.assertTrue(derivTest('h', 'hyi')) + # def test_Jvec_hzi_Hform(self): + # self.assertTrue(derivTest('h', 'hzi')) - def test_Jvec_hxr_Hform(self): - self.assertTrue(derivTest('h', 'jxr')) - def test_Jvec_hyr_Hform(self): - self.assertTrue(derivTest('h', 'jyr')) - def test_Jvec_hzr_Hform(self): - self.assertTrue(derivTest('h', 'jzr')) - def test_Jvec_hxi_Hform(self): - self.assertTrue(derivTest('h', 'jxi')) - def test_Jvec_hyi_Hform(self): - self.assertTrue(derivTest('h', 'jyi')) - def test_Jvec_hzi_Hform(self): - self.assertTrue(derivTest('h', 'jzi')) + # def test_Jvec_hxr_Hform(self): + # self.assertTrue(derivTest('h', 'jxr')) + # def test_Jvec_hyr_Hform(self): + # self.assertTrue(derivTest('h', 'jyr')) + # def test_Jvec_hzr_Hform(self): + # self.assertTrue(derivTest('h', 'jzr')) + # def test_Jvec_hxi_Hform(self): + # self.assertTrue(derivTest('h', 'jxi')) + # def test_Jvec_hyi_Hform(self): + # self.assertTrue(derivTest('h', 'jyi')) + # def test_Jvec_hzi_Hform(self): + # self.assertTrue(derivTest('h', 'jzi')) - def test_Jtvec_adjointTest_jxr_Jform(self): - self.assertTrue(adjointTest('j', 'jxr')) - def test_Jtvec_adjointTest_jyr_Jform(self): - self.assertTrue(adjointTest('j', 'jyr')) - def test_Jtvec_adjointTest_jzr_Jform(self): - self.assertTrue(adjointTest('j', 'jzr')) - def test_Jtvec_adjointTest_jxi_Jform(self): - self.assertTrue(adjointTest('j', 'jxi')) - def test_Jtvec_adjointTest_jyi_Jform(self): - self.assertTrue(adjointTest('j', 'jyi')) - def test_Jtvec_adjointTest_jzi_Jform(self): - self.assertTrue(adjointTest('j', 'jzi')) + # def test_Jtvec_adjointTest_jxr_Jform(self): + # self.assertTrue(adjointTest('j', 'jxr')) + # def test_Jtvec_adjointTest_jyr_Jform(self): + # self.assertTrue(adjointTest('j', 'jyr')) + # def test_Jtvec_adjointTest_jzr_Jform(self): + # self.assertTrue(adjointTest('j', 'jzr')) + # def test_Jtvec_adjointTest_jxi_Jform(self): + # self.assertTrue(adjointTest('j', 'jxi')) + # def test_Jtvec_adjointTest_jyi_Jform(self): + # self.assertTrue(adjointTest('j', 'jyi')) + # def test_Jtvec_adjointTest_jzi_Jform(self): + # self.assertTrue(adjointTest('j', 'jzi')) - def test_Jtvec_adjointTest_hxr_Jform(self): - self.assertTrue(adjointTest('j', 'hxr')) - def test_Jtvec_adjointTest_hyr_Jform(self): - self.assertTrue(adjointTest('j', 'hyr')) - def test_Jtvec_adjointTest_hzr_Jform(self): - self.assertTrue(adjointTest('j', 'hzr')) - def test_Jtvec_adjointTest_hxi_Jform(self): - self.assertTrue(adjointTest('j', 'hxi')) - def test_Jtvec_adjointTest_hyi_Jform(self): - self.assertTrue(adjointTest('j', 'hyi')) - def test_Jtvec_adjointTest_hzi_Jform(self): - self.assertTrue(adjointTest('j', 'hzi')) + # def test_Jtvec_adjointTest_hxr_Jform(self): + # self.assertTrue(adjointTest('j', 'hxr')) + # def test_Jtvec_adjointTest_hyr_Jform(self): + # self.assertTrue(adjointTest('j', 'hyr')) + # def test_Jtvec_adjointTest_hzr_Jform(self): + # self.assertTrue(adjointTest('j', 'hzr')) + # def test_Jtvec_adjointTest_hxi_Jform(self): + # self.assertTrue(adjointTest('j', 'hxi')) + # def test_Jtvec_adjointTest_hyi_Jform(self): + # self.assertTrue(adjointTest('j', 'hyi')) + # def test_Jtvec_adjointTest_hzi_Jform(self): + # 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_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_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_EB_CrossCheck_exr_Eform(self): + self.assertTrue(crossCheckTest('e', 'exr')) + def test_EB_CrossCheck_eyr_Eform(self): + self.assertTrue(crossCheckTest('e', 'eyr')) + def test_EB_CrossCheck_ezr_Eform(self): + self.assertTrue(crossCheckTest('e', 'ezr')) + def test_EB_CrossCheck_exi_Eform(self): + self.assertTrue(crossCheckTest('e', 'exi')) + def test_EB_CrossCheck_eyi_Eform(self): + self.assertTrue(crossCheckTest('e', 'eyi')) + def test_EB_CrossCheck_ezi_Eform(self): + self.assertTrue(crossCheckTest('e', 'ezi')) + + def test_EB_CrossCheck_bxr_Eform(self): + self.assertTrue(crossCheckTest('e', 'bxr')) + def test_EB_CrossCheck_byr_Eform(self): + self.assertTrue(crossCheckTest('e', 'byr')) + def test_EB_CrossCheck_bzr_Eform(self): + self.assertTrue(crossCheckTest('e', 'bzr')) + def test_EB_CrossCheck_bxi_Eform(self): + self.assertTrue(crossCheckTest('e', 'bxi')) + def test_EB_CrossCheck_byi_Eform(self): + self.assertTrue(crossCheckTest('e', 'byi')) + def test_EB_CrossCheck_bzi_Eform(self): + self.assertTrue(crossCheckTest('e', 'bzi')) + + def test_HJ_CrossCheck_jxr_Jform(self): + self.assertTrue(crossCheckTest('j', 'jxr')) + def test_HJ_CrossCheck_jyr_Jform(self): + self.assertTrue(crossCheckTest('j', 'jyr')) + def test_HJ_CrossCheck_jzr_Jform(self): + self.assertTrue(crossCheckTest('j', 'jzr')) + def test_HJ_CrossCheck_jxi_Jform(self): + self.assertTrue(crossCheckTest('j', 'jxi')) + def test_HJ_CrossCheck_jyi_Jform(self): + self.assertTrue(crossCheckTest('j', 'jyi')) + def test_HJ_CrossCheck_jzi_Jform(self): + self.assertTrue(crossCheckTest('j', 'jzi')) + + def test_HJ_CrossCheck_hxr_Jform(self): + self.assertTrue(crossCheckTest('j', 'hxr')) + def test_HJ_CrossCheck_hyr_Jform(self): + self.assertTrue(crossCheckTest('j', 'hyr')) + def test_HJ_CrossCheck_hzr_Jform(self): + self.assertTrue(crossCheckTest('j', 'hzr')) + def test_HJ_CrossCheck_hxi_Jform(self): + self.assertTrue(crossCheckTest('j', 'hxi')) + def test_HJ_CrossCheck_hyi_Jform(self): + self.assertTrue(crossCheckTest('j', 'hyi')) + def test_HJ_CrossCheck_hzi_Jform(self): + self.assertTrue(crossCheckTest('j', 'hzi')) if __name__ == '__main__': unittest.main() From 5c19da60cc792745e47aa08c1c625226b7efcd36 Mon Sep 17 00:00:00 2001 From: Lindsey Date: Fri, 27 Feb 2015 12:13:01 -0800 Subject: [PATCH 16/25] fixed MeMuI, should be invmat=True, as opposed to 1/mu --- simpegEM/Base.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/simpegEM/Base.py b/simpegEM/Base.py index 38fca561..459a4e1e 100644 --- a/simpegEM/Base.py +++ b/simpegEM/Base.py @@ -52,7 +52,7 @@ class BaseEMProblem(Problem.BaseProblem): def MeMuI(self): # Assuming isotropic mu! if getattr(self, '_MeMuI', None) is None: - self._MeMuI = self.mesh.getEdgeInnerProduct(1/self.mu) + self._MeMuI = self.mesh.getEdgeInnerProduct(self.mu, invMat=True) return self._MeMuI @property From 2780a4cea7d3366f41cb0113376cf26a9a128047 Mon Sep 17 00:00:00 2001 From: Lindsey Date: Fri, 27 Feb 2015 13:55:18 -0800 Subject: [PATCH 17/25] HJ formulation implemented. EB formlation will fail the cross-check test... maybe a problem in the source definition? --- simpegEM/FDEM/FDEM.py | 42 +--- simpegEM/Tests/test_FDEM.py | 395 ++++++++++++++++++------------------ 2 files changed, 203 insertions(+), 234 deletions(-) diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py index 4638936b..6f74eb9c 100644 --- a/simpegEM/FDEM/FDEM.py +++ b/simpegEM/FDEM/FDEM.py @@ -453,13 +453,13 @@ class ProblemFDEM_j(BaseFDEMProblem): if fieldType == 'j': return j elif fieldType == 'h': - mui = self.MeMuI + MeMuI = self.MeMuI C = self.mesh.edgeCurl MfSigi = self.MfSigmai if not adjoint: - h = -(1./(1j*omega(freq))) * mui * ( C.T * ( MfSigi * j ) ) + h = -(1./(1j*omega(freq))) * MeMuI * ( C.T * ( MfSigi * j ) ) else: - h = -(1./(1j*omega(freq))) * MfSigi * ( C * ( mui.T * j ) ) + h = -(1./(1j*omega(freq))) * MfSigi.T * ( C * ( MeMuI.T * j ) ) return h raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) @@ -618,8 +618,7 @@ class ProblemFDEM_h(BaseFDEMProblem): # MfSigi = self.MfSigmai return C.T*h - # return -1j * omega(freq) * MeMu.T * (MfSigmaiinv * (CTinv * h)) - return C*h #- j_s # -iomega(freq) inv(MfSigmai) inv(C.T) MeMu + return C*h - j_s elif fieldType == 'h': return h raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) @@ -630,41 +629,12 @@ class ProblemFDEM_h(BaseFDEMProblem): if fieldType == 'j': C = self.mesh.edgeCurl - # MeMu = self.MeMu - # MfSigi = self.MfSigmai - - # if adjoint: - # MfSigiTCinv = self.Solver(MfSigi.T*C, **self.solverOpts) - # MeMu.T * (Cinv * (MfSigmaiTinv * h)) - # v1 = MeMu.T *( Cinv *(MfSigmaiTinv * v)) - # pt1 = -1j * omega(freq) * self.calcFieldsDeriv(h,freq,'h',v,adjoint=True) - - # pt2 = 1j * omega(freq) * dsig_dm.T * ( dsigi_dsig.T * ( dMf_dsigi.T * ( MfSigiTCinv * v) ) ) - # return pt1 + pt2 - - # CTMfSigiinv = self.Solver(C.T*MfSigi, **self.solverOpts) - # hDeriv = self.calcFieldsDeriv(h,freq,'h',v,adjoint=False) - - # pt1 = -1j * omega(freq) * (CTMfSigiinv * (MeMu * hDeriv)) - - # sig = self.curModel.transform - # sigi = 1/sig - # dsig_dm = self.curModel.transformDeriv - # dsigi_dsig = -Utils.sdiag(sigi)**2 - - # v1 = CTMfSigiinv * (MeMu * h) - - # dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(v1) - - # pt2 = 1j * omega(freq) * (CTMfSigiinv * (dMf_dsigi * (dsigi_dsig * (dsig_dm * v)))) - - # return pt1+pt2 + j_s = self.getjs(freq) if adjoint: dh = self.calcFieldsDeriv(h,freq,'h',C.T*v,adjoint=True) return dh dh = self.calcFieldsDeriv(h,freq,'h',v) - return C*dh - raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) + return C*dh - j_s elif fieldType == 'h': if adjoint: diff --git a/simpegEM/Tests/test_FDEM.py b/simpegEM/Tests/test_FDEM.py index ffbed6a6..a59cd470 100644 --- a/simpegEM/Tests/test_FDEM.py +++ b/simpegEM/Tests/test_FDEM.py @@ -5,10 +5,10 @@ import sys from scipy.constants import mu_0 TOL = 1e-4 -FLR = 1e-15 # "zero", so if residual below this --> pass regardless of order +FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order CONDUCTIVITY = 1e1 MU = mu_0 -freq = 1 +freq = 1e-1 addrandoms = True # important to addrandoms if testing HJ formulation with VMD source! (or else jz ~ 0) def getProblem(fdemType, comp): @@ -22,14 +22,13 @@ def getProblem(fdemType, comp): mapping = Maps.ExpMap(mesh) - x = np.linspace(-30,30,6) - XYZ = Utils.ndgrid(x,x,np.r_[0]) + x = np.array([np.linspace(-30,-15,3),np.linspace(15,30,3)]) #don't sample right by the transmitter + XYZ = Utils.ndgrid(x,x,np.r_[0.]) Rx0 = EM.FDEM.RxFDEM(XYZ, comp) Tx0 = EM.FDEM.TxFDEM(np.r_[0.,0.,0.], 'VMD', freq, [Rx0]) survey = EM.FDEM.SurveyFDEM([Tx0]) - print fdemType if fdemType == 'e': prb = EM.FDEM.ProblemFDEM_e(mesh, mapping=mapping) @@ -99,7 +98,7 @@ def crossCheckTest(fdemType, comp): prb1 = getProblem(fdemType, comp) mesh = prb1.mesh - print '%s formulation - %s' % (fdemType, comp) + print 'Cross Checking Forward: %s formulation - %s' % (fdemType, comp) m = np.log(np.ones(mesh.nC)*CONDUCTIVITY) mu = np.log(np.ones(mesh.nC)*MU) @@ -145,212 +144,212 @@ def crossCheckTest(fdemType, comp): class FDEM_DerivTests(unittest.TestCase): - # def test_Jvec_exr_Eform(self): - # self.assertTrue(derivTest('e', 'exr')) - # def test_Jvec_exr_Bform(self): - # self.assertTrue(derivTest('b', 'exr')) - # def test_Jvec_eyr_Eform(self): - # self.assertTrue(derivTest('e', 'eyr')) - # def test_Jvec_eyr_Bform(self): - # self.assertTrue(derivTest('b', 'eyr')) - # def test_Jvec_ezr_Eform(self): - # self.assertTrue(derivTest('e', 'ezr')) - # def test_Jvec_ezr_Bform(self): - # self.assertTrue(derivTest('b', 'ezr')) - # def test_Jvec_exi_Eform(self): - # self.assertTrue(derivTest('e', 'exi')) - # def test_Jvec_exi_Bform(self): - # self.assertTrue(derivTest('b', 'exi')) - # def test_Jvec_eyi_Eform(self): - # self.assertTrue(derivTest('e', 'eyi')) - # def test_Jvec_eyi_Bform(self): - # self.assertTrue(derivTest('b', 'eyi')) - # def test_Jvec_ezi_Eform(self): - # self.assertTrue(derivTest('e', 'ezi')) - # def test_Jvec_ezi_Bform(self): - # self.assertTrue(derivTest('b', 'ezi')) + def test_Jvec_exr_Eform(self): + self.assertTrue(derivTest('e', 'exr')) + def test_Jvec_exr_Bform(self): + self.assertTrue(derivTest('b', 'exr')) + def test_Jvec_eyr_Eform(self): + self.assertTrue(derivTest('e', 'eyr')) + def test_Jvec_eyr_Bform(self): + self.assertTrue(derivTest('b', 'eyr')) + def test_Jvec_ezr_Eform(self): + self.assertTrue(derivTest('e', 'ezr')) + def test_Jvec_ezr_Bform(self): + self.assertTrue(derivTest('b', 'ezr')) + def test_Jvec_exi_Eform(self): + self.assertTrue(derivTest('e', 'exi')) + def test_Jvec_exi_Bform(self): + self.assertTrue(derivTest('b', 'exi')) + def test_Jvec_eyi_Eform(self): + self.assertTrue(derivTest('e', 'eyi')) + def test_Jvec_eyi_Bform(self): + self.assertTrue(derivTest('b', 'eyi')) + def test_Jvec_ezi_Eform(self): + self.assertTrue(derivTest('e', 'ezi')) + def test_Jvec_ezi_Bform(self): + self.assertTrue(derivTest('b', 'ezi')) - # def test_Jvec_bxr_Eform(self): - # self.assertTrue(derivTest('e', 'bxr')) - # def test_Jvec_bxr_Bform(self): - # self.assertTrue(derivTest('b', 'bxr')) - # def test_Jvec_byr_Eform(self): - # self.assertTrue(derivTest('e', 'byr')) - # def test_Jvec_byr_Bform(self): - # self.assertTrue(derivTest('b', 'byr')) - # def test_Jvec_bzr_Eform(self): - # self.assertTrue(derivTest('e', 'bzr')) - # def test_Jvec_bzr_Bform(self): - # self.assertTrue(derivTest('b', 'bzr')) - # def test_Jvec_bxi_Eform(self): - # self.assertTrue(derivTest('e', 'bxi')) - # def test_Jvec_bxi_Bform(self): - # self.assertTrue(derivTest('b', 'bxi')) - # def test_Jvec_byi_Eform(self): - # self.assertTrue(derivTest('e', 'byi')) - # def test_Jvec_byi_Bform(self): - # self.assertTrue(derivTest('b', 'byi')) - # def test_Jvec_bzi_Eform(self): - # self.assertTrue(derivTest('e', 'bzi')) - # def test_Jvec_bzi_Bform(self): - # self.assertTrue(derivTest('b', 'bzi')) + def test_Jvec_bxr_Eform(self): + self.assertTrue(derivTest('e', 'bxr')) + def test_Jvec_bxr_Bform(self): + self.assertTrue(derivTest('b', 'bxr')) + def test_Jvec_byr_Eform(self): + self.assertTrue(derivTest('e', 'byr')) + def test_Jvec_byr_Bform(self): + self.assertTrue(derivTest('b', 'byr')) + def test_Jvec_bzr_Eform(self): + self.assertTrue(derivTest('e', 'bzr')) + def test_Jvec_bzr_Bform(self): + self.assertTrue(derivTest('b', 'bzr')) + def test_Jvec_bxi_Eform(self): + self.assertTrue(derivTest('e', 'bxi')) + def test_Jvec_bxi_Bform(self): + self.assertTrue(derivTest('b', 'bxi')) + def test_Jvec_byi_Eform(self): + self.assertTrue(derivTest('e', 'byi')) + def test_Jvec_byi_Bform(self): + self.assertTrue(derivTest('b', 'byi')) + def test_Jvec_bzi_Eform(self): + self.assertTrue(derivTest('e', 'bzi')) + def test_Jvec_bzi_Bform(self): + self.assertTrue(derivTest('b', 'bzi')) - # def test_Jtvec_adjointTest_exr_Eform(self): - # self.assertTrue(adjointTest('e', 'exr')) - # def test_Jtvec_adjointTest_exr_Bform(self): - # self.assertTrue(adjointTest('b', 'exr')) - # def test_Jtvec_adjointTest_eyr_Eform(self): - # self.assertTrue(adjointTest('e', 'eyr')) - # def test_Jtvec_adjointTest_eyr_Bform(self): - # self.assertTrue(adjointTest('b', 'eyr')) - # def test_Jtvec_adjointTest_ezr_Eform(self): - # self.assertTrue(adjointTest('e', 'ezr')) - # def test_Jtvec_adjointTest_ezr_Bform(self): - # self.assertTrue(adjointTest('b', 'ezr')) - # def test_Jtvec_adjointTest_exi_Eform(self): - # self.assertTrue(adjointTest('e', 'exi')) - # def test_Jtvec_adjointTest_exi_Bform(self): - # self.assertTrue(adjointTest('b', 'exi')) - # def test_Jtvec_adjointTest_eyi_Eform(self): - # self.assertTrue(adjointTest('e', 'eyi')) - # def test_Jtvec_adjointTest_eyi_Bform(self): - # self.assertTrue(adjointTest('b', 'eyi')) - # def test_Jtvec_adjointTest_ezi_Eform(self): - # self.assertTrue(adjointTest('e', 'ezi')) - # def test_Jtvec_adjointTest_ezi_Bform(self): - # self.assertTrue(adjointTest('b', 'ezi')) + def test_Jtvec_adjointTest_exr_Eform(self): + self.assertTrue(adjointTest('e', 'exr')) + def test_Jtvec_adjointTest_exr_Bform(self): + self.assertTrue(adjointTest('b', 'exr')) + def test_Jtvec_adjointTest_eyr_Eform(self): + self.assertTrue(adjointTest('e', 'eyr')) + def test_Jtvec_adjointTest_eyr_Bform(self): + self.assertTrue(adjointTest('b', 'eyr')) + def test_Jtvec_adjointTest_ezr_Eform(self): + self.assertTrue(adjointTest('e', 'ezr')) + def test_Jtvec_adjointTest_ezr_Bform(self): + self.assertTrue(adjointTest('b', 'ezr')) + def test_Jtvec_adjointTest_exi_Eform(self): + self.assertTrue(adjointTest('e', 'exi')) + def test_Jtvec_adjointTest_exi_Bform(self): + self.assertTrue(adjointTest('b', 'exi')) + def test_Jtvec_adjointTest_eyi_Eform(self): + self.assertTrue(adjointTest('e', 'eyi')) + def test_Jtvec_adjointTest_eyi_Bform(self): + self.assertTrue(adjointTest('b', 'eyi')) + def test_Jtvec_adjointTest_ezi_Eform(self): + self.assertTrue(adjointTest('e', 'ezi')) + def test_Jtvec_adjointTest_ezi_Bform(self): + self.assertTrue(adjointTest('b', 'ezi')) - # def test_Jtvec_adjointTest_bxr_Eform(self): - # self.assertTrue(adjointTest('e', 'bxr')) - # def test_Jtvec_adjointTest_bxr_Bform(self): - # self.assertTrue(adjointTest('b', 'bxr')) - # def test_Jtvec_adjointTest_byr_Eform(self): - # self.assertTrue(adjointTest('e', 'byr')) - # def test_Jtvec_adjointTest_byr_Bform(self): - # self.assertTrue(adjointTest('b', 'byr')) - # def test_Jtvec_adjointTest_bzr_Eform(self): - # self.assertTrue(adjointTest('e', 'bzr')) - # def test_Jtvec_adjointTest_bzr_Bform(self): - # self.assertTrue(adjointTest('b', 'bzr')) - # def test_Jtvec_adjointTest_bxi_Eform(self): - # self.assertTrue(adjointTest('e', 'bxi')) - # def test_Jtvec_adjointTest_bxi_Bform(self): - # self.assertTrue(adjointTest('b', 'bxi')) - # def test_Jtvec_adjointTest_byi_Eform(self): - # self.assertTrue(adjointTest('e', 'byi')) - # def test_Jtvec_adjointTest_byi_Bform(self): - # self.assertTrue(adjointTest('b', 'byi')) - # def test_Jtvec_adjointTest_bzi_Eform(self): - # self.assertTrue(adjointTest('e', 'bzi')) - # def test_Jtvec_adjointTest_bzi_Bform(self): - # self.assertTrue(adjointTest('b', 'bzi')) + def test_Jtvec_adjointTest_bxr_Eform(self): + self.assertTrue(adjointTest('e', 'bxr')) + def test_Jtvec_adjointTest_bxr_Bform(self): + self.assertTrue(adjointTest('b', 'bxr')) + def test_Jtvec_adjointTest_byr_Eform(self): + self.assertTrue(adjointTest('e', 'byr')) + def test_Jtvec_adjointTest_byr_Bform(self): + self.assertTrue(adjointTest('b', 'byr')) + def test_Jtvec_adjointTest_bzr_Eform(self): + self.assertTrue(adjointTest('e', 'bzr')) + def test_Jtvec_adjointTest_bzr_Bform(self): + self.assertTrue(adjointTest('b', 'bzr')) + def test_Jtvec_adjointTest_bxi_Eform(self): + self.assertTrue(adjointTest('e', 'bxi')) + def test_Jtvec_adjointTest_bxi_Bform(self): + self.assertTrue(adjointTest('b', 'bxi')) + def test_Jtvec_adjointTest_byi_Eform(self): + self.assertTrue(adjointTest('e', 'byi')) + def test_Jtvec_adjointTest_byi_Bform(self): + self.assertTrue(adjointTest('b', 'byi')) + def test_Jtvec_adjointTest_bzi_Eform(self): + self.assertTrue(adjointTest('e', 'bzi')) + def test_Jtvec_adjointTest_bzi_Bform(self): + self.assertTrue(adjointTest('b', 'bzi')) - # def test_Jvec_jxr_Jform(self): - # self.assertTrue(derivTest('j', 'jxr')) - # def test_Jvec_jyr_Jform(self): - # self.assertTrue(derivTest('j', 'jyr')) - # def test_Jvec_jzr_Jform(self): - # self.assertTrue(derivTest('j', 'jzr')) - # def test_Jvec_jxi_Jform(self): - # self.assertTrue(derivTest('j', 'jxi')) - # def test_Jvec_jyi_Jform(self): - # self.assertTrue(derivTest('j', 'jyi')) - # def test_Jvec_jzi_Jform(self): - # self.assertTrue(derivTest('j', 'jzi')) + def test_Jvec_jxr_Jform(self): + self.assertTrue(derivTest('j', 'jxr')) + def test_Jvec_jyr_Jform(self): + self.assertTrue(derivTest('j', 'jyr')) + def test_Jvec_jzr_Jform(self): + self.assertTrue(derivTest('j', 'jzr')) + def test_Jvec_jxi_Jform(self): + self.assertTrue(derivTest('j', 'jxi')) + def test_Jvec_jyi_Jform(self): + self.assertTrue(derivTest('j', 'jyi')) + def test_Jvec_jzi_Jform(self): + self.assertTrue(derivTest('j', 'jzi')) - # def test_Jvec_hxr_Jform(self): - # self.assertTrue(derivTest('j', 'hxr')) - # def test_Jvec_hyr_Jform(self): - # self.assertTrue(derivTest('j', 'hyr')) - # def test_Jvec_hzr_Jform(self): - # self.assertTrue(derivTest('j', 'hzr')) - # def test_Jvec_hxi_Jform(self): - # self.assertTrue(derivTest('j', 'hxi')) - # def test_Jvec_hyi_Jform(self): - # self.assertTrue(derivTest('j', 'hyi')) - # def test_Jvec_hzi_Jform(self): - # self.assertTrue(derivTest('j', 'hzi')) + def test_Jvec_hxr_Jform(self): + self.assertTrue(derivTest('j', 'hxr')) + def test_Jvec_hyr_Jform(self): + self.assertTrue(derivTest('j', 'hyr')) + def test_Jvec_hzr_Jform(self): + self.assertTrue(derivTest('j', 'hzr')) + def test_Jvec_hxi_Jform(self): + self.assertTrue(derivTest('j', 'hxi')) + def test_Jvec_hyi_Jform(self): + self.assertTrue(derivTest('j', 'hyi')) + def test_Jvec_hzi_Jform(self): + self.assertTrue(derivTest('j', 'hzi')) - # def test_Jvec_hxr_Hform(self): - # self.assertTrue(derivTest('h', 'hxr')) - # def test_Jvec_hyr_Hform(self): - # self.assertTrue(derivTest('h', 'hyr')) - # def test_Jvec_hzr_Hform(self): - # self.assertTrue(derivTest('h', 'hzr')) - # def test_Jvec_hxi_Hform(self): - # self.assertTrue(derivTest('h', 'hxi')) - # def test_Jvec_hyi_Hform(self): - # self.assertTrue(derivTest('h', 'hyi')) - # def test_Jvec_hzi_Hform(self): - # self.assertTrue(derivTest('h', 'hzi')) + def test_Jvec_hxr_Hform(self): + self.assertTrue(derivTest('h', 'hxr')) + def test_Jvec_hyr_Hform(self): + self.assertTrue(derivTest('h', 'hyr')) + def test_Jvec_hzr_Hform(self): + self.assertTrue(derivTest('h', 'hzr')) + def test_Jvec_hxi_Hform(self): + self.assertTrue(derivTest('h', 'hxi')) + def test_Jvec_hyi_Hform(self): + self.assertTrue(derivTest('h', 'hyi')) + def test_Jvec_hzi_Hform(self): + self.assertTrue(derivTest('h', 'hzi')) - # def test_Jvec_hxr_Hform(self): - # self.assertTrue(derivTest('h', 'jxr')) - # def test_Jvec_hyr_Hform(self): - # self.assertTrue(derivTest('h', 'jyr')) - # def test_Jvec_hzr_Hform(self): - # self.assertTrue(derivTest('h', 'jzr')) - # def test_Jvec_hxi_Hform(self): - # self.assertTrue(derivTest('h', 'jxi')) - # def test_Jvec_hyi_Hform(self): - # self.assertTrue(derivTest('h', 'jyi')) - # def test_Jvec_hzi_Hform(self): - # self.assertTrue(derivTest('h', 'jzi')) + def test_Jvec_hxr_Hform(self): + self.assertTrue(derivTest('h', 'jxr')) + def test_Jvec_hyr_Hform(self): + self.assertTrue(derivTest('h', 'jyr')) + def test_Jvec_hzr_Hform(self): + self.assertTrue(derivTest('h', 'jzr')) + def test_Jvec_hxi_Hform(self): + self.assertTrue(derivTest('h', 'jxi')) + def test_Jvec_hyi_Hform(self): + self.assertTrue(derivTest('h', 'jyi')) + def test_Jvec_hzi_Hform(self): + self.assertTrue(derivTest('h', 'jzi')) - # def test_Jtvec_adjointTest_jxr_Jform(self): - # self.assertTrue(adjointTest('j', 'jxr')) - # def test_Jtvec_adjointTest_jyr_Jform(self): - # self.assertTrue(adjointTest('j', 'jyr')) - # def test_Jtvec_adjointTest_jzr_Jform(self): - # self.assertTrue(adjointTest('j', 'jzr')) - # def test_Jtvec_adjointTest_jxi_Jform(self): - # self.assertTrue(adjointTest('j', 'jxi')) - # def test_Jtvec_adjointTest_jyi_Jform(self): - # self.assertTrue(adjointTest('j', 'jyi')) - # def test_Jtvec_adjointTest_jzi_Jform(self): - # self.assertTrue(adjointTest('j', 'jzi')) + def test_Jtvec_adjointTest_jxr_Jform(self): + self.assertTrue(adjointTest('j', 'jxr')) + def test_Jtvec_adjointTest_jyr_Jform(self): + self.assertTrue(adjointTest('j', 'jyr')) + def test_Jtvec_adjointTest_jzr_Jform(self): + self.assertTrue(adjointTest('j', 'jzr')) + def test_Jtvec_adjointTest_jxi_Jform(self): + self.assertTrue(adjointTest('j', 'jxi')) + def test_Jtvec_adjointTest_jyi_Jform(self): + self.assertTrue(adjointTest('j', 'jyi')) + def test_Jtvec_adjointTest_jzi_Jform(self): + self.assertTrue(adjointTest('j', 'jzi')) - # def test_Jtvec_adjointTest_hxr_Jform(self): - # self.assertTrue(adjointTest('j', 'hxr')) - # def test_Jtvec_adjointTest_hyr_Jform(self): - # self.assertTrue(adjointTest('j', 'hyr')) - # def test_Jtvec_adjointTest_hzr_Jform(self): - # self.assertTrue(adjointTest('j', 'hzr')) - # def test_Jtvec_adjointTest_hxi_Jform(self): - # self.assertTrue(adjointTest('j', 'hxi')) - # def test_Jtvec_adjointTest_hyi_Jform(self): - # self.assertTrue(adjointTest('j', 'hyi')) - # def test_Jtvec_adjointTest_hzi_Jform(self): - # self.assertTrue(adjointTest('j', 'hzi')) + def test_Jtvec_adjointTest_hxr_Jform(self): + self.assertTrue(adjointTest('j', 'hxr')) + def test_Jtvec_adjointTest_hyr_Jform(self): + self.assertTrue(adjointTest('j', 'hyr')) + def test_Jtvec_adjointTest_hzr_Jform(self): + self.assertTrue(adjointTest('j', 'hzr')) + def test_Jtvec_adjointTest_hxi_Jform(self): + self.assertTrue(adjointTest('j', 'hxi')) + def test_Jtvec_adjointTest_hyi_Jform(self): + self.assertTrue(adjointTest('j', 'hyi')) + def test_Jtvec_adjointTest_hzi_Jform(self): + 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_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_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_EB_CrossCheck_exr_Eform(self): From 84fedcf2a7af453c3ecab9981c6991ebffbd1467 Mon Sep 17 00:00:00 2001 From: Lindsey Date: Fri, 27 Feb 2015 15:25:13 -0800 Subject: [PATCH 18/25] seperated out calculation of j_s from calculation of rhs in j formulation --- simpegEM/FDEM/FDEM.py | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py index 6f74eb9c..7bc9abd5 100644 --- a/simpegEM/FDEM/FDEM.py +++ b/simpegEM/FDEM/FDEM.py @@ -415,12 +415,11 @@ class ProblemFDEM_j(BaseFDEMProblem): return C * ( MeMuI * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) ) ) - - def getRHS(self, freq): + def getjs(self,freq): """ :param float freq: Frequency :rtype: numpy.ndarray (nE, nTx) - :return: RHS + :return: j_s """ Txs = self.survey.getTransmitters(freq) rhs = range(len(Txs)) @@ -445,7 +444,15 @@ class ProblemFDEM_j(BaseFDEMProblem): MeMuI = self.MeMuI C = self.mesh.edgeCurl - j_s = C*MeMuI*C.T*a + return C*MeMuI*C.T*a + + def getRHS(self, freq): + """ + :param float freq: Frequency + :rtype: numpy.ndarray (nE, nTx) + :return: RHS + """ + j_s = self.getjs(freq) return -1j*omega(freq)*j_s def calcFields(self, sol, freq, fieldType, adjoint=False): From 62b4edcf6d32cb886406f6d7bc6a7bde1f025595 Mon Sep 17 00:00:00 2001 From: Lindsey Date: Fri, 27 Feb 2015 18:43:47 -0800 Subject: [PATCH 19/25] - changed H implementation to primary secondary for a dipole source (much faster) - removed extra MeMui matrix floating around in the B implementation - changed tests so that we do not cross-check the z-components, as this doesn't make sense for the primary secondary approach if we don't add back the primary --- simpegEM/FDEM/FDEM.py | 95 ++++++++++++++++++------------------- simpegEM/Tests/test_FDEM.py | 10 ++-- 2 files changed, 52 insertions(+), 53 deletions(-) diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py index 7bc9abd5..78bb8a0b 100644 --- a/simpegEM/FDEM/FDEM.py +++ b/simpegEM/FDEM/FDEM.py @@ -223,8 +223,9 @@ class ProblemFDEM_b(BaseFDEMProblem): mui = self.MfMui sigI = self.MeSigmaI C = self.mesh.edgeCurl + iomega = 1j * omega(freq) * sp.eye(self.mesh.nF) - return mui*C*sigI*C.T*mui + 1j*omega(freq)*mui + return C*sigI*C.T*mui + iomega def getADeriv(self, freq, u, v, adjoint=False): @@ -301,7 +302,7 @@ class ProblemFDEM_b(BaseFDEMProblem): C = self.mesh.edgeCurl b_0 = C*a - return -1j*omega(freq)*mui*b_0 + return -1j*omega(freq)*b_0 def calcFields(self, sol, freq, fieldType, adjoint=False): b = sol @@ -492,7 +493,7 @@ class ProblemFDEM_j(BaseFDEMProblem): -# Solving for h! - NOTE: WE ARE GOING TO NEED dRHS / dm ! +# Solving for h! - using primary- secondary approach class ProblemFDEM_h(BaseFDEMProblem): """ Using the H-J formulation of Maxwell's equations @@ -583,7 +584,7 @@ class ProblemFDEM_h(BaseFDEMProblem): MeMuI = self.MeMuI C = self.mesh.edgeCurl - return C*MeMuI*C.T*a + return MeMuI*C.T*a #C*MeMuI*C.T*a def getRHS(self, freq): """ @@ -591,64 +592,62 @@ class ProblemFDEM_h(BaseFDEMProblem): :rtype: numpy.ndarray (nE, nTx) :return: RHS """ + MeMu = self.MeMu MfSigi = self.MfSigmai C = self.mesh.edgeCurl - j_s = self.getjs(freq) - return C.T*MfSigi*j_s + Hp = self.getjs(freq) + return -1j*omega(freq)*MeMu*Hp #C.T*MfSigi*j_s - def getRHSDeriv(self, freq, v, adjoint=False): - """ - :param float freq: Frequency - :rtype: numpy.ndarray (nE, nTx) - :return: RHSDeriv - """ - C = self.mesh.edgeCurl - sig = self.curModel.transform - sigi = 1/sig - j_s = self.getjs(freq) - dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(j_s) - dsig_dm = self.curModel.transformDeriv - dsigi_dsig = -Utils.sdiag(sigi)**2 # only works for diagonal matrices + # def getRHSDeriv(self, freq, v, adjoint=False): + # """ + # :param float freq: Frequency + # :rtype: numpy.ndarray (nE, nTx) + # :return: RHSDeriv + # """ + # C = self.mesh.edgeCurl + # sig = self.curModel.transform + # sigi = 1/sig + # j_s = self.getjs(freq) + # dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(j_s) + # dsig_dm = self.curModel.transformDeriv + # dsigi_dsig = -Utils.sdiag(sigi)**2 # only works for diagonal matrices - if adjoint: - return dsig_dm.T * dsigi_dsig.T * dMf_dsigi.T * C * v - return C.T * dMf_dsigi * dsigi_dsig * dsig_dm * v + # if adjoint: + # return dsig_dm.T * dsigi_dsig.T * dMf_dsigi.T * C * v + # return C.T * dMf_dsigi * dsigi_dsig * dsig_dm * v def calcFields(self, sol, freq, fieldType, adjoint=False): h = sol if fieldType == 'j': - # NotImplementedError('fieldType "%s" is not implemented.' % fieldType) C = self.mesh.edgeCurl - j_s = self.getjs(freq) - if adjoint: - # MeMuI = self.MeMuI - # MfSigi = self.MfSigmai - + # j_s = self.getjs(freq) + if adjoint: return C.T*h - return C*h - j_s + return C*h #- j_s elif fieldType == 'h': return h raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) def calcFieldsDeriv(self, sol, freq, fieldType, v, adjoint=False): - h = sol - A = self.getA(freq) + return None + # h = sol + # A = self.getA(freq) - if fieldType == 'j': - C = self.mesh.edgeCurl - j_s = self.getjs(freq) - if adjoint: - dh = self.calcFieldsDeriv(h,freq,'h',C.T*v,adjoint=True) - return dh - dh = self.calcFieldsDeriv(h,freq,'h',v) - return C*dh - j_s + # if fieldType == 'j': + # C = self.mesh.edgeCurl + # j_s = self.getjs(freq) + # if adjoint: + # dh = self.calcFieldsDeriv(h,freq,'h',C.T*v,adjoint=True) + # return dh + # dh = self.calcFieldsDeriv(h,freq,'h',v) + # return C*dh - j_s - elif fieldType == 'h': - if adjoint: - ATinv = self.Solver(A.T, **self.solverOpts) - ATinvv = ATinv*v - return self.getRHSDeriv(freq,ATinvv,adjoint=True) - dRHSh = self.getRHSDeriv(freq,v,adjoint) - Ainv = self.Solver(A, **self.solverOpts) - return Ainv*dRHSh - raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) \ No newline at end of file + # elif fieldType == 'h': + # if adjoint: + # ATinv = self.Solver(A.T, **self.solverOpts) + # ATinvv = ATinv*v + # return self.getRHSDeriv(freq,ATinvv,adjoint=True) + # dRHSh = self.getRHSDeriv(freq,v,adjoint) + # Ainv = self.Solver(A, **self.solverOpts) + # return Ainv*dRHSh + # raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) \ No newline at end of file diff --git a/simpegEM/Tests/test_FDEM.py b/simpegEM/Tests/test_FDEM.py index a59cd470..4c2ecefd 100644 --- a/simpegEM/Tests/test_FDEM.py +++ b/simpegEM/Tests/test_FDEM.py @@ -323,7 +323,7 @@ class FDEM_DerivTests(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')) @@ -369,8 +369,8 @@ class FDEM_DerivTests(unittest.TestCase): self.assertTrue(crossCheckTest('e', 'bxr')) def test_EB_CrossCheck_byr_Eform(self): self.assertTrue(crossCheckTest('e', 'byr')) - def test_EB_CrossCheck_bzr_Eform(self): - self.assertTrue(crossCheckTest('e', 'bzr')) + # def test_EB_CrossCheck_bzr_Eform(self): + # self.assertTrue(crossCheckTest('e', 'bzr')) # Doesn't make sense to test this for p-s approach def test_EB_CrossCheck_bxi_Eform(self): self.assertTrue(crossCheckTest('e', 'bxi')) def test_EB_CrossCheck_byi_Eform(self): @@ -395,8 +395,8 @@ class FDEM_DerivTests(unittest.TestCase): self.assertTrue(crossCheckTest('j', 'hxr')) def test_HJ_CrossCheck_hyr_Jform(self): self.assertTrue(crossCheckTest('j', 'hyr')) - def test_HJ_CrossCheck_hzr_Jform(self): - self.assertTrue(crossCheckTest('j', 'hzr')) + # def test_HJ_CrossCheck_hzr_Jform(self): + # self.assertTrue(crossCheckTest('j', 'hzr')) # Doesn't make sense to test this for p-s approach def test_HJ_CrossCheck_hxi_Jform(self): self.assertTrue(crossCheckTest('j', 'hxi')) def test_HJ_CrossCheck_hyi_Jform(self): From 35678c587d02bc466598adc060aef76b4606e279 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sat, 28 Feb 2015 08:20:58 -0800 Subject: [PATCH 20/25] Removed extra mui in bformulation deriv terms --- simpegEM/FDEM/FDEM.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py index 78bb8a0b..8f549034 100644 --- a/simpegEM/FDEM/FDEM.py +++ b/simpegEM/FDEM/FDEM.py @@ -240,9 +240,10 @@ class ProblemFDEM_b(BaseFDEMProblem): dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig)(vec) if adjoint: - return dsig_dm.T * ( dMe_dsig.T * ( dMeSigmaI_dI.T * ( C.T * ( mui.T * v ) ) ) ) + return dsig_dm.T * ( dMe_dsig.T * ( dMeSigmaI_dI.T * ( C.T * v ) ) ) + + return C * ( dMeSigmaI_dI * ( dMe_dsig * ( dsig_dm * v ) ) ) - return mui * ( C * ( dMeSigmaI_dI * ( dMe_dsig * ( dsig_dm * v ) ) ) ) def getRHS(self, freq): """ From fa0fd7f23f5e756607c4770f553258a2093ea8f6 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sat, 28 Feb 2015 10:18:58 -0800 Subject: [PATCH 21/25] broke out calculation of source term from rhs so that you can do prb.getSource --- simpegEM/FDEM/FDEM.py | 386 +++++++++++++++++++++--------------- simpegEM/Tests/test_FDEM.py | 2 +- 2 files changed, 224 insertions(+), 164 deletions(-) diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py index 8f549034..e816be19 100644 --- a/simpegEM/FDEM/FDEM.py +++ b/simpegEM/FDEM/FDEM.py @@ -8,6 +8,102 @@ def omega(freq): """Change frequency to angular frequency, omega""" return 2.*np.pi*freq +def getSource(self,freq): + """ + :param float freq: Frequency + :rtype: numpy.ndarray (nE, nTx) + :return: RHS + """ + Txs = self.survey.getTransmitters(freq) + rhs = range(len(Txs)) + + solType = self.solType + + if solType == 'e' or solType == 'b': + gridEJx = self.mesh.gridEx + gridEJy = self.mesh.gridEy + gridEJz = self.mesh.gridEz + nEJ = self.mesh.nE + + gridBHx = self.mesh.gridFx + gridBHy = self.mesh.gridFy + gridBHz = self.mesh.gridFz + nBH = self.mesh.nF + + + C = self.mesh.edgeCurl + mui = self.MfMui + + elif solType == 'h' or solType == 'j': + gridEJx = self.mesh.gridFx + gridEJy = self.mesh.gridFy + gridEJz = self.mesh.gridFz + nEJ = self.mesh.nF + + gridBHx = self.mesh.gridEx + gridBHy = self.mesh.gridEy + gridBHz = self.mesh.gridEz + nBH = self.mesh.nE + + C = self.mesh.edgeCurl.T + mui = self.MeMuI + + else: + NotImplementedError('Only E or F sources') + + for i, tx in enumerate(Txs): + if self.mesh._meshType is 'CYL': + if self.mesh.isSymmetric: + if tx.txType == 'VMD': + SRC = Sources.MagneticDipoleVectorPotential(tx.loc, gridEJy, 'y') + elif tx.txType =='CircularLoop': + SRC = Sources.MagneticLoopVectorPotential(tx.loc, gridEJy, 'y', tx.radius) + else: + raise NotImplementedError('Only VMD and CircularLoop') + else: + raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!') + + elif self.mesh._meshType is 'TENSOR': + + if tx.txType == 'VMD': + src = Sources.MagneticDipoleVectorPotential + SRCx = src(tx.loc, gridEJx, 'x') + SRCy = src(tx.loc, gridEJy, 'y') + SRCz = src(tx.loc, gridEJz, 'z') + + elif tx.txType == 'VMD_B': + src = Sources.MagneticDipoleFields + SRCx = src(tx.loc, gridBHx, 'x') + SRCy = src(tx.loc, gridBHy, 'y') + SRCz = src(tx.loc, gridBHz, 'z') + + elif tx.txType == 'CircularLoop': + src = Sources.MagneticLoopVectorPotential + SRCx = src(tx.loc, gridEJx, 'x', tx.radius) + SRCy = src(tx.loc, gridEJy, 'y', tx.radius) + SRCz = src(tx.loc, gridEJz, 'z', tx.radius) + else: + + raise NotImplemented('%s txType is not implemented' % tx.txType) + SRC = np.concatenate((SRCx, SRCy, SRCz)) + + else: + raise Exception('Unknown mesh for VMD') + + rhs[i] = SRC + + # b-forumlation + if tx.txType == 'VMD_B': + b_0 = np.concatenate(rhs).reshape((nBH, len(Txs)), order='E') + else: + a = np.concatenate(rhs).reshape((nEJ, len(Txs)), order='F') + b_0 = C*a + + if solType == 'b' or solType == 'h': + return b_0 + elif solType == 'e' or solType == 'j': + return C.T*mui*b_0 + class BaseFDEMProblem(BaseEMProblem): """ We start by looking at Maxwell's equations in the electric field \\(\\vec{E}\\) and the magnetic flux density \\(\\vec{B}\\): @@ -104,6 +200,9 @@ class BaseFDEMProblem(BaseEMProblem): return Jtv + def getSource(self,freq): + return self.getSource(freq) + ########################################################################################## ################################ E-B Formulation ######################################### ########################################################################################## @@ -159,29 +258,29 @@ class ProblemFDEM_e(BaseFDEMProblem): :rtype: numpy.ndarray (nE, nTx) :return: RHS """ - Txs = self.survey.getTransmitters(freq) - rhs = range(len(Txs)) - for i, tx in enumerate(Txs): - if tx.txType == 'VMD': - src = Sources.MagneticDipoleVectorPotential - SRCx = src(tx.loc, self.mesh.gridEx, 'x') - SRCy = src(tx.loc, self.mesh.gridEy, 'y') - SRCz = src(tx.loc, self.mesh.gridEz, 'z') + # Txs = self.survey.getTransmitters(freq) + # rhs = range(len(Txs)) + # for i, tx in enumerate(Txs): + # if tx.txType == 'VMD': + # src = Sources.MagneticDipoleVectorPotential + # SRCx = src(tx.loc, self.mesh.gridEx, 'x') + # SRCy = src(tx.loc, self.mesh.gridEy, 'y') + # SRCz = src(tx.loc, self.mesh.gridEz, 'z') - elif tx.txType == 'CircularLoop': - src = Sources.MagneticLoopVectorPotential - SRCx = src(tx.loc, self.mesh.gridEx, 'x', tx.radius) - SRCy = src(tx.loc, self.mesh.gridEy, 'y', tx.radius) - SRCz = src(tx.loc, self.mesh.gridEz, 'z', tx.radius) - else: - raise NotImplemented('%s txType is not implemented' % tx.txType) - rhs[i] = np.concatenate((SRCx, SRCy, SRCz)) + # elif tx.txType == 'CircularLoop': + # src = Sources.MagneticLoopVectorPotential + # SRCx = src(tx.loc, self.mesh.gridEx, 'x', tx.radius) + # SRCy = src(tx.loc, self.mesh.gridEy, 'y', tx.radius) + # SRCz = src(tx.loc, self.mesh.gridEz, 'z', tx.radius) + # else: + # raise NotImplemented('%s txType is not implemented' % tx.txType) + # rhs[i] = np.concatenate((SRCx, SRCy, SRCz)) - a = np.concatenate(rhs).reshape((self.mesh.nE, len(Txs)), order='F') - mui = self.MfMui - C = self.mesh.edgeCurl + # a = np.concatenate(rhs).reshape((self.mesh.nE, len(Txs)), order='F') + # mui = self.MfMui + # C = self.mesh.edgeCurl - j_s = C.T*mui*C*a + j_s = getSource(self,freq) #C.T*mui*C*a return -1j*omega(freq)*j_s def calcFields(self, sol, freq, fieldType, adjoint=False): @@ -251,57 +350,57 @@ class ProblemFDEM_b(BaseFDEMProblem): :rtype: numpy.ndarray (nE, nTx) :return: RHS """ - Txs = self.survey.getTransmitters(freq) - rhs = range(len(Txs)) - for i, tx in enumerate(Txs): + # Txs = self.survey.getTransmitters(freq) + # rhs = range(len(Txs)) + # for i, tx in enumerate(Txs): - if self.mesh._meshType is 'CYL': - if self.mesh.isSymmetric: - if tx.txType == 'VMD': - SRC = Sources.MagneticDipoleVectorPotential(tx.loc, self.mesh.gridEy, 'y') - elif tx.txType =='CircularLoop': - SRC = Sources.MagneticLoopVectorPotential(tx.loc, self.mesh.gridEy, 'y', tx.radius) - else: - raise NotImplementedError('Only VMD and CircularLoop') - else: - raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!') + # if self.mesh._meshType is 'CYL': + # if self.mesh.isSymmetric: + # if tx.txType == 'VMD': + # SRC = Sources.MagneticDipoleVectorPotential(tx.loc, self.mesh.gridEy, 'y') + # elif tx.txType =='CircularLoop': + # SRC = Sources.MagneticLoopVectorPotential(tx.loc, self.mesh.gridEy, 'y', tx.radius) + # else: + # raise NotImplementedError('Only VMD and CircularLoop') + # else: + # raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!') - elif self.mesh._meshType is 'TENSOR': + # elif self.mesh._meshType is 'TENSOR': - if tx.txType == 'VMD': - src = Sources.MagneticDipoleVectorPotential - SRCx = src(tx.loc, self.mesh.gridEx, 'x') - SRCy = src(tx.loc, self.mesh.gridEy, 'y') - SRCz = src(tx.loc, self.mesh.gridEz, 'z') + # if tx.txType == 'VMD': + # src = Sources.MagneticDipoleVectorPotential + # SRCx = src(tx.loc, self.mesh.gridEx, 'x') + # SRCy = src(tx.loc, self.mesh.gridEy, 'y') + # SRCz = src(tx.loc, self.mesh.gridEz, 'z') - elif tx.txType == 'VMD_B': - src = Sources.MagneticDipoleFields - SRCx = src(tx.loc, self.mesh.gridFx, 'x') - SRCy = src(tx.loc, self.mesh.gridFy, 'y') - SRCz = src(tx.loc, self.mesh.gridFz, 'z') + # elif tx.txType == 'VMD_B': + # src = Sources.MagneticDipoleFields + # SRCx = src(tx.loc, self.mesh.gridFx, 'x') + # SRCy = src(tx.loc, self.mesh.gridFy, 'y') + # SRCz = src(tx.loc, self.mesh.gridFz, 'z') - elif tx.txType == 'CircularLoop': - src = Sources.MagneticLoopVectorPotential - SRCx = src(tx.loc, self.mesh.gridEx, 'x', tx.radius) - SRCy = src(tx.loc, self.mesh.gridEy, 'y', tx.radius) - SRCz = src(tx.loc, self.mesh.gridEz, 'z', tx.radius) - else: + # elif tx.txType == 'CircularLoop': + # src = Sources.MagneticLoopVectorPotential + # SRCx = src(tx.loc, self.mesh.gridEx, 'x', tx.radius) + # SRCy = src(tx.loc, self.mesh.gridEy, 'y', tx.radius) + # SRCz = src(tx.loc, self.mesh.gridEz, 'z', tx.radius) + # else: - raise NotImplemented('%s txType is not implemented' % tx.txType) - SRC = np.concatenate((SRCx, SRCy, SRCz)) + # raise NotImplemented('%s txType is not implemented' % tx.txType) + # SRC = np.concatenate((SRCx, SRCy, SRCz)) - else: - raise Exception('Unknown mesh for VMD') + # else: + # raise Exception('Unknown mesh for VMD') - rhs[i] = SRC + # rhs[i] = SRC - mui = self.MfMui - if tx.txType == 'VMD_B': - b_0 = np.concatenate(rhs).reshape((self.mesh.nF, len(Txs)), order='F') - else: - a = np.concatenate(rhs).reshape((self.mesh.nE, len(Txs)), order='F') - C = self.mesh.edgeCurl - b_0 = C*a + # mui = self.MfMui + # if tx.txType == 'VMD_B': + # b_0 = np.concatenate(rhs).reshape((self.mesh.nF, len(Txs)), order='F') + # else: + # a = np.concatenate(rhs).reshape((self.mesh.nE, len(Txs)), order='F') + # C = self.mesh.edgeCurl + b_0 = getSource(self,freq) #C*a return -1j*omega(freq)*b_0 @@ -340,7 +439,6 @@ class ProblemFDEM_b(BaseFDEMProblem): raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) - ########################################################################################## ################################ H-J Formulation ######################################### ########################################################################################## @@ -417,36 +515,36 @@ class ProblemFDEM_j(BaseFDEMProblem): return C * ( MeMuI * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) ) ) - def getjs(self,freq): - """ - :param float freq: Frequency - :rtype: numpy.ndarray (nE, nTx) - :return: j_s - """ - Txs = self.survey.getTransmitters(freq) - rhs = range(len(Txs)) - for i, tx in enumerate(Txs): - if tx.txType == 'VMD': - src = Sources.MagneticDipoleVectorPotential - SRCx = src(tx.loc, self.mesh.gridFx, 'x') - SRCy = src(tx.loc, self.mesh.gridFy, 'y') - SRCz = src(tx.loc, self.mesh.gridFz, 'z') + # def getjs(self,freq): + # """ + # :param float freq: Frequency + # :rtype: numpy.ndarray (nE, nTx) + # :return: j_s + # """ + # Txs = self.survey.getTransmitters(freq) + # rhs = range(len(Txs)) + # for i, tx in enumerate(Txs): + # if tx.txType == 'VMD': + # src = Sources.MagneticDipoleVectorPotential + # SRCx = src(tx.loc, self.mesh.gridFx, 'x') + # SRCy = src(tx.loc, self.mesh.gridFy, 'y') + # SRCz = src(tx.loc, self.mesh.gridFz, 'z') - elif tx.txType == 'CircularLoop': - src = Sources.MagneticLoopVectorPotential - SRCx = src(tx.loc, self.mesh.gridFx, 'x', tx.radius) - SRCy = src(tx.loc, self.mesh.gridFy, 'y', tx.radius) - SRCz = src(tx.loc, self.mesh.gridFz, 'z', tx.radius) - else: - raise NotImplemented('%s txType is not implemented' % tx.txType) - rhs[i] = np.concatenate((SRCx, SRCy, SRCz)) + # elif tx.txType == 'CircularLoop': + # src = Sources.MagneticLoopVectorPotential + # SRCx = src(tx.loc, self.mesh.gridFx, 'x', tx.radius) + # SRCy = src(tx.loc, self.mesh.gridFy, 'y', tx.radius) + # SRCz = src(tx.loc, self.mesh.gridFz, 'z', tx.radius) + # else: + # raise NotImplemented('%s txType is not implemented' % tx.txType) + # rhs[i] = np.concatenate((SRCx, SRCy, SRCz)) - a = np.concatenate(rhs).reshape((self.mesh.nF, len(Txs)), order='F') - a = Utils.mkvc(a) - MeMuI = self.MeMuI - C = self.mesh.edgeCurl + # a = np.concatenate(rhs).reshape((self.mesh.nF, len(Txs)), order='F') + # a = Utils.mkvc(a) + # MeMuI = self.MeMuI + # C = self.mesh.edgeCurl - return C*MeMuI*C.T*a + # return C*MeMuI*C.T*a def getRHS(self, freq): """ @@ -454,7 +552,7 @@ class ProblemFDEM_j(BaseFDEMProblem): :rtype: numpy.ndarray (nE, nTx) :return: RHS """ - j_s = self.getjs(freq) + j_s = getSource(self,freq) return -1j*omega(freq)*j_s def calcFields(self, sol, freq, fieldType, adjoint=False): @@ -555,37 +653,37 @@ class ProblemFDEM_h(BaseFDEMProblem): return (C.T * (dMf_dsigi * (dsigi_dsig * (dsig_dm * v)))) - def getjs(self,freq): - """ - :param float freq: Frequency - :rtype: numpy.ndarray (nE, nTx) - :return: j_s - """ - Txs = self.survey.getTransmitters(freq) - rhs = range(len(Txs)) - for i, tx in enumerate(Txs): - if tx.txType == 'VMD': - src = Sources.MagneticDipoleVectorPotential - SRCx = src(tx.loc, self.mesh.gridFx, 'x') - SRCy = src(tx.loc, self.mesh.gridFy, 'y') - SRCz = src(tx.loc, self.mesh.gridFz, 'z') + # def getjs(self,freq): + # """ + # :param float freq: Frequency + # :rtype: numpy.ndarray (nE, nTx) + # :return: j_s + # """ + # Txs = self.survey.getTransmitters(freq) + # rhs = range(len(Txs)) + # for i, tx in enumerate(Txs): + # if tx.txType == 'VMD': + # src = Sources.MagneticDipoleVectorPotential + # SRCx = src(tx.loc, self.mesh.gridFx, 'x') + # SRCy = src(tx.loc, self.mesh.gridFy, 'y') + # SRCz = src(tx.loc, self.mesh.gridFz, 'z') - elif tx.txType == 'CircularLoop': - src = Sources.MagneticLoopVectorPotential - SRCx = src(tx.loc, self.mesh.gridFx, 'x', tx.radius) - SRCy = src(tx.loc, self.mesh.gridFy, 'y', tx.radius) - SRCz = src(tx.loc, self.mesh.gridFz, 'z', tx.radius) - else: - raise NotImplemented('%s txType is not implemented' % tx.txType) - rhs[i] = np.concatenate((SRCx, SRCy, SRCz)) + # elif tx.txType == 'CircularLoop': + # src = Sources.MagneticLoopVectorPotential + # SRCx = src(tx.loc, self.mesh.gridFx, 'x', tx.radius) + # SRCy = src(tx.loc, self.mesh.gridFy, 'y', tx.radius) + # SRCz = src(tx.loc, self.mesh.gridFz, 'z', tx.radius) + # else: + # raise NotImplemented('%s txType is not implemented' % tx.txType) + # rhs[i] = np.concatenate((SRCx, SRCy, SRCz)) - a = np.concatenate(rhs).reshape((self.mesh.nF, len(Txs)), order='F') - a = Utils.mkvc(a) + # a = np.concatenate(rhs).reshape((self.mesh.nF, len(Txs)), order='F') + # a = Utils.mkvc(a) - MeMuI = self.MeMuI - C = self.mesh.edgeCurl + # MeMuI = self.MeMuI + # C = self.mesh.edgeCurl - return MeMuI*C.T*a #C*MeMuI*C.T*a + # return MeMuI*C.T*a #C*MeMuI*C.T*a def getRHS(self, freq): """ @@ -593,29 +691,12 @@ class ProblemFDEM_h(BaseFDEMProblem): :rtype: numpy.ndarray (nE, nTx) :return: RHS """ - MeMu = self.MeMu - MfSigi = self.MfSigmai - C = self.mesh.edgeCurl - Hp = self.getjs(freq) - return -1j*omega(freq)*MeMu*Hp #C.T*MfSigi*j_s - - # def getRHSDeriv(self, freq, v, adjoint=False): - # """ - # :param float freq: Frequency - # :rtype: numpy.ndarray (nE, nTx) - # :return: RHSDeriv - # """ - # C = self.mesh.edgeCurl - # sig = self.curModel.transform - # sigi = 1/sig - # j_s = self.getjs(freq) - # dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(j_s) - # dsig_dm = self.curModel.transformDeriv - # dsigi_dsig = -Utils.sdiag(sigi)**2 # only works for diagonal matrices - - # if adjoint: - # return dsig_dm.T * dsigi_dsig.T * dMf_dsigi.T * C * v - # return C.T * dMf_dsigi * dsigi_dsig * dsig_dm * v + # MeMu = self.MeMu + # MfSigi = self.MfSigmai + # C = self.mesh.edgeCurl + # Hp = self.getjs(freq) + b_0 = getSource(self,freq) + return -1j*omega(freq)*b_0 #C.T*MfSigi*j_s def calcFields(self, sol, freq, fieldType, adjoint=False): h = sol @@ -630,25 +711,4 @@ class ProblemFDEM_h(BaseFDEMProblem): raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) def calcFieldsDeriv(self, sol, freq, fieldType, v, adjoint=False): - return None - # h = sol - # A = self.getA(freq) - - # if fieldType == 'j': - # C = self.mesh.edgeCurl - # j_s = self.getjs(freq) - # if adjoint: - # dh = self.calcFieldsDeriv(h,freq,'h',C.T*v,adjoint=True) - # return dh - # dh = self.calcFieldsDeriv(h,freq,'h',v) - # return C*dh - j_s - - # elif fieldType == 'h': - # if adjoint: - # ATinv = self.Solver(A.T, **self.solverOpts) - # ATinvv = ATinv*v - # return self.getRHSDeriv(freq,ATinvv,adjoint=True) - # dRHSh = self.getRHSDeriv(freq,v,adjoint) - # Ainv = self.Solver(A, **self.solverOpts) - # return Ainv*dRHSh - # raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) \ No newline at end of file + return None \ No newline at end of file diff --git a/simpegEM/Tests/test_FDEM.py b/simpegEM/Tests/test_FDEM.py index 4c2ecefd..3a1c4c92 100644 --- a/simpegEM/Tests/test_FDEM.py +++ b/simpegEM/Tests/test_FDEM.py @@ -9,7 +9,7 @@ FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order CONDUCTIVITY = 1e1 MU = mu_0 freq = 1e-1 -addrandoms = True # important to addrandoms if testing HJ formulation with VMD source! (or else jz ~ 0) +addrandoms = True def getProblem(fdemType, comp): cs = 5. From a697833b9f0c79ad5fab10f2ec44172386561153 Mon Sep 17 00:00:00 2001 From: Lindsey Date: Mon, 2 Mar 2015 14:49:06 -0800 Subject: [PATCH 22/25] changed todo for constant mu to anisotropic mu --- simpegEM/Base.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/simpegEM/Base.py b/simpegEM/Base.py index 459a4e1e..2b005540 100644 --- a/simpegEM/Base.py +++ b/simpegEM/Base.py @@ -50,14 +50,14 @@ class BaseEMProblem(Problem.BaseProblem): @property def MeMuI(self): - # Assuming isotropic mu! + # TODO: Assuming isotropic mu if getattr(self, '_MeMuI', None) is None: self._MeMuI = self.mesh.getEdgeInnerProduct(self.mu, invMat=True) return self._MeMuI @property def MeMu(self): - #TODO: assuming constant mu + #TODO: Assuming isotropic mu if getattr(self, '_MeMu', None) is None: self._MeMu = self.mesh.getEdgeInnerProduct(self.mu) return self._MeMu From ed61178b452dbfd3809dadd5ba8d529ffe61a142 Mon Sep 17 00:00:00 2001 From: Lindsey Date: Mon, 2 Mar 2015 15:07:50 -0800 Subject: [PATCH 23/25] removed comment blocks of old code from get RHS in each formulation --- simpegEM/FDEM/FDEM.py | 133 +----------------------------------------- 1 file changed, 1 insertion(+), 132 deletions(-) diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py index e816be19..81b0d3ab 100644 --- a/simpegEM/FDEM/FDEM.py +++ b/simpegEM/FDEM/FDEM.py @@ -258,27 +258,6 @@ class ProblemFDEM_e(BaseFDEMProblem): :rtype: numpy.ndarray (nE, nTx) :return: RHS """ - # Txs = self.survey.getTransmitters(freq) - # rhs = range(len(Txs)) - # for i, tx in enumerate(Txs): - # if tx.txType == 'VMD': - # src = Sources.MagneticDipoleVectorPotential - # SRCx = src(tx.loc, self.mesh.gridEx, 'x') - # SRCy = src(tx.loc, self.mesh.gridEy, 'y') - # SRCz = src(tx.loc, self.mesh.gridEz, 'z') - - # elif tx.txType == 'CircularLoop': - # src = Sources.MagneticLoopVectorPotential - # SRCx = src(tx.loc, self.mesh.gridEx, 'x', tx.radius) - # SRCy = src(tx.loc, self.mesh.gridEy, 'y', tx.radius) - # SRCz = src(tx.loc, self.mesh.gridEz, 'z', tx.radius) - # else: - # raise NotImplemented('%s txType is not implemented' % tx.txType) - # rhs[i] = np.concatenate((SRCx, SRCy, SRCz)) - - # a = np.concatenate(rhs).reshape((self.mesh.nE, len(Txs)), order='F') - # mui = self.MfMui - # C = self.mesh.edgeCurl j_s = getSource(self,freq) #C.T*mui*C*a return -1j*omega(freq)*j_s @@ -350,57 +329,8 @@ class ProblemFDEM_b(BaseFDEMProblem): :rtype: numpy.ndarray (nE, nTx) :return: RHS """ - # Txs = self.survey.getTransmitters(freq) - # rhs = range(len(Txs)) - # for i, tx in enumerate(Txs): - - # if self.mesh._meshType is 'CYL': - # if self.mesh.isSymmetric: - # if tx.txType == 'VMD': - # SRC = Sources.MagneticDipoleVectorPotential(tx.loc, self.mesh.gridEy, 'y') - # elif tx.txType =='CircularLoop': - # SRC = Sources.MagneticLoopVectorPotential(tx.loc, self.mesh.gridEy, 'y', tx.radius) - # else: - # raise NotImplementedError('Only VMD and CircularLoop') - # else: - # raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!') - # elif self.mesh._meshType is 'TENSOR': - - # if tx.txType == 'VMD': - # src = Sources.MagneticDipoleVectorPotential - # SRCx = src(tx.loc, self.mesh.gridEx, 'x') - # SRCy = src(tx.loc, self.mesh.gridEy, 'y') - # SRCz = src(tx.loc, self.mesh.gridEz, 'z') - - # elif tx.txType == 'VMD_B': - # src = Sources.MagneticDipoleFields - # SRCx = src(tx.loc, self.mesh.gridFx, 'x') - # SRCy = src(tx.loc, self.mesh.gridFy, 'y') - # SRCz = src(tx.loc, self.mesh.gridFz, 'z') - - # elif tx.txType == 'CircularLoop': - # src = Sources.MagneticLoopVectorPotential - # SRCx = src(tx.loc, self.mesh.gridEx, 'x', tx.radius) - # SRCy = src(tx.loc, self.mesh.gridEy, 'y', tx.radius) - # SRCz = src(tx.loc, self.mesh.gridEz, 'z', tx.radius) - # else: - - # raise NotImplemented('%s txType is not implemented' % tx.txType) - # SRC = np.concatenate((SRCx, SRCy, SRCz)) - - # else: - # raise Exception('Unknown mesh for VMD') - - # rhs[i] = SRC - - # mui = self.MfMui - # if tx.txType == 'VMD_B': - # b_0 = np.concatenate(rhs).reshape((self.mesh.nF, len(Txs)), order='F') - # else: - # a = np.concatenate(rhs).reshape((self.mesh.nE, len(Txs)), order='F') - # C = self.mesh.edgeCurl - b_0 = getSource(self,freq) #C*a + b_0 = getSource(self,freq) return -1j*omega(freq)*b_0 @@ -515,36 +445,6 @@ class ProblemFDEM_j(BaseFDEMProblem): return C * ( MeMuI * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) ) ) - # def getjs(self,freq): - # """ - # :param float freq: Frequency - # :rtype: numpy.ndarray (nE, nTx) - # :return: j_s - # """ - # Txs = self.survey.getTransmitters(freq) - # rhs = range(len(Txs)) - # for i, tx in enumerate(Txs): - # if tx.txType == 'VMD': - # src = Sources.MagneticDipoleVectorPotential - # SRCx = src(tx.loc, self.mesh.gridFx, 'x') - # SRCy = src(tx.loc, self.mesh.gridFy, 'y') - # SRCz = src(tx.loc, self.mesh.gridFz, 'z') - - # elif tx.txType == 'CircularLoop': - # src = Sources.MagneticLoopVectorPotential - # SRCx = src(tx.loc, self.mesh.gridFx, 'x', tx.radius) - # SRCy = src(tx.loc, self.mesh.gridFy, 'y', tx.radius) - # SRCz = src(tx.loc, self.mesh.gridFz, 'z', tx.radius) - # else: - # raise NotImplemented('%s txType is not implemented' % tx.txType) - # rhs[i] = np.concatenate((SRCx, SRCy, SRCz)) - - # a = np.concatenate(rhs).reshape((self.mesh.nF, len(Txs)), order='F') - # a = Utils.mkvc(a) - # MeMuI = self.MeMuI - # C = self.mesh.edgeCurl - - # return C*MeMuI*C.T*a def getRHS(self, freq): """ @@ -653,37 +553,6 @@ class ProblemFDEM_h(BaseFDEMProblem): return (C.T * (dMf_dsigi * (dsigi_dsig * (dsig_dm * v)))) - # def getjs(self,freq): - # """ - # :param float freq: Frequency - # :rtype: numpy.ndarray (nE, nTx) - # :return: j_s - # """ - # Txs = self.survey.getTransmitters(freq) - # rhs = range(len(Txs)) - # for i, tx in enumerate(Txs): - # if tx.txType == 'VMD': - # src = Sources.MagneticDipoleVectorPotential - # SRCx = src(tx.loc, self.mesh.gridFx, 'x') - # SRCy = src(tx.loc, self.mesh.gridFy, 'y') - # SRCz = src(tx.loc, self.mesh.gridFz, 'z') - - # elif tx.txType == 'CircularLoop': - # src = Sources.MagneticLoopVectorPotential - # SRCx = src(tx.loc, self.mesh.gridFx, 'x', tx.radius) - # SRCy = src(tx.loc, self.mesh.gridFy, 'y', tx.radius) - # SRCz = src(tx.loc, self.mesh.gridFz, 'z', tx.radius) - # else: - # raise NotImplemented('%s txType is not implemented' % tx.txType) - # rhs[i] = np.concatenate((SRCx, SRCy, SRCz)) - - # a = np.concatenate(rhs).reshape((self.mesh.nF, len(Txs)), order='F') - # a = Utils.mkvc(a) - - # MeMuI = self.MeMuI - # C = self.mesh.edgeCurl - - # return MeMuI*C.T*a #C*MeMuI*C.T*a def getRHS(self, freq): """ From c480279ae510d994500329da4b6977f02d51b897 Mon Sep 17 00:00:00 2001 From: Lindsey Date: Mon, 2 Mar 2015 15:43:17 -0800 Subject: [PATCH 24/25] added _makeASymmetric option to j, b implementations --- simpegEM/Base.py | 8 ++++++++ simpegEM/FDEM/FDEM.py | 28 +++++++++++++++++++++++----- 2 files changed, 31 insertions(+), 5 deletions(-) diff --git a/simpegEM/Base.py b/simpegEM/Base.py index 2b005540..d1e6bb08 100644 --- a/simpegEM/Base.py +++ b/simpegEM/Base.py @@ -17,6 +17,14 @@ class BaseEMProblem(Problem.BaseProblem): verbose = False + #################################################### + # Make A Symmetric + #################################################### + @property + def _makeASymmetric(self): + if getattr(self, '__makeASymmetric', None) is None: + self.__makeASymmetric = True + return self.__makeASymmetric #################################################### # Mu Model diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py index 81b0d3ab..9943665d 100644 --- a/simpegEM/FDEM/FDEM.py +++ b/simpegEM/FDEM/FDEM.py @@ -242,6 +242,7 @@ class ProblemFDEM_e(BaseFDEMProblem): return C.T*mui*C + 1j*omega(freq)*sig + def getADeriv(self, freq, u, v, adjoint=False): sig = self.curModel.transform dsig_dm = self.curModel.transformDeriv @@ -259,7 +260,7 @@ class ProblemFDEM_e(BaseFDEMProblem): :return: RHS """ - j_s = getSource(self,freq) #C.T*mui*C*a + j_s = getSource(self,freq) return -1j*omega(freq)*j_s def calcFields(self, sol, freq, fieldType, adjoint=False): @@ -303,7 +304,11 @@ class ProblemFDEM_b(BaseFDEMProblem): C = self.mesh.edgeCurl iomega = 1j * omega(freq) * sp.eye(self.mesh.nF) - return C*sigI*C.T*mui + iomega + A = C*sigI*C.T*mui + iomega + + if self._makeASymmetric is True: + return mui.T*A + return A def getADeriv(self, freq, u, v, adjoint=False): @@ -332,7 +337,11 @@ class ProblemFDEM_b(BaseFDEMProblem): b_0 = getSource(self,freq) - return -1j*omega(freq)*b_0 + rhs = -1j*omega(freq)*b_0 + if self._makeASymmetric is True: + mui = self.MfMui + return mui.T*rhs + return rhs def calcFields(self, sol, freq, fieldType, adjoint=False): b = sol @@ -421,7 +430,11 @@ class ProblemFDEM_j(BaseFDEMProblem): C = self.mesh.edgeCurl iomega = 1j * omega(freq) * sp.eye(self.mesh.nF) - return C * MeMuI * C.T * MfSigi + iomega + A = C * MeMuI * C.T * MfSigi + iomega + + if self._makeASymmetric is True: + return MfSigi.T*A + return A def getADeriv(self, freq, u, v, adjoint=False): @@ -453,7 +466,12 @@ class ProblemFDEM_j(BaseFDEMProblem): :return: RHS """ j_s = getSource(self,freq) - return -1j*omega(freq)*j_s + rhs = -1j*omega(freq)*j_s + + if self._makeASymmetric is True: + MfSigi = self.MfSigmai + return MfSigi.T*rhs + return rhs def calcFields(self, sol, freq, fieldType, adjoint=False): j = sol From 65842379b5d7d67632453019c886609d0fc86484 Mon Sep 17 00:00:00 2001 From: Lindsey Date: Tue, 3 Mar 2015 10:02:19 -0800 Subject: [PATCH 25/25] Added self._makeASymmetric in b and j formulations (e and h are symmetric already) --- simpegEM/FDEM/FDEM.py | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py index 9943665d..1a45d139 100644 --- a/simpegEM/FDEM/FDEM.py +++ b/simpegEM/FDEM/FDEM.py @@ -323,8 +323,12 @@ class ProblemFDEM_b(BaseFDEMProblem): dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig)(vec) if adjoint: + if self._makeASymmetric is True: + v = mui * v return dsig_dm.T * ( dMe_dsig.T * ( dMeSigmaI_dI.T * ( C.T * v ) ) ) + if self._makeASymmetric is True: + return mui.T * ( C * ( dMeSigmaI_dI * ( dMe_dsig * ( dsig_dm * v ) ) ) ) return C * ( dMeSigmaI_dI * ( dMe_dsig * ( dsig_dm * v ) ) ) @@ -446,6 +450,7 @@ class ProblemFDEM_j(BaseFDEMProblem): """ MeMuI = self.MeMuI + MfSigi = self.MfSigmai C = self.mesh.edgeCurl sig = self.curModel.transform sigi = 1/sig @@ -454,8 +459,12 @@ class ProblemFDEM_j(BaseFDEMProblem): dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(u) if adjoint: + if self._makeASymmetric is True: + v = MfSigi * v return dsig_dm.T * ( dsigi_dsig.T *( dMf_dsigi.T * ( C * ( MeMuI.T * ( C.T * v ) ) ) ) ) + if self._makeASymmetric is True: + return MfSigi.T * ( C * ( MeMuI * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) ) ) ) return C * ( MeMuI * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) ) ) @@ -578,21 +587,16 @@ class ProblemFDEM_h(BaseFDEMProblem): :rtype: numpy.ndarray (nE, nTx) :return: RHS """ - # MeMu = self.MeMu - # MfSigi = self.MfSigmai - # C = self.mesh.edgeCurl - # Hp = self.getjs(freq) b_0 = getSource(self,freq) - return -1j*omega(freq)*b_0 #C.T*MfSigi*j_s + return -1j*omega(freq)*b_0 def calcFields(self, sol, freq, fieldType, adjoint=False): h = sol if fieldType == 'j': C = self.mesh.edgeCurl - # j_s = self.getjs(freq) if adjoint: return C.T*h - return C*h #- j_s + return C*h elif fieldType == 'h': return h raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)