wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-08-08 11:26:56 +08:00
Code Issues Packages Projects Releases Wiki Activity

Adjoint solve not yet working.

Browse Source
This commit is contained in:
rowanc1
2014-04-27 23:24:49 -07:00
parent df7676f14c
commit 3514ad56d3
5 changed files with 85 additions and 91 deletions
Download Patch File Download Diff File Expand all files Collapse all files
docs/api_TDEM_derivation.rst
+16 -28
View File
@@ -307,38 +307,21 @@ Multiplying \\(\\mathbf{J}^\\top\\) onto a vector can be broken into three steps
\end{array}
\right]
For the last time-step \\(t=N\\):
For the all time-steps (going backwards in time):
.. math::
\begin{align}
\frac{1}{\delta t} \MfMui \vec{y}_{b}^{(N)} + \MfMui \dcurl \vec{y}_{e}^{(N)} = \vec{p}_b^{(N)} \\
\dcurl^\top \MfMui \vec{y}_b^{(N)} - \MeSig \vec{y}_e^{(N)} = \vec{p}_e^{(N)}
\end{align}
A \vec{y}^{(t)} + B \vec{y}^{(t+1)} = \vec{p}^{(t)}
.. math::
\begin{align}
\left( \MfMui \dcurl \MeSig^{-1} \dcurl^\top \MfMui + \frac{1}{\delta t} \MfMui \right) \vec{y}_{b}^{(N)} = \MfMui \dcurl \MeSig^{-1} \vec{p}_e^{(N)} + \vec{p}_b^{(N)} \\
\vec{y}_e^{(N)} = \MeSig^{-1} \dcurl^\top \MfMui \vec{y}_b^{(N)} - \MeSig^{-1} \vec{p}_e^{(N)}
\end{align}
For the rest of the time-steps (going backwards in time)
.. math::
A \vec{y}^{(t-1)} + B \vec{y}^{(t)} = \vec{p}^{(t-1)}
.. math::
\begin{align}
\frac{1}{\delta t} \MfMui\vec{y}_{b}^{(t-1)} + \MfMui\dcurl \vec{y}_{e}^{(t-1)}
- \frac{1}{\delta t} \MfMui \vec{y}_{b}^{(t)}
= \vec{p}_b^{(t-1)} \\
\dcurl^\top \MfMui \vec{y}_b^{(t-1)} - \MeSig \vec{y}_e^{(t-1)} = \vec{p}_e^{(t-1)}
\frac{1}{\delta t} \MfMui\vec{y}_{b}^{(t)} + \MfMui\dcurl \vec{y}_{e}^{(t)}
- \frac{1}{\delta t} \MfMui \vec{y}_{b}^{(t+1)}
= \vec{p}_b^{(t)} \\
\dcurl^\top \MfMui \vec{y}_b^{(t)} - \MeSig \vec{y}_e^{(t)} = \vec{p}_e^{(t)}
\end{align}
and
@@ -346,8 +329,13 @@ and
.. math::
\begin{align}
\left( \MfMui \dcurl \MeSig^{-1} \dcurl^\top \MfMui + \frac{1}{\delta t} \MfMui \right) \vec{y}_{b}^{(t-1)} =
\frac{1}{\delta t} \MfMui \vec{y}_b^{(t)}
+ \MfMui \dcurl \MeSig^{-1} \vec{p}_e^{(t-1)} + \vec{p}_b^{(t-1)} \\
\vec{y}_e^{(t-1)} = \MeSig^{-1} \dcurl^\top \MfMui \vec{y}_b^{(t-1)} - \MeSig^{-1} \vec{p}_e^{(t-1)}
\left( \MfMui \dcurl \MeSig^{-1} \dcurl^\top \MfMui + \frac{1}{\delta t} \MfMui \right) \vec{y}_{b}^{(t)} =
\frac{1}{\delta t} \MfMui \vec{y}_b^{(t+1)}
+ \MfMui \dcurl \MeSig^{-1} \vec{p}_e^{(t)} + \vec{p}_b^{(t)} \\
\vec{y}_e^{(t)} = \MeSig^{-1} \dcurl^\top \MfMui \vec{y}_b^{(t)} - \MeSig^{-1} \vec{p}_e^{(t)}
\end{align}
.. note::
For the last time step, \\\(t=N\\\), the term: \\\(\\frac{1}{\\delta t} \\MfMui \\vec{y}_b^{(N+1)}\\\) is zero.
simpegEM/TDEM/BaseTDEM.py
+3 -5
View File
@@ -57,8 +57,7 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
return F
def adjoint(self, m, RHS, CalcFields, F=None):
if F is None:
F = FieldsTDEM(self.mesh, self.survey)
F = F or FieldsTDEM(self.mesh, self.survey)
dtFact = None
for tInd, dt in reversed(list(enumerate(self.timeSteps))):
@@ -66,13 +65,12 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
dtFact = dt
A = self.getA(tInd)
# print 'Factoring... (dt = ' + str(dt) + ')'
Asolve = Solver(A, options=self.solverOpts)
Asolve = self.Solver(A, **self.solverOpts)
# print 'Done'
rhs = RHS(tInd, F)
sol = Asolve.solve(rhs)
if sol.ndim == 1:
sol.shape = (sol.size,1)
newFields = CalcFields(sol, self.solType, tInd)
F[:,:,tInd] = newFields
F[:,:,tInd+1] = CalcFields(sol, self.solType, tInd)
return F
simpegEM/TDEM/SurveyTDEM.py
+1 -2
View File
@@ -39,8 +39,7 @@ class RxTDEM(Survey.BaseTimeRx):
if not adjoint:
return P * Utils.mkvc(v[tx, self.projField, :])
elif adjoint:
Ptv = P.T * v[tx, self]
return Ptv
return P.T * v[tx, self]
class FieldsTDEM(Survey.TimeFields):
simpegEM/TDEM/TDEM_b.py
+16 -7
View File
@@ -64,7 +64,7 @@ class ProblemTDEM_b(BaseTDEMProblem):
p = self.survey.projectFieldsDeriv(u, v=v, adjoint=True)
y = self.solveAht(m, p)
w = self.Gtvec(m, y, u)
return w
return - w
def Gvec(self, m, vec, u=None):
"""
@@ -138,18 +138,27 @@ class ProblemTDEM_b(BaseTDEMProblem):
def solveAht(self, m, p):
def AhtRHS(tInd, u):
rhs = self.MfMui*self.mesh.edgeCurl*self.MeSigmaI*p[:,'e',tInd] + p[:,'b',tInd]
nTx, nF = self.survey.nTx, self.mesh.nF
rhs = np.zeros(nF if nTx == 1 else (nF, nTx))
if 'e' in p:
rhs += self.MfMui*self.mesh.edgeCurl*self.MeSigmaI*p[:,'e',tInd]
if 'b' in p:
rhs += p[:,'b',tInd]
if tInd == self.nT-1:
return rhs
dt = self.timeSteps[tInd+1]
dt = self.timeSteps[tInd]
return rhs + 1.0/dt*self.MfMui*u[:,'b',tInd+1]
def AhtCalcFields(sol, solType, tInd):
b = sol
y_b = sol
if self.survey.nTx == 1:
b = mkvc(b)
e = self.MeSigmaI*self.mesh.edgeCurl.T*self.MfMui*b - self.MeSigmaI*p[:,'e',tInd]
return {'b':b, 'e':e}
y_b = mkvc(y_b)
y_e = self.MeSigmaI*self.mesh.edgeCurl.T*self.MfMui*y_b
if 'e' in p:
y_e += - self.MeSigmaI*p[:,'e',tInd]
return {'b':y_b, 'e':y_e}
self.curModel = m
return self.adjoint(m, AhtRHS, AhtCalcFields)
simpegEM/Tests/test_TDEM_b_DerivAdjoint.py
+49 -49
View File
@@ -168,22 +168,22 @@ class TDEM_bDerivTests(unittest.TestCase):
# passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20)
# self.assertTrue(passed)
def test_Deriv_J(self):
# def test_Deriv_J(self):
prb = self.prb
prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)]
mesh = self.mesh
sigma = self.sigma
# prb = self.prb
# prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)]
# mesh = self.mesh
# sigma = self.sigma
# d_sig = 0.8*sigma #np.random.rand(mesh.nCz)
d_sig = 10*np.random.rand(prb.mapping.nP)
# # d_sig = 0.8*sigma #np.random.rand(mesh.nCz)
# d_sig = 10*np.random.rand(prb.mapping.nP)
derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(sigma, mx)]
print '\n'
print 'test_Deriv_J'
passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=4, eps=1e-20)
self.assertTrue(passed)
# derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(sigma, mx)]
# print '\n'
# print 'test_Deriv_J'
# passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=4, eps=1e-20)
# self.assertTrue(passed)
# def test_projectAdjoint(self):
# prb = self.prb
@@ -223,48 +223,48 @@ class TDEM_bDerivTests(unittest.TestCase):
# V2 = f1.tovec().dot(prb.AhtVec(sigma, f2).tovec())
# self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6)
# def test_solveAhtVsAhtVec(self):
# prb = self.prb
# mesh = self.mesh
# sigma = np.random.rand(prb.mapping.nP)
def test_solveAhtVsAhtVec(self):
prb = self.prb
mesh = self.mesh
sigma = np.random.rand(prb.mapping.nP)
# f1 = EM.TDEM.FieldsTDEM(mesh, 1, prb.nT, 'b')
# for i in range(prb.nT):
# f1.set_b(np.random.rand(mesh.nF, 1), i)
# f1.set_e(np.random.rand(mesh.nE, 1), i)
f1 = EM.TDEM.FieldsTDEM(mesh,prb.survey)
for i in range(prb.nT):
f1[:,'b',i] = np.random.rand(mesh.nF, 1)
f1[:,'e',i] = np.random.rand(mesh.nE, 1)
# f2 = prb.solveAht(sigma, f1)
# f3 = prb.AhtVec(sigma, f2)
f2 = prb.solveAht(sigma, f1)
f3 = prb.AhtVec(sigma, f2)
# if plotIt:
# import matplotlib.pyplot as plt
# plt.plot(f3.tovec())
# plt.plot(f1.tovec())
# plt.show()
# V1 = np.linalg.norm(f3.tovec()-f1.tovec())
# V2 = np.linalg.norm(f1.tovec())
# print V1, V2
# print 'I am gunna fail this one: boo. :('
# self.assertLess(V1/V2, 1e-6)
if plotIt:
import matplotlib.pyplot as plt
plt.plot(f3.tovec())
plt.plot(f1.tovec())
plt.show()
V1 = np.linalg.norm(f3.tovec()-f1.tovec())
V2 = np.linalg.norm(f1.tovec())
print V1, V2
print 'I am gunna fail this one: boo. :('
self.assertLess(V1/V2, 1e-6)
# def test_adjointsolveAhVssolveAht(self):
# prb = self.prb
# mesh = self.mesh
# sigma = self.sigma
def test_adjointsolveAhVssolveAht(self):
prb = self.prb
mesh = self.mesh
sigma = self.sigma
# f1 = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
# for i in range(1,prb.nT+1):
# f1[:,'b',i] = np.random.rand(mesh.nF, 1)
# f1[:,'e',i] = np.random.rand(mesh.nE, 1)
f1 = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
for i in range(1,prb.nT+1):
f1[:,'b',i] = np.random.rand(mesh.nF, 1)
f1[:,'e',i] = np.random.rand(mesh.nE, 1)
# f2 = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
# for i in range(1,prb.nT+1):
# f2[:,'b',i] = np.random.rand(mesh.nF, 1)
# f2[:,'e',i] = np.random.rand(mesh.nE, 1)
f2 = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
for i in range(1,prb.nT+1):
f2[:,'b',i] = np.random.rand(mesh.nF, 1)
f2[:,'e',i] = np.random.rand(mesh.nE, 1)
# V1 = f2.tovec().dot(prb.solveAh(sigma, f1).tovec())
# V2 = f1.tovec().dot(prb.solveAht(sigma, f2).tovec())
# self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6)
V1 = f2.tovec().dot(prb.solveAh(sigma, f1).tovec())
V2 = f1.tovec().dot(prb.solveAht(sigma, f2).tovec())
self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6)
# def test_adjointGvecVsGtvec(self):
# mesh = self.mesh
@@ -274,12 +274,12 @@ class TDEM_bDerivTests(unittest.TestCase):
# sigma = np.random.rand(prb.mapping.nP)
# u = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
# for i in range(prb.nT):
# for i in range(1,prb.nT+1):
# u[:,'b',i] = np.random.rand(mesh.nF, 1)
# u[:,'e',i] = np.random.rand(mesh.nE, 1)
# v = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
# for i in range(prb.nT):
# for i in range(1,prb.nT+1):
# v[:,'b',i] = np.random.rand(mesh.nF, 1)
# v[:,'e',i] = np.random.rand(mesh.nE, 1)