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

Added Gtvec method to TDEM. Renamed G to Gvec. Added tests.

Browse Source
This commit is contained in:
Dave Marchant
2014-02-18 20:34:34 -08:00
parent 4cd25fde82
commit edc241eec2
2 changed files with 36 additions and 6 deletions
Download Patch File Download Diff File Expand all files Collapse all files
simpegEM/TDEM/TDEM_b.py
+14 -3
View File
@@ -1,5 +1,6 @@
from BaseTDEM import ProblemBaseTDEM
from FieldsTDEM import FieldsTDEM
from SimPEG.Utils import mkvc
import numpy as np
class ProblemTDEM_b(ProblemBaseTDEM):
@@ -43,15 +44,15 @@ class ProblemTDEM_b(ProblemBaseTDEM):
def Jvec(self, m, v, u=None):
if u is None:
u = self.fields(m)
p = self.G(m, v, u)
p = self.Gvec(m, v, u)
y = self.solveAh(m, p)
return self.data.dpred(m, u=y)
def G(self, m, v, u=None):
def Gvec(self, m, v, u=None):
if u is None:
u = self.fields(m)
p = FieldsTDEM(self.mesh, 1, self.times.size, 'b')
c = self.mesh.getEdgeMassDeriv()*self.model.transformDeriv(m)*v
c = self.mesh.getEdgeMassDeriv()*self.model.transformDeriv(None)*v
for i in range(self.times.size):
ei = u.get_e(i)
pVal = np.empty_like(ei)
@@ -62,6 +63,16 @@ class ProblemTDEM_b(ProblemBaseTDEM):
p.set_b(np.zeros((self.mesh.nF,1)), i)
return p
def Gtvec(self, m, v, u=None):
if u is None:
u = self.fields(m)
tmp = np.zeros((self.mesh.nE,self.data.nTx))
for i in range(self.nTimes):
tmp += v.get_e(i)*u.get_e(i)
p = -mkvc(self.model.transformDeriv(None).T*self.mesh.getEdgeMassDeriv().T*tmp)
return p
def solveAh(self, m, p):
def AhRHS(tInd, u):
rhs = self.MfMui*self.mesh.edgeCurl*self.MeSigmaI*p.get_e(tInd) + p.get_b(tInd)
simpegEM/Tests/test_forward_EMproblem.py
+22 -3
View File
@@ -180,7 +180,7 @@ class TDEM_bDerivTests(unittest.TestCase):
h = 1.
a = np.linalg.norm(self.prb.AhVec(sigma+h*dm, f).fieldVec() - self.prb.AhVec(sigma, f).fieldVec())
b = np.linalg.norm(self.prb.AhVec(sigma+h*dm, f).fieldVec() - self.prb.AhVec(sigma, f).fieldVec() - h*self.prb.G(sigma, dm, u=f).fieldVec())
b = np.linalg.norm(self.prb.AhVec(sigma+h*dm, f).fieldVec() - self.prb.AhVec(sigma, f).fieldVec() - h*self.prb.Gvec(sigma, dm, u=f).fieldVec())
# Assuming that the gradient is exact to machine precision
self.assertTrue(b<1e-16)
@@ -202,7 +202,7 @@ class TDEM_bDerivTests(unittest.TestCase):
for i, h in enumerate(hv):
f = prb.fields(sigma)
fstep = prb.fields(sigma + h*d_sig)
dcdm = prb.G(sigma, h*d_sig, u=f) # TODO: make negative!?!?
dcdm = prb.Gvec(sigma, h*d_sig, u=f) # TODO: make negative!?!?
dudm = prb.solveAh(sigma, dcdm)
linear = np.linalg.norm(f.fieldVec() - fstep.fieldVec())
@@ -305,7 +305,26 @@ class TDEM_bDerivTests(unittest.TestCase):
V2 = f1.fieldVec().dot(prb.solveAht(sigma, f2).fieldVec())
self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-8)
def test_adjointGvecVsGtvec(self):
mesh = self.mesh
prb = self.prb
m = np.random.rand(mesh.nCz)
sigma = np.random.rand(mesh.nCz)
u = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nTimes, 'b')
for i in range(u.nTimes):
u.set_b(np.random.rand(mesh.nF, 1), i)
u.set_e(np.random.rand(mesh.nE, 1), i)
v = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nTimes, 'b')
for i in range(v.nTimes):
v.set_b(np.random.rand(mesh.nF, 1), i)
v.set_e(np.random.rand(mesh.nE, 1), i)
V1 = m.dot(prb.Gtvec(sigma, v, u))
V2 = v.fieldVec().dot(prb.Gvec(sigma, m, u).fieldVec())
self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-8)
if __name__ == '__main__':