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

working jvec

Browse Source
This commit is contained in:
rowanc1
2014-04-27 22:17:55 -07:00
parent d1f23081d6
commit df7676f14c
3 changed files with 41 additions and 52 deletions
Download Patch File Download Diff File Expand all files Collapse all files
docs/api_TDEM_derivation.rst
+6 -18
View File
@@ -255,25 +255,8 @@ Multiplying **J** onto a vector can be broken into three steps
\vec{p}_e^{(n)} = - \diag{\e^{(n)}} \Ace \diag{V} m
\end{align}
First time step
.. math::
\begin{align}
\frac{1}{\delta t} \MfMui \vec{y}_{b}^{(1)} + \MfMui \dcurl \vec{y}_{e}^{(1)} = \vec{p}_b^{(1)} \\
\dcurl^\top \MfMui \vec{y}_b^{(1)} - \MeSig \vec{y}_e^{(1)} = \vec{p}_e^{(1)}
\end{align}
.. math::
\begin{align}
\left( \MfMui \dcurl \MeSig^{-1} \dcurl^\top \MfMui + \frac{1}{\delta t} \MfMui \right) \vec{y}_{b}^{(1)} = \MfMui \dcurl \MeSig^{-1} \vec{p}_e^{(1)} + \vec{p}_b^{(1)} \\
\vec{y}_e^{(1)} = \MeSig^{-1} \dcurl^\top \MfMui \vec{y}_b^{(1)} - \MeSig^{-1} \vec{p}_e^{(1)}
\end{align}
Remaining time steps:
For all time steps:
.. math::
@@ -295,6 +278,11 @@ and
\vec{y}_e^{(t+1)} = \MeSig^{-1} \dcurl^\top \MfMui \vec{y}_b^{(t+1)} - \MeSig^{-1} \vec{p}_e^{(t+1)}
\end{align}
.. note::
For the first time step, \\\(t=0\\\), the term: \\\(\\frac{1}{\\delta t} \\MfMui \\vec{y}_b^{(0)}\\\) is zero.
Implementing **J** transpose times a vector
simpegEM/TDEM/TDEM_b.py
+8 -7
View File
@@ -52,7 +52,7 @@ class ProblemTDEM_b(BaseTDEMProblem):
p = self.Gvec(m, v, u)
y = self.solveAh(m, p)
Jv = self.survey.projectFieldsDeriv(u, v=y)
return mkvc(Jv)
return - mkvc(Jv)
def Jtvec(self, m, v, u=None):
if u is None:
@@ -117,19 +117,20 @@ class ProblemTDEM_b(BaseTDEMProblem):
return p
def solveAh(self, m, p):
def AhRHS(tInd, u):
def AhRHS(tInd, y):
rhs = self.MfMui*self.mesh.edgeCurl*self.MeSigmaI*p[:,'e',tInd+1] + p[:,'b',tInd+1]
if tInd == 0:
return rhs
dt = self.timeSteps[tInd]
return rhs + 1.0/dt*self.MfMui*u[:,'b',tInd]
return rhs + 1.0/dt*self.MfMui*y[:,'b',tInd]
def AhCalcFields(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+1]
return {'b':b, 'e':e}
y_b = mkvc(y_b)
y_e = self.MeSigmaI*self.mesh.edgeCurl.T*self.MfMui*y_b - self.MeSigmaI*p[:,'e',tInd+1]
return {'b':y_b, 'e':y_e}
self.curModel = m
return self.forward(m, AhRHS, AhCalcFields)
simpegEM/Tests/test_TDEM_b_DerivAdjoint.py
+27 -27
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
@@ -247,24 +247,24 @@ class TDEM_bDerivTests(unittest.TestCase):
# 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