diff --git a/SimPEG/EM/Static/DC/ProblemDC.py b/SimPEG/EM/Static/DC/ProblemDC.py index ff3fff0f..a31731a1 100644 --- a/SimPEG/EM/Static/DC/ProblemDC.py +++ b/SimPEG/EM/Static/DC/ProblemDC.py @@ -44,22 +44,47 @@ class BaseDCProblem(BaseEMProblem): df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None) df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False) Jv[src, rx] = rx.evalDeriv(src, self.mesh, f, df_dm_v) - Ainv.clean() return Utils.mkvc(Jv) def Jtvec(self, m, v, f=None): - raise NotImplementedError + if f is None: + f = self.fields(m) + + self.curModel = m + + # Ensure v is a data object. + if not isinstance(v, self.dataPair): + v = self.dataPair(self.survey, v) + + Jtv = np.zeros(m.size) + AT = self.getA().T + ATinv = self.Solver(AT, **self.solverOpts) + + for src in self.survey.srcList: + u_src = f[src, self._solutionType] + for rx in src.rxList: + PTv = rx.evalDeriv(src, self.mesh, f, v[src, rx], adjoint=True) # wrt f, need possibility wrt m + df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None) + df_duT, df_dmT = df_duTFun(src, None, PTv, adjoint=True) + + ATinvdf_duT = ATinv * df_duT + + dA_dmT = self.getADeriv(u_src, ATinvdf_duT, adjoint=True) + dRHS_dmT = self.getRHSDeriv(src, ATinvdf_duT, adjoint=True) + du_dmT = -dA_dmT + dRHS_dmT + Jtv += df_dmT + du_dmT + + return Utils.mkvc(Jtv) def getSourceTerm(self): """ takes concept of source and turns it into a matrix """ """ - Evaluates the sources for a given frequency and puts them in matrix form + Evaluates the sources, and puts them in matrix form - :param float freq: Frequency :rtype: (numpy.ndarray, numpy.ndarray) - :return: s_m, s_e (nE or nF, nSrc) + :return: q (nC or nN, nSrc) """ Srcs = self.survey.srcList @@ -128,8 +153,10 @@ class Problem3D_N(BaseDCProblem): """ Derivative of the right hand side with respect to the model """ - qDeriv = src.evalDeriv(self, adjoint=adjoint) - return qDeriv + # TODO: add qDeriv for RHS depending on m + # qDeriv = src.evalDeriv(self, adjoint=adjoint) + # return qDeriv + return Zero() class Problem3D_CC(BaseDCProblem): @@ -169,7 +196,7 @@ class Problem3D_CC(BaseDCProblem): if adjoint: # if self._makeASymmetric is True: # v = V * v - return D * MfRhoIDeriv(D * v) + return( MfRhoIDeriv( D.T * u ).T) * ( D.T * v) # I think we should deprecate this for DC problem. # if self._makeASymmetric is True: diff --git a/SimPEG/EM/Static/DC/RxDC.py b/SimPEG/EM/Static/DC/RxDC.py index 7692c880..7d50392e 100644 --- a/SimPEG/EM/Static/DC/RxDC.py +++ b/SimPEG/EM/Static/DC/RxDC.py @@ -38,7 +38,13 @@ class BaseRx(SimPEG.Survey.BaseRx): def evalDeriv(self, src, mesh, f, v, adjoint=False): P = self.getP(mesh, self.projGLoc(f)) - return P*v + if not adjoint: + return P*v + elif adjoint: + return P.T*v + + + # DC.Rx.Dipole(locs) class Dipole(BaseRx): diff --git a/tests/em/static/test_DC.py b/tests/em/static/test_DC.py index 7de2ea45..6c777acf 100644 --- a/tests/em/static/test_DC.py +++ b/tests/em/static/test_DC.py @@ -47,31 +47,21 @@ class DCProblemTests(unittest.TestCase): passed = Tests.checkDerivative(derChk, self.m0, plotIt=False) self.assertTrue(passed) - # def test_adjoint(self): - # # Adjoint Test - # u = np.random.rand(self.mesh.nC*self.survey.nSrc) - # v = np.random.rand(self.mesh.nC) - # w = np.random.rand(self.survey.dobs.shape[0]) - # wtJv = w.dot(self.p.Jvec(self.m0, v)) - # vtJtw = v.dot(self.p.Jtvec(self.m0, w)) - # passed = np.abs(wtJv - vtJtw) < 1e-10 - # print 'Adjoint Test', np.abs(wtJv - vtJtw), passed - # self.assertTrue(passed) - - # def test_dataObj(self): - # derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)] - # passed = Tests.checkDerivative(derChk, self.m0, plotIt=False) - # self.assertTrue(passed) - - - # def test_massMatrices(self): - # Gu = np.random.rand(self.mesh.nF) - # def derChk(m): - # self.p.curModel = m - # return [self.p.Msig * Gu, self.p.dMdsig(Gu)] - # passed = Tests.checkDerivative(derChk, self.m0, plotIt=False) - # self.assertTrue(passed) + def test_adjoint(self): + # Adjoint Test + u = np.random.rand(self.mesh.nC*self.survey.nSrc) + v = np.random.rand(self.mesh.nC) + w = np.random.rand(self.survey.dobs.shape[0]) + wtJv = w.dot(self.p.Jvec(self.m0, v)) + vtJtw = v.dot(self.p.Jtvec(self.m0, w)) + passed = np.abs(wtJv - vtJtw) < 1e-10 + print 'Adjoint Test', np.abs(wtJv - vtJtw), passed + self.assertTrue(passed) + def test_dataObj(self): + derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)] + passed = Tests.checkDerivative(derChk, self.m0, plotIt=False) + self.assertTrue(passed) if __name__ == '__main__': unittest.main()