From a1ecef07096e0a13eaf3281ded961487415ff915 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sat, 12 Mar 2016 15:13:17 -0800 Subject: [PATCH] first shot through of passing Jtvec for TDEM problem (code will need to be cleaned up, but it passes!) --- SimPEG/EM/TDEM/SurveyTDEM.py | 24 +++--- SimPEG/EM/TDEM/TDEM.py | 94 +++++++++++++++++------ tests/em/tdem/test_TDEM_b_DerivAdjoint.py | 39 +++++++--- 3 files changed, 110 insertions(+), 47 deletions(-) diff --git a/SimPEG/EM/TDEM/SurveyTDEM.py b/SimPEG/EM/TDEM/SurveyTDEM.py index 90f26d55..07783f29 100644 --- a/SimPEG/EM/TDEM/SurveyTDEM.py +++ b/SimPEG/EM/TDEM/SurveyTDEM.py @@ -1,8 +1,8 @@ import SimPEG from SimPEG import np, Utils -from SimPEG.Utils import Zero, Identity -from scipy.constants import mu_0 -from SimPEG.EM.Utils import * +from SimPEG.Utils import Zero, Identity +from scipy.constants import mu_0 +from SimPEG.EM.Utils import * #################################################### # Receivers @@ -60,13 +60,15 @@ class Rx(SimPEG.Survey.BaseTimeRx): u_part = Utils.mkvc(u[src, self.projField, :]) return P*u_part - def evalDeriv(self, src, mesh, timeMesh, df_dm, adjoint=False): + def evalDeriv(self, src, mesh, timeMesh, v, adjoint=False): P = self.getP(mesh, timeMesh) if not adjoint: - return P * Utils.mkvc(df_dm[src, self.projField+'Deriv', :]) + return P * v #Utils.mkvc(v[src, self.projField+'Deriv', :]) elif adjoint: - return P.T * df_dm[src, self] + # dP_dF_T = P.T * v #[src, self] + # newshape = (len(dP_dF_T)/timeMesh.nN, timeMesh.nN ) + return P.T * v #np.reshape(dP_dF_T, newshape, order='F') #################################################### # Sources @@ -83,10 +85,10 @@ class BaseWaveform(object): ), "Waveform object must be an instance of a %s BaseWaveform class."%(pair.__name__) def eval(self, time): - raise NotImplementedError + raise NotImplementedError def evalDeriv(self, time): - raise NotImplementedError # needed for E-formulation + raise NotImplementedError # needed for E-formulation class StepOffWaveform(BaseWaveform): @@ -136,7 +138,7 @@ class BaseSrc(SimPEG.Survey.BaseSrc): def __init__(self, rxList, waveform = None ): self.waveform = waveform or StepOffWaveform() - SimPEG.Survey.BaseSrc.__init__(self, rxList) + SimPEG.Survey.BaseSrc.__init__(self, rxList) def bInitial(self, prob): @@ -148,7 +150,7 @@ class BaseSrc(SimPEG.Survey.BaseSrc): def eval(self, prob, time): S_m = self.S_m(prob, time) S_e = self.S_e(prob, time) - return S_m, S_e + return S_m, S_e def evalDeriv(self, prob, time, v=None, adjoint=False): if v is not None: @@ -170,7 +172,7 @@ class BaseSrc(SimPEG.Survey.BaseSrc): class MagDipole(BaseSrc): - def __init__(self, rxList, waveform=None, loc=None, orientation='Z', moment=1., mu=mu_0): + def __init__(self, rxList, waveform=None, loc=None, orientation='Z', moment=1., mu=mu_0): self.loc = loc self.orientation = orientation diff --git a/SimPEG/EM/TDEM/TDEM.py b/SimPEG/EM/TDEM/TDEM.py index 457494be..bc4a0c5e 100644 --- a/SimPEG/EM/TDEM/TDEM.py +++ b/SimPEG/EM/TDEM/TDEM.py @@ -63,6 +63,17 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem): def Jvec(self, m, v, u=None): + """ + Jvec computes the sensitivity times a vector + + .. math:: + \mathbf{J} \mathbf{v} = \\frac{d\mathbf{P}}{d\mathbf{F}} \left( \\frac{d\mathbf{F}}{d\mathbf{u}} \\frac{d\mathbf{u}}{d\mathbf{m}} + \\frac{\partial\mathbf{F}}{\partial\mathbf{m}} \\right) \mathbf{v} + + where + + .. math:: + \mathbf{A} \\frac{d\mathbf{u}}{d\mathbf{m}} + \\frac{d\mathbf{A}(\mathbf{u})}{d\mathbf{m}} = \\frac{d \mathbf{RHS}}{d \mathbf{m}} + """ if u is None: u = self.fields(m) @@ -81,7 +92,7 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem): Adiaginv = None - for tInd, dt in zip(range(self.nT+1), self.timeSteps): + for tInd, dt in zip(range(self.nT), self.timeSteps): if Adiaginv is not None and (tInd > 0 and dt != self.timeSteps[tInd - 1]):# keep factors if dt is the same as previous step b/c A will be the same Adiaginv.clean() Adiaginv = None @@ -115,7 +126,7 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem): for src in self.survey.srcList: for rx in src.rxList: - Jv[src,rx] = rx.evalDeriv(src, self.mesh, self.timeMesh, df_dm_v) + Jv[src,rx] = rx.evalDeriv(src, self.mesh, self.timeMesh, Utils.mkvc(df_dm_v[src,'%sDeriv'%rx.projField,:])) Adiaginv.clean() return Utils.mkvc(Jv) @@ -140,61 +151,94 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem): PT_v = Fields_Derivs(self.mesh, self.survey) #PT_v is a fields object + # TODO: This will only work for b formulation right now b/c of the mesh.nF df_duT_v = np.zeros((self.mesh.nF,self.nT+1)) - ATinv_df_duT_v = np.zeros((self.mesh.nF,self.nT)) + ATinv_df_duT_v = np.zeros((self.mesh.nF,self.nT+1)) + JTv = np.zeros(m.shape) + + + # TODO : this is pretty ugly + + # Loop over sources and receivers to create a fields object: PT_v for src in self.survey.srcList: - # for rx in src.rxList: + # initialize empty fields derivs for projField in set([rx.projField for rx in src.rxList]): - PT_v[src,'%sDeriv'%projField, :] = rx.evalDeriv(src, self.mesh, self.timeMesh, v, adjoint = True) # All the fields for a given src, reciever. + PT_v[src,'%sDeriv'%projField, :] = np.zeros_like(u[src, '%s'%projField, : ]) + # loop over recievers and sum contributions to fields object + for rx in src.rxList: + # for projField in set([rx.projField for rx in src.rxList]): + curPT_v = rx.evalDeriv(src, self.mesh, self.timeMesh, Utils.mkvc(v[src,rx]), adjoint=True) + PT_v[src,'%sDeriv'%rx.projField, :] += np.reshape(curPT_v,(len(curPT_v)/self.timeMesh.nN, self.timeMesh.nN), order='F') # All the fields for a given src, reciever. + + # print np.linalg.norm(PT_v[src,'bDeriv',:]) + # for src in self.survey.srcList: + # initialize empty fields derivs + for projField in set([rx.projField for rx in src.rxList]): df_duTFun = getattr(u, '_%sDeriv'%projField, None) - # TODO: don't need to recompute df_dmT_v every time... only need it once - df_duT_v_cur, df_dmT_v = df_duTFun(None, src, None, PT_v[src,'%sDeriv'%projField,:], adjoint=True) # this seems odd + df_duT_v_cur, df_dmT_v = df_duTFun(None, src, None, PT_v[src,'%sDeriv'%projField,:], adjoint=True) + JTv = JTv + df_dmT_v df_duT_v = df_duT_v + df_duT_v_cur - AdiagTinv = None - JTv = df_dmT_v + # print np.linalg.norm(df_duT_v) - for tInd in reversed(range(self.nT)): #enumerate(reversed(list(self.timeSteps))): - if AdiagTinv is not None and (tInd < self.nT and self.timeSteps[tInd] != self.timeSteps[tInd + 1]):# keep factors if dt is the same as previous step b/c A will be the same + AdiagTinv = None + + + # for tInd in reversed(range(self.nT)): #enumerate(reversed(list(self.timeSteps))): + for tInd in reversed(range(self.nT)) : # reversed(self.timeSteps)): + if AdiagTinv is not None: # and (tInd <= self.nT and dt != self.timeSteps[tInd]): + # (tInd < self.nT and self.timeSteps[tInd] != self.timeSteps[tInd + 1]):# keep factors if dt is the same as previous step b/c A will be the same AdiagTinv.clean() AdiagTinv = None + if AdiagTinv is None: Adiag = self.getAdiag(tInd) AdiagTinv = self.Solver(Adiag.T, **self.solverOpts) - Asubdiag = self.getAsubdiag(tInd) - # solve against df_duT_v - - if tInd < self.nT: - # print Utils.mkvc(AdiagTinv * df_duT_v[:,tInd],2).shape, ATinv_df_duT_v[:,tInd].shape - ATinv_df_duT_v[:,tInd] = AdiagTinv * df_duT_v[:,tInd] + if tInd >= self.nT-1: + ATinv_df_duT_v[:,tInd+1] = AdiagTinv * df_duT_v[:,tInd+1] else: - ATinv_df_duT_v[:,tInd] = AdiagTinv * (df_duT_v[:,tInd+1] - Asubdiag.T * df_duT_v[:,tInd]) + Asubdiag = self.getAsubdiag(tInd+1) + ATinv_df_duT_v[:,tInd+1] = AdiagTinv * (df_duT_v[:,tInd+1] - Asubdiag.T * ATinv_df_duT_v[:,tInd+2]) - un_src = u[src,ftype,tInd] - dAT_dm_v = self.getAdiagDeriv(tInd, un_src, ATinv_df_duT_v[:,tInd], adjoint=True) # cell centered on time mesh + # for src in self.survey.srcList: + un_src = u[src,ftype,tInd+1] + dAT_dm_v = self.getAdiagDeriv(None, un_src, ATinv_df_duT_v[:,tInd+1], adjoint=True) # cell centered on time mesh - dRHST_dm_v = self.getRHSDeriv(tInd, src, ATinv_df_duT_v[:,tInd], adjoint=True) # on nodes of time mesh + dRHST_dm_v = self.getRHSDeriv(tInd+1, src, ATinv_df_duT_v[:,tInd+1], adjoint=True) # on nodes of time mesh # dAsubdiag_dm_v = 0 JTv = JTv + (-dAT_dm_v + dRHST_dm_v) + # JTv = JTv + + + # tInd = 0 + # un_src = u[src,ftype,tInd] + # # dAT_dm_v = self.getAdiagDeriv(None, un_src, self.getInitialFieldsDeriv(), adjoint=True) + # Asubdiag = self.getAsubdiag(tInd) + # ATinv_df_duT_v[:,tInd] = AdiagTinv * (- Asubdiag.T * df_duT_v[:,tInd]) + # # - self.getAsubdiag(tInd).T * df_duT_v[:,tInd] + # dAT_dm_v = self.getAdiagDeriv(None, un_src, ATinv_df_duT_v[:,tInd], adjoint=True) # cell centered on time mesh + + # dRHST_dm_v = self.getRHSDeriv(tInd, src, ATinv_df_duT_v[:,tInd], adjoint=True) + + # JTv = JTv + (- dAT_dm_v + dRHST_dm_v) + + + return JTv - - - # for i, src in enumerate(self.survey.srcList): # un_src = u[src,ftype,tInd+1] # fields for this source at tInd @@ -402,7 +446,7 @@ class Problem_b(BaseTDEMProblem): MfMui = self.MfMui _, S_e = src.eval(tInd, self) - S_mDeriv, S_eDeriv = src.evalDeriv(self.times[tInd], self, adjoint=adjoint) # I think this is tInd+1 ? + S_mDeriv, S_eDeriv = src.evalDeriv(self.times[tInd], self, adjoint=adjoint) if adjoint: if self._makeASymmetric is True: diff --git a/tests/em/tdem/test_TDEM_b_DerivAdjoint.py b/tests/em/tdem/test_TDEM_b_DerivAdjoint.py index 0f358499..76bb9008 100644 --- a/tests/em/tdem/test_TDEM_b_DerivAdjoint.py +++ b/tests/em/tdem/test_TDEM_b_DerivAdjoint.py @@ -4,7 +4,7 @@ from SimPEG import EM plotIt = False -testDeriv = False +testDeriv = True testAdjoint = True tol = 1e-6 @@ -30,6 +30,7 @@ def setUp(rxcomp='bz'): prb = EM.TDEM.Problem_b(mesh, mapping=mapping) prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)] + # prb.timeSteps = [(1e-05, 10), (1e-05, 50), (1e-05, 50) ] #, (2.5e-4, 10)] try: from pymatsolver import MumpsSolver @@ -48,7 +49,7 @@ class TDEM_bDerivTests(unittest.TestCase): - def test_ADeriv(self): + def test_Deriv_Pieces(self): prb, m0, mesh = setUp() tInd = 0 @@ -63,19 +64,34 @@ class TDEM_bDerivTests(unittest.TestCase): return Av, ADeriv_dm - # def A_adjointTest(): + 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)/np.abs(V1) < tol + print 'AdjointTest', V1, V2, passed + self.assertTrue(passed) + + # def P_adjointTest(): + # print '\n Testing P_adjoint' # m = np.random.rand(prb.mapping.nP) - # d = np.random.rand(prb.survey.nD) # v = np.random.rand(prb.mesh.nF) + # d = np.random.rand(prb.survey.nD) + # prb.curModel = m0 - # 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) + # for src in prb.survey.srcList: + # for rx in src.rxList: - # Tests.checkDerivative(AderivTest, m0, plotIt=False, num=4, eps=1e-20) + print '\n Testing ADeriv' + Tests.checkDerivative(AderivTest, m0, plotIt=False, num=4, eps=1e-20) + A_adjointTest() @@ -100,10 +116,11 @@ class TDEM_bDerivTests(unittest.TestCase): if testAdjoint: def test_adjointJvecVsJtvec(self): + print '\n Adjoint Testing Jvec, Jtvec' prb, m0, mesh = setUp() m = np.random.rand(prb.mapping.nP) - d = np.random.rand(prb.survey.nD) + d = np.random.randn(prb.survey.nD) V1 = d.dot(prb.Jvec(m0, m)) V2 = m.dot(prb.Jtvec(m0, d))