From 5cf0acd153577706e5c1a4b218fcdd91e23d0bd6 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sun, 6 Mar 2016 16:06:14 -0800 Subject: [PATCH] TDEM bderiv from b formulation working --- SimPEG/EM/TDEM/TDEM.py | 18 ++++------ tests/em/tdem/test_TDEM_b_DerivAdjoint.py | 42 +++++++++++------------ 2 files changed, 27 insertions(+), 33 deletions(-) diff --git a/SimPEG/EM/TDEM/TDEM.py b/SimPEG/EM/TDEM/TDEM.py index a3dbac07..8f0a8759 100644 --- a/SimPEG/EM/TDEM/TDEM.py +++ b/SimPEG/EM/TDEM/TDEM.py @@ -68,9 +68,12 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem): self.curModel = m Jv = self.dataPair(self.survey) + + # mat to store previous time-step's solution deriv times a vector for each source dun_dm_v = self.getInitialFieldsDeriv(v) # can over-write this at each timestep - df_dm = Fields_Derivs(self.mesh, self.survey) + df_dm_v = Fields_Derivs(self.mesh, self.survey) # store the field derivs we need to project to calc full deriv + Ainv = None for tInd, dt in enumerate(self.timeSteps): @@ -84,12 +87,12 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem): for i, src in enumerate(self.survey.srcList): # compute next du_dm_v for next timestep - u_src = u[src,ftype,tInd] + u_src = u[src,ftype,tInd+1] rhs_v = self.getJRHS(tInd, src, u_src, v, dun_dm_v[:,i]) for rx in src.rxList: df_dmFun = getattr(u, '_%sDeriv'%rx.projField, None) - df_dm[src, '%sDeriv'%rx.projField , tInd] = df_dmFun(src, dun_dm_v, v) + df_dm_v[src, '%sDeriv'%rx.projField , tInd] = df_dmFun(src, dun_dm_v[:,i], v) # over-write with this time-steps (if not on last timestep) if tInd != len(self.timeSteps): @@ -97,18 +100,11 @@ 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) + Jv[src,rx] = rx.evalDeriv(src, self.mesh, self.timeMesh, df_dm_v) Ainv.clean() return Utils.mkvc(Jv) - - - # if tInd == 0: - # dbn_dm_v = get - # else: - - # rhs_v = self.getJRHS(tInd, m, u, v, dbn_dm_v) def getJRHS(self, tInd, src, u, v, dbn_dm_v, adjoint = False): diff --git a/tests/em/tdem/test_TDEM_b_DerivAdjoint.py b/tests/em/tdem/test_TDEM_b_DerivAdjoint.py index ca311958..988cec87 100644 --- a/tests/em/tdem/test_TDEM_b_DerivAdjoint.py +++ b/tests/em/tdem/test_TDEM_b_DerivAdjoint.py @@ -11,12 +11,12 @@ class TDEM_bDerivTests(unittest.TestCase): cs = 5. ncx = 20 - ncy = 6 + ncy = 10 npad = 20 hx = [(cs,ncx), (cs,npad,1.3)] hy = [(cs,npad,-1.3), (cs,ncy), (cs,npad,1.3)] mesh = Mesh.CylMesh([hx,1,hy], '00C') - +# active = mesh.vectorCCz<0. activeMap = Maps.InjectActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz) mapping = Maps.ExpMap(mesh) * Maps.SurjectVertical1D(mesh) * activeMap @@ -193,31 +193,29 @@ class TDEM_bDerivTests(unittest.TestCase): Tests.checkDerivative(AderivTest, self.m, plotIt=False, num=4, eps=1e-20) - def test_Fields_Deriv(self): - prb = self.prb - tInd = 0 + # def test_Fields_Deriv(self): + # prb = self.prb + # tInd = 10 - v = np.random.rand(self.mesh.nF) + # 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) + # def FieldsDerivs(m): - sol = Ainv * RHS - Ainv.clean() + # sol = prb.fields(m)[:,'bSolution',tInd] - prb.curModel = self.m - f = prb.fields(self.m) - deriv = lambda dm: f._bDeriv(prb.survey.srcList[0], np.zeros_like(sol), dm) + # prb.curModel = self.m + # f = prb.fields(self.m) - return sol, deriv + # df_dm_v = EM.TDEM.FieldsTDEM.Fields_Derivs(mesh, survey) + # for i in range(tInd): + # Ainv = prb.Solver(prb.getA(tInd)) + # df_dm_v = Ainv * - Tests.checkDerivative(FieldsDerivs, self.m, plotIt=False, num=4, eps=1e-20) + # deriv = lambda dm: f._bDeriv(prb.survey.srcList[0], f[:,'bSolution',tInd], dm) + + # return sol, deriv + + # Tests.checkDerivative(FieldsDerivs, self.m, plotIt=False, num=4, eps=1e-20) def test_Deriv_J(self): @@ -233,7 +231,7 @@ class TDEM_bDerivTests(unittest.TestCase): 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, num=5, eps=1e-20) + Tests.checkDerivative(derChk, m, plotIt=False, num=2, eps=1e-20) # def test_projectAdjoint(self): # prb = self.prb