Tested and passing HJ formulation, solving for J. Only works with diagonal anisotropy

This commit is contained in:
Lindsey
2015-02-25 16:28:54 -08:00
parent 2a6f0ab440
commit 57380ea2d8
3 changed files with 189 additions and 168 deletions
+1 -1
View File
@@ -7,7 +7,7 @@ class BaseEMProblem(Problem.BaseProblem):
Problem.BaseProblem.__init__(self, mesh, **kwargs)
solType = None
storeTheseFields = ['e', 'b', 'j', 'h']
storeTheseFields = ['e', 'b']
surveyPair = Survey.BaseSurvey
dataPair = Survey.Data
+29 -8
View File
@@ -105,6 +105,9 @@ class BaseFDEMProblem(BaseEMProblem):
return Jtv
##########################################################################################
################################ E-B Formulation #########################################
##########################################################################################
class ProblemFDEM_e(BaseFDEMProblem):
"""
@@ -335,6 +338,10 @@ class ProblemFDEM_b(BaseFDEMProblem):
return None
raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
##########################################################################################
################################ H-J Formulation #########################################
##########################################################################################
class ProblemFDEM_j(BaseFDEMProblem):
"""
@@ -360,6 +367,7 @@ class ProblemFDEM_j(BaseFDEMProblem):
def __init__(self, model, **kwargs):
BaseFDEMProblem.__init__(self, model, **kwargs)
BaseFDEMProblem.storeTheseFields = ['j','h']
def getA(self, freq):
"""
@@ -377,17 +385,18 @@ class ProblemFDEM_j(BaseFDEMProblem):
def getADeriv(self, freq, u, v, adjoint=False):
mui = self.MeMui
MeMui = self.MeMui
C = self.mesh.edgeCurl
sig = self.curModel.transform
sigi = 1/sig
dsig_dm = self.curModel.transformDeriv
dMf_dsig = self.mesh.getFaceInnerProductDeriv(sig)(u)
dMfSigi_di = - self.MfSigmai**2
dsigi_dsig = -Utils.sdiag(sigi)**2
dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(u)
if adjoint:
return dsig_dm.T * ( dMf_dsig.T * ( dMfSigi_di.T * ( C * ( mui.T * ( C.T * v ) ) ) ) )
return dsig_dm.T * ( dsigi_dsig.T *( dMf_dsigi.T * ( C * ( MeMui.T * ( C.T * v ) ) ) ) )
return C * ( mui * ( C.T * ( dMfSigi_di * ( dMf_dsig * ( dsig_dm * v ) ) ) ) )
return C * ( MeMui * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) ) )
def getRHS(self, freq):
@@ -428,10 +437,11 @@ class ProblemFDEM_j(BaseFDEMProblem):
elif fieldType == 'h':
mui = self.MeMui
C = self.mesh.edgeCurl
MfSigi = self.MfSigmai
if not adjoint:
h = -(1./(1j*omega(freq))) * mui * ( C.T * j )
h = -(1./(1j*omega(freq))) * mui * ( C.T * ( MfSigi * j ) )
else:
h = -(1./(1j*omega(freq))) * C * ( mui.T * j )
h = -(1./(1j*omega(freq))) * MfSigi * ( C * ( mui.T * j ) )
return h
raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
@@ -440,5 +450,16 @@ class ProblemFDEM_j(BaseFDEMProblem):
if fieldType == 'j':
return None
elif fieldType == 'h':
return None
MeMui = self.MeMui
C = self.mesh.edgeCurl
sig = self.curModel.transform
sigi = 1/sig
dsig_dm = self.curModel.transformDeriv
dsigi_dsig = -Utils.sdiag(sigi)**2
dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(j)
sigi = self.MfSigmai
if not adjoint:
return -(1./(1j*omega(freq))) * MeMui * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) )
else:
return -(1./(1j*omega(freq))) * dsig_dm.T * ( dsigi_dsig.T * ( dMf_dsigi.T * ( C * ( MeMui.T * v ) ) ) )
raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
+159 -159
View File
@@ -5,9 +5,10 @@ import sys
from scipy.constants import mu_0
TOL = 1e-4
CONDUCTIVITY = 1e3
FLR = 1e-15 # "zero", so if residual below this --> pass regardless of order
CONDUCTIVITY = 1e1
MU = mu_0
addrandoms = True
addrandoms = True # important to addrandoms if testing HJ formulation with VMD source! (or else jz ~ 0)
def getProblem(fdemType, comp):
cs = 5.
@@ -23,7 +24,7 @@ def getProblem(fdemType, comp):
x = np.linspace(-30,30,6)
XYZ = Utils.ndgrid(x,x,np.r_[0])
Rx0 = EM.FDEM.RxFDEM(XYZ, comp)
Tx0 = EM.FDEM.TxFDEM(np.r_[4.,2.,2.], 'VMD', 1e-2, [Rx0])
Tx0 = EM.FDEM.TxFDEM(np.r_[0.,0.,0.], 'VMD', 1, [Rx0])
survey = EM.FDEM.SurveyFDEM([Tx0])
@@ -39,11 +40,11 @@ def getProblem(fdemType, comp):
raise NotImplementedError()
prb.pair(survey)
try:
from pymatsolver import MumpsSolver
prb.Solver = MumpsSolver
except ImportError, e:
pass
# try:
# from pymatsolver import MumpsSolver
# prb.Solver = MumpsSolver
# except ImportError, e:
# pass
return prb
@@ -55,19 +56,20 @@ def adjointTest(fdemType, comp):
mu = np.log(np.ones(prb.mesh.nC)*MU)
if addrandoms is True:
m = m + np.random.randn(prb.mesh.nC)*CONDUCTIVITY*1e-3
mu = mu + np.random.randn(prb.mesh.nC)*MU*1e-3
m = m + np.random.randn(prb.mesh.nC)*CONDUCTIVITY*1e-1
mu = mu + np.random.randn(prb.mesh.nC)*MU*1e-1
prb.mu = mu
survey = prb.survey
u = prb.fields(m)
v = np.random.rand(survey.nD)
w = np.random.rand(prb.mapping.nP)
u = prb.fields(m)
vJw = v.dot(prb.Jvec(m, w, u=u))
wJtv = w.dot(prb.Jtvec(m, v, u=u))
tol = TOL*(10**int(np.log10(np.abs(vJw))))
tol = np.max([TOL*(10**int(np.log10(np.abs(vJw)))),FLR])
print vJw, wJtv, vJw - wJtv, tol, np.abs(vJw - wJtv) < tol
return np.abs(vJw - wJtv) < tol
@@ -78,174 +80,172 @@ def derivTest(fdemType, comp):
mu = np.log(np.ones(prb.mesh.nC)*MU)
if addrandoms is True:
x0 = x0 + np.random.randn(prb.mesh.nC)*CONDUCTIVITY*1e-3
mu = mu + np.random.randn(prb.mesh.nC)*MU*1e-3
x0 = x0 + np.random.randn(prb.mesh.nC)*CONDUCTIVITY*1e-1
mu = mu + np.random.randn(prb.mesh.nC)*MU*1e-1
prb.mu = mu
survey = prb.survey
def fun(x):
return survey.dpred(x), lambda x: prb.Jvec(x0, x)
return Tests.checkDerivative(fun, x0, num=3, plotIt=False, eps=1e-25)
return Tests.checkDerivative(fun, x0, num=3, plotIt=False, eps=FLR)
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_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_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_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'))
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'))
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__':