From 350818d80256c4dfc6a363ae6645288a2e7819b7 Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Sun, 1 May 2016 13:21:48 -0700 Subject: [PATCH] add analytic test for 2D dc problems. --- tests/em/static/test_DC_2D_analytic.py | 71 ++++++++++++++ tests/em/static/test_DC_2D_jvecjtvecadj.py | 102 ++++++++++----------- tests/em/static/test_DC_analytic.py | 48 +++++----- 3 files changed, 145 insertions(+), 76 deletions(-) create mode 100644 tests/em/static/test_DC_2D_analytic.py diff --git a/tests/em/static/test_DC_2D_analytic.py b/tests/em/static/test_DC_2D_analytic.py new file mode 100644 index 00000000..53d494e2 --- /dev/null +++ b/tests/em/static/test_DC_2D_analytic.py @@ -0,0 +1,71 @@ +import unittest +from SimPEG import Mesh, Utils, EM, Maps, np +import SimPEG.EM.Static.DC as DC + +class DCProblemAnalyticTests(unittest.TestCase): + + def setUp(self): + + cs = 25. + hx = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)] + hy = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)] + hz = [(cs,7, -1.3),(cs,20)] + mesh = Mesh.TensorMesh([hx, hy, hz],x0="CCN") + sigma = np.ones(mesh.nC)*1e-2 + + x = mesh.vectorCCx[(mesh.vectorCCx>-155.)&(mesh.vectorCCx<155.)] + y = mesh.vectorCCx[(mesh.vectorCCy>-155.)&(mesh.vectorCCy<155.)] + Aloc = np.r_[-200., 0., 0.] + Bloc = np.r_[200., 0., 0.] + M = Utils.ndgrid(x-25.,y, np.r_[0.]) + N = Utils.ndgrid(x+25.,y, np.r_[0.]) + phiA = EM.Analytics.DCAnalyticHalf(Aloc, [M,N], 1e-2, flag="halfspace") + phiB = EM.Analytics.DCAnalyticHalf(Bloc, [M,N], 1e-2, flag="halfspace") + data_anal = phiA-phiB + + rx = DC.Rx.Dipole(M, N) + src = DC.Src.Dipole([rx], Aloc, Bloc) + survey = DC.Survey([src]) + + self.survey = survey + self.mesh = mesh + self.sigma = sigma + self.data_anal = data_anal + + try: + from pymatsolver import MumpsSolver + self.Solver = MumpsSolver + except ImportError, e: + self.Solver = SolverLU + + def test_Problem3D_N(self): + problem = DC.Problem3D_N(self.mesh) + problem.Solver = self.Solver + problem.pair(self.survey) + data = self.survey.dpred(self.sigma) + err= np.linalg.norm(data-self.data_anal)/np.linalg.norm(self.data_anal) + if err < 0.2: + passed = True + print ">> DC analytic test for Problem3D_N is passed" + else: + passed = False + print ">> DC analytic test for Problem3D_N is failed" + self.assertTrue(passed) + + def test_Problem3D_CC(self): + problem = DC.Problem3D_CC(self.mesh) + problem.Solver = self.Solver + problem.pair(self.survey) + data = self.survey.dpred(self.sigma) + err= np.linalg.norm(data-self.data_anal)/np.linalg.norm(self.data_anal) + if err < 0.2: + passed = True + print ">> DC analytic test for Problem3D_CC is passed" + else: + passed = False + print ">> DC analytic test for Problem3D_CC is failed" + self.assertTrue(passed) + +if __name__ == '__main__': + unittest.main() + diff --git a/tests/em/static/test_DC_2D_jvecjtvecadj.py b/tests/em/static/test_DC_2D_jvecjtvecadj.py index 9412614e..b96529cf 100644 --- a/tests/em/static/test_DC_2D_jvecjtvecadj.py +++ b/tests/em/static/test_DC_2D_jvecjtvecadj.py @@ -3,65 +3,65 @@ from SimPEG import * import SimPEG.EM.Static.DC as DC -# class DCProblem_2DTestsCC(unittest.TestCase): +class DCProblem_2DTestsCC(unittest.TestCase): -# def setUp(self): + def setUp(self): -# cs = 12.5 -# hx = [(cs,7, -1.3),(cs,61),(cs,7, 1.3)] -# hy = [(cs,7, -1.3),(cs,20)] -# mesh = Mesh.TensorMesh([hx, hy],x0="CN") -# x = np.linspace(-135, 250., 20) -# M = Utils.ndgrid(x-12.5, np.r_[0.]) -# N = Utils.ndgrid(x+12.5, np.r_[0.]) -# A0loc = np.r_[-150, 0.] -# A1loc = np.r_[-130, 0.] -# rxloc = [np.c_[M, np.zeros(20)], np.c_[N, np.zeros(20)]] -# rx = DC.Rx.Dipole_ky(M, N) -# src0 = DC.Src.Pole([rx], A0loc) -# src1 = DC.Src.Pole([rx], A1loc) -# survey = DC.Survey_ky([src0, src1]) -# problem = DC.Problem2D_CC(mesh, mapping=[('rho', Maps.IdentityMap(mesh))]) -# problem.pair(survey) + cs = 12.5 + hx = [(cs,7, -1.3),(cs,61),(cs,7, 1.3)] + hy = [(cs,7, -1.3),(cs,20)] + mesh = Mesh.TensorMesh([hx, hy],x0="CN") + x = np.linspace(-135, 250., 20) + M = Utils.ndgrid(x-12.5, np.r_[0.]) + N = Utils.ndgrid(x+12.5, np.r_[0.]) + A0loc = np.r_[-150, 0.] + A1loc = np.r_[-130, 0.] + rxloc = [np.c_[M, np.zeros(20)], np.c_[N, np.zeros(20)]] + rx = DC.Rx.Dipole_ky(M, N) + src0 = DC.Src.Pole([rx], A0loc) + src1 = DC.Src.Pole([rx], A1loc) + survey = DC.Survey_ky([src0, src1]) + problem = DC.Problem2D_CC(mesh, mapping=[('rho', Maps.IdentityMap(mesh))]) + problem.pair(survey) -# mSynth = np.ones(mesh.nC)*1. -# survey.makeSyntheticData(mSynth) + mSynth = np.ones(mesh.nC)*1. + survey.makeSyntheticData(mSynth) -# # Now set up the problem to do some minimization -# dmis = DataMisfit.l2_DataMisfit(survey) -# reg = Regularization.Tikhonov(mesh) -# opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6) -# invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e0) -# inv = Inversion.BaseInversion(invProb) + # Now set up the problem to do some minimization + dmis = DataMisfit.l2_DataMisfit(survey) + reg = Regularization.Tikhonov(mesh) + opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6) + invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e0) + inv = Inversion.BaseInversion(invProb) -# self.inv = inv -# self.reg = reg -# self.p = problem -# self.mesh = mesh -# self.m0 = mSynth -# self.survey = survey -# self.dmis = dmis + self.inv = inv + self.reg = reg + self.p = problem + self.mesh = mesh + self.m0 = mSynth + self.survey = survey + self.dmis = dmis -# def test_misfit(self): -# derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)] -# passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3) -# self.assertTrue(passed) + def test_misfit(self): + derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)] + passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3) + self.assertTrue(passed) -# def test_adjoint(self): -# # Adjoint Test -# u = np.random.rand(self.mesh.nC*self.survey.nSrc) -# v = np.random.rand(self.mesh.nC) -# w = np.random.rand(self.survey.dobs.shape[0]) -# wtJv = w.dot(self.p.Jvec(self.m0, v)) -# vtJtw = v.dot(self.p.Jtvec(self.m0, w)) -# passed = np.abs(wtJv - vtJtw) < 1e-10 -# print 'Adjoint Test', np.abs(wtJv - vtJtw), passed -# self.assertTrue(passed) + def test_adjoint(self): + # Adjoint Test + u = np.random.rand(self.mesh.nC*self.survey.nSrc) + v = np.random.rand(self.mesh.nC) + w = np.random.rand(self.survey.dobs.shape[0]) + wtJv = w.dot(self.p.Jvec(self.m0, v)) + vtJtw = v.dot(self.p.Jtvec(self.m0, w)) + passed = np.abs(wtJv - vtJtw) < 1e-10 + print 'Adjoint Test', np.abs(wtJv - vtJtw), passed + self.assertTrue(passed) -# def test_dataObj(self): -# derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)] -# passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3) -# self.assertTrue(passed) + def test_dataObj(self): + derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)] + passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3) + self.assertTrue(passed) class DCProblemTestsN(unittest.TestCase): diff --git a/tests/em/static/test_DC_analytic.py b/tests/em/static/test_DC_analytic.py index 53d494e2..e0b8b611 100644 --- a/tests/em/static/test_DC_analytic.py +++ b/tests/em/static/test_DC_analytic.py @@ -6,26 +6,23 @@ class DCProblemAnalyticTests(unittest.TestCase): def setUp(self): - cs = 25. - hx = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)] - hy = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)] - hz = [(cs,7, -1.3),(cs,20)] - mesh = Mesh.TensorMesh([hx, hy, hz],x0="CCN") - sigma = np.ones(mesh.nC)*1e-2 + cs = 12.5 + hx = [(cs,7, -1.3),(cs,61),(cs,7, 1.3)] + hy = [(cs,7, -1.3),(cs,20)] + mesh = Mesh.TensorMesh([hx, hy],x0="CN") + sighalf = 1e-2 + sigma = np.ones(mesh.nC)*sighalf + x = np.linspace(-135, 250., 20) + M = Utils.ndgrid(x-12.5, np.r_[0.]) + N = Utils.ndgrid(x+12.5, np.r_[0.]) + A0loc = np.r_[-150, 0.] + A1loc = np.r_[-130, 0.] + rxloc = [np.c_[M, np.zeros(20)], np.c_[N, np.zeros(20)]] + data_anal = EM.Analytics.DCAnalyticHalf(np.r_[A0loc, 0.], rxloc, sighalf, flag="halfspace") - x = mesh.vectorCCx[(mesh.vectorCCx>-155.)&(mesh.vectorCCx<155.)] - y = mesh.vectorCCx[(mesh.vectorCCy>-155.)&(mesh.vectorCCy<155.)] - Aloc = np.r_[-200., 0., 0.] - Bloc = np.r_[200., 0., 0.] - M = Utils.ndgrid(x-25.,y, np.r_[0.]) - N = Utils.ndgrid(x+25.,y, np.r_[0.]) - phiA = EM.Analytics.DCAnalyticHalf(Aloc, [M,N], 1e-2, flag="halfspace") - phiB = EM.Analytics.DCAnalyticHalf(Bloc, [M,N], 1e-2, flag="halfspace") - data_anal = phiA-phiB - - rx = DC.Rx.Dipole(M, N) - src = DC.Src.Dipole([rx], Aloc, Bloc) - survey = DC.Survey([src]) + rx = DC.Rx.Dipole_ky(M, N) + src0 = DC.Src.Pole([rx], A0loc) + survey = DC.Survey_ky([src0]) self.survey = survey self.mesh = mesh @@ -39,12 +36,13 @@ class DCProblemAnalyticTests(unittest.TestCase): self.Solver = SolverLU def test_Problem3D_N(self): - problem = DC.Problem3D_N(self.mesh) + + problem = DC.Problem2D_N(self.mesh) problem.Solver = self.Solver problem.pair(self.survey) data = self.survey.dpred(self.sigma) - err= np.linalg.norm(data-self.data_anal)/np.linalg.norm(self.data_anal) - if err < 0.2: + err= np.linalg.norm((data-self.data_anal)/self.data_anal)**2 / self.data_anal.size + if err < 0.05: passed = True print ">> DC analytic test for Problem3D_N is passed" else: @@ -53,12 +51,12 @@ class DCProblemAnalyticTests(unittest.TestCase): self.assertTrue(passed) def test_Problem3D_CC(self): - problem = DC.Problem3D_CC(self.mesh) + problem = DC.Problem2D_CC(self.mesh) problem.Solver = self.Solver problem.pair(self.survey) data = self.survey.dpred(self.sigma) - err= np.linalg.norm(data-self.data_anal)/np.linalg.norm(self.data_anal) - if err < 0.2: + err= np.linalg.norm((data-self.data_anal)/self.data_anal)**2 / self.data_anal.size + if err < 0.05: passed = True print ">> DC analytic test for Problem3D_CC is passed" else: