diff --git a/docs/api_TDEM_derivation.rst b/docs/api_TDEM_derivation.rst index 364070e5..0d8bfcea 100644 --- a/docs/api_TDEM_derivation.rst +++ b/docs/api_TDEM_derivation.rst @@ -255,25 +255,8 @@ Multiplying **J** onto a vector can be broken into three steps \vec{p}_e^{(n)} = - \diag{\e^{(n)}} \Ace \diag{V} m \end{align} -First time step -.. math:: - - \begin{align} - \frac{1}{\delta t} \MfMui \vec{y}_{b}^{(1)} + \MfMui \dcurl \vec{y}_{e}^{(1)} = \vec{p}_b^{(1)} \\ - \dcurl^\top \MfMui \vec{y}_b^{(1)} - \MeSig \vec{y}_e^{(1)} = \vec{p}_e^{(1)} - \end{align} - - -.. math:: - - \begin{align} - \left( \MfMui \dcurl \MeSig^{-1} \dcurl^\top \MfMui + \frac{1}{\delta t} \MfMui \right) \vec{y}_{b}^{(1)} = \MfMui \dcurl \MeSig^{-1} \vec{p}_e^{(1)} + \vec{p}_b^{(1)} \\ - \vec{y}_e^{(1)} = \MeSig^{-1} \dcurl^\top \MfMui \vec{y}_b^{(1)} - \MeSig^{-1} \vec{p}_e^{(1)} - \end{align} - - -Remaining time steps: +For all time steps: .. math:: @@ -295,6 +278,11 @@ and \vec{y}_e^{(t+1)} = \MeSig^{-1} \dcurl^\top \MfMui \vec{y}_b^{(t+1)} - \MeSig^{-1} \vec{p}_e^{(t+1)} \end{align} +.. note:: + + For the first time step, \\\(t=0\\\), the term: \\\(\\frac{1}{\\delta t} \\MfMui \\vec{y}_b^{(0)}\\\) is zero. + + Implementing **J** transpose times a vector diff --git a/simpegEM/TDEM/TDEM_b.py b/simpegEM/TDEM/TDEM_b.py index 2999f01e..e60b9ec9 100644 --- a/simpegEM/TDEM/TDEM_b.py +++ b/simpegEM/TDEM/TDEM_b.py @@ -52,7 +52,7 @@ class ProblemTDEM_b(BaseTDEMProblem): p = self.Gvec(m, v, u) y = self.solveAh(m, p) Jv = self.survey.projectFieldsDeriv(u, v=y) - return mkvc(Jv) + return - mkvc(Jv) def Jtvec(self, m, v, u=None): if u is None: @@ -117,19 +117,20 @@ class ProblemTDEM_b(BaseTDEMProblem): return p def solveAh(self, m, p): - def AhRHS(tInd, u): + + def AhRHS(tInd, y): rhs = self.MfMui*self.mesh.edgeCurl*self.MeSigmaI*p[:,'e',tInd+1] + p[:,'b',tInd+1] if tInd == 0: return rhs dt = self.timeSteps[tInd] - return rhs + 1.0/dt*self.MfMui*u[:,'b',tInd] + return rhs + 1.0/dt*self.MfMui*y[:,'b',tInd] def AhCalcFields(sol, solType, tInd): - b = sol + y_b = sol if self.survey.nTx == 1: - b = mkvc(b) - e = self.MeSigmaI*self.mesh.edgeCurl.T*self.MfMui*b - self.MeSigmaI*p[:,'e',tInd+1] - return {'b':b, 'e':e} + y_b = mkvc(y_b) + y_e = self.MeSigmaI*self.mesh.edgeCurl.T*self.MfMui*y_b - self.MeSigmaI*p[:,'e',tInd+1] + return {'b':y_b, 'e':y_e} self.curModel = m return self.forward(m, AhRHS, AhCalcFields) diff --git a/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py b/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py index a7b8ea0f..9f505ce3 100644 --- a/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py +++ b/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py @@ -168,22 +168,22 @@ class TDEM_bDerivTests(unittest.TestCase): # passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20) # self.assertTrue(passed) - # def test_Deriv_J(self): + def test_Deriv_J(self): - # prb = self.prb - # prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)] - # mesh = self.mesh - # sigma = self.sigma + prb = self.prb + prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)] + mesh = self.mesh + sigma = self.sigma - # # d_sig = 0.8*sigma #np.random.rand(mesh.nCz) - # d_sig = 10*np.random.rand(prb.mapping.nP) + # d_sig = 0.8*sigma #np.random.rand(mesh.nCz) + d_sig = 10*np.random.rand(prb.mapping.nP) - # derChk = lambda m: [prb.survey.dpred(m), lambda mx: -prb.Jvec(sigma, mx)] - # print '\n' - # print 'test_Deriv_J' - # passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=4, eps=1e-20) - # self.assertTrue(passed) + derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(sigma, mx)] + print '\n' + print 'test_Deriv_J' + passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=4, eps=1e-20) + self.assertTrue(passed) # def test_projectAdjoint(self): # prb = self.prb @@ -247,24 +247,24 @@ class TDEM_bDerivTests(unittest.TestCase): # 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 + # 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) + # 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) + # 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()) - self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6) + # V1 = f2.tovec().dot(prb.solveAh(sigma, f1).tovec()) + # V2 = f1.tovec().dot(prb.solveAht(sigma, f2).tovec()) + # self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6) # def test_adjointGvecVsGtvec(self): # mesh = self.mesh