diff --git a/SimPEG/EM/TDEM/FieldsTDEM.py b/SimPEG/EM/TDEM/FieldsTDEM.py index a15d93c0..c5f5d2d7 100644 --- a/SimPEG/EM/TDEM/FieldsTDEM.py +++ b/SimPEG/EM/TDEM/FieldsTDEM.py @@ -31,6 +31,8 @@ class Fields(SimPEG.Problem.TimeFields): knownFields = {} dtype = float + + class Fields_Derivs(Fields): knownFields = { 'bDeriv': 'F', @@ -65,7 +67,7 @@ class Fields_b(Fields): def _bDeriv(self, tInd, src, dun_dm_v, v, adjoint=False): if adjoint is True: - raise NotImplementedError + return self._bDeriv_u(tInd, src, v, adjoint), self._bDeriv_m(tInd, src, v, adjoint) return self._bDeriv_u(tInd, src, dun_dm_v) + self._bDeriv_m(tInd, src, v) def _e(self, bSolution, srcList, tInd): diff --git a/SimPEG/EM/TDEM/SurveyTDEM.py b/SimPEG/EM/TDEM/SurveyTDEM.py index 79041121..90f26d55 100644 --- a/SimPEG/EM/TDEM/SurveyTDEM.py +++ b/SimPEG/EM/TDEM/SurveyTDEM.py @@ -252,24 +252,25 @@ class Survey(SimPEG.Survey.BaseSurvey): return data def evalDeriv(self, u, v=None, adjoint=False): - assert v is not None, 'v to multiply must be provided.' + raise Exception('Use Receivers to project fields deriv.') + # assert v is not None, 'v to multiply must be provided.' - if not adjoint: - data = SimPEG.Survey.Data(self) - for src in self.srcList: - for rx in src.rxList: - data[src, rx] = rx.evalDeriv(src, self.mesh, self.prob.timeMesh, u, v) - return data - else: - f = FieldsTDEM(self.mesh, self) - for src in self.srcList: - for rx in src.rxList: - Ptv = rx.evalDeriv(src, self.mesh, self.prob.timeMesh, u, v, adjoint=True) - Ptv = Ptv.reshape((-1, self.prob.timeMesh.nN), order='F') - if rx.projField not in f: # first time we are projecting - f[src, rx.projField, :] = Ptv - else: # there are already fields, so let's add to them! - f[src, rx.projField, :] += Ptv - return f + # if not adjoint: + # data = SimPEG.Survey.Data(self) + # for src in self.srcList: + # for rx in src.rxList: + # data[src, rx] = rx.evalDeriv(src, self.mesh, self.prob.timeMesh, u, v) + # return data + # else: + # f = FieldsTDEM(self.mesh, self) + # for src in self.srcList: + # for rx in src.rxList: + # Ptv = rx.evalDeriv(src, self.mesh, self.prob.timeMesh, u, v, adjoint=True) + # Ptv = Ptv.reshape((-1, self.prob.timeMesh.nN), order='F') + # if rx.projField not in f: # first time we are projecting + # f[src, rx.projField, :] = Ptv + # else: # there are already fields, so let's add to them! + # f[src, rx.projField, :] += Ptv + # return f diff --git a/SimPEG/EM/TDEM/TDEM.py b/SimPEG/EM/TDEM/TDEM.py index ff4f8cb8..5984b898 100644 --- a/SimPEG/EM/TDEM/TDEM.py +++ b/SimPEG/EM/TDEM/TDEM.py @@ -16,8 +16,6 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem): Problem.BaseTimeProblem.__init__(self, mesh, mapping=mapping, **kwargs) - # _FieldsForward_pair = FieldsTDEM #: used for the forward calculation only - def fields(self, m): """ Solve the forward problem for the fields. @@ -49,12 +47,15 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem): if self.verbose: print 'Done' rhs = self.getRHS(tInd, F) + if self.verbose: print ' Solving... (tInd = %d)'%tInd sol = Ainv * rhs if self.verbose: print ' Done...' + if sol.ndim == 1: sol.shape = (sol.size,1) F[:,self._fieldType+'Solution',tInd+1] = sol + Ainv.clean() return F @@ -105,7 +106,60 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem): Ainv.clean() return Utils.mkvc(Jv) - + + def Jtvec(self, m, v, u=None): + + if u is None: + u = self.fields(m) + + self.curModel = m + ftype = self._fieldType + 'Solution' # the thing we solved for + + # Ensure v is a data object. + if not isinstance(v, self.dataPair): + v = self.dataPair(self.survey, v) + + # TODO: make this general + if self._fieldType is 'b': + dun_dmT_v = np.zeros((len(m), self.mesh.nF)) + + # df_dm_v = Fields_Derivs(self.mesh, self.survey) + JTv = np.zeros(m.size) + PT_v = Fields_Derivs(self.mesh, self.survey) #PT_v is a fields object + + for src in self.survey.srcList: + for rx in src.rxList: + PT_v[src,'%sDeriv'%rx.projField, :] = rx.evalDeriv(src, self.mesh, self.timeMesh, v, adjoint = True) # All the fields for a given src, reciever. + + + ATinv = None + + for tInd, dt in enumerate(reversed(list(self.timeSteps))): + if ATinv is not None and (tInd < self.nT and dt != self.timeSteps[tInd - 1]):# keep factors if dt is the same as previous step b/c A will be the same + ATinv.clean() + ATinv = None + + if ATinv is None: + A = self.getA(tInd) + ATinv = self.Solver(A.T, **self.solverOpts) + + for i, src in enumerate(self.survey.srcList): + + u_src = u[src,ftype,tInd] # fields for this source at tInd + + for rx in src.rxList: + + df_duTFun = getattr(u, '_%sDeriv'%rx.projField, None) + df_duT_v, df_dmT_v = df_duTFun(tInd, src, None, PT_v[src,'%sDeriv'%rx.projField,tInd-1], adjoint=True) + + ATinv_df_duT_v = ATinv * df_duT_v + rhsT_v = self.getJRHS(tInd, src, u_src, ATinv_df_duT_v, dun_dmT_v[:,i], adjoint = True) + + JTv += rhsT_v + df_dmT_v + + return Utils.mkvc(JTv) + + def getJRHS(self, tInd, src, u, v, dbn_dm_v, adjoint = False): @@ -115,9 +169,7 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem): b = - dA_dm + dRHS_dm return b - - def Jtvec(self, m, v, u=None): - raise NotImplementedError + def getSourceTerm(self, tInd): @@ -245,7 +297,7 @@ class Problem_b(BaseTDEMProblem): if adjoint: if self._makeASymmetric is True: v = MfMui * v - return MfMui.T * ( C * ( MeSigmaIDeriv.T * ( C.T * v ) ) ) + return MeSigmaIDeriv(C.T * ( MfMui * u )).T * ( C.T * v ) ADeriv = ( C * ( MeSigmaIDeriv(C.T * ( MfMui * u )) * v ) ) if self._makeASymmetric is True: @@ -279,27 +331,24 @@ class Problem_b(BaseTDEMProblem): MfMui = self.MfMui _, S_e = src.eval(tInd+1, self) # I think this is tInd+1 ? - S_mDeriv_v, S_eDeriv_v = src.evalDeriv(self.times[tInd+1], self, v=v, adjoint=adjoint) # I think this is tInd+1 ? - - # B_n = np.c_[[F[src,'b',tInd] for src in self.survey.srcList]].T - # if B_n.shape[0] is not 1: - # raise NotImplementedError('getRHS not implemented for this shape of B_n') + S_mDeriv, S_eDeriv = src.evalDeriv(self.times[tInd+1], self, adjoint=adjoint) # I think this is tInd+1 ? if adjoint: - raise NotImplementedError + if self._makeASymmetric is True: + v = self.MfMui * v + if isinstance(S_e, Utils.Zero): + MeSigmaIDerivT_v = Utils.Zero() + else: + MeSigmaIDerivT_v = MeSigmaIDeriv(S_e).T * v + RHSDeriv = MeSigmaIDerivT_v + S_eDeriv( MeSigmaI.T * ( C.T * v ) ) + S_mDeriv(v) + dbn_dm_v / dt #this will be given the transposed version + return RHSDeriv - - if isinstance(S_e,Utils.Zero): + if isinstance(S_e, Utils.Zero): MeSigmaIDeriv_v = Utils.Zero() else: MeSigmaIDeriv_v = MeSigmaIDeriv(S_e) * v - # if isinstance(S_eDeriv, Utils.Zero): - # MeSigmaI_S_eDeriv_v = Utils.Zero() - # else: - # MeSigmaI_S_eDeriv_v = MeSigmaI * S_eDeriv(v) - - RHSDeriv = (C * (MeSigmaIDeriv_v + MeSigmaI * S_eDeriv_v) + S_mDeriv_v) + dbn_dm_v / dt + RHSDeriv = (C * (MeSigmaIDeriv_v + MeSigmaI * S_eDeriv(v) + S_mDeriv(v))) + dbn_dm_v / dt if self._makeASymmetric is True: return self.MfMui.T * RHSDeriv diff --git a/tests/em/tdem/test_TDEM_b_DerivAdjoint.py b/tests/em/tdem/test_TDEM_b_DerivAdjoint.py index cdf10e05..6c534159 100644 --- a/tests/em/tdem/test_TDEM_b_DerivAdjoint.py +++ b/tests/em/tdem/test_TDEM_b_DerivAdjoint.py @@ -3,6 +3,10 @@ from SimPEG import * from SimPEG import EM plotIt = False + +testDeriv = False +testAdjoint = True + tol = 1e-6 def setUp(rxcomp='bz'): @@ -70,132 +74,29 @@ class TDEM_bDerivTests(unittest.TestCase): print 'test_Jvec_%s' %(rxcomp) Tests.checkDerivative(derChk, m, plotIt=False, num=2, eps=1e-20) - def test_Jvec_b_bx(self): - self.JvecTest('bx') + if testDeriv: + def test_Jvec_b_bx(self): + self.JvecTest('bx') - def test_Jvec_b_bz(self): - self.JvecTest('bz') + def test_Jvec_b_bz(self): + self.JvecTest('bz') - def test_Jvec_b_ey(self): - self.JvecTest('ey') - + def test_Jvec_b_ey(self): + self.JvecTest('ey') + - # def test_projectAdjoint(self): - # prb = self.prb - # survey = prb.survey - # mesh = self.mesh + if testAdjoint: + def test_adjointJvecVsJtvec(self): + prb, m0, mesh = setUp() - # # Generate random fields and data - # f = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey) - # for i in range(prb.nT): - # f[:,'b',i] = np.random.rand(mesh.nF, 1) - # f[:,'e',i] = np.random.rand(mesh.nE, 1) - # d_vec = np.random.rand(survey.nD) - # d = Survey.Data(survey,v=d_vec) - - # # Check that d.T*Q*f = f.T*Q.T*d - # V1 = d_vec.dot(survey.evalDeriv(None, v=f).tovec()) - # V2 = f.tovec().dot(survey.evalDeriv(None, v=d, adjoint=True).tovec()) + m = np.random.rand(prb.mapping.nP) + d = np.random.rand(prb.survey.nD) - # self.assertTrue((V1-V2)/np.abs(V1) < tol) - - # def test_adjointAhVsAht(self): - # prb = self.prb - # mesh = self.mesh - # sigma = self.sigma - - # f1 = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey) - # for i in range(1,prb.nT+1): - # f1[:,'b',i] = np.random.rand(mesh.nF, 1) - # f1[:,'e',i] = np.random.rand(mesh.nE, 1) - - # f2 = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey) - # for i in range(1,prb.nT+1): - # f2[:,'b',i] = np.random.rand(mesh.nF, 1) - # f2[:,'e',i] = np.random.rand(mesh.nE, 1) - - # V1 = f2.tovec().dot(prb._AhVec(sigma, f1).tovec()) - # V2 = f1.tovec().dot(prb._AhtVec(sigma, f2).tovec()) - # self.assertTrue(np.abs(V1-V2)/np.abs(V1) < tol) - - # # def test_solveAhtVsAhtVec(self): - # # prb = self.prb - # # mesh = self.mesh - # # sigma = np.random.rand(prb.mapping.nP) - - # # f1 = EM.TDEM.FieldsTDEM(mesh,prb.survey) - # # for i in range(1,prb.nT+1): - # # f1[:,'b',i] = np.random.rand(mesh.nF, 1) - # # f1[:,'e',i] = np.random.rand(mesh.nE, 1) - - # # f2 = prb.solveAht(sigma, f1) - # # f3 = prb._AhtVec(sigma, f2) - - # # if True: - # # import matplotlib.pyplot as plt - # # plt.plot(f3.tovec(),'b') - # # plt.plot(f1.tovec(),'r') - # # plt.show() - # # V1 = np.linalg.norm(f3.tovec()-f1.tovec()) - # # V2 = np.linalg.norm(f1.tovec()) - # # print 'AhtVsAhtVec', V1, V2, f1.tovec() - # # print 'I am gunna fail this one: boo. :(' - # # self.assertLess(V1/V2, 1e-6) - - # # def test_adjointsolveAhVssolveAht(self): - # # prb = self.prb - # # mesh = self.mesh - # # sigma = self.sigma - - # # f1 = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey) - # # for i in range(1,prb.nT+1): - # # f1[:,'b',i] = np.random.rand(mesh.nF, 1) - # # f1[:,'e',i] = np.random.rand(mesh.nE, 1) - - # # f2 = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey) - # # for i in range(1,prb.nT+1): - # # f2[:,'b',i] = np.random.rand(mesh.nF, 1) - # # f2[:,'e',i] = np.random.rand(mesh.nE, 1) - - # # V1 = f2.tovec().dot(prb.solveAh(sigma, f1).tovec()) - # # V2 = f1.tovec().dot(prb.solveAht(sigma, f2).tovec()) - # # print V1, V2 - # # self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6) - - # def test_adjointGvecVsGtvec(self): - # mesh = self.mesh - # prb = self.prb - - # m = np.random.rand(prb.mapping.nP) - # sigma = np.random.rand(prb.mapping.nP) - - # u = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey) - # for i in range(1,prb.nT+1): - # u[:,'b',i] = np.random.rand(mesh.nF, 1) - # u[:,'e',i] = np.random.rand(mesh.nE, 1) - - # v = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey) - # for i in range(1,prb.nT+1): - # v[:,'b',i] = np.random.rand(mesh.nF, 1) - # v[:,'e',i] = np.random.rand(mesh.nE, 1) - - # V1 = m.dot(prb.Gtvec(sigma, v, u)) - # V2 = v.tovec().dot(prb.Gvec(sigma, m, u).tovec()) - # self.assertTrue(np.abs(V1-V2)/np.abs(V1) < tol) - - # def test_adjointJvecVsJtvec(self): - # mesh = self.mesh - # prb = self.prb - # sigma = self.sigma - - # m = np.random.rand(prb.mapping.nP) - # d = np.random.rand(prb.survey.nD) - - # V1 = d.dot(prb.Jvec(sigma, m)) - # V2 = m.dot(prb.Jtvec(sigma, d)) - # passed = np.abs(V1-V2)/np.abs(V1) < tol - # print 'AdjointTest', V1, V2, passed - # self.assertTrue(passed) + V1 = d.dot(prb.Jvec(m0, m)) + V2 = m.dot(prb.Jtvec(m0, d)) + passed = np.abs(V1-V2)/np.abs(V1) < tol + print 'AdjointTest', V1, V2, passed + self.assertTrue(passed)