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

Update current changes ... about u -> f

Browse Source
This commit is contained in:
seogi_macbook
2016-06-01 23:45:52 -07:00
parent 8e6536b2ed
commit 6f8f065643
1 changed files with 19 additions and 19 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/EM/TDEM/TDEM.py
+19 -19
View File
@@ -62,7 +62,7 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
return F
def Jvec(self, m, v, u=None):
def Jvec(self, m, v, f=None):
"""
Jvec computes the sensitivity times a vector
@@ -75,8 +75,8 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
\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)
if f is None:
f = self.fields(m)
ftype = self._fieldType + 'Solution' # the thing we solved for
self.curModel = m
@@ -115,17 +115,17 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
# here, we are lagging by a timestep, so filling in as we go
for projField in set([rx.projField for rx in src.rxList]):
df_dmFun = getattr(u, '_%sDeriv'%projField, None)
df_dmFun = getattr(f, '_%sDeriv'%projField, None)
# df_dm_v is dense, but we only need the times at (rx.P.T * ones > 0)
# This should be called rx.footprint
df_dm_v[src, '%sDeriv'%projField , tInd] = df_dmFun(tInd, src, dun_dm_v[:,i], v)
un_src = u[src,ftype,tInd+1]
un_src = f[src,ftype,tInd+1]
dA_dm_v = self.getAdiagDeriv(tInd, un_src, v) # cell centered on time mesh
dRHS_dm_v = self.getRHSDeriv(tInd+1, src, v) # on nodes of time mesh
dAsubdiag_dm_v = self.getAsubdiagDeriv(tInd, u[src,ftype,tInd], v)
dAsubdiag_dm_v = self.getAsubdiagDeriv(tInd, f[src,ftype,tInd], v)
JRHS = dRHS_dm_v - dAsubdiag_dm_v - dA_dm_v
@@ -141,7 +141,7 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
return Utils.mkvc(Jv)
def Jtvec(self, m, v, u=None):
def Jtvec(self, m, v, f=None):
"""
Jvec computes the adjoint of the sensitivity times a vector
@@ -155,8 +155,8 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
\\frac{d\mathbf{u}}{d\mathbf{m}} ^\\top \mathbf{A}^\\top + \\frac{d\mathbf{A}(\mathbf{u})}{d\mathbf{m}} ^ \\top = \\frac{d \mathbf{RHS}}{d \mathbf{m}} ^ \\top
"""
if u is None:
u = self.fields(m)
if f is None:
f = self.fields(m)
self.curModel = m
ftype = self._fieldType + 'Solution' # the thing we solved for
@@ -166,7 +166,7 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
v = self.dataPair(self.survey, v)
df_duT_v = Fields_Derivs(self.mesh, self.survey)
ATinv_df_duT_v = np.zeros((len(self.survey.srcList), len(u[self.survey.srcList[0],ftype,0]))) # same size as fields at a single timestep
ATinv_df_duT_v = np.zeros((len(self.survey.srcList), len(f[self.survey.srcList[0],ftype,0]))) # same size as fields at a single timestep
JTv = np.zeros(m.shape)
@@ -174,20 +174,20 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
for src in self.survey.srcList:
PT_v = Fields_Derivs(self.mesh, self.survey) # initialize storage for PT_v (don't need to preserve over sources)
# initialize size
df_duT_v[src, '%sDeriv'%self._fieldType, :] = np.zeros_like(u[src, self._fieldType, :])
df_duT_v[src, '%sDeriv'%self._fieldType, :] = np.zeros_like(f[src, self._fieldType, :])
for rx in src.rxList:
PT_v[src,'%sDeriv'%rx.projField,:] = rx.evalDeriv(src, self.mesh, self.timeMesh, Utils.mkvc(v[src,rx]), adjoint=True) # this is +=
# PT_v = np.reshape(curPT_v,(len(curPT_v)/self.timeMesh.nN, self.timeMesh.nN), order='F')
df_duTFun = getattr(u, '_%sDeriv'%rx.projField, None)
df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None)
for tInd in range(self.nT+1):
cur = df_duTFun(tInd, src, None, Utils.mkvc(PT_v[src,'%sDeriv'%rx.projField,tInd]), adjoint=True)
df_duT_v[src, '%sDeriv'%self._fieldType, tInd] = df_duT_v[src, '%sDeriv'%self._fieldType, tInd] + Utils.mkvc(cur[0],2)
JTv = JTv + cur[1]
JTv = cur[1] + JTv
del PT_v # no longer need this
@@ -227,10 +227,10 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
# (Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1]) - Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:]))
if tInd < self.nT - 1:
dAsubdiagT_dm_v = self.getAsubdiagDeriv(tInd+1, u[src,ftype,tInd+1], ATinv_df_duT_v[isrc,:], adjoint = True)
dAsubdiagT_dm_v = self.getAsubdiagDeriv(tInd+1, f[src,ftype,tInd+1], ATinv_df_duT_v[isrc,:], adjoint = True)
if tInd > -1:
un_src = u[src,ftype,tInd+1]
un_src = f[src,ftype,tInd+1]
dAT_dm_v = self.getAdiagDeriv(tInd, un_src, ATinv_df_duT_v[isrc,:], adjoint=True) # cell centered on time mesh
dRHST_dm_v = self.getRHSDeriv(tInd+1, src, ATinv_df_duT_v[isrc,:], adjoint=True) # on nodes of time mesh
@@ -240,7 +240,7 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
# print np.linalg.norm(self.getInitialFieldsDeriv(src, df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1], adjoint=True))
# print np.linalg.norm(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1])
# vec = - Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:]) + Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1])
# dAsubdiagT_dm_v = self.getAsubdiagDeriv(tInd+1, u[src,ftype,tInd+1], Utils.mkvc(ATinv_df_duT_v[isrc,:]), adjoint = True)
# dAsubdiagT_dm_v = self.getAsubdiagDeriv(tInd+1, f[src,ftype,tInd+1], Utils.mkvc(ATinv_df_duT_v[isrc,:]), adjoint = True)
dRHST_dm_v = Utils.mkvc(self.getInitialFieldsDeriv(src, Utils.mkvc(ATinv_df_duT_v[isrc,:]) , adjoint=True))
JTv = JTv + Utils.mkvc( -dAsubdiagT_dm_v + dRHST_dm_v) #
@@ -253,7 +253,7 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
# JTv = JTv + Utils.mkvc(dRHST_dm_v0 + dRHST_dm_v1)
# print 'here'
# inFields = self.getInitialFieldsDeriv(u[src,ftype,tInd+1], Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1]), adjoint=True)
# inFields = self.getInitialFieldsDeriv(f[src,ftype,tInd+1], Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1]), adjoint=True)
# # - Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:]), adjoint=True)
# print inFields.shape
# JTv = JTv + inFields
@@ -268,11 +268,11 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
# for src in self.survey.srcList:
# for projField in set(rx.projField):
# v = AdiagTinv * (Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,0]) - Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:]))
# JTv = JTv - Utils.mkvc(self.getAdiagDeriv(0, u[src, ftype, tInd], v, adjoint = True))
# JTv = JTv - Utils.mkvc(self.getAdiagDeriv(0, f[src, ftype, tInd], v, adjoint = True))
# # JTv = JTv + self.getInitialFieldsDeriv(Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,0] - Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:])), adjoint=True)
return Utils.mkvc(JTv)
return Utils.mkvc(JTv).astype(float)