mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-25 13:30:06 +08:00
Adjoint test working.
This commit is contained in:
rowanc1
+13
-5
@@ -137,19 +137,27 @@ class ProblemTDEM_b(BaseTDEMProblem):
|
||||
|
||||
def solveAht(self, m, p):
|
||||
|
||||
def AhtRHS(tInd, u):
|
||||
# fLoc 0 1 2 3
|
||||
# |-----|-----|-----|
|
||||
# tInd 0 1 2 / /
|
||||
# / __/
|
||||
# 2 (tInd=2 uses fields 3 and would use 4 but it doesn't exist)
|
||||
# / __/
|
||||
# 1 (tInd=1 uses fields 2 and 3)
|
||||
|
||||
def AhtRHS(tInd, y):
|
||||
nTx, nF = self.survey.nTx, self.mesh.nF
|
||||
rhs = np.zeros(nF if nTx == 1 else (nF, nTx))
|
||||
|
||||
if 'e' in p:
|
||||
rhs += self.MfMui*self.mesh.edgeCurl*self.MeSigmaI*p[:,'e',tInd]
|
||||
rhs += self.MfMui*self.mesh.edgeCurl*self.MeSigmaI*p[:,'e',tInd+1]
|
||||
if 'b' in p:
|
||||
rhs += p[:,'b',tInd]
|
||||
rhs += p[:,'b',tInd+1]
|
||||
|
||||
if tInd == self.nT-1:
|
||||
return rhs
|
||||
dt = self.timeSteps[tInd]
|
||||
return rhs + 1.0/dt*self.MfMui*u[:,'b',tInd+1]
|
||||
dt = self.timeSteps[tInd+1]
|
||||
return rhs + 1.0/dt*self.MfMui*y[:,'b',tInd+2]
|
||||
|
||||
def AhtCalcFields(sol, solType, tInd):
|
||||
y_b = sol
|
||||
|
||||
@@ -28,7 +28,9 @@ class TDEM_bDerivTests(unittest.TestCase):
|
||||
survey = EM.TDEM.SurveyTDEM([tx])
|
||||
|
||||
self.prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping)
|
||||
# self.prb.timeSteps = [1e-5]
|
||||
self.prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)]
|
||||
# self.prb.timeSteps = [(1e-05, 100)]
|
||||
|
||||
self.sigma = np.ones(mesh.nCz)*1e-8
|
||||
self.sigma[mesh.vectorCCz<0] = 1e-1
|
||||
@@ -37,217 +39,174 @@ class TDEM_bDerivTests(unittest.TestCase):
|
||||
self.prb.pair(survey)
|
||||
self.mesh = mesh
|
||||
|
||||
# def test_AhVec(self):
|
||||
# """
|
||||
# Test that fields and AhVec produce consistent results
|
||||
# """
|
||||
def test_AhVec(self):
|
||||
"""
|
||||
Test that fields and AhVec produce consistent results
|
||||
"""
|
||||
|
||||
# prb = self.prb
|
||||
# sigma = self.sigma
|
||||
|
||||
# u = prb.fields(sigma)
|
||||
# Ahu = prb.AhVec(sigma, u)
|
||||
|
||||
# V1 = Ahu[:,'b',1]
|
||||
# V2 = 1./prb.timeSteps[0]*prb.MfMui*u[:,'b',0]
|
||||
# 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)
|
||||
|
||||
# for i in range(2,prb.nT):
|
||||
|
||||
# dt = prb.timeSteps[i]
|
||||
|
||||
# V1 = Ahu[:,'b',i]
|
||||
# V2 = 1.0/dt*prb.MfMui*u[:,'b', i-1]
|
||||
# # print np.linalg.norm(V1), np.linalg.norm(V2)
|
||||
# self.assertLess(np.linalg.norm(V1)/np.linalg.norm(V2), 1.e-6)
|
||||
|
||||
# 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)
|
||||
|
||||
# def test_AhVecVSMat_OneTS(self):
|
||||
|
||||
# prb = self.prb
|
||||
# prb.timeSteps = [1e-05]
|
||||
# sigma = self.sigma
|
||||
# prb.curModel = sigma
|
||||
|
||||
# dt = prb.timeSteps[0]
|
||||
# a11 = 1/dt*prb.MfMui*sp.eye(prb.mesh.nF)
|
||||
# a12 = prb.MfMui*prb.mesh.edgeCurl
|
||||
# a21 = prb.mesh.edgeCurl.T*prb.MfMui
|
||||
# a22 = -prb.MeSigma
|
||||
# A = sp.bmat([[a11,a12],[a21,a22]])
|
||||
|
||||
# f = prb.fields(sigma)
|
||||
# u1 = A*f.tovec()
|
||||
# u2 = prb.AhVec(sigma,f).tovec()
|
||||
|
||||
# self.assertTrue(np.linalg.norm(u1-u2)/np.linalg.norm(u1)<1e-12)
|
||||
|
||||
# def test_solveAhVSMat_OneTS(self):
|
||||
# prb = self.prb
|
||||
|
||||
# prb.timeSteps = [1e-05]
|
||||
|
||||
# sigma = self.sigma
|
||||
# prb.curModel = sigma
|
||||
|
||||
# dt = prb.timeSteps[0]
|
||||
# a11 = 1.0/dt*prb.MfMui*sp.eye(prb.mesh.nF)
|
||||
# a12 = prb.MfMui*prb.mesh.edgeCurl
|
||||
# a21 = prb.mesh.edgeCurl.T*prb.MfMui
|
||||
# a22 = -prb.MeSigma
|
||||
# A = sp.bmat([[a11,a12],[a21,a22]])
|
||||
|
||||
# f = prb.fields(sigma)
|
||||
# f[:,:,0] = {'e':0,'b':0}
|
||||
# f[:,'b',1] = 0
|
||||
# f[:,'e',1] = np.random.rand(prb.mesh.nE,1)
|
||||
|
||||
# self.assertTrue(np.all(np.r_[f[:,'b',1],f[:,'e',1]] == f.tovec()))
|
||||
|
||||
# u1 = prb.solveAh(sigma,f).tovec().flatten()
|
||||
# u2 = sp.linalg.spsolve(A.tocsr(),f.tovec())
|
||||
|
||||
# self.assertLess(np.linalg.norm(u1-u2),1e-8)
|
||||
|
||||
# def test_solveAhVsAhVec(self):
|
||||
|
||||
# prb = self.prb
|
||||
# mesh = self.prb.mesh
|
||||
# sigma = self.sigma
|
||||
# self.prb.curModel = sigma
|
||||
|
||||
# f = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
|
||||
# f[:,'b',:] = 0.0
|
||||
# for i in range(prb.nT):
|
||||
# f[:,'e', i] = np.random.rand(mesh.nE, 1)
|
||||
|
||||
# Ahf = prb.AhVec(sigma, f)
|
||||
# f_test = prb.solveAh(sigma, Ahf)
|
||||
|
||||
# u1 = f.tovec()
|
||||
# u2 = f_test.tovec()
|
||||
# self.assertTrue(np.linalg.norm(u1-u2)<1e-8)
|
||||
|
||||
# def test_DerivG(self):
|
||||
# """
|
||||
# Test the derivative of c with respect to sigma
|
||||
# """
|
||||
|
||||
# # Random model and perturbation
|
||||
# sigma = np.random.rand(self.prb.mapping.nP)
|
||||
|
||||
# f = self.prb.fields(sigma)
|
||||
# dm = 1000*np.random.rand(self.prb.mapping.nP)
|
||||
# h = 0.01
|
||||
|
||||
# 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)
|
||||
|
||||
# def test_Deriv_dUdM(self):
|
||||
|
||||
# prb = self.prb
|
||||
# prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)]
|
||||
# mesh = self.mesh
|
||||
# sigma = self.sigma
|
||||
|
||||
# dm = 10*np.random.rand(prb.mapping.nP)
|
||||
# f = prb.fields(sigma)
|
||||
|
||||
# derChk = lambda m: [self.prb.fields(m).tovec(), lambda mx: -prb.solveAh(sigma, prb.Gvec(sigma, mx, u=f)).tovec()]
|
||||
# print '\n'
|
||||
# print 'test_Deriv_dUdM'
|
||||
# passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20)
|
||||
# self.assertTrue(passed)
|
||||
|
||||
# def test_Deriv_J(self):
|
||||
|
||||
# prb = self.prb
|
||||
# prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)]
|
||||
# mesh = self.mesh
|
||||
# sigma = self.sigma
|
||||
|
||||
# # d_sig = 0.8*sigma #np.random.rand(mesh.nCz)
|
||||
# d_sig = 10*np.random.rand(prb.mapping.nP)
|
||||
|
||||
|
||||
# derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(sigma, mx)]
|
||||
# print '\n'
|
||||
# print 'test_Deriv_J'
|
||||
# passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=4, eps=1e-20)
|
||||
# self.assertTrue(passed)
|
||||
|
||||
# def test_projectAdjoint(self):
|
||||
# prb = self.prb
|
||||
# survey = prb.survey
|
||||
# mesh = self.mesh
|
||||
|
||||
# # Generate random fields and data
|
||||
# f = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
|
||||
# for i in range(prb.nT):
|
||||
# f[:,'b',i] = np.random.rand(mesh.nF, 1)
|
||||
# f[:,'e',i] = np.random.rand(mesh.nE, 1)
|
||||
# d_vec = np.random.rand(survey.nD, survey.nTx).flatten()
|
||||
# d = Survey.Data(survey,v=d_vec)
|
||||
|
||||
# # 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())
|
||||
|
||||
# self.assertLess((V1-V2)/np.abs(V1), 1e-6)
|
||||
|
||||
# def test_adjointAhVsAht(self):
|
||||
# prb = self.prb
|
||||
# mesh = self.mesh
|
||||
# sigma = self.sigma
|
||||
|
||||
# f1 = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
|
||||
# for i in range(1,prb.nT+1):
|
||||
# f1[:,'b',i] = np.random.rand(mesh.nF, 1)
|
||||
# f1[:,'e',i] = np.random.rand(mesh.nE, 1)
|
||||
|
||||
# f2 = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
|
||||
# for i in range(1,prb.nT+1):
|
||||
# f2[:,'b',i] = np.random.rand(mesh.nF, 1)
|
||||
# f2[:,'e',i] = np.random.rand(mesh.nE, 1)
|
||||
|
||||
# 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)
|
||||
|
||||
def test_solveAhtVsAhtVec(self):
|
||||
prb = self.prb
|
||||
mesh = self.mesh
|
||||
sigma = np.random.rand(prb.mapping.nP)
|
||||
sigma = self.sigma
|
||||
|
||||
f1 = EM.TDEM.FieldsTDEM(mesh,prb.survey)
|
||||
u = prb.fields(sigma)
|
||||
Ahu = prb.AhVec(sigma, u)
|
||||
|
||||
V1 = Ahu[:,'b',1]
|
||||
V2 = 1./prb.timeSteps[0]*prb.MfMui*u[:,'b',0]
|
||||
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)
|
||||
|
||||
for i in range(2,prb.nT):
|
||||
|
||||
dt = prb.timeSteps[i]
|
||||
|
||||
V1 = Ahu[:,'b',i]
|
||||
V2 = 1.0/dt*prb.MfMui*u[:,'b', i-1]
|
||||
# print np.linalg.norm(V1), np.linalg.norm(V2)
|
||||
self.assertLess(np.linalg.norm(V1)/np.linalg.norm(V2), 1.e-6)
|
||||
|
||||
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)
|
||||
|
||||
def test_AhVecVSMat_OneTS(self):
|
||||
|
||||
prb = self.prb
|
||||
prb.timeSteps = [1e-05]
|
||||
sigma = self.sigma
|
||||
prb.curModel = sigma
|
||||
|
||||
dt = prb.timeSteps[0]
|
||||
a11 = 1/dt*prb.MfMui*sp.eye(prb.mesh.nF)
|
||||
a12 = prb.MfMui*prb.mesh.edgeCurl
|
||||
a21 = prb.mesh.edgeCurl.T*prb.MfMui
|
||||
a22 = -prb.MeSigma
|
||||
A = sp.bmat([[a11,a12],[a21,a22]])
|
||||
|
||||
f = prb.fields(sigma)
|
||||
u1 = A*f.tovec()
|
||||
u2 = prb.AhVec(sigma,f).tovec()
|
||||
|
||||
self.assertTrue(np.linalg.norm(u1-u2)/np.linalg.norm(u1)<1e-12)
|
||||
|
||||
def test_solveAhVSMat_OneTS(self):
|
||||
prb = self.prb
|
||||
|
||||
prb.timeSteps = [1e-05]
|
||||
|
||||
sigma = self.sigma
|
||||
prb.curModel = sigma
|
||||
|
||||
dt = prb.timeSteps[0]
|
||||
a11 = 1.0/dt*prb.MfMui*sp.eye(prb.mesh.nF)
|
||||
a12 = prb.MfMui*prb.mesh.edgeCurl
|
||||
a21 = prb.mesh.edgeCurl.T*prb.MfMui
|
||||
a22 = -prb.MeSigma
|
||||
A = sp.bmat([[a11,a12],[a21,a22]])
|
||||
|
||||
f = prb.fields(sigma)
|
||||
f[:,:,0] = {'e':0,'b':0}
|
||||
f[:,'b',1] = 0
|
||||
f[:,'e',1] = np.random.rand(prb.mesh.nE,1)
|
||||
|
||||
self.assertTrue(np.all(np.r_[f[:,'b',1],f[:,'e',1]] == f.tovec()))
|
||||
|
||||
u1 = prb.solveAh(sigma,f).tovec().flatten()
|
||||
u2 = sp.linalg.spsolve(A.tocsr(),f.tovec())
|
||||
|
||||
self.assertLess(np.linalg.norm(u1-u2),1e-8)
|
||||
|
||||
def test_solveAhVsAhVec(self):
|
||||
|
||||
prb = self.prb
|
||||
mesh = self.prb.mesh
|
||||
sigma = self.sigma
|
||||
self.prb.curModel = sigma
|
||||
|
||||
f = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
|
||||
f[:,'b',:] = 0.0
|
||||
for i in range(prb.nT):
|
||||
f1[:,'b',i] = np.random.rand(mesh.nF, 1)
|
||||
f1[:,'e',i] = np.random.rand(mesh.nE, 1)
|
||||
f[:,'e', i] = np.random.rand(mesh.nE, 1)
|
||||
|
||||
f2 = prb.solveAht(sigma, f1)
|
||||
f3 = prb.AhtVec(sigma, f2)
|
||||
Ahf = prb.AhVec(sigma, f)
|
||||
f_test = prb.solveAh(sigma, Ahf)
|
||||
|
||||
if plotIt:
|
||||
import matplotlib.pyplot as plt
|
||||
plt.plot(f3.tovec())
|
||||
plt.plot(f1.tovec())
|
||||
plt.show()
|
||||
V1 = np.linalg.norm(f3.tovec()-f1.tovec())
|
||||
V2 = np.linalg.norm(f1.tovec())
|
||||
print V1, V2
|
||||
print 'I am gunna fail this one: boo. :('
|
||||
self.assertLess(V1/V2, 1e-6)
|
||||
u1 = f.tovec()
|
||||
u2 = f_test.tovec()
|
||||
self.assertTrue(np.linalg.norm(u1-u2)<1e-8)
|
||||
|
||||
def test_adjointsolveAhVssolveAht(self):
|
||||
def test_DerivG(self):
|
||||
"""
|
||||
Test the derivative of c with respect to sigma
|
||||
"""
|
||||
|
||||
# Random model and perturbation
|
||||
sigma = np.random.rand(self.prb.mapping.nP)
|
||||
|
||||
f = self.prb.fields(sigma)
|
||||
dm = 1000*np.random.rand(self.prb.mapping.nP)
|
||||
h = 0.01
|
||||
|
||||
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)
|
||||
|
||||
def test_Deriv_dUdM(self):
|
||||
|
||||
prb = self.prb
|
||||
prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)]
|
||||
mesh = self.mesh
|
||||
sigma = self.sigma
|
||||
|
||||
dm = 10*np.random.rand(prb.mapping.nP)
|
||||
f = prb.fields(sigma)
|
||||
|
||||
derChk = lambda m: [self.prb.fields(m).tovec(), lambda mx: -prb.solveAh(sigma, prb.Gvec(sigma, mx, u=f)).tovec()]
|
||||
print '\n'
|
||||
print 'test_Deriv_dUdM'
|
||||
passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20)
|
||||
self.assertTrue(passed)
|
||||
|
||||
def test_Deriv_J(self):
|
||||
|
||||
prb = self.prb
|
||||
prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)]
|
||||
mesh = self.mesh
|
||||
sigma = self.sigma
|
||||
|
||||
# d_sig = 0.8*sigma #np.random.rand(mesh.nCz)
|
||||
d_sig = 10*np.random.rand(prb.mapping.nP)
|
||||
|
||||
|
||||
derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(sigma, mx)]
|
||||
print '\n'
|
||||
print 'test_Deriv_J'
|
||||
passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=4, eps=1e-20)
|
||||
self.assertTrue(passed)
|
||||
|
||||
def test_projectAdjoint(self):
|
||||
prb = self.prb
|
||||
survey = prb.survey
|
||||
mesh = self.mesh
|
||||
|
||||
# Generate random fields and data
|
||||
f = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
|
||||
for i in range(prb.nT):
|
||||
f[:,'b',i] = np.random.rand(mesh.nF, 1)
|
||||
f[:,'e',i] = np.random.rand(mesh.nE, 1)
|
||||
d_vec = np.random.rand(survey.nD, survey.nTx).flatten()
|
||||
d = Survey.Data(survey,v=d_vec)
|
||||
|
||||
# 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())
|
||||
|
||||
self.assertLess((V1-V2)/np.abs(V1), 1e-6)
|
||||
|
||||
def test_adjointAhVsAht(self):
|
||||
prb = self.prb
|
||||
mesh = self.mesh
|
||||
sigma = self.sigma
|
||||
@@ -262,43 +221,88 @@ class TDEM_bDerivTests(unittest.TestCase):
|
||||
f2[:,'b',i] = np.random.rand(mesh.nF, 1)
|
||||
f2[:,'e',i] = np.random.rand(mesh.nE, 1)
|
||||
|
||||
V1 = f2.tovec().dot(prb.solveAh(sigma, f1).tovec())
|
||||
V2 = f1.tovec().dot(prb.solveAht(sigma, f2).tovec())
|
||||
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)
|
||||
|
||||
# def test_adjointGvecVsGtvec(self):
|
||||
# mesh = self.mesh
|
||||
# def test_solveAhtVsAhtVec(self):
|
||||
# prb = self.prb
|
||||
|
||||
# m = np.random.rand(prb.mapping.nP)
|
||||
# mesh = self.mesh
|
||||
# sigma = np.random.rand(prb.mapping.nP)
|
||||
|
||||
# u = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
|
||||
# f1 = EM.TDEM.FieldsTDEM(mesh,prb.survey)
|
||||
# for i in range(1,prb.nT+1):
|
||||
# u[:,'b',i] = np.random.rand(mesh.nF, 1)
|
||||
# u[:,'e',i] = np.random.rand(mesh.nE, 1)
|
||||
# f1[:,'b',i] = np.random.rand(mesh.nF, 1)
|
||||
# f1[:,'e',i] = np.random.rand(mesh.nE, 1)
|
||||
|
||||
# v = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
|
||||
# for i in range(1,prb.nT+1):
|
||||
# v[:,'b',i] = np.random.rand(mesh.nF, 1)
|
||||
# v[:,'e',i] = np.random.rand(mesh.nE, 1)
|
||||
# f2 = prb.solveAht(sigma, f1)
|
||||
# f3 = prb.AhtVec(sigma, f2)
|
||||
|
||||
# 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)
|
||||
# if True:
|
||||
# import matplotlib.pyplot as plt
|
||||
# plt.plot(f3.tovec(),'b')
|
||||
# plt.plot(f1.tovec(),'r')
|
||||
# plt.show()
|
||||
# V1 = np.linalg.norm(f3.tovec()-f1.tovec())
|
||||
# V2 = np.linalg.norm(f1.tovec())
|
||||
# print 'AhtVsAhtVec', V1, V2, f1.tovec()
|
||||
# print 'I am gunna fail this one: boo. :('
|
||||
# self.assertLess(V1/V2, 1e-6)
|
||||
|
||||
# def test_adjointJvecVsJtvec(self):
|
||||
# mesh = self.mesh
|
||||
# def test_adjointsolveAhVssolveAht(self):
|
||||
# prb = self.prb
|
||||
# mesh = self.mesh
|
||||
# sigma = self.sigma
|
||||
|
||||
# m = np.random.rand(prb.mapping.nP)
|
||||
# d = np.random.rand(prb.survey.nD)
|
||||
# f1 = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
|
||||
# for i in range(1,prb.nT+1):
|
||||
# f1[:,'b',i] = np.random.rand(mesh.nF, 1)
|
||||
# f1[:,'e',i] = np.random.rand(mesh.nE, 1)
|
||||
|
||||
# V1 = d.dot(prb.Jvec(sigma, m))
|
||||
# V2 = m.dot(prb.Jtvec(sigma, d))
|
||||
# f2 = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
|
||||
# for i in range(1,prb.nT+1):
|
||||
# f2[:,'b',i] = np.random.rand(mesh.nF, 1)
|
||||
# f2[:,'e',i] = np.random.rand(mesh.nE, 1)
|
||||
|
||||
# V1 = f2.tovec().dot(prb.solveAh(sigma, f1).tovec())
|
||||
# V2 = f1.tovec().dot(prb.solveAht(sigma, f2).tovec())
|
||||
# print V1, V2
|
||||
# self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6)
|
||||
|
||||
def test_adjointGvecVsGtvec(self):
|
||||
mesh = self.mesh
|
||||
prb = self.prb
|
||||
|
||||
m = np.random.rand(prb.mapping.nP)
|
||||
sigma = np.random.rand(prb.mapping.nP)
|
||||
|
||||
u = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
|
||||
for i in range(1,prb.nT+1):
|
||||
u[:,'b',i] = np.random.rand(mesh.nF, 1)
|
||||
u[:,'e',i] = np.random.rand(mesh.nE, 1)
|
||||
|
||||
v = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
|
||||
for i in range(1,prb.nT+1):
|
||||
v[:,'b',i] = np.random.rand(mesh.nF, 1)
|
||||
v[:,'e',i] = np.random.rand(mesh.nE, 1)
|
||||
|
||||
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)
|
||||
|
||||
def test_adjointJvecVsJtvec(self):
|
||||
mesh = self.mesh
|
||||
prb = self.prb
|
||||
sigma = self.sigma
|
||||
|
||||
m = np.random.rand(prb.mapping.nP)
|
||||
d = np.random.rand(prb.survey.nD)
|
||||
|
||||
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)
|
||||
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
|
||||
Reference in New Issue
Block a user