mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-17 11:32:59 +08:00
jtvec runs, fails
This commit is contained in:
Lindsey Heagy
@@ -31,6 +31,8 @@ class Fields(SimPEG.Problem.TimeFields):
|
||||
knownFields = {}
|
||||
dtype = float
|
||||
|
||||
|
||||
|
||||
class Fields_Derivs(Fields):
|
||||
knownFields = {
|
||||
'bDeriv': 'F',
|
||||
@@ -65,7 +67,7 @@ class Fields_b(Fields):
|
||||
|
||||
def _bDeriv(self, tInd, src, dun_dm_v, v, adjoint=False):
|
||||
if adjoint is True:
|
||||
raise NotImplementedError
|
||||
return self._bDeriv_u(tInd, src, v, adjoint), self._bDeriv_m(tInd, src, v, adjoint)
|
||||
return self._bDeriv_u(tInd, src, dun_dm_v) + self._bDeriv_m(tInd, src, v)
|
||||
|
||||
def _e(self, bSolution, srcList, tInd):
|
||||
|
||||
@@ -252,24 +252,25 @@ class Survey(SimPEG.Survey.BaseSurvey):
|
||||
return data
|
||||
|
||||
def evalDeriv(self, u, v=None, adjoint=False):
|
||||
assert v is not None, 'v to multiply must be provided.'
|
||||
raise Exception('Use Receivers to project fields deriv.')
|
||||
# assert v is not None, 'v to multiply must be provided.'
|
||||
|
||||
if not adjoint:
|
||||
data = SimPEG.Survey.Data(self)
|
||||
for src in self.srcList:
|
||||
for rx in src.rxList:
|
||||
data[src, rx] = rx.evalDeriv(src, self.mesh, self.prob.timeMesh, u, v)
|
||||
return data
|
||||
else:
|
||||
f = FieldsTDEM(self.mesh, self)
|
||||
for src in self.srcList:
|
||||
for rx in src.rxList:
|
||||
Ptv = rx.evalDeriv(src, self.mesh, self.prob.timeMesh, u, v, adjoint=True)
|
||||
Ptv = Ptv.reshape((-1, self.prob.timeMesh.nN), order='F')
|
||||
if rx.projField not in f: # first time we are projecting
|
||||
f[src, rx.projField, :] = Ptv
|
||||
else: # there are already fields, so let's add to them!
|
||||
f[src, rx.projField, :] += Ptv
|
||||
return f
|
||||
# if not adjoint:
|
||||
# data = SimPEG.Survey.Data(self)
|
||||
# for src in self.srcList:
|
||||
# for rx in src.rxList:
|
||||
# data[src, rx] = rx.evalDeriv(src, self.mesh, self.prob.timeMesh, u, v)
|
||||
# return data
|
||||
# else:
|
||||
# f = FieldsTDEM(self.mesh, self)
|
||||
# for src in self.srcList:
|
||||
# for rx in src.rxList:
|
||||
# Ptv = rx.evalDeriv(src, self.mesh, self.prob.timeMesh, u, v, adjoint=True)
|
||||
# Ptv = Ptv.reshape((-1, self.prob.timeMesh.nN), order='F')
|
||||
# if rx.projField not in f: # first time we are projecting
|
||||
# f[src, rx.projField, :] = Ptv
|
||||
# else: # there are already fields, so let's add to them!
|
||||
# f[src, rx.projField, :] += Ptv
|
||||
# return f
|
||||
|
||||
|
||||
|
||||
+70
-21
@@ -16,8 +16,6 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
|
||||
Problem.BaseTimeProblem.__init__(self, mesh, mapping=mapping, **kwargs)
|
||||
|
||||
|
||||
# _FieldsForward_pair = FieldsTDEM #: used for the forward calculation only
|
||||
|
||||
def fields(self, m):
|
||||
"""
|
||||
Solve the forward problem for the fields.
|
||||
@@ -49,12 +47,15 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
|
||||
if self.verbose: print 'Done'
|
||||
|
||||
rhs = self.getRHS(tInd, F)
|
||||
|
||||
if self.verbose: print ' Solving... (tInd = %d)'%tInd
|
||||
sol = Ainv * rhs
|
||||
if self.verbose: print ' Done...'
|
||||
|
||||
if sol.ndim == 1:
|
||||
sol.shape = (sol.size,1)
|
||||
F[:,self._fieldType+'Solution',tInd+1] = sol
|
||||
|
||||
Ainv.clean()
|
||||
return F
|
||||
|
||||
@@ -105,7 +106,60 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
|
||||
Ainv.clean()
|
||||
return Utils.mkvc(Jv)
|
||||
|
||||
|
||||
|
||||
def Jtvec(self, m, v, u=None):
|
||||
|
||||
if u is None:
|
||||
u = self.fields(m)
|
||||
|
||||
self.curModel = m
|
||||
ftype = self._fieldType + 'Solution' # the thing we solved for
|
||||
|
||||
# Ensure v is a data object.
|
||||
if not isinstance(v, self.dataPair):
|
||||
v = self.dataPair(self.survey, v)
|
||||
|
||||
# TODO: make this general
|
||||
if self._fieldType is 'b':
|
||||
dun_dmT_v = np.zeros((len(m), self.mesh.nF))
|
||||
|
||||
# df_dm_v = Fields_Derivs(self.mesh, self.survey)
|
||||
JTv = np.zeros(m.size)
|
||||
PT_v = Fields_Derivs(self.mesh, self.survey) #PT_v is a fields object
|
||||
|
||||
for src in self.survey.srcList:
|
||||
for rx in src.rxList:
|
||||
PT_v[src,'%sDeriv'%rx.projField, :] = rx.evalDeriv(src, self.mesh, self.timeMesh, v, adjoint = True) # All the fields for a given src, reciever.
|
||||
|
||||
|
||||
ATinv = None
|
||||
|
||||
for tInd, dt in enumerate(reversed(list(self.timeSteps))):
|
||||
if ATinv is not None and (tInd < self.nT and dt != self.timeSteps[tInd - 1]):# keep factors if dt is the same as previous step b/c A will be the same
|
||||
ATinv.clean()
|
||||
ATinv = None
|
||||
|
||||
if ATinv is None:
|
||||
A = self.getA(tInd)
|
||||
ATinv = self.Solver(A.T, **self.solverOpts)
|
||||
|
||||
for i, src in enumerate(self.survey.srcList):
|
||||
|
||||
u_src = u[src,ftype,tInd] # fields for this source at tInd
|
||||
|
||||
for rx in src.rxList:
|
||||
|
||||
df_duTFun = getattr(u, '_%sDeriv'%rx.projField, None)
|
||||
df_duT_v, df_dmT_v = df_duTFun(tInd, src, None, PT_v[src,'%sDeriv'%rx.projField,tInd-1], adjoint=True)
|
||||
|
||||
ATinv_df_duT_v = ATinv * df_duT_v
|
||||
rhsT_v = self.getJRHS(tInd, src, u_src, ATinv_df_duT_v, dun_dmT_v[:,i], adjoint = True)
|
||||
|
||||
JTv += rhsT_v + df_dmT_v
|
||||
|
||||
return Utils.mkvc(JTv)
|
||||
|
||||
|
||||
|
||||
def getJRHS(self, tInd, src, u, v, dbn_dm_v, adjoint = False):
|
||||
|
||||
@@ -115,9 +169,7 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
|
||||
b = - dA_dm + dRHS_dm
|
||||
|
||||
return b
|
||||
|
||||
def Jtvec(self, m, v, u=None):
|
||||
raise NotImplementedError
|
||||
|
||||
|
||||
def getSourceTerm(self, tInd):
|
||||
|
||||
@@ -245,7 +297,7 @@ class Problem_b(BaseTDEMProblem):
|
||||
if adjoint:
|
||||
if self._makeASymmetric is True:
|
||||
v = MfMui * v
|
||||
return MfMui.T * ( C * ( MeSigmaIDeriv.T * ( C.T * v ) ) )
|
||||
return MeSigmaIDeriv(C.T * ( MfMui * u )).T * ( C.T * v )
|
||||
|
||||
ADeriv = ( C * ( MeSigmaIDeriv(C.T * ( MfMui * u )) * v ) )
|
||||
if self._makeASymmetric is True:
|
||||
@@ -279,27 +331,24 @@ class Problem_b(BaseTDEMProblem):
|
||||
MfMui = self.MfMui
|
||||
|
||||
_, S_e = src.eval(tInd+1, self) # I think this is tInd+1 ?
|
||||
S_mDeriv_v, S_eDeriv_v = src.evalDeriv(self.times[tInd+1], self, v=v, adjoint=adjoint) # I think this is tInd+1 ?
|
||||
|
||||
# B_n = np.c_[[F[src,'b',tInd] for src in self.survey.srcList]].T
|
||||
# if B_n.shape[0] is not 1:
|
||||
# raise NotImplementedError('getRHS not implemented for this shape of B_n')
|
||||
S_mDeriv, S_eDeriv = src.evalDeriv(self.times[tInd+1], self, adjoint=adjoint) # I think this is tInd+1 ?
|
||||
|
||||
if adjoint:
|
||||
raise NotImplementedError
|
||||
if self._makeASymmetric is True:
|
||||
v = self.MfMui * v
|
||||
if isinstance(S_e, Utils.Zero):
|
||||
MeSigmaIDerivT_v = Utils.Zero()
|
||||
else:
|
||||
MeSigmaIDerivT_v = MeSigmaIDeriv(S_e).T * v
|
||||
RHSDeriv = MeSigmaIDerivT_v + S_eDeriv( MeSigmaI.T * ( C.T * v ) ) + S_mDeriv(v) + dbn_dm_v / dt #this will be given the transposed version
|
||||
return RHSDeriv
|
||||
|
||||
|
||||
if isinstance(S_e,Utils.Zero):
|
||||
if isinstance(S_e, Utils.Zero):
|
||||
MeSigmaIDeriv_v = Utils.Zero()
|
||||
else:
|
||||
MeSigmaIDeriv_v = MeSigmaIDeriv(S_e) * v
|
||||
|
||||
# if isinstance(S_eDeriv, Utils.Zero):
|
||||
# MeSigmaI_S_eDeriv_v = Utils.Zero()
|
||||
# else:
|
||||
# MeSigmaI_S_eDeriv_v = MeSigmaI * S_eDeriv(v)
|
||||
|
||||
RHSDeriv = (C * (MeSigmaIDeriv_v + MeSigmaI * S_eDeriv_v) + S_mDeriv_v) + dbn_dm_v / dt
|
||||
RHSDeriv = (C * (MeSigmaIDeriv_v + MeSigmaI * S_eDeriv(v) + S_mDeriv(v))) + dbn_dm_v / dt
|
||||
|
||||
if self._makeASymmetric is True:
|
||||
return self.MfMui.T * RHSDeriv
|
||||
|
||||
@@ -3,6 +3,10 @@ from SimPEG import *
|
||||
from SimPEG import EM
|
||||
|
||||
plotIt = False
|
||||
|
||||
testDeriv = False
|
||||
testAdjoint = True
|
||||
|
||||
tol = 1e-6
|
||||
|
||||
def setUp(rxcomp='bz'):
|
||||
@@ -70,132 +74,29 @@ class TDEM_bDerivTests(unittest.TestCase):
|
||||
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')
|
||||
if testDeriv:
|
||||
def test_Jvec_b_bx(self):
|
||||
self.JvecTest('bx')
|
||||
|
||||
def test_Jvec_b_bz(self):
|
||||
self.JvecTest('bz')
|
||||
def test_Jvec_b_bz(self):
|
||||
self.JvecTest('bz')
|
||||
|
||||
def test_Jvec_b_ey(self):
|
||||
self.JvecTest('ey')
|
||||
|
||||
def test_Jvec_b_ey(self):
|
||||
self.JvecTest('ey')
|
||||
|
||||
|
||||
# def test_projectAdjoint(self):
|
||||
# prb = self.prb
|
||||
# survey = prb.survey
|
||||
# mesh = self.mesh
|
||||
if testAdjoint:
|
||||
def test_adjointJvecVsJtvec(self):
|
||||
prb, m0, mesh = setUp()
|
||||
|
||||
# # 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.evalDeriv(None, v=f).tovec())
|
||||
# V2 = f.tovec().dot(survey.evalDeriv(None, v=d, adjoint=True).tovec())
|
||||
m = np.random.rand(prb.mapping.nP)
|
||||
d = np.random.rand(prb.survey.nD)
|
||||
|
||||
# self.assertTrue((V1-V2)/np.abs(V1) < tol)
|
||||
|
||||
# def test_adjointAhVsAht(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)
|
||||
|
||||
# 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._AhVec(sigma, f1).tovec())
|
||||
# V2 = f1.tovec().dot(prb._AhtVec(sigma, f2).tovec())
|
||||
# self.assertTrue(np.abs(V1-V2)/np.abs(V1) < tol)
|
||||
|
||||
# # def test_solveAhtVsAhtVec(self):
|
||||
# # prb = self.prb
|
||||
# # mesh = self.mesh
|
||||
# # sigma = np.random.rand(prb.mapping.nP)
|
||||
|
||||
# # f1 = EM.TDEM.FieldsTDEM(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 = prb.solveAht(sigma, f1)
|
||||
# # f3 = prb._AhtVec(sigma, f2)
|
||||
|
||||
# # if True:
|
||||
# # import matplotlib.pyplot as plt
|
||||
# # plt.plot(f3.tovec(),'b')
|
||||
# # plt.plot(f1.tovec(),'r')
|
||||
# # plt.show()
|
||||
# # V1 = np.linalg.norm(f3.tovec()-f1.tovec())
|
||||
# # V2 = np.linalg.norm(f1.tovec())
|
||||
# # print 'AhtVsAhtVec', V1, V2, f1.tovec()
|
||||
# # 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
|
||||
|
||||
# # 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)
|
||||
|
||||
# # V1 = f2.tovec().dot(prb.solveAh(sigma, f1).tovec())
|
||||
# # V2 = f1.tovec().dot(prb.solveAht(sigma, f2).tovec())
|
||||
# # print V1, V2
|
||||
# # self.assertLess(np.abs(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, 1)
|
||||
# u[:,'e',i] = np.random.rand(mesh.nE, 1)
|
||||
|
||||
# v = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
|
||||
# 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)
|
||||
|
||||
# V1 = m.dot(prb.Gtvec(sigma, v, u))
|
||||
# V2 = v.tovec().dot(prb.Gvec(sigma, m, u).tovec())
|
||||
# self.assertTrue(np.abs(V1-V2)/np.abs(V1) < tol)
|
||||
|
||||
# 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))
|
||||
# passed = np.abs(V1-V2)/np.abs(V1) < tol
|
||||
# print 'AdjointTest', V1, V2, passed
|
||||
# self.assertTrue(passed)
|
||||
V1 = d.dot(prb.Jvec(m0, m))
|
||||
V2 = m.dot(prb.Jtvec(m0, d))
|
||||
passed = np.abs(V1-V2)/np.abs(V1) < tol
|
||||
print 'AdjointTest', V1, V2, passed
|
||||
self.assertTrue(passed)
|
||||
|
||||
|
||||
|
||||
|
||||
Reference in New Issue
Block a user