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

sketch of derivs

Browse Source
This commit is contained in:
Lindsey Heagy
2016-02-24 15:24:56 -08:00
parent 617241ad4e
commit 8bd027c2b2
5 changed files with 101 additions and 18 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/EM/TDEM/SurveyTDEM.py
+4 -1
View File
@@ -118,7 +118,10 @@ class BaseSrc(SimPEG.Survey.BaseSrc):
SimPEG.Survey.BaseSrc.__init__(self, rxList)
def bInitial(self):
def bInitial(self, prob):
return Zero()
def bInitialDeriv(self, prob, v=None, adjoint=False):
return Zero()
def eval(self, prob, time):
SimPEG/EM/TDEM/TDEM.py
+87 -11
View File
@@ -34,7 +34,7 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
# set initial fields
for i, src in enumerate(self.survey.srcList):
F[src,'bSolution',0] = src.bInitial(self)
F[src,self._fieldType+'Solution',0] = src.bInitial(self) # TODO: this will only work for b formulation
# timestep to solve forward
Ainv = None
@@ -59,8 +59,84 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
Ainv.clean()
return F
# @Utils.timeIt
# def diagsJ(self, tInd, m, u, v, 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
# Adiag = self.getA(tInd)
# Asub = Utils.speye(len(u))
# if self._makeASymmetric:
# Asub = self.MfMui.T * Asub #TODO: this will only work for E-B formulation
# dA_dm = self.getADeriv(tInd, u, v, adjoint)
# dRHS_dm = self.getRHSDeriv(self, tInd, src, v, adjoint)
# b = - dA_dm + dRHS_dm
# return Adiag, Asub, b
def Jvec(self, m, v, u=None):
# raise NotImplementedError
if u is None:
u = self.fields(m)
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)
# b = - dA_dm + dRHS_dm
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
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
Adiaginv.clean()
Adiaginv = None
if Adiaginv is None:
Adiag = self.getA(tInd)
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[]
# Adiag, Asub, b = self.diagsJ(0, m, u, v, adjoint=False)
# AdiagInv
# Jv[]
def Jtvec(self, m, v, u=None):
raise NotImplementedError
@@ -156,7 +232,7 @@ class Problem_b(BaseTDEMProblem):
return MfMui.T * A
return A
def getADeriv(self, freq, u, v, adjoint=False):
def getADeriv(self, tInd, u, v, adjoint=False):
dt = self.timeSteps[tInd]
C = self.mesh.edgeCurl
MeSigmaI = self.MeSigmaIDeriv
@@ -191,9 +267,7 @@ class Problem_b(BaseTDEMProblem):
return MfMui.T * rhs
return rhs
def getRHSDeriv(self, tInd, F, src, v, adjoint=False):
raise NotImplementedError
def getRHSDeriv(self, tInd, src, v, adjoint=False):
dt = self.timeSteps[tInd]
C = self.mesh.edgeCurl
@@ -201,14 +275,16 @@ class Problem_b(BaseTDEMProblem):
MeSigmaIDeriv = self.MeSigmaIDeriv
MfMui = self.MfMui
S_m, S_e = src.eval(tInd+1, self) # I think this is tInd+1 ?
S_m, S_e = src.evalDeriv(tInd+1, self, adjoint=adjoint) # I think this is tInd+1 ?
_, 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 ?
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')
# 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')
# return B_n + dt * (C * (MeSigmaIDeriv * S_e) + S_m)
if adjoint:
raise NotImplementedError
return dt * (C * (MeSigmaIDeriv * S_e + MeSigmaI * S_eDeriv) + S_mDeriv)
tests/em/tdem/test_TDEM_b_DerivAdjoint.py
+5 -1
View File
@@ -40,11 +40,15 @@ class TDEM_bDerivTests(unittest.TestCase):
self.sigma = np.ones(mesh.nCz)*1e-8
self.sigma[mesh.vectorCCz<0] = 1e-1
self.sigma = np.log(self.sigma[active])
self.m = np.log(self.sigma[active])
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
tests/em/tdem/test_TDEM_b_MultiSrc_DerivAdjoint.py
+2 -2
View File
@@ -149,5 +149,5 @@ class TDEM_bDerivTests(unittest.TestCase):
if __name__ == '__main__':
unittest.main()
# if __name__ == '__main__':
# unittest.main()
tests/em/tdem/test_TDEM_combos.py
+3 -3
View File
@@ -60,7 +60,7 @@ def dotestAdjoint(prb, mesh, sigma):
print 'AdjointTest', V1, V2
return np.abs(V1-V2)/np.abs(V1), 1e-6
class TDEM_bDerivTests(unittest.TestCase):
# class TDEM_bDerivTests(unittest.TestCase):
# def test_Jvec_bx(self): self.assertTrue(dotestJvec(*getProb(rxTypes='bx')))
# def test_Adjoint_bx(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='bx')))
@@ -90,5 +90,5 @@ class TDEM_bDerivTests(unittest.TestCase):
# def test_Adjoint_eybzdbxdt(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='ey,bz,dbxdt')))
if __name__ == '__main__':
unittest.main()
# if __name__ == '__main__':
# unittest.main()