mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-07 10:40:46 +08:00
TDEM bderiv from b formulation working
This commit is contained in:
+7
-11
@@ -68,9 +68,12 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
|
||||
self.curModel = m
|
||||
|
||||
Jv = self.dataPair(self.survey)
|
||||
|
||||
# mat to store previous time-step's solution deriv times a vector for each source
|
||||
dun_dm_v = self.getInitialFieldsDeriv(v) # can over-write this at each timestep
|
||||
|
||||
df_dm = Fields_Derivs(self.mesh, self.survey)
|
||||
df_dm_v = Fields_Derivs(self.mesh, self.survey) # store the field derivs we need to project to calc full deriv
|
||||
|
||||
Ainv = None
|
||||
|
||||
for tInd, dt in enumerate(self.timeSteps):
|
||||
@@ -84,12 +87,12 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
|
||||
|
||||
for i, src in enumerate(self.survey.srcList):
|
||||
# compute next du_dm_v for next timestep
|
||||
u_src = u[src,ftype,tInd]
|
||||
u_src = u[src,ftype,tInd+1]
|
||||
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)
|
||||
df_dm_v[src, '%sDeriv'%rx.projField , tInd] = df_dmFun(src, dun_dm_v[:,i], v)
|
||||
|
||||
# over-write with this time-steps (if not on last timestep)
|
||||
if tInd != len(self.timeSteps):
|
||||
@@ -97,18 +100,11 @@ 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)
|
||||
Jv[src,rx] = rx.evalDeriv(src, self.mesh, self.timeMesh, df_dm_v)
|
||||
|
||||
Ainv.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):
|
||||
|
||||
@@ -11,12 +11,12 @@ class TDEM_bDerivTests(unittest.TestCase):
|
||||
|
||||
cs = 5.
|
||||
ncx = 20
|
||||
ncy = 6
|
||||
ncy = 10
|
||||
npad = 20
|
||||
hx = [(cs,ncx), (cs,npad,1.3)]
|
||||
hy = [(cs,npad,-1.3), (cs,ncy), (cs,npad,1.3)]
|
||||
mesh = Mesh.CylMesh([hx,1,hy], '00C')
|
||||
|
||||
#
|
||||
active = mesh.vectorCCz<0.
|
||||
activeMap = Maps.InjectActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
|
||||
mapping = Maps.ExpMap(mesh) * Maps.SurjectVertical1D(mesh) * activeMap
|
||||
@@ -193,31 +193,29 @@ class TDEM_bDerivTests(unittest.TestCase):
|
||||
|
||||
Tests.checkDerivative(AderivTest, self.m, plotIt=False, num=4, eps=1e-20)
|
||||
|
||||
def test_Fields_Deriv(self):
|
||||
prb = self.prb
|
||||
tInd = 0
|
||||
# def test_Fields_Deriv(self):
|
||||
# prb = self.prb
|
||||
# tInd = 10
|
||||
|
||||
v = np.random.rand(self.mesh.nF)
|
||||
# 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)
|
||||
# def FieldsDerivs(m):
|
||||
|
||||
sol = Ainv * RHS
|
||||
Ainv.clean()
|
||||
# sol = prb.fields(m)[:,'bSolution',tInd]
|
||||
|
||||
prb.curModel = self.m
|
||||
f = prb.fields(self.m)
|
||||
deriv = lambda dm: f._bDeriv(prb.survey.srcList[0], np.zeros_like(sol), dm)
|
||||
# prb.curModel = self.m
|
||||
# f = prb.fields(self.m)
|
||||
|
||||
return sol, deriv
|
||||
# df_dm_v = EM.TDEM.FieldsTDEM.Fields_Derivs(mesh, survey)
|
||||
# for i in range(tInd):
|
||||
# Ainv = prb.Solver(prb.getA(tInd))
|
||||
# df_dm_v = Ainv *
|
||||
|
||||
Tests.checkDerivative(FieldsDerivs, self.m, plotIt=False, num=4, eps=1e-20)
|
||||
# deriv = lambda dm: f._bDeriv(prb.survey.srcList[0], f[:,'bSolution',tInd], dm)
|
||||
|
||||
# return sol, deriv
|
||||
|
||||
# Tests.checkDerivative(FieldsDerivs, self.m, plotIt=False, num=4, eps=1e-20)
|
||||
|
||||
def test_Deriv_J(self):
|
||||
|
||||
@@ -233,7 +231,7 @@ class TDEM_bDerivTests(unittest.TestCase):
|
||||
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, num=5, eps=1e-20)
|
||||
Tests.checkDerivative(derChk, m, plotIt=False, num=2, eps=1e-20)
|
||||
|
||||
# def test_projectAdjoint(self):
|
||||
# prb = self.prb
|
||||
|
||||
Reference in New Issue
Block a user