From cd94fea61d3eb95afef8c158f4df95245ca3c83d Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 2 Jun 2015 19:23:30 -0700 Subject: [PATCH] start of Jvec for Problem_e --- simpegEM/FDEM/FDEM.py | 61 +++++++++++++++++++++++++++---------- simpegEM/FDEM/FieldsFDEM.py | 32 ++++++++++++++----- simpegEM/FDEM/SurveyFDEM.py | 2 +- simpegEM/Tests/test_FDEM.py | 54 ++++++++++++++++---------------- 4 files changed, 98 insertions(+), 51 deletions(-) diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py index d93a3509..83fc9ee5 100644 --- a/simpegEM/FDEM/FDEM.py +++ b/simpegEM/FDEM/FDEM.py @@ -44,21 +44,39 @@ class BaseFDEMProblem(BaseEMProblem): for freq in self.survey.freqs: A = self.getA(freq) - Ainv = self.Solver(A, **self.solverOpts) + dF_duI = self.Solver(A, **self.solverOpts) - for src in self.survey.getSource(freq): - u_src = u[src, self.solType] - w = self.getADeriv(freq, u_src, v) - Ainvw = Ainv * w + for src in self.survey.getSrcByFreq(freq): + u_src = u[src, self._fieldType] + dF_dm = self.getADeriv(freq, u_src, v) + dRHS_dm = self.getRHSDeriv(src, v) + if dRHS_dm is None: + du_dm = dF_duI * ( - dF_dm ) + else: + du_dm = dF_duI * ( - dF_dm + dRHS_dm ) for rx in src.rxList: - fAinvw = self.calcFields(Ainvw, freq, rx.projField) + dAl_duFun = getattr(u, '_%sDeriv_u'%rx.projField, None) + dAl_du = dAl_duFun(src, du_dm, adjoint=False) + if dAl_du is not None: + du_dm = dAl_du + + dAl_dmFun = getattr(u, '_%sDeriv_m'%rx.projField, None) + dAl_dm = dAl_dmFun(src, v, adjoint=False) + if dAl_dm is not None: + du_dm += dAl_dm + P = lambda v: rx.projectFieldsDeriv(src, self.mesh, u, v) - Jv[src, rx] = - P(fAinvw) + Jv[src, rx] = P(du_dm) - df_dm = self.calcFieldsDeriv(u_src, freq, rx.projField, v) - if df_dm is not None: - Jv[src, rx] += P(df_dm) + # fAinvw = self.calcFields(Ainvw, freq, rx.projField) + # P = lambda v: rx.projectFieldsDeriv(src, self.mesh, u, v) + + # Jv[src, rx] = - P(fAinvw) + + # df_dm = self.calcFieldsDeriv(u_src, freq, rx.projField, v) + # if df_dm is not None: + # Jv[src, rx] += P(df_dm) return Utils.mkvc(Jv) @@ -201,9 +219,20 @@ class ProblemFDEM_e(BaseFDEMProblem): return RHS - def getRHSDeriv(self, freq, u, v, adjoint=False): - raise NotImplementedError('getRHSDeriv not implemented yet') - return None + def getRHSDeriv(self, src, v, adjoint=False): + S_mDeriv, S_eDeriv = src.evalDeriv(self, v, adjoint) + if adjoint: + # evalDeriv(MfMui.T* C * v, adjoint = True) + raise Exception('Not implemented') + + if S_mDeriv is not None and S_eDeriv is not None: + return C.T * (MfMui * S_mDeriv) -1j*omega(freq)*S_eDeriv + elif S_mDeriv is not None: + return C.T * (MfMui * S_mDeriv) + elif S_eDeriv is not None: + return -1j*omega(freq)*S_eDeriv + else: + return None class ProblemFDEM_b(BaseFDEMProblem): @@ -274,7 +303,7 @@ class ProblemFDEM_b(BaseFDEMProblem): return RHS - def getRHSDeriv(self, freq, u, v, adjoint=False): + def getRHSDeriv(self, freq, v, adjoint=False): raise NotImplementedError('getRHSDeriv not implemented yet') return None @@ -385,7 +414,7 @@ class ProblemFDEM_j(BaseFDEMProblem): return RHS - def getRHSDeriv(self, freq, u, v, adjoint=False): + def getRHSDeriv(self, freq, v, adjoint=False): raise NotImplementedError('getRHSDeriv not implemented yet') return None @@ -466,7 +495,7 @@ class ProblemFDEM_h(BaseFDEMProblem): return RHS - def getRHSDeriv(self, freq, u, v, adjoint=False): + def getRHSDeriv(self, freq, v, adjoint=False): raise NotImplementedError('getRHSDeriv not implemented yet') return None diff --git a/simpegEM/FDEM/FieldsFDEM.py b/simpegEM/FDEM/FieldsFDEM.py index 9f6b95c7..8844b3ce 100644 --- a/simpegEM/FDEM/FieldsFDEM.py +++ b/simpegEM/FDEM/FieldsFDEM.py @@ -39,6 +39,12 @@ class FieldsFDEM_e(FieldsFDEM): def _e(self, eSolution, srcList): return self._ePrimary(eSolution,srcList) + self._eSecondary(eSolution,srcList) + def _eDeriv_u(self, src, v, adjoint = False): + return None + + def _eDeriv_m(self, src, v, adjoint = False): + return None + def _bPrimary(self, eSolution, srcList): bPrimary = np.zeros([self._edgeCurl.shape[0],eSolution.shape[1]],dtype = complex) for i, src in enumerate(srcList): @@ -57,16 +63,28 @@ class FieldsFDEM_e(FieldsFDEM): b[:,i] += 1./(1j*omega(src.freq)) * S_m return b + def _bSecondaryDeriv_u(self, src, v, adjoint = False): + C = self._edgeCurl + if adjoint: + return - 1./(1j*omega(src.freq)) * (C.T * v) + return - 1./(1j*omega(src.freq)) * (C * v) + + def _bSecondaryDeriv_m(self, src, v, adjoint = False): + S_mDeriv, _ = src.evalDeriv(self.survey.prob, v, adjoint) + if S_mDeriv is not None: + return 1./(1j * omega(src.freq)) * S_mDeriv + return None + def _b(self, eSolution, srcList): return self._bPrimary(eSolution, srcList) + self._bSecondary(eSolution, srcList) - def _bDeriv(self, e, srcList, v, adjoint=False): - raise NotImplementedError('Fields Derivs Not Implemented Yet') - # S_mDeriv,_ = src.getSourceDeriv(self.survey.prob, v, adjoint) - # if S_mDeriv is None: - # return None - # else: - # return 1./(1j*omega(src.freq)) * S_mDeriv + def _bDeriv_u(self, src, v, adjoint=False): + # Primary does not depend on u + return self._bSecondaryDeriv_u(src, v, adjoint) + + def _bDeriv_m(self, src, v, adjoint=False): + # Assuming the primary does not depend on the model + return self._bSecondaryDeriv_m(src, v, adjoint) class FieldsFDEM_b(FieldsFDEM): diff --git a/simpegEM/FDEM/SurveyFDEM.py b/simpegEM/FDEM/SurveyFDEM.py index e2df08c9..84a278c1 100644 --- a/simpegEM/FDEM/SurveyFDEM.py +++ b/simpegEM/FDEM/SurveyFDEM.py @@ -100,7 +100,7 @@ class SrcFDEM(Survey.BaseSrc): S_e = self.S_e(prob) return S_m, S_e - def evalDeriv(self, prob, v, adjoint=None): + def evalDeriv(self, prob, v, adjoint=False): return self.S_mDeriv(prob,v,adjoint), self.S_eDeriv(prob,v,adjoint) def bPrimary(self,prob): diff --git a/simpegEM/Tests/test_FDEM.py b/simpegEM/Tests/test_FDEM.py index f6f71982..fe9e69ed 100644 --- a/simpegEM/Tests/test_FDEM.py +++ b/simpegEM/Tests/test_FDEM.py @@ -5,11 +5,11 @@ import sys from scipy.constants import mu_0 import copy -testDerivs = False -testCrossCheck = True +testDerivs = True +testCrossCheck = False testAdjoint = False testEB = True -testHJ = True +testHJ = False verbose = False @@ -164,54 +164,54 @@ class FDEM_DerivTests(unittest.TestCase): if testEB: 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_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_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_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_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_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_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_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_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_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_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_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_bzi_Bform(self): + # self.assertTrue(derivTest('b', 'bzi')) if testHJ: def test_Jvec_jxr_Jform(self):