mirror of
https://github.com/wassname/simpeg.git
synced 2026-06-29 21:09:14 +08:00
adjoint hooked up for b formulation
This commit is contained in:
@@ -86,15 +86,12 @@ class Fields_b(Fields):
|
||||
_, S_e = src.eval(self.survey.prob, self.survey.prob.times[tInd])
|
||||
bSolution = self[[src],'bSolution',tInd]
|
||||
|
||||
if adjoint:
|
||||
v = self.MeSigmaI.T * v
|
||||
|
||||
_, S_eDeriv_v = src.evalDeriv(self.survey.prob.times[tInd], self, v=v, adjoint=adjoint)
|
||||
_, S_eDeriv = src.evalDeriv(self.survey.prob.times[tInd], self, adjoint=adjoint)
|
||||
|
||||
if adjoint is True:
|
||||
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 ) ).T * v - S_eDeriv(self.MeSigmaI.T * v)
|
||||
|
||||
return self.MeSigmaIDeriv(-S_e + self.edgeCurl.T * ( self.MfMui * bSolution)) * v - self.MeSigmaI * 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:
|
||||
|
||||
+39
-12
@@ -168,7 +168,8 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
|
||||
df_duT_v[src, '%sDeriv'%self._fieldType, :] = np.zeros_like(u[src, self._fieldType, :])
|
||||
|
||||
for rx in src.rxList:
|
||||
PT_v[src,'%sDeriv'%rx.projField,:] = rx.evalDeriv(src, self.mesh, self.timeMesh, Utils.mkvc(v[src,rx]), adjoint=True)
|
||||
PT_v[src,'%sDeriv'%rx.projField,:] = rx.evalDeriv(src, self.mesh, self.timeMesh, Utils.mkvc(v[src,rx]), adjoint=True) # this is +=
|
||||
|
||||
# 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)
|
||||
@@ -187,13 +188,14 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
|
||||
AdiagTinv = None
|
||||
|
||||
# Do the back-solve through time
|
||||
for tInd in reversed(range(self.nT)):
|
||||
for tIndP in reversed(range(self.nT + 1)):
|
||||
tInd = tIndP - 1
|
||||
if AdiagTinv is not None and (tInd <= self.nT and self.timeSteps[tInd] != self.timeSteps[tInd+1]): # if the previous timestep is the same --> no need to refactor the matrix
|
||||
AdiagTinv.clean()
|
||||
AdiagTinv = None
|
||||
|
||||
# refactor if we need to
|
||||
if AdiagTinv is None:
|
||||
if AdiagTinv is None and tInd > -1:
|
||||
Adiag = self.getAdiag(tInd)
|
||||
AdiagTinv = self.Solver(Adiag.T, **self.solverOpts)
|
||||
|
||||
@@ -204,25 +206,46 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
|
||||
for isrc, src in enumerate(self.survey.srcList):
|
||||
# solve against df_duT_v
|
||||
if tInd >= self.nT-1:
|
||||
# last timestep (first to be solved)
|
||||
ATinv_df_duT_v[isrc,:] = AdiagTinv * df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1]
|
||||
else:
|
||||
elif tInd > -1:
|
||||
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,:]))
|
||||
|
||||
un_src = u[src,ftype,tInd+1]
|
||||
dAT_dm_v = self.getAdiagDeriv(None, un_src, ATinv_df_duT_v[isrc,:], adjoint=True) # cell centered on time mesh
|
||||
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'
|
||||
|
||||
dRHST_dm_v = self.getRHSDeriv(tInd+1, src, ATinv_df_duT_v[isrc,:], adjoint=True) # on nodes of time mesh
|
||||
# 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_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)
|
||||
# JTv = JTv + Utils.mkvc(dRHST_dm_v0 + dRHST_dm_v1)
|
||||
|
||||
# print 'here'
|
||||
# inFields = self.getInitialFieldsDeriv(u[src,ftype,tInd+1], Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1]), adjoint=True)
|
||||
# # - Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:]), adjoint=True)
|
||||
# print inFields.shape
|
||||
# JTv = JTv + inFields
|
||||
# dAsubdiag_dm_v = 0
|
||||
|
||||
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:
|
||||
for projField in set(rx.projField):
|
||||
JTv = JTv + self.getInitialFieldsDeriv(Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,0]), adjoint=True)
|
||||
# Asubdiag = self.getAsubdiag(0)
|
||||
# for src in self.survey.srcList:
|
||||
# for projField in set(rx.projField):
|
||||
# v = AdiagTinv * (Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,0]) - Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:]))
|
||||
# JTv = JTv - Utils.mkvc(self.getAdiagDeriv(0, u[src, ftype, tInd], v, adjoint = True))
|
||||
# # JTv = JTv + self.getInitialFieldsDeriv(Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,0] - Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:])), adjoint=True)
|
||||
|
||||
|
||||
return Utils.mkvc(JTv)
|
||||
@@ -274,7 +297,10 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
|
||||
ifieldsDeriv[:,i] = ifieldsDeriv[:,i] + getattr(src, '%sInitialDeriv'%self._fieldType, None)(self,v,adjoint)
|
||||
|
||||
if adjoint is True:
|
||||
ifieldsDeriv = ifieldsDeriv.sum()
|
||||
ifieldsDeriv = np.zeros_like(self.curModel)
|
||||
print v.shape
|
||||
# ifieldsDeriv = self.getAdiagDeriv(None, u, v, adjoint)
|
||||
# ifieldsDeriv = ifieldsDeriv.sum()
|
||||
|
||||
return ifieldsDeriv
|
||||
|
||||
@@ -337,6 +363,7 @@ class Problem_b(BaseTDEMProblem):
|
||||
(\mathbf{I} + \mathbf{dt} \mathbf{C} \mathbf{M_{\sigma}^e}^{-1} \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f})
|
||||
|
||||
"""
|
||||
assert tInd >= 0 and tInd < self.nT
|
||||
|
||||
dt = self.timeSteps[tInd]
|
||||
C = self.mesh.edgeCurl
|
||||
|
||||
@@ -23,7 +23,7 @@ def setUp(rxcomp='bz'):
|
||||
mapping = Maps.ExpMap(mesh) * Maps.SurjectVertical1D(mesh) * activeMap
|
||||
|
||||
rxOffset = 10.
|
||||
rx = EM.TDEM.Rx(np.array([[rxOffset, 0., -1e-2]]), np.logspace(-4,-3, 20), rxcomp)
|
||||
rx = EM.TDEM.Rx(np.array([[rxOffset, 0., -1e-2]]), np.logspace(-4,-3, 20), rxcomp) #,]
|
||||
src = EM.TDEM.SurveyTDEM.MagDipole([rx], loc=np.array([0., 0., 0.]))
|
||||
|
||||
survey = EM.TDEM.Survey([src])
|
||||
@@ -45,13 +45,14 @@ def setUp(rxcomp='bz'):
|
||||
|
||||
return prb, m, mesh
|
||||
|
||||
|
||||
class TDEM_DerivTests(unittest.TestCase):
|
||||
|
||||
# ====== TEST A ========== #
|
||||
|
||||
def test_Deriv_Pieces(self):
|
||||
def test_AderivTest(self):
|
||||
prb, m0, mesh = setUp()
|
||||
tInd = 0
|
||||
tInd = 2
|
||||
|
||||
v = np.random.rand(mesh.nF)
|
||||
|
||||
@@ -64,36 +65,78 @@ class TDEM_DerivTests(unittest.TestCase):
|
||||
|
||||
return Av, ADeriv_dm
|
||||
|
||||
def A_adjointTest():
|
||||
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)
|
||||
prb.curModel = m0
|
||||
|
||||
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) < TOL * (np.abs(V1) + np.abs(V2))/2.
|
||||
print 'AdjointTest', V1, V2, passed
|
||||
self.assertTrue(passed)
|
||||
|
||||
print '\n Testing ADeriv'
|
||||
Tests.checkDerivative(AderivTest, m0, plotIt=False, num=4, eps=1e-20)
|
||||
A_adjointTest()
|
||||
|
||||
def test_A_adjointTest(self):
|
||||
prb, m0, mesh = setUp()
|
||||
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)
|
||||
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
|
||||
self.assertTrue(passed)
|
||||
|
||||
# ====== TEST Fields Deriv Pieces ========== #
|
||||
|
||||
def test_eDeriv_m_adjoint(self):
|
||||
prb, m0, mesh = setUp()
|
||||
tInd = 0
|
||||
|
||||
v = np.random.rand(mesh.nF)
|
||||
|
||||
print '\n Testing eDeriv_m Adjoint'
|
||||
|
||||
prb, m0, mesh = setUp()
|
||||
f = prb.fields(m0)
|
||||
|
||||
m = np.random.rand(prb.mapping.nP)
|
||||
e = np.random.randn(prb.mesh.nE)
|
||||
V1 = e.dot(f._eDeriv_m(1, prb.survey.srcList[0], m))
|
||||
V2 = m.dot(f._eDeriv_m(1, prb.survey.srcList[0], e, adjoint=True))
|
||||
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)
|
||||
|
||||
def test_eDeriv_u_adjoint(self):
|
||||
print '\n Testing eDeriv_u Adjoint'
|
||||
|
||||
prb, m0, mesh = setUp()
|
||||
f = prb.fields(m0)
|
||||
|
||||
b = np.random.rand(prb.mesh.nF)
|
||||
e = np.random.randn(prb.mesh.nE)
|
||||
V1 = e.dot(f._eDeriv_u(1, prb.survey.srcList[0], b))
|
||||
V2 = b.dot(f._eDeriv_u(1, prb.survey.srcList[0], e, adjoint=True))
|
||||
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)
|
||||
|
||||
|
||||
# ====== TEST Jvec ========== #
|
||||
|
||||
def JvecTest(self, rxcomp):
|
||||
prb, m, mesh = setUp(rxcomp)
|
||||
|
||||
derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(m, mx)]
|
||||
print '\n'
|
||||
print 'test_Jvec_%s' %(rxcomp)
|
||||
Tests.checkDerivative(derChk, m, plotIt=False, num=2, eps=1e-20)
|
||||
|
||||
if testDeriv:
|
||||
|
||||
def JvecTest(self, rxcomp):
|
||||
prb, m, mesh = setUp(rxcomp)
|
||||
|
||||
derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(m, mx)]
|
||||
print '\n'
|
||||
print 'test_Jvec_%s' %(rxcomp)
|
||||
Tests.checkDerivative(derChk, m, plotIt=False, num=2, eps=1e-20)
|
||||
|
||||
def test_Jvec_b_bx(self):
|
||||
self.JvecTest('bx')
|
||||
|
||||
@@ -103,33 +146,38 @@ class TDEM_DerivTests(unittest.TestCase):
|
||||
def test_Jvec_b_ey(self):
|
||||
self.JvecTest('ey')
|
||||
|
||||
else:
|
||||
pass
|
||||
|
||||
|
||||
# ====== TEST Jtvec ========== #
|
||||
|
||||
def adjointJvecVsJtvecTest(self, rxcomp='bz'):
|
||||
|
||||
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))
|
||||
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:
|
||||
def test_Jvec_b_bx(self):
|
||||
self.adjointJvecVsJtvecTest('bx')
|
||||
|
||||
def test_Jvec_b_bz(self):
|
||||
self.adjointJvecVsJtvecTest('bz')
|
||||
def JvecVsJtvecTest(self, rxcomp='bz'):
|
||||
|
||||
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))
|
||||
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)
|
||||
|
||||
def test_Jvec_adjoint_b_bx(self):
|
||||
self.JvecVsJtvecTest('bx')
|
||||
|
||||
def test_Jvec_adjoint_b_bz(self):
|
||||
self.JvecVsJtvecTest('bz')
|
||||
|
||||
def test_Jvec_adjoint_b_ey(self):
|
||||
self.JvecVsJtvecTest('ey')
|
||||
|
||||
def test_Jvec_b_ey(self):
|
||||
self.adjointJvecVsJtvecTest('ey')
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
|
||||
@@ -6,11 +6,9 @@ plotIt = False
|
||||
testDeriv = True
|
||||
testAdjoint = True
|
||||
|
||||
TOL = 1e-6
|
||||
TOL = 1e-5
|
||||
|
||||
class TDEM_bDerivTests(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
def setUp(self, rxcomp='bz'):
|
||||
|
||||
cs = 5.
|
||||
ncx = 20
|
||||
@@ -25,46 +23,58 @@ class TDEM_bDerivTests(unittest.TestCase):
|
||||
mapping = Maps.ExpMap(mesh) * Maps.SurjectVertical1D(mesh) * activeMap
|
||||
|
||||
rxOffset = 40.
|
||||
rx = EM.TDEM.Rx(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20), 'bz')
|
||||
rx = EM.TDEM.Rx(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20), rxcomp)
|
||||
src = EM.TDEM.SurveyTDEM.MagDipole( [rx], loc=np.array([0., 0., 0.]))
|
||||
rx2 = EM.TDEM.Rx(np.array([[rxOffset-10, 0., 0.]]), np.logspace(-5,-4, 25), 'bz')
|
||||
rx2 = EM.TDEM.Rx(np.array([[rxOffset-10, 0., 0.]]), np.logspace(-5,-4, 25), rxcomp)
|
||||
src2 = EM.TDEM.SurveyTDEM.MagDipole( [rx2], loc=np.array([0., 0., 0.]))
|
||||
|
||||
survey = EM.TDEM.Survey([src,src2])
|
||||
|
||||
self.prb = EM.TDEM.Problem_b(mesh, mapping=mapping)
|
||||
# self.prb.timeSteps = [1e-5]
|
||||
self.prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)]
|
||||
# self.prb.timeSteps = [(1e-05, 100)]
|
||||
prb = EM.TDEM.Problem_b(mesh, mapping=mapping)
|
||||
# prb.timeSteps = [1e-5]
|
||||
prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)]
|
||||
# prb.timeSteps = [(1e-05, 100)]
|
||||
|
||||
try:
|
||||
from pymatsolver import MumpsSolver
|
||||
self.prb.Solver = MumpsSolver
|
||||
prb.Solver = MumpsSolver
|
||||
except ImportError, e:
|
||||
self.prb.Solver = SolverLU
|
||||
prb.Solver = SolverLU
|
||||
|
||||
self.m = np.log(1e-1)*np.ones(self.prb.mapping.nP) + 1e-2*np.random.randn(self.prb.mapping.nP)
|
||||
m = np.log(1e-1)*np.ones(prb.mapping.nP) + 1e-2*np.random.randn(prb.mapping.nP)
|
||||
|
||||
prb.pair(survey)
|
||||
|
||||
return mesh, prb, m
|
||||
|
||||
class TDEM_bDerivTests(unittest.TestCase):
|
||||
|
||||
self.prb.pair(survey)
|
||||
self.mesh = mesh
|
||||
|
||||
if testDeriv:
|
||||
def test_Deriv_J(self):
|
||||
def Deriv_J(self, rxcomp='bz'):
|
||||
|
||||
mesh, prb, m0 = setUp(rxcomp)
|
||||
|
||||
prb = self.prb
|
||||
prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)]
|
||||
mesh = self.mesh
|
||||
|
||||
derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(self.m, mx)]
|
||||
derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(m0, mx)]
|
||||
print '\n'
|
||||
print 'test_Deriv_J'
|
||||
Tests.checkDerivative(derChk, self.m, plotIt=False, num=3, eps=1e-20)
|
||||
print 'test_Deriv_J %s'%rxcomp
|
||||
Tests.checkDerivative(derChk, m0, plotIt=False, num=3, eps=1e-20)
|
||||
|
||||
def test_Jvec_bx(self):
|
||||
self.Deriv_J('bx')
|
||||
|
||||
def test_Jvec_bz(self):
|
||||
self.Deriv_J('bz')
|
||||
|
||||
def test_Jvec_ey(self):
|
||||
self.Deriv_J('ey')
|
||||
|
||||
if testAdjoint:
|
||||
def test_adjointJvecVsJtvec(self):
|
||||
mesh = self.mesh
|
||||
prb = self.prb
|
||||
m0 = self.m
|
||||
def adjointJvecVsJtvec(self, rxcomp='bz'):
|
||||
print ' \n Testing Adjoint %s' %rxcomp
|
||||
mesh, prb, m0 = setUp(rxcomp)
|
||||
|
||||
m = np.random.rand(prb.mapping.nP)
|
||||
d = np.random.rand(prb.survey.nD)
|
||||
@@ -77,6 +87,15 @@ class TDEM_bDerivTests(unittest.TestCase):
|
||||
print ' ', V1, V2, np.abs(V1-V2), tol, passed
|
||||
self.assertTrue(passed)
|
||||
|
||||
def test_JvecVsJtvec_bx(self):
|
||||
self.adjointJvecVsJtvec('bx')
|
||||
|
||||
def test_JvecVsJtvec_bz(self):
|
||||
self.adjointJvecVsJtvec('bz')
|
||||
|
||||
def test_JvecVsJtvec_ey(self):
|
||||
self.adjointJvecVsJtvec('ey')
|
||||
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
|
||||
Reference in New Issue
Block a user