TDEM tests running and passing individually

This commit is contained in:
Lindsey Heagy
2015-07-01 00:21:54 -07:00
parent 265a1189e8
commit 4fa5e9426d
9 changed files with 66 additions and 57 deletions
+1 -1
View File
@@ -36,7 +36,7 @@ if plotIt:
rxOffset=1e-3
rx = EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 30]]), np.logspace(-5,-3, 31), 'bz')
src = EM.TDEM.SrcTDEM(np.array([0., 0., 80]), 'VMD_MVP', [rx])
src = EM.TDEM.SrcTDEM_VMD_MVP([rx], np.array([0., 0., 80]))
survey = EM.TDEM.SurveyTDEM([src])
prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping)
+4 -4
View File
@@ -14,20 +14,20 @@ class FieldsTDEM(Problem.TimeFields):
def tovec(self):
nSrc, nF, nE = self.survey.nSrc, self.mesh.nF, self.mesh.nE
u = np.empty(0 if nSrc == 1 else (0, nSrc))
u = np.empty((0,nSrc)) #((0,1) if nSrc == 1 else (0, nSrc))
for i in range(self.survey.prob.nT):
if 'b' in self:
b = self[:,'b',i+1]
else:
b = np.zeros(nF if nSrc == 1 else (nF, nSrc))
b = np.zeros((nF,nSrc)) # if nSrc == 1 else (nF, nSrc))
if 'e' in self:
e = self[:,'e',i+1]
else:
e = np.zeros(nE if nSrc == 1 else (nE, nSrc))
e = np.zeros((nE,nSrc)) # if nSrc == 1 else (nE, nSrc))
u = np.concatenate((u, b, e))
return Utils.mkvc(u)
return Utils.mkvc(u,nSrc)
class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
+15 -10
View File
@@ -74,7 +74,9 @@ class ProblemTDEM_b(BaseTDEMProblem):
def getRHS(self, tInd, F):
dt = self.timeSteps[tInd]
B_n = np.c_[[F[src,'b',tInd] for src in self.survey.srcList]].T
RHS = (1.0/dt)*self.MfMui*mkvc(B_n)
if B_n.shape[0] is not 1:
raise NotImplementedError('getRHS not implemented for this shape of B_n')
RHS = (1.0/dt)*self.MfMui*B_n[0,:,:] #TODO: This is a hack
return RHS
####################################################
@@ -106,15 +108,17 @@ class ProblemTDEM_b(BaseTDEMProblem):
# fake initial 'e' fields
p[:, 'e', 0] = 0.0
dMdsig = self.mesh.getEdgeInnerProductDeriv(self.curModel.transform)
dsigdm_x_v = self.curModel.transformDeriv*vec
dMdsig = self.MeSigmaDeriv
# self.mesh.getEdgeInnerProductDeriv(self.curModel.transform)
# dsigdm_x_v = self.curModel.sigmaDeriv*vec
# dsigdm_x_v = self.curModel.transformDeriv*vec
for i in range(1,self.nT+1):
# TODO: G[1] may be dependent on the model
# for a galvanic source (deriv of the dc problem)
#
# Do multiplication for all src in self.survey.srcList
for src in self.survey.srcList:
p[src, 'e', i] = - dMdsig(u[src,'e',i]) * dsigdm_x_v
p[src, 'e', i] = - dMdsig(u[src,'e',i]) * vec
return p
def Gtvec(self, m, vec, u=None):
@@ -130,8 +134,9 @@ class ProblemTDEM_b(BaseTDEMProblem):
if u is None:
u = self.fields(m)
self.curModel = m
dMdsig = self.mesh.getEdgeInnerProductDeriv(self.curModel.transform)
dsigdm = self.curModel.transformDeriv
# dMdsig = self.mesh.getEdgeInnerProductDeriv(self.curModel.transform)
# dsigdm = self.curModel.transformDeriv
MeSigmaDeriv = self.MeSigmaDeriv
nSrc = self.survey.nSrc
VUs = None
@@ -139,11 +144,11 @@ class ProblemTDEM_b(BaseTDEMProblem):
for i in range(1,self.nT+1):
vu = None
for src in self.survey.srcList:
vusrc = dMdsig(u[src,'e',i]).T * vec[src,'e',i]
vusrc = MeSigmaDeriv(u[src,'e',i]).T * vec[src,'e',i]
vu = vusrc if vu is None else vu + vusrc
VUs = vu if VUs is None else VUs + vu
p = -dsigdm.T*VUs
return p
# p = -dsigdm.T*VUs
return -VUs
def solveAh(self, m, p):
"""
@@ -242,7 +247,7 @@ class ProblemTDEM_b(BaseTDEMProblem):
def AhtRHS(tInd, y):
nSrc, nF = self.survey.nSrc, self.mesh.nF
rhs = np.zeros(nF if nSrc == 1 else (nF, nSrc))
rhs = np.zeros((nF,1) if nSrc == 1 else (nF, nSrc))
if 'e' in p:
rhs += self.MfMui*(self.mesh.edgeCurl*(self.MeSigmaI*p[:,'e',tInd+1]))
+5 -6
View File
@@ -108,11 +108,11 @@ def adjointTest(fdemType, comp):
prb = getProblem(fdemType, comp)
print 'Adjoint %s formulation - %s' % (fdemType, comp)
m = np.log(np.ones(prb.mesh.nC)*CONDUCTIVITY)
m = np.log(np.ones(prb.mapping.nP)*CONDUCTIVITY)
mu = np.ones(prb.mesh.nC)*MU
if addrandoms is True:
m = m + np.random.randn(prb.mesh.nC)*np.log(CONDUCTIVITY)*1e-1
m = m + np.random.randn(prb.mapping.nP)*np.log(CONDUCTIVITY)*1e-1
mu = mu + np.random.randn(prb.mesh.nC)*MU*1e-1
survey = prb.survey
@@ -121,7 +121,6 @@ def adjointTest(fdemType, comp):
u = prb.fields(m)
v = np.random.rand(survey.nD)
# print prb.PropMap.PropModel.nP
w = np.random.rand(prb.mesh.nC)
vJw = v.dot(prb.Jvec(m, w, u))
@@ -135,12 +134,12 @@ def derivTest(fdemType, comp):
prb = getProblem(fdemType, comp)
print '%s formulation - %s' % (fdemType, comp)
x0 = np.log(np.ones(prb.mesh.nC)*CONDUCTIVITY)
x0 = np.log(np.ones(prb.mapping.nP)*CONDUCTIVITY)
mu = np.log(np.ones(prb.mesh.nC)*MU)
if addrandoms is True:
x0 = x0 + np.random.randn(prb.mesh.nC)*np.log(CONDUCTIVITY)*1e-1
mu = mu + np.random.randn(prb.mesh.nC)*MU*1e-1
x0 = x0 + np.random.randn(prb.mapping.nP)*np.log(CONDUCTIVITY)*1e-1
mu = mu + np.random.randn(prb.mapping.nP)*MU*1e-1
prb.PropMap.PropModel.mu = mu
prb.PropMap.PropModel.mui = 1./mu
+1
View File
@@ -76,6 +76,7 @@ class FDEM_analyticTests(unittest.TestCase):
orderMag = 1.6
passed = np.abs(np.mean(diff - np.log10(np.abs(mu_0*np.imag(an))))) > orderMag
self.assertTrue(passed)
return passed
if __name__ == '__main__':
+17 -14
View File
@@ -3,6 +3,7 @@ from SimPEG import *
import simpegEM as EM
plotIt = False
tol = 1e-6
class TDEM_bDerivTests(unittest.TestCase):
@@ -22,7 +23,7 @@ class TDEM_bDerivTests(unittest.TestCase):
rxOffset = 40.
rx = EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20), 'bz')
src = EM.TDEM.SrcTDEM(np.array([0., 0., 0.]), 'VMD_MVP', [rx])
src = EM.TDEM.SrcTDEM_VMD_MVP([rx], loc=np.array([0., 0., 0.]))
survey = EM.TDEM.SurveyTDEM([src])
@@ -60,7 +61,7 @@ class TDEM_bDerivTests(unittest.TestCase):
self.assertLess(np.linalg.norm(V1-V2)/np.linalg.norm(V2), 1.e-6)
V1 = Ahu[:,'e',1]
self.assertLess(np.linalg.norm(V1), 1.e-6)
return np.linalg.norm(V1) < 1.e-6
for i in range(2,prb.nT):
@@ -74,7 +75,7 @@ class TDEM_bDerivTests(unittest.TestCase):
V1 = Ahu[:,'e',i]
V2 = prb.MeSigma*u[:,'e',i]
# print np.linalg.norm(V1), np.linalg.norm(V2)
self.assertLess(np.linalg.norm(V1)/np.linalg.norm(V2), 1.e-6)
return np.linalg.norm(V1)/np.linalg.norm(V2), 1.e-6
def test_AhVecVSMat_OneTS(self):
@@ -94,7 +95,7 @@ class TDEM_bDerivTests(unittest.TestCase):
u1 = A*f.tovec()
u2 = prb._AhVec(sigma,f).tovec()
self.assertTrue(np.linalg.norm(u1-u2)/np.linalg.norm(u1)<1e-12)
return np.linalg.norm(u1-u2)/np.linalg.norm(u1)<1e-12
def test_solveAhVSMat_OneTS(self):
prb = self.prb
@@ -120,7 +121,7 @@ class TDEM_bDerivTests(unittest.TestCase):
u1 = prb.solveAh(sigma,f).tovec().flatten()
u2 = sp.linalg.spsolve(A.tocsr(),f.tovec())
self.assertLess(np.linalg.norm(u1-u2),1e-8)
return np.linalg.norm(u1-u2),1e-8
def test_solveAhVsAhVec(self):
@@ -139,7 +140,7 @@ class TDEM_bDerivTests(unittest.TestCase):
u1 = f.tovec()
u2 = f_test.tovec()
self.assertTrue(np.linalg.norm(u1-u2)<1e-8)
return np.linalg.norm(u1-u2)<1e-8
def test_DerivG(self):
"""
@@ -156,7 +157,7 @@ class TDEM_bDerivTests(unittest.TestCase):
derChk = lambda m: [self.prb._AhVec(m, f).tovec(), lambda mx: self.prb.Gvec(sigma, mx, u=f).tovec()]
print '\ntest_DerivG'
passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20)
self.assertTrue(passed)
return passed
def test_Deriv_dUdM(self):
@@ -172,7 +173,7 @@ class TDEM_bDerivTests(unittest.TestCase):
print '\n'
print 'test_Deriv_dUdM'
passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20)
self.assertTrue(passed)
return passed
def test_Deriv_J(self):
@@ -189,7 +190,7 @@ class TDEM_bDerivTests(unittest.TestCase):
print '\n'
print 'test_Deriv_J'
passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=4, eps=1e-20)
self.assertTrue(passed)
return passed
def test_projectAdjoint(self):
prb = self.prb
@@ -208,7 +209,7 @@ class TDEM_bDerivTests(unittest.TestCase):
V1 = d_vec.dot(survey.projectFieldsDeriv(None, v=f).tovec())
V2 = f.tovec().dot(survey.projectFieldsDeriv(None, v=d, adjoint=True).tovec())
self.assertLess((V1-V2)/np.abs(V1), 1e-6)
return (V1-V2)/np.abs(V1) < tol
def test_adjointAhVsAht(self):
prb = self.prb
@@ -227,7 +228,7 @@ class TDEM_bDerivTests(unittest.TestCase):
V1 = f2.tovec().dot(prb._AhVec(sigma, f1).tovec())
V2 = f1.tovec().dot(prb._AhtVec(sigma, f2).tovec())
self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6)
return np.abs(V1-V2)/np.abs(V1) < tol
# def test_solveAhtVsAhtVec(self):
# prb = self.prb
@@ -292,7 +293,7 @@ class TDEM_bDerivTests(unittest.TestCase):
V1 = m.dot(prb.Gtvec(sigma, v, u))
V2 = v.tovec().dot(prb.Gvec(sigma, m, u).tovec())
self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6)
return np.abs(V1-V2)/np.abs(V1) < tol
def test_adjointJvecVsJtvec(self):
mesh = self.mesh
@@ -304,8 +305,10 @@ class TDEM_bDerivTests(unittest.TestCase):
V1 = d.dot(prb.Jvec(sigma, m))
V2 = m.dot(prb.Jtvec(sigma, d))
print 'AdjointTest', V1, V2
self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6)
passed = np.abs(V1-V2)/np.abs(V1) < tol
print 'AdjointTest', V1, V2, passed
return passed
@@ -22,9 +22,9 @@ class TDEM_bDerivTests(unittest.TestCase):
rxOffset = 40.
rx = EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20), 'bz')
src = EM.TDEM.SrcTDEM(np.array([0., 0., 0.]), 'VMD_MVP', [rx])
src = EM.TDEM.SrcTDEM_VMD_MVP( [rx], loc=np.array([0., 0., 0.]))
rx2 = EM.TDEM.RxTDEM(np.array([[rxOffset-10, 0., 0.]]), np.logspace(-5,-4, 25), 'bz')
src2 = EM.TDEM.SrcTDEM(np.array([0., 0., 0.]), 'VMD_MVP', [rx2])
src2 = EM.TDEM.SrcTDEM_VMD_MVP( [rx2], loc=np.array([0., 0., 0.]))
survey = EM.TDEM.SurveyTDEM([src,src2])
@@ -61,7 +61,7 @@ class TDEM_bDerivTests(unittest.TestCase):
derChk = lambda m: [self.prb._AhVec(m, f).tovec(), lambda mx: self.prb.Gvec(sigma, mx, u=f).tovec()]
print '\ntest_DerivG'
passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20)
self.assertTrue(passed)
return passed
def test_Deriv_dUdM(self):
@@ -77,7 +77,7 @@ class TDEM_bDerivTests(unittest.TestCase):
print '\n'
print 'test_Deriv_dUdM'
passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20)
self.assertTrue(passed)
return passed
def test_Deriv_J(self):
@@ -94,7 +94,7 @@ class TDEM_bDerivTests(unittest.TestCase):
print '\n'
print 'test_Deriv_J'
passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=4, eps=1e-20)
self.assertTrue(passed)
return passed
def test_projectAdjoint(self):
prb = self.prb
@@ -112,9 +112,9 @@ class TDEM_bDerivTests(unittest.TestCase):
# Check that d.T*Q*f = f.T*Q.T*d
V1 = d_vec.dot(survey.projectFieldsDeriv(None, v=f).tovec())
V2 = f.tovec().dot(survey.projectFieldsDeriv(None, v=d, adjoint=True).tovec())
V2 = np.sum((f.tovec())*(survey.projectFieldsDeriv(None, v=d, adjoint=True).tovec()))
self.assertLess((V1-V2)/np.abs(V1), 1e-6)
return (V1-V2)/np.abs(V1) < 1e-6
def test_adjointGvecVsGtvec(self):
mesh = self.mesh
@@ -134,8 +134,8 @@ class TDEM_bDerivTests(unittest.TestCase):
v[:,'e',i] = np.random.rand(mesh.nE, 2)
V1 = m.dot(prb.Gtvec(sigma, v, u))
V2 = v.tovec().dot(prb.Gvec(sigma, m, u).tovec())
self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6)
V2 = np.sum(v.tovec()*prb.Gvec(sigma, m, u).tovec())
return np.abs(V1-V2)/np.abs(V1) <1e-6
def test_adjointJvecVsJtvec(self):
mesh = self.mesh
@@ -148,7 +148,7 @@ class TDEM_bDerivTests(unittest.TestCase):
V1 = d.dot(prb.Jvec(sigma, m))
V2 = m.dot(prb.Jtvec(sigma, d))
print 'AdjointTest', V1, V2
self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6)
return np.abs(V1-V2)/np.abs(V1) < 1e-6
+1 -1
View File
@@ -22,7 +22,7 @@ def getProb(meshType='CYL',rxTypes='bx,bz',nSrc=1):
srcs = []
for ii in range(nSrc):
rxs = [EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20 + ii), rxType) for rxType in rxTypes.split(',')]
srcs += [EM.TDEM.SrcTDEM(np.array([0., 0., 0.]), 'VMD_MVP', rxs)]
srcs += [EM.TDEM.SrcTDEM_VMD_MVP(rxs,np.array([0., 0., 0.]))]
survey = EM.TDEM.SurveyTDEM(srcs)
+12 -11
View File
@@ -28,7 +28,8 @@ def halfSpaceProblemAnaDiff(meshType, sig_half=1e-2, rxOffset=50., bounds=[1e-5,
mapping = Maps.ExpMap(mesh) * Maps.Vertical1DMap(mesh) * actMap
rx = EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-5,-4, 21), 'bz')
src = EM.TDEM.SrcTDEM(np.array([0., 0., 0.]), 'VMD_MVP', [rx])
src = EM.TDEM.SrcTDEM_VMD_MVP([rx], loc=np.array([0., 0., 0.]))
# src = EM.TDEM.SrcTDEM([rx], loc=np.array([0., 0., 0.]))
survey = EM.TDEM.SurveyTDEM([src])
prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping)
@@ -60,26 +61,26 @@ def halfSpaceProblemAnaDiff(meshType, sig_half=1e-2, rxOffset=50., bounds=[1e-5,
class TDEM_bTests(unittest.TestCase):
def test_analitic_p2_CYL_50m(self):
def test_analytic_p2_CYL_50m(self):
self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=50., sig_half=1e+2) < 0.01)
def test_analitic_p1_CYL_50m(self):
def test_analytic_p1_CYL_50m(self):
self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=50., sig_half=1e+1) < 0.01)
def test_analitic_p0_CYL_50m(self):
def test_analytic_p0_CYL_50m(self):
self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=50., sig_half=1e+0) < 0.01)
def test_analitic_m1_CYL_50m(self):
def test_analytic_m1_CYL_50m(self):
self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=50., sig_half=1e-1) < 0.01)
def test_analitic_m2_CYL_50m(self):
def test_analytic_m2_CYL_50m(self):
self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=50., sig_half=1e-2) < 0.01)
def test_analitic_m3_CYL_50m(self):
def test_analytic_m3_CYL_50m(self):
self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=50., sig_half=1e-3) < 0.02)
def test_analitic_p0_CYL_1m(self):
def test_analytic_p0_CYL_1m(self):
self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=1.0, sig_half=1e+0) < 0.01)
def test_analitic_m1_CYL_1m(self):
def test_analytic_m1_CYL_1m(self):
self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=1.0, sig_half=1e-1) < 0.01)
def test_analitic_m2_CYL_1m(self):
def test_analytic_m2_CYL_1m(self):
self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=1.0, sig_half=1e-2) < 0.01)
def test_analitic_m3_CYL_1m(self):
def test_analytic_m3_CYL_1m(self):
self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=1.0, sig_half=1e-3) < 0.02)