From 8d4e001301ca7a97db6d779a518c9e4ccc87f14d Mon Sep 17 00:00:00 2001 From: Lindsey Date: Thu, 26 Feb 2015 13:51:42 -0800 Subject: [PATCH] HJ formulation solving for h, with h data --- simpegEM/FDEM/FDEM.py | 66 +++++++++++++++++++++++++++---------- simpegEM/Tests/test_FDEM.py | 48 +++++++++++++-------------- 2 files changed, 73 insertions(+), 41 deletions(-) diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py index 5343c45d..56edd731 100644 --- a/simpegEM/FDEM/FDEM.py +++ b/simpegEM/FDEM/FDEM.py @@ -538,7 +538,7 @@ class ProblemFDEM_h(BaseFDEMProblem): if adjoint: 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): @@ -606,28 +606,60 @@ class ProblemFDEM_h(BaseFDEMProblem): h = sol if fieldType == 'j': NotImplementedError('fieldType "%s" is not implemented.' % fieldType) - C = self.mesh.edgeCurl - j_s = self.getjs(freq) - if adjoint: - NotImplementedError('fieldType "%s" is not implemented.' % fieldType) - # return C.T*h # TODO: THIS IS WRONG - NotImplementedError('fieldType "%s" is not implemented.' % fieldType) - return C*h # - j_s + # C = self.mesh.edgeCurl + # j_s = self.getjs(freq) + # if adjoint: + # MeMu = self.MeMu + # MfSigi = self.MfSigmai + # MfSigmaiinv = self.Solver(MfSigi.T, **self.solverOpts) + # 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': return h raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) def calcFieldsDeriv(self, sol, freq, fieldType, v, adjoint=False): h = sol + A = self.getA(freq) + if fieldType == 'j': - NotImplementedError('fieldType "%s" is not implemented.' % fieldType) - C = self.mesh.edgeCurl - drhs = self.getRHSDeriv(freq,v,adjoint) - if adjoint: - NotImplementedError('fieldType "%s" is not implemented.' % fieldType) - # return -C.T * drhs - NotImplementedError('fieldType "%s" is not implemented.' % fieldType) - # return -C * drhs + # C = self.mesh.edgeCurl + # MeMu = self.MeMu + # MfSigi = self.MfSigmai + # MfSigmaiinv = self.Solver(MfSigi, **self.solverOpts) + # CTinv = self.Solver(C.T, **self.solverOpts) + # hDeriv = self.calcFieldsDeriv(h,freq,'h',v,adjoint) + + # 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': - 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) \ No newline at end of file diff --git a/simpegEM/Tests/test_FDEM.py b/simpegEM/Tests/test_FDEM.py index 7894dbab..453cd756 100644 --- a/simpegEM/Tests/test_FDEM.py +++ b/simpegEM/Tests/test_FDEM.py @@ -224,18 +224,18 @@ class FDEM_DerivTests(unittest.TestCase): def test_Jvec_hzi_Jform(self): self.assertTrue(derivTest('j', 'hzi')) - # def test_Jvec_hxr_Hform(self): - # self.assertTrue(derivTest('h', 'hxr')) - # def test_Jvec_hyr_Hform(self): - # self.assertTrue(derivTest('h', 'hyr')) - # def test_Jvec_hzr_Hform(self): - # self.assertTrue(derivTest('h', 'hzr')) - # def test_Jvec_hxi_Hform(self): - # self.assertTrue(derivTest('h', 'hxi')) - # def test_Jvec_hyi_Hform(self): - # self.assertTrue(derivTest('h', 'hyi')) - # def test_Jvec_hzi_Hform(self): - # self.assertTrue(derivTest('h', 'hzi')) + def test_Jvec_hxr_Hform(self): + self.assertTrue(derivTest('h', 'hxr')) + def test_Jvec_hyr_Hform(self): + self.assertTrue(derivTest('h', 'hyr')) + def test_Jvec_hzr_Hform(self): + self.assertTrue(derivTest('h', 'hzr')) + def test_Jvec_hxi_Hform(self): + self.assertTrue(derivTest('h', 'hxi')) + def test_Jvec_hyi_Hform(self): + self.assertTrue(derivTest('h', 'hyi')) + def test_Jvec_hzi_Hform(self): + self.assertTrue(derivTest('h', 'hzi')) # def test_Jvec_hxr_Hform(self): # self.assertTrue(derivTest('h', 'jxr')) @@ -276,18 +276,18 @@ class FDEM_DerivTests(unittest.TestCase): def test_Jtvec_adjointTest_hzi_Jform(self): self.assertTrue(adjointTest('j', 'hzi')) - # def test_Jtvec_adjointTest_hxr_Hform(self): - # self.assertTrue(adjointTest('h', 'hxr')) - # def test_Jtvec_adjointTest_hyr_Hform(self): - # self.assertTrue(adjointTest('h', 'hyr')) - # def test_Jtvec_adjointTest_hzr_Hform(self): - # self.assertTrue(adjointTest('h', 'hzr')) - # def test_Jtvec_adjointTest_hxi_Hform(self): - # self.assertTrue(adjointTest('h', 'hxi')) - # def test_Jtvec_adjointTest_hyi_Hform(self): - # self.assertTrue(adjointTest('h', 'hyi')) - # def test_Jtvec_adjointTest_hzi_Hform(self): - # self.assertTrue(adjointTest('h', 'hzi')) + def test_Jtvec_adjointTest_hxr_Hform(self): + self.assertTrue(adjointTest('h', 'hxr')) + def test_Jtvec_adjointTest_hyr_Hform(self): + self.assertTrue(adjointTest('h', 'hyr')) + def test_Jtvec_adjointTest_hzr_Hform(self): + self.assertTrue(adjointTest('h', 'hzr')) + def test_Jtvec_adjointTest_hxi_Hform(self): + self.assertTrue(adjointTest('h', 'hxi')) + def test_Jtvec_adjointTest_hyi_Hform(self): + self.assertTrue(adjointTest('h', 'hyi')) + def test_Jtvec_adjointTest_hzi_Hform(self): + self.assertTrue(adjointTest('h', 'hzi')) # def test_Jtvec_adjointTest_hxr_Hform(self): # self.assertTrue(adjointTest('h', 'jxr'))