diff --git a/docs/api_TDEM_derivation.rst b/docs/api_TDEM_derivation.rst index 0d8bfcea..af3fc2fc 100644 --- a/docs/api_TDEM_derivation.rst +++ b/docs/api_TDEM_derivation.rst @@ -307,38 +307,21 @@ Multiplying \\(\\mathbf{J}^\\top\\) onto a vector can be broken into three steps \end{array} \right] -For the last time-step \\(t=N\\): +For the all time-steps (going backwards in time): + .. math:: - \begin{align} - \frac{1}{\delta t} \MfMui \vec{y}_{b}^{(N)} + \MfMui \dcurl \vec{y}_{e}^{(N)} = \vec{p}_b^{(N)} \\ - \dcurl^\top \MfMui \vec{y}_b^{(N)} - \MeSig \vec{y}_e^{(N)} = \vec{p}_e^{(N)} - \end{align} + A \vec{y}^{(t)} + B \vec{y}^{(t+1)} = \vec{p}^{(t)} .. math:: \begin{align} - \left( \MfMui \dcurl \MeSig^{-1} \dcurl^\top \MfMui + \frac{1}{\delta t} \MfMui \right) \vec{y}_{b}^{(N)} = \MfMui \dcurl \MeSig^{-1} \vec{p}_e^{(N)} + \vec{p}_b^{(N)} \\ - \vec{y}_e^{(N)} = \MeSig^{-1} \dcurl^\top \MfMui \vec{y}_b^{(N)} - \MeSig^{-1} \vec{p}_e^{(N)} - \end{align} - -For the rest of the time-steps (going backwards in time) - - -.. math:: - - A \vec{y}^{(t-1)} + B \vec{y}^{(t)} = \vec{p}^{(t-1)} - - -.. math:: - - \begin{align} - \frac{1}{\delta t} \MfMui\vec{y}_{b}^{(t-1)} + \MfMui\dcurl \vec{y}_{e}^{(t-1)} - - \frac{1}{\delta t} \MfMui \vec{y}_{b}^{(t)} - = \vec{p}_b^{(t-1)} \\ - \dcurl^\top \MfMui \vec{y}_b^{(t-1)} - \MeSig \vec{y}_e^{(t-1)} = \vec{p}_e^{(t-1)} + \frac{1}{\delta t} \MfMui\vec{y}_{b}^{(t)} + \MfMui\dcurl \vec{y}_{e}^{(t)} + - \frac{1}{\delta t} \MfMui \vec{y}_{b}^{(t+1)} + = \vec{p}_b^{(t)} \\ + \dcurl^\top \MfMui \vec{y}_b^{(t)} - \MeSig \vec{y}_e^{(t)} = \vec{p}_e^{(t)} \end{align} and @@ -346,8 +329,13 @@ and .. math:: \begin{align} - \left( \MfMui \dcurl \MeSig^{-1} \dcurl^\top \MfMui + \frac{1}{\delta t} \MfMui \right) \vec{y}_{b}^{(t-1)} = - \frac{1}{\delta t} \MfMui \vec{y}_b^{(t)} - + \MfMui \dcurl \MeSig^{-1} \vec{p}_e^{(t-1)} + \vec{p}_b^{(t-1)} \\ - \vec{y}_e^{(t-1)} = \MeSig^{-1} \dcurl^\top \MfMui \vec{y}_b^{(t-1)} - \MeSig^{-1} \vec{p}_e^{(t-1)} + \left( \MfMui \dcurl \MeSig^{-1} \dcurl^\top \MfMui + \frac{1}{\delta t} \MfMui \right) \vec{y}_{b}^{(t)} = + \frac{1}{\delta t} \MfMui \vec{y}_b^{(t+1)} + + \MfMui \dcurl \MeSig^{-1} \vec{p}_e^{(t)} + \vec{p}_b^{(t)} \\ + \vec{y}_e^{(t)} = \MeSig^{-1} \dcurl^\top \MfMui \vec{y}_b^{(t)} - \MeSig^{-1} \vec{p}_e^{(t)} \end{align} + + +.. note:: + + For the last time step, \\\(t=N\\\), the term: \\\(\\frac{1}{\\delta t} \\MfMui \\vec{y}_b^{(N+1)}\\\) is zero. diff --git a/simpegEM/TDEM/BaseTDEM.py b/simpegEM/TDEM/BaseTDEM.py index 12911847..cacf6eaa 100644 --- a/simpegEM/TDEM/BaseTDEM.py +++ b/simpegEM/TDEM/BaseTDEM.py @@ -57,8 +57,7 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem): return F def adjoint(self, m, RHS, CalcFields, F=None): - if F is None: - F = FieldsTDEM(self.mesh, self.survey) + F = F or FieldsTDEM(self.mesh, self.survey) dtFact = None for tInd, dt in reversed(list(enumerate(self.timeSteps))): @@ -66,13 +65,12 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem): dtFact = dt A = self.getA(tInd) # print 'Factoring... (dt = ' + str(dt) + ')' - Asolve = Solver(A, options=self.solverOpts) + Asolve = self.Solver(A, **self.solverOpts) # print 'Done' rhs = RHS(tInd, F) sol = Asolve.solve(rhs) if sol.ndim == 1: sol.shape = (sol.size,1) - newFields = CalcFields(sol, self.solType, tInd) - F[:,:,tInd] = newFields + F[:,:,tInd+1] = CalcFields(sol, self.solType, tInd) return F diff --git a/simpegEM/TDEM/SurveyTDEM.py b/simpegEM/TDEM/SurveyTDEM.py index 3bc5fde6..722846f5 100644 --- a/simpegEM/TDEM/SurveyTDEM.py +++ b/simpegEM/TDEM/SurveyTDEM.py @@ -39,8 +39,7 @@ class RxTDEM(Survey.BaseTimeRx): if not adjoint: return P * Utils.mkvc(v[tx, self.projField, :]) elif adjoint: - Ptv = P.T * v[tx, self] - return Ptv + return P.T * v[tx, self] class FieldsTDEM(Survey.TimeFields): diff --git a/simpegEM/TDEM/TDEM_b.py b/simpegEM/TDEM/TDEM_b.py index e60b9ec9..95304962 100644 --- a/simpegEM/TDEM/TDEM_b.py +++ b/simpegEM/TDEM/TDEM_b.py @@ -64,7 +64,7 @@ class ProblemTDEM_b(BaseTDEMProblem): p = self.survey.projectFieldsDeriv(u, v=v, adjoint=True) y = self.solveAht(m, p) w = self.Gtvec(m, y, u) - return w + return - w def Gvec(self, m, vec, u=None): """ @@ -138,18 +138,27 @@ class ProblemTDEM_b(BaseTDEMProblem): def solveAht(self, m, p): def AhtRHS(tInd, u): - rhs = self.MfMui*self.mesh.edgeCurl*self.MeSigmaI*p[:,'e',tInd] + p[:,'b',tInd] + nTx, nF = self.survey.nTx, self.mesh.nF + rhs = np.zeros(nF if nTx == 1 else (nF, nTx)) + + if 'e' in p: + rhs += self.MfMui*self.mesh.edgeCurl*self.MeSigmaI*p[:,'e',tInd] + if 'b' in p: + rhs += p[:,'b',tInd] + if tInd == self.nT-1: return rhs - dt = self.timeSteps[tInd+1] + dt = self.timeSteps[tInd] return rhs + 1.0/dt*self.MfMui*u[:,'b',tInd+1] def AhtCalcFields(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] - return {'b':b, 'e':e} + y_b = mkvc(y_b) + y_e = self.MeSigmaI*self.mesh.edgeCurl.T*self.MfMui*y_b + if 'e' in p: + y_e += - self.MeSigmaI*p[:,'e',tInd] + return {'b':y_b, 'e':y_e} self.curModel = m return self.adjoint(m, AhtRHS, AhtCalcFields) diff --git a/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py b/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py index 9f505ce3..7b32f484 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 @@ -223,48 +223,48 @@ class TDEM_bDerivTests(unittest.TestCase): # V2 = f1.tovec().dot(prb.AhtVec(sigma, f2).tovec()) # self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6) - # def test_solveAhtVsAhtVec(self): - # prb = self.prb - # mesh = self.mesh - # sigma = np.random.rand(prb.mapping.nP) + def test_solveAhtVsAhtVec(self): + prb = self.prb + mesh = self.mesh + sigma = np.random.rand(prb.mapping.nP) - # f1 = EM.TDEM.FieldsTDEM(mesh, 1, prb.nT, 'b') - # for i in range(prb.nT): - # f1.set_b(np.random.rand(mesh.nF, 1), i) - # f1.set_e(np.random.rand(mesh.nE, 1), i) + f1 = EM.TDEM.FieldsTDEM(mesh,prb.survey) + for i in range(prb.nT): + 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) + f2 = prb.solveAht(sigma, f1) + f3 = prb.AhtVec(sigma, f2) - # if plotIt: - # import matplotlib.pyplot as plt - # plt.plot(f3.tovec()) - # plt.plot(f1.tovec()) - # plt.show() - # V1 = np.linalg.norm(f3.tovec()-f1.tovec()) - # V2 = np.linalg.norm(f1.tovec()) - # print V1, V2 - # print 'I am gunna fail this one: boo. :(' - # self.assertLess(V1/V2, 1e-6) + if plotIt: + import matplotlib.pyplot as plt + plt.plot(f3.tovec()) + plt.plot(f1.tovec()) + plt.show() + V1 = np.linalg.norm(f3.tovec()-f1.tovec()) + V2 = np.linalg.norm(f1.tovec()) + print V1, V2 + 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 @@ -274,12 +274,12 @@ class TDEM_bDerivTests(unittest.TestCase): # sigma = np.random.rand(prb.mapping.nP) # u = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey) - # for i in range(prb.nT): + # 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(prb.nT): + # 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)