wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-06-28 17:03:04 +08:00
Code Issues Packages Projects Releases Wiki Activity

Modifications for inversion

Browse Source
This commit is contained in:
seogi
2014-07-04 14:34:39 -07:00
parent 83a7a7bc12
commit d43caf092b
7 changed files with 2412 additions and 5275 deletions
Download Patch File Download Diff File Expand all files Collapse all files
simpegPF/BaseMag.py
+10 -17
View File
@@ -121,31 +121,24 @@ class BaseMagMap(Maps.IdentityMap):
def __init__(self, mesh, **kwargs):
Maps.IdentityMap.__init__(self, mesh)
def transform(self, m):
def _transform(self, m):
return mu_0*(1 + m)
def transformDeriv(self, m):
def deriv(self, m):
return mu_0*sp.identity(self.nP)
class BaseDepthMap(Maps.IdentityMap):
"""BaseDepthMagModel"""
class WeightMap(Maps.IdentityMap):
"""Weighted Map for distributed parameters"""
def __init__(self, mesh, **kwargs):
def __init__(self, mesh, weight, **kwargs):
Maps.IdentityMap.__init__(self, mesh)
self.mesh = mesh
self.active_ind = kwargs['active_ind']
self.c = kwargs['c']
self.weight = weight
def transform(self, m):
weight = abs(self.mesh.gridCC[:,2])**self.c
weight = weight/weight.max()
weight[~self.active_ind] = 1.
return m*weight
def _transform(self, m):
return m*self.weight
def transformDeriv(self, m):
weight = abs(self.mesh.gridCC[:,2])**self.c
weight = weight/weight.max()
weight[~self.active_ind] = 1.
return Utils.sdiag(weight)
def deriv(self, m):
return Utils.sdiag(self.weight)
simpegPF/Magnetics.py
+17 -9
View File
@@ -33,7 +33,7 @@ class MagneticsDiffSecondary(Problem.BaseProblem):
def MfMu0(self): return self._MfMu0
def makeMassMatrices(self, m):
mu = self.mapping.transform(m)
mu = self.mapping*m
self._MfMui = self.mesh.getFaceInnerProduct(1./mu)/self.mesh.dim
# self._MfMui = self.mesh.getFaceInnerProduct(1./mu)
#TODO: this will break if tensor mu
@@ -61,13 +61,16 @@ class MagneticsDiffSecondary(Problem.BaseProblem):
Dface = self.mesh.faceDiv
Mc = Utils.sdiag(self.mesh.vol)
mu = self.mapping.transform(m)
mu = self.mapping*m
chi = mu/mu_0-1
#temporary fix
Bbc, Bbc_const = CongruousMagBC(self.mesh, self.survey.B0, chi)
self.Bbc = Bbc
self.Bbc_const = Bbc_const
return self._Div*self.MfMuI*self.MfMu0*B0 - self._Div*B0 + Mc*Dface*self._Pout.T*Bbc
# return self._Div*self.MfMuI*self.MfMu0*B0 - self._Div*B0 + Mc*Dface*self._Pout.T*Bbc
return self._Div*self.MfMuI*self.MfMu0*B0 - self._Div*B0
def getA(self, m):
"""
@@ -184,8 +187,8 @@ class MagneticsDiffSecondary(Problem.BaseProblem):
u = self.fields(m)
B, u = u['B'], u['u']
mu = self.mapping.transform(m)
dmudm = self.mapping.transformDeriv(m)
mu = self.mapping*(m)
dmudm = self.mapping.deriv(m)
dchidmu = Utils.sdiag(1/mu_0*np.ones(self.mesh.nC))
vol = self.mesh.vol
@@ -207,7 +210,9 @@ class MagneticsDiffSecondary(Problem.BaseProblem):
dCdm_RHS1 = Div * (Utils.sdiag( self.MfMu0*B0 ) * dMfMuI)
temp1 = (Dface*(self._Pout.T*self.Bbc_const*self.Bbc))
dCdm_RHS2v = (Utils.sdiag(vol)*temp1)*np.inner(vol, dchidmu*dmudm*v)
dCdm_RHSv = dCdm_RHS1*(dmudm*v) + dCdm_RHS2v
#temporary fix
# dCdm_RHSv = dCdm_RHS1*(dmudm*v) + dCdm_RHS2v
dCdm_RHSv = dCdm_RHS1*(dmudm*v)
dCdm_v = dCdm_A*v - dCdm_RHSv
m1 = sp.linalg.interface.aslinearoperator(Utils.sdiag(1/dCdu.diagonal()))
@@ -261,8 +266,8 @@ class MagneticsDiffSecondary(Problem.BaseProblem):
u = self.fields(m)
B, u = u['B'], u['u']
mu = self.mapping.transform(m)
dmudm = self.mapping.transformDeriv(m)
mu = self.mapping*(m)
dmudm = self.mapping.deriv(m)
dchidmu = Utils.sdiag(1/mu_0*np.ones(self.mesh.nC))
vol = self.mesh.vol
@@ -306,7 +311,10 @@ class MagneticsDiffSecondary(Problem.BaseProblem):
dCdm_RHS2tsol = (dmudm.T*dchidmu.T*vol)*np.inner(temp2, temp1sol)
# dCdm_RHSv = dCdm_RHS1*(dmudm*v) + dCdm_RHS2v
dCdm_RHStsol = dCdm_RHS1tsol - dCdm_RHS2tsol
#temporary fix
# dCdm_RHStsol = dCdm_RHS1tsol - dCdm_RHS2tsol
dCdm_RHStsol = dCdm_RHS1tsol
# dCdm_RHSv = dCdm_RHS1*(dmudm*v) + dCdm_RHS2v
# dCdm_v = dCdm_A*v - dCdm_RHSv
simpegPF/Tests/test_sensitivity_PFproblem.py
+91 -91
View File
@@ -57,7 +57,7 @@ class MagSensProblemTests(unittest.TestCase):
def test_mass(self):
print '\n >>Derivative for MfMuI works.'
mu = self.model.transform(self.chi)
mu = self.model*self.chi
def MfmuI(mu):
chi = mu/mu_0-1
@@ -94,11 +94,11 @@ class MagSensProblemTests(unittest.TestCase):
aveF2CC = self.prob.mesh.aveF2CC
def Cm_A(chi):
dmudm = self.model.transformDeriv(chi)
dmudm = self.model.deriv(chi)
u = self.u
# chi = mu/mu_0-1
self.prob.makeMassMatrices(chi)
mu = self.model.transform(self.chi)
mu = self.model*(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)
@@ -109,10 +109,10 @@ class MagSensProblemTests(unittest.TestCase):
def dCdm_A(chi, v):
dmudm = self.model.transformDeriv(chi)
dmudm = self.model.deriv(chi)
u = self.u
self.prob.makeMassMatrices(chi)
mu = self.model.transform(self.chi)
mu = self.model*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)
@@ -132,7 +132,7 @@ class MagSensProblemTests(unittest.TestCase):
print '\n >>Derivative for Cm_RHS.'
u = self.u
Div = self.prob._Div
mu = self.model.transform(self.chi)
mu = self.model*self.chi
vol = self.prob.mesh.vol
Mc = Utils.sdiag(vol)
aveF2CC = self.prob.mesh.aveF2CC
@@ -142,7 +142,7 @@ class MagSensProblemTests(unittest.TestCase):
def Cm_RHS(chi):
self.prob.makeMassMatrices(chi)
dmudm = self.model.transformDeriv(chi)
dmudm = self.model.deriv(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()
@@ -158,7 +158,7 @@ class MagSensProblemTests(unittest.TestCase):
self.prob.makeMassMatrices(chi)
dmudm = self.model.transformDeriv(chi)
dmudm = self.model.deriv(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()
@@ -176,105 +176,105 @@ class MagSensProblemTests(unittest.TestCase):
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*(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.deriv(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*(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.deriv(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)
simpegPF/notebooks/Jacobian.ipynb
-402
View File
File diff suppressed because one or more lines are too long
simpegPF/notebooks/Linear Inversion.ipynb
-546
View File
File diff suppressed because one or more lines are too long
simpegPF/notebooks/MagInversion.ipynb
+2294 -3873
View File
File diff suppressed because one or more lines are too long
simpegPF/notebooks/Play_analytic_fun.ipynb
-337
View File
File diff suppressed because one or more lines are too long