From cbefdd2e84483a7724dc660e45e9d9fedc45dcfb Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Tue, 15 Apr 2014 14:59:25 -0700 Subject: [PATCH 1/3] derivatives not yet working. --- simpegPF/BaseMag.py | 12 +- simpegPF/Magnetics.py | 22 +- simpegPF/Tests/test_forward_PFproblem.py | 4 +- simpegPF/Tests/test_sensitivity_PFproblem.py | 381 +++++++++---------- 4 files changed, 207 insertions(+), 212 deletions(-) diff --git a/simpegPF/BaseMag.py b/simpegPF/BaseMag.py index 5be1198c..9cef01ec 100644 --- a/simpegPF/BaseMag.py +++ b/simpegPF/BaseMag.py @@ -1,4 +1,4 @@ -from SimPEG import Model, Survey, Utils, np, sp +from SimPEG import Maps, Survey, Utils, np, sp from scipy.constants import mu_0 @@ -115,11 +115,11 @@ class MagSurveyBx(object): return bfx -class BaseMagModel(Model.BaseModel): - """BaseMagModel""" +class BaseMagMap(Maps.IdentityMap): + """BaseMagMap""" def __init__(self, mesh, **kwargs): - Model.BaseModel.__init__(self, mesh) + Maps.IdentityMap.__init__(self, mesh) def transform(self, m): @@ -129,11 +129,11 @@ class BaseMagModel(Model.BaseModel): return mu_0*sp.identity(self.nP) -class BaseDepthModel(Model.BaseModel): +class BaseDepthMap(Maps.IdentityMap): """BaseDepthMagModel""" def __init__(self, mesh, **kwargs): - Model.BaseModel.__init__(self, mesh) + Maps.IdentityMap.__init__(self, mesh) self.mesh = mesh self.active_ind = kwargs['active_ind'] self.c = kwargs['c'] diff --git a/simpegPF/Magnetics.py b/simpegPF/Magnetics.py index 51b7932a..26c67185 100644 --- a/simpegPF/Magnetics.py +++ b/simpegPF/Magnetics.py @@ -11,10 +11,10 @@ class MagneticsDiffSecondary(Problem.BaseProblem): """ surveyPair = BaseMag.BaseMagSurvey - modelPair = BaseMag.BaseMagModel + modelPair = BaseMag.BaseMagMap - def __init__(self, model, **kwargs): - Problem.BaseProblem.__init__(self, model, **kwargs) + def __init__(self, model, mapping=None, **kwargs): + Problem.BaseProblem.__init__(self, model, mapping=mapping, **kwargs) Pbc, Pin, self._Pout = \ self.mesh.getBCProjWF('neumann', discretization='CC') @@ -33,12 +33,12 @@ class MagneticsDiffSecondary(Problem.BaseProblem): def MfMu0(self): return self._MfMu0 def makeMassMatrices(self, m): - mu = self.model.transform(m) - self._MfMui = self.mesh.getFaceMass(1./mu) + mu = self.mapping.transform(m) + self._MfMui = self.mesh.getFaceInnerProduct(1./mu) # 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.getFaceMass(1/mu_0) + self._MfMu0 = self.mesh.getFaceInnerProduct(1./mu_0) # self._MfMu0 = self.mesh.getFaceInnerProduct(1/mu_0) @Utils.requires('survey') @@ -61,7 +61,7 @@ class MagneticsDiffSecondary(Problem.BaseProblem): Dface = self.mesh.faceDiv Mc = Utils.sdiag(self.mesh.vol) - mu = self.model.transform(m) + mu = self.mapping.transform(m) chi = mu/mu_0-1 Bbc, Bbc_const = CongruousMagBC(self.mesh, self.survey.B0, chi) @@ -184,8 +184,8 @@ class MagneticsDiffSecondary(Problem.BaseProblem): u = self.fields(m) B, u = u['B'], u['u'] - mu = self.model.transform(m) - dmudm = self.model.transformDeriv(m) + mu = self.mapping.transform(m) + dmudm = self.mapping.transformDeriv(m) dchidmu = Utils.sdiag(1/mu_0*np.ones(self.mesh.nC)) vol = self.mesh.vol @@ -261,8 +261,8 @@ class MagneticsDiffSecondary(Problem.BaseProblem): u = self.fields(m) B, u = u['B'], u['u'] - mu = self.model.transform(m) - dmudm = self.model.transformDeriv(m) + mu = self.mapping.transform(m) + dmudm = self.mapping.transformDeriv(m) dchidmu = Utils.sdiag(1/mu_0*np.ones(self.mesh.nC)) vol = self.mesh.vol diff --git a/simpegPF/Tests/test_forward_PFproblem.py b/simpegPF/Tests/test_forward_PFproblem.py index 6d8d5bcf..18419fda 100644 --- a/simpegPF/Tests/test_forward_PFproblem.py +++ b/simpegPF/Tests/test_forward_PFproblem.py @@ -19,8 +19,8 @@ class MagFwdProblemTests(unittest.TestCase): chi = np.ones(M.nC)*chibkg sph_ind = PF.MagAnalytics.spheremodel(M, 0., 0., 0., 100) chi[sph_ind] = chiblk - model = PF.BaseMag.BaseMagModel(M) - prob = PF.Magnetics.MagneticsDiffSecondary(model) + model = PF.BaseMag.BaseMagMap(M) + prob = PF.Magnetics.MagneticsDiffSecondary(M, mapping=model) self.prob = prob self.M = M self.chi = chi diff --git a/simpegPF/Tests/test_sensitivity_PFproblem.py b/simpegPF/Tests/test_sensitivity_PFproblem.py index 6be19f43..2d256997 100644 --- a/simpegPF/Tests/test_sensitivity_PFproblem.py +++ b/simpegPF/Tests/test_sensitivity_PFproblem.py @@ -26,7 +26,7 @@ class MagSensProblemTests(unittest.TestCase): b0 = PF.MagAnalytics.IDTtoxyz(Inc, Dec, Btot) sph_ind = PF.MagAnalytics.spheremodel(M, 0., 0., 0., 100) chi[sph_ind] = chiblk - model = PF.BaseMag.BaseMagModel(M) + model = PF.BaseMag.BaseMagMap(M) survey = BaseMag.BaseMagSurvey() @@ -40,7 +40,7 @@ class MagSensProblemTests(unittest.TestCase): - prob = PF.Magnetics.MagneticsDiffSecondary(model) + prob = PF.Magnetics.MagneticsDiffSecondary(M, mapping=model) prob.pair(survey) dpre = survey.dpred(chi) @@ -56,261 +56,256 @@ 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) def test_Jvec(self): print ">> Derivative test for Jvec" - mu = self.model.transform(self.chi) 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 - a = self.prob.Jvec(self.chi, d_chi) - - derChk = lambda m: [self.survey.dpred(m), lambda mx: self.prob.Jvec(self.chi, mx)] + derChk = lambda m: (self.survey.dpred(m), lambda v: self.prob.Jvec(m, v)) # 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_Jtvec(self): - print ">> Derivative test for Jtvec" - mu = self.model.transform(self.chi) - dobs = self.survey.dpred(self.chi) + # def test_Jtvec(self): + # print ">> Derivative test for Jtvec" + # dobs = self.survey.dpred(self.chi) - def misfit (m, dobs): - 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 .. --; - derChk = lambda m: misfit(m, dobs) - passed = Tests.checkDerivative(derChk, 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() From 1d5c3ced397e564ca7d008dd58aca38decc01830 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Tue, 15 Apr 2014 15:07:00 -0700 Subject: [PATCH 2/3] 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() From 95f0cb997278b28398eb2800a5477a76f776c2bd Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Tue, 15 Apr 2014 15:41:44 -0700 Subject: [PATCH 3/3] docs updates --- README.md | 2 +- docs/index.rst | 9 +-------- 2 files changed, 2 insertions(+), 9 deletions(-) diff --git a/README.md b/README.md index 1015642a..26ea05de 100644 --- a/README.md +++ b/README.md @@ -16,7 +16,7 @@ Tests: [https://travis-ci.org/simpeg/simpegpf](https://travis-ci.org/simpeg/simpegpf) Build Status: -[![Build Status](https://travis-ci.org/simpeg/simpegpf.png)](https://travis-ci.org/simpeg/simpegpf) +[![Build Status](https://travis-ci.org/simpeg/simpegpf.svg?branch=master)](https://travis-ci.org/simpeg/simpegpf) Bugs & Issues: [https://github.com/simpeg/simpegpf/issues](https://github.com/simpeg/simpegpf/issues) diff --git a/docs/index.rst b/docs/index.rst index bcd70591..e4a0a372 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -20,22 +20,15 @@ Potential fields api_PF - - Testing SimPEG ============== * Master Branch - .. image:: https://travis-ci.org/simpeg/simpegpf.png?branch=master + .. image:: https://travis-ci.org/simpeg/simpegpf.svg?branch=master :target: https://travis-ci.org/simpeg/simpegpf :alt: Master Branch :align: center -* Develop Branch - .. image:: https://travis-ci.org/simpeg/simpegpf.png?branch=develop - :target: https://travis-ci.org/simpeg/simpegpf - :alt: Develop Branch - :align: center Project Index & Search