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

Added AhtVec method & test.

Browse Source
This commit is contained in:
Dave Marchant
2014-02-18 15:49:54 -08:00
parent 1128fbe39b
commit 7bb70103aa
2 changed files with 39 additions and 3 deletions
Download Patch File Download Diff File Expand all files Collapse all files
simpegEM/TDEM/TDEM_b.py
+20 -1
View File
@@ -93,7 +93,7 @@ class ProblemTDEM_b(ProblemBaseTDEM):
f = FieldsTDEM(self.mesh, 1, self.times.size, 'b')
f.set_b(b, 0)
f.set_e(e, 0)
for i in range(1,self.times.size):
for i in range(1,self.nTimes):
dt = self.getDt(i)
b = 1/dt*self.MfMui*u.get_b(i) + self.MfMui*self.mesh.edgeCurl*u.get_e(i) - 1/dt*self.MfMui*u.get_b(i-1)
e = self.mesh.edgeCurl.T*self.MfMui*u.get_b(i) - self.MeSigma*u.get_e(i)
@@ -101,6 +101,25 @@ class ProblemTDEM_b(ProblemBaseTDEM):
f.set_e(e, i)
return f
def AhtVec(self, m, u=None):
if u is None:
u = self.fields(m)
self.makeMassMatrices(m)
f = FieldsTDEM(self.mesh, 1, self.times.size, 'b')
for i in range(self.nTimes-1):
b = 1/self.getDt(i)*self.MfMui*u.get_b(i) + self.MfMui*self.mesh.edgeCurl*u.get_e(i) - 1/self.getDt(i+1)*self.MfMui*u.get_b(i+1)
e = self.mesh.edgeCurl.T*self.MfMui*u.get_b(i) - self.MeSigma*u.get_e(i)
f.set_b(b, i)
f.set_e(e, i)
N = self.nTimes - 1
b = 1/self.getDt(N)*self.MfMui*u.get_b(N) + self.MfMui*self.mesh.edgeCurl*u.get_e(N)
e = self.mesh.edgeCurl.T*self.MfMui*u.get_b(N) - self.MeSigma*u.get_e(N)
f.set_b(b, N)
f.set_e(e, N)
return f
if __name__ == '__main__':
from SimPEG import *
import simpegEM as EM
simpegEM/Tests/test_forward_EMproblem.py
+19 -2
View File
@@ -69,7 +69,6 @@ class TDEM_bDerivTests(unittest.TestCase):
self.prb.pair(self.dat)
self.mesh = mesh
def test_AhVec(self):
"""
Test that fields and AhVec produce consistent results
@@ -233,7 +232,6 @@ class TDEM_bDerivTests(unittest.TestCase):
passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=2, eps=1e-20)
self.assertTrue(passed)
def test_projectAdjoint(self):
prb = self.prb
dat = self.dat
@@ -252,6 +250,25 @@ class TDEM_bDerivTests(unittest.TestCase):
self.assertTrue((V1-V2)/V1<1e-12)
def test_adjointAhVsAht(self):
prb = self.prb
mesh = self.mesh
sigma = self.sigma
f1 = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nTimes, 'b')
for i in range(f1.nTimes):
f1.set_b(np.random.rand(mesh.nF, 1), i)
f1.set_e(np.random.rand(mesh.nE, 1), i)
f2 = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nTimes, 'b')
for i in range(f2.nTimes):
f2.set_b(np.random.rand(mesh.nF, 1), i)
f2.set_e(np.random.rand(mesh.nE, 1), i)
V1 = f2.fieldVec().dot(prb.AhVec(sigma, f1).fieldVec())
V2 = f1.fieldVec().dot(prb.AhtVec(sigma, f2).fieldVec())
self.assertLess(np.abs(V1-V2)/V1, 1e-12)
if __name__ == '__main__':
unittest.main()