mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-08 05:22:40 +08:00
bx, bz running and passing, ey failing adjoint --> I think the initial fields are not being taken care of correctly in the deriv
This commit is contained in:
Lindsey Heagy
@@ -7,7 +7,7 @@ from SimPEG.Utils import Zero, Identity
|
||||
|
||||
class Fields(SimPEG.Problem.TimeFields):
|
||||
"""
|
||||
|
||||
|
||||
Fancy Field Storage for a TDEM survey. Only one field type is stored for
|
||||
each problem, the rest are computed. The fields obejct acts like an array and is indexed by
|
||||
|
||||
@@ -29,7 +29,7 @@ class Fields(SimPEG.Problem.TimeFields):
|
||||
"""
|
||||
|
||||
knownFields = {}
|
||||
dtype = float
|
||||
dtype = float
|
||||
|
||||
|
||||
|
||||
@@ -40,7 +40,7 @@ class Fields_Derivs(Fields):
|
||||
'hDeriv': 'E',
|
||||
'jDeriv': 'F'
|
||||
}
|
||||
|
||||
|
||||
|
||||
class Fields_b(Fields):
|
||||
"""Fancy Field Storage for a TDEM survey."""
|
||||
@@ -65,7 +65,7 @@ class Fields_b(Fields):
|
||||
def _bDeriv_m(self, tInd, src, v, adjoint=False):
|
||||
return Zero()
|
||||
|
||||
def _bDeriv(self, tInd, src, dun_dm_v, v, adjoint=False):
|
||||
def _bDeriv(self, tInd, src, dun_dm_v, v, adjoint=False):
|
||||
if adjoint is True:
|
||||
return self._bDeriv_u(tInd, src, v, adjoint), self._bDeriv_m(tInd, src, v, adjoint)
|
||||
return self._bDeriv_u(tInd, src, dun_dm_v) + self._bDeriv_m(tInd, src, v)
|
||||
@@ -73,27 +73,30 @@ class Fields_b(Fields):
|
||||
def _e(self, bSolution, srcList, tInd):
|
||||
e = self.MeSigmaI * ( self.edgeCurl.T * ( self.MfMui * bSolution ) )
|
||||
for i, src in enumerate(srcList):
|
||||
_, S_e = src.eval(self.survey.prob, self.survey.prob.times[tInd])
|
||||
_, S_e = src.eval(self.survey.prob, self.survey.prob.times[tInd])
|
||||
e[:,i] = e[:,i] - self.MeSigmaI * S_e
|
||||
return e
|
||||
return e
|
||||
|
||||
def _eDeriv_u(self, tInd, src, dun_dm_v, adjoint = False):
|
||||
if adjoint is True:
|
||||
raise NotImplementedError
|
||||
return self.MfMui.T * ( self.edgeCurl * ( self.MeSigmaI.T * dun_dm_v ) )
|
||||
return self.MeSigmaI * ( self.edgeCurl.T * ( self.MfMui * dun_dm_v ) )
|
||||
|
||||
def _eDeriv_m(self, tInd, src, v, adjoint = False):
|
||||
if adjoint is True:
|
||||
raise NotImplementedError
|
||||
|
||||
_, S_e = src.eval(self.survey.prob, self.survey.prob.times[tInd])
|
||||
bSolution = self[[src],'bSolution',tInd]
|
||||
|
||||
_, S_e = src.eval(self.survey.prob, self.survey.prob.times[tInd])
|
||||
_, S_eDeriv = src.evalDeriv(self.survey.prob.times[tInd], self, v=v)
|
||||
if adjoint:
|
||||
v = self.MeSigmaI.T * v
|
||||
|
||||
return self.MeSigmaIDeriv(self.edgeCurl.T * ( self.MfMui * bSolution) ) * v - self.MeSigmaIDeriv(S_e) * v - self.MeSigmaI * S_eDeriv
|
||||
_, S_eDeriv_v = src.evalDeriv(self.survey.prob.times[tInd], self, v=v, adjoint=adjoint)
|
||||
|
||||
def _eDeriv(self, tInd, src, dun_dm_v, v, adjoint=False):
|
||||
if adjoint is True:
|
||||
raise NotImplementedError
|
||||
return self.MeSigmaIDeriv(-S_e + self.edgeCurl.T * ( self.MfMui * bSolution ) ).T * v - S_eDeriv_v
|
||||
|
||||
return self.MeSigmaIDeriv(-S_e + self.edgeCurl.T * ( self.MfMui * bSolution)) * v - self.MeSigmaI * S_eDeriv_v
|
||||
|
||||
def _eDeriv(self, tInd, src, dun_dm_v, v, adjoint=False):
|
||||
if adjoint is True:
|
||||
return self._eDeriv_u(tInd, src, v, adjoint), self._eDeriv_m(tInd, src, v, adjoint)
|
||||
return self._eDeriv_u(tInd, src, dun_dm_v) + self._eDeriv_m(tInd, src, v)
|
||||
|
||||
+16
-8
@@ -163,20 +163,25 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
|
||||
|
||||
# Loop over sources and receivers to create a fields object: PT_v, df_duT_v, df_dmT_v
|
||||
for src in self.survey.srcList:
|
||||
|
||||
PT_v = Fields_Derivs(self.mesh, self.survey) # initialize storage for PT_v (don't need to preserve over sources)
|
||||
# initialize size
|
||||
df_duT_v[src, '%sDeriv'%self._fieldType, :]= np.zeros_like(u[src, self._fieldType, :])
|
||||
df_duT_v[src, '%sDeriv'%self._fieldType, :] = np.zeros_like(u[src, self._fieldType, :])
|
||||
|
||||
for rx in src.rxList:
|
||||
curPT_v = rx.evalDeriv(src, self.mesh, self.timeMesh, Utils.mkvc(v[src,rx]), adjoint=True)
|
||||
PT_v = np.reshape(curPT_v,(len(curPT_v)/self.timeMesh.nN, self.timeMesh.nN), order='F')
|
||||
PT_v[src,'%sDeriv'%rx.projField,:] = rx.evalDeriv(src, self.mesh, self.timeMesh, Utils.mkvc(v[src,rx]), adjoint=True)
|
||||
# PT_v = np.reshape(curPT_v,(len(curPT_v)/self.timeMesh.nN, self.timeMesh.nN), order='F')
|
||||
|
||||
df_duTFun = getattr(u, '_%sDeriv'%rx.projField, None)
|
||||
df_duT_v_cur, df_dmT_v = df_duTFun(None, src, None, PT_v, adjoint=True)
|
||||
|
||||
JTv = JTv + df_dmT_v
|
||||
for tInd in range(self.nT+1):
|
||||
cur = df_duTFun(tInd, src, None, Utils.mkvc(PT_v[src,'%sDeriv'%rx.projField,tInd]), adjoint=True)
|
||||
|
||||
df_duT_v[src, '%sDeriv'%self._fieldType, :] += df_duT_v_cur
|
||||
df_duT_v[src, '%sDeriv'%self._fieldType, tInd] = df_duT_v[src, '%sDeriv'%self._fieldType, tInd] + Utils.mkvc(cur[0],2)
|
||||
JTv = JTv + cur[1]
|
||||
|
||||
|
||||
|
||||
del PT_v # no longer need this
|
||||
|
||||
|
||||
AdiagTinv = None
|
||||
@@ -212,9 +217,12 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
|
||||
JTv = JTv + Utils.mkvc(-dAT_dm_v + dRHST_dm_v)
|
||||
|
||||
|
||||
# Missing the 0 step
|
||||
|
||||
# adding du_dm^T * dF_du^T * P^T vfor time 0 (no dRHS_dm_v at time 0)
|
||||
for src in self.survey.srcList:
|
||||
JTv = JTv + self.getInitialFieldsDeriv(Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,0]), adjoint=True)
|
||||
for projField in set(rx.projField):
|
||||
JTv = JTv + self.getInitialFieldsDeriv(Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,0]), adjoint=True)
|
||||
|
||||
|
||||
return Utils.mkvc(JTv)
|
||||
|
||||
@@ -4,15 +4,15 @@ from SimPEG import EM
|
||||
|
||||
plotIt = False
|
||||
|
||||
testDeriv = False
|
||||
testDeriv = True
|
||||
testAdjoint = True
|
||||
|
||||
tol = 1e-6
|
||||
TOL = 1e-5
|
||||
|
||||
def setUp(rxcomp='bz'):
|
||||
cs = 5.
|
||||
ncx = 20
|
||||
ncy = 10
|
||||
ncy = 15
|
||||
npad = 20
|
||||
hx = [(cs,ncx), (cs,npad,1.3)]
|
||||
hy = [(cs,npad,-1.3), (cs,ncy), (cs,npad,1.3)]
|
||||
@@ -74,7 +74,7 @@ class TDEM_DerivTests(unittest.TestCase):
|
||||
tInd = 0 # 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)/np.abs(V1) < tol
|
||||
passed = np.abs(V1-V2) < TOL * (np.abs(V1) + np.abs(V2))/2.
|
||||
print 'AdjointTest', V1, V2, passed
|
||||
self.assertTrue(passed)
|
||||
|
||||
@@ -107,16 +107,18 @@ class TDEM_DerivTests(unittest.TestCase):
|
||||
# ====== TEST Jtvec ========== #
|
||||
|
||||
def adjointJvecVsJtvecTest(self, rxcomp='bz'):
|
||||
print '\n Adjoint Testing Jvec, Jtvec %s' %(rxcomp)
|
||||
prb, m0, mesh = setUp(rxcomp)
|
||||
|
||||
print '\nAdjoint Testing Jvec, Jtvec %s' %(rxcomp)
|
||||
|
||||
prb, m0, mesh = setUp(rxcomp)
|
||||
m = np.random.rand(prb.mapping.nP)
|
||||
d = np.random.randn(prb.survey.nD)
|
||||
|
||||
V1 = d.dot(prb.Jvec(m0, m))
|
||||
V2 = m.dot(prb.Jtvec(m0, d))
|
||||
passed = np.abs(V1-V2)/np.abs(V1) < tol
|
||||
print 'AdjointTest', V1, V2, passed
|
||||
tol = TOL * (np.abs(V1) + np.abs(V2)) / 2.
|
||||
passed = np.abs(V1-V2) < tol
|
||||
|
||||
print ' ', V1, V2, np.abs(V1-V2), tol, passed
|
||||
self.assertTrue(passed)
|
||||
|
||||
if testAdjoint:
|
||||
|
||||
@@ -6,6 +6,7 @@ plotIt = False
|
||||
testDeriv = True
|
||||
testAdjoint = True
|
||||
|
||||
TOL = 1e-6
|
||||
|
||||
class TDEM_bDerivTests(unittest.TestCase):
|
||||
|
||||
@@ -70,8 +71,11 @@ class TDEM_bDerivTests(unittest.TestCase):
|
||||
|
||||
V1 = d.dot(prb.Jvec(m0, m))
|
||||
V2 = m.dot(prb.Jtvec(m0, d))
|
||||
print 'AdjointTest', V1, V2
|
||||
self.assertTrue(np.abs(V1-V2)/np.abs(V1) < 1e-6)
|
||||
|
||||
tol = TOL * (np.abs(V1) + np.abs(V2)) / 2.
|
||||
passed = np.abs(V1-V2) < tol
|
||||
print ' ', V1, V2, np.abs(V1-V2), tol, passed
|
||||
self.assertTrue(passed)
|
||||
|
||||
|
||||
|
||||
|
||||
Reference in New Issue
Block a user