From edc241eec2463232847066a53f6caea120de1a85 Mon Sep 17 00:00:00 2001 From: Dave Marchant Date: Tue, 18 Feb 2014 20:34:34 -0800 Subject: [PATCH] Added Gtvec method to TDEM. Renamed G to Gvec. Added tests. --- simpegEM/TDEM/TDEM_b.py | 17 +++++++++++++--- simpegEM/Tests/test_forward_EMproblem.py | 25 +++++++++++++++++++++--- 2 files changed, 36 insertions(+), 6 deletions(-) diff --git a/simpegEM/TDEM/TDEM_b.py b/simpegEM/TDEM/TDEM_b.py index 35aaf5c9..613daab9 100644 --- a/simpegEM/TDEM/TDEM_b.py +++ b/simpegEM/TDEM/TDEM_b.py @@ -1,5 +1,6 @@ from BaseTDEM import ProblemBaseTDEM from FieldsTDEM import FieldsTDEM +from SimPEG.Utils import mkvc import numpy as np class ProblemTDEM_b(ProblemBaseTDEM): @@ -43,15 +44,15 @@ class ProblemTDEM_b(ProblemBaseTDEM): def Jvec(self, m, v, u=None): if u is None: u = self.fields(m) - p = self.G(m, v, u) + p = self.Gvec(m, v, u) y = self.solveAh(m, p) return self.data.dpred(m, u=y) - def G(self, m, v, u=None): + def Gvec(self, m, v, u=None): if u is None: u = self.fields(m) p = FieldsTDEM(self.mesh, 1, self.times.size, 'b') - c = self.mesh.getEdgeMassDeriv()*self.model.transformDeriv(m)*v + c = self.mesh.getEdgeMassDeriv()*self.model.transformDeriv(None)*v for i in range(self.times.size): ei = u.get_e(i) pVal = np.empty_like(ei) @@ -62,6 +63,16 @@ class ProblemTDEM_b(ProblemBaseTDEM): p.set_b(np.zeros((self.mesh.nF,1)), i) return p + def Gtvec(self, m, v, u=None): + if u is None: + u = self.fields(m) + tmp = np.zeros((self.mesh.nE,self.data.nTx)) + for i in range(self.nTimes): + tmp += v.get_e(i)*u.get_e(i) + p = -mkvc(self.model.transformDeriv(None).T*self.mesh.getEdgeMassDeriv().T*tmp) + return p + + def solveAh(self, m, p): def AhRHS(tInd, u): rhs = self.MfMui*self.mesh.edgeCurl*self.MeSigmaI*p.get_e(tInd) + p.get_b(tInd) diff --git a/simpegEM/Tests/test_forward_EMproblem.py b/simpegEM/Tests/test_forward_EMproblem.py index 95dcb72f..4de5d920 100644 --- a/simpegEM/Tests/test_forward_EMproblem.py +++ b/simpegEM/Tests/test_forward_EMproblem.py @@ -180,7 +180,7 @@ class TDEM_bDerivTests(unittest.TestCase): h = 1. a = np.linalg.norm(self.prb.AhVec(sigma+h*dm, f).fieldVec() - self.prb.AhVec(sigma, f).fieldVec()) - b = np.linalg.norm(self.prb.AhVec(sigma+h*dm, f).fieldVec() - self.prb.AhVec(sigma, f).fieldVec() - h*self.prb.G(sigma, dm, u=f).fieldVec()) + b = np.linalg.norm(self.prb.AhVec(sigma+h*dm, f).fieldVec() - self.prb.AhVec(sigma, f).fieldVec() - h*self.prb.Gvec(sigma, dm, u=f).fieldVec()) # Assuming that the gradient is exact to machine precision self.assertTrue(b<1e-16) @@ -202,7 +202,7 @@ class TDEM_bDerivTests(unittest.TestCase): for i, h in enumerate(hv): f = prb.fields(sigma) fstep = prb.fields(sigma + h*d_sig) - dcdm = prb.G(sigma, h*d_sig, u=f) # TODO: make negative!?!? + dcdm = prb.Gvec(sigma, h*d_sig, u=f) # TODO: make negative!?!? dudm = prb.solveAh(sigma, dcdm) linear = np.linalg.norm(f.fieldVec() - fstep.fieldVec()) @@ -305,7 +305,26 @@ class TDEM_bDerivTests(unittest.TestCase): V2 = f1.fieldVec().dot(prb.solveAht(sigma, f2).fieldVec()) self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-8) - + def test_adjointGvecVsGtvec(self): + mesh = self.mesh + prb = self.prb + + m = np.random.rand(mesh.nCz) + sigma = np.random.rand(mesh.nCz) + + u = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nTimes, 'b') + for i in range(u.nTimes): + u.set_b(np.random.rand(mesh.nF, 1), i) + u.set_e(np.random.rand(mesh.nE, 1), i) + + v = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nTimes, 'b') + for i in range(v.nTimes): + v.set_b(np.random.rand(mesh.nF, 1), i) + v.set_e(np.random.rand(mesh.nE, 1), i) + + V1 = m.dot(prb.Gtvec(sigma, v, u)) + V2 = v.fieldVec().dot(prb.Gvec(sigma, m, u).fieldVec()) + self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-8) if __name__ == '__main__':