mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-13 11:12:36 +08:00
Tested J and Jt for b and e formulations. Generalized code so it is easy to reuse. New receiver types.
This commit is contained in:
rowanc1
+88
-115
@@ -20,6 +20,7 @@ class BaseProblemFDEM(Problem.BaseProblem):
|
||||
def __init__(self, model, **kwargs):
|
||||
Problem.BaseProblem.__init__(self, model, **kwargs)
|
||||
|
||||
solType = None
|
||||
storeTheseFields = ['e', 'b']
|
||||
|
||||
surveyPair = SurveyFDEM
|
||||
@@ -49,7 +50,7 @@ class BaseProblemFDEM(Problem.BaseProblem):
|
||||
def MeSigma(self):
|
||||
#TODO: hardcoded to sigma as the model
|
||||
if getattr(self, '_MeSigma', None) is None:
|
||||
sigma = self.currentTransformedModel
|
||||
sigma = self.curTModel
|
||||
self._MeSigma = self.mesh.getEdgeInnerProduct(sigma)
|
||||
return self._MeSigma
|
||||
|
||||
@@ -60,10 +61,22 @@ class BaseProblemFDEM(Problem.BaseProblem):
|
||||
self._MeSigmaI = Utils.sdiag(1/self.MeSigma.diagonal())
|
||||
return self._MeSigmaI
|
||||
|
||||
currentTransformedModel = Utils.dependentProperty('_currentTransformedModel', None, ['_MeSigma', '_MeSigmaI'], 'Sets the current model, and removes dependent mass matrices.')
|
||||
curModel = Utils.dependentProperty('_curModel', None, ['_MeSigma', '_MeSigmaI', '_curTModel', '_curTModelDeriv'], 'Sets the current model, and removes dependent mass matrices.')
|
||||
|
||||
@property
|
||||
def curTModel(self):
|
||||
if getattr(self, '_curTModel', None) is None:
|
||||
self._curTModel = self.model.transform(self.curModel)
|
||||
return self._curTModel
|
||||
|
||||
@property
|
||||
def curTModelDeriv(self):
|
||||
if getattr(self, '_curTModelDeriv', None) is None:
|
||||
self._curTModelDeriv = self.model.transformDeriv(self.curModel)
|
||||
return self._curTModelDeriv
|
||||
|
||||
def fields(self, m):
|
||||
self.currentTransformedModel = self.model.transform(m)
|
||||
self.curModel = m
|
||||
F = self.forward(m, self.getRHS, self.calcFields)
|
||||
return F
|
||||
|
||||
@@ -80,10 +93,54 @@ class BaseProblemFDEM(Problem.BaseProblem):
|
||||
|
||||
return F
|
||||
|
||||
def Jvec(self, m, v, u=None):
|
||||
if u is None:
|
||||
u = self.fields(m)
|
||||
|
||||
self.curModel = m
|
||||
|
||||
Jv = self.dataPair(self.survey)
|
||||
|
||||
for freq in self.survey.freqs:
|
||||
A = self.getA(freq)
|
||||
solver = self.Solver(A, **self.solverOpts)
|
||||
|
||||
for tx in self.survey.getTransmitters(freq):
|
||||
w = self.getADeriv(freq, u[tx, self.solType], v)
|
||||
Ainvw = solver.solve(w)
|
||||
P = tx.projectFieldsDeriv(self.mesh, u)
|
||||
Jv[tx] = -P*Ainvw
|
||||
|
||||
return Utils.mkvc(Jv)
|
||||
|
||||
def Jtvec(self, m, v, u=None):
|
||||
if u is None:
|
||||
u = self.fields(m)
|
||||
|
||||
self.curModel = m
|
||||
|
||||
# Ensure v is a data object.
|
||||
if not isinstance(v, self.dataPair):
|
||||
v = self.dataPair(self.survey, v)
|
||||
|
||||
Jtv = np.zeros(self.model.nP, dtype=complex)
|
||||
|
||||
for freq in self.survey.freqs:
|
||||
AT = self.getA(freq).T
|
||||
solver = self.Solver(AT, **self.solverOpts)
|
||||
|
||||
for tx in self.survey.getTransmitters(freq):
|
||||
P = tx.projectFieldsDeriv(self.mesh, u)
|
||||
w = solver.solve( - P.T * v[tx])
|
||||
Jtv += self.getADeriv(freq, u[tx, self.solType], w, adjoint=True)
|
||||
|
||||
return Jtv
|
||||
|
||||
class ProblemFDEM_e(BaseProblemFDEM):
|
||||
"""
|
||||
Solving for e!
|
||||
"""
|
||||
solType = 'e'
|
||||
|
||||
def __init__(self, model, **kwargs):
|
||||
BaseProblemFDEM.__init__(self, model, **kwargs)
|
||||
@@ -100,6 +157,16 @@ class ProblemFDEM_e(BaseProblemFDEM):
|
||||
|
||||
return C.T*mui*C + 1j*omega(freq)*sig
|
||||
|
||||
def getADeriv(self, freq, u, v, adjoint=False):
|
||||
sig = self.curTModel
|
||||
dsig_dm = self.curTModelDeriv
|
||||
dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig, v=u)
|
||||
|
||||
if adjoint:
|
||||
return 1j * omega(freq) * ( dsig_dm.T * ( dMe_dsig.T * v ) )
|
||||
|
||||
return 1j * omega(freq) * ( dMe_dsig * ( dsig_dm * v ) )
|
||||
|
||||
def getRHS(self, freq):
|
||||
"""
|
||||
:param float freq: Frequency
|
||||
@@ -139,65 +206,12 @@ class ProblemFDEM_e(BaseProblemFDEM):
|
||||
|
||||
return fDict
|
||||
|
||||
def Jvec(self, m, v, u=None):
|
||||
if u is None:
|
||||
u = self.fields(m)
|
||||
|
||||
sig = self.model.transform(m)
|
||||
self.currentTransformedModel = sig
|
||||
|
||||
Jv = self.dataPair(self.survey)
|
||||
dsig_dm = self.model.transformDeriv(m)
|
||||
|
||||
for i, freq in enumerate(self.survey.freqs):
|
||||
e = u[freq, 'e']
|
||||
A = self.getA(freq)
|
||||
solver = self.Solver(A, **self.solverOpts)
|
||||
|
||||
for txi, tx in enumerate(self.survey.getTransmitters(freq)):
|
||||
dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig, v=e[:,txi])
|
||||
|
||||
P = tx.projectFieldsDeriv(self.mesh, u)
|
||||
b = 1j*omega(freq) * ( dMe_dsig * ( dsig_dm * v ) )
|
||||
Ainvb = solver.solve(b)
|
||||
Jv[tx] = -P*Ainvb
|
||||
|
||||
return Utils.mkvc(Jv)
|
||||
|
||||
|
||||
def Jtvec(self, m, v, u=None):
|
||||
if u is None:
|
||||
u = self.fields(m)
|
||||
|
||||
sig = self.model.transform(m)
|
||||
self.currentTransformedModel = sig
|
||||
|
||||
# Ensure v is a data object.
|
||||
if not isinstance(v, self.dataPair):
|
||||
v = self.dataPair(self.survey, v)
|
||||
|
||||
Jtv = np.zeros(self.model.nP, dtype=complex)
|
||||
|
||||
dsig_dm = self.model.transformDeriv(m)
|
||||
|
||||
for i, freq in enumerate(self.survey.freqs):
|
||||
e = u[freq, 'e']
|
||||
AT = self.getA(freq).T
|
||||
solver = self.Solver(AT, **self.solverOpts)
|
||||
|
||||
for txi, tx in enumerate(self.survey.getTransmitters(freq)):
|
||||
dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig, v=e[:,txi])
|
||||
|
||||
P = tx.projectFieldsDeriv(self.mesh, u)
|
||||
w = solver.solve(P.T * v[tx])
|
||||
Jtv += - 1j*omega(freq) * ( dsig_dm.T * ( dMe_dsig.T * w ) )
|
||||
|
||||
return Jtv
|
||||
|
||||
class ProblemFDEM_b(BaseProblemFDEM):
|
||||
"""
|
||||
Solving for b!
|
||||
"""
|
||||
solType = 'b'
|
||||
|
||||
def __init__(self, model, **kwargs):
|
||||
BaseProblemFDEM.__init__(self, model, **kwargs)
|
||||
@@ -214,6 +228,23 @@ class ProblemFDEM_b(BaseProblemFDEM):
|
||||
|
||||
return mui*C*sigI*C.T*mui + 1j*omega(freq)*mui
|
||||
|
||||
def getADeriv(self, freq, u, v, adjoint=False):
|
||||
|
||||
mui = self.MfMui
|
||||
C = self.mesh.edgeCurl
|
||||
sig = self.curTModel
|
||||
dsig_dm = self.curTModelDeriv
|
||||
#TODO: This only works if diagonal (no tensors)...
|
||||
dMeSigmaI_dI = - self.MeSigmaI**2
|
||||
|
||||
vec = (C.T*(mui*u))
|
||||
dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig, v=vec)
|
||||
|
||||
if adjoint:
|
||||
return dsig_dm.T * ( dMe_dsig.T * ( dMeSigmaI_dI.T * ( C.T * ( mui.T * v ) ) ) )
|
||||
|
||||
return mui * ( C * ( dMeSigmaI_dI * ( dMe_dsig * ( dsig_dm * v ) ) ) )
|
||||
|
||||
def getRHS(self, freq):
|
||||
"""
|
||||
:param float freq: Frequency
|
||||
@@ -253,61 +284,3 @@ class ProblemFDEM_b(BaseProblemFDEM):
|
||||
|
||||
return fDict
|
||||
|
||||
def Jvec(self, m, v, u=None):
|
||||
if u is None:
|
||||
u = self.fields(m)
|
||||
|
||||
raise NotImplemented('')
|
||||
|
||||
# sig = self.model.transform(m)
|
||||
# self.currentTransformedModel = sig
|
||||
|
||||
# Jv = self.dataPair(self.survey)
|
||||
# dsig_dm = self.model.transformDeriv(m)
|
||||
|
||||
# for i, freq in enumerate(self.survey.freqs):
|
||||
# e = u[freq, 'e']
|
||||
# A = self.getA(freq)
|
||||
# solver = self.Solver(A, **self.solverOpts)
|
||||
|
||||
# for txi, tx in enumerate(self.survey.getTransmitters(freq)):
|
||||
# dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig, v=e[:,txi])
|
||||
|
||||
# P = tx.projectFieldsDeriv(self.mesh, u)
|
||||
# b = 1j*omega(freq) * ( dMe_dsig * ( dsig_dm * v ) )
|
||||
# Ainvb = solver.solve(b)
|
||||
# Jv[tx] = -P*Ainvb
|
||||
|
||||
# return Utils.mkvc(Jv)
|
||||
|
||||
|
||||
def Jtvec(self, m, v, u=None):
|
||||
if u is None:
|
||||
u = self.fields(m)
|
||||
|
||||
# Ensure v is a data object.
|
||||
if not isinstance(v, self.dataPair):
|
||||
v = self.dataPair(self.survey, v)
|
||||
|
||||
raise NotImplemented('')
|
||||
|
||||
# sig = self.model.transform(m)
|
||||
# self.currentTransformedModel = sig
|
||||
|
||||
# Jtv = np.zeros(self.model.nP, dtype=complex)
|
||||
|
||||
# dsig_dm = self.model.transformDeriv(m)
|
||||
|
||||
# for i, freq in enumerate(self.survey.freqs):
|
||||
# e = u[freq, 'e']
|
||||
# AT = self.getA(freq).T
|
||||
# solver = self.Solver(AT, **self.solverOpts)
|
||||
|
||||
# for txi, tx in enumerate(self.survey.getTransmitters(freq)):
|
||||
# dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig, v=e[:,txi])
|
||||
|
||||
# P = tx.projectFieldsDeriv(self.mesh, u)
|
||||
# w = solver.solve(P.T * v[tx])
|
||||
# Jtv += - 1j*omega(freq) * ( dsig_dm.T * ( dMe_dsig.T * w ) )
|
||||
|
||||
# return Jtv
|
||||
|
||||
@@ -2,7 +2,15 @@ from SimPEG import Survey, Utils, np, sp
|
||||
|
||||
class RxListFDEM(Survey.BaseRxList):
|
||||
|
||||
knownRxTypes = ['Ex', 'Ey', 'Ez']
|
||||
knownRxTypes = {
|
||||
'Ex':'Ex',
|
||||
'Ey':'Ey',
|
||||
'Ez':'Ez',
|
||||
|
||||
'Bx':'Fx',
|
||||
'By':'Fy',
|
||||
'Bz':'Fz',
|
||||
}
|
||||
|
||||
def __init__(self, locs, rxType):
|
||||
Survey.BaseRxList.__init__(self, locs, rxType)
|
||||
@@ -11,7 +19,7 @@ class RxListFDEM(Survey.BaseRxList):
|
||||
|
||||
def getP(self, mesh):
|
||||
if mesh not in self._Ps:
|
||||
self._Ps[mesh] = mesh.getInterpolationMat(self.locs, self.rxType)
|
||||
self._Ps[mesh] = mesh.getInterpolationMat(self.locs, self.knownRxTypes[self.rxType])
|
||||
return self._Ps[mesh]
|
||||
|
||||
|
||||
@@ -37,6 +45,8 @@ class TxFDEM(Survey.BaseTx):
|
||||
|
||||
if self.rxList.rxType in ['Ex', 'Ey', 'Ez']:
|
||||
u_part = u[self, 'e']
|
||||
elif self.rxList.rxType in ['Bx', 'By', 'Bz']:
|
||||
u_part = u[self, 'b']
|
||||
else:
|
||||
raise NotImplemented('Unknown receiver type.')
|
||||
|
||||
|
||||
+76
-30
@@ -4,44 +4,90 @@ import simpegEM as EM
|
||||
|
||||
TOL = 1e-10
|
||||
|
||||
class FDEM_bDerivTests(unittest.TestCase):
|
||||
def getProblem(fdemType):
|
||||
cs = 5.
|
||||
ncx, ncy, ncz = 2, 2, 2
|
||||
npad = 3
|
||||
hx = Utils.meshTensors(((npad,cs), (ncx,cs), (npad,cs)))
|
||||
hy = Utils.meshTensors(((npad,cs), (ncy,cs), (npad,cs)))
|
||||
hz = Utils.meshTensors(((npad,cs), (ncz,cs), (npad,cs)))
|
||||
mesh = Mesh.TensorMesh([hx,hy,hz])
|
||||
|
||||
model = Model.LogModel(mesh)
|
||||
|
||||
x = np.linspace(5,10,3)
|
||||
XYZ = Utils.ndgrid(x,x,np.r_[0])
|
||||
if fdemType == 'e':
|
||||
rxList = EM.FDEM.RxListFDEM(XYZ, 'Ex')
|
||||
elif fdemType == 'b':
|
||||
rxList = EM.FDEM.RxListFDEM(XYZ, 'Bx')
|
||||
else:
|
||||
raise NotImplementedError()
|
||||
|
||||
Tx0 = EM.FDEM.TxFDEM(np.r_[4.,2.,2.], 'VMD', 1e-2, rxList)
|
||||
|
||||
x = np.linspace(5,10,3)
|
||||
XYZ = Utils.ndgrid(x,x,np.r_[0])
|
||||
if fdemType == 'e':
|
||||
rxList = EM.FDEM.RxListFDEM(XYZ, 'Ey')
|
||||
elif fdemType == 'b':
|
||||
rxList = EM.FDEM.RxListFDEM(XYZ, 'By')
|
||||
else:
|
||||
raise NotImplementedError()
|
||||
|
||||
Tx1 = EM.FDEM.TxFDEM(np.r_[4.,2.,2.], 'VMD', 1e-4, rxList)
|
||||
|
||||
survey = EM.FDEM.SurveyFDEM([Tx0, Tx1])
|
||||
|
||||
if fdemType == 'e':
|
||||
prb = EM.FDEM.ProblemFDEM_e(model)
|
||||
elif fdemType == 'b':
|
||||
prb = EM.FDEM.ProblemFDEM_b(model)
|
||||
else:
|
||||
raise NotImplementedError()
|
||||
prb.pair(survey)
|
||||
|
||||
return prb
|
||||
|
||||
class FDEM_DerivTests_e(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
|
||||
cs = 5.
|
||||
ncx, ncy, ncz = 2, 2, 2
|
||||
npad = 3
|
||||
hx = Utils.meshTensors(((npad,cs), (ncx,cs), (npad,cs)))
|
||||
hy = Utils.meshTensors(((npad,cs), (ncy,cs), (npad,cs)))
|
||||
hz = Utils.meshTensors(((npad,cs), (ncz,cs), (npad,cs)))
|
||||
mesh = Mesh.TensorMesh([hx,hy,hz])
|
||||
|
||||
model = Model.LogModel(mesh)
|
||||
|
||||
x = np.linspace(5,10,3)
|
||||
XYZ = Utils.ndgrid(x,x,np.r_[0])
|
||||
rxList = EM.FDEM.RxListFDEM(XYZ, 'Ex')
|
||||
Tx0 = EM.FDEM.TxFDEM(np.r_[4.,2.,2.], 'VMD', 1e-2, rxList)
|
||||
|
||||
x = np.linspace(5,10,3)
|
||||
XYZ = Utils.ndgrid(x,x,np.r_[0])
|
||||
rxList = EM.FDEM.RxListFDEM(XYZ, 'Ey')
|
||||
Tx1 = EM.FDEM.TxFDEM(np.r_[4.,2.,2.], 'VMD', 1e-4, rxList)
|
||||
|
||||
survey = EM.FDEM.SurveyFDEM([Tx0, Tx1])
|
||||
|
||||
prb = EM.FDEM.ProblemFDEM_e(model)
|
||||
prb.pair(survey)
|
||||
|
||||
self.sigma = np.log(np.ones(mesh.nC)*1e-3)
|
||||
prb = getProblem('e')
|
||||
self.prb = prb
|
||||
self.survey = survey
|
||||
self.sigma = np.log(np.ones(prb.mesh.nC)*1e-3)
|
||||
self.survey = prb.survey
|
||||
|
||||
def test_Jvec(self):
|
||||
x0 = self.sigma
|
||||
def fun(x):
|
||||
return self.survey.dpred(x), lambda x: self.prb.Jvec(x0, x)
|
||||
passed = Tests.checkDerivative(fun, x0, num=3, plotIt=False, eps=1e-18)
|
||||
passed = Tests.checkDerivative(fun, x0, num=3, plotIt=False, eps=1e-25)
|
||||
self.assertTrue(passed)
|
||||
|
||||
def test_Jtvec_adjointTest(self):
|
||||
v = np.random.rand(self.survey.nD)
|
||||
w = np.random.rand(self.prb.model.nP)
|
||||
|
||||
m = self.sigma
|
||||
u = self.prb.fields(m)
|
||||
vJw = v.dot(self.prb.Jvec(m, w, u=u))
|
||||
wJtv = w.dot(self.prb.Jtvec(m, v, u=u))
|
||||
self.assertTrue(vJw - wJtv < TOL)
|
||||
|
||||
|
||||
class FDEM_DerivTests_b(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
prb = getProblem('e')
|
||||
self.prb = prb
|
||||
self.sigma = np.log(np.ones(prb.mesh.nC)*1e-3)
|
||||
self.survey = prb.survey
|
||||
|
||||
def test_Jvec(self):
|
||||
x0 = self.sigma
|
||||
def fun(x):
|
||||
return self.survey.dpred(x), lambda x: self.prb.Jvec(x0, x)
|
||||
passed = Tests.checkDerivative(fun, x0, num=3, plotIt=False, eps=1e-25)
|
||||
self.assertTrue(passed)
|
||||
|
||||
def test_Jtvec_adjointTest(self):
|
||||
|
||||
Reference in New Issue
Block a user