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

Updates to derivs for multiple Txs

Browse Source
This commit is contained in:
rowanc1
2014-04-28 12:34:44 -07:00
parent 4c15267c0c
commit 8e5c93accb
4 changed files with 165 additions and 6 deletions
Download Patch File Download Diff File Expand all files Collapse all files
simpegEM/TDEM/SurveyTDEM.py
+2 -2
View File
@@ -60,8 +60,8 @@ class FieldsTDEM(Survey.TimeFields):
e = self[:,'e',i+1]
else:
e = np.zeros(nE if nTx == 1 else (nE, nTx))
u = np.r_[u, b, e]
return u
u = np.concatenate((u, b, e))
return Utils.mkvc(u)
class TxTDEM(Survey.BaseTx):
rxPair = RxTDEM
simpegEM/TDEM/TDEM_b.py
+12 -3
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 - mkvc(w)
def Gvec(self, m, vec, u=None):
"""
@@ -107,10 +107,13 @@ class ProblemTDEM_b(BaseTDEMProblem):
if u is None:
u = self.fields(m)
nTx, nE = self.survey.nTx, self.mesh.nE
tmp = np.zeros(nE if nTx == 1 else (nE,nTx))
tmp = np.zeros(nE)
# Here we can do internal multiplications of Gt*v and then multiply by MsigDeriv.T in one go.
for i in range(1,self.nT+1):
tmp += vec[:,'e',i]*u[:,'e',i]
vu = vec[:,'e',i]*u[:,'e',i]
if nTx > 1:
vu = vu.sum(axis=1)
tmp += vu
curModel = self.mapping.transform(m)
p = -mkvc(self.mapping.transformDeriv(m).T*self.mesh.getEdgeInnerProductDeriv(curModel).T*tmp)
@@ -137,6 +140,12 @@ class ProblemTDEM_b(BaseTDEMProblem):
def solveAht(self, m, p):
# Mini Example:
#
# nT = 3, len(times) == 4, fields stored in F[:,:,1:4]
#
# 0 is held for initial conditions (this shifts the storage by +1)
# ^
# fLoc 0 1 2 3
# |-----|-----|-----|
# tInd 0 1 2 / /
simpegEM/Tests/test_TDEM_b_DerivAdjoint.py
+1 -1
View File
@@ -197,7 +197,7 @@ class TDEM_bDerivTests(unittest.TestCase):
for i in range(prb.nT):
f[:,'b',i] = np.random.rand(mesh.nF, 1)
f[:,'e',i] = np.random.rand(mesh.nE, 1)
d_vec = np.random.rand(survey.nD, survey.nTx).flatten()
d_vec = np.random.rand(survey.nD)
d = Survey.Data(survey,v=d_vec)
# Check that d.T*Q*f = f.T*Q.T*d
simpegEM/Tests/test_TDEM_b_MultiTx_DerivAdjoint.py
+150
View File
@@ -0,0 +1,150 @@
import unittest
from SimPEG import *
import simpegEM as EM
plotIt = False
class TDEM_bDerivTests(unittest.TestCase):
def setUp(self):
cs = 5.
ncx = 20
ncy = 6
npad = 20
hx = [(cs,ncx), (cs,npad,1.3)]
hy = [(cs,npad,-1.3), (cs,ncy), (cs,npad,1.3)]
mesh = Mesh.CylMesh([hx,1,hy], '00C')
active = mesh.vectorCCz<0.
activeMap = Maps.ActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
mapping = Maps.ComboMap(mesh,
[Maps.ExpMap, Maps.Vertical1DMap, activeMap])
rxOffset = 40.
rx = EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20), 'bz')
tx = EM.TDEM.TxTDEM(np.array([0., 0., 0.]), 'VMD_MVP', [rx])
rx2 = EM.TDEM.RxTDEM(np.array([[rxOffset-10, 0., 0.]]), np.logspace(-5,-4, 25), 'bz')
tx2 = EM.TDEM.TxTDEM(np.array([0., 0., 0.]), 'VMD_MVP', [rx2])
survey = EM.TDEM.SurveyTDEM([tx,tx2])
self.prb = EM.TDEM.ProblemTDEM_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)]
self.sigma = np.ones(mesh.nCz)*1e-8
self.sigma[mesh.vectorCCz<0] = 1e-1
self.sigma = np.log(self.sigma[active])
self.prb.pair(survey)
self.mesh = mesh
def test_DerivG(self):
"""
Test the derivative of c with respect to sigma
"""
# Random model and perturbation
sigma = np.random.rand(self.prb.mapping.nP)
f = self.prb.fields(sigma)
dm = 1000*np.random.rand(self.prb.mapping.nP)
h = 0.01
derChk = lambda m: [self.prb.AhVec(m, f).tovec(), lambda mx: self.prb.Gvec(sigma, mx, u=f).tovec()]
print '\ntest_DerivG'
passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20)
self.assertTrue(passed)
def test_Deriv_dUdM(self):
prb = self.prb
prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)]
mesh = self.mesh
sigma = self.sigma
dm = 10*np.random.rand(prb.mapping.nP)
f = prb.fields(sigma)
derChk = lambda m: [self.prb.fields(m).tovec(), lambda mx: -prb.solveAh(sigma, prb.Gvec(sigma, mx, u=f)).tovec()]
print '\n'
print 'test_Deriv_dUdM'
passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20)
self.assertTrue(passed)
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
# 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)
def test_projectAdjoint(self):
prb = self.prb
survey = prb.survey
mesh = self.mesh
# Generate random fields and data
f = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
for i in range(prb.nT):
f[:,'b',i] = np.random.rand(mesh.nF, 1)
f[:,'e',i] = np.random.rand(mesh.nE, 1)
d_vec = np.random.rand(survey.nD)
d = Survey.Data(survey,v=d_vec)
# Check that d.T*Q*f = f.T*Q.T*d
V1 = d_vec.dot(survey.projectFieldsDeriv(None, v=f).tovec())
V2 = f.tovec().dot(survey.projectFieldsDeriv(None, v=d, adjoint=True).tovec())
self.assertLess((V1-V2)/np.abs(V1), 1e-6)
def test_adjointGvecVsGtvec(self):
mesh = self.mesh
prb = self.prb
m = np.random.rand(prb.mapping.nP)
sigma = np.random.rand(prb.mapping.nP)
u = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
for i in range(1,prb.nT+1):
u[:,'b',i] = np.random.rand(mesh.nF, 2)
u[:,'e',i] = np.random.rand(mesh.nE, 2)
v = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
for i in range(1,prb.nT+1):
v[:,'b',i] = np.random.rand(mesh.nF, 2)
v[:,'e',i] = np.random.rand(mesh.nE, 2)
V1 = m.dot(prb.Gtvec(sigma, v, u))
V2 = v.tovec().dot(prb.Gvec(sigma, m, u).tovec())
self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6)
def test_adjointJvecVsJtvec(self):
mesh = self.mesh
prb = self.prb
sigma = self.sigma
m = np.random.rand(prb.mapping.nP)
d = np.random.rand(prb.survey.nD)
V1 = d.dot(prb.Jvec(sigma, m))
V2 = m.dot(prb.Jtvec(sigma, d))
print 'AdjointTest', V1, V2
self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6)
if __name__ == '__main__':
unittest.main()