From 57380ea2d8ed53ef9929844adf66b79a447e3523 Mon Sep 17 00:00:00 2001 From: Lindsey Date: Wed, 25 Feb 2015 16:28:54 -0800 Subject: [PATCH] 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__':