From 1d5c3ced397e564ca7d008dd58aca38decc01830 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Tue, 15 Apr 2014 15:07:00 -0700 Subject: [PATCH] Fix derivatives. --- simpegPF/Magnetics.py | 4 +- simpegPF/Tests/test_sensitivity_PFproblem.py | 370 +++++++++---------- 2 files changed, 187 insertions(+), 187 deletions(-) diff --git a/simpegPF/Magnetics.py b/simpegPF/Magnetics.py index 26c67185..a19a6c82 100644 --- a/simpegPF/Magnetics.py +++ b/simpegPF/Magnetics.py @@ -34,11 +34,11 @@ class MagneticsDiffSecondary(Problem.BaseProblem): def makeMassMatrices(self, m): mu = self.mapping.transform(m) - self._MfMui = self.mesh.getFaceInnerProduct(1./mu) + self._MfMui = self.mesh.getFaceInnerProduct(1./mu)/self.mesh.dim # self._MfMui = self.mesh.getFaceInnerProduct(1./mu) #TODO: this will break if tensor mu self._MfMuI = Utils.sdiag(1./self._MfMui.diagonal()) - self._MfMu0 = self.mesh.getFaceInnerProduct(1./mu_0) + self._MfMu0 = self.mesh.getFaceInnerProduct(1./mu_0)/self.mesh.dim # self._MfMu0 = self.mesh.getFaceInnerProduct(1/mu_0) @Utils.requires('survey') diff --git a/simpegPF/Tests/test_sensitivity_PFproblem.py b/simpegPF/Tests/test_sensitivity_PFproblem.py index 2d256997..63676579 100644 --- a/simpegPF/Tests/test_sensitivity_PFproblem.py +++ b/simpegPF/Tests/test_sensitivity_PFproblem.py @@ -56,226 +56,226 @@ class MagSensProblemTests(unittest.TestCase): - # def test_mass(self): - # print '\n >>Derivative for MfMuI works.' - # mu = self.model.transform(self.chi) - # def MfmuI(mu): + def test_mass(self): + print '\n >>Derivative for MfMuI works.' + mu = self.model.transform(self.chi) + def MfmuI(mu): - # chi = mu/mu_0-1 - # self.prob.makeMassMatrices(chi) - # vol = self.prob.mesh.vol - # aveF2CC = self.prob.mesh.aveF2CC - # MfMuI = self.prob.MfMuI.diagonal() + chi = mu/mu_0-1 + self.prob.makeMassMatrices(chi) + vol = self.prob.mesh.vol + aveF2CC = self.prob.mesh.aveF2CC + MfMuI = self.prob.MfMuI.diagonal() - # return MfMuI + return MfMuI - # def dMfmuI(mu, v): + def dMfmuI(mu, v): - # chi = mu/mu_0-1 - # self.prob.makeMassMatrices(chi) - # vol = self.prob.mesh.vol - # aveF2CC = self.prob.mesh.aveF2CC - # MfMuI = self.prob.MfMuI.diagonal() - # dMfMuI = Utils.sdiag(MfMuI**2)*aveF2CC.T*Utils.sdiag(vol*1./mu**2) + chi = mu/mu_0-1 + self.prob.makeMassMatrices(chi) + vol = self.prob.mesh.vol + aveF2CC = self.prob.mesh.aveF2CC + MfMuI = self.prob.MfMuI.diagonal() + dMfMuI = Utils.sdiag(MfMuI**2)*aveF2CC.T*Utils.sdiag(vol*1./mu**2) - # return dMfMuI*v + return dMfMuI*v - # d_mu = mu*0.8 - # derChk = lambda m: [MfmuI(m), lambda mx: dMfmuI(self.chi, mx)] - # passed = Tests.checkDerivative(derChk, mu, num=4, dx = d_mu, plotIt=False) + d_mu = mu*0.8 + derChk = lambda m: [MfmuI(m), lambda mx: dMfmuI(self.chi, mx)] + passed = Tests.checkDerivative(derChk, mu, num=4, dx = d_mu, plotIt=False) - # self.assertTrue(passed) + self.assertTrue(passed) - # def test_dCdm_Av(self): - # print '\n >>Derivative for Cm_A.' - # Div = self.prob._Div - # vol = self.prob.mesh.vol - # aveF2CC = self.prob.mesh.aveF2CC + def test_dCdm_Av(self): + print '\n >>Derivative for Cm_A.' + Div = self.prob._Div + vol = self.prob.mesh.vol + aveF2CC = self.prob.mesh.aveF2CC - # def Cm_A(chi): - # dmudm = self.model.transformDeriv(chi) - # u = self.u - # # chi = mu/mu_0-1 - # self.prob.makeMassMatrices(chi) - # mu = self.model.transform(self.chi) - # A = self.prob.getA(self.chi) - # MfMuIvec = 1/self.prob.MfMui.diagonal() - # dMfMuI = Utils.sdiag(MfMuIvec**2)*aveF2CC.T*Utils.sdiag(vol*1./mu**2) + def Cm_A(chi): + dmudm = self.model.transformDeriv(chi) + u = self.u + # chi = mu/mu_0-1 + self.prob.makeMassMatrices(chi) + mu = self.model.transform(self.chi) + A = self.prob.getA(self.chi) + MfMuIvec = 1/self.prob.MfMui.diagonal() + dMfMuI = Utils.sdiag(MfMuIvec**2)*aveF2CC.T*Utils.sdiag(vol*1./mu**2) - # Cm_A = A*u + Cm_A = A*u - # return Cm_A + return Cm_A - # def dCdm_A(chi, v): + def dCdm_A(chi, v): - # dmudm = self.model.transformDeriv(chi) - # u = self.u - # self.prob.makeMassMatrices(chi) - # mu = self.model.transform(self.chi) - # A = self.prob.getA(self.chi) - # MfMuIvec = 1/self.prob.MfMui.diagonal() - # dMfMuI = Utils.sdiag(MfMuIvec**2)*aveF2CC.T*Utils.sdiag(vol*1./mu**2) + dmudm = self.model.transformDeriv(chi) + u = self.u + self.prob.makeMassMatrices(chi) + mu = self.model.transform(self.chi) + A = self.prob.getA(self.chi) + MfMuIvec = 1/self.prob.MfMui.diagonal() + dMfMuI = Utils.sdiag(MfMuIvec**2)*aveF2CC.T*Utils.sdiag(vol*1./mu**2) - # Cm_A = A*u - # dCdm_A = Div * ( Utils.sdiag( Div.T * u )* dMfMuI *dmudm ) + Cm_A = A*u + dCdm_A = Div * ( Utils.sdiag( Div.T * u )* dMfMuI *dmudm ) - # return dCdm_A*v + return dCdm_A*v - # d_chi = self.chi*0.8 - # derChk = lambda m: [Cm_A(m), lambda mx: dCdm_A(self.chi, mx)] - # passed = Tests.checkDerivative(derChk, self.chi, num=4, dx = d_chi, plotIt=False) - # self.assertTrue(passed) + d_chi = self.chi*0.8 + derChk = lambda m: [Cm_A(m), lambda mx: dCdm_A(self.chi, mx)] + passed = Tests.checkDerivative(derChk, self.chi, num=4, dx = d_chi, plotIt=False) + self.assertTrue(passed) - # def test_dCdmu_RHS(self): - # print '\n >>Derivative for Cm_RHS.' - # u = self.u - # Div = self.prob._Div - # mu = self.model.transform(self.chi) - # vol = self.prob.mesh.vol - # Mc = Utils.sdiag(vol) - # aveF2CC = self.prob.mesh.aveF2CC - # B0 = self.prob.getB0() - # Dface = self.prob.mesh.faceDiv + def test_dCdmu_RHS(self): + print '\n >>Derivative for Cm_RHS.' + u = self.u + Div = self.prob._Div + mu = self.model.transform(self.chi) + vol = self.prob.mesh.vol + Mc = Utils.sdiag(vol) + aveF2CC = self.prob.mesh.aveF2CC + B0 = self.prob.getB0() + Dface = self.prob.mesh.faceDiv - # def Cm_RHS(chi): + def Cm_RHS(chi): - # self.prob.makeMassMatrices(chi) - # dmudm = self.model.transformDeriv(chi) - # dchidmu = Utils.sdiag(1/(dmudm.diagonal())) - # Bbc, Bbc_const = PF.MagAnalytics.CongruousMagBC(self.prob.mesh, self.survey.B0, chi) - # MfMuIvec = 1/self.prob.MfMui.diagonal() - # dMfMuI = Utils.sdiag(MfMuIvec**2)*aveF2CC.T*Utils.sdiag(vol*1./mu**2) - # RHS1 = Div*self.prob.MfMuI*self.prob.MfMu0*B0 - # RHS2 = Mc*Dface*self.prob._Pout.T*Bbc - # RHS = RHS1 + RHS2 + Div*B0 + self.prob.makeMassMatrices(chi) + dmudm = self.model.transformDeriv(chi) + dchidmu = Utils.sdiag(1/(dmudm.diagonal())) + Bbc, Bbc_const = PF.MagAnalytics.CongruousMagBC(self.prob.mesh, self.survey.B0, chi) + MfMuIvec = 1/self.prob.MfMui.diagonal() + dMfMuI = Utils.sdiag(MfMuIvec**2)*aveF2CC.T*Utils.sdiag(vol*1./mu**2) + RHS1 = Div*self.prob.MfMuI*self.prob.MfMu0*B0 + RHS2 = Mc*Dface*self.prob._Pout.T*Bbc + RHS = RHS1 + RHS2 + Div*B0 - # return RHS + return RHS - # def dCdm_RHS(chi, v): + def dCdm_RHS(chi, v): - # self.prob.makeMassMatrices(chi) - # dmudm = self.model.transformDeriv(chi) - # dmdmu = Utils.sdiag(1/(dmudm.diagonal())) - # Bbc, Bbc_const = PF.MagAnalytics.CongruousMagBC(self.prob.mesh, self.survey.B0, chi) - # MfMuIvec = 1/self.prob.MfMui.diagonal() - # dMfMuI = Utils.sdiag(MfMuIvec**2)*aveF2CC.T*Utils.sdiag(vol*1./mu**2) - # dCdm_RHS1 = Div * (Utils.sdiag( self.prob.MfMu0*B0 ) * dMfMuI) - # temp1 = (Dface*(self.prob._Pout.T*Bbc_const*Bbc)) - # dCdm_RHS2v = (Utils.sdiag(vol)*temp1)*np.inner(vol, v) - # dCdm_RHSv = dCdm_RHS1*(dmudm*v) + dCdm_RHS2v + self.prob.makeMassMatrices(chi) + dmudm = self.model.transformDeriv(chi) + dmdmu = Utils.sdiag(1/(dmudm.diagonal())) + Bbc, Bbc_const = PF.MagAnalytics.CongruousMagBC(self.prob.mesh, self.survey.B0, chi) + MfMuIvec = 1/self.prob.MfMui.diagonal() + dMfMuI = Utils.sdiag(MfMuIvec**2)*aveF2CC.T*Utils.sdiag(vol*1./mu**2) + dCdm_RHS1 = Div * (Utils.sdiag( self.prob.MfMu0*B0 ) * dMfMuI) + temp1 = (Dface*(self.prob._Pout.T*Bbc_const*Bbc)) + dCdm_RHS2v = (Utils.sdiag(vol)*temp1)*np.inner(vol, v) + dCdm_RHSv = dCdm_RHS1*(dmudm*v) + dCdm_RHS2v - # return dCdm_RHSv + return dCdm_RHSv - # d_chi = self.chi*0.8 - # derChk = lambda m: [Cm_RHS(m), lambda mx: dCdm_RHS(self.chi, mx)] - # passed = Tests.checkDerivative(derChk, self.chi, num=4, dx = d_chi, plotIt=False) - # self.assertTrue(passed) + d_chi = self.chi*0.8 + derChk = lambda m: [Cm_RHS(m), lambda mx: dCdm_RHS(self.chi, mx)] + passed = Tests.checkDerivative(derChk, self.chi, num=4, dx = d_chi, plotIt=False) + self.assertTrue(passed) - # def test_dudm(self): - # print ">> Derivative test for dudm" - # u = self.u - # Div = self.prob._Div - # mu = self.model.transform(self.chi) - # vol = self.prob.mesh.vol - # Mc = Utils.sdiag(vol) - # aveF2CC = self.prob.mesh.aveF2CC - # B0 = self.prob.getB0() - # Dface = self.prob.mesh.faceDiv + def test_dudm(self): + print ">> Derivative test for dudm" + u = self.u + Div = self.prob._Div + mu = self.model.transform(self.chi) + vol = self.prob.mesh.vol + Mc = Utils.sdiag(vol) + aveF2CC = self.prob.mesh.aveF2CC + B0 = self.prob.getB0() + Dface = self.prob.mesh.faceDiv - # def ufun(chi): - # u = self.prob.fields(chi)['u'] - # return u + def ufun(chi): + u = self.prob.fields(chi)['u'] + return u - # def dudm(chi, v): + def dudm(chi, v): - # chi = mu/mu_0-1 - # self.prob.makeMassMatrices(chi) - # u = self.u - # dmudm = self.model.transformDeriv(chi) - # dmdmu = Utils.sdiag(1/(dmudm.diagonal())) - # Bbc, Bbc_const = PF.MagAnalytics.CongruousMagBC(self.prob.mesh, self.survey.B0, chi) - # MfMuIvec = 1/self.prob.MfMui.diagonal() - # dMfMuI = Utils.sdiag(MfMuIvec**2)*aveF2CC.T*Utils.sdiag(vol*1./mu**2) - # dCdu = self.prob.getA(chi) - # dCdm_A = Div * ( Utils.sdiag( Div.T * u )* dMfMuI *dmudm ) - # dCdm_RHS1 = Div * (Utils.sdiag( self.prob.MfMu0*B0 ) * dMfMuI) - # temp1 = (Dface*(self.prob._Pout.T*Bbc_const*Bbc)) - # dCdm_RHS2v = (Utils.sdiag(vol)*temp1)*np.inner(vol, v) - # dCdm_RHSv = dCdm_RHS1*(dmudm*v) + dCdm_RHS2v - # dCdm_v = dCdm_A*v - dCdm_RHSv - # m1 = sp.linalg.interface.aslinearoperator(Utils.sdiag(1/dCdu.diagonal())) - # sol, info = sp.linalg.bicgstab(dCdu, dCdm_v, tol=1e-8, maxiter=1000, M=m1) + chi = mu/mu_0-1 + self.prob.makeMassMatrices(chi) + u = self.u + dmudm = self.model.transformDeriv(chi) + dmdmu = Utils.sdiag(1/(dmudm.diagonal())) + Bbc, Bbc_const = PF.MagAnalytics.CongruousMagBC(self.prob.mesh, self.survey.B0, chi) + MfMuIvec = 1/self.prob.MfMui.diagonal() + dMfMuI = Utils.sdiag(MfMuIvec**2)*aveF2CC.T*Utils.sdiag(vol*1./mu**2) + dCdu = self.prob.getA(chi) + dCdm_A = Div * ( Utils.sdiag( Div.T * u )* dMfMuI *dmudm ) + dCdm_RHS1 = Div * (Utils.sdiag( self.prob.MfMu0*B0 ) * dMfMuI) + temp1 = (Dface*(self.prob._Pout.T*Bbc_const*Bbc)) + dCdm_RHS2v = (Utils.sdiag(vol)*temp1)*np.inner(vol, v) + dCdm_RHSv = dCdm_RHS1*(dmudm*v) + dCdm_RHS2v + dCdm_v = dCdm_A*v - dCdm_RHSv + m1 = sp.linalg.interface.aslinearoperator(Utils.sdiag(1/dCdu.diagonal())) + sol, info = sp.linalg.bicgstab(dCdu, dCdm_v, tol=1e-8, maxiter=1000, M=m1) - # dudm = -sol + dudm = -sol - # return dudm + return dudm - # d_chi = 10.0*self.chi #np.random.rand(mesh.nCz) - # d_sph_ind = PF.MagAnalytics.spheremodel(self.prob.mesh, 0., 0., -50., 50) - # d_chi[d_sph_ind] = 0.1 + d_chi = 10.0*self.chi #np.random.rand(mesh.nCz) + d_sph_ind = PF.MagAnalytics.spheremodel(self.prob.mesh, 0., 0., -50., 50) + d_chi[d_sph_ind] = 0.1 - # derChk = lambda m: [ufun(m), lambda mx: dudm(self.chi, mx)] - # # TODO: I am not sure why the order get worse as step decreases .. --; - # passed = Tests.checkDerivative(derChk, self.chi, num=2, dx = d_chi, plotIt=False) - # self.assertTrue(passed) + derChk = lambda m: [ufun(m), lambda mx: dudm(self.chi, mx)] + # TODO: I am not sure why the order get worse as step decreases .. --; + passed = Tests.checkDerivative(derChk, self.chi, num=2, dx = d_chi, plotIt=False) + self.assertTrue(passed) - # def test_dBdm(self): - # print ">> Derivative test for dBdm" - # u = self.u - # Div = self.prob._Div - # mu = self.model.transform(self.chi) - # vol = self.prob.mesh.vol - # Mc = Utils.sdiag(vol) - # aveF2CC = self.prob.mesh.aveF2CC - # B0 = self.prob.getB0() - # Dface = self.prob.mesh.faceDiv + def test_dBdm(self): + print ">> Derivative test for dBdm" + u = self.u + Div = self.prob._Div + mu = self.model.transform(self.chi) + vol = self.prob.mesh.vol + Mc = Utils.sdiag(vol) + aveF2CC = self.prob.mesh.aveF2CC + B0 = self.prob.getB0() + Dface = self.prob.mesh.faceDiv - # def Bfun(chi): - # B = self.prob.fields(chi)['B'] - # return B + def Bfun(chi): + B = self.prob.fields(chi)['B'] + return B - # def dBdm(chi, v): + def dBdm(chi, v): - # chi = mu/mu_0-1 - # self.prob.makeMassMatrices(chi) - # u = self.u - # dmudm = self.model.transformDeriv(chi) - # dmdmu = Utils.sdiag(1/(dmudm.diagonal())) - # Bbc, Bbc_const = PF.MagAnalytics.CongruousMagBC(self.prob.mesh, self.survey.B0, chi) - # MfMuIvec = 1/self.prob.MfMui.diagonal() - # dMfMuI = Utils.sdiag(MfMuIvec**2)*aveF2CC.T*Utils.sdiag(vol*1./mu**2) - # dCdu = self.prob.getA(chi) - # dCdm_A = Div * ( Utils.sdiag( Div.T * u )* dMfMuI *dmudm ) - # dCdm_RHS1 = Div * (Utils.sdiag( self.prob.MfMu0*B0 ) * dMfMuI) - # temp1 = (Dface*(self.prob._Pout.T*Bbc_const*Bbc)) - # dCdm_RHS2v = (Utils.sdiag(vol)*temp1)*np.inner(vol, v) - # dCdm_RHSv = dCdm_RHS1*(dmudm*v) + dCdm_RHS2v - # dCdm_v = dCdm_A*v - dCdm_RHSv - # m1 = sp.linalg.interface.aslinearoperator(Utils.sdiag(1/dCdu.diagonal())) - # sol, info = sp.linalg.bicgstab(dCdu, dCdm_v, tol=1e-8, maxiter=1000, M=m1) + chi = mu/mu_0-1 + self.prob.makeMassMatrices(chi) + u = self.u + dmudm = self.model.transformDeriv(chi) + dmdmu = Utils.sdiag(1/(dmudm.diagonal())) + Bbc, Bbc_const = PF.MagAnalytics.CongruousMagBC(self.prob.mesh, self.survey.B0, chi) + MfMuIvec = 1/self.prob.MfMui.diagonal() + dMfMuI = Utils.sdiag(MfMuIvec**2)*aveF2CC.T*Utils.sdiag(vol*1./mu**2) + dCdu = self.prob.getA(chi) + dCdm_A = Div * ( Utils.sdiag( Div.T * u )* dMfMuI *dmudm ) + dCdm_RHS1 = Div * (Utils.sdiag( self.prob.MfMu0*B0 ) * dMfMuI) + temp1 = (Dface*(self.prob._Pout.T*Bbc_const*Bbc)) + dCdm_RHS2v = (Utils.sdiag(vol)*temp1)*np.inner(vol, v) + dCdm_RHSv = dCdm_RHS1*(dmudm*v) + dCdm_RHS2v + dCdm_v = dCdm_A*v - dCdm_RHSv + m1 = sp.linalg.interface.aslinearoperator(Utils.sdiag(1/dCdu.diagonal())) + sol, info = sp.linalg.bicgstab(dCdu, dCdm_v, tol=1e-8, maxiter=1000, M=m1) - # dudm = -sol - # dBdmv = ( Utils.sdiag(self.prob.MfMu0*B0)*(dMfMuI * (dmudm*v)) \ - # - Utils.sdiag(Div.T*u)*(dMfMuI* (dmudm*v)) \ - # - self.prob.MfMuI*(Div.T* (dudm)) ) + dudm = -sol + dBdmv = ( Utils.sdiag(self.prob.MfMu0*B0)*(dMfMuI * (dmudm*v)) \ + - Utils.sdiag(Div.T*u)*(dMfMuI* (dmudm*v)) \ + - self.prob.MfMuI*(Div.T* (dudm)) ) - # return dBdmv + return dBdmv - # d_chi = 10.0*self.chi #np.random.rand(mesh.nCz) - # d_sph_ind = PF.MagAnalytics.spheremodel(self.prob.mesh, 0., 0., -50., 50) - # d_chi[d_sph_ind] = 0.1 + d_chi = 10.0*self.chi #np.random.rand(mesh.nCz) + d_sph_ind = PF.MagAnalytics.spheremodel(self.prob.mesh, 0., 0., -50., 50) + d_chi[d_sph_ind] = 0.1 - # derChk = lambda m: [Bfun(m), lambda mx: dBdm(self.chi, mx)] - # # TODO: I am not sure why the order get worse as step decreases .. --; - # passed = Tests.checkDerivative(derChk, self.chi, num=2, dx = d_chi, plotIt=False) - # self.assertTrue(passed) + derChk = lambda m: [Bfun(m), lambda mx: dBdm(self.chi, mx)] + # TODO: I am not sure why the order get worse as step decreases .. --; + passed = Tests.checkDerivative(derChk, self.chi, num=2, dx = d_chi, plotIt=False) + self.assertTrue(passed) @@ -291,21 +291,21 @@ class MagSensProblemTests(unittest.TestCase): passed = Tests.checkDerivative(derChk, self.chi, num=2, dx = d_chi, plotIt=False) self.assertTrue(passed) - # def test_Jtvec(self): - # print ">> Derivative test for Jtvec" - # dobs = self.survey.dpred(self.chi) + def test_Jtvec(self): + print ">> Derivative test for Jtvec" + dobs = self.survey.dpred(self.chi) - # def misfit(m): - # dpre = self.survey.dpred(m) - # misfit = 0.5*np.linalg.norm(dpre-dobs)**2 - # residual = dpre-dobs - # dmisfit = self.prob.Jtvec(self.chi, residual) + def misfit(m): + dpre = self.survey.dpred(m) + misfit = 0.5*np.linalg.norm(dpre-dobs)**2 + residual = dpre-dobs + dmisfit = self.prob.Jtvec(self.chi, residual) - # return misfit, dmisfit + return misfit, dmisfit - # # TODO: I am not sure why the order get worse as step decreases .. --; - # passed = Tests.checkDerivative(misfit, self.chi, num=4, plotIt=False) - # self.assertTrue(passed) + # TODO: I am not sure why the order get worse as step decreases .. --; + passed = Tests.checkDerivative(misfit, self.chi, num=4, plotIt=False) + self.assertTrue(passed) if __name__ == '__main__': unittest.main()