From fe91312917cfde73169950aacd5f119a4ac8fd0c Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sun, 13 Mar 2016 14:02:32 -0700 Subject: [PATCH] bx, bz running and passing, ey failing adjoint --> I think the initial fields are not being taken care of correctly in the deriv --- SimPEG/EM/TDEM/FieldsTDEM.py | 33 ++++++++++--------- SimPEG/EM/TDEM/TDEM.py | 24 +++++++++----- tests/em/tdem/test_TDEM_b_DerivAdjoint.py | 20 ++++++----- .../tdem/test_TDEM_b_MultiSrc_DerivAdjoint.py | 8 +++-- 4 files changed, 51 insertions(+), 34 deletions(-) diff --git a/SimPEG/EM/TDEM/FieldsTDEM.py b/SimPEG/EM/TDEM/FieldsTDEM.py index 0f4e78a5..6129cf01 100644 --- a/SimPEG/EM/TDEM/FieldsTDEM.py +++ b/SimPEG/EM/TDEM/FieldsTDEM.py @@ -7,7 +7,7 @@ from SimPEG.Utils import Zero, Identity class Fields(SimPEG.Problem.TimeFields): """ - + Fancy Field Storage for a TDEM survey. Only one field type is stored for each problem, the rest are computed. The fields obejct acts like an array and is indexed by @@ -29,7 +29,7 @@ class Fields(SimPEG.Problem.TimeFields): """ knownFields = {} - dtype = float + dtype = float @@ -40,7 +40,7 @@ class Fields_Derivs(Fields): 'hDeriv': 'E', 'jDeriv': 'F' } - + class Fields_b(Fields): """Fancy Field Storage for a TDEM survey.""" @@ -65,7 +65,7 @@ class Fields_b(Fields): def _bDeriv_m(self, tInd, src, v, adjoint=False): return Zero() - def _bDeriv(self, tInd, src, dun_dm_v, v, adjoint=False): + def _bDeriv(self, tInd, src, dun_dm_v, v, adjoint=False): if adjoint is True: return self._bDeriv_u(tInd, src, v, adjoint), self._bDeriv_m(tInd, src, v, adjoint) return self._bDeriv_u(tInd, src, dun_dm_v) + self._bDeriv_m(tInd, src, v) @@ -73,27 +73,30 @@ class Fields_b(Fields): def _e(self, bSolution, srcList, tInd): e = self.MeSigmaI * ( self.edgeCurl.T * ( self.MfMui * bSolution ) ) for i, src in enumerate(srcList): - _, S_e = src.eval(self.survey.prob, self.survey.prob.times[tInd]) + _, S_e = src.eval(self.survey.prob, self.survey.prob.times[tInd]) e[:,i] = e[:,i] - self.MeSigmaI * S_e - return e + return e def _eDeriv_u(self, tInd, src, dun_dm_v, adjoint = False): if adjoint is True: - raise NotImplementedError + return self.MfMui.T * ( self.edgeCurl * ( self.MeSigmaI.T * dun_dm_v ) ) return self.MeSigmaI * ( self.edgeCurl.T * ( self.MfMui * dun_dm_v ) ) def _eDeriv_m(self, tInd, src, v, adjoint = False): - if adjoint is True: - raise NotImplementedError - + _, S_e = src.eval(self.survey.prob, self.survey.prob.times[tInd]) bSolution = self[[src],'bSolution',tInd] - _, S_e = src.eval(self.survey.prob, self.survey.prob.times[tInd]) - _, S_eDeriv = src.evalDeriv(self.survey.prob.times[tInd], self, v=v) + if adjoint: + v = self.MeSigmaI.T * v - return self.MeSigmaIDeriv(self.edgeCurl.T * ( self.MfMui * bSolution) ) * v - self.MeSigmaIDeriv(S_e) * v - self.MeSigmaI * S_eDeriv + _, S_eDeriv_v = src.evalDeriv(self.survey.prob.times[tInd], self, v=v, adjoint=adjoint) - def _eDeriv(self, tInd, src, dun_dm_v, v, adjoint=False): if adjoint is True: - raise NotImplementedError + 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)) * v - self.MeSigmaI * S_eDeriv_v + + def _eDeriv(self, tInd, src, dun_dm_v, v, adjoint=False): + if adjoint is True: + return self._eDeriv_u(tInd, src, v, adjoint), self._eDeriv_m(tInd, src, v, adjoint) return self._eDeriv_u(tInd, src, dun_dm_v) + self._eDeriv_m(tInd, src, v) diff --git a/SimPEG/EM/TDEM/TDEM.py b/SimPEG/EM/TDEM/TDEM.py index bccc9ed3..78fd9287 100644 --- a/SimPEG/EM/TDEM/TDEM.py +++ b/SimPEG/EM/TDEM/TDEM.py @@ -163,20 +163,25 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem): # Loop over sources and receivers to create a fields object: PT_v, df_duT_v, df_dmT_v for src in self.survey.srcList: - + PT_v = Fields_Derivs(self.mesh, self.survey) # initialize storage for PT_v (don't need to preserve over sources) # initialize size - df_duT_v[src, '%sDeriv'%self._fieldType, :]= np.zeros_like(u[src, self._fieldType, :]) + df_duT_v[src, '%sDeriv'%self._fieldType, :] = np.zeros_like(u[src, self._fieldType, :]) for rx in src.rxList: - curPT_v = rx.evalDeriv(src, self.mesh, self.timeMesh, Utils.mkvc(v[src,rx]), adjoint=True) - PT_v = np.reshape(curPT_v,(len(curPT_v)/self.timeMesh.nN, self.timeMesh.nN), order='F') + PT_v[src,'%sDeriv'%rx.projField,:] = rx.evalDeriv(src, self.mesh, self.timeMesh, Utils.mkvc(v[src,rx]), adjoint=True) + # 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) - df_duT_v_cur, df_dmT_v = df_duTFun(None, src, None, PT_v, adjoint=True) - JTv = JTv + df_dmT_v + for tInd in range(self.nT+1): + cur = df_duTFun(tInd, src, None, Utils.mkvc(PT_v[src,'%sDeriv'%rx.projField,tInd]), adjoint=True) - df_duT_v[src, '%sDeriv'%self._fieldType, :] += df_duT_v_cur + df_duT_v[src, '%sDeriv'%self._fieldType, tInd] = df_duT_v[src, '%sDeriv'%self._fieldType, tInd] + Utils.mkvc(cur[0],2) + JTv = JTv + cur[1] + + + + del PT_v # no longer need this AdiagTinv = None @@ -212,9 +217,12 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem): 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: - JTv = JTv + self.getInitialFieldsDeriv(Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,0]), adjoint=True) + for projField in set(rx.projField): + JTv = JTv + self.getInitialFieldsDeriv(Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,0]), adjoint=True) return Utils.mkvc(JTv) diff --git a/tests/em/tdem/test_TDEM_b_DerivAdjoint.py b/tests/em/tdem/test_TDEM_b_DerivAdjoint.py index c5a431df..a4c58aed 100644 --- a/tests/em/tdem/test_TDEM_b_DerivAdjoint.py +++ b/tests/em/tdem/test_TDEM_b_DerivAdjoint.py @@ -4,15 +4,15 @@ from SimPEG import EM plotIt = False -testDeriv = False +testDeriv = True testAdjoint = True -tol = 1e-6 +TOL = 1e-5 def setUp(rxcomp='bz'): cs = 5. ncx = 20 - ncy = 10 + ncy = 15 npad = 20 hx = [(cs,ncx), (cs,npad,1.3)] hy = [(cs,npad,-1.3), (cs,ncy), (cs,npad,1.3)] @@ -74,7 +74,7 @@ class TDEM_DerivTests(unittest.TestCase): 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 + passed = np.abs(V1-V2) < TOL * (np.abs(V1) + np.abs(V2))/2. print 'AdjointTest', V1, V2, passed self.assertTrue(passed) @@ -107,16 +107,18 @@ class TDEM_DerivTests(unittest.TestCase): # ====== TEST Jtvec ========== # def adjointJvecVsJtvecTest(self, rxcomp='bz'): - print '\n Adjoint Testing Jvec, Jtvec %s' %(rxcomp) - prb, m0, mesh = setUp(rxcomp) + 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)) - passed = np.abs(V1-V2)/np.abs(V1) < tol - print 'AdjointTest', V1, V2, passed + 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: diff --git a/tests/em/tdem/test_TDEM_b_MultiSrc_DerivAdjoint.py b/tests/em/tdem/test_TDEM_b_MultiSrc_DerivAdjoint.py index e5886baa..7d0fc724 100644 --- a/tests/em/tdem/test_TDEM_b_MultiSrc_DerivAdjoint.py +++ b/tests/em/tdem/test_TDEM_b_MultiSrc_DerivAdjoint.py @@ -6,6 +6,7 @@ plotIt = False testDeriv = True testAdjoint = True +TOL = 1e-6 class TDEM_bDerivTests(unittest.TestCase): @@ -70,8 +71,11 @@ class TDEM_bDerivTests(unittest.TestCase): V1 = d.dot(prb.Jvec(m0, m)) V2 = m.dot(prb.Jtvec(m0, d)) - print 'AdjointTest', V1, V2 - self.assertTrue(np.abs(V1-V2)/np.abs(V1) < 1e-6) + + 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)