mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-11 06:23:06 +08:00
e hooked up with Jvec
This commit is contained in:
Lindsey Heagy
@@ -49,35 +49,49 @@ class Fields_b(Fields):
|
||||
}
|
||||
|
||||
def startup(self):
|
||||
self.MeSigmaI = self.survey.prob.MeSigmaI
|
||||
self.edgeCurl = self.survey.prob.mesh.edgeCurl
|
||||
self.MfMui = self.survey.prob.MfMui
|
||||
self.MeSigmaI = self.survey.prob.MeSigmaI
|
||||
self.MeSigmaIDeriv = self.survey.prob.MeSigmaIDeriv
|
||||
self.edgeCurl = self.survey.prob.mesh.edgeCurl
|
||||
self.MfMui = self.survey.prob.MfMui
|
||||
|
||||
def _b(self, bSolution, srcList, tInd):
|
||||
return bSolution
|
||||
|
||||
def _bDeriv_u(self, src, dun_dm_v, adjoint = False):
|
||||
def _bDeriv_u(self, tInd, src, dun_dm_v, adjoint = False):
|
||||
return Identity()*dun_dm_v
|
||||
|
||||
def _bDeriv_m(self, src, v, adjoint = False):
|
||||
def _bDeriv_m(self, tInd, src, v, adjoint = False):
|
||||
return Zero()
|
||||
|
||||
def _bDeriv(self, src, dun_dm_v, v, adjoint=False):
|
||||
def _bDeriv(self, tInd, 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)
|
||||
return self._bDeriv_u(tInd, src, dun_dm_v) + self._bDeriv_m(tInd, src, v)
|
||||
|
||||
def _e(self, bSolution, srcList, tInd):
|
||||
e = self.MeSigmaI * ( self.edgeCurl.T * ( self.MfMui * bSolution ) )
|
||||
for i, src in enumerate(srcList):
|
||||
_, S_e = src.eval(self.prob, tInd)
|
||||
e[:,i,tInd] = e[:,i,tInd] - self.MeSigmaI * S_e
|
||||
_, S_e = src.eval(self.survey.prob, self.survey.prob.times[tInd])
|
||||
e[:,i] = e[:,i] - self.MeSigmaI * S_e
|
||||
return e
|
||||
|
||||
def _eDeriv_u(self, src, dun_dm_v, adjoint = False):
|
||||
raise NotImplementedError
|
||||
def _eDeriv_u(self, tInd, src, dun_dm_v, adjoint = False):
|
||||
if adjoint is True:
|
||||
raise NotImplementedError
|
||||
return self.MeSigmaI * ( self.edgeCurl.T * ( self.MfMui * dun_dm_v ) )
|
||||
|
||||
def _eDeriv_m(self, src, v, adjoint = False):
|
||||
raise NotImplementedError
|
||||
def _eDeriv_m(self, tInd, src, v, adjoint = False):
|
||||
if adjoint is True:
|
||||
raise NotImplementedError
|
||||
|
||||
bSolution = self[[src],'bSolution',tInd]
|
||||
|
||||
_, S_e = src.eval(self.survey.prob, self.survey.prob.times[tInd])
|
||||
_, S_eDeriv = src.evalDeriv(self.survey.prob.times[tInd], self, v=v)
|
||||
|
||||
return self.MeSigmaIDeriv(self.edgeCurl.T * ( self.MfMui * bSolution) ) * v - self.MeSigmaIDeriv(S_e) * v - self.MeSigmaI * S_eDeriv
|
||||
|
||||
def _eDeriv(self, tInd, src, dun_dm_v, v, adjoint=False):
|
||||
if adjoint is True:
|
||||
raise NotImplementedError
|
||||
return self._eDeriv_u(tInd, src, dun_dm_v) + self._eDeriv_m(tInd, src, v)
|
||||
|
||||
@@ -92,7 +92,7 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
|
||||
|
||||
for rx in src.rxList:
|
||||
df_dmFun = getattr(u, '_%sDeriv'%rx.projField, None)
|
||||
df_dm_v[src, '%sDeriv'%rx.projField , tInd] = df_dmFun(src, dun_dm_v[:,i], v)
|
||||
df_dm_v[src, '%sDeriv'%rx.projField , tInd] = df_dmFun(tInd, src, dun_dm_v[:,i], v)
|
||||
|
||||
# over-write with this time-steps (if not on last timestep)
|
||||
if tInd != len(self.timeSteps):
|
||||
|
||||
@@ -5,234 +5,81 @@ from SimPEG import EM
|
||||
plotIt = False
|
||||
tol = 1e-6
|
||||
|
||||
def setUp(rxcomp='bz'):
|
||||
cs = 5.
|
||||
ncx = 20
|
||||
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
|
||||
|
||||
rxOffset = 10.
|
||||
rx = EM.TDEM.Rx(np.array([[rxOffset, 0., -1e-2]]), np.logspace(-4,-3, 20), rxcomp)
|
||||
src = EM.TDEM.SurveyTDEM.MagDipole([rx], loc=np.array([0., 0., 0.]))
|
||||
|
||||
survey = EM.TDEM.Survey([src])
|
||||
|
||||
prb = EM.TDEM.Problem_b(mesh, mapping=mapping)
|
||||
prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)]
|
||||
|
||||
try:
|
||||
from pymatsolver import MumpsSolver
|
||||
prb.Solver = MumpsSolver
|
||||
except ImportError, e:
|
||||
prb.Solver = SolverLU
|
||||
|
||||
m = np.log(1e-1)*np.ones(prb.mapping.nP) + 1e-2*np.random.randn(prb.mapping.nP)
|
||||
|
||||
prb.pair(survey)
|
||||
mesh = mesh
|
||||
|
||||
return prb, m, mesh
|
||||
|
||||
class TDEM_bDerivTests(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
|
||||
cs = 5.
|
||||
ncx = 20
|
||||
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
|
||||
|
||||
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])
|
||||
|
||||
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:
|
||||
from pymatsolver import MumpsSolver
|
||||
self.prb.Solver = MumpsSolver
|
||||
except ImportError, e:
|
||||
self.prb.Solver = SolverLU
|
||||
|
||||
# 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
|
||||
|
||||
# def test_AhVec(self):
|
||||
# """
|
||||
# Test that fields and AhVec produce consistent results
|
||||
# """
|
||||
|
||||
# prb = self.prb
|
||||
# sigma = self.sigma
|
||||
|
||||
# u = prb.fields(sigma)
|
||||
# Ahu = prb._AhVec(sigma, u)
|
||||
|
||||
# V1 = Ahu[:,'b',1]
|
||||
# V2 = 1./prb.timeSteps[0]*prb.MfMui*u[:,'b',0]
|
||||
# self.assertLess(np.linalg.norm(V1-V2)/np.linalg.norm(V2), 1.e-6)
|
||||
|
||||
# V1 = Ahu[:,'e',1]
|
||||
# return np.linalg.norm(V1) < 1.e-6
|
||||
|
||||
# for i in range(2,prb.nT):
|
||||
|
||||
# dt = prb.timeSteps[i]
|
||||
|
||||
# V1 = Ahu[:,'b',i]
|
||||
# V2 = 1.0/dt*prb.MfMui*u[:,'b', i-1]
|
||||
# # print np.linalg.norm(V1), np.linalg.norm(V2)
|
||||
# self.assertLess(np.linalg.norm(V1)/np.linalg.norm(V2), 1.e-6)
|
||||
|
||||
# V1 = Ahu[:,'e',i]
|
||||
# V2 = prb.MeSigma*u[:,'e',i]
|
||||
# # print np.linalg.norm(V1), np.linalg.norm(V2)
|
||||
# return np.linalg.norm(V1)/np.linalg.norm(V2), 1.e-6
|
||||
|
||||
# def test_AhVecVSMat_OneTS(self):
|
||||
|
||||
# prb = self.prb
|
||||
# prb.timeSteps = [1e-05]
|
||||
# sigma = self.sigma
|
||||
# prb.curModel = sigma
|
||||
|
||||
# dt = prb.timeSteps[0]
|
||||
# a11 = 1/dt*prb.MfMui*sp.identity(prb.mesh.nF)
|
||||
# a12 = prb.MfMui*prb.mesh.edgeCurl
|
||||
# a21 = prb.mesh.edgeCurl.T*prb.MfMui
|
||||
# a22 = -prb.MeSigma
|
||||
# A = sp.bmat([[a11,a12],[a21,a22]])
|
||||
|
||||
# f = prb.fields(sigma)
|
||||
# u1 = A*f.tovec()
|
||||
# u2 = prb._AhVec(sigma,f).tovec()
|
||||
|
||||
# self.assertTrue(np.linalg.norm(u1-u2)/np.linalg.norm(u1)<1e-12)
|
||||
|
||||
# def test_solveAhVSMat_OneTS(self):
|
||||
# prb = self.prb
|
||||
|
||||
# prb.timeSteps = [1e-05]
|
||||
|
||||
# sigma = self.sigma
|
||||
# prb.curModel = sigma
|
||||
|
||||
# dt = prb.timeSteps[0]
|
||||
# a11 = 1.0/dt*prb.MfMui*sp.identity(prb.mesh.nF)
|
||||
# a12 = prb.MfMui*prb.mesh.edgeCurl
|
||||
# a21 = prb.mesh.edgeCurl.T*prb.MfMui
|
||||
# a22 = -prb.MeSigma
|
||||
# A = sp.bmat([[a11,a12],[a21,a22]])
|
||||
|
||||
# f = prb.fields(sigma)
|
||||
# f[:,:,0] = {'b':0}
|
||||
# f[:,'b',1] = 0
|
||||
|
||||
# self.assertTrue(np.all(np.r_[f[:,'b',1],f[:,'e',1]] == f.tovec()))
|
||||
|
||||
# u1 = prb.solveAh(sigma,f).tovec().flatten()
|
||||
# u2 = sp.linalg.spsolve(A.tocsr(),f.tovec())
|
||||
|
||||
# self.assertTrue(np.linalg.norm(u1-u2)<1e-8)
|
||||
|
||||
# def test_solveAhVsAhVec(self):
|
||||
|
||||
# prb = self.prb
|
||||
# mesh = self.prb.mesh
|
||||
# sigma = self.sigma
|
||||
# self.prb.curModel = sigma
|
||||
|
||||
# f = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
|
||||
# f[:,'b',:] = 0.0
|
||||
# for i in range(prb.nT):
|
||||
# f[:,'e', i] = np.random.rand(mesh.nE, 1)
|
||||
|
||||
# Ahf = prb._AhVec(sigma, f)
|
||||
# f_test = prb.solveAh(sigma, Ahf)
|
||||
|
||||
# u1 = f.tovec()
|
||||
# u2 = f_test.tovec()
|
||||
# self.assertTrue(np.linalg.norm(u1-u2)<1e-8)
|
||||
|
||||
# def test_DerivG(self):
|
||||
# """
|
||||
# Test the derivative of c with respect to sigma
|
||||
# """
|
||||
|
||||
# # Random model and perturbation
|
||||
# sigma = np.random.rand(self.prb.mapping.nP)
|
||||
|
||||
# f = self.prb.fields(sigma)
|
||||
# dm = 1000*np.random.rand(self.prb.mapping.nP)
|
||||
# h = 0.01
|
||||
|
||||
# derChk = lambda m: [self.prb._AhVec(m, f).tovec(), lambda mx: self.prb.Gvec(sigma, mx, u=f).tovec()]
|
||||
# print '\ntest_DerivG'
|
||||
# passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20)
|
||||
# return passed
|
||||
|
||||
# def test_Deriv_dUdM(self):
|
||||
|
||||
# prb = self.prb
|
||||
# prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)]
|
||||
# mesh = self.mesh
|
||||
# sigma = self.sigma
|
||||
|
||||
# dm = 10*np.random.rand(prb.mapping.nP)
|
||||
# f = prb.fields(sigma)
|
||||
|
||||
# derChk = lambda m: [self.prb.fields(m).tovec(), lambda mx: -prb.solveAh(sigma, prb.Gvec(sigma, mx, u=f)).tovec()]
|
||||
# print '\n'
|
||||
# print 'test_Deriv_dUdM'
|
||||
# Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20)
|
||||
|
||||
|
||||
def test_ADeriv(self):
|
||||
prb = self.prb
|
||||
prb, m0, mesh = setUp()
|
||||
tInd = 0
|
||||
|
||||
v = np.random.rand(self.mesh.nF)
|
||||
v = np.random.rand(mesh.nF)
|
||||
|
||||
def AderivTest(m):
|
||||
prb.curModel = m
|
||||
A = prb.getA(tInd)
|
||||
Av = A*v
|
||||
prb.curModel = self.m
|
||||
prb.curModel = m0
|
||||
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)
|
||||
Tests.checkDerivative(AderivTest, m0, plotIt=False, num=4, eps=1e-20)
|
||||
|
||||
# def test_Fields_Deriv(self):
|
||||
# prb = self.prb
|
||||
# tInd = 10
|
||||
|
||||
# v = np.random.rand(self.mesh.nF)
|
||||
|
||||
# def FieldsDerivs(m):
|
||||
|
||||
# sol = prb.fields(m)[:,'bSolution',tInd]
|
||||
|
||||
# prb.curModel = self.m
|
||||
# f = prb.fields(self.m)
|
||||
|
||||
# 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 *
|
||||
|
||||
# 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):
|
||||
|
||||
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_m = 0.1*np.random.randn(prb.mapping.nP)
|
||||
|
||||
def JvecTest(self, rxcomp):
|
||||
prb, m, mesh = setUp(rxcomp)
|
||||
|
||||
derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(m, mx)]
|
||||
print '\n'
|
||||
print 'test_Deriv_J'
|
||||
print 'test_Jvec_%s' %(rxcomp)
|
||||
Tests.checkDerivative(derChk, m, plotIt=False, num=2, eps=1e-20)
|
||||
|
||||
def test_Jvec_b_bx(self):
|
||||
self.JvecTest('bx')
|
||||
|
||||
def test_Jvec_b_bz(self):
|
||||
self.JvecTest('bz')
|
||||
|
||||
def test_Jvec_b_ey(self):
|
||||
self.JvecTest('ey')
|
||||
|
||||
|
||||
# def test_projectAdjoint(self):
|
||||
# prb = self.prb
|
||||
# survey = prb.survey
|
||||
|
||||
Reference in New Issue
Block a user