From d102ec9c00952f6e49f8b169080626d9cfd7a5fa Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Fri, 6 Nov 2015 10:00:30 -0800 Subject: [PATCH 01/14] Zero, Identity for E formulation --- SimPEG/EM/FDEM/FDEM.py | 64 +++++++++++++++++++----------------- SimPEG/EM/FDEM/FieldsFDEM.py | 15 +++++---- SimPEG/EM/FDEM/SurveyFDEM.py | 5 +-- 3 files changed, 44 insertions(+), 40 deletions(-) diff --git a/SimPEG/EM/FDEM/FDEM.py b/SimPEG/EM/FDEM/FDEM.py index 6d3df44f..260907f5 100644 --- a/SimPEG/EM/FDEM/FDEM.py +++ b/SimPEG/EM/FDEM/FDEM.py @@ -73,26 +73,29 @@ class BaseFDEMProblem(BaseEMProblem): ftype = self._fieldType + 'Solution' u_src = f[src, ftype] dA_dm = self.getADeriv_m(freq, u_src, v) - dRHS_dm = self.getRHSDeriv_m(src, v) - if dRHS_dm is None: - du_dm = dA_duI * ( - dA_dm ) - else: - du_dm = dA_duI * ( - dA_dm + dRHS_dm ) + dRHS_dm = self.getRHSDeriv_m(freq, src, v) + # if dRHS_dm is None: + # du_dm = dA_duI * ( - dA_dm ) + # else: + du_dm = dA_duI * ( - dA_dm + dRHS_dm ) + + for rx in src.rxList: # df_duFun = u.deriv_u(rx.fieldsUsed, m) df_duFun = getattr(f, '_%sDeriv_u'%rx.projField, None) df_du = df_duFun(src, du_dm, adjoint=False) - if df_du is not None: - du_dm = df_du + # if df_du is not None: + + du_dm = df_du df_dmFun = getattr(f, '_%sDeriv_m'%rx.projField, None) df_dm = df_dmFun(src, v, adjoint=False) - if df_dm is not None: - du_dm += df_dm + # if df_dm is not None: + + du_dm += df_dm P = lambda v: rx.projectFieldsDeriv(src, self.mesh, f, v) # wrt u, also have wrt m - Jv[src, rx] = P(du_dm) return Utils.mkvc(Jv) @@ -249,12 +252,11 @@ class ProblemFDEM_e(BaseFDEMProblem): C = self.mesh.edgeCurl MfMui = self.MfMui - # RHS = C.T * (MfMui * S_m) -1j * omega(freq) * Me * S_e RHS = C.T * (MfMui * S_m) -1j * omega(freq) * S_e return RHS - def getRHSDeriv_m(self, src, v, adjoint=False): + def getRHSDeriv_m(self, freq, src, v, adjoint=False): C = self.mesh.edgeCurl MfMui = self.MfMui S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint) @@ -263,25 +265,25 @@ class ProblemFDEM_e(BaseFDEMProblem): dRHS = MfMui * (C * v) S_mDerivv = S_mDeriv(dRHS) S_eDerivv = S_eDeriv(v) - if S_mDerivv is not None and S_eDerivv is not None: - return S_mDerivv - 1j * omega(freq) * S_eDerivv - elif S_mDerivv is not None: - return S_mDerivv - elif S_eDerivv is not None: - return - 1j * omega(freq) * S_eDerivv - else: - return None + # if S_mDerivv is not None and S_eDerivv is not None: + return S_mDerivv - 1j * omega(freq) * S_eDerivv + # elif S_mDerivv is not None: + # return S_mDerivv + # elif S_eDerivv is not None: + # return - 1j * omega(freq) * S_eDerivv + # else: + # return None else: S_mDerivv, S_eDerivv = S_mDeriv(v), S_eDeriv(v) - if S_mDerivv is not None and S_eDerivv is not None: - return C.T * (MfMui * S_mDerivv) -1j * omega(freq) * S_eDerivv - elif S_mDerivv is not None: - return C.T * (MfMui * S_mDerivv) - elif S_eDerivv is not None: - return -1j * omega(freq) * S_eDerivv - else: - return None + # if S_mDerivv is not None and S_eDerivv is not None: + return C.T * (MfMui * S_mDerivv) -1j * omega(freq) * S_eDerivv + # elif S_mDerivv is not None: + # return C.T * (MfMui * S_mDerivv) + # elif S_eDerivv is not None: + # return -1j * omega(freq) * S_eDerivv + # else: + # return None class ProblemFDEM_b(BaseFDEMProblem): @@ -372,7 +374,7 @@ class ProblemFDEM_b(BaseFDEMProblem): return RHS - def getRHSDeriv_m(self, src, v, adjoint=False): + def getRHSDeriv_m(self, freq, src, v, adjoint=False): C = self.mesh.edgeCurl S_m, S_e = src.eval(self) MfMui = self.MfMui @@ -519,7 +521,7 @@ class ProblemFDEM_j(BaseFDEMProblem): return RHS - def getRHSDeriv_m(self, src, v, adjoint=False): + def getRHSDeriv_m(self, freq, src, v, adjoint=False): C = self.mesh.edgeCurl MeMuI = self.MeMuI S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint) @@ -627,7 +629,7 @@ class ProblemFDEM_h(BaseFDEMProblem): return RHS - def getRHSDeriv_m(self, src, v, adjoint=False): + def getRHSDeriv_m(self, freq, src, v, adjoint=False): _, S_e = src.eval(self) C = self.mesh.edgeCurl MfRho = self.MfRho diff --git a/SimPEG/EM/FDEM/FieldsFDEM.py b/SimPEG/EM/FDEM/FieldsFDEM.py index bb786bd1..103ba821 100644 --- a/SimPEG/EM/FDEM/FieldsFDEM.py +++ b/SimPEG/EM/FDEM/FieldsFDEM.py @@ -1,5 +1,6 @@ from SimPEG import Survey, Problem, Utils, np, sp from SimPEG.EM.Utils.EMUtils import omega +from SimPEG.Utils import Zero, Identity class FieldsFDEM(Problem.Fields): @@ -40,11 +41,11 @@ class FieldsFDEM_e(FieldsFDEM): return self._ePrimary(eSolution,srcList) + self._eSecondary(eSolution,srcList) def _eDeriv_u(self, src, v, adjoint = False): - return None + return Identity()*v def _eDeriv_m(self, src, v, adjoint = False): # assuming primary does not depend on the model - return None + return Zero() def _bPrimary(self, eSolution, srcList): bPrimary = np.zeros([self._edgeCurl.shape[0],eSolution.shape[1]],dtype = complex) @@ -73,9 +74,9 @@ class FieldsFDEM_e(FieldsFDEM): def _bSecondaryDeriv_m(self, src, v, adjoint = False): S_mDeriv, _ = src.evalDeriv(self.prob, adjoint) S_mDeriv = S_mDeriv(v) - if S_mDeriv is not None: - return 1./(1j * omega(src.freq)) * S_mDeriv - return None + # 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) @@ -126,11 +127,11 @@ class FieldsFDEM_b(FieldsFDEM): return self._bPrimary(bSolution, srcList) + self._bSecondary(bSolution, srcList) def _bDeriv_u(self, src, v, adjoint=False): - return None + return Identity()*v def _bDeriv_m(self, src, v, adjoint=False): # assuming primary does not depend on the model - return None + return Zero() def _ePrimary(self, bSolution, srcList): ePrimary = np.zeros([self._edgeCurl.shape[1],bSolution.shape[1]],dtype = complex) diff --git a/SimPEG/EM/FDEM/SurveyFDEM.py b/SimPEG/EM/FDEM/SurveyFDEM.py index dbda9f80..5db248cd 100644 --- a/SimPEG/EM/FDEM/SurveyFDEM.py +++ b/SimPEG/EM/FDEM/SurveyFDEM.py @@ -2,6 +2,7 @@ from SimPEG import Survey, Problem, Utils, np, sp from SimPEG.EM.Utils import SrcUtils from SimPEG.EM.Utils.EMUtils import omega, e_from_j, j_from_e, b_from_h, h_from_b from scipy.constants import mu_0 +from SimPEG.Utils import Zero, Identity #################################################### # Receivers @@ -123,10 +124,10 @@ class SrcFDEM(Survey.BaseSrc): return None def S_mDeriv(self, prob, v, adjoint = False): - return None + return Zero() def S_eDeriv(self, prob, v, adjoint = False): - return None + return Zero() class SrcFDEM_RawVec_e(SrcFDEM): From 4553d6db518513eb4cee97ccb6e637e68ea6b65e Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Fri, 6 Nov 2015 10:05:44 -0800 Subject: [PATCH 02/14] b formulation Zero Identity --- SimPEG/EM/FDEM/FDEM.py | 58 ++++++++++++++++++------------------ SimPEG/EM/FDEM/FieldsFDEM.py | 8 ++--- 2 files changed, 33 insertions(+), 33 deletions(-) diff --git a/SimPEG/EM/FDEM/FDEM.py b/SimPEG/EM/FDEM/FDEM.py index 260907f5..a20f1dc7 100644 --- a/SimPEG/EM/FDEM/FDEM.py +++ b/SimPEG/EM/FDEM/FDEM.py @@ -383,40 +383,40 @@ class ProblemFDEM_b(BaseFDEMProblem): if self._makeASymmetric and adjoint: v = self.MfMui * v - if S_e is not None: - MeSigmaIDeriv = self.MeSigmaIDeriv(S_e) - if not adjoint: - RHSderiv = C * (MeSigmaIDeriv * v) - elif adjoint: - RHSderiv = MeSigmaIDeriv.T * (C.T * v) - else: - RHSderiv = None + # if S_e is not None: + MeSigmaIDeriv = self.MeSigmaIDeriv(S_e) + if not adjoint: + RHSderiv = C * (MeSigmaIDeriv * v) + elif adjoint: + RHSderiv = MeSigmaIDeriv.T * (C.T * v) + # else: + # RHSderiv = None S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint) S_mDeriv, S_eDeriv = S_mDeriv(v), S_eDeriv(v) - if S_mDeriv is not None and S_eDeriv is not None: - if not adjoint: - SrcDeriv = S_mDeriv + C * (self.MeSigmaI * S_eDeriv) - elif adjoint: - SrcDeriv = S_mDeriv + Self.MeSigmaI.T * ( C.T * S_eDeriv) - elif S_mDeriv is not None: - SrcDeriv = S_mDeriv - elif S_eDeriv is not None: - if not adjoint: - SrcDeriv = C * (self.MeSigmaI * S_eDeriv) - elif adjoint: - SrcDeriv = self.MeSigmaI.T * ( C.T * S_eDeriv) - else: - SrcDeriv = None + # if S_mDeriv is not None and S_eDeriv is not None: + if not adjoint: + SrcDeriv = S_mDeriv + C * (self.MeSigmaI * S_eDeriv) + elif adjoint: + SrcDeriv = S_mDeriv + Self.MeSigmaI.T * ( C.T * S_eDeriv) + # elif S_mDeriv is not None: + # SrcDeriv = S_mDeriv + # elif S_eDeriv is not None: + # if not adjoint: + # SrcDeriv = C * (self.MeSigmaI * S_eDeriv) + # elif adjoint: + # SrcDeriv = self.MeSigmaI.T * ( C.T * S_eDeriv) + # else: + # SrcDeriv = None - if RHSderiv is not None and SrcDeriv is not None: - RHSderiv += SrcDeriv - elif SrcDeriv is not None: - RHSderiv = SrcDeriv + # if RHSderiv is not None and SrcDeriv is not None: + RHSderiv += SrcDeriv + # elif SrcDeriv is not None: + # RHSderiv = SrcDeriv - if RHSderiv is not None: - if self._makeASymmetric is True and not adjoint: - return MfMui.T * RHSderiv + # if RHSderiv is not None: + if self._makeASymmetric is True and not adjoint: + return MfMui.T * RHSderiv return RHSderiv diff --git a/SimPEG/EM/FDEM/FieldsFDEM.py b/SimPEG/EM/FDEM/FieldsFDEM.py index 103ba821..00b0d155 100644 --- a/SimPEG/EM/FDEM/FieldsFDEM.py +++ b/SimPEG/EM/FDEM/FieldsFDEM.py @@ -145,8 +145,8 @@ class FieldsFDEM_b(FieldsFDEM): e = self._MeSigmaI * ( self._edgeCurl.T * ( self._MfMui * bSolution)) for i,src in enumerate(srcList): _,S_e = src.eval(self.prob) - if S_e is not None: - e[:,i] += -self._MeSigmaI * S_e + # if S_e is not None: + e[:,i] += -self._MeSigmaI * S_e return e def _eSecondaryDeriv_u(self, src, v, adjoint=False): @@ -164,8 +164,8 @@ class FieldsFDEM_b(FieldsFDEM): Me = Me.T w = self._edgeCurl.T * (self._MfMui * bSolution) - if S_e is not None: - w += -Utils.mkvc(Me * S_e,2) + # if S_e is not None: + w += -Utils.mkvc(Me * S_e,2) if not adjoint: de_dm = self._MeSigmaIDeriv(w) * v From 0af1f598916754a586567e0a1e9cf57f616999f1 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Fri, 6 Nov 2015 10:58:23 -0800 Subject: [PATCH 03/14] h,j jvec now with Zero and identity classes --- SimPEG/EM/FDEM/FDEM.py | 75 +++++++++++++++++++----------------- SimPEG/EM/FDEM/FieldsFDEM.py | 22 +++++------ SimPEG/EM/FDEM/SurveyFDEM.py | 4 +- 3 files changed, 53 insertions(+), 48 deletions(-) diff --git a/SimPEG/EM/FDEM/FDEM.py b/SimPEG/EM/FDEM/FDEM.py index a20f1dc7..f0fd080a 100644 --- a/SimPEG/EM/FDEM/FDEM.py +++ b/SimPEG/EM/FDEM/FDEM.py @@ -4,6 +4,7 @@ from SurveyFDEM import SurveyFDEM from FieldsFDEM import FieldsFDEM, FieldsFDEM_e, FieldsFDEM_b, FieldsFDEM_h, FieldsFDEM_j from SimPEG.EM.Base import BaseEMProblem from SimPEG.EM.Utils.EMUtils import omega +from SimPEG.Utils import Zero, Identity class BaseFDEMProblem(BaseEMProblem): @@ -73,7 +74,7 @@ class BaseFDEMProblem(BaseEMProblem): ftype = self._fieldType + 'Solution' u_src = f[src, ftype] dA_dm = self.getADeriv_m(freq, u_src, v) - dRHS_dm = self.getRHSDeriv_m(freq, src, v) + dRHS_dm = self.getRHSDeriv_m(freq, src, v) # if dRHS_dm is None: # du_dm = dA_duI * ( - dA_dm ) # else: @@ -532,25 +533,25 @@ class ProblemFDEM_j(BaseFDEMProblem): v = MfRho*v S_mDerivv = S_mDeriv(MeMuI.T * (C.T * v)) S_eDerivv = S_eDeriv(v) - if S_mDerivv is not None and S_eDerivv is not None: - return S_mDerivv - 1j * omega(freq) * S_eDerivv - elif S_mDerivv is not None: - return S_mDerivv - elif S_eDerivv is not None: - return - 1j * omega(freq) * S_eDerivv - else: - return None + # if S_mDerivv is not None and S_eDerivv is not None: + return S_mDerivv - 1j * omega(freq) * S_eDerivv + # elif S_mDerivv is not None: + # return S_mDerivv + # elif S_eDerivv is not None: + # return - 1j * omega(freq) * S_eDerivv + # else: + # return None else: S_mDerivv, S_eDerivv = S_mDeriv(v), S_eDeriv(v) - if S_mDerivv is not None and S_eDerivv is not None: - RHSDeriv = C * (MeMuI * S_mDerivv) - 1j * omega(freq) * S_eDerivv - elif S_mDerivv is not None: - RHSDeriv = C * (MeMuI * S_mDerivv) - elif S_eDerivv is not None: - RHSDeriv = - 1j * omega(freq) * S_eDerivv - else: - return None + # if S_mDerivv is not None and S_eDerivv is not None: + RHSDeriv = C * (MeMuI * S_mDerivv) - 1j * omega(freq) * S_eDerivv + # elif S_mDerivv is not None: + # RHSDeriv = C * (MeMuI * S_mDerivv) + # elif S_eDerivv is not None: + # RHSDeriv = - 1j * omega(freq) * S_eDerivv + # else: + # return None if self._makeASymmetric: MfRho = self.MfRho @@ -634,30 +635,34 @@ class ProblemFDEM_h(BaseFDEMProblem): C = self.mesh.edgeCurl MfRho = self.MfRho - RHSDeriv = None + # RHSDeriv = Zero() - if S_e is not None: - MfRhoDeriv = self.MfRhoDeriv(S_e) - if not adjoint: - RHSDeriv = C.T * (MfRhoDeriv * v) - elif adjoint: - RHSDeriv = MfRhoDeriv.T * (C * v) + # if S_e is not None: + MfRhoDeriv = self.MfRhoDeriv(S_e) + if not adjoint: + RHSDeriv = C.T * (MfRhoDeriv * v) + elif adjoint: + RHSDeriv = MfRhoDeriv.T * (C * v) + + # print S_e, RHSDeriv S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint) S_mDeriv = S_mDeriv(v) S_eDeriv = S_eDeriv(v) - if S_mDeriv is not None: - if RHSDeriv is not None: - RHSDeriv += S_mDeriv(v) - else: - RHSDeriv = S_mDeriv(v) - if S_eDeriv is not None: - if RHSDeriv is not None: - RHSDeriv += C.T * (MfRho * S_e) - else: - RHSDeriv = C.T * (MfRho * S_e) + # if S_mDeriv is not None: + # if RHSDeriv is not None: + RHSDeriv += S_mDeriv - return RHSDeriv + # else: + # RHSDeriv = S_mDeriv(v) + # if S_eDeriv is not None: + # if RHSDeriv is not None: + return RHSDeriv + C.T * (MfRho * S_eDeriv) + # print RHSDeriv + # # else: + # # RHSDeriv = C.T * (MfRho * S_e) + # print RHSDeriv + # return RHSDeriv diff --git a/SimPEG/EM/FDEM/FieldsFDEM.py b/SimPEG/EM/FDEM/FieldsFDEM.py index 00b0d155..32f04d58 100644 --- a/SimPEG/EM/FDEM/FieldsFDEM.py +++ b/SimPEG/EM/FDEM/FieldsFDEM.py @@ -228,11 +228,11 @@ class FieldsFDEM_j(FieldsFDEM): return self._jPrimary(jSolution, srcList) + self._jSecondary(jSolution, srcList) def _jDeriv_u(self, src, v, adjoint=False): - return None + return Identity()*v def _jDeriv_m(self, src, v, adjoint=False): # assuming primary does not depend on the model - return None + return Zero() def _hPrimary(self, jSolution, srcList): hPrimary = np.zeros([self._edgeCurl.shape[1],jSolution.shape[1]],dtype = complex) @@ -274,12 +274,12 @@ class FieldsFDEM_j(FieldsFDEM): if not adjoint: S_mDeriv = S_mDeriv(v) - if S_mDeriv is not None: - hDeriv_m += 1./(1j*omega(src.freq)) * MeMuI * (Me * S_mDeriv) + # if S_mDeriv is not None: + hDeriv_m += 1./(1j*omega(src.freq)) * MeMuI * (Me * S_mDeriv) elif adjoint: S_mDeriv = S_mDeriv(Me.T * (MeMuI.T * v)) - if S_mDeriv is not None: - hDeriv_m += 1./(1j*omega(src.freq)) * S_mDeriv + # if S_mDeriv is not None: + hDeriv_m += 1./(1j*omega(src.freq)) * S_mDeriv return hDeriv_m @@ -329,11 +329,11 @@ class FieldsFDEM_h(FieldsFDEM): return self._hPrimary(hSolution, srcList) + self._hSecondary(hSolution, srcList) def _hDeriv_u(self, src, v, adjoint=False): - return None + return Identity()*v def _hDeriv_m(self, src, v, adjoint=False): # assuming primary does not depend on the model - return None + return Zero() def _jPrimary(self, hSolution, srcList): jPrimary = np.zeros([self._edgeCurl.shape[0], hSolution.shape[1]], dtype = complex) @@ -360,9 +360,9 @@ class FieldsFDEM_h(FieldsFDEM): def _jSecondaryDeriv_m(self, src, v, adjoint=False): _,S_eDeriv = src.evalDeriv(self.prob, adjoint) S_eDeriv = S_eDeriv(v) - if S_eDeriv is not None: - return -S_eDeriv - return None + # if S_eDeriv is not None: + return -S_eDeriv + # return None def _j(self, hSolution, srcList): return self._jPrimary(hSolution, srcList) + self._jSecondary(hSolution, srcList) diff --git a/SimPEG/EM/FDEM/SurveyFDEM.py b/SimPEG/EM/FDEM/SurveyFDEM.py index 5db248cd..a962fe7e 100644 --- a/SimPEG/EM/FDEM/SurveyFDEM.py +++ b/SimPEG/EM/FDEM/SurveyFDEM.py @@ -118,10 +118,10 @@ class SrcFDEM(Survey.BaseSrc): return None def S_m(self, prob): - return None + return Zero() def S_e(self, prob): - return None + return Zero() def S_mDeriv(self, prob, v, adjoint = False): return Zero() From 639a6e7a06f9346c6662fd83e0707d3a12d4cd7f Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sun, 8 Nov 2015 14:09:18 -0800 Subject: [PATCH 04/14] Adjoint for all 4 formulations with Zero Identity --- SimPEG/EM/FDEM/FDEM.py | 32 ++++++++++++++++---------------- SimPEG/EM/FDEM/SurveyFDEM.py | 4 ++-- 2 files changed, 18 insertions(+), 18 deletions(-) diff --git a/SimPEG/EM/FDEM/FDEM.py b/SimPEG/EM/FDEM/FDEM.py index f0fd080a..37c79236 100644 --- a/SimPEG/EM/FDEM/FDEM.py +++ b/SimPEG/EM/FDEM/FDEM.py @@ -130,30 +130,30 @@ class BaseFDEMProblem(BaseEMProblem): df_duTFun = getattr(f, '_%sDeriv_u'%rx.projField, None) df_duT = df_duTFun(src, PTv, adjoint=True) - if df_duT is not None: - dA_duIT = ATinv * df_duT - else: - dA_duIT = ATinv * PTv + # if df_duT is not None: + dA_duIT = ATinv * df_duT + # else: + # dA_duIT = ATinv * PTv dA_dmT = self.getADeriv_m(freq, u_src, dA_duIT, adjoint=True) - dRHS_dmT = self.getRHSDeriv_m(src, dA_duIT, adjoint=True) + dRHS_dmT = self.getRHSDeriv_m(freq,src, dA_duIT, adjoint=True) - if dRHS_dmT is None: - du_dmT = - dA_dmT - else: - du_dmT = -dA_dmT + dRHS_dmT + # if dRHS_dmT is None: + # du_dmT = - dA_dmT + # else: + du_dmT = -dA_dmT + dRHS_dmT df_dmFun = getattr(f, '_%sDeriv_m'%rx.projField, None) dfT_dm = df_dmFun(src, PTv, adjoint=True) - if dfT_dm is not None: - du_dmT += dfT_dm + # if dfT_dm is not None: + du_dmT += dfT_dm real_or_imag = rx.projComp - if real_or_imag == 'real': - Jtv += du_dmT.real - elif real_or_imag == 'imag': - Jtv += - du_dmT.real + if real_or_imag is 'real': + Jtv += np.array(du_dmT,dtype=complex).real + elif real_or_imag is 'imag': + Jtv += - np.array(du_dmT,dtype=complex).real else: raise Exception('Must be real or imag') @@ -399,7 +399,7 @@ class ProblemFDEM_b(BaseFDEMProblem): if not adjoint: SrcDeriv = S_mDeriv + C * (self.MeSigmaI * S_eDeriv) elif adjoint: - SrcDeriv = S_mDeriv + Self.MeSigmaI.T * ( C.T * S_eDeriv) + SrcDeriv = S_mDeriv + self.MeSigmaI.T * ( C.T * S_eDeriv) # elif S_mDeriv is not None: # SrcDeriv = S_mDeriv # elif S_eDeriv is not None: diff --git a/SimPEG/EM/FDEM/SurveyFDEM.py b/SimPEG/EM/FDEM/SurveyFDEM.py index a962fe7e..5db248cd 100644 --- a/SimPEG/EM/FDEM/SurveyFDEM.py +++ b/SimPEG/EM/FDEM/SurveyFDEM.py @@ -118,10 +118,10 @@ class SrcFDEM(Survey.BaseSrc): return None def S_m(self, prob): - return Zero() + return None def S_e(self, prob): - return Zero() + return None def S_mDeriv(self, prob, v, adjoint = False): return Zero() From 297ea916b49d70317bfaf3de37ad96556aa43819 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sun, 8 Nov 2015 14:58:28 -0800 Subject: [PATCH 05/14] light notation clean up --- SimPEG/EM/FDEM/FDEM.py | 63 ++++++++++-------------------------------- 1 file changed, 14 insertions(+), 49 deletions(-) diff --git a/SimPEG/EM/FDEM/FDEM.py b/SimPEG/EM/FDEM/FDEM.py index 37c79236..349b3bdb 100644 --- a/SimPEG/EM/FDEM/FDEM.py +++ b/SimPEG/EM/FDEM/FDEM.py @@ -67,37 +67,28 @@ class BaseFDEMProblem(BaseEMProblem): Jv = self.dataPair(self.survey) for freq in self.survey.freqs: - dA_du = self.getA(freq) # - dA_duI = self.Solver(dA_du, **self.solverOpts) + A = self.getA(freq) # + Ainv = self.Solver(A, **self.solverOpts) for src in self.survey.getSrcByFreq(freq): ftype = self._fieldType + 'Solution' u_src = f[src, ftype] dA_dm = self.getADeriv_m(freq, u_src, v) dRHS_dm = self.getRHSDeriv_m(freq, src, v) - # if dRHS_dm is None: - # du_dm = dA_duI * ( - dA_dm ) - # else: - du_dm = dA_duI * ( - dA_dm + dRHS_dm ) + du_dm = Ainv * ( - dA_dm + dRHS_dm ) - for rx in src.rxList: - # df_duFun = u.deriv_u(rx.fieldsUsed, m) df_duFun = getattr(f, '_%sDeriv_u'%rx.projField, None) - df_du = df_duFun(src, du_dm, adjoint=False) - # if df_du is not None: - - du_dm = df_du + df_dudu_dm = df_duFun(src, du_dm, adjoint=False) df_dmFun = getattr(f, '_%sDeriv_m'%rx.projField, None) df_dm = df_dmFun(src, v, adjoint=False) - # if df_dm is not None: - - du_dm += df_dm + + Df_Dm = np.array(df_dudu_dm + df_dm,dtype=complex) P = lambda v: rx.projectFieldsDeriv(src, self.mesh, f, v) # wrt u, also have wrt m - Jv[src, rx] = P(du_dm) + Jv[src, rx] = P(Df_Dm) return Utils.mkvc(Jv) @@ -130,23 +121,16 @@ class BaseFDEMProblem(BaseEMProblem): df_duTFun = getattr(f, '_%sDeriv_u'%rx.projField, None) df_duT = df_duTFun(src, PTv, adjoint=True) - # if df_duT is not None: - dA_duIT = ATinv * df_duT - # else: - # dA_duIT = ATinv * PTv + + ATinvdf_duT = ATinv * df_duT - dA_dmT = self.getADeriv_m(freq, u_src, dA_duIT, adjoint=True) - - dRHS_dmT = self.getRHSDeriv_m(freq,src, dA_duIT, adjoint=True) - - # if dRHS_dmT is None: - # du_dmT = - dA_dmT - # else: + dA_dmT = self.getADeriv_m(freq, u_src, ATinvdf_duT, adjoint=True) + dRHS_dmT = self.getRHSDeriv_m(freq,src, ATinvdf_duT, adjoint=True) du_dmT = -dA_dmT + dRHS_dmT df_dmFun = getattr(f, '_%sDeriv_m'%rx.projField, None) dfT_dm = df_dmFun(src, PTv, adjoint=True) - # if dfT_dm is not None: + du_dmT += dfT_dm real_or_imag = rx.projComp @@ -379,43 +363,27 @@ class ProblemFDEM_b(BaseFDEMProblem): C = self.mesh.edgeCurl S_m, S_e = src.eval(self) MfMui = self.MfMui - # Me = self.Me if self._makeASymmetric and adjoint: v = self.MfMui * v - # if S_e is not None: MeSigmaIDeriv = self.MeSigmaIDeriv(S_e) if not adjoint: RHSderiv = C * (MeSigmaIDeriv * v) elif adjoint: RHSderiv = MeSigmaIDeriv.T * (C.T * v) - # else: - # RHSderiv = None S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint) S_mDeriv, S_eDeriv = S_mDeriv(v), S_eDeriv(v) - # if S_mDeriv is not None and S_eDeriv is not None: + if not adjoint: SrcDeriv = S_mDeriv + C * (self.MeSigmaI * S_eDeriv) elif adjoint: SrcDeriv = S_mDeriv + self.MeSigmaI.T * ( C.T * S_eDeriv) - # elif S_mDeriv is not None: - # SrcDeriv = S_mDeriv - # elif S_eDeriv is not None: - # if not adjoint: - # SrcDeriv = C * (self.MeSigmaI * S_eDeriv) - # elif adjoint: - # SrcDeriv = self.MeSigmaI.T * ( C.T * S_eDeriv) - # else: - # SrcDeriv = None - # if RHSderiv is not None and SrcDeriv is not None: + print RHSderiv, SrcDeriv RHSderiv += SrcDeriv - # elif SrcDeriv is not None: - # RHSderiv = SrcDeriv - # if RHSderiv is not None: if self._makeASymmetric is True and not adjoint: return MfMui.T * RHSderiv @@ -635,9 +603,6 @@ class ProblemFDEM_h(BaseFDEMProblem): C = self.mesh.edgeCurl MfRho = self.MfRho - # RHSDeriv = Zero() - - # if S_e is not None: MfRhoDeriv = self.MfRhoDeriv(S_e) if not adjoint: RHSDeriv = C.T * (MfRhoDeriv * v) From 915a2f84464cfd7ee33e0c1e2352ff5aecfd5345 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sun, 8 Nov 2015 15:18:19 -0800 Subject: [PATCH 06/14] cleaned out commented-out code --- SimPEG/EM/FDEM/FDEM.py | 86 +++++++----------------------------------- 1 file changed, 13 insertions(+), 73 deletions(-) diff --git a/SimPEG/EM/FDEM/FDEM.py b/SimPEG/EM/FDEM/FDEM.py index 349b3bdb..559f4f7c 100644 --- a/SimPEG/EM/FDEM/FDEM.py +++ b/SimPEG/EM/FDEM/FDEM.py @@ -248,27 +248,10 @@ class ProblemFDEM_e(BaseFDEMProblem): if adjoint: dRHS = MfMui * (C * v) - S_mDerivv = S_mDeriv(dRHS) - S_eDerivv = S_eDeriv(v) - # if S_mDerivv is not None and S_eDerivv is not None: - return S_mDerivv - 1j * omega(freq) * S_eDerivv - # elif S_mDerivv is not None: - # return S_mDerivv - # elif S_eDerivv is not None: - # return - 1j * omega(freq) * S_eDerivv - # else: - # return None - else: - S_mDerivv, S_eDerivv = S_mDeriv(v), S_eDeriv(v) + return S_mDeriv(dRHS) - 1j * omega(freq) * S_eDeriv(v) - # if S_mDerivv is not None and S_eDerivv is not None: - return C.T * (MfMui * S_mDerivv) -1j * omega(freq) * S_eDerivv - # elif S_mDerivv is not None: - # return C.T * (MfMui * S_mDerivv) - # elif S_eDerivv is not None: - # return -1j * omega(freq) * S_eDerivv - # else: - # return None + else: + return C.T * (MfMui * S_mDeriv(v)) -1j * omega(freq) * S_eDeriv(v) class ProblemFDEM_b(BaseFDEMProblem): @@ -349,7 +332,6 @@ class ProblemFDEM_b(BaseFDEMProblem): S_m, S_e = self.getSourceTerm(freq) C = self.mesh.edgeCurl MeSigmaI = self.MeSigmaI - # Me = self.Me RHS = S_m + C * ( MeSigmaI * S_e ) @@ -368,26 +350,19 @@ class ProblemFDEM_b(BaseFDEMProblem): v = self.MfMui * v MeSigmaIDeriv = self.MeSigmaIDeriv(S_e) + S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint) + if not adjoint: RHSderiv = C * (MeSigmaIDeriv * v) + SrcDeriv = S_mDeriv(v) + C * (self.MeSigmaI * S_eDeriv(v)) elif adjoint: RHSderiv = MeSigmaIDeriv.T * (C.T * v) - - S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint) - S_mDeriv, S_eDeriv = S_mDeriv(v), S_eDeriv(v) - - if not adjoint: - SrcDeriv = S_mDeriv + C * (self.MeSigmaI * S_eDeriv) - elif adjoint: - SrcDeriv = S_mDeriv + self.MeSigmaI.T * ( C.T * S_eDeriv) - - print RHSderiv, SrcDeriv - RHSderiv += SrcDeriv + SrcDeriv = S_mDeriv(v) + self.MeSigmaI.T * (C.T * S_eDeriv(v)) if self._makeASymmetric is True and not adjoint: - return MfMui.T * RHSderiv + return MfMui.T * (RHSderiv + SrcDeriv) - return RHSderiv + return RHSderiv + SrcDeriv @@ -499,27 +474,10 @@ class ProblemFDEM_j(BaseFDEMProblem): if self._makeASymmetric: MfRho = self.MfRho v = MfRho*v - S_mDerivv = S_mDeriv(MeMuI.T * (C.T * v)) - S_eDerivv = S_eDeriv(v) - # if S_mDerivv is not None and S_eDerivv is not None: - return S_mDerivv - 1j * omega(freq) * S_eDerivv - # elif S_mDerivv is not None: - # return S_mDerivv - # elif S_eDerivv is not None: - # return - 1j * omega(freq) * S_eDerivv - # else: - # return None - else: - S_mDerivv, S_eDerivv = S_mDeriv(v), S_eDeriv(v) + return S_mDeriv(MeMuI.T * (C.T * v)) - 1j * omega(freq) * S_eDeriv(v) - # if S_mDerivv is not None and S_eDerivv is not None: - RHSDeriv = C * (MeMuI * S_mDerivv) - 1j * omega(freq) * S_eDerivv - # elif S_mDerivv is not None: - # RHSDeriv = C * (MeMuI * S_mDerivv) - # elif S_eDerivv is not None: - # RHSDeriv = - 1j * omega(freq) * S_eDerivv - # else: - # return None + else: + RHSDeriv = C * (MeMuI * S_mDeriv(v)) - 1j * omega(freq) * S_eDeriv(v) if self._makeASymmetric: MfRho = self.MfRho @@ -609,25 +567,7 @@ class ProblemFDEM_h(BaseFDEMProblem): elif adjoint: RHSDeriv = MfRhoDeriv.T * (C * v) - # print S_e, RHSDeriv - S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint) - S_mDeriv = S_mDeriv(v) - S_eDeriv = S_eDeriv(v) - - # if S_mDeriv is not None: - # if RHSDeriv is not None: - RHSDeriv += S_mDeriv - - # else: - # RHSDeriv = S_mDeriv(v) - # if S_eDeriv is not None: - # if RHSDeriv is not None: - return RHSDeriv + C.T * (MfRho * S_eDeriv) - # print RHSDeriv - # # else: - # # RHSDeriv = C.T * (MfRho * S_e) - # print RHSDeriv - # return RHSDeriv + return RHSDeriv + S_mDeriv(v) + C.T * (MfRho * S_eDeriv(v)) From 554f507ba754144d2d12f9fa0b1703304edf5004 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Mon, 9 Nov 2015 09:15:02 -0800 Subject: [PATCH 07/14] switch order of addition to satisfy Zero class --- SimPEG/EM/FDEM/FDEM.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SimPEG/EM/FDEM/FDEM.py b/SimPEG/EM/FDEM/FDEM.py index 559f4f7c..0f3d5a09 100644 --- a/SimPEG/EM/FDEM/FDEM.py +++ b/SimPEG/EM/FDEM/FDEM.py @@ -360,7 +360,7 @@ class ProblemFDEM_b(BaseFDEMProblem): SrcDeriv = S_mDeriv(v) + self.MeSigmaI.T * (C.T * S_eDeriv(v)) if self._makeASymmetric is True and not adjoint: - return MfMui.T * (RHSderiv + SrcDeriv) + return MfMui.T * (SrcDeriv + RHSderiv) return RHSderiv + SrcDeriv From ff68400320cf1660eb9c94b3e9c2077390a1aef1 Mon Sep 17 00:00:00 2001 From: Brendan Smithyman Date: Tue, 10 Nov 2015 12:26:55 -0500 Subject: [PATCH 08/14] Be compatible with outside metaclasses, for now Making this change would allow use of metaclasses in a more straightforward way for code that inherits from both SimPEG classes and external classes. In the future, this may need to be addressed in more detail. --- SimPEG/Utils/codeutils.py | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/SimPEG/Utils/codeutils.py b/SimPEG/Utils/codeutils.py index 0ba57b2e..4a9a28a7 100644 --- a/SimPEG/Utils/codeutils.py +++ b/SimPEG/Utils/codeutils.py @@ -3,10 +3,7 @@ import time import numpy as np from functools import wraps - -class SimPEGMetaClass(type): - def __new__(cls, name, bases, attrs): - return super(SimPEGMetaClass, cls).__new__(cls, name, bases, attrs) +SimPEGMetaClass = type def memProfileWrapper(towrap, *funNames): """ From 815311bfecf92d65dfc7d5cbd7d96d821843675f Mon Sep 17 00:00:00 2001 From: Lindsey Date: Mon, 23 Nov 2015 18:03:49 -0800 Subject: [PATCH 09/14] cleaned up sources --- SimPEG/EM/FDEM/FieldsFDEM.py | 6 ++-- SimPEG/EM/FDEM/SrcFDEM.py | 36 ++-------------------- tests/em/fdem/forward/test_FDEM_forward.py | 2 +- 3 files changed, 5 insertions(+), 39 deletions(-) diff --git a/SimPEG/EM/FDEM/FieldsFDEM.py b/SimPEG/EM/FDEM/FieldsFDEM.py index be766bcc..fa6a6485 100644 --- a/SimPEG/EM/FDEM/FieldsFDEM.py +++ b/SimPEG/EM/FDEM/FieldsFDEM.py @@ -53,8 +53,7 @@ class Fields_e(Fields): bPrimary = np.zeros([self._edgeCurl.shape[0],eSolution.shape[1]],dtype = complex) for i, src in enumerate(srcList): bp = src.bPrimary(self.prob) - if bp is not None: - bPrimary[:,i] += bp + bPrimary[:,i] += bp return bPrimary def _bSecondary(self, eSolution, srcList): @@ -147,7 +146,6 @@ class Fields_b(Fields): e = self._MeSigmaI * ( self._edgeCurl.T * ( self._MfMui * bSolution)) for i,src in enumerate(srcList): _,S_e = src.eval(self.prob) - # if S_e is not None: e[:,i] += -self._MeSigmaI * S_e return e @@ -316,7 +314,7 @@ class Fields_h(Fields): for i, src in enumerate(srcList): hp = src.hPrimary(self.prob) hPrimary[:,i] += hp - return hPrimary + return hPrimary def _hSecondary(self, hSolution, srcList): return hSolution diff --git a/SimPEG/EM/FDEM/SrcFDEM.py b/SimPEG/EM/FDEM/SrcFDEM.py index c4303731..b29768ac 100644 --- a/SimPEG/EM/FDEM/SrcFDEM.py +++ b/SimPEG/EM/FDEM/SrcFDEM.py @@ -52,30 +52,14 @@ class RawVec_e(BaseSrc): :param rxList: receiver list """ - def __init__(self, rxList, freq, S_e, ePrimary=None, bPrimary=None, hPrimary=None, jPrimary=None): + def __init__(self, rxList, freq, S_e): #, ePrimary=None, bPrimary=None, hPrimary=None, jPrimary=None): self._S_e = np.array(S_e,dtype=complex) - self._ePrimary = ePrimary - self._bPrimary = bPrimary - self._hPrimary = hPrimary - self._jPrimary = jPrimary self.freq = float(freq) BaseSrc.__init__(self, rxList) def S_e(self, prob): return self._S_e - def ePrimary(self, prob): - return self._ePrimary - - def bPrimary(self, prob): - return self._bPrimary - - def hPrimary(self, prob): - return self._hPrimary - - def jPrimary(self, prob): - return self._jPrimary - class RawVec_m(BaseSrc): """ @@ -86,32 +70,16 @@ class RawVec_m(BaseSrc): :param rxList: receiver list """ - def __init__(self, rxList, freq, S_m, integrate = True, ePrimary=None, bPrimary=None, hPrimary=None, jPrimary=None): + def __init__(self, rxList, freq, S_m, integrate = True): #ePrimary=Zero(), bPrimary=Zero(), hPrimary=Zero(), jPrimary=Zero()): self._S_m = np.array(S_m,dtype=complex) self.freq = float(freq) self.integrate = integrate - self._ePrimary = np.array(ePrimary,dtype=complex) - self._bPrimary = np.array(bPrimary,dtype=complex) - self._hPrimary = np.array(hPrimary,dtype=complex) - self._jPrimary = np.array(jPrimary,dtype=complex) BaseSrc.__init__(self, rxList) def S_m(self, prob): return self._S_m - def ePrimary(self, prob): - return self._ePrimary - - def bPrimary(self, prob): - return self._bPrimary - - def hPrimary(self, prob): - return self._hPrimary - - def jPrimary(self, prob): - return self._jPrimary - class RawVec(BaseSrc): """ diff --git a/tests/em/fdem/forward/test_FDEM_forward.py b/tests/em/fdem/forward/test_FDEM_forward.py index ce6aed4e..437f3708 100644 --- a/tests/em/fdem/forward/test_FDEM_forward.py +++ b/tests/em/fdem/forward/test_FDEM_forward.py @@ -6,7 +6,7 @@ from scipy.constants import mu_0 from SimPEG.EM.Utils.testingUtils import getFDEMProblem testEB = True -testHJ = False +testHJ = True verbose = False From f3fb1e64812c38bf79eba65dd654ba1ed84efae3 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 24 Nov 2015 10:05:28 -0800 Subject: [PATCH 10/14] added an __iadd__, __isub__ to zero class --- SimPEG/Utils/matutils.py | 2 ++ tests/utils/test_Zero.py | 7 +++++++ 2 files changed, 9 insertions(+) diff --git a/SimPEG/Utils/matutils.py b/SimPEG/Utils/matutils.py index ee58bf86..b38bb4a1 100644 --- a/SimPEG/Utils/matutils.py +++ b/SimPEG/Utils/matutils.py @@ -399,8 +399,10 @@ def diagEst(matFun, n, k=None, approach='Probing'): class Zero(object): def __add__(self, v):return v def __radd__(self, v):return v + def __iadd__(self, v):return v def __sub__(self, v):return -v def __rsub__(self, v):return v + def __isub__(self, v):return v def __mul__(self, v):return self def __rmul__(self, v):return self def __div__(self, v): return self diff --git a/tests/utils/test_Zero.py b/tests/utils/test_Zero.py index be7c9bbe..594de6a6 100644 --- a/tests/utils/test_Zero.py +++ b/tests/utils/test_Zero.py @@ -20,6 +20,13 @@ class Tests(unittest.TestCase): assert 3*z == 0 assert z*3 == 0 assert z/3 == 0 + + a = 1 + a += z + assert a == 1 + a = 1 + a += z + assert a == 1 self.assertRaises(ZeroDivisionError, lambda:3/z) def test_mat_zero(self): From 9b847c0e32c099d87733eb54c6a04ae9f8591bec Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 24 Nov 2015 10:06:15 -0800 Subject: [PATCH 11/14] replaced all numpy (a += b) with (a = a + b) because it plays more nicely with the Zero class --- SimPEG/EM/FDEM/FDEM.py | 4 ++-- SimPEG/EM/FDEM/FieldsFDEM.py | 36 +++++++++++++++--------------------- 2 files changed, 17 insertions(+), 23 deletions(-) diff --git a/SimPEG/EM/FDEM/FDEM.py b/SimPEG/EM/FDEM/FDEM.py index c91d6c59..eec4bab9 100644 --- a/SimPEG/EM/FDEM/FDEM.py +++ b/SimPEG/EM/FDEM/FDEM.py @@ -161,8 +161,8 @@ class BaseFDEMProblem(BaseEMProblem): for i, src in enumerate(Srcs): smi, sei = src.eval(self) - S_m[:,i] += smi - S_e[:,i] += sei + S_m[:,i] = S_m[:,i] + smi + S_e[:,i] = S_e[:,i] + sei return S_m, S_e diff --git a/SimPEG/EM/FDEM/FieldsFDEM.py b/SimPEG/EM/FDEM/FieldsFDEM.py index fa6a6485..f90dbfef 100644 --- a/SimPEG/EM/FDEM/FieldsFDEM.py +++ b/SimPEG/EM/FDEM/FieldsFDEM.py @@ -33,7 +33,7 @@ class Fields_e(Fields): ePrimary = np.zeros_like(eSolution) for i, src in enumerate(srcList): ep = src.ePrimary(self.prob) - ePrimary[:,i] += ep + ePrimary[:,i] = ePrimary[:,i] + ep return ePrimary def _eSecondary(self, eSolution, srcList): @@ -53,7 +53,7 @@ class Fields_e(Fields): bPrimary = np.zeros([self._edgeCurl.shape[0],eSolution.shape[1]],dtype = complex) for i, src in enumerate(srcList): bp = src.bPrimary(self.prob) - bPrimary[:,i] += bp + bPrimary[:,i] = bPrimary[:,i] + bp return bPrimary def _bSecondary(self, eSolution, srcList): @@ -62,8 +62,7 @@ class Fields_e(Fields): for i, src in enumerate(srcList): b[:,i] *= - 1./(1j*omega(src.freq)) S_m, _ = src.eval(self.prob) - if S_m is not None: - b[:,i] += 1./(1j*omega(src.freq)) * S_m + b[:,i] = b[:,i]+ 1./(1j*omega(src.freq)) * S_m return b def _bSecondaryDeriv_u(self, src, v, adjoint = False): @@ -75,9 +74,7 @@ class Fields_e(Fields): def _bSecondaryDeriv_m(self, src, v, adjoint = False): S_mDeriv, _ = src.evalDeriv(self.prob, adjoint) S_mDeriv = S_mDeriv(v) - # 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) @@ -117,8 +114,7 @@ class Fields_b(Fields): bPrimary = np.zeros_like(bSolution) for i, src in enumerate(srcList): bp = src.bPrimary(self.prob) - # if bp is not None: - bPrimary[:,i] += bp + bPrimary[:,i] = bPrimary[:,i] + bp return bPrimary def _bSecondary(self, bSolution, srcList): @@ -138,15 +134,14 @@ class Fields_b(Fields): ePrimary = np.zeros([self._edgeCurl.shape[1],bSolution.shape[1]],dtype = complex) for i,src in enumerate(srcList): ep = src.ePrimary(self.prob) - if ep is not None: - ePrimary[:,i] += ep + ePrimary[:,i] = ePrimary[:,i] + ep return ePrimary def _eSecondary(self, bSolution, srcList): e = self._MeSigmaI * ( self._edgeCurl.T * ( self._MfMui * bSolution)) for i,src in enumerate(srcList): _,S_e = src.eval(self.prob) - e[:,i] += -self._MeSigmaI * S_e + e[:,i] = e[:,i]+ -self._MeSigmaI * S_e return e def _eSecondaryDeriv_u(self, src, v, adjoint=False): @@ -164,7 +159,6 @@ class Fields_b(Fields): Me = Me.T w = self._edgeCurl.T * (self._MfMui * bSolution) - # if S_e is not None: w += -Utils.mkvc(Me * S_e,2) if not adjoint: @@ -176,7 +170,7 @@ class Fields_b(Fields): Se_Deriv = S_eDeriv(v) if Se_Deriv is not None: - de_dm += -self._MeSigmaI * Se_Deriv + de_dm = de_dm-self._MeSigmaI * Se_Deriv return de_dm @@ -217,7 +211,7 @@ class Fields_j(Fields): jPrimary = np.zeros_like(jSolution,dtype = complex) for i, src in enumerate(srcList): jp = src.jPrimary(self.prob) - jPrimary[:,i] += jp + jPrimary[:,i] = jPrimary[:,i] + jp return jPrimary def _jSecondary(self, jSolution, srcList): @@ -237,7 +231,7 @@ class Fields_j(Fields): hPrimary = np.zeros([self._edgeCurl.shape[1],jSolution.shape[1]],dtype = complex) for i, src in enumerate(srcList): hp = src.hPrimary(self.prob) - hPrimary[:,i] += hp + hPrimary[:,i] = hPrimary[:,i] + hp return hPrimary def _hSecondary(self, jSolution, srcList): @@ -245,7 +239,7 @@ class Fields_j(Fields): for i, src in enumerate(srcList): h[:,i] *= -1./(1j*omega(src.freq)) S_m,_ = src.eval(self.prob) - h[:,i] += 1./(1j*omega(src.freq)) * self._MeMuI * (S_m) + h[:,i] = h[:,i]+ 1./(1j*omega(src.freq)) * self._MeMuI * (S_m) return h def _hSecondaryDeriv_u(self, src, v, adjoint=False): @@ -271,10 +265,10 @@ class Fields_j(Fields): if not adjoint: S_mDeriv = S_mDeriv(v) - hDeriv_m += 1./(1j*omega(src.freq)) * MeMuI * (Me * S_mDeriv) + hDeriv_m = hDeriv_m + 1./(1j*omega(src.freq)) * MeMuI * (Me * S_mDeriv) elif adjoint: S_mDeriv = S_mDeriv(Me.T * (MeMuI.T * v)) - hDeriv_m += 1./(1j*omega(src.freq)) * S_mDeriv + hDeriv_m = hDeriv_m + 1./(1j*omega(src.freq)) * S_mDeriv return hDeriv_m @@ -313,7 +307,7 @@ class Fields_h(Fields): hPrimary = np.zeros_like(hSolution,dtype = complex) for i, src in enumerate(srcList): hp = src.hPrimary(self.prob) - hPrimary[:,i] += hp + hPrimary[:,i] = hPrimary[:,i] + hp return hPrimary def _hSecondary(self, hSolution, srcList): @@ -333,14 +327,14 @@ class Fields_h(Fields): jPrimary = np.zeros([self._edgeCurl.shape[0], hSolution.shape[1]], dtype = complex) for i, src in enumerate(srcList): jp = src.jPrimary(self.prob) - jPrimary[:,i] += jp + jPrimary[:,i] = jPrimary[:,i] + jp return jPrimary def _jSecondary(self, hSolution, srcList): j = self._edgeCurl*hSolution for i, src in enumerate(srcList): _,S_e = src.eval(self.prob) - j[:,i] += -S_e + j[:,i] = j[:,i]+ -S_e return j def _jSecondaryDeriv_u(self, src, v, adjoint=False): From 71a5c8b2c1743cf8a892c3d3d59404c46bc4c947 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 24 Nov 2015 11:01:51 -0800 Subject: [PATCH 12/14] fix sizes in Fields object (make an array not a vector) --- SimPEG/EM/FDEM/FieldsFDEM.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/SimPEG/EM/FDEM/FieldsFDEM.py b/SimPEG/EM/FDEM/FieldsFDEM.py index f90dbfef..8f6fafe9 100644 --- a/SimPEG/EM/FDEM/FieldsFDEM.py +++ b/SimPEG/EM/FDEM/FieldsFDEM.py @@ -159,7 +159,7 @@ class Fields_b(Fields): Me = Me.T w = self._edgeCurl.T * (self._MfMui * bSolution) - w += -Utils.mkvc(Me * S_e,2) + w = w - Utils.mkvc(Me * S_e,2) if not adjoint: de_dm = self._MeSigmaIDeriv(w) * v @@ -169,8 +169,7 @@ class Fields_b(Fields): _, S_eDeriv = src.evalDeriv(self.prob, adjoint) Se_Deriv = S_eDeriv(v) - if Se_Deriv is not None: - de_dm = de_dm-self._MeSigmaI * Se_Deriv + de_dm = de_dm - self._MeSigmaI * Se_Deriv return de_dm From 5bd9209d6ab3aaa93ca4d59d4a552f312093a1dc Mon Sep 17 00:00:00 2001 From: Lindsey Date: Tue, 24 Nov 2015 14:38:58 -0800 Subject: [PATCH 13/14] type cast b in solver so that it is not an object if we add zero --- SimPEG/Utils/SolverUtils.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/SimPEG/Utils/SolverUtils.py b/SimPEG/Utils/SolverUtils.py index 279d2b06..26ff3e2a 100644 --- a/SimPEG/Utils/SolverUtils.py +++ b/SimPEG/Utils/SolverUtils.py @@ -37,12 +37,20 @@ def SolverWrapD(fun, factorize=True, checkAccuracy=True, accuracyTol=1e-6): if len(b.shape) == 1 or b.shape[1] == 1: b = b.flatten() # Just one RHS + + if b.dtype is np.dtype('O'): + b = b.astype(type(b[0])) + if factorize: X = self.solver.solve(b, **self.kwargs) else: X = fun(self.A, b, **self.kwargs) else: # Multiple RHSs + if b.dtype is np.dtype('O'): + b = b.astype(type(b[0,0])) + X = np.empty_like(b) + for i in range(b.shape[1]): if factorize: X[:,i] = self.solver.solve(b[:,i]) From 37f56f25f362a9546e49eedeceae8414d5d5b8d6 Mon Sep 17 00:00:00 2001 From: Lindsey Date: Tue, 24 Nov 2015 16:15:57 -0800 Subject: [PATCH 14/14] cleaned out Zero, Identity import from FDEM.py --- SimPEG/EM/FDEM/FDEM.py | 1 - 1 file changed, 1 deletion(-) diff --git a/SimPEG/EM/FDEM/FDEM.py b/SimPEG/EM/FDEM/FDEM.py index eec4bab9..f2167fd8 100644 --- a/SimPEG/EM/FDEM/FDEM.py +++ b/SimPEG/EM/FDEM/FDEM.py @@ -3,7 +3,6 @@ from scipy.constants import mu_0 from SurveyFDEM import Survey as SurveyFDEM from FieldsFDEM import Fields, Fields_e, Fields_b, Fields_h, Fields_j from SimPEG.EM.Base import BaseEMProblem -from SimPEG.Utils import Zero, Identity from SimPEG.EM.Utils import omega