wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-09-12 12:51:28 +08:00
Code Issues Packages Projects Releases Wiki Activity

HJ formulation solving for h, with h data

Browse Source
This commit is contained in:
Lindsey
2015-02-26 13:51:42 -08:00
parent 5e323591b9
commit 8d4e001301
2 changed files with 73 additions and 41 deletions
Download Patch File Download Diff File Expand all files Collapse all files
simpegEM/FDEM/FDEM.py
+49 -17
View File
@@ -538,7 +538,7 @@ class ProblemFDEM_h(BaseFDEMProblem):
if adjoint: if adjoint:
return (dsig_dm.T * (dsigi_dsig.T * (dMf_dsigi.T * (C * v)))) return (dsig_dm.T * (dsigi_dsig.T * (dMf_dsigi.T * (C * v))))
return (C.T * (dMf_dsigi * (dsigi_dsig * (dsig_dm * v)))) return (C.T * (dMf_dsigi * (dsigi_dsig * (dsig_dm * v))))
def getjs(self,freq): def getjs(self,freq):
@@ -606,28 +606,60 @@ class ProblemFDEM_h(BaseFDEMProblem):
h = sol h = sol
if fieldType == 'j': if fieldType == 'j':
NotImplementedError('fieldType "%s" is not implemented.' % fieldType) NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
C = self.mesh.edgeCurl # C = self.mesh.edgeCurl
j_s = self.getjs(freq) # j_s = self.getjs(freq)
if adjoint: # if adjoint:
NotImplementedError('fieldType "%s" is not implemented.' % fieldType) # MeMu = self.MeMu
# return C.T*h # TODO: THIS IS WRONG # MfSigi = self.MfSigmai
NotImplementedError('fieldType "%s" is not implemented.' % fieldType) # MfSigmaiinv = self.Solver(MfSigi.T, **self.solverOpts)
return C*h # - j_s # Cinv = self.Solver(C, **self.solverOpts)
# return -1j * omega(freq) * MeMu.T * (MfSigmaiinv * (CTinv * h))
# return C*h - j_s # -iomega(freq) inv(MfSigmai) inv(C.T) MeMu
elif fieldType == 'h': elif fieldType == 'h':
return h return h
raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
def calcFieldsDeriv(self, sol, freq, fieldType, v, adjoint=False): def calcFieldsDeriv(self, sol, freq, fieldType, v, adjoint=False):
h = sol h = sol
A = self.getA(freq)
if fieldType == 'j': if fieldType == 'j':
NotImplementedError('fieldType "%s" is not implemented.' % fieldType) # C = self.mesh.edgeCurl
C = self.mesh.edgeCurl # MeMu = self.MeMu
drhs = self.getRHSDeriv(freq,v,adjoint) # MfSigi = self.MfSigmai
if adjoint: # MfSigmaiinv = self.Solver(MfSigi, **self.solverOpts)
NotImplementedError('fieldType "%s" is not implemented.' % fieldType) # CTinv = self.Solver(C.T, **self.solverOpts)
# return -C.T * drhs # hDeriv = self.calcFieldsDeriv(h,freq,'h',v,adjoint)
NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
# return -C * drhs # pt1 = -1j * omega(freq) * (MfSigmaiinv * (CTinv * (MeMu * hDeriv)))
# sig = self.curModel.transform
# sigi = 1/sig
# dsig_dm = self.curModel.transformDeriv
# dsigi_dsig = -Utils.sdiag(sigi)**2
# v1 = MfSigmaiinv *(CTinv * (MeMu * h))
# dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(v1)
# pt2 = 1j * omega(freq) * (MfSigmaiinv * (dMf_dsigi * v1))
# return pt1+pt2
# if adjoint:
# MfSigmaiTinv = self.Solver(MfSigi.T, **self.solverOpts)
# Cinv = self.Solver(C, **self.solverOpts)
# MeMu.T * (Cinv * (MfSigmaiTinv * h))
# pt1 = -1j * omega(freq) * hDeriv.T * MeMu.T * Cinv *(MfSigmaiTinv * v)
raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
elif fieldType == 'h': elif fieldType == 'h':
return -self.getRHSDeriv(freq,v,adjoint) if adjoint:
ATinv = self.Solver(A.T, **self.solverOpts)
ATinvv = ATinv*v
return self.getRHSDeriv(freq,ATinvv,adjoint=True)
dRHSh = self.getRHSDeriv(freq,v,adjoint)
Ainv = self.Solver(A, **self.solverOpts)
return Ainv*dRHSh
raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
simpegEM/Tests/test_FDEM.py
+24 -24
View File
@@ -224,18 +224,18 @@ class FDEM_DerivTests(unittest.TestCase):
def test_Jvec_hzi_Jform(self): def test_Jvec_hzi_Jform(self):
self.assertTrue(derivTest('j', 'hzi')) self.assertTrue(derivTest('j', 'hzi'))
# def test_Jvec_hxr_Hform(self): def test_Jvec_hxr_Hform(self):
# self.assertTrue(derivTest('h', 'hxr')) self.assertTrue(derivTest('h', 'hxr'))
# def test_Jvec_hyr_Hform(self): def test_Jvec_hyr_Hform(self):
# self.assertTrue(derivTest('h', 'hyr')) self.assertTrue(derivTest('h', 'hyr'))
# def test_Jvec_hzr_Hform(self): def test_Jvec_hzr_Hform(self):
# self.assertTrue(derivTest('h', 'hzr')) self.assertTrue(derivTest('h', 'hzr'))
# def test_Jvec_hxi_Hform(self): def test_Jvec_hxi_Hform(self):
# self.assertTrue(derivTest('h', 'hxi')) self.assertTrue(derivTest('h', 'hxi'))
# def test_Jvec_hyi_Hform(self): def test_Jvec_hyi_Hform(self):
# self.assertTrue(derivTest('h', 'hyi')) self.assertTrue(derivTest('h', 'hyi'))
# def test_Jvec_hzi_Hform(self): def test_Jvec_hzi_Hform(self):
# self.assertTrue(derivTest('h', 'hzi')) self.assertTrue(derivTest('h', 'hzi'))
# def test_Jvec_hxr_Hform(self): # def test_Jvec_hxr_Hform(self):
# self.assertTrue(derivTest('h', 'jxr')) # self.assertTrue(derivTest('h', 'jxr'))
@@ -276,18 +276,18 @@ class FDEM_DerivTests(unittest.TestCase):
def test_Jtvec_adjointTest_hzi_Jform(self): def test_Jtvec_adjointTest_hzi_Jform(self):
self.assertTrue(adjointTest('j', 'hzi')) self.assertTrue(adjointTest('j', 'hzi'))
# def test_Jtvec_adjointTest_hxr_Hform(self): def test_Jtvec_adjointTest_hxr_Hform(self):
# self.assertTrue(adjointTest('h', 'hxr')) self.assertTrue(adjointTest('h', 'hxr'))
# def test_Jtvec_adjointTest_hyr_Hform(self): def test_Jtvec_adjointTest_hyr_Hform(self):
# self.assertTrue(adjointTest('h', 'hyr')) self.assertTrue(adjointTest('h', 'hyr'))
# def test_Jtvec_adjointTest_hzr_Hform(self): def test_Jtvec_adjointTest_hzr_Hform(self):
# self.assertTrue(adjointTest('h', 'hzr')) self.assertTrue(adjointTest('h', 'hzr'))
# def test_Jtvec_adjointTest_hxi_Hform(self): def test_Jtvec_adjointTest_hxi_Hform(self):
# self.assertTrue(adjointTest('h', 'hxi')) self.assertTrue(adjointTest('h', 'hxi'))
# def test_Jtvec_adjointTest_hyi_Hform(self): def test_Jtvec_adjointTest_hyi_Hform(self):
# self.assertTrue(adjointTest('h', 'hyi')) self.assertTrue(adjointTest('h', 'hyi'))
# def test_Jtvec_adjointTest_hzi_Hform(self): def test_Jtvec_adjointTest_hzi_Hform(self):
# self.assertTrue(adjointTest('h', 'hzi')) self.assertTrue(adjointTest('h', 'hzi'))
# def test_Jtvec_adjointTest_hxr_Hform(self): # def test_Jtvec_adjointTest_hxr_Hform(self):
# self.assertTrue(adjointTest('h', 'jxr')) # self.assertTrue(adjointTest('h', 'jxr'))