From 664adb04acd65b3557b419ac76ba863d01a38cc7 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sun, 6 Mar 2016 15:10:13 -0800 Subject: [PATCH] light notation cleanup in Jvec, testing ADeriv --> passes, Jvec is still first order --- SimPEG/EM/TDEM/TDEM.py | 26 ++++------ tests/em/tdem/test_TDEM_b_DerivAdjoint.py | 63 +++++++++++++++++++---- 2 files changed, 64 insertions(+), 25 deletions(-) diff --git a/SimPEG/EM/TDEM/TDEM.py b/SimPEG/EM/TDEM/TDEM.py index 261b877c..a3dbac07 100644 --- a/SimPEG/EM/TDEM/TDEM.py +++ b/SimPEG/EM/TDEM/TDEM.py @@ -61,8 +61,6 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem): def Jvec(self, m, v, u=None): - # raise NotImplementedError - if u is None: u = self.fields(m) @@ -70,22 +68,19 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem): self.curModel = m Jv = self.dataPair(self.survey) - dun_dm_v = self.getInitialFieldsDeriv(v) # can over-write this at each timestep df_dm = Fields_Derivs(self.mesh, self.survey) - # np.zeros((dun_dm_v.shape[0], self.nT, self.survey.nSrc)) - - Adiag, Asub = None, None + Ainv = None for tInd, dt in enumerate(self.timeSteps): - if Adiag 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() - Adiag, Asub = None, None + if Ainv 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 + Ainv.clean() + Ainv = None - if Adiag is None: - Adiag, Asub = self.getJdiags(tInd, adjoint = False) - Adiaginv = self.Solver(Adiag) + if Ainv is None: + A = self.getA(tInd) + Ainv = self.Solver(A, **self.solverOpts) for i, src in enumerate(self.survey.srcList): # compute next du_dm_v for next timestep @@ -98,13 +93,13 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem): # over-write with this time-steps (if not on last timestep) if tInd != len(self.timeSteps): - dun_dm_v[:,i] = Adiaginv * rhs_v + dun_dm_v[:,i] = Ainv * rhs_v for src in self.survey.srcList: for rx in src.rxList: Jv[src,rx] = rx.evalDeriv(src, self.mesh, self.timeMesh, df_dm) - Adiaginv.clean() + Ainv.clean() return Utils.mkvc(Jv) @@ -248,9 +243,8 @@ class Problem_b(BaseTDEMProblem): def getADeriv(self, tInd, u, v, adjoint=False): C = self.mesh.edgeCurl - MeSigmaIDeriv = lambda u: self.MeSigmaIDeriv(u) + MeSigmaIDeriv = lambda x: self.MeSigmaIDeriv(x) MfMui = self.MfMui - I = Utils.speye(self.mesh.nF) 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 4dd443a8..ca311958 100644 --- a/tests/em/tdem/test_TDEM_b_DerivAdjoint.py +++ b/tests/em/tdem/test_TDEM_b_DerivAdjoint.py @@ -21,8 +21,8 @@ class TDEM_bDerivTests(unittest.TestCase): activeMap = Maps.InjectActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz) mapping = Maps.ExpMap(mesh) * Maps.SurjectVertical1D(mesh) * activeMap - rxOffset = 40. - rx = EM.TDEM.Rx(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20), 'bz') + rxOffset = 10. + rx = EM.TDEM.Rx(np.array([[rxOffset, 0., -1e-2]]), np.logspace(-4,-3, 20), 'bz') src = EM.TDEM.SurveyTDEM.MagDipole([rx], loc=np.array([0., 0., 0.])) survey = EM.TDEM.Survey([src]) @@ -30,6 +30,7 @@ class TDEM_bDerivTests(unittest.TestCase): self.prb = EM.TDEM.Problem_b(mesh, mapping=mapping) # self.prb.timeSteps = [1e-5] self.prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)] + # self.prb.__makeASymmetric = False # self.prb.timeSteps = [(1e-05, 100)] try: @@ -38,9 +39,10 @@ class TDEM_bDerivTests(unittest.TestCase): except ImportError, e: self.prb.Solver = SolverLU - self.sigma = np.ones(mesh.nCz)*1e-8 - self.sigma[mesh.vectorCCz<0] = 1e-1 - self.m = np.log(self.sigma[active]) + # self.sigma = np.ones(mesh.nCz)*1e-8 + # self.sigma[active] = 1e-1 + # self.sigma[active] += 1e-2*np.random.rand(len(active)) + self.m = np.log(1e-1)*np.ones(self.prb.mapping.nP) + 1e-2*np.random.randn(self.prb.mapping.nP) self.prb.pair(survey) self.mesh = mesh @@ -174,21 +176,64 @@ class TDEM_bDerivTests(unittest.TestCase): # print 'test_Deriv_dUdM' # Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20) + def test_ADeriv(self): + prb = self.prb + tInd = 0 + + v = np.random.rand(self.mesh.nF) + + def AderivTest(m): + prb.curModel = m + A = prb.getA(tInd) + Av = A*v + prb.curModel = self.m + ADeriv_dm = lambda dm: prb.getADeriv(tInd, v, dm) + + return Av, ADeriv_dm + + Tests.checkDerivative(AderivTest, self.m, plotIt=False, num=4, eps=1e-20) + + def test_Fields_Deriv(self): + prb = self.prb + tInd = 0 + + v = np.random.rand(self.mesh.nF) + + def FieldsDerivs(m): + prb.curModel = m + F = prb.fieldsPair(self.mesh, self.prb.survey) + # set initial fields + F[:,prb._fieldType+'Solution',0] = prb.getInitialFields() + A = prb.getA(tInd) + Ainv = prb.Solver(A) + RHS = prb.getRHS(tInd, F) + + sol = Ainv * RHS + Ainv.clean() + + prb.curModel = self.m + f = prb.fields(self.m) + deriv = lambda dm: f._bDeriv(prb.survey.srcList[0], np.zeros_like(sol), dm) + + return sol, deriv + + Tests.checkDerivative(FieldsDerivs, self.m, plotIt=False, num=4, eps=1e-20) + def test_Deriv_J(self): prb = self.prb - prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)] + # prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)] mesh = self.mesh m = self.m # d_sig = 0.8*sigma #np.random.rand(mesh.nCz) - d_m = np.random.rand(prb.mapping.nP) + # d_m = 0.1*np.random.randn(prb.mapping.nP) derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(m, mx)] print '\n' print 'test_Deriv_J' - Tests.checkDerivative(derChk, m, plotIt=False, dx=d_m, num=4, eps=1e-20) + Tests.checkDerivative(derChk, m, plotIt=False, num=5, eps=1e-20) # def test_projectAdjoint(self): # prb = self.prb @@ -202,7 +247,7 @@ class TDEM_bDerivTests(unittest.TestCase): # f[:,'e',i] = np.random.rand(mesh.nE, 1) # d_vec = np.random.rand(survey.nD) # d = Survey.Data(survey,v=d_vec) - + # # Check that d.T*Q*f = f.T*Q.T*d # V1 = d_vec.dot(survey.evalDeriv(None, v=f).tovec()) # V2 = f.tovec().dot(survey.evalDeriv(None, v=d, adjoint=True).tovec())