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

first shot through of passing Jtvec for TDEM problem (code will need to be cleaned up, but it passes!)

Browse Source
This commit is contained in:
Lindsey Heagy
2016-03-12 15:13:17 -08:00
parent fb66acea11
commit a1ecef0709
3 changed files with 110 additions and 47 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/EM/TDEM/SurveyTDEM.py
+13 -11
View File
@@ -1,8 +1,8 @@
import SimPEG
from SimPEG import np, Utils
from SimPEG.Utils import Zero, Identity
from scipy.constants import mu_0
from SimPEG.EM.Utils import *
from SimPEG.Utils import Zero, Identity
from scipy.constants import mu_0
from SimPEG.EM.Utils import *
####################################################
# Receivers
@@ -60,13 +60,15 @@ class Rx(SimPEG.Survey.BaseTimeRx):
u_part = Utils.mkvc(u[src, self.projField, :])
return P*u_part
def evalDeriv(self, src, mesh, timeMesh, df_dm, adjoint=False):
def evalDeriv(self, src, mesh, timeMesh, v, adjoint=False):
P = self.getP(mesh, timeMesh)
if not adjoint:
return P * Utils.mkvc(df_dm[src, self.projField+'Deriv', :])
return P * v #Utils.mkvc(v[src, self.projField+'Deriv', :])
elif adjoint:
return P.T * df_dm[src, self]
# dP_dF_T = P.T * v #[src, self]
# newshape = (len(dP_dF_T)/timeMesh.nN, timeMesh.nN )
return P.T * v #np.reshape(dP_dF_T, newshape, order='F')
####################################################
# Sources
@@ -83,10 +85,10 @@ class BaseWaveform(object):
), "Waveform object must be an instance of a %s BaseWaveform class."%(pair.__name__)
def eval(self, time):
raise NotImplementedError
raise NotImplementedError
def evalDeriv(self, time):
raise NotImplementedError # needed for E-formulation
raise NotImplementedError # needed for E-formulation
class StepOffWaveform(BaseWaveform):
@@ -136,7 +138,7 @@ class BaseSrc(SimPEG.Survey.BaseSrc):
def __init__(self, rxList, waveform = None ):
self.waveform = waveform or StepOffWaveform()
SimPEG.Survey.BaseSrc.__init__(self, rxList)
SimPEG.Survey.BaseSrc.__init__(self, rxList)
def bInitial(self, prob):
@@ -148,7 +150,7 @@ class BaseSrc(SimPEG.Survey.BaseSrc):
def eval(self, prob, time):
S_m = self.S_m(prob, time)
S_e = self.S_e(prob, time)
return S_m, S_e
return S_m, S_e
def evalDeriv(self, prob, time, v=None, adjoint=False):
if v is not None:
@@ -170,7 +172,7 @@ class BaseSrc(SimPEG.Survey.BaseSrc):
class MagDipole(BaseSrc):
def __init__(self, rxList, waveform=None, loc=None, orientation='Z', moment=1., mu=mu_0):
def __init__(self, rxList, waveform=None, loc=None, orientation='Z', moment=1., mu=mu_0):
self.loc = loc
self.orientation = orientation
SimPEG/EM/TDEM/TDEM.py
+69 -25
View File
@@ -63,6 +63,17 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
def Jvec(self, m, v, u=None):
"""
Jvec computes the sensitivity times a vector
.. math::
\mathbf{J} \mathbf{v} = \\frac{d\mathbf{P}}{d\mathbf{F}} \left( \\frac{d\mathbf{F}}{d\mathbf{u}} \\frac{d\mathbf{u}}{d\mathbf{m}} + \\frac{\partial\mathbf{F}}{\partial\mathbf{m}} \\right) \mathbf{v}
where
.. math::
\mathbf{A} \\frac{d\mathbf{u}}{d\mathbf{m}} + \\frac{d\mathbf{A}(\mathbf{u})}{d\mathbf{m}} = \\frac{d \mathbf{RHS}}{d \mathbf{m}}
"""
if u is None:
u = self.fields(m)
@@ -81,7 +92,7 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
Adiaginv = None
for tInd, dt in zip(range(self.nT+1), self.timeSteps):
for tInd, dt in zip(range(self.nT), self.timeSteps):
if Adiaginv is not None and (tInd > 0 and dt != self.timeSteps[tInd - 1]):# keep factors if dt is the same as previous step b/c A will be the same
Adiaginv.clean()
Adiaginv = None
@@ -115,7 +126,7 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
for src in self.survey.srcList:
for rx in src.rxList:
Jv[src,rx] = rx.evalDeriv(src, self.mesh, self.timeMesh, df_dm_v)
Jv[src,rx] = rx.evalDeriv(src, self.mesh, self.timeMesh, Utils.mkvc(df_dm_v[src,'%sDeriv'%rx.projField,:]))
Adiaginv.clean()
return Utils.mkvc(Jv)
@@ -140,61 +151,94 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
PT_v = Fields_Derivs(self.mesh, self.survey) #PT_v is a fields object
# TODO: This will only work for b formulation right now b/c of the mesh.nF
df_duT_v = np.zeros((self.mesh.nF,self.nT+1))
ATinv_df_duT_v = np.zeros((self.mesh.nF,self.nT))
ATinv_df_duT_v = np.zeros((self.mesh.nF,self.nT+1))
JTv = np.zeros(m.shape)
# TODO : this is pretty ugly
# Loop over sources and receivers to create a fields object: PT_v
for src in self.survey.srcList:
# for rx in src.rxList:
# initialize empty fields derivs
for projField in set([rx.projField for rx in src.rxList]):
PT_v[src,'%sDeriv'%projField, :] = rx.evalDeriv(src, self.mesh, self.timeMesh, v, adjoint = True) # All the fields for a given src, reciever.
PT_v[src,'%sDeriv'%projField, :] = np.zeros_like(u[src, '%s'%projField, : ])
# loop over recievers and sum contributions to fields object
for rx in src.rxList:
# for projField in set([rx.projField for rx in src.rxList]):
curPT_v = rx.evalDeriv(src, self.mesh, self.timeMesh, Utils.mkvc(v[src,rx]), adjoint=True)
PT_v[src,'%sDeriv'%rx.projField, :] += np.reshape(curPT_v,(len(curPT_v)/self.timeMesh.nN, self.timeMesh.nN), order='F') # All the fields for a given src, reciever.
# print np.linalg.norm(PT_v[src,'bDeriv',:])
# for src in self.survey.srcList:
# initialize empty fields derivs
for projField in set([rx.projField for rx in src.rxList]):
df_duTFun = getattr(u, '_%sDeriv'%projField, None)
# TODO: don't need to recompute df_dmT_v every time... only need it once
df_duT_v_cur, df_dmT_v = df_duTFun(None, src, None, PT_v[src,'%sDeriv'%projField,:], adjoint=True) # this seems odd
df_duT_v_cur, df_dmT_v = df_duTFun(None, src, None, PT_v[src,'%sDeriv'%projField,:], adjoint=True)
JTv = JTv + df_dmT_v
df_duT_v = df_duT_v + df_duT_v_cur
AdiagTinv = None
JTv = df_dmT_v
# print np.linalg.norm(df_duT_v)
for tInd in reversed(range(self.nT)): #enumerate(reversed(list(self.timeSteps))):
if AdiagTinv is not None and (tInd < self.nT and self.timeSteps[tInd] != self.timeSteps[tInd + 1]):# keep factors if dt is the same as previous step b/c A will be the same
AdiagTinv = None
# for tInd in reversed(range(self.nT)): #enumerate(reversed(list(self.timeSteps))):
for tInd in reversed(range(self.nT)) : # reversed(self.timeSteps)):
if AdiagTinv is not None: # and (tInd <= self.nT and dt != self.timeSteps[tInd]):
# (tInd < self.nT and self.timeSteps[tInd] != self.timeSteps[tInd + 1]):# keep factors if dt is the same as previous step b/c A will be the same
AdiagTinv.clean()
AdiagTinv = None
if AdiagTinv is None:
Adiag = self.getAdiag(tInd)
AdiagTinv = self.Solver(Adiag.T, **self.solverOpts)
Asubdiag = self.getAsubdiag(tInd)
# solve against df_duT_v
if tInd < self.nT:
# print Utils.mkvc(AdiagTinv * df_duT_v[:,tInd],2).shape, ATinv_df_duT_v[:,tInd].shape
ATinv_df_duT_v[:,tInd] = AdiagTinv * df_duT_v[:,tInd]
if tInd >= self.nT-1:
ATinv_df_duT_v[:,tInd+1] = AdiagTinv * df_duT_v[:,tInd+1]
else:
ATinv_df_duT_v[:,tInd] = AdiagTinv * (df_duT_v[:,tInd+1] - Asubdiag.T * df_duT_v[:,tInd])
Asubdiag = self.getAsubdiag(tInd+1)
ATinv_df_duT_v[:,tInd+1] = AdiagTinv * (df_duT_v[:,tInd+1] - Asubdiag.T * ATinv_df_duT_v[:,tInd+2])
un_src = u[src,ftype,tInd]
dAT_dm_v = self.getAdiagDeriv(tInd, un_src, ATinv_df_duT_v[:,tInd], adjoint=True) # cell centered on time mesh
# for src in self.survey.srcList:
un_src = u[src,ftype,tInd+1]
dAT_dm_v = self.getAdiagDeriv(None, un_src, ATinv_df_duT_v[:,tInd+1], adjoint=True) # cell centered on time mesh
dRHST_dm_v = self.getRHSDeriv(tInd, src, ATinv_df_duT_v[:,tInd], adjoint=True) # on nodes of time mesh
dRHST_dm_v = self.getRHSDeriv(tInd+1, src, ATinv_df_duT_v[:,tInd+1], adjoint=True) # on nodes of time mesh
# dAsubdiag_dm_v = 0
JTv = JTv + (-dAT_dm_v + dRHST_dm_v)
# JTv = JTv +
# tInd = 0
# un_src = u[src,ftype,tInd]
# # dAT_dm_v = self.getAdiagDeriv(None, un_src, self.getInitialFieldsDeriv(), adjoint=True)
# Asubdiag = self.getAsubdiag(tInd)
# ATinv_df_duT_v[:,tInd] = AdiagTinv * (- Asubdiag.T * df_duT_v[:,tInd])
# # - self.getAsubdiag(tInd).T * df_duT_v[:,tInd]
# dAT_dm_v = self.getAdiagDeriv(None, un_src, ATinv_df_duT_v[:,tInd], adjoint=True) # cell centered on time mesh
# dRHST_dm_v = self.getRHSDeriv(tInd, src, ATinv_df_duT_v[:,tInd], adjoint=True)
# JTv = JTv + (- dAT_dm_v + dRHST_dm_v)
return JTv
# for i, src in enumerate(self.survey.srcList):
# un_src = u[src,ftype,tInd+1] # fields for this source at tInd
@@ -402,7 +446,7 @@ class Problem_b(BaseTDEMProblem):
MfMui = self.MfMui
_, S_e = src.eval(tInd, self)
S_mDeriv, S_eDeriv = src.evalDeriv(self.times[tInd], self, adjoint=adjoint) # I think this is tInd+1 ?
S_mDeriv, S_eDeriv = src.evalDeriv(self.times[tInd], self, adjoint=adjoint)
if adjoint:
if self._makeASymmetric is True:
tests/em/tdem/test_TDEM_b_DerivAdjoint.py
+28 -11
View File
@@ -4,7 +4,7 @@ from SimPEG import EM
plotIt = False
testDeriv = False
testDeriv = True
testAdjoint = True
tol = 1e-6
@@ -30,6 +30,7 @@ def setUp(rxcomp='bz'):
prb = EM.TDEM.Problem_b(mesh, mapping=mapping)
prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)]
# prb.timeSteps = [(1e-05, 10), (1e-05, 50), (1e-05, 50) ] #, (2.5e-4, 10)]
try:
from pymatsolver import MumpsSolver
@@ -48,7 +49,7 @@ class TDEM_bDerivTests(unittest.TestCase):
def test_ADeriv(self):
def test_Deriv_Pieces(self):
prb, m0, mesh = setUp()
tInd = 0
@@ -63,19 +64,34 @@ class TDEM_bDerivTests(unittest.TestCase):
return Av, ADeriv_dm
# def A_adjointTest():
def A_adjointTest():
print '\n Testing A_adjoint'
m = np.random.rand(prb.mapping.nP)
v = np.random.rand(prb.mesh.nF)
u = np.random.rand(prb.mesh.nF)
prb.curModel = m0
tInd = 0 # not actually used
V1 = v.dot(prb.getAdiagDeriv(tInd, u, m))
V2 = m.dot(prb.getAdiagDeriv(tInd, u, v, adjoint=True))
passed = np.abs(V1-V2)/np.abs(V1) < tol
print 'AdjointTest', V1, V2, passed
self.assertTrue(passed)
# def P_adjointTest():
# print '\n Testing P_adjoint'
# m = np.random.rand(prb.mapping.nP)
# d = np.random.rand(prb.survey.nD)
# v = np.random.rand(prb.mesh.nF)
# d = np.random.rand(prb.survey.nD)
# prb.curModel = m0
# 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)
# for src in prb.survey.srcList:
# for rx in src.rxList:
# Tests.checkDerivative(AderivTest, m0, plotIt=False, num=4, eps=1e-20)
print '\n Testing ADeriv'
Tests.checkDerivative(AderivTest, m0, plotIt=False, num=4, eps=1e-20)
A_adjointTest()
@@ -100,10 +116,11 @@ class TDEM_bDerivTests(unittest.TestCase):
if testAdjoint:
def test_adjointJvecVsJtvec(self):
print '\n Adjoint Testing Jvec, Jtvec'
prb, m0, mesh = setUp()
m = np.random.rand(prb.mapping.nP)
d = np.random.rand(prb.survey.nD)
d = np.random.randn(prb.survey.nD)
V1 = d.dot(prb.Jvec(m0, m))
V2 = m.dot(prb.Jtvec(m0, d))