e,b,h,j with deriv and adjoint from e formulation

This commit is contained in:
Lindsey Heagy
2016-02-20 11:05:25 -08:00
parent 4c3c2c361c
commit ce49249664
4 changed files with 144 additions and 152 deletions
+63 -71
View File
@@ -178,7 +178,11 @@ class Fields_e(Fields):
'bPrimary' : ['eSolution','F','_bPrimary'],
'bSecondary' : ['eSolution','F','_bSecondary'],
'j' : ['eSolution','CCV','_j'],
'jPrimary' : ['eSolution','CCV','_jPrimary'],
'jSecondary' : ['eSolution','CCV','_jSecondary'],
'h' : ['eSolution','CCV','_h'],
'hPrimary' : ['eSolution','CCV','_hPrimary'],
'hSecondary' : ['eSolution','CCV','_hSecondary'],
}
def __init__(self, mesh, survey, **kwargs):
@@ -192,6 +196,7 @@ class Fields_e(Fields):
self._nC = self.survey.prob.mesh.nC
self._MeSigma = self.survey.prob.MeSigma
self._MeSigmaDeriv = self.survey.prob.MeSigmaDeriv
self._MfMui = self.survey.prob.MfMui
def _GLoc(self, fieldType):
if fieldType == 'e':
@@ -203,6 +208,7 @@ class Fields_e(Fields):
else:
raise Exception('Field type must be e, b, h, j')
def _ePrimary(self, eSolution, srcList):
"""
Primary electric field from source
@@ -293,9 +299,9 @@ class Fields_e(Fields):
b[:,i] = b[:,i]+ 1./(1j*omega(src.freq)) * S_m
return b
def _bSecondaryDeriv_u(self, src, du_dm_v, adjoint = False):
def _bDeriv_u(self, src, du_dm_v, adjoint = False):
"""
Derivative of the secondary magnetic flux density with respect to the thing we solved for
Derivative of the magnetic flux density with respect to the thing we solved for
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray du_dm_v: vector to take product with
@@ -310,94 +316,80 @@ class Fields_e(Fields):
return - 1./(1j*omega(src.freq)) * (C * du_dm_v)
def _bSecondaryDeriv_m(self, src, v, adjoint = False):
def _bDeriv_m(self, src, v, adjoint = False):
"""
Derivative of the secondary magnetic flux density with respect to the inversion model.
Derivative of the magnetic flux density with respect to the inversion model.
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
:return: product of the secondary magnetic flux density derivative with respect to the inversion model with a vector
:return: product of the magnetic flux density derivative with respect to the inversion model with a vector
"""
S_mDeriv, _ = src.evalDeriv(self.prob, v, adjoint)
return 1./(1j * omega(src.freq)) * S_mDeriv
def _bDeriv_u(self, src, du_dm_v, adjoint=False):
"""
Partial derivative of the total magnetic flux density with respect to the thing we solved for
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray du_dm_v: vector to take product with
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
:return: product of the derivative of the magnetic flux density with respect to the field we solved for with a vector
"""
return self._bSecondaryDeriv_u(src, du_dm_v, adjoint)
def _bDeriv_m(self, src, v, adjoint=False):
"""
Partial derivative of the total magnetic flux density with respect to the inversion model.
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
:rtype: SimPEG.Utils.Zero
:return: product of the magnetic flux density derivative with respect to the inversion model with a vector
"""
# Assuming the primary does not depend on the model
return self._bSecondaryDeriv_m(src, eSolution, v, adjoint)
def _j(self, eSolution, srcList):
def _jSecondary(self, eSolution, srcList):
aveE2CCV = self._aveE2CCV
n = int(aveE2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy
Sigma = self._MeSigma
VI = sdiag(1./np.kron(np.ones(n), self.prob.mesh.vol))
return VI * (aveE2CCV * (Sigma *eSolution) )
n = int(aveE2CCV.shape[0] / self._nC) # number of components (instead of checking if cyl or not)
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
return VI * (aveE2CCV * (self._MeSigma * eSolution) )
def _jDeriv_u(self, src, v, adjoint = False):
eSolution = self._eSecondary
def _jPrimary(self, eSolution, srcList):
aveE2CCV = self._aveE2CCV
n = int(aveE2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy
Sigma = self._MeSigma
n = int(aveE2CCV.shape[0] / self._nC) # number of components (instead of checking if cyl or not)
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
return VI * (aveE2CCV * (self._MeSigma * self._ePrimary(eSolution, srcList)) )
VI = sdiag(1./np.kron(np.ones(n), self.prob.mesh.vol))
def _jDeriv_u(self, src, du_dm_v, adjoint = False):
n = int(self._aveE2CCV.shape[0] / self._nC) # number of components (instead of checking if cyl or not)
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
if not adjoint:
return VI * (aveE2CCV * (Sigma * (self._eDeriv_u(src, v, adjoint) ) ) )
return self._eDeriv_u(src, Sigma.T * (aveE2CCV.T * (VI.T * v) ), adjoint)
if adjoint:
return self._eDeriv_u(src, self._MeSigma.T * (self._aveE2CCV.T * (VI.T * du_dm_v) ), adjoint)
return VI * (self._aveE2CCV * (self._MeSigma * (self._eDeriv_u(src, du_dm_v, adjoint) ) ) )
def _jDeriv_m(self, src, v, adjoint = False):
eSolution = self._fields['eSolution']
e = self._e(self._fields['eSolution'], [src])
n = int(self._aveE2CCV.shape[0] / self._nC) #number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
aveE2CCV = self._aveE2CCV
Sigma = self._MeSigma
SigmaDeriv = self._MeSigmaDeriv
e = self._e(eSolution, [src])
n = int(aveE2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy
if adjoint:
return self._MeSigmaDeriv(e).T * (self._aveE2CCV.T * (VI.T * v)) + self._eDeriv_m(src, self._aveE2CCV.T * (VI.T * v), adjoint=adjoint)
return VI * (self._aveE2CCV * ( self._eDeriv_m(src, v, adjoint=adjoint) + self._MeSigmaDeriv(e) * v))
VI = sdiag(1./np.kron(np.ones(n), self.prob.mesh.vol))
if not adjoint:
return VI * (aveE2CCV * ( self._eDeriv_m(src, v, adjoint=adjoint) + SigmaDeriv(e) * v))
return SigmaDeriv(aveE2CCV.T * (VI.T * e), adjoint=adjoint) * v + self._eDeriv_m(src, aveE2CCV.T * (VI.T * v), adjoint=adjoint)
def _h(self, eSolution, srcList):
b = self._b(eSolution, srcList)
Mui = self.survey.prob.MfMui
aveF2CCV = self._aveF2CCV
n = int(aveF2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy
# Mui = sdiag(sp.kron(np.ones(n), mui))
VI = sdiag(1./np.kron(np.ones(n), self.prob.mesh.vol))
def _hPrimary(self, eSolution, srcList):
n = int(self._aveF2CCV.shape[0] / self._nC) # Number of Components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
return VI * (self._aveF2CCV * (self._MfMui * self._bPrimary(eSolution, srcList)))
return VI * (aveF2CCV * (Mui * b))
def _hSecondary(self, eSolution, srcList):
n = int(self._aveF2CCV.shape[0] / self._nC) # Number of Components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
return VI * (self._aveF2CCV * (self._MfMui * self._bSecondary(eSolution, srcList)))
def _hDeriv_u(self, src, du_dm_v, adjoint = False):
n = int(self._aveF2CCV.shape[0] / self._nC) # Number of Components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
if adjoint:
v = self._MfMui.T * (self._aveF2CCV.T * (VI.T * du_dm_v))
return self._bDeriv_u(src, v, adjoint=adjoint)
return VI * (self._aveF2CCV * (self._MfMui * self._bDeriv_u(src, du_dm_v, adjoint = adjoint)))
def _hDeriv_m(self, src, v, adjoint = False):
n = int(self._aveF2CCV.shape[0] / self._nC) # Number of Components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
if adjoint:
v = self._MfMui.T * (self._aveF2CCV.T * (VI.T * v))
return self._bDeriv_m(src, v, adjoint=adjoint)
return VI * (self._aveF2CCV * (self._MfMui * self._bDeriv_m(src, v, adjoint = adjoint)))
@@ -617,7 +609,7 @@ class Fields_b(Fields):
n = int(aveE2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy
# Sigma = sdiag(np.kron(np.ones(n), sigma))
Sigma = self.prob.MeSigma
VI = sdiag(1./np.kron(np.ones(n), self.prob.mesh.vol))
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
e = self._e(bSolution, srcList)
@@ -629,7 +621,7 @@ class Fields_b(Fields):
aveF2CCV = self._aveF2CCV
n = int(aveF2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy
# Mui = sdiag(sp.kron(np.ones(n), mui))
VI = sdiag(1./np.kron(np.ones(n), self.prob.mesh.vol))
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
return VI * (aveF2CCV * (Mui * b))
@@ -868,7 +860,7 @@ class Fields_j(Fields):
n = int(aveF2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy
Rho = self.prob.MfRho
VI = sdiag(1./np.kron(np.ones(n), self.prob.mesh.vol))
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
j = self._j(jSolution, srcList)
@@ -887,7 +879,7 @@ class Fields_j(Fields):
aveE2CCV = self._aveE2CCV
n = int(aveE2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy
# Mu = sdiag(sp.kron(np.ones(n), mu))
VI = sdiag(1./np.kron(np.ones(n), self.prob.mesh.vol))
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
return VI * (aveE2CCV * (Mu * h))
@@ -1088,7 +1080,7 @@ class Fields_h(Fields):
n = int(aveF2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy
Rho = self.prob.MfRho
VI = sdiag(1./np.kron(np.ones(n), self.prob.mesh.vol))
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
j = self._j(hSolution, srcList)
@@ -1107,6 +1099,6 @@ class Fields_h(Fields):
aveE2CCV = self._aveE2CCV
n = int(aveE2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy
# Mu = sdiag(sp.kron(np.ones(n), mu))
VI = sdiag(1./np.kron(np.ones(n), self.prob.mesh.vol))
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
return VI * (aveE2CCV * (Mu * h))
+1 -1
View File
@@ -104,7 +104,7 @@ class Rx(SimPEG.Survey.BaseRx):
# print self.knownRxTypes[self.rxType][:2], 'Deriv', projGLoc
# projGLoc += self.knownRxTypes[self.rxType][1]
P = self.getP(mesh)
P = self.getP(mesh, self.projGLoc(u))
if not adjoint:
Pv_complex = P * v
@@ -46,57 +46,57 @@ def adjointTest(fdemType, comp):
class FDEM_AdjointTests(unittest.TestCase):
if testE:
# def test_Jtvec_adjointTest_exr_Eform(self):
# self.assertTrue(adjointTest('e', 'exr'))
# def test_Jtvec_adjointTest_eyr_Eform(self):
# self.assertTrue(adjointTest('e', 'eyr'))
# def test_Jtvec_adjointTest_ezr_Eform(self):
# self.assertTrue(adjointTest('e', 'ezr'))
# def test_Jtvec_adjointTest_exi_Eform(self):
# self.assertTrue(adjointTest('e', 'exi'))
# def test_Jtvec_adjointTest_eyi_Eform(self):
# self.assertTrue(adjointTest('e', 'eyi'))
# def test_Jtvec_adjointTest_ezi_Eform(self):
# self.assertTrue(adjointTest('e', 'ezi'))
def test_Jtvec_adjointTest_exr_Eform(self):
self.assertTrue(adjointTest('e', 'exr'))
def test_Jtvec_adjointTest_eyr_Eform(self):
self.assertTrue(adjointTest('e', 'eyr'))
def test_Jtvec_adjointTest_ezr_Eform(self):
self.assertTrue(adjointTest('e', 'ezr'))
def test_Jtvec_adjointTest_exi_Eform(self):
self.assertTrue(adjointTest('e', 'exi'))
def test_Jtvec_adjointTest_eyi_Eform(self):
self.assertTrue(adjointTest('e', 'eyi'))
def test_Jtvec_adjointTest_ezi_Eform(self):
self.assertTrue(adjointTest('e', 'ezi'))
# def test_Jtvec_adjointTest_bxr_Eform(self):
# self.assertTrue(adjointTest('e', 'bxr'))
# def test_Jtvec_adjointTest_byr_Eform(self):
# self.assertTrue(adjointTest('e', 'byr'))
# def test_Jtvec_adjointTest_bzr_Eform(self):
# self.assertTrue(adjointTest('e', 'bzr'))
# def test_Jtvec_adjointTest_bxi_Eform(self):
# self.assertTrue(adjointTest('e', 'bxi'))
# def test_Jtvec_adjointTest_byi_Eform(self):
# self.assertTrue(adjointTest('e', 'byi'))
# def test_Jtvec_adjointTest_bzi_Eform(self):
# self.assertTrue(adjointTest('e', 'bzi'))
def test_Jtvec_adjointTest_bxr_Eform(self):
self.assertTrue(adjointTest('e', 'bxr'))
def test_Jtvec_adjointTest_byr_Eform(self):
self.assertTrue(adjointTest('e', 'byr'))
def test_Jtvec_adjointTest_bzr_Eform(self):
self.assertTrue(adjointTest('e', 'bzr'))
def test_Jtvec_adjointTest_bxi_Eform(self):
self.assertTrue(adjointTest('e', 'bxi'))
def test_Jtvec_adjointTest_byi_Eform(self):
self.assertTrue(adjointTest('e', 'byi'))
def test_Jtvec_adjointTest_bzi_Eform(self):
self.assertTrue(adjointTest('e', 'bzi'))
def test_Jtvec_adjointTest_exr_Eform(self):
self.assertTrue(adjointTest('e', 'jxr'))
# def test_Jtvec_adjointTest_eyr_Eform(self):
# self.assertTrue(adjointTest('e', 'jyr'))
# def test_Jtvec_adjointTest_ezr_Eform(self):
# self.assertTrue(adjointTest('e', 'jzr'))
# def test_Jtvec_adjointTest_exi_Eform(self):
# self.assertTrue(adjointTest('e', 'jxi'))
# def test_Jtvec_adjointTest_eyi_Eform(self):
# self.assertTrue(adjointTest('e', 'jyi'))
# def test_Jtvec_adjointTest_ezi_Eform(self):
# self.assertTrue(adjointTest('e', 'jzi'))
def test_Jtvec_adjointTest_eyr_Eform(self):
self.assertTrue(adjointTest('e', 'jyr'))
def test_Jtvec_adjointTest_ezr_Eform(self):
self.assertTrue(adjointTest('e', 'jzr'))
def test_Jtvec_adjointTest_exi_Eform(self):
self.assertTrue(adjointTest('e', 'jxi'))
def test_Jtvec_adjointTest_eyi_Eform(self):
self.assertTrue(adjointTest('e', 'jyi'))
def test_Jtvec_adjointTest_ezi_Eform(self):
self.assertTrue(adjointTest('e', 'jzi'))
# def test_Jtvec_adjointTest_bxr_Eform(self):
# self.assertTrue(adjointTest('e', 'hxr'))
# def test_Jtvec_adjointTest_byr_Eform(self):
# self.assertTrue(adjointTest('e', 'hyr'))
# def test_Jtvec_adjointTest_bzr_Eform(self):
# self.assertTrue(adjointTest('e', 'hzr'))
# def test_Jtvec_adjointTest_bxi_Eform(self):
# self.assertTrue(adjointTest('e', 'hxi'))
# def test_Jtvec_adjointTest_byi_Eform(self):
# self.assertTrue(adjointTest('e', 'hyi'))
# def test_Jtvec_adjointTest_bzi_Eform(self):
# self.assertTrue(adjointTest('e', 'hzi'))
def test_Jtvec_adjointTest_bxr_Eform(self):
self.assertTrue(adjointTest('e', 'hxr'))
def test_Jtvec_adjointTest_byr_Eform(self):
self.assertTrue(adjointTest('e', 'hyr'))
def test_Jtvec_adjointTest_bzr_Eform(self):
self.assertTrue(adjointTest('e', 'hzr'))
def test_Jtvec_adjointTest_bxi_Eform(self):
self.assertTrue(adjointTest('e', 'hxi'))
def test_Jtvec_adjointTest_byi_Eform(self):
self.assertTrue(adjointTest('e', 'hyi'))
def test_Jtvec_adjointTest_bzi_Eform(self):
self.assertTrue(adjointTest('e', 'hzi'))
if testB:
def test_Jtvec_adjointTest_exr_Bform(self):
@@ -84,45 +84,45 @@ class FDEM_DerivTests(unittest.TestCase):
def test_Jvec_ezi_Eform(self):
self.assertTrue(derivTest('e', 'jzi'))
# def test_Jvec_bxr_Eform(self):
# self.assertTrue(derivTest('e', 'hxr'))
# def test_Jvec_byr_Eform(self):
# self.assertTrue(derivTest('e', 'hyr'))
# def test_Jvec_bzr_Eform(self):
# self.assertTrue(derivTest('e', 'hzr'))
# def test_Jvec_bxi_Eform(self):
# self.assertTrue(derivTest('e', 'hxi'))
# def test_Jvec_byi_Eform(self):
# self.assertTrue(derivTest('e', 'hyi'))
# def test_Jvec_bzi_Eform(self):
# self.assertTrue(derivTest('e', 'hzi'))
def test_Jvec_bxr_Eform(self):
self.assertTrue(derivTest('e', 'hxr'))
def test_Jvec_byr_Eform(self):
self.assertTrue(derivTest('e', 'hyr'))
def test_Jvec_bzr_Eform(self):
self.assertTrue(derivTest('e', 'hzr'))
def test_Jvec_bxi_Eform(self):
self.assertTrue(derivTest('e', 'hxi'))
def test_Jvec_byi_Eform(self):
self.assertTrue(derivTest('e', 'hyi'))
def test_Jvec_bzi_Eform(self):
self.assertTrue(derivTest('e', 'hzi'))
if testB:
def test_Jvec_exr_Bform(self):
self.assertTrue(derivTest('b', 'exr'))
# def test_Jvec_eyr_Bform(self):
# self.assertTrue(derivTest('b', 'eyr'))
# def test_Jvec_ezr_Bform(self):
# self.assertTrue(derivTest('b', 'ezr'))
# def test_Jvec_exi_Bform(self):
# self.assertTrue(derivTest('b', 'exi'))
# def test_Jvec_eyi_Bform(self):
# self.assertTrue(derivTest('b', 'eyi'))
# def test_Jvec_ezi_Bform(self):
# self.assertTrue(derivTest('b', 'ezi'))
def test_Jvec_eyr_Bform(self):
self.assertTrue(derivTest('b', 'eyr'))
def test_Jvec_ezr_Bform(self):
self.assertTrue(derivTest('b', 'ezr'))
def test_Jvec_exi_Bform(self):
self.assertTrue(derivTest('b', 'exi'))
def test_Jvec_eyi_Bform(self):
self.assertTrue(derivTest('b', 'eyi'))
def test_Jvec_ezi_Bform(self):
self.assertTrue(derivTest('b', 'ezi'))
# def test_Jvec_bxr_Bform(self):
# self.assertTrue(derivTest('b', 'bxr'))
# def test_Jvec_byr_Bform(self):
# self.assertTrue(derivTest('b', 'byr'))
# def test_Jvec_bzr_Bform(self):
# self.assertTrue(derivTest('b', 'bzr'))
# def test_Jvec_bxi_Bform(self):
# self.assertTrue(derivTest('b', 'bxi'))
# def test_Jvec_byi_Bform(self):
# self.assertTrue(derivTest('b', 'byi'))
# def test_Jvec_bzi_Bform(self):
# self.assertTrue(derivTest('b', 'bzi'))
def test_Jvec_bxr_Bform(self):
self.assertTrue(derivTest('b', 'bxr'))
def test_Jvec_byr_Bform(self):
self.assertTrue(derivTest('b', 'byr'))
def test_Jvec_bzr_Bform(self):
self.assertTrue(derivTest('b', 'bzr'))
def test_Jvec_bxi_Bform(self):
self.assertTrue(derivTest('b', 'bxi'))
def test_Jvec_byi_Bform(self):
self.assertTrue(derivTest('b', 'byi'))
def test_Jvec_bzi_Bform(self):
self.assertTrue(derivTest('b', 'bzi'))
if testHJ:
def test_Jvec_jxr_Jform(self):