mirror of
https://github.com/wassname/simpeg.git
synced 2026-08-13 12:40:08 +08:00
resolved conflicts
This commit is contained in:
seogi
+47
-16
@@ -11,10 +11,14 @@ class BaseMagData(Data.BaseData):
|
||||
def __init__(self, **kwargs):
|
||||
Data.BaseData.__init__(self, **kwargs)
|
||||
|
||||
#TODO: change to inc, dec, intensity
|
||||
def setBackgroundField(self, x=1., y=0., z=0.):
|
||||
# Primary field in x-direction (background)
|
||||
self.B0 = np.r_[x,y,z]
|
||||
|
||||
def setBackgroundField(self, Inc, Dec, Btot):
|
||||
|
||||
Bx = Btot*np.cos(Inc/180.*np.pi)*np.sin(Dec/180.*np.pi)
|
||||
By = Btot*np.cos(Inc/180.*np.pi)*np.cos(Dec/180.*np.pi)
|
||||
Bz = -Btot*np.sin(Inc/180.*np.pi)
|
||||
|
||||
self.B0 = np.r_[Bx,By,Bz]
|
||||
|
||||
@property
|
||||
def Qfx(self):
|
||||
@@ -34,36 +38,63 @@ class BaseMagData(Data.BaseData):
|
||||
self._Qfz = self.prob.mesh.getInterpolationMat(self.rxLoc,'Fz')
|
||||
return self._Qfz
|
||||
|
||||
def projectFields(self, B):
|
||||
def projectFields(self, u):
|
||||
"""
|
||||
This function projects the fields onto the data space.
|
||||
|
||||
Esepcially, here for we use total magnetic intensity (TMI) data,
|
||||
which is common in practice.
|
||||
|
||||
First we project our B on to data location
|
||||
|
||||
.. math::
|
||||
d_\\text{pred} = \mathbf{P} u(m)
|
||||
|
||||
\mathbf{B}_{rec} = \mathbf{P} \mathbf{B}
|
||||
|
||||
then we take the dot product between B and b_0
|
||||
|
||||
.. math ::
|
||||
|
||||
\\text{TMI} = \\vec{B}_s \cdot \hat{B}_0
|
||||
|
||||
"""
|
||||
#TODO: There can be some different tyes of data like |B| or B
|
||||
|
||||
# bfx = self.Qfx*B
|
||||
# bfy = self.Qfy*B
|
||||
bfz = self.Qfz*B
|
||||
return bfz
|
||||
bfx = self.Qfx*u['B']
|
||||
bfy = self.Qfy*u['B']
|
||||
bfz = self.Qfz*u['B']
|
||||
|
||||
# return np.sqrt(bfx**2 + bfy**2 + bfz**2)
|
||||
|
||||
# Generate unit vector
|
||||
B0 = self.prob.data.B0
|
||||
Bot = np.sqrt(B0[0]**2+B0[1]**2+B0[2]**2)
|
||||
box = B0[0]/Bot
|
||||
boy = B0[1]/Bot
|
||||
boz = B0[2]/Bot
|
||||
|
||||
# return bfx*box + bfx*boy + bfx*boz
|
||||
return bfx*box + bfy*boy + bfz*boz
|
||||
|
||||
# return np.sqrt(bfx**2 + bfy**2 + bfz**2)
|
||||
|
||||
@Utils.count
|
||||
def projectFieldsDeriv(self, B):
|
||||
"""
|
||||
This function projects the fields onto the data space.
|
||||
|
||||
|
||||
.. math::
|
||||
|
||||
\\frac{\partial d_\\text{pred}}{\partial u} = \mathbf{P}
|
||||
\\frac{\partial d_\\text{pred}}{\partial \mathbf{B}} = \mathbf{P}
|
||||
|
||||
Esepcially, this function is for TMI data type
|
||||
|
||||
"""
|
||||
return self.Qfz
|
||||
# Generate unit vector
|
||||
B0 = self.prob.data.B0
|
||||
Bot = np.sqrt(B0[0]**2+B0[1]**2+B0[2]**2)
|
||||
box = B0[0]/Bot
|
||||
boy = B0[1]/Bot
|
||||
boz = B0[2]/Bot
|
||||
|
||||
return self.Qfx*box+self.Qfy*boy+self.Qfz*boz
|
||||
|
||||
|
||||
def projectFieldsAsVector(self, B):
|
||||
|
||||
+195
-92
@@ -41,6 +41,7 @@ class MagneticsDiffSecondary(Problem.BaseProblem):
|
||||
self._MfMu0 = self.mesh.getFaceMass(1/mu_0)
|
||||
# self._MfMu0 = self.mesh.getFaceInnerProduct(1/mu_0)
|
||||
|
||||
@Utils.requires('data')
|
||||
def getB0(self):
|
||||
b0 = self.data.B0
|
||||
B0 = np.r_[b0[0]*np.ones(self.mesh.nFx),
|
||||
@@ -49,17 +50,21 @@ class MagneticsDiffSecondary(Problem.BaseProblem):
|
||||
return B0
|
||||
|
||||
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}
|
||||
|
||||
"""
|
||||
B0 = self.getB0()
|
||||
Dface = self.mesh.faceDiv
|
||||
Mc = Utils.sdiag(self.mesh.vol)
|
||||
|
||||
chi = self.model.transform(m, asMu=False)
|
||||
Bbc, const = CongruousMagBC(self.mesh, self.data.B0, chi)
|
||||
self.Bbc_const = const
|
||||
Bbc, Bbc_const = CongruousMagBC(self.mesh, self.data.B0, chi)
|
||||
self.Bbc = Bbc
|
||||
#TODO: put congrous BC back in
|
||||
# return self._Div*self.MfMuI*self.MfMu0*B0 - self._Div*B0 #+ Mc*Dface*self._Pout.T*Bbc
|
||||
self.Bbc_const = Bbc_const
|
||||
return self._Div*self.MfMuI*self.MfMu0*B0 - self._Div*B0 + Mc*Dface*self._Pout.T*Bbc
|
||||
|
||||
def getA(self, m):
|
||||
@@ -68,41 +73,39 @@ class MagneticsDiffSecondary(Problem.BaseProblem):
|
||||
|
||||
The A matrix has the form:
|
||||
|
||||
.. math ::
|
||||
|
||||
\mathbf{A}\mathbf{u} = \mathbf{rhs}
|
||||
|
||||
\mathbf{A} = - \Div(\MfMui)^{-1}\Div^{T}
|
||||
|
||||
\mathbf{rhs} = - \Div(\MfMui)^{-1}\mathbf{M}^f_{\\frac{1}{\mu_0}}\mathbf{B}_0 + \Div\mathbf{B}_0-\diag(v)\mathbf{D} \mathbf{P}_{out}^T \mathbf{B}_{sBC}
|
||||
.. math ::
|
||||
|
||||
\mathbf{A} = \Div(\MfMui)^{-1}\Div^{T}
|
||||
|
||||
"""
|
||||
return self._Div*self.MfMuI*self._Div.T
|
||||
|
||||
|
||||
def fields(self, m):
|
||||
"""
|
||||
Return magnetic potential (u) and flux (B)
|
||||
u: defined on the cell center [nC x 1]
|
||||
B: defined on the cell center [nF x 1]
|
||||
|
||||
After we compute u, then we update B.
|
||||
|
||||
.. math ::
|
||||
|
||||
\mathbf{B}_s = (\MfMui)^{-1}\mathbf{M}^f_{\mu_0^{-1}}\mathbf{B}_0-\mathbf{B}_0 -(\MfMui)^{-1}\Div^T \mathbf{u}
|
||||
|
||||
"""
|
||||
self.makeMassMatrices(m)
|
||||
A = self.getA(m)
|
||||
rhs = self.getRHS(m)
|
||||
|
||||
m1 = sp.linalg.interface.aslinearoperator(Utils.sdiag(1/A.diagonal()))
|
||||
u, info = sp.linalg.bicgstab(A, rhs, tol=1e-6, maxiter=1000, M=m1)
|
||||
|
||||
B0 = self.getB0()
|
||||
|
||||
B = self.MfMuI*self.MfMu0*B0-B0-self.MfMuI*self._Div.T*u
|
||||
|
||||
#TODO: Create a mag fields object class.
|
||||
# F = self.getInitialFields()
|
||||
# e.g. {'B': B, 'u': u}
|
||||
|
||||
return {'B': B, 'u': u}
|
||||
|
||||
# return self.forward(m, self.getRHS, self.calcFields, F=F)
|
||||
|
||||
@Utils.timeIt
|
||||
def Jvec(self, m, v, u=None):
|
||||
def Jvec(self, m, v, fields=None):
|
||||
"""
|
||||
Computing Jacobian multiplied by vector
|
||||
|
||||
@@ -122,7 +125,7 @@ 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
|
||||
|
||||
|
||||
.. math ::
|
||||
|
||||
\\nabla_u \mathbf{C}(\mathbf{u}) = \mathbf{A}
|
||||
@@ -138,66 +141,177 @@ class MagneticsDiffSecondary(Problem.BaseProblem):
|
||||
\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[
|
||||
\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}}
|
||||
\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,
|
||||
|
||||
.. math ::
|
||||
|
||||
\\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
|
||||
|
||||
.. math ::
|
||||
|
||||
\\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
|
||||
|
||||
.. note ::
|
||||
|
||||
We only want to evalute
|
||||
|
||||
.. math ::
|
||||
|
||||
\mathbf{J}\mathbf{v} = \\frac{\delta \mathbf{P}\mathbf{B}} {\delta \mathbf{m}}\mathbf{v}
|
||||
|
||||
Since forming sensitivity matrix is very expensive in that this monster is "big" and "dense" matrix!!
|
||||
|
||||
|
||||
"""
|
||||
if u is None:
|
||||
u = self.fields(m)
|
||||
|
||||
#TODO: B, u = u['B'], u['u']
|
||||
|
||||
B, u = u['B'], u['u']
|
||||
if fields is None:
|
||||
fields = self.fields(m)
|
||||
|
||||
B, u = fields['B'], fields['u']
|
||||
mu = self.model.transform(m, asMu=True)
|
||||
dmudm = self.model.transform(m, asMu=True)
|
||||
|
||||
P = self.data.projectFieldsDeriv(u)
|
||||
|
||||
A = self.getA(m)
|
||||
dCdu = A
|
||||
# (Av_m)^-1
|
||||
# -(Av_m)^-2 * MfMu_dm * d/dm(1/mu(m))
|
||||
# (Av_m)^-2 * MfMu_dm * diag(mu(m)^-2) * mT_dm
|
||||
|
||||
#TODO: only works for diagonal MfMui
|
||||
# Some chain rule!
|
||||
|
||||
|
||||
# harm_dm = Utils.sdiag(self.MfMui.diagonal()**(-2))
|
||||
# MfMu_dm = self.mesh.getFaceMassDeriv()
|
||||
# dmuI_dm = Utils.sdiag(mu**(-2))
|
||||
# mT_dm = self.model.transformDeriv(m, asMu=True)
|
||||
|
||||
getFIPconst = 1./3
|
||||
MfMuIvec = 1/self.MfMui.diagonal()*getFIPconst
|
||||
dMfMuI = Utils.sdiag(MfMuIvec**2)*self.mesh.aveF2CC.T*Utils.sdiag(self.mesh.vol*1./mu**2)
|
||||
dmudm = self.model.transformDeriv(m, asMu=True)
|
||||
dmdmu = Utils.sdiag(1/(dmudm.diagonal()))
|
||||
|
||||
vol = self.mesh.vol
|
||||
Div = self._Div
|
||||
# lots-o-bracket for vector multiplication first!
|
||||
# MfMu_dmXv = harm_dm * ( MfMu_dm * ( dmuI_dm * ( mT_dm * v ) ) )
|
||||
#dCdm_A = D * ( Utils.sdiag( D.T * u ) * MfMu_dmXv )
|
||||
|
||||
dCdm_A = dmudm*Div * ( Utils.sdiag( Div.T * u * dMfMuI ) )
|
||||
|
||||
# rhs = D * MfMuI * MfMu0 * B0
|
||||
|
||||
Dface = self.mesh.faceDiv
|
||||
P = self.data.projectFieldsDeriv(B) # Projection matrix
|
||||
B0 = self.getB0()
|
||||
|
||||
#TODO: add congrous stuff
|
||||
dCdm_RHS = dmudm* Div * Utils.sdiag( self.MfMu0*B0 ) * dMfMuI - Utils.sdiag(self.mesh.vol)*self.mesh.faceDiv*self.Pout.T*self.Bbc*self.Bbc_const
|
||||
MfMuIvec = 1/self.MfMui.diagonal()
|
||||
dMfMuI = Utils.sdiag(MfMuIvec**2)*self.mesh.aveF2CC.T*Utils.sdiag(vol*1./mu**2)
|
||||
|
||||
# A = self._Div*self.MfMuI*self._Div.T
|
||||
# RHS = Div*MfMuI*MfMu0*B0 - Div*B0 + Mc*Dface*Pout.T*Bbc
|
||||
# C(m,u) = A*m-rhs
|
||||
# dudm = -(dCdu)^(-1)dCdm
|
||||
|
||||
# c(m,u) = A(m)u - rhs(m)
|
||||
dCdm = dCdm_A - dCdm_RHS
|
||||
dCdu = self.getA(m)
|
||||
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_v = dCdm_A*v - dCdm_RHSv
|
||||
|
||||
solve = Solver(dCdu)
|
||||
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)
|
||||
|
||||
#TODO: Multiply by the dP(u(m))/du
|
||||
# We transformed u in our fields object.
|
||||
# ( dBdu * + dBdm(u) )
|
||||
Jv = - P * solve.solve(dCdm)
|
||||
return Utils.mkvc(Jv)
|
||||
if info > 0:
|
||||
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
|
||||
dBdmv = ( Utils.sdiag(self.MfMu0*B0)*(dMfMuI * (dmudm*v)) \
|
||||
- Utils.sdiag(Div.T*u)*(dMfMuI* (dmudm*v)) \
|
||||
- self.MfMuI*(Div.T* (dudm)) )
|
||||
|
||||
return Utils.mkvc(P*dBdmv)
|
||||
|
||||
@Utils.timeIt
|
||||
def Jtvec(self, m, v, fields=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[
|
||||
\diag(\M^f_{\mu_{0}^{-1} } \mathbf{B}_0) \dMfMuI \\
|
||||
- \diag (\Div^T\mathbf{u})\dMfMuI
|
||||
\\right ]\\right]^{T}
|
||||
|
||||
- \left[\mathbf{P}_{deriv}(\MfMui)^{-1}\Div^T\\frac{\delta\mathbf{u}}{\delta \mathbf{m}} \\right]^{T}
|
||||
|
||||
where
|
||||
|
||||
.. math ::
|
||||
|
||||
\mathbf{P}_{derv} = \\frac{\partial \mathbf{P}}{\partial\mathbf{B}}
|
||||
|
||||
.. note ::
|
||||
|
||||
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)
|
||||
|
||||
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()))
|
||||
|
||||
vol = self.mesh.vol
|
||||
Div = self._Div
|
||||
Dface = self.mesh.faceDiv
|
||||
P = self.data.projectFieldsDeriv(B) # Projection matrix
|
||||
B0 = self.getB0()
|
||||
|
||||
MfMuIvec = 1/self.MfMui.diagonal()
|
||||
dMfMuI = Utils.sdiag(MfMuIvec**2)*self.mesh.aveF2CC.T*Utils.sdiag(vol*1./mu**2)
|
||||
|
||||
# A = self._Div*self.MfMuI*self._Div.T
|
||||
# RHS = Div*MfMuI*MfMu0*B0 - Div*B0 + Mc*Dface*Pout.T*Bbc
|
||||
# C(m,u) = A*m-rhs
|
||||
# dudm = -(dCdu)^(-1)dCdm
|
||||
|
||||
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)
|
||||
|
||||
if info > 0:
|
||||
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_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_RHStsol = dCdm_RHS1tsol - dCdm_RHS2tsol
|
||||
|
||||
# dCdm_RHSv = dCdm_RHS1*(dmudm*v) + dCdm_RHS2v
|
||||
# dCdm_v = dCdm_A*v - dCdm_RHSv
|
||||
|
||||
Ctv = dCdm_Atsol - dCdm_RHStsol
|
||||
|
||||
# B = self.MfMuI*self.MfMu0*B0-B0-self.MfMuI*self._Div.T*u
|
||||
# dBdm = d\mudm*dBd\mu
|
||||
# dPBdm^T*v = Atemp^T*P^T*v - Btemp^T*P^T*v - Ctv
|
||||
|
||||
Atemp = Utils.sdiag(self.MfMu0*B0)*(dMfMuI * (dmudm))
|
||||
Btemp = Utils.sdiag(Div.T*u)*(dMfMuI* (dmudm))
|
||||
Jtv = Atemp.T*(P.T*v) - Btemp.T*(P.T*v) - Ctv
|
||||
|
||||
return Utils.mkvc(Jtv)
|
||||
|
||||
|
||||
|
||||
@@ -218,7 +332,12 @@ if __name__ == '__main__':
|
||||
# mu = (1.+chi)*mu_0
|
||||
|
||||
data = BaseMag.BaseMagData()
|
||||
data.setBackgroundField(x=1., y=1., z=0.)
|
||||
Inc = 90.
|
||||
Dec = 0.
|
||||
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)
|
||||
@@ -230,30 +349,14 @@ if __name__ == '__main__':
|
||||
|
||||
prob.pair(data)
|
||||
|
||||
B = prob.fields(chi)
|
||||
# mesh.plotSlice(B, 'F', view='vec', showIt=True)
|
||||
dpred = data.dpred(chi)
|
||||
fig = plt.figure( figsize = (8,5) )
|
||||
ax = plt.subplot(111)
|
||||
dat = plt.imshow(np.reshape(dpred, (xr.size, yr.size), order='F'), extent=[min(xr), max(xr), min(yr), max(yr)])
|
||||
plt.colorbar(dat, ax = ax)
|
||||
plt.show()
|
||||
|
||||
dpred = data.dpred(chi, u=B)
|
||||
|
||||
|
||||
# ##################
|
||||
# # Test J
|
||||
# ##################
|
||||
|
||||
# d_chi = 0.8*chi #np.random.rand(mesh.nCz)
|
||||
# d_sph_ind = spheremodel(mesh, 0., 0., -100., 50)
|
||||
# d_chi[d_sph_ind] = 0.02
|
||||
|
||||
# from SimPEG.Tests import checkDerivative
|
||||
|
||||
# derChk = lambda m: [prob.data.dpred(m), lambda mx: -prob.Jvec(chi, mx)]
|
||||
# print '\n'
|
||||
# passed = checkDerivative(derChk, chi, plotIt=False, dx=d_chi, num=2)
|
||||
|
||||
|
||||
# # plt.pcolor(X, Y, dpred.reshape(X.shape, order='F'))
|
||||
|
||||
# # plt.show()
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -4,7 +4,7 @@ import matplotlib.pyplot as plt
|
||||
import simpegPF as PF
|
||||
|
||||
|
||||
class MagProblemTests(unittest.TestCase):
|
||||
class MagFwdProblemTests(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
|
||||
@@ -29,7 +29,13 @@ class MagProblemTests(unittest.TestCase):
|
||||
def test_anal_forward(self):
|
||||
|
||||
data = PF.BaseMag.BaseMagData()
|
||||
data.setBackgroundField(x=1., y=1., z=0.)
|
||||
|
||||
Inc = 90.
|
||||
Dec = 0.
|
||||
Btot = 51000
|
||||
|
||||
b0 = PF.MagAnalytics.IDTtoxyz(Inc, Dec, Btot)
|
||||
data.setBackgroundField(Inc, Dec, Btot)
|
||||
xr = np.linspace(-300, 300, 41)
|
||||
yr = np.linspace(-300, 300, 41)
|
||||
X, Y = np.meshgrid(xr, yr)
|
||||
@@ -38,9 +44,10 @@ class MagProblemTests(unittest.TestCase):
|
||||
data.rxLoc = rxLoc
|
||||
|
||||
self.prob.pair(data)
|
||||
B = self.prob.fields(self.chi)
|
||||
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,np.array([1.,1.,0.]),'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])
|
||||
@@ -50,6 +57,6 @@ class MagProblemTests(unittest.TestCase):
|
||||
else:
|
||||
print "Anaytic test is passed"
|
||||
pass
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
|
||||
@@ -0,0 +1,316 @@
|
||||
import unittest
|
||||
from SimPEG import *
|
||||
from simpegPF import BaseMag
|
||||
import matplotlib.pyplot as plt
|
||||
import simpegPF as PF
|
||||
from scipy.constants import mu_0
|
||||
|
||||
|
||||
class MagSensProblemTests(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
|
||||
hxind = ((5,25,1.3),(21, 25.),(5,25,1.3))
|
||||
hyind = ((5,25,1.3),(21, 25.),(5,25,1.3))
|
||||
hzind = ((5,25,1.3),(20, 25.),(5,25,1.3))
|
||||
hx, hy, hz = Utils.meshTensors(hxind, hyind, hzind)
|
||||
M = Mesh.TensorMesh([hx, hy, hz], [-hx.sum()/2,-hy.sum()/2,-hz.sum()/2])
|
||||
chibkg = 0.001
|
||||
chiblk = 0.01
|
||||
chi = np.ones(M.nC)*chibkg
|
||||
|
||||
Inc = 90.
|
||||
Dec = 0.
|
||||
Btot = 51000
|
||||
|
||||
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)
|
||||
|
||||
data = BaseMag.BaseMagData()
|
||||
|
||||
data.setBackgroundField(Inc, Dec, Btot)
|
||||
xr = np.linspace(-300, 300, 41)
|
||||
yr = np.linspace(-300, 300, 41)
|
||||
X, Y = np.meshgrid(xr, yr)
|
||||
Z = np.ones((xr.size, yr.size))*150
|
||||
rxLoc = np.c_[Utils.mkvc(X), Utils.mkvc(Y), Utils.mkvc(Z)]
|
||||
data.rxLoc = rxLoc
|
||||
|
||||
|
||||
|
||||
prob = PF.Magnetics.MagneticsDiffSecondary(M, model)
|
||||
prob.pair(data)
|
||||
dpre = data.dpred(chi)
|
||||
|
||||
fields = prob.fields(chi)
|
||||
self.u = fields['u']
|
||||
self.B = fields['B']
|
||||
|
||||
self.data = data
|
||||
self.model = model
|
||||
self.prob = prob
|
||||
self.M = M
|
||||
self.chi = chi
|
||||
|
||||
|
||||
|
||||
def test_mass(self):
|
||||
print '\n >>Derivative for MfMuI works.'
|
||||
mu = self.model.transform(self.chi, asMu=True)
|
||||
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()
|
||||
|
||||
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)
|
||||
|
||||
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
|
||||
vol = self.prob.mesh.vol
|
||||
aveF2CC = self.prob.mesh.aveF2CC
|
||||
|
||||
def Cm_A(chi):
|
||||
dmudm = self.model.transformDeriv(chi, asMu=True)
|
||||
u = self.u
|
||||
# chi = mu/mu_0-1
|
||||
self.prob.makeMassMatrices(chi)
|
||||
mu = self.model.transform(self.chi, asMu=True)
|
||||
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
|
||||
|
||||
return Cm_A
|
||||
|
||||
def dCdm_A(chi, v):
|
||||
|
||||
dmudm = self.model.transformDeriv(chi, asMu=True)
|
||||
u = self.u
|
||||
self.prob.makeMassMatrices(chi)
|
||||
mu = self.model.transform(self.chi, asMu=True)
|
||||
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 )
|
||||
|
||||
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)
|
||||
|
||||
|
||||
def test_dCdmu_RHS(self):
|
||||
print '\n >>Derivative for Cm_RHS.'
|
||||
u = self.u
|
||||
Div = self.prob._Div
|
||||
mu = self.model.transform(self.chi, asMu=True)
|
||||
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):
|
||||
|
||||
self.prob.makeMassMatrices(chi)
|
||||
dmudm = self.model.transformDeriv(chi, asMu=True)
|
||||
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
|
||||
|
||||
return RHS
|
||||
|
||||
|
||||
def dCdm_RHS(chi, v):
|
||||
|
||||
|
||||
self.prob.makeMassMatrices(chi)
|
||||
dmudm = self.model.transformDeriv(chi, asMu=True)
|
||||
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)
|
||||
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)
|
||||
|
||||
|
||||
def test_dudm(self):
|
||||
print ">> Derivative test for dudm"
|
||||
u = self.u
|
||||
Div = self.prob._Div
|
||||
mu = self.model.transform(self.chi, asMu=True)
|
||||
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 dudm(chi, v):
|
||||
|
||||
chi = mu/mu_0-1
|
||||
self.prob.makeMassMatrices(chi)
|
||||
u = self.u
|
||||
dmudm = self.model.transformDeriv(chi, asMu=True)
|
||||
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)
|
||||
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
|
||||
|
||||
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
|
||||
|
||||
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, asMu=True)
|
||||
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 dBdm(chi, v):
|
||||
|
||||
chi = mu/mu_0-1
|
||||
self.prob.makeMassMatrices(chi)
|
||||
u = self.u
|
||||
dmudm = self.model.transformDeriv(chi, asMu=True)
|
||||
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)
|
||||
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)) )
|
||||
|
||||
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
|
||||
|
||||
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)
|
||||
self.assertTrue(passed)
|
||||
|
||||
|
||||
|
||||
def test_Jvec(self):
|
||||
print ">> Derivative test for Jvec"
|
||||
mu = self.model.transform(self.chi, asMu=True)
|
||||
|
||||
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.data.dpred(m), lambda mx: self.prob.Jvec(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_Jtvec(self):
|
||||
print ">> Derivative test for Jtvec"
|
||||
mu = self.model.transform(self.chi, asMu=True)
|
||||
dobs = self.data.dpred(self.chi)
|
||||
|
||||
def misfit (m, dobs):
|
||||
dpre = self.data.dpred(m)
|
||||
misfit = 0.5*np.linalg.norm(dpre-dobs)**2
|
||||
residual = dpre-dobs
|
||||
dmisfit = self.prob.Jtvec(self.chi, residual)
|
||||
|
||||
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)
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
File diff suppressed because one or more lines are too long
Reference in New Issue
Block a user