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

tdem deriv runs but is first order at the moment

Browse Source
This commit is contained in:
Lindsey Heagy
2016-03-06 11:13:48 -08:00
parent 0bfc816ecc
commit 4f31e4e002
4 changed files with 144 additions and 84 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/EM/TDEM/FieldsTDEM.py
+17 -4
View File
@@ -31,7 +31,15 @@ class Fields(SimPEG.Problem.TimeFields):
knownFields = {}
dtype = float
class Fields_Derivs(Fields):
knownFields = {
'bDeriv': 'F',
'eDeriv': 'E',
'hDeriv': 'E',
'jDeriv': 'F'
}
class Fields_b(Fields):
"""Fancy Field Storage for a TDEM survey."""
knownFields = {'bSolution': 'F'}
@@ -48,12 +56,17 @@ class Fields_b(Fields):
def _b(self, bSolution, srcList, tInd):
return bSolution
def _bDeriv_u(self, src, du_dm_v, adjoint = False):
return Identity()*v
def _bDeriv_u(self, src, dun_dm_v, adjoint = False):
return Identity()*dun_dm_v
def _bDeriv_m(self, src, v, adjoint = False):
return Zero()
def _bDeriv(self, src, dun_dm_v, v, adjoint=False):
if adjoint is True:
raise NotImplementedError
return self._bDeriv_u(src, dun_dm_v) + self._bDeriv_m(src, v)
def _e(self, bSolution, srcList, tInd):
e = self.MeSigmaI * ( self.edgeCurl.T * ( self.MfMui * bSolution ) )
for i, src in enumerate(srcList):
@@ -61,7 +74,7 @@ class Fields_b(Fields):
e[:,i,tInd] = e[:,i,tInd] - self.MeSigmaI * S_e
return e
def _eDeriv_u(self, src, du_dm_v, adjoint = False):
def _eDeriv_u(self, src, dun_dm_v, adjoint = False):
raise NotImplementedError
def _eDeriv_m(self, src, v, adjoint = False):
SimPEG/EM/TDEM/SurveyTDEM.py
+10 -10
View File
@@ -60,13 +60,13 @@ class Rx(SimPEG.Survey.BaseTimeRx):
u_part = Utils.mkvc(u[src, self.projField, :])
return P*u_part
def evalDeriv(self, src, mesh, timeMesh, u, v, adjoint=False):
def evalDeriv(self, src, mesh, timeMesh, df_dm, adjoint=False):
P = self.getP(mesh, timeMesh)
if not adjoint:
return P * Utils.mkvc(v[src, self.projField, :])
return P * Utils.mkvc(df_dm[src, self.projField+'Deriv', :])
elif adjoint:
return P.T * v[src, self]
return P.T * df_dm[src, self]
####################################################
# Sources
@@ -121,7 +121,7 @@ class BaseSrc(SimPEG.Survey.BaseSrc):
rxPair = Rx
integrate = True
waveformPair = StepOffWaveform
waveformPair = BaseWaveform
@property
def waveform(self):
@@ -134,8 +134,8 @@ class BaseSrc(SimPEG.Survey.BaseSrc):
self._waveform = val
def __init__(self, rxList, waveform = None):
self.waveform = waveform
def __init__(self, rxList, waveform = None ):
self.waveform = waveform or StepOffWaveform()
SimPEG.Survey.BaseSrc.__init__(self, rxList)
@@ -162,11 +162,11 @@ class BaseSrc(SimPEG.Survey.BaseSrc):
def S_e(self, prob, time):
return Zero()
def S_mDeriv(self, prob, time):
raise NotImplementedError
def S_mDeriv(self, prob, time, v=None, adjoint=False):
return Zero()
def S_eDeriv(self, prob, time):
raise NotImplementedError
def S_eDeriv(self, prob, time, v=None, adjoint=False):
return Zero()
class MagDipole(BaseSrc):
SimPEG/EM/TDEM/TDEM.py
+102 -51
View File
@@ -33,8 +33,7 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
F = self.fieldsPair(self.mesh, self.survey)
# set initial fields
for i, src in enumerate(self.survey.srcList):
F[src,self._fieldType+'Solution',0] = src.bInitial(self) # TODO: this will only work for b formulation
F[:,self._fieldType+'Solution',0] = self.getInitialFields()
# timestep to solve forward
Ainv = None
@@ -67,50 +66,64 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
if u is None:
u = self.fields(m)
ftype = self._fieldType + 'Solution' # the thing we solved for
self.curModel = m
# def getb(tInd, src, u, v):
# dA_dm = self.getADeriv(tInd, u, v, adjoint=False)
# dRHS_dm = self.getRHSDeriv(tInd, src, v, adjoint=False)
Jv = self.dataPair(self.survey)
# b = - dA_dm + dRHS_dm
dun_dm_v = self.getInitialFieldsDeriv(v) # can over-write this at each timestep
df_dm = Fields_Derivs(self.mesh, self.survey)
# np.zeros((dun_dm_v.shape[0], self.nT, self.survey.nSrc))
Jv = self.dataPair(self.survey)
print Jv.shape
raise NotImplementedError
Asub = Utils.speye(len(u))
if self._makeASymmetric:
Asub = self.MfMui.T * Asub #TODO: this will only work for E-B formulation
Adiag, Asub = None, None
Adiaginv = None
if self._makeASymmetric:
Asub = self.MfMui.T * Asub #TODO: this will only work for E-B formulation
for tInd, dT in enumerate(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
for tInd, dt in enumerate(self.timeSteps):
if Adiag 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
Adiag, Asub = None, None
if Adiaginv is None:
Adiag = self.getA(tInd)
if Adiag is None:
Adiag, Asub = self.getJdiags(tInd, adjoint = False)
Adiaginv = self.Solver(Adiag)
for src in self.survey.srcList:
# just construction of RHS
dRHS_dm_v = self.getRHSDeriv(tInd, src, v, adjoint=False)
if tInd == 0:
dRHS_dm_v = dRHS_dm_v + src.bInitialDeriv(self, v, adjoint=False)
# else:
# db_dm = Jv[]
for i, src in enumerate(self.survey.srcList):
# compute next du_dm_v for next timestep
u_src = u[src,ftype,tInd]
rhs_v = self.getJRHS(tInd, src, u_src, v, dun_dm_v[:,i])
for rx in src.rxList:
df_dmFun = getattr(u, '_%sDeriv'%rx.projField, None)
df_dm[src, '%sDeriv'%rx.projField , tInd] = df_dmFun(src, dun_dm_v, v)
# Adiag, Asub, b = self.diagsJ(0, m, u, v, adjoint=False)
# AdiagInv
# Jv[]
# over-write with this time-steps (if not on last timestep)
if tInd != len(self.timeSteps):
dun_dm_v[:,i] = Adiaginv * rhs_v
for src in self.survey.srcList:
for rx in src.rxList:
Jv[src,rx] = rx.evalDeriv(src, self.mesh, self.timeMesh, df_dm)
Adiaginv.clean()
return Utils.mkvc(Jv)
# if tInd == 0:
# dbn_dm_v = get
# else:
# rhs_v = self.getJRHS(tInd, m, u, v, dbn_dm_v)
def getJRHS(self, tInd, src, u, v, dbn_dm_v, adjoint = False):
dA_dm = self.getADeriv(tInd, u, v, adjoint)
dRHS_dm = self.getRHSDeriv(tInd, src, v, dbn_dm_v, adjoint)
b = - dA_dm + dRHS_dm
return b
def Jtvec(self, m, v, u=None):
raise NotImplementedError
@@ -133,6 +146,33 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
return S_m, S_e
def getInitialFields(self):
Srcs = self.survey.srcList
if self._fieldType is 'b' or self._fieldType is 'j':
ifields = np.zeros((self.mesh.nF, len(Srcs)))
elif self._fieldType is 'e' or self._fieldType is 'h':
ifields = np.zeros((self.mesh.nE, len(Srcs)))
for i,src in enumerate(Srcs):
ifields[:,i] = ifields[:,i] + getattr(src, '%sInitial'%self._fieldType, None)(self)
return ifields
def getInitialFieldsDeriv(self, v, adjoint=False):
Srcs = self.survey.srcList
if self._fieldType is 'b' or self._fieldType is 'j':
ifieldsDeriv = np.zeros((self.mesh.nF, len(Srcs)))
elif self._fieldType is 'e' or self._fieldType is 'h':
ifieldsDeriv = np.zeros((self.mesh.nE, len(Srcs)))
for i,src in enumerate(Srcs):
ifieldsDeriv[:,i] = ifieldsDeriv[:,i] + getattr(src, '%sInitialDeriv'%self._fieldType, None)(self,v,adjoint)
return ifieldsDeriv
##########################################################################################
@@ -208,7 +248,7 @@ class Problem_b(BaseTDEMProblem):
def getADeriv(self, tInd, u, v, adjoint=False):
C = self.mesh.edgeCurl
MeSigmaI = lambda u: self.MeSigmaIDeriv(u)
MeSigmaIDeriv = lambda u: self.MeSigmaIDeriv(u)
MfMui = self.MfMui
I = Utils.speye(self.mesh.nF)
@@ -219,7 +259,7 @@ class Problem_b(BaseTDEMProblem):
ADeriv = ( C * ( MeSigmaIDeriv(C.T * ( MfMui * u )) * v ) )
if self._makeASymmetric is True:
return MeMui.T * ADeriv
return MfMui.T * ADeriv
return ADeriv
@@ -249,7 +289,7 @@ class Problem_b(BaseTDEMProblem):
MfMui = self.MfMui
_, S_e = src.eval(tInd+1, self) # I think this is tInd+1 ?
S_mDeriv, S_eDeriv = src.evalDeriv(tInd+1, self, adjoint=adjoint) # 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:
@@ -257,38 +297,49 @@ class Problem_b(BaseTDEMProblem):
if adjoint:
raise NotImplementedError
return (C * (MeSigmaIDeriv (S_e) * v + MeSigmaI * S_eDeriv) + S_mDeriv) + bdn_dm_v / dt
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
if self._makeASymmetric is True:
return self.MfMui.T * RHSDeriv
return RHSDeriv
@Utils.timeIt
def diagsJ(self, tInd, m, u, v, adjoint = False):
def getJdiags(self, tInd, adjoint = False):
# The matrix that we are computing has the form:
#
# - - - - - -
# - - - - - -
# | Adiag | | uderiv1 | | b1 |
# | Asub Adiag | | uderiv2 | | b2 |
# | Asub Adiag | | uderiv3 | = | b3 |
# | ... ... | | ... | | .. |
# | Asub Adiag | | uderivn | | bn |
# - - - - - -
# - - - - - -
if adjoint:
raise NotImplementedError
dt = self.timeSteps[tInd]
Adiag = self.getA(tInd)
Asub = Utils.speye(len(u))
Asub = - 1./dt * Utils.speye(self.mesh.nF)
if self._makeASymmetric:
Asub = self.MfMui.T * Asub
dA_dm = self.getADeriv(tInd, u, v, adjoint)
dRHS_dm = self.getRHSDeriv(self, tInd, src, v, dbn_dm_v, adjoint)
b = - dA_dm + dRHS_dm
return Adiag, Asub, b
return Adiag, Asub
tests/em/tdem/test_TDEM_b_DerivAdjoint.py
+15 -19
View File
@@ -22,12 +22,12 @@ class TDEM_bDerivTests(unittest.TestCase):
mapping = Maps.ExpMap(mesh) * Maps.SurjectVertical1D(mesh) * activeMap
rxOffset = 40.
rx = EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20), 'bz')
src = EM.TDEM.SrcTDEM_VMD_MVP([rx], loc=np.array([0., 0., 0.]))
rx = EM.TDEM.Rx(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20), 'bz')
src = EM.TDEM.SurveyTDEM.MagDipole([rx], loc=np.array([0., 0., 0.]))
survey = EM.TDEM.SurveyTDEM([src])
survey = EM.TDEM.Survey([src])
self.prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping)
self.prb = EM.TDEM.Problem_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)]
@@ -45,10 +45,6 @@ class TDEM_bDerivTests(unittest.TestCase):
self.prb.pair(survey)
self.mesh = mesh
def test_ADeriv(self):
prb = self.prb
m = self.m
# def test_AhVec(self):
# """
# Test that fields and AhVec produce consistent results
@@ -178,21 +174,21 @@ class TDEM_bDerivTests(unittest.TestCase):
# print 'test_Deriv_dUdM'
# Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20)
# 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
m = self.m
# # 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_m = 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'
# Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=4, eps=1e-20)
derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(m, mx)]
print '\n'
print 'test_Deriv_J'
Tests.checkDerivative(derChk, m, plotIt=False, dx=d_m, num=4, eps=1e-20)
# def test_projectAdjoint(self):
# prb = self.prb