From ef602eaab1f03dfbfb5bba8b55170f4968f4dca4 Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Thu, 28 Apr 2016 18:13:18 -0700 Subject: [PATCH] working Jtvec --- SimPEG/EM/Static/DC/ProblemDC_2D.py | 64 ++++++++++++++-------- tests/em/static/test_DC_2D_jvecjtvecadj.py | 32 +++++------ 2 files changed, 57 insertions(+), 39 deletions(-) diff --git a/SimPEG/EM/Static/DC/ProblemDC_2D.py b/SimPEG/EM/Static/DC/ProblemDC_2D.py index 8ebb9c67..30c10d97 100644 --- a/SimPEG/EM/Static/DC/ProblemDC_2D.py +++ b/SimPEG/EM/Static/DC/ProblemDC_2D.py @@ -70,38 +70,57 @@ class BaseDCProblem_2D(BaseEMProblem): Jv[src, rx] += Jv1_temp*dky[iky] /2.*np.cos(ky*y) Jv[src, rx] += Jv0[src, rx]*dky[iky]/2.*np.cos(ky*y) Jv0[src, rx] = Jv1_temp.copy() - JV[iky,isrc,:] = Jv1_temp.copy() return Utils.mkvc(Jv) - # def Jtvec(self, m, v, f=None): - # if f is None: - # f = self.fields(m) + def Jtvec(self, m, v, f=None): + if f is None: + f = self.fields(m) - # self.curModel = m + self.curModel = m - # # Ensure v is a data object. - # if not isinstance(v, self.dataPair): - # v = self.dataPair(self.survey, v) + # Ensure v is a data object. + if not isinstance(v, self.dataPair): + v = self.dataPair(self.survey, v) - # Jtv = np.zeros(m.size) - # AT = self.getA() + Jtv = np.zeros(m.size) + + # Assume y=0. + # This needs some thoughts to implement in general when src is dipole + dky = np.diff(self.kys) + dky = np.r_[dky[0], dky] + y = 0. - # for src in self.survey.srcList: - # u_src = f[src, self._solutionType] - # for rx in src.rxList: - # PTv = rx.evalDeriv(src, self.mesh, f, v[src, rx], adjoint=True) # wrt f, need possibility wrt m - # df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None) - # df_duT, df_dmT = df_duTFun(src, None, PTv, adjoint=True) + for src in self.survey.srcList: - # ATinvdf_duT = self.Ainv * df_duT + for rx in src.rxList: - # dA_dmT = self.getADeriv(u_src, ATinvdf_duT, adjoint=True) - # dRHS_dmT = self.getRHSDeriv(src, ATinvdf_duT, adjoint=True) - # du_dmT = -dA_dmT + dRHS_dmT - # Jtv += df_dmT + du_dmT + Jtv_temp1 = np.zeros(m.size) + Jtv_temp0 = np.zeros(m.size) - # return Utils.mkvc(Jtv) + for iky in range(self.nky): + u_src = f[src, self._solutionType, iky] + ky = self.kys[iky] + AT = self.getA(ky) + PTv = rx.evalDeriv(ky, src, self.mesh, f, v[src, rx], adjoint=True) # wrt f, need possibility wrt m + df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None) + df_duT, df_dmT = df_duTFun(iky, src, None, PTv, adjoint=True) + + ATinvdf_duT = self.Ainv[iky] * df_duT + + dA_dmT = self.getADeriv(ky, u_src, ATinvdf_duT, adjoint=True) + dRHS_dmT = self.getRHSDeriv(ky, src, ATinvdf_duT, adjoint=True) + du_dmT = -dA_dmT + dRHS_dmT + Jtv_temp1 = 1./np.pi*(df_dmT + du_dmT) + # Trapezoidal intergration + if iky==0: + #First assigment + Jtv += Jtv_temp1*dky[iky]*np.cos(ky*y) + else: + Jtv += Jtv_temp1*dky[iky]/2.*np.cos(ky*y) + Jtv += Jtv_temp0*dky[iky]/2.*np.cos(ky*y) + Jtv_temp0 = Jtv_temp1.copy() + return Utils.mkvc(Jtv) def getSourceTerm(self, ky): """ @@ -167,7 +186,6 @@ class Problem2D_CC(BaseDCProblem_2D): rho = self.curModel.rho if adjoint: return(MfRhoIDeriv( G * u ).T) * ( D.T * v) + ky**2*Utils.sdiag(u.flatten()*vol*(-1./rho**2))*v - return D * ((MfRhoIDeriv( G * u )) * v) + ky**2*Utils.sdiag(u.flatten()*vol*(-1./rho**2))*v def getRHS(self, ky): diff --git a/tests/em/static/test_DC_2D_jvecjtvecadj.py b/tests/em/static/test_DC_2D_jvecjtvecadj.py index ad7198e9..6beee640 100644 --- a/tests/em/static/test_DC_2D_jvecjtvecadj.py +++ b/tests/em/static/test_DC_2D_jvecjtvecadj.py @@ -24,14 +24,14 @@ class DCProblem_2DTestsCC(unittest.TestCase): problem = DC.Problem2D_CC(mesh, mapping=[('rho', Maps.IdentityMap(mesh))]) problem.pair(survey) - mSynth = np.ones(mesh.nC) + mSynth = np.ones(mesh.nC)*1. survey.makeSyntheticData(mSynth) # Now set up the problem to do some minimization dmis = DataMisfit.l2_DataMisfit(survey) reg = Regularization.Tikhonov(mesh) opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6) - invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4) + invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e0) inv = Inversion.BaseInversion(invProb) self.inv = inv @@ -47,21 +47,21 @@ class DCProblem_2DTestsCC(unittest.TestCase): passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3) self.assertTrue(passed) - # def test_adjoint(self): - # # Adjoint Test - # u = np.random.rand(self.mesh.nC*self.survey.nSrc) - # v = np.random.rand(self.mesh.nC) - # w = np.random.rand(self.survey.dobs.shape[0]) - # wtJv = w.dot(self.p.Jvec(self.m0, v)) - # vtJtw = v.dot(self.p.Jtvec(self.m0, w)) - # passed = np.abs(wtJv - vtJtw) < 1e-10 - # print 'Adjoint Test', np.abs(wtJv - vtJtw), passed - # self.assertTrue(passed) + def test_adjoint(self): + # Adjoint Test + u = np.random.rand(self.mesh.nC*self.survey.nSrc) + v = np.random.rand(self.mesh.nC) + w = np.random.rand(self.survey.dobs.shape[0]) + wtJv = w.dot(self.p.Jvec(self.m0, v)) + vtJtw = v.dot(self.p.Jtvec(self.m0, w)) + passed = np.abs(wtJv - vtJtw) < 1e-10 + print 'Adjoint Test', np.abs(wtJv - vtJtw), passed + self.assertTrue(passed) - # def test_dataObj(self): - # derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)] - # passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3) - # self.assertTrue(passed) + def test_dataObj(self): + derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)] + passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3) + self.assertTrue(passed) # class DCProblemTestsN(unittest.TestCase):