diff --git a/SimPEG/EM/Static/IP/ProblemIP.py b/SimPEG/EM/Static/IP/ProblemIP.py index d35bf09e..dcf12d6b 100644 --- a/SimPEG/EM/Static/IP/ProblemIP.py +++ b/SimPEG/EM/Static/IP/ProblemIP.py @@ -13,8 +13,6 @@ class IPPropMap(Maps.PropMap): (\\(\\eta\\)) is the default inversion property """ eta = Maps.Property("Electrical Chargeability", defaultInvProp = True) - # sigma = Maps.Property("Electrical Conductivity", defaultVal=mu_0, propertyLink=('rho',Maps.ReciprocalMap)) - # rho = Maps.Property("Electrical Resistivity", propertyLink=('sigma', Maps.ReciprocalMap)) class BaseIPProblem(BaseEMProblem): @@ -92,9 +90,10 @@ class BaseIPProblem(BaseEMProblem): dRHS_dmT = self.getRHSDeriv(src, ATinvdf_duT, adjoint=True) du_dmT = -dA_dmT + dRHS_dmT Jtv += df_dmT + du_dmT - + # Conductivity ((d u / d log sigma).T) if self._formulation is 'EB': return -Utils.mkvc(Jtv) + # Conductivity ((d u / d log rho).T) if self._formulation is 'HJ': return Utils.mkvc(Jtv) diff --git a/tests/em/static/test_DC_2D_analytic.py b/tests/em/static/test_DC_2D_analytic.py index 53d494e2..a6e3e6ab 100644 --- a/tests/em/static/test_DC_2D_analytic.py +++ b/tests/em/static/test_DC_2D_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 - - 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 + 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") rx = DC.Rx.Dipole(M, N) - src = DC.Src.Dipole([rx], Aloc, Bloc) - survey = DC.Survey([src]) + src0 = DC.Src.Pole([rx], A0loc) + survey = DC.Survey([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: diff --git a/tests/em/static/test_DC_analytic.py b/tests/em/static/test_DC_analytic.py index e0b8b611..53d494e2 100644 --- a/tests/em/static/test_DC_analytic.py +++ b/tests/em/static/test_DC_analytic.py @@ -6,23 +6,26 @@ class DCProblemAnalyticTests(unittest.TestCase): 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") - 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") + 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 - rx = DC.Rx.Dipole_ky(M, N) - src0 = DC.Src.Pole([rx], A0loc) - survey = DC.Survey_ky([src0]) + 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 @@ -36,13 +39,12 @@ class DCProblemAnalyticTests(unittest.TestCase): self.Solver = SolverLU def test_Problem3D_N(self): - - problem = DC.Problem2D_N(self.mesh) + 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)/self.data_anal)**2 / self.data_anal.size - if err < 0.05: + 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: @@ -51,12 +53,12 @@ class DCProblemAnalyticTests(unittest.TestCase): self.assertTrue(passed) def test_Problem3D_CC(self): - problem = DC.Problem2D_CC(self.mesh) + 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)/self.data_anal)**2 / self.data_anal.size - if err < 0.05: + 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: diff --git a/tests/em/static/test_IP_fwd.py b/tests/em/static/test_IP_fwd.py new file mode 100644 index 00000000..9643a0bb --- /dev/null +++ b/tests/em/static/test_IP_fwd.py @@ -0,0 +1,77 @@ +import unittest +from SimPEG import Mesh, Utils, EM, Maps, np +import SimPEG.EM.Static.DC as DC + +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)]] + + 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]) + + self.surveyDC = surveyDC + self.surveyIP = surveyIP + + 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.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 + if err < 0.05: + 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.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 + if err < 0.05: + 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() +