From 354e57f24ed9ef4e571b2522f50deba60503d238 Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Mon, 2 May 2016 16:47:16 -0700 Subject: [PATCH] Problem3D_CC and _N for IP are all tested a) fwd b) jvec, jtvec c) adjoint --- tests/em/static/test_DC_2D_analytic.py | 2 +- tests/em/static/test_IP_fwd.py | 85 +++++++++------- tests/em/static/test_IP_jvecjtvecadj.py | 126 ++++++++++++++++++++++++ 3 files changed, 179 insertions(+), 34 deletions(-) create mode 100644 tests/em/static/test_IP_jvecjtvecadj.py diff --git a/tests/em/static/test_DC_2D_analytic.py b/tests/em/static/test_DC_2D_analytic.py index a6e3e6ab..7a8fe1b6 100644 --- a/tests/em/static/test_DC_2D_analytic.py +++ b/tests/em/static/test_DC_2D_analytic.py @@ -22,7 +22,7 @@ class DCProblemAnalyticTests(unittest.TestCase): rx = DC.Rx.Dipole(M, N) src0 = DC.Src.Pole([rx], A0loc) - survey = DC.Survey([src0]) + survey = DC.Survey_ky([src0]) self.survey = survey self.mesh = mesh diff --git a/tests/em/static/test_IP_fwd.py b/tests/em/static/test_IP_fwd.py index 9643a0bb..98bb7e90 100644 --- a/tests/em/static/test_IP_fwd.py +++ b/tests/em/static/test_IP_fwd.py @@ -1,41 +1,43 @@ import unittest from SimPEG import Mesh, Utils, EM, Maps, np import SimPEG.EM.Static.DC as DC +import SimPEG.EM.Static.IP as IP class IPProblemAnalyticTests(unittest.TestCase): def setUp(self): cs = 12.5 - hx = [(cs,2, -1.3),(cs,61),(cs,2, 1.3)] - hy = [(cs,2, -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)]] + npad=2 + hx = [(cs,npad, -1.3),(cs,21),(cs,npad, 1.3)] + hy = [(cs,npad, -1.3),(cs,21),(cs,npad, 1.3)] + hz = [(cs,npad, -1.3),(cs,20)] + mesh = Mesh.TensorMesh([hx, hy, hz],x0="CCN") + x = mesh.vectorCCx[(mesh.vectorCCx>-80.)&(mesh.vectorCCx<80.)] + y = mesh.vectorCCx[(mesh.vectorCCy>-80.)&(mesh.vectorCCy<80.)] + Aloc = np.r_[-100., 0., 0.] + Bloc = np.r_[100., 0., 0.] + M = Utils.ndgrid(x-12.5,y, np.r_[0.]) + N = Utils.ndgrid(x+12.5,y, np.r_[0.]) + radius = 50. + xc = np.r_[0., 0., -100] blkind = Utils.ModelBuilder.getIndicesSphere(xc, radius, mesh.gridCC) sigmaInf = np.ones(mesh.nC)*1e-2 eta = np.zeros(mesh.nC) eta[blkind] = 0.1 sigma0 = sigmaInf*(1.-eta) - rx = DC.Rx.Dipole_ky(M, N) - src0 = DC.Src.Pole([rx], A0loc) - surveyDC = DC.Survey([src0]) - surveyIP = DC.Survey_ky([src0]) + rx = DC.Rx.Dipole(M, N) + src = DC.Src.Dipole([rx], Aloc, Bloc) + surveyDC = DC.Survey([src]) self.surveyDC = surveyDC - self.surveyIP = surveyIP - self.mesh = mesh - self.sigma = sigma - self.data_anal = data_anal + self.sigmaInf = sigmaInf + self.sigma0 = sigma0 + self.src = src + self.eta = eta try: from pymatsolver import MumpsSolver @@ -45,31 +47,48 @@ class IPProblemAnalyticTests(unittest.TestCase): def test_Problem3D_N(self): - 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)/self.data_anal)**2 / self.data_anal.size + problemDC = DC.Problem3D_N(self.mesh) + problemDC.Solver = self.Solver + problemDC.pair(self.surveyDC) + data0 = self.surveyDC.dpred(self.sigma0) + finf = problemDC.fields(self.sigmaInf) + datainf = self.surveyDC.dpred(self.sigmaInf, f=finf) + problemIP = IP.Problem3D_N(self.mesh, sigma=self.sigmaInf, Ainv=problemDC.Ainv, f=finf) + problemIP.Solver = self.Solver + surveyIP = IP.Survey([self.src]) + problemIP.pair(surveyIP) + data_full = data0 - datainf + data = surveyIP.dpred(self.eta) + err= np.linalg.norm((data-data_full)/data_full)**2 / data_full.size if err < 0.05: passed = True - print ">> DC analytic test for Problem3D_N is passed" + print ">> IP forward test for Problem3D_N is passed" else: passed = False - print ">> DC analytic test for Problem3D_N is failed" + print ">> IP forward test for Problem3D_N is failed" self.assertTrue(passed) def test_Problem3D_CC(self): - 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)/self.data_anal)**2 / self.data_anal.size + + problemDC = DC.Problem3D_CC(self.mesh) + problemDC.Solver = self.Solver + problemDC.pair(self.surveyDC) + data0 = self.surveyDC.dpred(self.sigma0) + finf = problemDC.fields(self.sigmaInf) + datainf = self.surveyDC.dpred(self.sigmaInf, f=finf) + problemIP = IP.Problem3D_CC(self.mesh, rho=1./self.sigmaInf, Ainv=problemDC.Ainv, f=finf) + problemIP.Solver = self.Solver + surveyIP = IP.Survey([self.src]) + problemIP.pair(surveyIP) + data_full = data0 - datainf + data = surveyIP.dpred(self.eta) + err= np.linalg.norm((data-data_full)/data_full)**2 / data_full.size if err < 0.05: passed = True - print ">> DC analytic test for Problem3D_CC is passed" + print ">> IP forward test for Problem3D_CC is passed" else: passed = False - print ">> DC analytic test for Problem3D_CC is failed" + print ">> IP forward test for Problem3D_CC is failed" self.assertTrue(passed) if __name__ == '__main__': diff --git a/tests/em/static/test_IP_jvecjtvecadj.py b/tests/em/static/test_IP_jvecjtvecadj.py new file mode 100644 index 00000000..7a455784 --- /dev/null +++ b/tests/em/static/test_IP_jvecjtvecadj.py @@ -0,0 +1,126 @@ +import unittest +from SimPEG import * +import SimPEG.EM.Static.DC as DC +import SimPEG.EM.Static.IP as IP + + +class IPProblemTestsCC(unittest.TestCase): + + def setUp(self): + + aSpacing=2.5 + nElecs=5 + + surveySize = nElecs*aSpacing - aSpacing + cs = surveySize/nElecs/4 + + mesh = Mesh.TensorMesh([ + [(cs,10, -1.3),(cs,surveySize/cs),(cs,10, 1.3)], + [(cs,3, -1.3),(cs,3,1.3)], + # [(cs,5, -1.3),(cs,10)] + ],'CN') + + srcList = DC.Utils.WennerSrcList(nElecs, aSpacing, in2D=True) + survey = IP.Survey(srcList) + sigma = np.ones(mesh.nC) + problem = IP.Problem3D_CC(mesh, rho=1./sigma) + problem.pair(survey) + mSynth = np.ones(mesh.nC)*0.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=1e4) + 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 + + 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_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 IPProblemTestsN(unittest.TestCase): + + def setUp(self): + + aSpacing=2.5 + nElecs=5 + + surveySize = nElecs*aSpacing - aSpacing + cs = surveySize/nElecs/4 + + mesh = Mesh.TensorMesh([ + [(cs,10, -1.3),(cs,surveySize/cs),(cs,10, 1.3)], + [(cs,3, -1.3),(cs,3,1.3)], + # [(cs,5, -1.3),(cs,10)] + ],'CN') + + srcList = DC.Utils.WennerSrcList(nElecs, aSpacing, in2D=True) + survey = IP.Survey(srcList) + sigma = np.ones(mesh.nC) + problem = IP.Problem3D_N(mesh, sigma=sigma) + problem.pair(survey) + mSynth = np.ones(mesh.nC)*0.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=1e4) + 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 + + 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-8 + 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) + +if __name__ == '__main__': + unittest.main()