From 1d2eac62a3902e9098e07eb5492228bab450c60f Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sun, 6 Mar 2016 21:42:30 -0800 Subject: [PATCH] e hooked up with Jvec --- SimPEG/EM/TDEM/FieldsTDEM.py | 40 ++-- SimPEG/EM/TDEM/TDEM.py | 2 +- tests/em/tdem/test_TDEM_b_DerivAdjoint.py | 259 +++++----------------- 3 files changed, 81 insertions(+), 220 deletions(-) diff --git a/SimPEG/EM/TDEM/FieldsTDEM.py b/SimPEG/EM/TDEM/FieldsTDEM.py index 5ba3abf3..a15d93c0 100644 --- a/SimPEG/EM/TDEM/FieldsTDEM.py +++ b/SimPEG/EM/TDEM/FieldsTDEM.py @@ -49,35 +49,49 @@ class Fields_b(Fields): } def startup(self): - self.MeSigmaI = self.survey.prob.MeSigmaI - self.edgeCurl = self.survey.prob.mesh.edgeCurl - self.MfMui = self.survey.prob.MfMui + self.MeSigmaI = self.survey.prob.MeSigmaI + self.MeSigmaIDeriv = self.survey.prob.MeSigmaIDeriv + self.edgeCurl = self.survey.prob.mesh.edgeCurl + self.MfMui = self.survey.prob.MfMui def _b(self, bSolution, srcList, tInd): return bSolution - def _bDeriv_u(self, src, dun_dm_v, adjoint = False): + def _bDeriv_u(self, tInd, src, dun_dm_v, adjoint = False): return Identity()*dun_dm_v - def _bDeriv_m(self, src, v, adjoint = False): + def _bDeriv_m(self, tInd, src, v, adjoint = False): return Zero() - def _bDeriv(self, src, dun_dm_v, v, adjoint=False): + def _bDeriv(self, tInd, src, dun_dm_v, v, adjoint=False): if adjoint is True: raise NotImplementedError - return self._bDeriv_u(src, dun_dm_v) + self._bDeriv_m(src, v) + return self._bDeriv_u(tInd, src, dun_dm_v) + self._bDeriv_m(tInd, src, v) 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.prob, tInd) - e[:,i,tInd] = e[:,i,tInd] - self.MeSigmaI * S_e + _, S_e = src.eval(self.survey.prob, self.survey.prob.times[tInd]) + e[:,i] = e[:,i] - self.MeSigmaI * S_e return e - def _eDeriv_u(self, src, dun_dm_v, adjoint = False): - raise NotImplementedError + def _eDeriv_u(self, tInd, src, dun_dm_v, adjoint = False): + if adjoint is True: + raise NotImplementedError + return self.MeSigmaI * ( self.edgeCurl.T * ( self.MfMui * dun_dm_v ) ) - def _eDeriv_m(self, src, v, adjoint = False): - raise NotImplementedError + def _eDeriv_m(self, tInd, src, v, adjoint = False): + if adjoint is True: + raise NotImplementedError + 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) + + return self.MeSigmaIDeriv(self.edgeCurl.T * ( self.MfMui * bSolution) ) * v - self.MeSigmaIDeriv(S_e) * v - self.MeSigmaI * S_eDeriv + + def _eDeriv(self, tInd, src, dun_dm_v, v, adjoint=False): + if adjoint is True: + raise NotImplementedError + 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 8f0a8759..ff4f8cb8 100644 --- a/SimPEG/EM/TDEM/TDEM.py +++ b/SimPEG/EM/TDEM/TDEM.py @@ -92,7 +92,7 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem): for rx in src.rxList: df_dmFun = getattr(u, '_%sDeriv'%rx.projField, None) - df_dm_v[src, '%sDeriv'%rx.projField , tInd] = df_dmFun(src, dun_dm_v[:,i], v) + df_dm_v[src, '%sDeriv'%rx.projField , tInd] = df_dmFun(tInd, src, dun_dm_v[:,i], v) # over-write with this time-steps (if not on last timestep) if tInd != len(self.timeSteps): diff --git a/tests/em/tdem/test_TDEM_b_DerivAdjoint.py b/tests/em/tdem/test_TDEM_b_DerivAdjoint.py index 988cec87..cdf10e05 100644 --- a/tests/em/tdem/test_TDEM_b_DerivAdjoint.py +++ b/tests/em/tdem/test_TDEM_b_DerivAdjoint.py @@ -5,234 +5,81 @@ from SimPEG import EM plotIt = False tol = 1e-6 +def setUp(rxcomp='bz'): + cs = 5. + ncx = 20 + 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 + + rxOffset = 10. + rx = EM.TDEM.Rx(np.array([[rxOffset, 0., -1e-2]]), np.logspace(-4,-3, 20), rxcomp) + src = EM.TDEM.SurveyTDEM.MagDipole([rx], loc=np.array([0., 0., 0.])) + + survey = EM.TDEM.Survey([src]) + + prb = EM.TDEM.Problem_b(mesh, mapping=mapping) + prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)] + + try: + from pymatsolver import MumpsSolver + prb.Solver = MumpsSolver + except ImportError, e: + prb.Solver = SolverLU + + m = np.log(1e-1)*np.ones(prb.mapping.nP) + 1e-2*np.random.randn(prb.mapping.nP) + + prb.pair(survey) + mesh = mesh + + return prb, m, mesh + class TDEM_bDerivTests(unittest.TestCase): - def setUp(self): - - cs = 5. - ncx = 20 - 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 - - 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]) - - 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: - from pymatsolver import MumpsSolver - self.prb.Solver = MumpsSolver - except ImportError, e: - self.prb.Solver = SolverLU - - # 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 - - # def test_AhVec(self): - # """ - # Test that fields and AhVec produce consistent results - # """ - - # prb = self.prb - # sigma = self.sigma - - # u = prb.fields(sigma) - # Ahu = prb._AhVec(sigma, u) - - # V1 = Ahu[:,'b',1] - # V2 = 1./prb.timeSteps[0]*prb.MfMui*u[:,'b',0] - # self.assertLess(np.linalg.norm(V1-V2)/np.linalg.norm(V2), 1.e-6) - - # V1 = Ahu[:,'e',1] - # return np.linalg.norm(V1) < 1.e-6 - - # for i in range(2,prb.nT): - - # dt = prb.timeSteps[i] - - # V1 = Ahu[:,'b',i] - # V2 = 1.0/dt*prb.MfMui*u[:,'b', i-1] - # # print np.linalg.norm(V1), np.linalg.norm(V2) - # self.assertLess(np.linalg.norm(V1)/np.linalg.norm(V2), 1.e-6) - - # V1 = Ahu[:,'e',i] - # V2 = prb.MeSigma*u[:,'e',i] - # # print np.linalg.norm(V1), np.linalg.norm(V2) - # return np.linalg.norm(V1)/np.linalg.norm(V2), 1.e-6 - - # def test_AhVecVSMat_OneTS(self): - - # prb = self.prb - # prb.timeSteps = [1e-05] - # sigma = self.sigma - # prb.curModel = sigma - - # dt = prb.timeSteps[0] - # a11 = 1/dt*prb.MfMui*sp.identity(prb.mesh.nF) - # a12 = prb.MfMui*prb.mesh.edgeCurl - # a21 = prb.mesh.edgeCurl.T*prb.MfMui - # a22 = -prb.MeSigma - # A = sp.bmat([[a11,a12],[a21,a22]]) - - # f = prb.fields(sigma) - # u1 = A*f.tovec() - # u2 = prb._AhVec(sigma,f).tovec() - - # self.assertTrue(np.linalg.norm(u1-u2)/np.linalg.norm(u1)<1e-12) - - # def test_solveAhVSMat_OneTS(self): - # prb = self.prb - - # prb.timeSteps = [1e-05] - - # sigma = self.sigma - # prb.curModel = sigma - - # dt = prb.timeSteps[0] - # a11 = 1.0/dt*prb.MfMui*sp.identity(prb.mesh.nF) - # a12 = prb.MfMui*prb.mesh.edgeCurl - # a21 = prb.mesh.edgeCurl.T*prb.MfMui - # a22 = -prb.MeSigma - # A = sp.bmat([[a11,a12],[a21,a22]]) - - # f = prb.fields(sigma) - # f[:,:,0] = {'b':0} - # f[:,'b',1] = 0 - - # self.assertTrue(np.all(np.r_[f[:,'b',1],f[:,'e',1]] == f.tovec())) - - # u1 = prb.solveAh(sigma,f).tovec().flatten() - # u2 = sp.linalg.spsolve(A.tocsr(),f.tovec()) - - # self.assertTrue(np.linalg.norm(u1-u2)<1e-8) - - # def test_solveAhVsAhVec(self): - - # prb = self.prb - # mesh = self.prb.mesh - # sigma = self.sigma - # self.prb.curModel = sigma - - # f = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey) - # f[:,'b',:] = 0.0 - # for i in range(prb.nT): - # f[:,'e', i] = np.random.rand(mesh.nE, 1) - - # Ahf = prb._AhVec(sigma, f) - # f_test = prb.solveAh(sigma, Ahf) - - # u1 = f.tovec() - # u2 = f_test.tovec() - # self.assertTrue(np.linalg.norm(u1-u2)<1e-8) - - # def test_DerivG(self): - # """ - # Test the derivative of c with respect to sigma - # """ - - # # Random model and perturbation - # sigma = np.random.rand(self.prb.mapping.nP) - - # f = self.prb.fields(sigma) - # dm = 1000*np.random.rand(self.prb.mapping.nP) - # h = 0.01 - - # derChk = lambda m: [self.prb._AhVec(m, f).tovec(), lambda mx: self.prb.Gvec(sigma, mx, u=f).tovec()] - # print '\ntest_DerivG' - # passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20) - # return passed - - # def test_Deriv_dUdM(self): - - # prb = self.prb - # prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)] - # mesh = self.mesh - # sigma = self.sigma - - # dm = 10*np.random.rand(prb.mapping.nP) - # f = prb.fields(sigma) - - # derChk = lambda m: [self.prb.fields(m).tovec(), lambda mx: -prb.solveAh(sigma, prb.Gvec(sigma, mx, u=f)).tovec()] - # print '\n' - # print 'test_Deriv_dUdM' - # Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20) + def test_ADeriv(self): - prb = self.prb + prb, m0, mesh = setUp() tInd = 0 - v = np.random.rand(self.mesh.nF) + v = np.random.rand(mesh.nF) def AderivTest(m): prb.curModel = m A = prb.getA(tInd) Av = A*v - prb.curModel = self.m + prb.curModel = m0 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) + Tests.checkDerivative(AderivTest, m0, plotIt=False, num=4, eps=1e-20) - # def test_Fields_Deriv(self): - # prb = self.prb - # tInd = 10 - - # v = np.random.rand(self.mesh.nF) - - # def FieldsDerivs(m): - - # sol = prb.fields(m)[:,'bSolution',tInd] - - # prb.curModel = self.m - # f = prb.fields(self.m) - - # 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 * - - # 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): - - prb = self.prb - # 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 = 0.1*np.random.randn(prb.mapping.nP) + def JvecTest(self, rxcomp): + prb, m, mesh = setUp(rxcomp) derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(m, mx)] print '\n' - print 'test_Deriv_J' + print 'test_Jvec_%s' %(rxcomp) Tests.checkDerivative(derChk, m, plotIt=False, num=2, eps=1e-20) + def test_Jvec_b_bx(self): + self.JvecTest('bx') + + def test_Jvec_b_bz(self): + self.JvecTest('bz') + + def test_Jvec_b_ey(self): + self.JvecTest('ey') + + # def test_projectAdjoint(self): # prb = self.prb # survey = prb.survey