Start of HJ formulation: FAILING TESTS right now

This commit is contained in:
Lindsey
2015-02-24 18:00:34 -08:00
parent 4e69ee1f1f
commit 2a6f0ab440
4 changed files with 234 additions and 121 deletions
+13 -10
View File
@@ -7,7 +7,7 @@ class BaseEMProblem(Problem.BaseProblem):
Problem.BaseProblem.__init__(self, mesh, **kwargs)
solType = None
storeTheseFields = ['e', 'b']
storeTheseFields = ['e', 'b', 'j', 'h']
surveyPair = Survey.BaseSurvey
dataPair = Survey.Data
@@ -30,6 +30,8 @@ class BaseEMProblem(Problem.BaseProblem):
def mu(self, value):
if getattr(self, '_MfMui', None) is not None:
del self._MfMui
if getattr(self, '_MeMui', None) is not None:
del self._MeMui
self._mu = value
@@ -45,11 +47,11 @@ class BaseEMProblem(Problem.BaseProblem):
return self._MfMui
@property
def MeMuI(self):
def MeMui(self):
#TODO: assuming constant mu
if getattr(self, '_MeMuI', None) is None:
self._MeMuI = self.mesh.getEdgeInnerProduct(1/mu_0)
return self._MeMuI
if getattr(self, '_MeMui', None) is None:
self._MeMui = self.mesh.getEdgeInnerProduct(1/mu_0)
return self._MeMui
@property
def Me(self):
@@ -74,14 +76,15 @@ class BaseEMProblem(Problem.BaseProblem):
return self._MeSigmaI
@property
def MfSigmaI(self):
def MfSigmai(self):
#TODO: hardcoded to sigma as the model
if getattr(self, '_MfSigmaI', None) is None:
#TODO: hardcoded to sigma diagonal
if getattr(self, '_MfSigmai', None) is None:
sigma = self.curModel.transform
self._MfSigmaI = self.mesh.getFaceInnerProduct(sigma, invMat=True)
return self._MfSigmaI
self._MfSigmai = self.mesh.getFaceInnerProduct(1/sigma)
return self._MfSigmai
deleteTheseOnModelUpdate = ['_MeSigma', '_MeSigmaI','_MfSigmaI']
deleteTheseOnModelUpdate = ['_MeSigma', '_MeSigmaI','_MfSigmai']
def fields(self, m):
self.curModel = m
+52 -14
View File
@@ -350,6 +350,10 @@ class ProblemFDEM_j(BaseFDEMProblem):
.. math::
\\nabla \\times ( \\mu^{-1} \\nabla \\times \\sigma^{-1} \\vec{J} ) + i\\omega \\vec{J} = - i\\omega\\vec{J_s}
.. note::
This implementation does not yet work with full anisotropy!!
"""
solType = 'j'
@@ -364,43 +368,77 @@ class ProblemFDEM_j(BaseFDEMProblem):
:return: A
"""
muI = self.MeMuI
sigI = self.MfSigmaI
mui = self.MeMui
sigi = self.MfSigmai
C = self.mesh.edgeCurl
iomega = 1j * omega(freq) * sp.eye(self.mesh.nF)
return C * muI * C.T * sigI + 1j * omega(freq)
return C * mui * C.T * sigi + iomega
def getADeriv(self, freq, u, v, adjoint=False):
muI = self.MeMuI
mui = self.MeMui
C = self.mesh.edgeCurl
sig = self.curModel.transform
dsig_dm = self.curModel.transformDeriv
#TODO: This only works if diagonal (no tensors)...
dMfSigmaI_dI = - self.MfSigmaI**2
dMf_dsig = self.mesh.getFaceInnerProductDeriv(sig)(u)
dMfSigi_di = - self.MfSigmai**2
if adjoint:
raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
return dsig_dm.T * ( dMf_dsig.T * ( dMfSigi_di.T * ( C * ( mui.T * ( C.T * v ) ) ) ) )
return C * ( muI * ( C.T * ( dMfSigmaI_dI * ( dMf_dsig * ( dsig_dm * v ) ) ) ) )
return C * ( mui * ( C.T * ( dMfSigi_di * ( dMf_dsig * ( dsig_dm * v ) ) ) ) )
raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
def getRHS(self, freq):
"""
:param float freq: Frequency
:rtype: numpy.ndarray (nE, nTx)
:return: RHS
"""
raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
Txs = self.survey.getTransmitters(freq)
rhs = range(len(Txs))
for i, tx in enumerate(Txs):
if tx.txType == 'VMD':
src = Sources.MagneticDipoleVectorPotential
SRCx = src(tx.loc, self.mesh.gridFx, 'x')
SRCy = src(tx.loc, self.mesh.gridFy, 'y')
SRCz = src(tx.loc, self.mesh.gridFz, 'z')
elif tx.txType == 'CircularLoop':
src = Sources.MagneticLoopVectorPotential
SRCx = src(tx.loc, self.mesh.gridFx, 'x', tx.radius)
SRCy = src(tx.loc, self.mesh.gridFy, 'y', tx.radius)
SRCz = src(tx.loc, self.mesh.gridFz, 'z', tx.radius)
else:
raise NotImplemented('%s txType is not implemented' % tx.txType)
rhs[i] = np.concatenate((SRCx, SRCy, SRCz))
a = np.concatenate(rhs).reshape((self.mesh.nF, len(Txs)), order='F')
mui = self.MeMui
C = self.mesh.edgeCurl
j_s = C*mui*C.T*a
return -1j*omega(freq)*j_s
def calcFields(self, sol, freq, fieldType, adjoint=False):
j = sol
if fieldType == 'j':
return j
elif fieldType == 'h':
mui = self.MeMui
C = self.mesh.edgeCurl
if not adjoint:
h = -(1./(1j*omega(freq))) * mui * ( C.T * j )
else:
h = -(1./(1j*omega(freq))) * C * ( mui.T * j )
return h
raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
def calcFieldsDeriv(self, sol, freq, fieldType, v, adjoint=False):
j = sol
if fieldType == 'j':
return None
elif fieldType == 'h':
return None
raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
+15 -1
View File
@@ -16,6 +16,20 @@ class RxFDEM(Survey.BaseRx):
'bxi':['b', 'Fx', 'imag'],
'byi':['b', 'Fy', 'imag'],
'bzi':['b', 'Fz', 'imag'],
'jxr':['j', 'Fx', 'real'],
'jyr':['j', 'Fy', 'real'],
'jzr':['j', 'Fz', 'real'],
'jxi':['j', 'Fx', 'imag'],
'jyi':['j', 'Fy', 'imag'],
'jzi':['j', 'Fz', 'imag'],
'hxr':['h', 'Ex', 'real'],
'hyr':['h', 'Ey', 'real'],
'hzr':['h', 'Ez', 'real'],
'hxi':['h', 'Ex', 'imag'],
'hyi':['h', 'Ey', 'imag'],
'hzi':['h', 'Ez', 'imag'],
}
radius = None
@@ -84,7 +98,7 @@ class TxFDEM(Survey.BaseTx):
class FieldsFDEM(Problem.Fields):
"""Fancy Field Storage for a FDEM survey."""
knownFields = {'b': 'F', 'e': 'E'}
knownFields = {'b': 'F', 'e': 'E', 'j': 'F', 'h': 'E'}
dtype = complex
+154 -96
View File
@@ -27,10 +27,14 @@ def getProblem(fdemType, comp):
survey = EM.FDEM.SurveyFDEM([Tx0])
print fdemType
if fdemType == 'e':
prb = EM.FDEM.ProblemFDEM_e(mesh, mapping=mapping)
elif fdemType == 'b':
prb = EM.FDEM.ProblemFDEM_b(mesh, mapping=mapping)
elif fdemType == 'j':
prb = EM.FDEM.ProblemFDEM_j(mesh, mapping=mapping)
else:
raise NotImplementedError()
prb.pair(survey)
@@ -86,109 +90,163 @@ def derivTest(fdemType, comp):
class FDEM_DerivTests(unittest.TestCase):
def test_Jvec_exr_Eform(self):
self.assertTrue(derivTest('e', 'exr'))
def test_Jvec_exr_Bform(self):
self.assertTrue(derivTest('b', 'exr'))
def test_Jvec_eyr_Eform(self):
self.assertTrue(derivTest('e', 'eyr'))
def test_Jvec_eyr_Bform(self):
self.assertTrue(derivTest('b', 'eyr'))
def test_Jvec_ezr_Eform(self):
self.assertTrue(derivTest('e', 'ezr'))
def test_Jvec_ezr_Bform(self):
self.assertTrue(derivTest('b', 'ezr'))
def test_Jvec_exi_Eform(self):
self.assertTrue(derivTest('e', 'exi'))
def test_Jvec_exi_Bform(self):
self.assertTrue(derivTest('b', 'exi'))
def test_Jvec_eyi_Eform(self):
self.assertTrue(derivTest('e', 'eyi'))
def test_Jvec_eyi_Bform(self):
self.assertTrue(derivTest('b', 'eyi'))
def test_Jvec_ezi_Eform(self):
self.assertTrue(derivTest('e', 'ezi'))
def test_Jvec_ezi_Bform(self):
self.assertTrue(derivTest('b', 'ezi'))
# def test_Jvec_exr_Eform(self):
# self.assertTrue(derivTest('e', 'exr'))
# def test_Jvec_exr_Bform(self):
# self.assertTrue(derivTest('b', 'exr'))
# def test_Jvec_eyr_Eform(self):
# self.assertTrue(derivTest('e', 'eyr'))
# def test_Jvec_eyr_Bform(self):
# self.assertTrue(derivTest('b', 'eyr'))
# def test_Jvec_ezr_Eform(self):
# self.assertTrue(derivTest('e', 'ezr'))
# def test_Jvec_ezr_Bform(self):
# self.assertTrue(derivTest('b', 'ezr'))
# def test_Jvec_exi_Eform(self):
# self.assertTrue(derivTest('e', 'exi'))
# def test_Jvec_exi_Bform(self):
# self.assertTrue(derivTest('b', 'exi'))
# def test_Jvec_eyi_Eform(self):
# self.assertTrue(derivTest('e', 'eyi'))
# def test_Jvec_eyi_Bform(self):
# self.assertTrue(derivTest('b', 'eyi'))
# def test_Jvec_ezi_Eform(self):
# self.assertTrue(derivTest('e', 'ezi'))
# def test_Jvec_ezi_Bform(self):
# self.assertTrue(derivTest('b', 'ezi'))
def test_Jvec_bxr_Eform(self):
self.assertTrue(derivTest('e', 'bxr'))
def test_Jvec_bxr_Bform(self):
self.assertTrue(derivTest('b', 'bxr'))
def test_Jvec_byr_Eform(self):
self.assertTrue(derivTest('e', 'byr'))
def test_Jvec_byr_Bform(self):
self.assertTrue(derivTest('b', 'byr'))
def test_Jvec_bzr_Eform(self):
self.assertTrue(derivTest('e', 'bzr'))
def test_Jvec_bzr_Bform(self):
self.assertTrue(derivTest('b', 'bzr'))
def test_Jvec_bxi_Eform(self):
self.assertTrue(derivTest('e', 'bxi'))
def test_Jvec_bxi_Bform(self):
self.assertTrue(derivTest('b', 'bxi'))
def test_Jvec_byi_Eform(self):
self.assertTrue(derivTest('e', 'byi'))
def test_Jvec_byi_Bform(self):
self.assertTrue(derivTest('b', 'byi'))
def test_Jvec_bzi_Eform(self):
self.assertTrue(derivTest('e', 'bzi'))
def test_Jvec_bzi_Bform(self):
self.assertTrue(derivTest('b', 'bzi'))
# def test_Jvec_bxr_Eform(self):
# self.assertTrue(derivTest('e', 'bxr'))
# def test_Jvec_bxr_Bform(self):
# self.assertTrue(derivTest('b', 'bxr'))
# def test_Jvec_byr_Eform(self):
# self.assertTrue(derivTest('e', 'byr'))
# def test_Jvec_byr_Bform(self):
# self.assertTrue(derivTest('b', 'byr'))
# def test_Jvec_bzr_Eform(self):
# self.assertTrue(derivTest('e', 'bzr'))
# def test_Jvec_bzr_Bform(self):
# self.assertTrue(derivTest('b', 'bzr'))
# def test_Jvec_bxi_Eform(self):
# self.assertTrue(derivTest('e', 'bxi'))
# def test_Jvec_bxi_Bform(self):
# self.assertTrue(derivTest('b', 'bxi'))
# def test_Jvec_byi_Eform(self):
# self.assertTrue(derivTest('e', 'byi'))
# def test_Jvec_byi_Bform(self):
# self.assertTrue(derivTest('b', 'byi'))
# def test_Jvec_bzi_Eform(self):
# self.assertTrue(derivTest('e', 'bzi'))
# def test_Jvec_bzi_Bform(self):
# self.assertTrue(derivTest('b', 'bzi'))
def test_Jtvec_adjointTest_exr_Eform(self):
self.assertTrue(adjointTest('e', 'exr'))
def test_Jtvec_adjointTest_exr_Bform(self):
self.assertTrue(adjointTest('b', 'exr'))
def test_Jtvec_adjointTest_eyr_Eform(self):
self.assertTrue(adjointTest('e', 'eyr'))
def test_Jtvec_adjointTest_eyr_Bform(self):
self.assertTrue(adjointTest('b', 'eyr'))
def test_Jtvec_adjointTest_ezr_Eform(self):
self.assertTrue(adjointTest('e', 'ezr'))
def test_Jtvec_adjointTest_ezr_Bform(self):
self.assertTrue(adjointTest('b', 'ezr'))
def test_Jtvec_adjointTest_exi_Eform(self):
self.assertTrue(adjointTest('e', 'exi'))
def test_Jtvec_adjointTest_exi_Bform(self):
self.assertTrue(adjointTest('b', 'exi'))
def test_Jtvec_adjointTest_eyi_Eform(self):
self.assertTrue(adjointTest('e', 'eyi'))
def test_Jtvec_adjointTest_eyi_Bform(self):
self.assertTrue(adjointTest('b', 'eyi'))
def test_Jtvec_adjointTest_ezi_Eform(self):
self.assertTrue(adjointTest('e', 'ezi'))
def test_Jtvec_adjointTest_ezi_Bform(self):
self.assertTrue(adjointTest('b', 'ezi'))
# def test_Jtvec_adjointTest_exr_Eform(self):
# self.assertTrue(adjointTest('e', 'exr'))
# def test_Jtvec_adjointTest_exr_Bform(self):
# self.assertTrue(adjointTest('b', 'exr'))
# def test_Jtvec_adjointTest_eyr_Eform(self):
# self.assertTrue(adjointTest('e', 'eyr'))
# def test_Jtvec_adjointTest_eyr_Bform(self):
# self.assertTrue(adjointTest('b', 'eyr'))
# def test_Jtvec_adjointTest_ezr_Eform(self):
# self.assertTrue(adjointTest('e', 'ezr'))
# def test_Jtvec_adjointTest_ezr_Bform(self):
# self.assertTrue(adjointTest('b', 'ezr'))
# def test_Jtvec_adjointTest_exi_Eform(self):
# self.assertTrue(adjointTest('e', 'exi'))
# def test_Jtvec_adjointTest_exi_Bform(self):
# self.assertTrue(adjointTest('b', 'exi'))
# def test_Jtvec_adjointTest_eyi_Eform(self):
# self.assertTrue(adjointTest('e', 'eyi'))
# def test_Jtvec_adjointTest_eyi_Bform(self):
# self.assertTrue(adjointTest('b', 'eyi'))
# def test_Jtvec_adjointTest_ezi_Eform(self):
# self.assertTrue(adjointTest('e', 'ezi'))
# def test_Jtvec_adjointTest_ezi_Bform(self):
# self.assertTrue(adjointTest('b', 'ezi'))
def test_Jtvec_adjointTest_bxr_Eform(self):
self.assertTrue(adjointTest('e', 'bxr'))
def test_Jtvec_adjointTest_bxr_Bform(self):
self.assertTrue(adjointTest('b', 'bxr'))
def test_Jtvec_adjointTest_byr_Eform(self):
self.assertTrue(adjointTest('e', 'byr'))
def test_Jtvec_adjointTest_byr_Bform(self):
self.assertTrue(adjointTest('b', 'byr'))
def test_Jtvec_adjointTest_bzr_Eform(self):
self.assertTrue(adjointTest('e', 'bzr'))
def test_Jtvec_adjointTest_bzr_Bform(self):
self.assertTrue(adjointTest('b', 'bzr'))
def test_Jtvec_adjointTest_bxi_Eform(self):
self.assertTrue(adjointTest('e', 'bxi'))
def test_Jtvec_adjointTest_bxi_Bform(self):
self.assertTrue(adjointTest('b', 'bxi'))
def test_Jtvec_adjointTest_byi_Eform(self):
self.assertTrue(adjointTest('e', 'byi'))
def test_Jtvec_adjointTest_byi_Bform(self):
self.assertTrue(adjointTest('b', 'byi'))
def test_Jtvec_adjointTest_bzi_Eform(self):
self.assertTrue(adjointTest('e', 'bzi'))
def test_Jtvec_adjointTest_bzi_Bform(self):
self.assertTrue(adjointTest('b', 'bzi'))
# def test_Jtvec_adjointTest_bxr_Eform(self):
# self.assertTrue(adjointTest('e', 'bxr'))
# def test_Jtvec_adjointTest_bxr_Bform(self):
# self.assertTrue(adjointTest('b', 'bxr'))
# def test_Jtvec_adjointTest_byr_Eform(self):
# self.assertTrue(adjointTest('e', 'byr'))
# def test_Jtvec_adjointTest_byr_Bform(self):
# self.assertTrue(adjointTest('b', 'byr'))
# def test_Jtvec_adjointTest_bzr_Eform(self):
# self.assertTrue(adjointTest('e', 'bzr'))
# def test_Jtvec_adjointTest_bzr_Bform(self):
# self.assertTrue(adjointTest('b', 'bzr'))
# def test_Jtvec_adjointTest_bxi_Eform(self):
# self.assertTrue(adjointTest('e', 'bxi'))
# def test_Jtvec_adjointTest_bxi_Bform(self):
# self.assertTrue(adjointTest('b', 'bxi'))
# def test_Jtvec_adjointTest_byi_Eform(self):
# self.assertTrue(adjointTest('e', 'byi'))
# def test_Jtvec_adjointTest_byi_Bform(self):
# self.assertTrue(adjointTest('b', 'byi'))
# def test_Jtvec_adjointTest_bzi_Eform(self):
# self.assertTrue(adjointTest('e', 'bzi'))
# def test_Jtvec_adjointTest_bzi_Bform(self):
# self.assertTrue(adjointTest('b', 'bzi'))
# def test_Jvec_jxr_Jform(self):
# self.assertTrue(derivTest('j', 'jxr'))
# def test_Jvec_jyr_Jform(self):
# self.assertTrue(derivTest('j', 'jyr'))
# def test_Jvec_jzr_Jform(self):
# self.assertTrue(derivTest('j', 'jzr'))
# def test_Jvec_jxi_Jform(self):
# self.assertTrue(derivTest('j', 'jxi'))
# def test_Jvec_jyi_Jform(self):
# self.assertTrue(derivTest('j', 'jyi'))
# def test_Jvec_jzi_Jform(self):
# self.assertTrue(derivTest('j', 'jzi'))
# def test_Jvec_hxr_Jform(self):
# self.assertTrue(derivTest('j', 'hxr'))
# def test_Jvec_hyr_Jform(self):
# self.assertTrue(derivTest('j', 'hyr'))
# def test_Jvec_hzr_Jform(self):
# self.assertTrue(derivTest('j', 'hzr'))
# def test_Jvec_hxi_Jform(self):
# self.assertTrue(derivTest('j', 'hxi'))
# def test_Jvec_hyi_Jform(self):
# self.assertTrue(derivTest('j', 'hyi'))
# def test_Jvec_hzi_Jform(self):
# self.assertTrue(derivTest('j', 'hzi'))
# def test_Jtvec_adjointTest_jxr_Jform(self):
# self.assertTrue(adjointTest('j', 'jxr'))
# def test_Jtvec_adjointTest_jyr_Jform(self):
# self.assertTrue(adjointTest('j', 'jyr'))
# def test_Jtvec_adjointTest_jzr_Jform(self):
# self.assertTrue(adjointTest('j', 'jzr'))
# def test_Jtvec_adjointTest_jxi_Jform(self):
# self.assertTrue(adjointTest('j', 'jxi'))
# def test_Jtvec_adjointTest_jyi_Jform(self):
# self.assertTrue(adjointTest('j', 'jyi'))
# def test_Jtvec_adjointTest_jzi_Jform(self):
# self.assertTrue(adjointTest('j', 'jzi'))R
def test_Jtvec_adjointTest_hxr_Jform(self):
self.assertTrue(adjointTest('j', 'hxr'))
# def test_Jtvec_adjointTest_hyr_Jform(self):
# self.assertTrue(adjointTest('j', 'hyr'))
# def test_Jtvec_adjointTest_hzr_Jform(self):
# self.assertTrue(adjointTest('j', 'hzr'))
# def test_Jtvec_adjointTest_hxi_Jform(self):
# self.assertTrue(adjointTest('j', 'hxi'))
# def test_Jtvec_adjointTest_hyi_Jform(self):
# self.assertTrue(adjointTest('j', 'hyi'))
# def test_Jtvec_adjointTest_hzi_Jform(self):
# self.assertTrue(adjointTest('j', 'hzi'))
if __name__ == '__main__':
unittest.main()