diff --git a/SimPEG/EM/TDEM/FieldsTDEM.py b/SimPEG/EM/TDEM/FieldsTDEM.py index 6129cf01..033dd68f 100644 --- a/SimPEG/EM/TDEM/FieldsTDEM.py +++ b/SimPEG/EM/TDEM/FieldsTDEM.py @@ -86,15 +86,12 @@ class Fields_b(Fields): _, S_e = src.eval(self.survey.prob, self.survey.prob.times[tInd]) bSolution = self[[src],'bSolution',tInd] - if adjoint: - v = self.MeSigmaI.T * v - - _, S_eDeriv_v = src.evalDeriv(self.survey.prob.times[tInd], self, v=v, adjoint=adjoint) + _, S_eDeriv = src.evalDeriv(self.survey.prob.times[tInd], self, adjoint=adjoint) if adjoint is True: - return self.MeSigmaIDeriv(-S_e + self.edgeCurl.T * ( self.MfMui * bSolution ) ).T * v - S_eDeriv_v + return self.MeSigmaIDeriv(-S_e + self.edgeCurl.T * ( self.MfMui * bSolution ) ).T * v - S_eDeriv(self.MeSigmaI.T * v) - return self.MeSigmaIDeriv(-S_e + self.edgeCurl.T * ( self.MfMui * bSolution)) * v - self.MeSigmaI * S_eDeriv_v + return self.MeSigmaIDeriv(-S_e + self.edgeCurl.T * ( self.MfMui * bSolution)) * v - self.MeSigmaI * S_eDeriv(v) def _eDeriv(self, tInd, src, dun_dm_v, v, adjoint=False): if adjoint is True: diff --git a/SimPEG/EM/TDEM/TDEM.py b/SimPEG/EM/TDEM/TDEM.py index 78fd9287..c3fc824d 100644 --- a/SimPEG/EM/TDEM/TDEM.py +++ b/SimPEG/EM/TDEM/TDEM.py @@ -168,7 +168,8 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem): df_duT_v[src, '%sDeriv'%self._fieldType, :] = np.zeros_like(u[src, self._fieldType, :]) for rx in src.rxList: - PT_v[src,'%sDeriv'%rx.projField,:] = rx.evalDeriv(src, self.mesh, self.timeMesh, Utils.mkvc(v[src,rx]), adjoint=True) + PT_v[src,'%sDeriv'%rx.projField,:] = rx.evalDeriv(src, self.mesh, self.timeMesh, Utils.mkvc(v[src,rx]), adjoint=True) # this is += + # PT_v = np.reshape(curPT_v,(len(curPT_v)/self.timeMesh.nN, self.timeMesh.nN), order='F') df_duTFun = getattr(u, '_%sDeriv'%rx.projField, None) @@ -187,13 +188,14 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem): AdiagTinv = None # Do the back-solve through time - for tInd in reversed(range(self.nT)): + for tIndP in reversed(range(self.nT + 1)): + tInd = tIndP - 1 if AdiagTinv is not None and (tInd <= self.nT and self.timeSteps[tInd] != self.timeSteps[tInd+1]): # if the previous timestep is the same --> no need to refactor the matrix AdiagTinv.clean() AdiagTinv = None # refactor if we need to - if AdiagTinv is None: + if AdiagTinv is None and tInd > -1: Adiag = self.getAdiag(tInd) AdiagTinv = self.Solver(Adiag.T, **self.solverOpts) @@ -204,25 +206,46 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem): for isrc, src in enumerate(self.survey.srcList): # solve against df_duT_v if tInd >= self.nT-1: + # last timestep (first to be solved) ATinv_df_duT_v[isrc,:] = AdiagTinv * df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1] - else: + elif tInd > -1: ATinv_df_duT_v[isrc,:] = AdiagTinv * (Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1]) - Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:])) + else: + # AdiagTinv = I + ATinv_df_duT_v[isrc,:] = (Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1]) - Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:])) - un_src = u[src,ftype,tInd+1] - dAT_dm_v = self.getAdiagDeriv(None, un_src, ATinv_df_duT_v[isrc,:], adjoint=True) # cell centered on time mesh + if tInd > -1: + un_src = u[src,ftype,tInd+1] + dAT_dm_v = self.getAdiagDeriv(tInd, un_src, ATinv_df_duT_v[isrc,:], adjoint=True) # cell centered on time mesh + dRHST_dm_v = self.getRHSDeriv(tInd+1, src, ATinv_df_duT_v[isrc,:], adjoint=True) # on nodes of time mesh + JTv = JTv + Utils.mkvc(-dAT_dm_v + dRHST_dm_v) + # else: + # print 'here' - dRHST_dm_v = self.getRHSDeriv(tInd+1, src, ATinv_df_duT_v[isrc,:], adjoint=True) # on nodes of time mesh + # un_src = u[src,ftype,tInd+1] + # # dAT_dm_v = self.getAdiagDeriv(tInd, un_src, ATinv_df_duT_v[isrc,:], adjoint=True) # cell centered on time mesh + # dRHST_dm_v0 = self.getRHSDeriv(tInd+1, src, ATinv_df_duT_v[isrc,:], adjoint=True) # on nodes of time mesh + # dRHST_dm_v1 = self.getInitialFieldsDeriv( Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1]), adjoint=True) + # JTv = JTv + Utils.mkvc(dRHST_dm_v0 + dRHST_dm_v1) + + # print 'here' + # inFields = self.getInitialFieldsDeriv(u[src,ftype,tInd+1], Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1]), adjoint=True) + # # - Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:]), adjoint=True) + # print inFields.shape + # JTv = JTv + inFields # dAsubdiag_dm_v = 0 - JTv = JTv + Utils.mkvc(-dAT_dm_v + dRHST_dm_v) # Missing the 0 step # adding du_dm^T * dF_du^T * P^T vfor time 0 (no dRHS_dm_v at time 0) - for src in self.survey.srcList: - for projField in set(rx.projField): - JTv = JTv + self.getInitialFieldsDeriv(Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,0]), adjoint=True) + # Asubdiag = self.getAsubdiag(0) + # for src in self.survey.srcList: + # for projField in set(rx.projField): + # v = AdiagTinv * (Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,0]) - Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:])) + # JTv = JTv - Utils.mkvc(self.getAdiagDeriv(0, u[src, ftype, tInd], v, adjoint = True)) + # # JTv = JTv + self.getInitialFieldsDeriv(Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,0] - Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:])), adjoint=True) return Utils.mkvc(JTv) @@ -274,7 +297,10 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem): ifieldsDeriv[:,i] = ifieldsDeriv[:,i] + getattr(src, '%sInitialDeriv'%self._fieldType, None)(self,v,adjoint) if adjoint is True: - ifieldsDeriv = ifieldsDeriv.sum() + ifieldsDeriv = np.zeros_like(self.curModel) + print v.shape + # ifieldsDeriv = self.getAdiagDeriv(None, u, v, adjoint) + # ifieldsDeriv = ifieldsDeriv.sum() return ifieldsDeriv @@ -337,6 +363,7 @@ class Problem_b(BaseTDEMProblem): (\mathbf{I} + \mathbf{dt} \mathbf{C} \mathbf{M_{\sigma}^e}^{-1} \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f}) """ + assert tInd >= 0 and tInd < self.nT dt = self.timeSteps[tInd] C = self.mesh.edgeCurl diff --git a/tests/em/tdem/test_TDEM_b_DerivAdjoint.py b/tests/em/tdem/test_TDEM_b_DerivAdjoint.py index a4c58aed..5796228a 100644 --- a/tests/em/tdem/test_TDEM_b_DerivAdjoint.py +++ b/tests/em/tdem/test_TDEM_b_DerivAdjoint.py @@ -23,7 +23,7 @@ def setUp(rxcomp='bz'): mapping = Maps.ExpMap(mesh) * Maps.SurjectVertical1D(mesh) * activeMap rxOffset = 10. - rx = EM.TDEM.Rx(np.array([[rxOffset, 0., -1e-2]]), np.logspace(-4,-3, 20), rxcomp) + rx = EM.TDEM.Rx(np.array([[rxOffset, 0., -1e-2]]), np.logspace(-4,-3, 20), rxcomp) #,] src = EM.TDEM.SurveyTDEM.MagDipole([rx], loc=np.array([0., 0., 0.])) survey = EM.TDEM.Survey([src]) @@ -45,13 +45,14 @@ def setUp(rxcomp='bz'): return prb, m, mesh + class TDEM_DerivTests(unittest.TestCase): # ====== TEST A ========== # - def test_Deriv_Pieces(self): + def test_AderivTest(self): prb, m0, mesh = setUp() - tInd = 0 + tInd = 2 v = np.random.rand(mesh.nF) @@ -64,36 +65,78 @@ class TDEM_DerivTests(unittest.TestCase): return Av, ADeriv_dm - def A_adjointTest(): - print '\n Testing A_adjoint' - m = np.random.rand(prb.mapping.nP) - v = np.random.rand(prb.mesh.nF) - u = np.random.rand(prb.mesh.nF) - prb.curModel = m0 - - tInd = 0 # not actually used - V1 = v.dot(prb.getAdiagDeriv(tInd, u, m)) - V2 = m.dot(prb.getAdiagDeriv(tInd, u, v, adjoint=True)) - passed = np.abs(V1-V2) < TOL * (np.abs(V1) + np.abs(V2))/2. - print 'AdjointTest', V1, V2, passed - self.assertTrue(passed) - print '\n Testing ADeriv' Tests.checkDerivative(AderivTest, m0, plotIt=False, num=4, eps=1e-20) - A_adjointTest() + + def test_A_adjointTest(self): + prb, m0, mesh = setUp() + tInd = 2 + + print '\n Testing A_adjoint' + m = np.random.rand(prb.mapping.nP) + v = np.random.rand(prb.mesh.nF) + u = np.random.rand(prb.mesh.nF) + prb.curModel = m0 + + tInd = 2 # not actually used + V1 = v.dot(prb.getAdiagDeriv(tInd, u, m)) + V2 = m.dot(prb.getAdiagDeriv(tInd, u, v, adjoint=True)) + passed = np.abs(V1-V2) < TOL * (np.abs(V1) + np.abs(V2))/2. + print 'AdjointTest', V1, V2, passed + self.assertTrue(passed) + +# ====== TEST Fields Deriv Pieces ========== # + + def test_eDeriv_m_adjoint(self): + prb, m0, mesh = setUp() + tInd = 0 + + v = np.random.rand(mesh.nF) + + print '\n Testing eDeriv_m Adjoint' + + prb, m0, mesh = setUp() + f = prb.fields(m0) + + m = np.random.rand(prb.mapping.nP) + e = np.random.randn(prb.mesh.nE) + V1 = e.dot(f._eDeriv_m(1, prb.survey.srcList[0], m)) + V2 = m.dot(f._eDeriv_m(1, prb.survey.srcList[0], e, adjoint=True)) + tol = TOL * (np.abs(V1) + np.abs(V2)) / 2. + passed = np.abs(V1-V2) < tol + + print ' ', V1, V2, np.abs(V1-V2), tol, passed + self.assertTrue(passed) + + def test_eDeriv_u_adjoint(self): + print '\n Testing eDeriv_u Adjoint' + + prb, m0, mesh = setUp() + f = prb.fields(m0) + + b = np.random.rand(prb.mesh.nF) + e = np.random.randn(prb.mesh.nE) + V1 = e.dot(f._eDeriv_u(1, prb.survey.srcList[0], b)) + V2 = b.dot(f._eDeriv_u(1, prb.survey.srcList[0], e, adjoint=True)) + tol = TOL * (np.abs(V1) + np.abs(V2)) / 2. + passed = np.abs(V1-V2) < tol + + print ' ', V1, V2, np.abs(V1-V2), tol, passed + self.assertTrue(passed) # ====== TEST Jvec ========== # - def JvecTest(self, rxcomp): - prb, m, mesh = setUp(rxcomp) - - derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(m, mx)] - print '\n' - print 'test_Jvec_%s' %(rxcomp) - Tests.checkDerivative(derChk, m, plotIt=False, num=2, eps=1e-20) - if testDeriv: + + def JvecTest(self, rxcomp): + prb, m, mesh = setUp(rxcomp) + + derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(m, mx)] + print '\n' + 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') @@ -103,33 +146,38 @@ class TDEM_DerivTests(unittest.TestCase): def test_Jvec_b_ey(self): self.JvecTest('ey') + else: + pass + # ====== TEST Jtvec ========== # - def adjointJvecVsJtvecTest(self, rxcomp='bz'): - - print '\nAdjoint Testing Jvec, Jtvec %s' %(rxcomp) - - prb, m0, mesh = setUp(rxcomp) - m = np.random.rand(prb.mapping.nP) - d = np.random.randn(prb.survey.nD) - V1 = d.dot(prb.Jvec(m0, m)) - V2 = m.dot(prb.Jtvec(m0, d)) - tol = TOL * (np.abs(V1) + np.abs(V2)) / 2. - passed = np.abs(V1-V2) < tol - - print ' ', V1, V2, np.abs(V1-V2), tol, passed - self.assertTrue(passed) - if testAdjoint: - def test_Jvec_b_bx(self): - self.adjointJvecVsJtvecTest('bx') - def test_Jvec_b_bz(self): - self.adjointJvecVsJtvecTest('bz') + def JvecVsJtvecTest(self, rxcomp='bz'): + + print '\nAdjoint Testing Jvec, Jtvec %s' %(rxcomp) + + prb, m0, mesh = setUp(rxcomp) + m = np.random.rand(prb.mapping.nP) + d = np.random.randn(prb.survey.nD) + V1 = d.dot(prb.Jvec(m0, m)) + V2 = m.dot(prb.Jtvec(m0, d)) + tol = TOL * (np.abs(V1) + np.abs(V2)) / 2. + passed = np.abs(V1-V2) < tol + + print ' ', V1, V2, np.abs(V1-V2), tol, passed + self.assertTrue(passed) + + def test_Jvec_adjoint_b_bx(self): + self.JvecVsJtvecTest('bx') + + def test_Jvec_adjoint_b_bz(self): + self.JvecVsJtvecTest('bz') + + def test_Jvec_adjoint_b_ey(self): + self.JvecVsJtvecTest('ey') - def test_Jvec_b_ey(self): - self.adjointJvecVsJtvecTest('ey') if __name__ == '__main__': diff --git a/tests/em/tdem/test_TDEM_b_MultiSrc_DerivAdjoint.py b/tests/em/tdem/test_TDEM_b_MultiSrc_DerivAdjoint.py index 7d0fc724..2d84de42 100644 --- a/tests/em/tdem/test_TDEM_b_MultiSrc_DerivAdjoint.py +++ b/tests/em/tdem/test_TDEM_b_MultiSrc_DerivAdjoint.py @@ -6,11 +6,9 @@ plotIt = False testDeriv = True testAdjoint = True -TOL = 1e-6 +TOL = 1e-5 -class TDEM_bDerivTests(unittest.TestCase): - - def setUp(self): +def setUp(self, rxcomp='bz'): cs = 5. ncx = 20 @@ -25,46 +23,58 @@ class TDEM_bDerivTests(unittest.TestCase): mapping = Maps.ExpMap(mesh) * Maps.SurjectVertical1D(mesh) * activeMap rxOffset = 40. - rx = EM.TDEM.Rx(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20), 'bz') + rx = EM.TDEM.Rx(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20), rxcomp) src = EM.TDEM.SurveyTDEM.MagDipole( [rx], loc=np.array([0., 0., 0.])) - rx2 = EM.TDEM.Rx(np.array([[rxOffset-10, 0., 0.]]), np.logspace(-5,-4, 25), 'bz') + rx2 = EM.TDEM.Rx(np.array([[rxOffset-10, 0., 0.]]), np.logspace(-5,-4, 25), rxcomp) src2 = EM.TDEM.SurveyTDEM.MagDipole( [rx2], loc=np.array([0., 0., 0.])) survey = EM.TDEM.Survey([src,src2]) - self.prb = EM.TDEM.Problem_b(mesh, mapping=mapping) - # self.prb.timeSteps = [1e-5] - self.prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)] - # self.prb.timeSteps = [(1e-05, 100)] + prb = EM.TDEM.Problem_b(mesh, mapping=mapping) + # prb.timeSteps = [1e-5] + prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)] + # prb.timeSteps = [(1e-05, 100)] try: from pymatsolver import MumpsSolver - self.prb.Solver = MumpsSolver + prb.Solver = MumpsSolver except ImportError, e: - self.prb.Solver = SolverLU + prb.Solver = SolverLU - self.m = np.log(1e-1)*np.ones(self.prb.mapping.nP) + 1e-2*np.random.randn(self.prb.mapping.nP) + m = np.log(1e-1)*np.ones(prb.mapping.nP) + 1e-2*np.random.randn(prb.mapping.nP) + + prb.pair(survey) + + return mesh, prb, m + +class TDEM_bDerivTests(unittest.TestCase): - self.prb.pair(survey) - self.mesh = mesh if testDeriv: - def test_Deriv_J(self): + def Deriv_J(self, rxcomp='bz'): + + mesh, prb, m0 = setUp(rxcomp) - prb = self.prb prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)] - mesh = self.mesh - derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(self.m, mx)] + derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(m0, mx)] print '\n' - print 'test_Deriv_J' - Tests.checkDerivative(derChk, self.m, plotIt=False, num=3, eps=1e-20) + print 'test_Deriv_J %s'%rxcomp + Tests.checkDerivative(derChk, m0, plotIt=False, num=3, eps=1e-20) + + def test_Jvec_bx(self): + self.Deriv_J('bx') + + def test_Jvec_bz(self): + self.Deriv_J('bz') + + def test_Jvec_ey(self): + self.Deriv_J('ey') if testAdjoint: - def test_adjointJvecVsJtvec(self): - mesh = self.mesh - prb = self.prb - m0 = self.m + def adjointJvecVsJtvec(self, rxcomp='bz'): + print ' \n Testing Adjoint %s' %rxcomp + mesh, prb, m0 = setUp(rxcomp) m = np.random.rand(prb.mapping.nP) d = np.random.rand(prb.survey.nD) @@ -77,6 +87,15 @@ class TDEM_bDerivTests(unittest.TestCase): print ' ', V1, V2, np.abs(V1-V2), tol, passed self.assertTrue(passed) + def test_JvecVsJtvec_bx(self): + self.adjointJvecVsJtvec('bx') + + def test_JvecVsJtvec_bz(self): + self.adjointJvecVsJtvec('bz') + + def test_JvecVsJtvec_ey(self): + self.adjointJvecVsJtvec('ey') + if __name__ == '__main__':