From 2d8bbdce451afd122c57a43d14f72c49ebb1599a Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sun, 20 Mar 2016 11:42:31 -0700 Subject: [PATCH] working to debug JTv for e-formulation... still some work to do --- SimPEG/EM/FDEM/FDEM.py | 72 +++++++++++------------ SimPEG/EM/TDEM/SurveyTDEM.py | 5 +- SimPEG/EM/TDEM/TDEM.py | 64 +++++++++++++------- tests/em/tdem/test_TDEM_b_DerivAdjoint.py | 47 +++++++++++---- 4 files changed, 116 insertions(+), 72 deletions(-) diff --git a/SimPEG/EM/FDEM/FDEM.py b/SimPEG/EM/FDEM/FDEM.py index 603d322c..e8cc9ba2 100644 --- a/SimPEG/EM/FDEM/FDEM.py +++ b/SimPEG/EM/FDEM/FDEM.py @@ -18,9 +18,9 @@ class BaseFDEMProblem(BaseEMProblem): {\mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{b} - \mathbf{M_{\sigma}^e} \mathbf{e} = \mathbf{s_e}} if using the E-B formulation (:code:`Problem_e` - or :code:`Problem_b`). Note that in this case, :math:`\mathbf{s_e}` is an integrated quantity. + or :code:`Problem_b`). Note that in this case, :math:`\mathbf{s_e}` is an integrated quantity. - If we write Maxwell's equations in terms of + If we write Maxwell's equations in terms of \\\(\\\mathbf{h}\\\) and current density \\\(\\\mathbf{j}\\\) .. math :: @@ -28,7 +28,7 @@ class BaseFDEMProblem(BaseEMProblem): \mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{j} + i \omega \mathbf{M_{\mu}^e} \mathbf{h} = \mathbf{s_m} \\\\ \mathbf{C} \mathbf{h} - \mathbf{j} = \mathbf{s_e} - if using the H-J formulation (:code:`Problem_j` or :code:`Problem_h`). Note that here, :math:`\mathbf{s_m}` is an integrated quantity. + if using the H-J formulation (:code:`Problem_j` or :code:`Problem_h`). Note that here, :math:`\mathbf{s_m}` is an integrated quantity. The problem performs the elimination so that we are solving the system for \\\(\\\mathbf{e},\\\mathbf{b},\\\mathbf{j} \\\) or \\\(\\\mathbf{h}\\\) """ @@ -39,7 +39,7 @@ class BaseFDEMProblem(BaseEMProblem): def fields(self, m=None): """ Solve the forward problem for the fields. - + :param numpy.array m: inversion model (nP,) :rtype numpy.array: :return F: forward solution @@ -65,9 +65,9 @@ class BaseFDEMProblem(BaseEMProblem): :param numpy.array m: inversion model (nP,) :param numpy.array v: vector which we take sensitivity product with (nP,) - :param SimPEG.EM.FDEM.Fields u: fields object + :param SimPEG.EM.FDEM.Fields u: fields object :rtype numpy.array: - :return: Jv (ndata,) + :return: Jv (ndata,) """ if u is None: @@ -85,13 +85,13 @@ class BaseFDEMProblem(BaseEMProblem): ftype = self._fieldType + 'Solution' u_src = u[src, ftype] dA_dm_v = self.getADeriv(freq, u_src, v) - dRHS_dm_v = self.getRHSDeriv(freq, src, v) + dRHS_dm_v = self.getRHSDeriv(freq, src, v) du_dm_v = Ainv * ( - dA_dm_v + dRHS_dm_v ) - + for rx in src.rxList: df_dmFun = getattr(u, '_%sDeriv'%rx.projField, None) df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False) - df_dm_v = np.array(df_dm_v,dtype=complex) + df_dm_v = np.array(df_dm_v, dtype=complex) Jv[src, rx] = rx.evalDeriv(src, self.mesh, u, df_dm_v) Ainv.clean() return Utils.mkvc(Jv) @@ -102,9 +102,9 @@ class BaseFDEMProblem(BaseEMProblem): :param numpy.array m: inversion model (nP,) :param numpy.array v: vector which we take adjoint product with (nP,) - :param SimPEG.EM.FDEM.Fields u: fields object + :param SimPEG.EM.FDEM.Fields u: fields object :rtype numpy.array: - :return: Jv (ndata,) + :return: Jv (ndata,) """ if u is None: @@ -148,7 +148,7 @@ class BaseFDEMProblem(BaseEMProblem): Jtv += - np.array(df_dmT,dtype=complex).real else: raise Exception('Must be real or imag') - + ATinv.clean() return Utils.mkvc(Jtv) @@ -158,7 +158,7 @@ class BaseFDEMProblem(BaseEMProblem): Evaluates the sources for a given frequency and puts them in matrix form :param float freq: Frequency - :rtype: (numpy.ndarray, numpy.ndarray) + :rtype: (numpy.ndarray, numpy.ndarray) :return: S_m, S_e (nE or nF, nSrc) """ Srcs = self.survey.getSrcByFreq(freq) @@ -211,7 +211,7 @@ class Problem_e(BaseFDEMProblem): def getA(self, freq): """ System matrix - + .. math :: \mathbf{A} = \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{C} + i \omega \mathbf{M^e_{\sigma}} @@ -234,12 +234,12 @@ class Problem_e(BaseFDEMProblem): .. math :: \\frac{\mathbf{A}(\mathbf{m}) \mathbf{v}}{d \mathbf{m}} = i \omega \\frac{d \mathbf{M^e_{\sigma}}\mathbf{v} }{d\mathbf{m}} - :param float freq: frequency - :param numpy.ndarray u: solution vector (nE,) + :param float freq: frequency + :param numpy.ndarray u: solution vector (nE,) :param numpy.ndarray v: vector to take prodct with (nP,) or (nD,) for adjoint :param bool adjoint: adjoint? :rtype: numpy.ndarray - :return: derivative of the system matrix times a vector (nP,) or adjoint (nD,) + :return: derivative of the system matrix times a vector (nP,) or adjoint (nD,) """ dsig_dm = self.curModel.sigmaDeriv @@ -252,13 +252,13 @@ class Problem_e(BaseFDEMProblem): def getRHS(self, freq): """ - Right hand side for the system + Right hand side for the system .. math :: \mathbf{RHS} = \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f}\mathbf{s_m} -i\omega\mathbf{M_e}\mathbf{s_e} :param float freq: Frequency - :rtype: numpy.ndarray + :rtype: numpy.ndarray :return: RHS (nE, nSrc) """ @@ -270,7 +270,7 @@ class Problem_e(BaseFDEMProblem): def getRHSDeriv(self, freq, src, v, adjoint=False): """ - Derivative of the right hand side with respect to the model + Derivative of the right hand side with respect to the model :param float freq: frequency :param SimPEG.EM.FDEM.Src src: FDEM source @@ -350,12 +350,12 @@ class Problem_b(BaseFDEMProblem): .. math :: \\frac{\mathbf{A}(\mathbf{m}) \mathbf{v}}{d \mathbf{m}} = \mathbf{C} \\frac{\mathbf{M^e_{\sigma}} \mathbf{v}}{d\mathbf{m}} - :param float freq: frequency - :param numpy.ndarray u: solution vector (nF,) + :param float freq: frequency + :param numpy.ndarray u: solution vector (nF,) :param numpy.ndarray v: vector to take prodct with (nP,) or (nD,) for adjoint :param bool adjoint: adjoint? :rtype: numpy.ndarray - :return: derivative of the system matrix times a vector (nP,) or adjoint (nD,) + :return: derivative of the system matrix times a vector (nP,) or adjoint (nD,) """ MfMui = self.MfMui @@ -377,13 +377,13 @@ class Problem_b(BaseFDEMProblem): def getRHS(self, freq): """ - Right hand side for the system + Right hand side for the system .. math :: \mathbf{RHS} = \mathbf{s_m} + \mathbf{M^e_{\sigma}}^{-1}\mathbf{s_e} :param float freq: Frequency - :rtype: numpy.ndarray + :rtype: numpy.ndarray :return: RHS (nE, nSrc) """ @@ -501,12 +501,12 @@ class Problem_j(BaseFDEMProblem): \\frac{\mathbf{A(\sigma)} \mathbf{v}}{d \mathbf{m}} = \mathbf{C} \mathbf{M^e_{mu^{-1}}} \mathbf{C^{\\top}} \\frac{d \mathbf{M^f_{\sigma^{-1}}}\mathbf{v} }{d \mathbf{m}} - :param float freq: frequency - :param numpy.ndarray u: solution vector (nF,) + :param float freq: frequency + :param numpy.ndarray u: solution vector (nF,) :param numpy.ndarray v: vector to take prodct with (nP,) or (nD,) for adjoint :param bool adjoint: adjoint? :rtype: numpy.ndarray - :return: derivative of the system matrix times a vector (nP,) or adjoint (nD,) + :return: derivative of the system matrix times a vector (nP,) or adjoint (nD,) """ MeMuI = self.MeMuI @@ -526,7 +526,7 @@ class Problem_j(BaseFDEMProblem): def getRHS(self, freq): """ - Right hand side for the system + Right hand side for the system .. math :: @@ -550,7 +550,7 @@ class Problem_j(BaseFDEMProblem): def getRHSDeriv(self, freq, src, v, adjoint=False): """ - Derivative of the right hand side with respect to the model + Derivative of the right hand side with respect to the model :param float freq: frequency :param SimPEG.EM.FDEM.Src src: FDEM source @@ -630,12 +630,12 @@ class Problem_h(BaseFDEMProblem): .. math:: \\frac{\mathbf{A}(\mathbf{m}) \mathbf{v}}{d \mathbf{m}} = \mathbf{C}^{\\top}\\frac{d \mathbf{M^f_{\\rho}}\mathbf{v} }{d\mathbf{m}} - :param float freq: frequency - :param numpy.ndarray u: solution vector (nE,) + :param float freq: frequency + :param numpy.ndarray u: solution vector (nE,) :param numpy.ndarray v: vector to take prodct with (nP,) or (nD,) for adjoint :param bool adjoint: adjoint? :rtype: numpy.ndarray - :return: derivative of the system matrix times a vector (nP,) or adjoint (nD,) + :return: derivative of the system matrix times a vector (nP,) or adjoint (nD,) """ MeMu = self.MeMu @@ -648,14 +648,14 @@ class Problem_h(BaseFDEMProblem): def getRHS(self, freq): """ - Right hand side for the system + Right hand side for the system .. math :: \mathbf{RHS} = \mathbf{M^e} \mathbf{s_m} + \mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{s_e} :param float freq: Frequency - :rtype: numpy.ndarray + :rtype: numpy.ndarray :return: RHS (nE, nSrc) """ @@ -667,7 +667,7 @@ class Problem_h(BaseFDEMProblem): def getRHSDeriv(self, freq, src, v, adjoint=False): """ - Derivative of the right hand side with respect to the model + Derivative of the right hand side with respect to the model :param float freq: frequency :param SimPEG.EM.FDEM.Src src: FDEM source diff --git a/SimPEG/EM/TDEM/SurveyTDEM.py b/SimPEG/EM/TDEM/SurveyTDEM.py index 2e8d3ce7..eb534978 100644 --- a/SimPEG/EM/TDEM/SurveyTDEM.py +++ b/SimPEG/EM/TDEM/SurveyTDEM.py @@ -132,7 +132,6 @@ class BaseSrc(SimPEG.Survey.BaseSrc): @waveform.setter def waveform(self, val): if self.waveform is None: - print val val._assertMatchesPair(self.waveformPair) self._mapping = val else: @@ -252,8 +251,10 @@ class MagDipole(BaseSrc): C = prob.mesh.edgeCurl S_e = self.S_e(prob, prob.t0) + # S_e doesn't depend on the model + if adjoint: - raise NotImplementedError + return MeSigmaIDeriv( -S_e + C.T * ( MfMui * b ) ).T * v return MeSigmaIDeriv( -S_e + C.T * ( MfMui * b ) ) * v diff --git a/SimPEG/EM/TDEM/TDEM.py b/SimPEG/EM/TDEM/TDEM.py index bc7c6b4c..7c451dc6 100644 --- a/SimPEG/EM/TDEM/TDEM.py +++ b/SimPEG/EM/TDEM/TDEM.py @@ -85,8 +85,16 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem): # mat to store previous time-step's solution deriv times a vector for each source # size: nu x nSrc - dun_dm_v = self.getInitialFieldsDeriv(v) # can over-write this at each timestep + # this is a bit silly + + # if self._fieldType is 'b' or self._fieldType is 'j': + # ifields = np.zeros((self.mesh.nF, len(Srcs))) + # elif self._fieldType is 'e' or self._fieldType is 'h': + # ifields = np.zeros((self.mesh.nE, len(Srcs))) + + # for i, src in enumerate(self.survey.srcList): + dun_dm_v = np.hstack([Utils.mkvc(self.getInitialFieldsDeriv(src,v),2) for src in self.survey.srcList]) # can over-write this at each timestep # df_dm_v = Fields_Derivs(self.mesh, self.survey) # store the field derivs we need to project to calc full deriv @@ -117,10 +125,7 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem): dA_dm_v = self.getAdiagDeriv(tInd, un_src, v) # cell centered on time mesh dRHS_dm_v = self.getRHSDeriv(tInd+1, src, v) # on nodes of time mesh - # if tInd >= 1: dAsubdiag_dm_v = self.getAsubdiagDeriv(tInd, u[src,ftype,tInd], v) - # else: - # dAsubdiag_dm_v = Zero() JRHS = dRHS_dm_v - dAsubdiag_dm_v - dA_dm_v @@ -201,6 +206,8 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem): Adiag = self.getAdiag(tInd) AdiagTinv = self.Solver(Adiag.T, **self.solverOpts) + dAsubdiag_dm_v = Zero() + if tInd < self.nT - 1: Asubdiag = self.getAsubdiag(tInd+1) @@ -211,20 +218,35 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem): # last timestep (first to be solved) ATinv_df_duT_v[isrc,:] = AdiagTinv * df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1] elif tInd > -1: + # else: 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,:])) + 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,:]) + # - Utils.mkvc(Asubdiag.T * Utils.mkvc(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,:])) + + if tInd < self.nT - 1: + dAsubdiagT_dm_v = self.getAsubdiagDeriv(tInd+1, u[src,ftype,tInd+1], ATinv_df_duT_v[isrc,:], adjoint = True) 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' - # un_src = u[src,ftype,tInd+1] + JTv = JTv + Utils.mkvc(- dAT_dm_v - dAsubdiag_dm_v + dRHST_dm_v) + else: + # dA_dm_v = self.getInitialFieldsDeriv(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1], adjoint=True) + # print np.linalg.norm(self.getInitialFieldsDeriv(src, df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1], adjoint=True)) + # print np.linalg.norm(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1]) + # vec = - Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:]) + Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1]) + # dAsubdiagT_dm_v = self.getAsubdiagDeriv(tInd+1, u[src,ftype,tInd+1], Utils.mkvc(ATinv_df_duT_v[isrc,:]), adjoint = True) + dRHST_dm_v = Utils.mkvc(self.getInitialFieldsDeriv(src, Utils.mkvc(ATinv_df_duT_v[isrc,:]) , adjoint=True)) + + JTv = JTv + Utils.mkvc( -dAsubdiagT_dm_v + dRHST_dm_v) # + + + # # 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) @@ -286,21 +308,20 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem): return ifields - def getInitialFieldsDeriv(self, v, adjoint=False): + def getInitialFieldsDeriv(self, src, v, adjoint=False): - Srcs = self.survey.srcList + if adjoint is False: + if self._fieldType is 'b' or self._fieldType is 'j': + ifieldsDeriv = np.zeros(self.mesh.nF) + elif self._fieldType is 'e' or self._fieldType is 'h': + ifieldsDeriv = np.zeros(self.mesh.nE) - if self._fieldType is 'b' or self._fieldType is 'j': - ifieldsDeriv = np.zeros((self.mesh.nF, len(Srcs))) - elif self._fieldType is 'e' or self._fieldType is 'h': - ifieldsDeriv = np.zeros((self.mesh.nE, len(Srcs))) + elif adjoint is True: + ifieldsDeriv = np.zeros(self.mapping.nP) - for i,src in enumerate(Srcs): - ifieldsDeriv[:,i] = ifieldsDeriv[:,i] + getattr(src, '%sInitialDeriv'%self._fieldType, None)(self,v,adjoint) + ifieldsDeriv = Utils.mkvc(getattr(src, '%sInitialDeriv'%self._fieldType, None)(self,v,adjoint)) + ifieldsDeriv - if adjoint is True: - ifieldsDeriv = np.zeros_like(self.curModel) - print v.shape + # ifieldsDeriv = Utils.mkvc(getattr(src, '%sInitialDeriv'%self._fieldType, None)(self,v,adjoint)) + ifieldsDeriv # ifieldsDeriv = self.getAdiagDeriv(None, u, v, adjoint) # ifieldsDeriv = ifieldsDeriv.sum() @@ -485,6 +506,7 @@ class Problem_e(BaseTDEMProblem): def getAdiagDeriv(self, tInd, u, v, adjoint=False): + assert tInd >= 0 and tInd < self.nT dt = self.timeSteps[tInd] C = self.mesh.edgeCurl @@ -492,7 +514,7 @@ class Problem_e(BaseTDEMProblem): MeSigmaDeriv = self.MeSigmaDeriv(u) if adjoint: - raise 1./dt * MeSigmaDeriv.T * v + return 1./dt * MeSigmaDeriv.T * v return 1./dt * MeSigmaDeriv * v diff --git a/tests/em/tdem/test_TDEM_b_DerivAdjoint.py b/tests/em/tdem/test_TDEM_b_DerivAdjoint.py index a297e6cb..c3ac40ba 100644 --- a/tests/em/tdem/test_TDEM_b_DerivAdjoint.py +++ b/tests/em/tdem/test_TDEM_b_DerivAdjoint.py @@ -5,7 +5,7 @@ from SimPEG import EM plotIt = False testDeriv = True -testAdjoint = False +testAdjoint = True TOL = 1e-5 @@ -54,13 +54,16 @@ class TDEM_DerivTests(unittest.TestCase): # ====== TEST A ========== # - def test_AderivTest(self): - prb, m0, mesh = setUp() + def AderivTest(self, prbtype): + prb, m0, mesh = setUp(prbtype) tInd = 2 + if prbtype == 'b': + nu = mesh.nF + elif prbtype == 'e': + nu = mesh.nE + v = np.random.rand(nu) - v = np.random.rand(mesh.nF) - - def AderivTest(m): + def AderivFun(m): prb.curModel = m A = prb.getAdiag(tInd) Av = A*v @@ -69,26 +72,41 @@ class TDEM_DerivTests(unittest.TestCase): return Av, ADeriv_dm - print '\n Testing ADeriv' - Tests.checkDerivative(AderivTest, m0, plotIt=False, num=4, eps=1e-20) + print '\n Testing ADeriv %s'%(prbtype) + Tests.checkDerivative(AderivFun, m0, plotIt=False, num=4, eps=1e-20) - def test_A_adjointTest(self): - prb, m0, mesh = setUp() + def A_adjointTest(self,prbtype): + prb, m0, mesh = setUp(prbtype) 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) + if prbtype == 'b': + nu = prb.mesh.nF + elif prbtype == 'e': + nu = prb.mesh.nE + + v = np.random.rand(nu) + u = np.random.rand(nu) 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 + print 'AdjointTest %s'%(prbtype), V1, V2, passed self.assertTrue(passed) + def test_Aderiv_b(self): + self.AderivTest('b') + def test_Aderiv_e(self): + self.AderivTest('e') + + def test_Aadjoint_b(self): + self.A_adjointTest('b') + def test_Aadjoint_e(self): + self.A_adjointTest('e') + # ====== TEST Fields Deriv Pieces ========== # def test_eDeriv_m_adjoint(self): @@ -194,6 +212,9 @@ class TDEM_DerivTests(unittest.TestCase): def test_Jvec_adjoint_b_ey(self): self.JvecVsJtvecTest('b', 'ey') + def test_Jvec_adjoint_e_ey(self): + self.JvecVsJtvecTest('e', 'ey') + if __name__ == '__main__':