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

light notation cleanup in Jvec, testing ADeriv --> passes, Jvec is still first order

Browse Source
This commit is contained in:
Lindsey Heagy
2016-03-06 15:10:13 -08:00
parent 4f31e4e002
commit 664adb04ac
2 changed files with 64 additions and 25 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/EM/TDEM/TDEM.py
+10 -16
View File
@@ -61,8 +61,6 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
def Jvec(self, m, v, u=None):
# raise NotImplementedError
if u is None:
u = self.fields(m)
@@ -70,22 +68,19 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
self.curModel = m
Jv = self.dataPair(self.survey)
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))
Adiag, Asub = None, None
Ainv = None
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()
Adiag, Asub = None, None
if Ainv 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
Ainv.clean()
Ainv = None
if Adiag is None:
Adiag, Asub = self.getJdiags(tInd, adjoint = False)
Adiaginv = self.Solver(Adiag)
if Ainv is None:
A = self.getA(tInd)
Ainv = self.Solver(A, **self.solverOpts)
for i, src in enumerate(self.survey.srcList):
# compute next du_dm_v for next timestep
@@ -98,13 +93,13 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
# over-write with this time-steps (if not on last timestep)
if tInd != len(self.timeSteps):
dun_dm_v[:,i] = Adiaginv * rhs_v
dun_dm_v[:,i] = Ainv * 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()
Ainv.clean()
return Utils.mkvc(Jv)
@@ -248,9 +243,8 @@ class Problem_b(BaseTDEMProblem):
def getADeriv(self, tInd, u, v, adjoint=False):
C = self.mesh.edgeCurl
MeSigmaIDeriv = lambda u: self.MeSigmaIDeriv(u)
MeSigmaIDeriv = lambda x: self.MeSigmaIDeriv(x)
MfMui = self.MfMui
I = Utils.speye(self.mesh.nF)
if adjoint:
if self._makeASymmetric is True:
tests/em/tdem/test_TDEM_b_DerivAdjoint.py
+54 -9
View File
@@ -21,8 +21,8 @@ class TDEM_bDerivTests(unittest.TestCase):
activeMap = Maps.InjectActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
mapping = Maps.ExpMap(mesh) * Maps.SurjectVertical1D(mesh) * activeMap
rxOffset = 40.
rx = EM.TDEM.Rx(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20), 'bz')
rxOffset = 10.
rx = EM.TDEM.Rx(np.array([[rxOffset, 0., -1e-2]]), np.logspace(-4,-3, 20), 'bz')
src = EM.TDEM.SurveyTDEM.MagDipole([rx], loc=np.array([0., 0., 0.]))
survey = EM.TDEM.Survey([src])
@@ -30,6 +30,7 @@ class TDEM_bDerivTests(unittest.TestCase):
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.__makeASymmetric = False
# self.prb.timeSteps = [(1e-05, 100)]
try:
@@ -38,9 +39,10 @@ class TDEM_bDerivTests(unittest.TestCase):
except ImportError, e:
self.prb.Solver = SolverLU
self.sigma = np.ones(mesh.nCz)*1e-8
self.sigma[mesh.vectorCCz<0] = 1e-1
self.m = np.log(self.sigma[active])
# self.sigma = np.ones(mesh.nCz)*1e-8
# self.sigma[active] = 1e-1
# self.sigma[active] += 1e-2*np.random.rand(len(active))
self.m = np.log(1e-1)*np.ones(self.prb.mapping.nP) + 1e-2*np.random.randn(self.prb.mapping.nP)
self.prb.pair(survey)
self.mesh = mesh
@@ -174,21 +176,64 @@ class TDEM_bDerivTests(unittest.TestCase):
# print 'test_Deriv_dUdM'
# Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20)
def test_ADeriv(self):
prb = self.prb
tInd = 0
v = np.random.rand(self.mesh.nF)
def AderivTest(m):
prb.curModel = m
A = prb.getA(tInd)
Av = A*v
prb.curModel = self.m
ADeriv_dm = lambda dm: prb.getADeriv(tInd, v, dm)
return Av, ADeriv_dm
Tests.checkDerivative(AderivTest, self.m, plotIt=False, num=4, eps=1e-20)
def test_Fields_Deriv(self):
prb = self.prb
tInd = 0
v = np.random.rand(self.mesh.nF)
def FieldsDerivs(m):
prb.curModel = m
F = prb.fieldsPair(self.mesh, self.prb.survey)
# set initial fields
F[:,prb._fieldType+'Solution',0] = prb.getInitialFields()
A = prb.getA(tInd)
Ainv = prb.Solver(A)
RHS = prb.getRHS(tInd, F)
sol = Ainv * RHS
Ainv.clean()
prb.curModel = self.m
f = prb.fields(self.m)
deriv = lambda dm: f._bDeriv(prb.survey.srcList[0], np.zeros_like(sol), dm)
return sol, deriv
Tests.checkDerivative(FieldsDerivs, self.m, plotIt=False, num=4, eps=1e-20)
def test_Deriv_J(self):
prb = self.prb
prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)]
# 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_m = np.random.rand(prb.mapping.nP)
# d_m = 0.1*np.random.randn(prb.mapping.nP)
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)
Tests.checkDerivative(derChk, m, plotIt=False, num=5, eps=1e-20)
# def test_projectAdjoint(self):
# prb = self.prb
@@ -202,7 +247,7 @@ class TDEM_bDerivTests(unittest.TestCase):
# 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())