wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-08-01 12:50:05 +08:00
Code Issues Packages Projects Releases Wiki Activity

add analytic test for 2D dc problems.

Browse Source
This commit is contained in:
seogi_macbook
2016-05-01 13:21:48 -07:00
parent 4df6f340d3
commit 350818d802
3 changed files with 145 additions and 76 deletions
Download Patch File Download Diff File Expand all files Collapse all files
tests/em/static/test_DC_2D_analytic.py
+71
View File
@@ -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()
tests/em/static/test_DC_2D_jvecjtvecadj.py
+51 -51
View File
@@ -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):
tests/em/static/test_DC_analytic.py
+23 -25
View File
@@ -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: