mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-03 16:24:45 +08:00
Initial work Mag inversion
This commit is contained in:
SEOGI KANG
+29
-12
@@ -42,13 +42,13 @@ class BaseMagData(Data.BaseData):
|
||||
"""
|
||||
This function projects the fields onto the data space.
|
||||
|
||||
Esepcially, here for we use total magnetic intensity (TMI) data,
|
||||
which is common in practice.
|
||||
Esepcially, here for we use total magnetic intensity (TMI) data,
|
||||
which is common in practice.
|
||||
|
||||
First we project our B on to data location
|
||||
First we project our B on to data location
|
||||
|
||||
.. math::
|
||||
|
||||
|
||||
\mathbf{B}_{rec} = \mathbf{P} \mathbf{B}
|
||||
|
||||
then we take the dot product between B and b_0
|
||||
@@ -121,13 +121,30 @@ class BaseMagModel(Model.BaseModel):
|
||||
def __init__(self, mesh, **kwargs):
|
||||
Model.BaseModel.__init__(self, mesh)
|
||||
|
||||
def transform(self, m, asMu=True):
|
||||
if asMu:
|
||||
return mu_0*(1 + m)
|
||||
return m
|
||||
def transform(self, m):
|
||||
|
||||
def transformDeriv(self, m, asMu=True):
|
||||
if asMu:
|
||||
return mu_0*sp.identity(self.nP)
|
||||
return sp.identity(self.nP)
|
||||
return mu_0*(1 + m)
|
||||
|
||||
def transformDeriv(self, m):
|
||||
|
||||
return mu_0*sp.identity(self.nP)
|
||||
|
||||
class BaseDepthModel(Model.BaseModel):
|
||||
"""BaseDepthMagModel"""
|
||||
|
||||
def __init__(self, mesh, **kwargs):
|
||||
Model.BaseModel.__init__(self, mesh)
|
||||
self.mesh = mesh
|
||||
self.active_ind = kwargs['active_ind']
|
||||
|
||||
def transform(self, m):
|
||||
weight = abs(self.mesh.gridCC[:,2])**1.0
|
||||
weight = weight/weight.max()
|
||||
weight[~self.active_ind] = 1.
|
||||
return m*weight
|
||||
|
||||
def transformDeriv(self, m):
|
||||
weight = abs(self.mesh.gridCC[:,2])**1.0
|
||||
weight = weight/weight.max()
|
||||
weight[~self.active_ind] = 1.
|
||||
return Utils.sdiag(weight)
|
||||
|
||||
+102
-65
@@ -1,4 +1,4 @@
|
||||
from SimPEG import Mesh, Problem, Utils, np, sp, Tests
|
||||
from SimPEG import *
|
||||
import BaseMag
|
||||
from scipy.constants import mu_0
|
||||
from MagAnalytics import spheremodel, CongruousMagBC
|
||||
@@ -51,7 +51,7 @@ class MagneticsDiffSecondary(Problem.BaseProblem):
|
||||
|
||||
def getRHS(self, m):
|
||||
"""
|
||||
|
||||
|
||||
.. math ::
|
||||
|
||||
\mathbf{rhs} = \Div(\MfMui)^{-1}\mathbf{M}^f_{\mu_0^{-1}}\mathbf{B}_0 - \Div\mathbf{B}_0+\diag(v)\mathbf{D} \mathbf{P}_{out}^T \mathbf{B}_{sBC}
|
||||
@@ -61,7 +61,9 @@ class MagneticsDiffSecondary(Problem.BaseProblem):
|
||||
Dface = self.mesh.faceDiv
|
||||
Mc = Utils.sdiag(self.mesh.vol)
|
||||
|
||||
chi = self.model.transform(m, asMu=False)
|
||||
mu = self.model.transform(m)
|
||||
chi = mu/mu_0-1
|
||||
|
||||
Bbc, Bbc_const = CongruousMagBC(self.mesh, self.data.B0, chi)
|
||||
self.Bbc = Bbc
|
||||
self.Bbc_const = Bbc_const
|
||||
@@ -105,11 +107,11 @@ class MagneticsDiffSecondary(Problem.BaseProblem):
|
||||
return {'B': B, 'u': u}
|
||||
|
||||
@Utils.timeIt
|
||||
def Jvec(self, m, v, fields=None):
|
||||
def Jvec(self, m, v, u=None):
|
||||
"""
|
||||
Computing Jacobian multiplied by vector
|
||||
|
||||
By setting our problem as
|
||||
|
||||
By setting our problem as
|
||||
|
||||
.. math ::
|
||||
|
||||
@@ -124,51 +126,51 @@ class MagneticsDiffSecondary(Problem.BaseProblem):
|
||||
|
||||
\\frac{\delta \mathbf{u}}{\delta \mathbf{m}} = - [\\nabla_u \mathbf{C}(\mathbf{u})]^{-1}\\nabla_m \mathbf{C}(\mathbf{m})
|
||||
|
||||
With some linear algebra we can have
|
||||
|
||||
With some linear algebra we can have
|
||||
|
||||
.. math ::
|
||||
|
||||
\\nabla_u \mathbf{C}(\mathbf{u}) = \mathbf{A}
|
||||
|
||||
\\nabla_m \mathbf{C}(\mathbf{m}) =
|
||||
\\frac{\partial \mathbf{A}}{\partial \mathbf{m}}(\mathbf{m})\mathbf{u} - \\frac{\partial \mathbf{rhs}(\mathbf{m})}{\partial \mathbf{m}}
|
||||
\\nabla_m \mathbf{C}(\mathbf{m}) =
|
||||
\\frac{\partial \mathbf{A}}{\partial \mathbf{m}}(\mathbf{m})\mathbf{u} - \\frac{\partial \mathbf{rhs}(\mathbf{m})}{\partial \mathbf{m}}
|
||||
|
||||
.. math ::
|
||||
.. math ::
|
||||
|
||||
\\frac{\partial \mathbf{A}}{\partial \mathbf{m}}(\mathbf{m})\mathbf{u} =
|
||||
\\frac{\partial \mathbf{A}}{\partial \mathbf{m}}(\mathbf{m})\mathbf{u} =
|
||||
\\frac{\partial \mathbf{\mu}}{\partial \mathbf{m}} \left[\Div \diag (\Div^T \mathbf{u}) \dMfMuI \\right]
|
||||
|
||||
\dMfMuI = \diag(\MfMui)^{-1}_{vec} \mathbf{Av}_{F2CC}^T\diag(\mathbf{v})\diag(\\frac{1}{\mu^2})
|
||||
|
||||
\\frac{\partial \mathbf{rhs}(\mathbf{m})}{\partial \mathbf{m}} = \\frac{\partial \mathbf{\mu}}{\partial \mathbf{m}} \left[
|
||||
\\frac{\partial \mathbf{rhs}(\mathbf{m})}{\partial \mathbf{m}} = \\frac{\partial \mathbf{\mu}}{\partial \mathbf{m}} \left[
|
||||
\Div \diag(\M^f_{\mu_{0}^{-1}}\mathbf{B}_0) \dMfMuI \\right] - \diag(\mathbf{v})\mathbf{D} \mathbf{P}_{out}^T\\frac{\partial B_{sBC}}{\partial \mathbf{m}}
|
||||
|
||||
In the end,
|
||||
In the end,
|
||||
|
||||
.. math ::
|
||||
|
||||
\\frac{\delta \mathbf{u}}{\delta \mathbf{m}} =
|
||||
- [ \mathbf{A} ]^{-1}\left[ \\frac{\partial \mathbf{A}}{\partial \mathbf{m}}(\mathbf{m})\mathbf{u}
|
||||
\\frac{\delta \mathbf{u}}{\delta \mathbf{m}} =
|
||||
- [ \mathbf{A} ]^{-1}\left[ \\frac{\partial \mathbf{A}}{\partial \mathbf{m}}(\mathbf{m})\mathbf{u}
|
||||
- \\frac{\partial \mathbf{rhs}(\mathbf{m})}{\partial \mathbf{m}} \\right]
|
||||
|
||||
A little tricky point here is we are not interested in potential (u), but interested in magnetic flux (B).
|
||||
Thus, we need sensitivity for B. Now we take derivative of B w.r.t m and have
|
||||
|
||||
A little tricky point here is we are not interested in potential (u), but interested in magnetic flux (B).
|
||||
Thus, we need sensitivity for B. Now we take derivative of B w.r.t m and have
|
||||
|
||||
.. math ::
|
||||
|
||||
\\frac{\delta \mathbf{B}} {\delta \mathbf{m}} = \\frac{\partial \mathbf{\mu} } {\partial \mathbf{m} }
|
||||
\left[
|
||||
\\frac{\delta \mathbf{B}} {\delta \mathbf{m}} = \\frac{\partial \mathbf{\mu} } {\partial \mathbf{m} }
|
||||
\left[
|
||||
\diag(\M^f_{\mu_{0}^{-1} } \mathbf{B}_0) \dMfMuI \\
|
||||
- \diag (\Div^T\mathbf{u})\dMfMuI
|
||||
\\right ]
|
||||
|
||||
- (\MfMui)^{-1}\Div^T\\frac{\delta\mathbf{u}}{\delta \mathbf{m}}
|
||||
|
||||
Finally we evaluate the above, but we should remember that
|
||||
- (\MfMui)^{-1}\Div^T\\frac{\delta\mathbf{u}}{\delta \mathbf{m}}
|
||||
|
||||
.. note ::
|
||||
Finally we evaluate the above, but we should remember that
|
||||
|
||||
We only want to evalute
|
||||
.. note ::
|
||||
|
||||
We only want to evalute
|
||||
|
||||
.. math ::
|
||||
|
||||
@@ -178,13 +180,13 @@ class MagneticsDiffSecondary(Problem.BaseProblem):
|
||||
|
||||
|
||||
"""
|
||||
if fields is None:
|
||||
fields = self.fields(m)
|
||||
if u is None:
|
||||
u = self.fields(m)
|
||||
|
||||
B, u = fields['B'], fields['u']
|
||||
mu = self.model.transform(m, asMu=True)
|
||||
dmudm = self.model.transformDeriv(m, asMu=True)
|
||||
dmdmu = Utils.sdiag(1/(dmudm.diagonal()))
|
||||
B, u = u['B'], u['u']
|
||||
mu = self.model.transform(m)
|
||||
dmudm = self.model.transformDeriv(m)
|
||||
dchidmu = Utils.sdiag(1/mu_0*np.ones(self.mesh.nC))
|
||||
|
||||
vol = self.mesh.vol
|
||||
Div = self._Div
|
||||
@@ -204,42 +206,43 @@ class MagneticsDiffSecondary(Problem.BaseProblem):
|
||||
dCdm_A = Div * ( Utils.sdiag( Div.T * u )* dMfMuI *dmudm )
|
||||
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, v)
|
||||
dCdm_RHSv = dCdm_RHS1*(dmudm*v) + dCdm_RHS2v
|
||||
dCdm_RHS2v = (Utils.sdiag(vol)*temp1)*np.inner(vol, dchidmu*dmudm*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)
|
||||
m1 = sp.linalg.interface.aslinearoperator(Utils.sdiag(1/dCdu.diagonal()))
|
||||
sol, info = sp.linalg.bicgstab(dCdu, dCdm_v, tol=1e-6, maxiter=1000, M=m1)
|
||||
|
||||
if info > 0:
|
||||
raise Exception ("Iterative solver did not work well")
|
||||
print "Iterative solver did not work well (Jvec)"
|
||||
# raise Exception ("Iterative solver did not work well")
|
||||
|
||||
# B = self.MfMuI*self.MfMu0*B0-B0-self.MfMuI*self._Div.T*u
|
||||
# dBdm = d\mudm*dBd\mu
|
||||
|
||||
dudm = -sol
|
||||
dudm = -sol
|
||||
dBdmv = ( Utils.sdiag(self.MfMu0*B0)*(dMfMuI * (dmudm*v)) \
|
||||
- Utils.sdiag(Div.T*u)*(dMfMuI* (dmudm*v)) \
|
||||
- self.MfMuI*(Div.T* (dudm)) )
|
||||
|
||||
- self.MfMuI*(Div.T* (dudm)) )
|
||||
|
||||
return Utils.mkvc(P*dBdmv)
|
||||
|
||||
@Utils.timeIt
|
||||
def Jtvec(self, m, v, fields=None):
|
||||
def Jtvec(self, m, v, u=None):
|
||||
"""
|
||||
Computing Jacobian^T multiplied by vector.
|
||||
|
||||
.. math ::
|
||||
|
||||
(\\frac{\delta \mathbf{P}\mathbf{B}} {\delta \mathbf{m}})^{T} = \left[ \mathbf{P}_{deriv}\\frac{\partial \mathbf{\mu} } {\partial \mathbf{m} }
|
||||
\left[
|
||||
(\\frac{\delta \mathbf{P}\mathbf{B}} {\delta \mathbf{m}})^{T} = \left[ \mathbf{P}_{deriv}\\frac{\partial \mathbf{\mu} } {\partial \mathbf{m} }
|
||||
\left[
|
||||
\diag(\M^f_{\mu_{0}^{-1} } \mathbf{B}_0) \dMfMuI \\
|
||||
- \diag (\Div^T\mathbf{u})\dMfMuI
|
||||
\\right ]\\right]^{T}
|
||||
\\right ]\\right]^{T}
|
||||
|
||||
- \left[\mathbf{P}_{deriv}(\MfMui)^{-1}\Div^T\\frac{\delta\mathbf{u}}{\delta \mathbf{m}} \\right]^{T}
|
||||
|
||||
where
|
||||
where
|
||||
|
||||
.. math ::
|
||||
|
||||
@@ -247,20 +250,20 @@ class MagneticsDiffSecondary(Problem.BaseProblem):
|
||||
|
||||
.. note ::
|
||||
|
||||
Here we only want to compute
|
||||
Here we only want to compute
|
||||
|
||||
.. math ::
|
||||
|
||||
\mathbf{J}^{T}\mathbf{v} = (\\frac{\delta \mathbf{P}\mathbf{B}} {\delta \mathbf{m}})^{T} \mathbf{v}
|
||||
|
||||
"""
|
||||
if fields is None:
|
||||
fields = self.fields(m)
|
||||
if u is None:
|
||||
u = self.fields(m)
|
||||
|
||||
B, u = fields['B'], fields['u']
|
||||
mu = self.model.transform(m, asMu=True)
|
||||
dmudm = self.model.transformDeriv(m, asMu=True)
|
||||
dmdmu = Utils.sdiag(1/(dmudm.diagonal()))
|
||||
B, u = u['B'], u['u']
|
||||
mu = self.model.transform(m)
|
||||
dmudm = self.model.transformDeriv(m)
|
||||
dchidmu = Utils.sdiag(1/mu_0*np.ones(self.mesh.nC))
|
||||
|
||||
vol = self.mesh.vol
|
||||
Div = self._Div
|
||||
@@ -279,26 +282,33 @@ class MagneticsDiffSecondary(Problem.BaseProblem):
|
||||
dCdu = self.getA(m)
|
||||
s = Div * ( self.MfMuI.T * ( P.T*v ) )
|
||||
|
||||
m1 = sp.linalg.interface.aslinearoperator(Utils.sdiag(1/(dCdu.T).diagonal()))
|
||||
sol, info = sp.linalg.bicgstab(dCdu.T, s, tol=1e-8, maxiter=1000, M=m1)
|
||||
m1 = sp.linalg.interface.aslinearoperator(Utils.sdiag(1/(dCdu.T).diagonal()))
|
||||
sol, info = sp.linalg.bicgstab(dCdu.T, s, tol=1e-6, maxiter=1000, M=m1)
|
||||
|
||||
if info > 0:
|
||||
raise Exception ("Iterative solver did not work well")
|
||||
print "Iterative solver did not work well (Jtvec)"
|
||||
# raise Exception ("Iterative solver did not work well")
|
||||
|
||||
|
||||
# dCdm_A = Div * ( Utils.sdiag( Div.T * u )* dMfMuI *dmudm )
|
||||
dCdm_Atsol = ( dMfMuI.T*( Utils.sdiag( Div.T * u ) * (Div.T * dmudm)) ) * sol
|
||||
# dCdm_Atsol = ( dMfMuI.T*( Utils.sdiag( Div.T * u ) * (Div.T * dmudm)) ) * sol
|
||||
dCdm_Atsol = ( dmudm.T * dMfMuI.T*( Utils.sdiag( Div.T * u ) * Div.T ) ) * sol
|
||||
|
||||
# dCdm_RHS1 = Div * (Utils.sdiag( self.MfMu0*B0 ) * dMfMuI)
|
||||
# dCdm_RHS1tsol = (dMfMuI.T*( Utils.sdiag( self.MfMu0*B0 ) ) * Div.T * dmudm) * sol
|
||||
dCdm_RHS1tsol = ( dmudm.T * dMfMuI.T*( Utils.sdiag( self.MfMu0*B0 ) ) * Div.T ) * sol
|
||||
|
||||
|
||||
# dCdm_RHS1 = Div * (Utils.sdiag( self.MfMu0*B0 ) * dMfMuI)
|
||||
dCdm_RHS1tsol = (dMfMuI.T*( Utils.sdiag( self.MfMu0*B0 ) ) * Div.T * dmudm) * sol
|
||||
|
||||
# temp1 = (Dface*(self._Pout.T*self.Bbc_const*self.Bbc))
|
||||
# dCdm_RHS2v = (Utils.sdiag(vol)*temp1)*np.inner(vol, v)
|
||||
temp1sol = ( Dface.T*( Utils.sdiag(vol)*sol ) )
|
||||
temp2 = self.Bbc_const*(self._Pout.T*self.Bbc).T
|
||||
dCdm_RHS2tsol = vol*np.inner(temp2, temp1sol)
|
||||
# dCdm_RHS2v = (Utils.sdiag(vol)*temp1)*np.inner(vol, dchidmu*dmudm*v)
|
||||
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
|
||||
|
||||
# dCdm_RHSv = dCdm_RHS1*(dmudm*v) + dCdm_RHS2v
|
||||
# dCdm_RHSv = dCdm_RHS1*(dmudm*v) + dCdm_RHS2v
|
||||
# dCdm_v = dCdm_A*v - dCdm_RHSv
|
||||
|
||||
Ctv = dCdm_Atsol - dCdm_RHStsol
|
||||
@@ -311,7 +321,35 @@ class MagneticsDiffSecondary(Problem.BaseProblem):
|
||||
Btemp = Utils.sdiag(Div.T*u)*(dMfMuI* (dmudm))
|
||||
Jtv = Atemp.T*(P.T*v) - Btemp.T*(P.T*v) - Ctv
|
||||
|
||||
return Utils.mkvc(Jtv)
|
||||
return Utils.mkvc(Jtv)
|
||||
|
||||
def MagneticsDiffSecondaryInv(mesh, model, data, **kwargs):
|
||||
|
||||
"""
|
||||
Inversion module for MagneticsDiffSecondary
|
||||
|
||||
"""
|
||||
from SimPEG import Optimization, Regularization, Parameters, ObjFunction, Inversion
|
||||
prob = MagneticsDiffSecondary(mesh, model)
|
||||
|
||||
if prob.ispaired:
|
||||
prob.unpair()
|
||||
if data.ispaired:
|
||||
data.unpair()
|
||||
prob.pair(data)
|
||||
|
||||
# Create an optimization program
|
||||
opt = Optimization.InexactGaussNewton(maxIter=5)
|
||||
opt.bfgsH0 = Solver(sp.identity(model.nP),flag='D')
|
||||
# Create a regularization program
|
||||
reg = Regularization.Tikhonov(model)
|
||||
# Create an objective function
|
||||
beta = Parameters.BetaSchedule(beta0=1e0)
|
||||
obj = ObjFunction.BaseObjFunction(data, reg, beta=beta)
|
||||
# Create an inversion object
|
||||
inv = Inversion.BaseInversion(obj, opt)
|
||||
|
||||
return inv, reg
|
||||
|
||||
|
||||
|
||||
@@ -337,7 +375,7 @@ if __name__ == '__main__':
|
||||
Btot = 51000
|
||||
|
||||
data.setBackgroundField(Inc, Dec, Btot)
|
||||
|
||||
|
||||
xr = np.linspace(-300, 300, 41)
|
||||
yr = np.linspace(-300, 300, 41)
|
||||
X, Y = np.meshgrid(xr, yr)
|
||||
@@ -346,7 +384,6 @@ if __name__ == '__main__':
|
||||
data.rxLoc = rxLoc
|
||||
|
||||
prob = MagneticsDiffSecondary(mesh, model)
|
||||
|
||||
prob.pair(data)
|
||||
|
||||
dpred = data.dpred(chi)
|
||||
|
||||
@@ -30,11 +30,11 @@ class MagFwdProblemTests(unittest.TestCase):
|
||||
|
||||
data = PF.BaseMag.BaseMagData()
|
||||
|
||||
Inc = 90.
|
||||
Dec = 0.
|
||||
Inc = 45.
|
||||
Dec = 45.
|
||||
Btot = 51000
|
||||
|
||||
b0 = PF.MagAnalytics.IDTtoxyz(Inc, Dec, Btot)
|
||||
b0 = PF.MagAnalytics.IDTtoxyz(Inc, Dec, Btot)
|
||||
data.setBackgroundField(Inc, Dec, Btot)
|
||||
xr = np.linspace(-300, 300, 41)
|
||||
yr = np.linspace(-300, 300, 41)
|
||||
@@ -42,21 +42,25 @@ class MagFwdProblemTests(unittest.TestCase):
|
||||
Z = np.ones((xr.size, yr.size))*150
|
||||
rxLoc = np.c_[Utils.mkvc(X), Utils.mkvc(Y), Utils.mkvc(Z)]
|
||||
data.rxLoc = rxLoc
|
||||
|
||||
|
||||
self.prob.pair(data)
|
||||
u = self.prob.fields(self.chi)
|
||||
B = u['B']
|
||||
|
||||
bxa,bya,bza = PF.MagAnalytics.MagSphereAnalFunA(rxLoc[:,0],rxLoc[:,1],rxLoc[:,2],100.,0.,0.,0.,0.01, b0,'secondary')
|
||||
|
||||
bxa,bya,bza = PF.MagAnalytics.MagSphereAnalFunA(rxLoc[:,0],rxLoc[:,1],rxLoc[:,2],100.,0.,0.,0.,0.01, b0,'secondary')
|
||||
|
||||
dpred = data.projectFieldsAsVector(B)
|
||||
err = np.linalg.norm(dpred-np.r_[bxa, bya, bza])/np.linalg.norm(np.r_[bxa, bya, bza])
|
||||
|
||||
plt.plot(dpred)
|
||||
plt.plot(np.r_[bxa, bya, bza])
|
||||
plt.show()
|
||||
|
||||
if err > 0.05:
|
||||
raise Exception('Anaytic test is failed T.T')
|
||||
else:
|
||||
print "Anaytic test is passed"
|
||||
pass
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
|
||||
@@ -18,7 +18,7 @@ class MagSensProblemTests(unittest.TestCase):
|
||||
chibkg = 0.001
|
||||
chiblk = 0.01
|
||||
chi = np.ones(M.nC)*chibkg
|
||||
|
||||
|
||||
Inc = 90.
|
||||
Dec = 0.
|
||||
Btot = 51000
|
||||
@@ -58,7 +58,7 @@ class MagSensProblemTests(unittest.TestCase):
|
||||
|
||||
def test_mass(self):
|
||||
print '\n >>Derivative for MfMuI works.'
|
||||
mu = self.model.transform(self.chi, asMu=True)
|
||||
mu = self.model.transform(self.chi)
|
||||
def MfmuI(mu):
|
||||
|
||||
chi = mu/mu_0-1
|
||||
@@ -70,36 +70,36 @@ class MagSensProblemTests(unittest.TestCase):
|
||||
return MfMuI
|
||||
|
||||
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)
|
||||
dMfMuI = Utils.sdiag(MfMuI**2)*aveF2CC.T*Utils.sdiag(vol*1./mu**2)
|
||||
|
||||
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)
|
||||
|
||||
|
||||
|
||||
|
||||
self.assertTrue(passed)
|
||||
|
||||
|
||||
|
||||
def test_dCdm_Av(self):
|
||||
print '\n >>Derivative for Cm_A.'
|
||||
Div = self.prob._Div
|
||||
Div = self.prob._Div
|
||||
vol = self.prob.mesh.vol
|
||||
aveF2CC = self.prob.mesh.aveF2CC
|
||||
|
||||
def Cm_A(chi):
|
||||
dmudm = self.model.transformDeriv(chi, asMu=True)
|
||||
dmudm = self.model.transformDeriv(chi)
|
||||
u = self.u
|
||||
# chi = mu/mu_0-1
|
||||
self.prob.makeMassMatrices(chi)
|
||||
mu = self.model.transform(self.chi, asMu=True)
|
||||
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)
|
||||
@@ -109,11 +109,11 @@ class MagSensProblemTests(unittest.TestCase):
|
||||
return Cm_A
|
||||
|
||||
def dCdm_A(chi, v):
|
||||
|
||||
dmudm = self.model.transformDeriv(chi, asMu=True)
|
||||
|
||||
dmudm = self.model.transformDeriv(chi)
|
||||
u = self.u
|
||||
self.prob.makeMassMatrices(chi)
|
||||
mu = self.model.transform(self.chi, asMu=True)
|
||||
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)
|
||||
@@ -133,7 +133,7 @@ class MagSensProblemTests(unittest.TestCase):
|
||||
print '\n >>Derivative for Cm_RHS.'
|
||||
u = self.u
|
||||
Div = self.prob._Div
|
||||
mu = self.model.transform(self.chi, asMu=True)
|
||||
mu = self.model.transform(self.chi)
|
||||
vol = self.prob.mesh.vol
|
||||
Mc = Utils.sdiag(vol)
|
||||
aveF2CC = self.prob.mesh.aveF2CC
|
||||
@@ -143,45 +143,45 @@ class MagSensProblemTests(unittest.TestCase):
|
||||
def Cm_RHS(chi):
|
||||
|
||||
self.prob.makeMassMatrices(chi)
|
||||
dmudm = self.model.transformDeriv(chi, asMu=True)
|
||||
dmudm = self.model.transformDeriv(chi)
|
||||
dchidmu = Utils.sdiag(1/(dmudm.diagonal()))
|
||||
Bbc, Bbc_const = PF.MagAnalytics.CongruousMagBC(self.prob.mesh, self.data.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
|
||||
RHS1 = Div*self.prob.MfMuI*self.prob.MfMu0*B0
|
||||
RHS2 = Mc*Dface*self.prob._Pout.T*Bbc
|
||||
RHS = RHS1 + RHS2 + Div*B0
|
||||
|
||||
return RHS
|
||||
|
||||
|
||||
def dCdm_RHS(chi, v):
|
||||
|
||||
|
||||
|
||||
self.prob.makeMassMatrices(chi)
|
||||
dmudm = self.model.transformDeriv(chi, asMu=True)
|
||||
dmudm = self.model.transformDeriv(chi)
|
||||
dmdmu = Utils.sdiag(1/(dmudm.diagonal()))
|
||||
Bbc, Bbc_const = PF.MagAnalytics.CongruousMagBC(self.prob.mesh, self.data.B0, chi)
|
||||
MfMuIvec = 1/self.prob.MfMui.diagonal()
|
||||
dMfMuI = Utils.sdiag(MfMuIvec**2)*aveF2CC.T*Utils.sdiag(vol*1./mu**2)
|
||||
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
|
||||
|
||||
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)
|
||||
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, asMu=True)
|
||||
mu = self.model.transform(self.chi)
|
||||
vol = self.prob.mesh.vol
|
||||
Mc = Utils.sdiag(vol)
|
||||
aveF2CC = self.prob.mesh.aveF2CC
|
||||
@@ -193,11 +193,11 @@ class MagSensProblemTests(unittest.TestCase):
|
||||
return u
|
||||
|
||||
def dudm(chi, v):
|
||||
|
||||
|
||||
chi = mu/mu_0-1
|
||||
self.prob.makeMassMatrices(chi)
|
||||
self.prob.makeMassMatrices(chi)
|
||||
u = self.u
|
||||
dmudm = self.model.transformDeriv(chi, asMu=True)
|
||||
dmudm = self.model.transformDeriv(chi)
|
||||
dmdmu = Utils.sdiag(1/(dmudm.diagonal()))
|
||||
Bbc, Bbc_const = PF.MagAnalytics.CongruousMagBC(self.prob.mesh, self.data.B0, chi)
|
||||
MfMuIvec = 1/self.prob.MfMui.diagonal()
|
||||
@@ -207,21 +207,21 @@ class MagSensProblemTests(unittest.TestCase):
|
||||
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_RHSv = dCdm_RHS1*(dmudm*v) + dCdm_RHS2v
|
||||
dCdm_v = dCdm_A*v - dCdm_RHSv
|
||||
m1 = sp.linalg.interface.aslinearoperator(Utils.sdiag(1/dCdu.diagonal()))
|
||||
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
|
||||
|
||||
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[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 .. --;
|
||||
# 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)
|
||||
|
||||
@@ -230,7 +230,7 @@ class MagSensProblemTests(unittest.TestCase):
|
||||
print ">> Derivative test for dBdm"
|
||||
u = self.u
|
||||
Div = self.prob._Div
|
||||
mu = self.model.transform(self.chi, asMu=True)
|
||||
mu = self.model.transform(self.chi)
|
||||
vol = self.prob.mesh.vol
|
||||
Mc = Utils.sdiag(vol)
|
||||
aveF2CC = self.prob.mesh.aveF2CC
|
||||
@@ -242,11 +242,11 @@ class MagSensProblemTests(unittest.TestCase):
|
||||
return B
|
||||
|
||||
def dBdm(chi, v):
|
||||
|
||||
|
||||
chi = mu/mu_0-1
|
||||
self.prob.makeMassMatrices(chi)
|
||||
self.prob.makeMassMatrices(chi)
|
||||
u = self.u
|
||||
dmudm = self.model.transformDeriv(chi, asMu=True)
|
||||
dmudm = self.model.transformDeriv(chi)
|
||||
dmdmu = Utils.sdiag(1/(dmudm.diagonal()))
|
||||
Bbc, Bbc_const = PF.MagAnalytics.CongruousMagBC(self.prob.mesh, self.data.B0, chi)
|
||||
MfMuIvec = 1/self.prob.MfMui.diagonal()
|
||||
@@ -256,47 +256,47 @@ class MagSensProblemTests(unittest.TestCase):
|
||||
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_RHSv = dCdm_RHS1*(dmudm*v) + dCdm_RHS2v
|
||||
dCdm_v = dCdm_A*v - dCdm_RHSv
|
||||
m1 = sp.linalg.interface.aslinearoperator(Utils.sdiag(1/dCdu.diagonal()))
|
||||
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
|
||||
dBdmv = ( Utils.sdiag(self.prob.MfMu0*B0)*(dMfMuI * (dmudm*v)) \
|
||||
- Utils.sdiag(Div.T*u)*(dMfMuI* (dmudm*v)) \
|
||||
- self.prob.MfMuI*(Div.T* (dudm)) )
|
||||
- self.prob.MfMuI*(Div.T* (dudm)) )
|
||||
|
||||
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[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=3, dx = d_chi, plotIt=False)
|
||||
# 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, asMu=True)
|
||||
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
|
||||
d_chi[d_sph_ind] = 0.1
|
||||
|
||||
a = self.prob.Jvec(self.chi, d_chi)
|
||||
|
||||
derChk = lambda m: [self.data.dpred(m), lambda mx: self.prob.Jvec(self.chi, mx)]
|
||||
# TODO: I am not sure why the order get worse as step decreases .. --;
|
||||
# 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, asMu=True)
|
||||
mu = self.model.transform(self.chi)
|
||||
dobs = self.data.dpred(self.chi)
|
||||
|
||||
def misfit (m, dobs):
|
||||
@@ -307,10 +307,10 @@ class MagSensProblemTests(unittest.TestCase):
|
||||
|
||||
return misfit, dmisfit
|
||||
|
||||
# TODO: I am not sure why the order get worse as step decreases .. --;
|
||||
# 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)
|
||||
self.assertTrue(passed)
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
|
||||
File diff suppressed because one or more lines are too long
File diff suppressed because one or more lines are too long
@@ -1,142 +0,0 @@
|
||||
|
||||
from SimPEG.Utils.matutils import kron3, speye, sdiag
|
||||
from SimPEG import *
|
||||
import numpy as np
|
||||
import scipy.sparse as sp
|
||||
def ddxFaceDivBC(n, bc):
|
||||
|
||||
ij = (np.array([0, n-1]),np.array([0, 1]))
|
||||
vals = np.zeros(2)
|
||||
|
||||
# Set the first side
|
||||
if(bc[0] == 'dirichlet'):
|
||||
vals[0] = 0
|
||||
elif(bc[0] == 'neumann'):
|
||||
vals[0] = -1
|
||||
# Set the second side
|
||||
if(bc[1] == 'dirichlet'):
|
||||
vals[1] = 0
|
||||
elif(bc[1] == 'neumann'):
|
||||
vals[1] = 1
|
||||
D = sp.csr_matrix((vals, ij), shape=(n,2))
|
||||
return D
|
||||
|
||||
|
||||
def faceDivBC(mesh, BC, ind):
|
||||
"""
|
||||
The facd divergence boundary condtion matrix
|
||||
|
||||
.. math::
|
||||
|
||||
|
||||
|
||||
"""
|
||||
# The number of cell centers in each direction
|
||||
n = mesh.nCv
|
||||
# Compute faceDivergence operator on faces
|
||||
if(mesh.dim == 1):
|
||||
D = ddxFaceDivBC(n[0], BC[0])
|
||||
elif(mesh.dim == 2):
|
||||
D1 = sp.kron(speye(n[1]), ddxFaceDivBC(n[0]), BC[0])
|
||||
D2 = sp.kron(ddxFaceDivBC(n[1], BC[1]), speye(n[0]))
|
||||
D = sp.hstack((D1, D2), format="csr")
|
||||
elif(mesh.dim == 3):
|
||||
D1 = kron3(speye(n[2]), speye(n[1]), ddxFaceDivBC(n[0], BC[0]))
|
||||
D2 = kron3(speye(n[2]), ddxFaceDivBC(n[1], BC[1]), speye(n[0]))
|
||||
D3 = kron3(ddxFaceDivBC(n[2], BC[2]), speye(n[1]), speye(n[0]))
|
||||
D = sp.hstack((D1, D2, D3), format="csr")
|
||||
# Compute areas of cell faces & volumes
|
||||
S = mesh.area[ind]
|
||||
V = mesh.vol
|
||||
mesh._faceDiv = sdiag(1/V)*D*sdiag(S)
|
||||
|
||||
return mesh._faceDiv
|
||||
|
||||
|
||||
def faceBCind(mesh):
|
||||
"""
|
||||
Find indices of boundary faces in each direction
|
||||
|
||||
"""
|
||||
if(mesh.dim==1):
|
||||
indxd = (mesh.gridFx[:,0]==min(mesh.gridFx[:,0]))
|
||||
indxu = (mesh.gridFx[:,0]==max(mesh.gridFx[:,0]))
|
||||
return indxd, indxu
|
||||
elif(mesh.dim==1):
|
||||
indxd = (mesh.gridFx[:,0]==min(mesh.gridFx[:,0]))
|
||||
indxu = (mesh.gridFx[:,0]==max(mesh.gridFx[:,0]))
|
||||
indyd = (mesh.gridFy[:,1]==min(mesh.gridFy[:,1]))
|
||||
indyu = (mesh.gridFy[:,1]==max(mesh.gridFy[:,1]))
|
||||
return indxd, indxu, indyd, indyu
|
||||
elif(mesh.dim==3):
|
||||
indxd = (mesh.gridFx[:,0]==min(mesh.gridFx[:,0]))
|
||||
indxu = (mesh.gridFx[:,0]==max(mesh.gridFx[:,0]))
|
||||
indyd = (mesh.gridFy[:,1]==min(mesh.gridFy[:,1]))
|
||||
indyu = (mesh.gridFy[:,1]==max(mesh.gridFy[:,1]))
|
||||
indzd = (mesh.gridFz[:,2]==min(mesh.gridFz[:,2]))
|
||||
indzu = (mesh.gridFz[:,2]==max(mesh.gridFz[:,2]))
|
||||
return indxd, indxu, indyd, indyu, indzd, indzu
|
||||
|
||||
|
||||
def spheremodel(mesh, x0, y0, z0, r):
|
||||
"""
|
||||
Generate model indicies for sphere
|
||||
- (x0, y0, z0 ): is the center location of sphere
|
||||
- r: is the radius of the sphere
|
||||
- it returns logical indicies of cell-center model
|
||||
"""
|
||||
ind = np.sqrt((mesh.gridCC[:,0]-x0)**2+(mesh.gridCC[:,1]-y0)**2+(mesh.gridCC[:,2]-z0)**2 ) < r
|
||||
return ind
|
||||
|
||||
|
||||
|
||||
def MagSphereAnalFun(x, y, z, R, x0, y0, z0, mu1, mu2, H0, flag):
|
||||
"""
|
||||
Analytic function for Magnetics problem. The set up here is
|
||||
magnetic sphere in whole-space.
|
||||
- (x0,y0,z0)
|
||||
- (x0, y0, z0 ): is the center location of sphere
|
||||
- r: is the radius of the sphere
|
||||
|
||||
.. math::
|
||||
|
||||
\mathbf{H}^p = H_0\hat{x}
|
||||
|
||||
|
||||
"""
|
||||
if (~np.size(x)==np.size(y)==np.size(z)):
|
||||
print "Specify same size of x, y, z"
|
||||
return
|
||||
dim = x.shape
|
||||
x = Utils.mkvc(x)
|
||||
y = Utils.mkvc(y)
|
||||
z = Utils.mkvc(z)
|
||||
|
||||
ind = np.sqrt((x-x0)**2+(y-y0)**2+(z-z0)**2 ) < R
|
||||
r = Utils.mkvc(np.sqrt((x-x0)**2+(y-y0)**2+(z-z0)**2 ))
|
||||
Bx = np.zeros(x.size)
|
||||
By = np.zeros(x.size)
|
||||
Bz = np.zeros(x.size)
|
||||
|
||||
# Inside of the sphere
|
||||
rf2 = 3*mu1/(mu2+2*mu1)
|
||||
if (flag == 'total'):
|
||||
Bx[ind] = mu2*H0*(rf2)
|
||||
elif (flag == 'secondary'):
|
||||
Bx[ind] = mu2*H0*(rf2)-mu1*H0
|
||||
|
||||
By[ind] = 0.
|
||||
Bz[ind] = 0.
|
||||
# Outside of the sphere
|
||||
rf1 = (mu2-mu1)/(mu2+2*mu1)
|
||||
if (flag == 'total'):
|
||||
Bx[~ind] = mu1*(H0+H0/r[~ind]**5*(R**3)*rf1*(2*x[~ind]**2-y[~ind]**2-z[~ind]**2))
|
||||
elif (flag == 'secondary'):
|
||||
Bx[~ind] = mu1*(H0/r[~ind]**5*(R**3)*rf1*(2*x[~ind]**2-y[~ind]**2-z[~ind]**2))
|
||||
|
||||
By[~ind] = mu1*(H0/r[~ind]**5*(R**3)*rf1*(3*x[~ind]*y[~ind]))
|
||||
Bz[~ind] = mu1*(H0/r[~ind]**5*(R**3)*rf1*(3*x[~ind]*z[~ind]))
|
||||
|
||||
return np.reshape(Bx, x.shape, order='F'), np.reshape(By, x.shape, order='F'), np.reshape(Bz, x.shape, order='F')
|
||||
|
||||
|
||||
Reference in New Issue
Block a user