diff --git a/simpegEM/TDEM/TDEM_b.py b/simpegEM/TDEM/TDEM_b.py index 7f00485b..2eada637 100644 --- a/simpegEM/TDEM/TDEM_b.py +++ b/simpegEM/TDEM/TDEM_b.py @@ -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 diff --git a/simpegEM/Tests/test_forward_EMproblem.py b/simpegEM/Tests/test_forward_EMproblem.py index 1f5750d8..42ad57dd 100644 --- a/simpegEM/Tests/test_forward_EMproblem.py +++ b/simpegEM/Tests/test_forward_EMproblem.py @@ -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()