wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-07-02 12:32:09 +08:00
Code Issues Packages Projects Releases Wiki Activity

working to debug JTv for e-formulation... still some work to do

Browse Source
This commit is contained in:
Lindsey Heagy
2016-03-20 11:42:31 -07:00
parent 0be942730a
commit 2d8bbdce45
4 changed files with 116 additions and 72 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/EM/FDEM/FDEM.py
+36 -36
View File
@@ -18,9 +18,9 @@ class BaseFDEMProblem(BaseEMProblem):
{\mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{b} - \mathbf{M_{\sigma}^e} \mathbf{e} = \mathbf{s_e}}
if using the E-B formulation (:code:`Problem_e`
or :code:`Problem_b`). Note that in this case, :math:`\mathbf{s_e}` is an integrated quantity.
or :code:`Problem_b`). Note that in this case, :math:`\mathbf{s_e}` is an integrated quantity.
If we write Maxwell's equations in terms of
If we write Maxwell's equations in terms of
\\\(\\\mathbf{h}\\\) and current density \\\(\\\mathbf{j}\\\)
.. math ::
@@ -28,7 +28,7 @@ class BaseFDEMProblem(BaseEMProblem):
\mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{j} + i \omega \mathbf{M_{\mu}^e} \mathbf{h} = \mathbf{s_m} \\\\
\mathbf{C} \mathbf{h} - \mathbf{j} = \mathbf{s_e}
if using the H-J formulation (:code:`Problem_j` or :code:`Problem_h`). Note that here, :math:`\mathbf{s_m}` is an integrated quantity.
if using the H-J formulation (:code:`Problem_j` or :code:`Problem_h`). Note that here, :math:`\mathbf{s_m}` is an integrated quantity.
The problem performs the elimination so that we are solving the system for \\\(\\\mathbf{e},\\\mathbf{b},\\\mathbf{j} \\\) or \\\(\\\mathbf{h}\\\)
"""
@@ -39,7 +39,7 @@ class BaseFDEMProblem(BaseEMProblem):
def fields(self, m=None):
"""
Solve the forward problem for the fields.
:param numpy.array m: inversion model (nP,)
:rtype numpy.array:
:return F: forward solution
@@ -65,9 +65,9 @@ class BaseFDEMProblem(BaseEMProblem):
:param numpy.array m: inversion model (nP,)
:param numpy.array v: vector which we take sensitivity product with (nP,)
:param SimPEG.EM.FDEM.Fields u: fields object
:param SimPEG.EM.FDEM.Fields u: fields object
:rtype numpy.array:
:return: Jv (ndata,)
:return: Jv (ndata,)
"""
if u is None:
@@ -85,13 +85,13 @@ class BaseFDEMProblem(BaseEMProblem):
ftype = self._fieldType + 'Solution'
u_src = u[src, ftype]
dA_dm_v = self.getADeriv(freq, u_src, v)
dRHS_dm_v = self.getRHSDeriv(freq, src, v)
dRHS_dm_v = self.getRHSDeriv(freq, src, v)
du_dm_v = Ainv * ( - dA_dm_v + dRHS_dm_v )
for rx in src.rxList:
df_dmFun = getattr(u, '_%sDeriv'%rx.projField, None)
df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False)
df_dm_v = np.array(df_dm_v,dtype=complex)
df_dm_v = np.array(df_dm_v, dtype=complex)
Jv[src, rx] = rx.evalDeriv(src, self.mesh, u, df_dm_v)
Ainv.clean()
return Utils.mkvc(Jv)
@@ -102,9 +102,9 @@ class BaseFDEMProblem(BaseEMProblem):
:param numpy.array m: inversion model (nP,)
:param numpy.array v: vector which we take adjoint product with (nP,)
:param SimPEG.EM.FDEM.Fields u: fields object
:param SimPEG.EM.FDEM.Fields u: fields object
:rtype numpy.array:
:return: Jv (ndata,)
:return: Jv (ndata,)
"""
if u is None:
@@ -148,7 +148,7 @@ class BaseFDEMProblem(BaseEMProblem):
Jtv += - np.array(df_dmT,dtype=complex).real
else:
raise Exception('Must be real or imag')
ATinv.clean()
return Utils.mkvc(Jtv)
@@ -158,7 +158,7 @@ class BaseFDEMProblem(BaseEMProblem):
Evaluates the sources for a given frequency and puts them in matrix form
:param float freq: Frequency
:rtype: (numpy.ndarray, numpy.ndarray)
:rtype: (numpy.ndarray, numpy.ndarray)
:return: S_m, S_e (nE or nF, nSrc)
"""
Srcs = self.survey.getSrcByFreq(freq)
@@ -211,7 +211,7 @@ class Problem_e(BaseFDEMProblem):
def getA(self, freq):
"""
System matrix
.. math ::
\mathbf{A} = \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{C} + i \omega \mathbf{M^e_{\sigma}}
@@ -234,12 +234,12 @@ class Problem_e(BaseFDEMProblem):
.. math ::
\\frac{\mathbf{A}(\mathbf{m}) \mathbf{v}}{d \mathbf{m}} = i \omega \\frac{d \mathbf{M^e_{\sigma}}\mathbf{v} }{d\mathbf{m}}
:param float freq: frequency
:param numpy.ndarray u: solution vector (nE,)
:param float freq: frequency
:param numpy.ndarray u: solution vector (nE,)
:param numpy.ndarray v: vector to take prodct with (nP,) or (nD,) for adjoint
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
:return: derivative of the system matrix times a vector (nP,) or adjoint (nD,)
:return: derivative of the system matrix times a vector (nP,) or adjoint (nD,)
"""
dsig_dm = self.curModel.sigmaDeriv
@@ -252,13 +252,13 @@ class Problem_e(BaseFDEMProblem):
def getRHS(self, freq):
"""
Right hand side for the system
Right hand side for the system
.. math ::
\mathbf{RHS} = \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f}\mathbf{s_m} -i\omega\mathbf{M_e}\mathbf{s_e}
:param float freq: Frequency
:rtype: numpy.ndarray
:rtype: numpy.ndarray
:return: RHS (nE, nSrc)
"""
@@ -270,7 +270,7 @@ class Problem_e(BaseFDEMProblem):
def getRHSDeriv(self, freq, src, v, adjoint=False):
"""
Derivative of the right hand side with respect to the model
Derivative of the right hand side with respect to the model
:param float freq: frequency
:param SimPEG.EM.FDEM.Src src: FDEM source
@@ -350,12 +350,12 @@ class Problem_b(BaseFDEMProblem):
.. math ::
\\frac{\mathbf{A}(\mathbf{m}) \mathbf{v}}{d \mathbf{m}} = \mathbf{C} \\frac{\mathbf{M^e_{\sigma}} \mathbf{v}}{d\mathbf{m}}
:param float freq: frequency
:param numpy.ndarray u: solution vector (nF,)
:param float freq: frequency
:param numpy.ndarray u: solution vector (nF,)
:param numpy.ndarray v: vector to take prodct with (nP,) or (nD,) for adjoint
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
:return: derivative of the system matrix times a vector (nP,) or adjoint (nD,)
:return: derivative of the system matrix times a vector (nP,) or adjoint (nD,)
"""
MfMui = self.MfMui
@@ -377,13 +377,13 @@ class Problem_b(BaseFDEMProblem):
def getRHS(self, freq):
"""
Right hand side for the system
Right hand side for the system
.. math ::
\mathbf{RHS} = \mathbf{s_m} + \mathbf{M^e_{\sigma}}^{-1}\mathbf{s_e}
:param float freq: Frequency
:rtype: numpy.ndarray
:rtype: numpy.ndarray
:return: RHS (nE, nSrc)
"""
@@ -501,12 +501,12 @@ class Problem_j(BaseFDEMProblem):
\\frac{\mathbf{A(\sigma)} \mathbf{v}}{d \mathbf{m}} = \mathbf{C} \mathbf{M^e_{mu^{-1}}} \mathbf{C^{\\top}} \\frac{d \mathbf{M^f_{\sigma^{-1}}}\mathbf{v} }{d \mathbf{m}}
:param float freq: frequency
:param numpy.ndarray u: solution vector (nF,)
:param float freq: frequency
:param numpy.ndarray u: solution vector (nF,)
:param numpy.ndarray v: vector to take prodct with (nP,) or (nD,) for adjoint
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
:return: derivative of the system matrix times a vector (nP,) or adjoint (nD,)
:return: derivative of the system matrix times a vector (nP,) or adjoint (nD,)
"""
MeMuI = self.MeMuI
@@ -526,7 +526,7 @@ class Problem_j(BaseFDEMProblem):
def getRHS(self, freq):
"""
Right hand side for the system
Right hand side for the system
.. math ::
@@ -550,7 +550,7 @@ class Problem_j(BaseFDEMProblem):
def getRHSDeriv(self, freq, src, v, adjoint=False):
"""
Derivative of the right hand side with respect to the model
Derivative of the right hand side with respect to the model
:param float freq: frequency
:param SimPEG.EM.FDEM.Src src: FDEM source
@@ -630,12 +630,12 @@ class Problem_h(BaseFDEMProblem):
.. math::
\\frac{\mathbf{A}(\mathbf{m}) \mathbf{v}}{d \mathbf{m}} = \mathbf{C}^{\\top}\\frac{d \mathbf{M^f_{\\rho}}\mathbf{v} }{d\mathbf{m}}
:param float freq: frequency
:param numpy.ndarray u: solution vector (nE,)
:param float freq: frequency
:param numpy.ndarray u: solution vector (nE,)
:param numpy.ndarray v: vector to take prodct with (nP,) or (nD,) for adjoint
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
:return: derivative of the system matrix times a vector (nP,) or adjoint (nD,)
:return: derivative of the system matrix times a vector (nP,) or adjoint (nD,)
"""
MeMu = self.MeMu
@@ -648,14 +648,14 @@ class Problem_h(BaseFDEMProblem):
def getRHS(self, freq):
"""
Right hand side for the system
Right hand side for the system
.. math ::
\mathbf{RHS} = \mathbf{M^e} \mathbf{s_m} + \mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{s_e}
:param float freq: Frequency
:rtype: numpy.ndarray
:rtype: numpy.ndarray
:return: RHS (nE, nSrc)
"""
@@ -667,7 +667,7 @@ class Problem_h(BaseFDEMProblem):
def getRHSDeriv(self, freq, src, v, adjoint=False):
"""
Derivative of the right hand side with respect to the model
Derivative of the right hand side with respect to the model
:param float freq: frequency
:param SimPEG.EM.FDEM.Src src: FDEM source
SimPEG/EM/TDEM/SurveyTDEM.py
+3 -2
View File
@@ -132,7 +132,6 @@ class BaseSrc(SimPEG.Survey.BaseSrc):
@waveform.setter
def waveform(self, val):
if self.waveform is None:
print val
val._assertMatchesPair(self.waveformPair)
self._mapping = val
else:
@@ -252,8 +251,10 @@ class MagDipole(BaseSrc):
C = prob.mesh.edgeCurl
S_e = self.S_e(prob, prob.t0)
# S_e doesn't depend on the model
if adjoint:
raise NotImplementedError
return MeSigmaIDeriv( -S_e + C.T * ( MfMui * b ) ).T * v
return MeSigmaIDeriv( -S_e + C.T * ( MfMui * b ) ) * v
SimPEG/EM/TDEM/TDEM.py
+43 -21
View File
@@ -85,8 +85,16 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
# mat to store previous time-step's solution deriv times a vector for each source
# size: nu x nSrc
dun_dm_v = self.getInitialFieldsDeriv(v) # can over-write this at each timestep
# this is a bit silly
# if self._fieldType is 'b' or self._fieldType is 'j':
# ifields = np.zeros((self.mesh.nF, len(Srcs)))
# elif self._fieldType is 'e' or self._fieldType is 'h':
# ifields = np.zeros((self.mesh.nE, len(Srcs)))
# for i, src in enumerate(self.survey.srcList):
dun_dm_v = np.hstack([Utils.mkvc(self.getInitialFieldsDeriv(src,v),2) for src in self.survey.srcList]) # can over-write this at each timestep
#
df_dm_v = Fields_Derivs(self.mesh, self.survey) # store the field derivs we need to project to calc full deriv
@@ -117,10 +125,7 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
dA_dm_v = self.getAdiagDeriv(tInd, un_src, v) # cell centered on time mesh
dRHS_dm_v = self.getRHSDeriv(tInd+1, src, v) # on nodes of time mesh
# if tInd >= 1:
dAsubdiag_dm_v = self.getAsubdiagDeriv(tInd, u[src,ftype,tInd], v)
# else:
# dAsubdiag_dm_v = Zero()
JRHS = dRHS_dm_v - dAsubdiag_dm_v - dA_dm_v
@@ -201,6 +206,8 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
Adiag = self.getAdiag(tInd)
AdiagTinv = self.Solver(Adiag.T, **self.solverOpts)
dAsubdiag_dm_v = Zero()
if tInd < self.nT - 1:
Asubdiag = self.getAsubdiag(tInd+1)
@@ -211,20 +218,35 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
# last timestep (first to be solved)
ATinv_df_duT_v[isrc,:] = AdiagTinv * df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1]
elif tInd > -1:
# else:
ATinv_df_duT_v[isrc,:] = AdiagTinv * (Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1]) - Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:]))
else:
# AdiagTinv = I
ATinv_df_duT_v[isrc,:] = (Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1]) - Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:]))
ATinv_df_duT_v[isrc,:] = Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1]) - Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:])
# - Utils.mkvc(Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:]))
# (Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1]) - Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:]))
if tInd < self.nT - 1:
dAsubdiagT_dm_v = self.getAsubdiagDeriv(tInd+1, u[src,ftype,tInd+1], ATinv_df_duT_v[isrc,:], adjoint = True)
if tInd > -1:
un_src = u[src,ftype,tInd+1]
dAT_dm_v = self.getAdiagDeriv(tInd, un_src, ATinv_df_duT_v[isrc,:], adjoint=True) # cell centered on time mesh
dRHST_dm_v = self.getRHSDeriv(tInd+1, src, ATinv_df_duT_v[isrc,:], adjoint=True) # on nodes of time mesh
JTv = JTv + Utils.mkvc(-dAT_dm_v + dRHST_dm_v)
# else:
# print 'here'
# un_src = u[src,ftype,tInd+1]
JTv = JTv + Utils.mkvc(- dAT_dm_v - dAsubdiag_dm_v + dRHST_dm_v)
else:
# dA_dm_v = self.getInitialFieldsDeriv(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1], adjoint=True)
# print np.linalg.norm(self.getInitialFieldsDeriv(src, df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1], adjoint=True))
# print np.linalg.norm(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1])
# vec = - Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:]) + Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1])
# dAsubdiagT_dm_v = self.getAsubdiagDeriv(tInd+1, u[src,ftype,tInd+1], Utils.mkvc(ATinv_df_duT_v[isrc,:]), adjoint = True)
dRHST_dm_v = Utils.mkvc(self.getInitialFieldsDeriv(src, Utils.mkvc(ATinv_df_duT_v[isrc,:]) , adjoint=True))
JTv = JTv + Utils.mkvc( -dAsubdiagT_dm_v + dRHST_dm_v) #
# # dAT_dm_v = self.getAdiagDeriv(tInd, un_src, ATinv_df_duT_v[isrc,:], adjoint=True) # cell centered on time mesh
# dRHST_dm_v0 = self.getRHSDeriv(tInd+1, src, ATinv_df_duT_v[isrc,:], adjoint=True) # on nodes of time mesh
# dRHST_dm_v1 = self.getInitialFieldsDeriv( Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1]), adjoint=True)
@@ -286,21 +308,20 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
return ifields
def getInitialFieldsDeriv(self, v, adjoint=False):
def getInitialFieldsDeriv(self, src, v, adjoint=False):
Srcs = self.survey.srcList
if adjoint is False:
if self._fieldType is 'b' or self._fieldType is 'j':
ifieldsDeriv = np.zeros(self.mesh.nF)
elif self._fieldType is 'e' or self._fieldType is 'h':
ifieldsDeriv = np.zeros(self.mesh.nE)
if self._fieldType is 'b' or self._fieldType is 'j':
ifieldsDeriv = np.zeros((self.mesh.nF, len(Srcs)))
elif self._fieldType is 'e' or self._fieldType is 'h':
ifieldsDeriv = np.zeros((self.mesh.nE, len(Srcs)))
elif adjoint is True:
ifieldsDeriv = np.zeros(self.mapping.nP)
for i,src in enumerate(Srcs):
ifieldsDeriv[:,i] = ifieldsDeriv[:,i] + getattr(src, '%sInitialDeriv'%self._fieldType, None)(self,v,adjoint)
ifieldsDeriv = Utils.mkvc(getattr(src, '%sInitialDeriv'%self._fieldType, None)(self,v,adjoint)) + ifieldsDeriv
if adjoint is True:
ifieldsDeriv = np.zeros_like(self.curModel)
print v.shape
# ifieldsDeriv = Utils.mkvc(getattr(src, '%sInitialDeriv'%self._fieldType, None)(self,v,adjoint)) + ifieldsDeriv
# ifieldsDeriv = self.getAdiagDeriv(None, u, v, adjoint)
# ifieldsDeriv = ifieldsDeriv.sum()
@@ -485,6 +506,7 @@ class Problem_e(BaseTDEMProblem):
def getAdiagDeriv(self, tInd, u, v, adjoint=False):
assert tInd >= 0 and tInd < self.nT
dt = self.timeSteps[tInd]
C = self.mesh.edgeCurl
@@ -492,7 +514,7 @@ class Problem_e(BaseTDEMProblem):
MeSigmaDeriv = self.MeSigmaDeriv(u)
if adjoint:
raise 1./dt * MeSigmaDeriv.T * v
return 1./dt * MeSigmaDeriv.T * v
return 1./dt * MeSigmaDeriv * v
tests/em/tdem/test_TDEM_b_DerivAdjoint.py
+34 -13
View File
@@ -5,7 +5,7 @@ from SimPEG import EM
plotIt = False
testDeriv = True
testAdjoint = False
testAdjoint = True
TOL = 1e-5
@@ -54,13 +54,16 @@ class TDEM_DerivTests(unittest.TestCase):
# ====== TEST A ========== #
def test_AderivTest(self):
prb, m0, mesh = setUp()
def AderivTest(self, prbtype):
prb, m0, mesh = setUp(prbtype)
tInd = 2
if prbtype == 'b':
nu = mesh.nF
elif prbtype == 'e':
nu = mesh.nE
v = np.random.rand(nu)
v = np.random.rand(mesh.nF)
def AderivTest(m):
def AderivFun(m):
prb.curModel = m
A = prb.getAdiag(tInd)
Av = A*v
@@ -69,26 +72,41 @@ class TDEM_DerivTests(unittest.TestCase):
return Av, ADeriv_dm
print '\n Testing ADeriv'
Tests.checkDerivative(AderivTest, m0, plotIt=False, num=4, eps=1e-20)
print '\n Testing ADeriv %s'%(prbtype)
Tests.checkDerivative(AderivFun, m0, plotIt=False, num=4, eps=1e-20)
def test_A_adjointTest(self):
prb, m0, mesh = setUp()
def A_adjointTest(self,prbtype):
prb, m0, mesh = setUp(prbtype)
tInd = 2
print '\n Testing A_adjoint'
m = np.random.rand(prb.mapping.nP)
v = np.random.rand(prb.mesh.nF)
u = np.random.rand(prb.mesh.nF)
if prbtype == 'b':
nu = prb.mesh.nF
elif prbtype == 'e':
nu = prb.mesh.nE
v = np.random.rand(nu)
u = np.random.rand(nu)
prb.curModel = m0
tInd = 2 # not actually used
V1 = v.dot(prb.getAdiagDeriv(tInd, u, m))
V2 = m.dot(prb.getAdiagDeriv(tInd, u, v, adjoint=True))
passed = np.abs(V1-V2) < TOL * (np.abs(V1) + np.abs(V2))/2.
print 'AdjointTest', V1, V2, passed
print 'AdjointTest %s'%(prbtype), V1, V2, passed
self.assertTrue(passed)
def test_Aderiv_b(self):
self.AderivTest('b')
def test_Aderiv_e(self):
self.AderivTest('e')
def test_Aadjoint_b(self):
self.A_adjointTest('b')
def test_Aadjoint_e(self):
self.A_adjointTest('e')
# ====== TEST Fields Deriv Pieces ========== #
def test_eDeriv_m_adjoint(self):
@@ -194,6 +212,9 @@ class TDEM_DerivTests(unittest.TestCase):
def test_Jvec_adjoint_b_ey(self):
self.JvecVsJtvecTest('b', 'ey')
def test_Jvec_adjoint_e_ey(self):
self.JvecVsJtvecTest('e', 'ey')
if __name__ == '__main__':