diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py index af27bf09..afee51b1 100644 --- a/simpegEM/FDEM/FDEM.py +++ b/simpegEM/FDEM/FDEM.py @@ -97,22 +97,20 @@ class ProblemFDEM_e(Problem.BaseProblem): u = self.fields(m) Jv = self.dataPair(self.survey) + dsig_dm = self.model.transformDeriv(m) for i, freq in enumerate(self.survey.freqs): e = u[freq, 'e'] A = self.getA(freq) solver = Solver(A, options=self.solveOpts) - for tx in self.survey.getTransmitters(freq): - dMe_dsig = self.mesh.getEdgeInnerProductDeriv(m, v=e) - dsig_dm = self.model.transformDeriv(m) + for txi, tx in enumerate(self.survey.getTransmitters(freq)): + dMe_dsig = self.mesh.getEdgeInnerProductDeriv(m, v=e[:,txi]) + P = tx.projectFieldsDeriv(self.mesh, u) b = 1j*omega(freq) * ( dMe_dsig * ( dsig_dm * v ) ) - Ab = solver.solve(b) - - P = tx.projectFieldsDeriv(self.mesh, u) - - Jv[tx] = -P*Ab + Ainvb = solver.solve(b) + Jv[tx] = -P*Ainvb return Utils.mkvc(Jv) @@ -120,7 +118,28 @@ class ProblemFDEM_e(Problem.BaseProblem): def Jtvec(self, m, v, u=None): if u is None: u = self.fields(m) - raise NotImplementedError('Jtvec todo!') + + # Ensure v is a data object. + if not isinstance(v, self.dataPair): + v = self.dataPair(self.survey, v) + + Jtv = np.zeros(self.model.nP, dtype=complex) + + dsig_dm = self.model.transformDeriv(m) + + for i, freq in enumerate(self.survey.freqs): + e = u[freq, 'e'] + AT = self.getA(freq).T + solver = Solver(AT, options=self.solveOpts) + + for txi, tx in enumerate(self.survey.getTransmitters(freq)): + dMe_dsig = self.mesh.getEdgeInnerProductDeriv(m, v=e[:,txi]) + + P = tx.projectFieldsDeriv(self.mesh, u) + w = solver.solve(P.T * v[tx]) + Jtv += - 1j*omega(freq) * ( dsig_dm.T * ( dMe_dsig.T * w ) ) + + return Jtv diff --git a/simpegEM/Tests/test_FDEM.py b/simpegEM/Tests/test_FDEM.py index 836ef717..a864b66b 100644 --- a/simpegEM/Tests/test_FDEM.py +++ b/simpegEM/Tests/test_FDEM.py @@ -2,6 +2,8 @@ import unittest from SimPEG import * import simpegEM as EM +TOL = 1e-10 + class FDEM_bDerivTests(unittest.TestCase): def setUp(self): @@ -26,7 +28,12 @@ class FDEM_bDerivTests(unittest.TestCase): rxList = EM.FDEM.RxListFDEM(XYZ, 'Ex') Tx0 = EM.FDEM.TxFDEM(None, 'VMD', 3, rxList) - survey = EM.FDEM.SurveyFDEM([Tx0]) + x = np.linspace(5,10,3) + XYZ = Utils.ndgrid(x,x,np.r_[0]) + rxList = EM.FDEM.RxListFDEM(XYZ, 'Ey') + Tx1 = EM.FDEM.TxFDEM(None, 'VMD', 3, rxList) + + survey = EM.FDEM.SurveyFDEM([Tx0, Tx1]) prb = EM.FDEM.ProblemFDEM_e(model) prb.pair(survey) @@ -37,7 +44,7 @@ class FDEM_bDerivTests(unittest.TestCase): j_sx[4,4,4] = 1 j_s = np.r_[Utils.mkvc(j_sx),np.zeros(mesh.nEy+mesh.nEz)] - prb.j_s = j_s + prb.j_s = np.c_[j_s, j_s] f = prb.fields(sigma) self.sigma = sigma @@ -51,13 +58,15 @@ class FDEM_bDerivTests(unittest.TestCase): passed = Tests.checkDerivative(fun, x0, num=3, plotIt=False) self.assertTrue(passed) - # def test_Jtvec(self): - # v = self.survey.nD - # x0 = self.sigma - # def fun(x): - # return self.survey.dpred(x), lambda x: self.prb.Jvec(x0, x) - # passed = Tests.checkDerivative(fun, x0, num=3, plotIt=False) - # self.assertTrue(passed) + def test_Jtvec_adjointTest(self): + v = np.random.rand(self.survey.nD) + w = np.random.rand(self.prb.model.nP) + + m = self.sigma + u = self.prb.fields(m) + vJw = v.dot(self.prb.Jvec(m, w, u=u)) + wJtv = w.dot(self.prb.Jtvec(m, v, u=u)) + self.assertTrue(vJw - wJtv < TOL) if __name__ == '__main__':