From d102ec9c00952f6e49f8b169080626d9cfd7a5fa Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Fri, 6 Nov 2015 10:00:30 -0800 Subject: [PATCH 01/28] 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/28] 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/28] 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/28] 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/28] 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/28] 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/28] 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 815311bfecf92d65dfc7d5cbd7d96d821843675f Mon Sep 17 00:00:00 2001 From: Lindsey Date: Mon, 23 Nov 2015 18:03:49 -0800 Subject: [PATCH 08/28] 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 09/28] 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 10/28] 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 11/28] 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 12/28] 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 13/28] 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 From 83e29b74c721b6358867dc8e1f7552ac765c110a Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Tue, 24 Nov 2015 19:28:06 -0800 Subject: [PATCH 14/28] Seed the random flowTests --- tests/flow/test_Richards.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tests/flow/test_Richards.py b/tests/flow/test_Richards.py index 6579370d..d63a6210 100644 --- a/tests/flow/test_Richards.py +++ b/tests/flow/test_Richards.py @@ -12,6 +12,8 @@ except Exception, e: TOL = 1E-8 +np.random.seed(0) + class TestModels(unittest.TestCase): def test_BaseHaverkamp_Theta(self): From 56d5019b945d267a940e9aeb33efa0c95f3136ef Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Tue, 24 Nov 2015 22:09:50 -0800 Subject: [PATCH 15/28] Updates to examples and documentation. --- ...DCfwd.py => Forward_BasicDirectCurrent.py} | 4 ++- .../{Linear.py => Inversion_Linear.py} | 23 +++++++------ SimPEG/Examples/Mesh_QuadTree_Create.py | 18 +++++++++++ SimPEG/Examples/Mesh_ThreeMeshes.py | 24 ++++++++++++++ SimPEG/Examples/__init__.py | 10 +++++- SimPEG/Mesh/TreeMesh.py | 12 ++++--- SimPEG/Utils/__init__.py | 1 + SimPEG/Utils/codeutils.py | 32 +++++++++++++++++++ docs/api_Examples.rst | 11 +++++-- docs/api_Mesh.rst | 18 ++--------- docs/examples/Forward_BasicDirectCurrent.rst | 17 ++++++++++ docs/examples/Inversion_Linear.rst | 17 ++++++++++ docs/examples/Mesh_QuadTree_Create.rst | 17 ++++++++++ docs/examples/Mesh_ThreeMeshes.rst | 17 ++++++++++ tests/examples/test_examples.py | 21 +++++++----- 15 files changed, 201 insertions(+), 41 deletions(-) rename SimPEG/Examples/{DCfwd.py => Forward_BasicDirectCurrent.py} (98%) rename SimPEG/Examples/{Linear.py => Inversion_Linear.py} (78%) create mode 100644 SimPEG/Examples/Mesh_QuadTree_Create.py create mode 100644 SimPEG/Examples/Mesh_ThreeMeshes.py create mode 100644 docs/examples/Forward_BasicDirectCurrent.rst create mode 100644 docs/examples/Inversion_Linear.rst create mode 100644 docs/examples/Mesh_QuadTree_Create.rst create mode 100644 docs/examples/Mesh_ThreeMeshes.rst diff --git a/SimPEG/Examples/DCfwd.py b/SimPEG/Examples/Forward_BasicDirectCurrent.py similarity index 98% rename from SimPEG/Examples/DCfwd.py rename to SimPEG/Examples/Forward_BasicDirectCurrent.py index 33e7aad5..a3c9ff97 100644 --- a/SimPEG/Examples/DCfwd.py +++ b/SimPEG/Examples/Forward_BasicDirectCurrent.py @@ -12,7 +12,6 @@ def run(plotIt=True): tM = Mesh.TensorMesh(sz) # Curvilinear Mesh rM = Mesh.CurvilinearMesh(Utils.meshutils.exampleLrmGrid(sz,'rotate')) - # Step2: Direct Current (DC) operator def DCfun(mesh, pts): D = mesh.faceDiv @@ -39,6 +38,7 @@ def run(plotIt=True): phirM = AinvrM*rhsrM if not plotIt: return + #Step4: Making Figure fig, axes = plt.subplots(1,2,figsize=(12*1.2,4*1.2)) label = ["(a)", "(b)"] @@ -69,7 +69,9 @@ def run(plotIt=True): else: axes[i].set_ylabel(" ") axes[i].set_xlabel("x") + plt.show() if __name__ == '__main__': + Utils._makeExample(__file__) run() diff --git a/SimPEG/Examples/Linear.py b/SimPEG/Examples/Inversion_Linear.py similarity index 78% rename from SimPEG/Examples/Linear.py rename to SimPEG/Examples/Inversion_Linear.py index b065682b..2446bd58 100644 --- a/SimPEG/Examples/Linear.py +++ b/SimPEG/Examples/Inversion_Linear.py @@ -23,7 +23,9 @@ class LinearProblem(Problem.BaseProblem): return self.G.T.dot(v) -def run(N, plotIt=True): +def run(N=100, plotIt=True): + np.random.seed(1) + mesh = Mesh.TensorMesh([N]) nk = 20 @@ -52,7 +54,7 @@ def run(N, plotIt=True): reg = Regularization.Tikhonov(mesh) dmis = DataMisfit.l2_DataMisfit(survey) - opt = Optimization.InexactGaussNewton(maxIter=20) + opt = Optimization.InexactGaussNewton(maxIter=35) invProb = InvProblem.BaseInvProblem(dmis, reg, opt) beta = Directives.BetaSchedule() betaest = Directives.BetaEstimate_ByEig() @@ -63,16 +65,19 @@ def run(N, plotIt=True): if plotIt: import matplotlib.pyplot as plt - plt.figure(1) - for i in range(prob.G.shape[0]): - plt.plot(prob.G[i,:]) - plt.figure(2) - plt.plot(M.vectorCCx, survey.mtrue, 'b-') - plt.plot(M.vectorCCx, mrec, 'r-') + fig, axes = plt.subplots(1,2,figsize=(12*1.2,4*1.2)) + for i in range(prob.G.shape[0]): + axes[0].plot(prob.G[i,:]) + axes[0].set_title('Columns of matrix G') + + axes[1].plot(M.vectorCCx, survey.mtrue, 'b-') + axes[1].plot(M.vectorCCx, mrec, 'r-') + axes[1].legend(('True Model', 'Recovered Model')) plt.show() return prob, survey, mesh, mrec if __name__ == '__main__': - run(100) + Utils._makeExample(__file__) + run() diff --git a/SimPEG/Examples/Mesh_QuadTree_Create.py b/SimPEG/Examples/Mesh_QuadTree_Create.py new file mode 100644 index 00000000..f4417855 --- /dev/null +++ b/SimPEG/Examples/Mesh_QuadTree_Create.py @@ -0,0 +1,18 @@ +from SimPEG import * + +def run(plotIt=True): + from SimPEG import Mesh, np + M = Mesh.TreeMesh([32,32]) + M.refine(3) + def function(cell): + xyz = cell.center + for i in range(3): + if np.abs(np.sin(xyz[0]*np.pi*2)*0.5 + 0.5 - xyz[1]) < 0.2*i: + return 6-i + return 0 + M.refine(function); + if plotIt: M.plotGrid(showIt=True) + +if __name__ == '__main__': + Utils._makeExample(__file__) + run() diff --git a/SimPEG/Examples/Mesh_ThreeMeshes.py b/SimPEG/Examples/Mesh_ThreeMeshes.py new file mode 100644 index 00000000..2fc7fa8d --- /dev/null +++ b/SimPEG/Examples/Mesh_ThreeMeshes.py @@ -0,0 +1,24 @@ +from SimPEG import * + +def run(plotIt=True): + sz = [16,16] + tM = Mesh.TensorMesh(sz) + qM = Mesh.TreeMesh(sz) + qM.refine(lambda cell: 4 if np.sqrt(((np.r_[cell.center]-0.5)**2).sum()) < 0.4 else 3) + rM = Mesh.CurvilinearMesh(Utils.meshutils.exampleLrmGrid(sz,'rotate')) + + if plotIt: + import matplotlib.pyplot as plt + fig, axes = plt.subplots(1,3,figsize=(14,5)) + opts = {} + tM.plotGrid(ax=axes[0], **opts) + axes[0].set_title('TensorMesh') + qM.plotGrid(ax=axes[1], **opts) + axes[1].set_title('TreeMesh') + rM.plotGrid(ax=axes[2], **opts) + axes[2].set_title('CurvilinearMesh') + plt.show() + +if __name__ == '__main__': + Utils._makeExample(__file__) + run() diff --git a/SimPEG/Examples/__init__.py b/SimPEG/Examples/__init__.py index 3fec0b99..53b1f9a8 100644 --- a/SimPEG/Examples/__init__.py +++ b/SimPEG/Examples/__init__.py @@ -1 +1,9 @@ -import Linear, DCfwd +# This will import everything in the directory into this file +from os import path as p +from glob import glob +__all__ = [] +for x in glob(p.join(p.dirname(__file__), '*.py')): + if not p.basename(x).startswith('__'): + __import__(p.basename(x)[:-3], globals(), locals()) + __all__ += [p.basename(x)] +del glob, p, x diff --git a/SimPEG/Mesh/TreeMesh.py b/SimPEG/Mesh/TreeMesh.py index 1e6c91c0..69a44fdc 100644 --- a/SimPEG/Mesh/TreeMesh.py +++ b/SimPEG/Mesh/TreeMesh.py @@ -1960,7 +1960,7 @@ class TreeMesh(BaseTensorMesh, InnerProducts): def plotGrid(self, ax=None, showIt=False, grid=True, - cells=True, cellLine=False, + cells=False, cellLine=False, nodes=False, facesX=False, facesY=False, facesZ=False, edgesX=False, edgesY=False, edgesZ=False): @@ -2009,6 +2009,8 @@ class TreeMesh(BaseTensorMesh, InnerProducts): if facesY: ax.plot(self._gridFy[self._hangingFy.keys(),0], self._gridFy[self._hangingFy.keys(),1], 'gs', ms=10, mfc='none', mec='g') ax.plot(self._gridFy[:,0], self._gridFy[:,1], 'g^') + ax.set_xlabel('x1') + ax.set_ylabel('x2') elif self.dim == 3: if cells: ax.plot(self.gridCC[:,0], self.gridCC[:,1], 'r.', zs=self.gridCC[:,2]) @@ -2056,7 +2058,6 @@ class TreeMesh(BaseTensorMesh, InnerProducts): ind = [key, hf[0]] ax.plot(self._gridEx[ind,0], self._gridEx[ind,1], 'k:', zs=self._gridEx[ind,2]) - if edgesY: ax.plot(self._gridEy[:,0], self._gridEy[:,1], 'k<', zs=self._gridEy[:,2]) ax.plot(self._gridEy[self._hangingEy.keys(),0], self._gridEy[self._hangingEy.keys(),1], 'ks', ms=10, mfc='none', mec='k', zs=self._gridEy[self._hangingEy.keys(),2]) @@ -2072,7 +2073,10 @@ class TreeMesh(BaseTensorMesh, InnerProducts): for hf in self._hangingEz[key]: ind = [key, hf[0]] ax.plot(self._gridEz[ind,0], self._gridEz[ind,1], 'k:', zs=self._gridEz[ind,2]) - + ax.set_xlabel('x1') + ax.set_ylabel('x2') + ax.set_zlabel('x3') + ax.grid(True) if showIt:plt.show() def plotImage(self, I, ax=None, showIt=True, grid=False): @@ -2191,7 +2195,7 @@ class Cell(object): @property def center(self): if getattr(self, '_center', None) is None: - self._center = self.mesh._cellC(self._pointer) + self._center = np.array(self.mesh._cellC(self._pointer)) return self._center @property def h(self): return self.mesh._cellH(self._pointer) diff --git a/SimPEG/Utils/__init__.py b/SimPEG/Utils/__init__.py index 5280ae79..67e1b117 100644 --- a/SimPEG/Utils/__init__.py +++ b/SimPEG/Utils/__init__.py @@ -1,5 +1,6 @@ from matutils import * from codeutils import * +from codeutils import _makeExample from meshutils import exampleLrmGrid, meshTensor, closestPoints, readUBCTensorMesh, writeUBCTensorMesh, writeUBCTensorModel, readVTRFile, writeVTRFile from curvutils import volTetra, faceInfo, indexCube from interputils import interpmat diff --git a/SimPEG/Utils/codeutils.py b/SimPEG/Utils/codeutils.py index 4a9a28a7..fef34282 100644 --- a/SimPEG/Utils/codeutils.py +++ b/SimPEG/Utils/codeutils.py @@ -227,3 +227,35 @@ def requires(var): return requiresVarWrapper return requiresVar + +def _makeExample(filePath): + + import os + name = filePath.split(os.path.sep)[-1][:-3] + out = """.. _examples_%s: + +.. ------------------------------ .. +.. THIS FILE IS AUTO GENEREATED .. +.. ------------------------------ .. + +%s +%s + +.. plot:: + + from SimPEG import Examples + Examples.%s.run() + +.. literalinclude:: ../../SimPEG/Examples/%s.py + :language: python + :linenos: +"""%(name,name.replace('_',' '),'='*len(name),name,name) + + rst = os.path.sep.join((filePath.split(os.path.sep)[:-3] + ['docs', 'examples', name + '.rst'])) + + f = open(rst, 'w') + f.write(out) + f.close() + + + diff --git a/docs/api_Examples.rst b/docs/api_Examples.rst index 214dd8e1..68589fb3 100644 --- a/docs/api_Examples.rst +++ b/docs/api_Examples.rst @@ -3,8 +3,15 @@ Examples ******** -Forward problem -=============== +.. toctree:: + :maxdepth: 1 + :glob: + + examples/* + + +External Notebooks +================== * `Example 1: Direct Current `_ * `Example 2: Seismic-Acoustic `_ diff --git a/docs/api_Mesh.rst b/docs/api_Mesh.rst index 7bf398b1..a7f7abae 100644 --- a/docs/api_Mesh.rst +++ b/docs/api_Mesh.rst @@ -23,23 +23,9 @@ the implementations. .. plot:: - from SimPEG import Mesh, Utils, np - import matplotlib.pyplot as plt - sz = [10,10] - tM = Mesh.TensorMesh(sz) - qM = Mesh.TreeMesh(sz) - qM.refine(lambda X: 1 if np.sqrt(((X-0.5)**2).sum()) < 0.3 else 0) - rM = Mesh.CurvilinearMesh(Utils.meshutils.exampleLrmGrid(sz,'rotate')) + from SimPEG import Examples + Examples.Mesh_ThreeMeshes.run() - fig, axes = plt.subplots(1,3,figsize=(14,5)) - opts = {} - tM.plotGrid(ax=axes[0], **opts) - axes[0].set_title('TensorMesh') - qM.plotGrid(ax=axes[1], **opts) - axes[1].set_title('TreeMesh') - rM.plotGrid(ax=axes[2], **opts) - axes[2].set_title('CurvilinearMesh') - plt.show() Variable Locations and Terminology diff --git a/docs/examples/Forward_BasicDirectCurrent.rst b/docs/examples/Forward_BasicDirectCurrent.rst new file mode 100644 index 00000000..6fc816e7 --- /dev/null +++ b/docs/examples/Forward_BasicDirectCurrent.rst @@ -0,0 +1,17 @@ +.. _examples_Forward_BasicDirectCurrent: + +.. ------------------------------ .. +.. THIS FILE IS AUTO GENEREATED .. +.. ------------------------------ .. + +Forward BasicDirectCurrent +========================== + +.. plot:: + + from SimPEG import Examples + Examples.Forward_BasicDirectCurrent.run() + +.. literalinclude:: ../../SimPEG/Examples/Forward_BasicDirectCurrent.py + :language: python + :linenos: diff --git a/docs/examples/Inversion_Linear.rst b/docs/examples/Inversion_Linear.rst new file mode 100644 index 00000000..637315ae --- /dev/null +++ b/docs/examples/Inversion_Linear.rst @@ -0,0 +1,17 @@ +.. _examples_Inversion_Linear: + +.. ------------------------------ .. +.. THIS FILE IS AUTO GENEREATED .. +.. ------------------------------ .. + +Inversion Linear +================ + +.. plot:: + + from SimPEG import Examples + Examples.Inversion_Linear.run() + +.. literalinclude:: ../../SimPEG/Examples/Inversion_Linear.py + :language: python + :linenos: diff --git a/docs/examples/Mesh_QuadTree_Create.rst b/docs/examples/Mesh_QuadTree_Create.rst new file mode 100644 index 00000000..644593d5 --- /dev/null +++ b/docs/examples/Mesh_QuadTree_Create.rst @@ -0,0 +1,17 @@ +.. _examples_Mesh_QuadTree_Create: + +.. ------------------------------ .. +.. THIS FILE IS AUTO GENEREATED .. +.. ------------------------------ .. + +Mesh QuadTree Create +==================== + +.. plot:: + + from SimPEG import Examples + Examples.Mesh_QuadTree_Create.run() + +.. literalinclude:: ../../SimPEG/Examples/Mesh_QuadTree_Create.py + :language: python + :linenos: diff --git a/docs/examples/Mesh_ThreeMeshes.rst b/docs/examples/Mesh_ThreeMeshes.rst new file mode 100644 index 00000000..948398a7 --- /dev/null +++ b/docs/examples/Mesh_ThreeMeshes.rst @@ -0,0 +1,17 @@ +.. _examples_Mesh_ThreeMeshes: + +.. ------------------------------ .. +.. THIS FILE IS AUTO GENEREATED .. +.. ------------------------------ .. + +Mesh ThreeMeshes +================ + +.. plot:: + + from SimPEG import Examples + Examples.Mesh_ThreeMeshes.run() + +.. literalinclude:: ../../SimPEG/Examples/Mesh_ThreeMeshes.py + :language: python + :linenos: diff --git a/tests/examples/test_examples.py b/tests/examples/test_examples.py index 1fcf05f5..baddb6db 100644 --- a/tests/examples/test_examples.py +++ b/tests/examples/test_examples.py @@ -1,17 +1,22 @@ import unittest import sys -from SimPEG.Examples import Linear, DCfwd +from SimPEG import Examples import numpy as np -class TestLinear(unittest.TestCase): - def test_running(self): - Linear.run(100, plotIt=False) +def get_test(test): + def func(self): + print '\nTesting %s.run(plotIt=False)\n'%test + getattr(Examples, test).run(plotIt=False) self.assertTrue(True) + return func +attrs = dict() +tests = [_ for _ in dir(Examples) if not _.startswith('_')] +for test in tests: + attrs['test_'+test] = get_test(test) + + +TestExamples = type('TestExamples', (unittest.TestCase,), attrs) -class TestDCfwd(unittest.TestCase): - def test_running(self): - DCfwd.run(plotIt=False) - self.assertTrue(True) if __name__ == '__main__': unittest.main() From bcfe9040154812fe4ed995038dca865fda387ab3 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Wed, 25 Nov 2015 08:26:27 -0800 Subject: [PATCH 16/28] remove test from the name in examples. --- tests/examples/test_examples.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/examples/test_examples.py b/tests/examples/test_examples.py index baddb6db..48cbaa65 100644 --- a/tests/examples/test_examples.py +++ b/tests/examples/test_examples.py @@ -3,7 +3,7 @@ import sys from SimPEG import Examples import numpy as np -def get_test(test): +def get(test): def func(self): print '\nTesting %s.run(plotIt=False)\n'%test getattr(Examples, test).run(plotIt=False) @@ -12,8 +12,8 @@ def get_test(test): attrs = dict() tests = [_ for _ in dir(Examples) if not _.startswith('_')] for test in tests: - attrs['test_'+test] = get_test(test) - + attrs['test_'+test] = get(test) +del get, tests, _ TestExamples = type('TestExamples', (unittest.TestCase,), attrs) From 06c66f6392ce4bd2e4517833df970f99d1c92982 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Wed, 25 Nov 2015 13:27:35 -0800 Subject: [PATCH 17/28] Addresses #183: Import errors for tree mesh --- SimPEG/Mesh/TreeMesh.py | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/SimPEG/Mesh/TreeMesh.py b/SimPEG/Mesh/TreeMesh.py index 1e6c91c0..f7a438ce 100644 --- a/SimPEG/Mesh/TreeMesh.py +++ b/SimPEG/Mesh/TreeMesh.py @@ -96,7 +96,13 @@ from mpl_toolkits.mplot3d import Axes3D import matplotlib.colors as colors import matplotlib.cm as cmx -import TreeUtils +try: + import TreeUtils + _IMPORT_TREEUTILS = True +except Exception, e: + _IMPORT_TREEUTILS = False + + from InnerProducts import InnerProducts from TensorMesh import TensorMesh, BaseTensorMesh import time @@ -108,6 +114,8 @@ class TreeMesh(BaseTensorMesh, InnerProducts): _meshType = 'TREE' def __init__(self, h, x0=None, levels=None): + if not _IMPORT_TREEUTILS: + raise Exception('Could not import the Cython code to run the TreeMesh Try:.\n\npython setup.py build_ext --inplace') assert type(h) is list, 'h must be a list' assert len(h) in [2,3], "There is only support for TreeMesh in 2D or 3D." From be3667b4abee8f19cc4ae1a23663ec083bf056c9 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Wed, 25 Nov 2015 16:03:08 -0800 Subject: [PATCH 18/28] New Examples - Put all examples in same directory - Make a single test - use __init__.py to create the docs automatically --- .travis.yml | 1 - SimPEG/EM/Examples/__init__.py | 1 - .../EM_FDEM_1D_Inversion.py} | 7 +++ .../FLOW_Richards_1D_Celia1990.py} | 34 +++++++++++ SimPEG/Examples/Forward_BasicDirectCurrent.py | 1 - SimPEG/Examples/Inversion_Linear.py | 53 +++++++++------- SimPEG/Examples/Mesh_Basic_PlotImage.py | 46 ++++++++++++++ ...esh_ThreeMeshes.py => Mesh_Basic_Types.py} | 8 ++- ...ee_Create.py => Mesh_QuadTree_Creation.py} | 14 ++++- SimPEG/Examples/Mesh_QuadTree_HangingNodes.py | 32 ++++++++++ SimPEG/Examples/Mesh_Tensor_Creation.py | 35 +++++++++++ SimPEG/Examples/__init__.py | 60 ++++++++++++++++++- SimPEG/FLOW/Examples/__init__.py | 1 - SimPEG/Mesh/TreeMesh.py | 39 +++++++----- SimPEG/Utils/__init__.py | 1 - SimPEG/Utils/codeutils.py | 32 ---------- docs/examples/EM_FDEM_1D_Inversion.rst | 26 ++++++++ docs/examples/FLOW_Richards_1D_Celia1990.rst | 52 ++++++++++++++++ docs/examples/Forward_BasicDirectCurrent.rst | 10 +++- docs/examples/Inversion_Linear.rst | 19 ++++-- docs/examples/Mesh_Basic_PlotImage.rst | 27 +++++++++ docs/examples/Mesh_Basic_Types.rst | 26 ++++++++ docs/examples/Mesh_QuadTree_Create.rst | 17 ------ docs/examples/Mesh_QuadTree_Creation.rst | 31 ++++++++++ docs/examples/Mesh_QuadTree_HangingNodes.rst | 31 ++++++++++ docs/examples/Mesh_Tensor_Creation.rst | 43 +++++++++++++ docs/examples/Mesh_ThreeMeshes.rst | 17 ------ tests/em/examples/__init__.py | 11 ---- tests/em/examples/test_Examples.py | 10 ---- tests/flow/test_examples.py | 12 ---- 30 files changed, 543 insertions(+), 154 deletions(-) delete mode 100644 SimPEG/EM/Examples/__init__.py rename SimPEG/{EM/Examples/CylInversion.py => Examples/EM_FDEM_1D_Inversion.py} (96%) rename SimPEG/{FLOW/Examples/Celia1990.py => Examples/FLOW_Richards_1D_Celia1990.py} (54%) create mode 100644 SimPEG/Examples/Mesh_Basic_PlotImage.py rename SimPEG/Examples/{Mesh_ThreeMeshes.py => Mesh_Basic_Types.py} (82%) rename SimPEG/Examples/{Mesh_QuadTree_Create.py => Mesh_QuadTree_Creation.py} (50%) create mode 100644 SimPEG/Examples/Mesh_QuadTree_HangingNodes.py create mode 100644 SimPEG/Examples/Mesh_Tensor_Creation.py delete mode 100644 SimPEG/FLOW/Examples/__init__.py create mode 100644 docs/examples/EM_FDEM_1D_Inversion.rst create mode 100644 docs/examples/FLOW_Richards_1D_Celia1990.rst create mode 100644 docs/examples/Mesh_Basic_PlotImage.rst create mode 100644 docs/examples/Mesh_Basic_Types.rst delete mode 100644 docs/examples/Mesh_QuadTree_Create.rst create mode 100644 docs/examples/Mesh_QuadTree_Creation.rst create mode 100644 docs/examples/Mesh_QuadTree_HangingNodes.rst create mode 100644 docs/examples/Mesh_Tensor_Creation.rst delete mode 100644 docs/examples/Mesh_ThreeMeshes.rst delete mode 100644 tests/em/examples/__init__.py delete mode 100644 tests/em/examples/test_Examples.py delete mode 100644 tests/flow/test_examples.py diff --git a/.travis.yml b/.travis.yml index c2c672bf..b9826bab 100644 --- a/.travis.yml +++ b/.travis.yml @@ -5,7 +5,6 @@ python: sudo: false env: - - TEST_DIR=tests/em/examples - TEST_DIR=tests/em/fdem/forward - TEST_DIR=tests/em/fdem/inverse/derivs - TEST_DIR=tests/em/fdem/inverse/adjoint diff --git a/SimPEG/EM/Examples/__init__.py b/SimPEG/EM/Examples/__init__.py deleted file mode 100644 index eb36678d..00000000 --- a/SimPEG/EM/Examples/__init__.py +++ /dev/null @@ -1 +0,0 @@ -import CylInversion diff --git a/SimPEG/EM/Examples/CylInversion.py b/SimPEG/Examples/EM_FDEM_1D_Inversion.py similarity index 96% rename from SimPEG/EM/Examples/CylInversion.py rename to SimPEG/Examples/EM_FDEM_1D_Inversion.py index cfcfcfc1..fdf2750e 100644 --- a/SimPEG/EM/Examples/CylInversion.py +++ b/SimPEG/Examples/EM_FDEM_1D_Inversion.py @@ -4,6 +4,13 @@ from scipy.constants import mu_0 import matplotlib.pyplot as plt def run(plotIt=True): + """ + EM: FDEM: 1D: Inversion + ======================= + + Here we will create and run a FDEM 1D inversion. + + """ cs, ncx, ncz, npad = 5., 25, 15, 15 hx = [(cs,ncx), (cs,npad,1.3)] diff --git a/SimPEG/FLOW/Examples/Celia1990.py b/SimPEG/Examples/FLOW_Richards_1D_Celia1990.py similarity index 54% rename from SimPEG/FLOW/Examples/Celia1990.py rename to SimPEG/Examples/FLOW_Richards_1D_Celia1990.py index 24ae82a6..7c663e6b 100644 --- a/SimPEG/FLOW/Examples/Celia1990.py +++ b/SimPEG/Examples/FLOW_Richards_1D_Celia1990.py @@ -3,6 +3,39 @@ from SimPEG.FLOW import Richards import matplotlib.pyplot as plt def run(plotIt=True): + """ + FLOW: Richards: 1D: Celia1990 + ============================= + + There are two different forms of Richards equation that differ + on how they deal with the non-linearity in the time-stepping term. + + The most fundamental form, referred to as the + 'mixed'-form of Richards Equation Celia1990_ + + .. math:: + + \\frac{\partial \\theta(\psi)}{\partial t} - \\nabla \cdot k(\psi) \\nabla \psi - \\frac{\partial k(\psi)}{\partial z} = 0 + \quad \psi \in \Omega + + where \\\\(\\\\theta\\\\) is water content, and \\\\(\\\\psi\\\\) is pressure head. + This formulation of Richards equation is called the + 'mixed'-form because the equation is parameterized in \\\\(\\\\psi\\\\) + but the time-stepping is in terms of \\\\(\\\\theta\\\\). + + As noted in Celia1990_ the 'head'-based form of Richards + equation can be written in the continuous form as: + + .. math:: + + \\frac{\partial \\theta}{\partial \psi}\\frac{\partial \psi}{\partial t} - \\nabla \cdot k(\psi) \\nabla \psi - \\frac{\partial k(\psi)}{\partial z} = 0 \quad \psi \in \Omega + + However, it can be shown that this does not conserve mass in the discrete formulation. + + Here we reproduce the results from Celia1990_ demonstrating the head-based formulation and the mixed-formulation. + + .. _Celia1990: http://www.webpages.uidaho.edu/ch/papers/Celia.pdf + """ M = Mesh.TensorMesh([np.ones(40)]) M.setCellGradBC('dirichlet') params = Richards.Empirical.HaverkampParams().celia1990 @@ -47,6 +80,7 @@ def run(plotIt=True): plt.xlabel('Depth, cm') plt.ylabel('Pressure Head, cm') plt.legend(('$\Delta t$ = 10 sec','$\Delta t$ = 30 sec','$\Delta t$ = 120 sec')) + plt.show() if __name__ == '__main__': run() diff --git a/SimPEG/Examples/Forward_BasicDirectCurrent.py b/SimPEG/Examples/Forward_BasicDirectCurrent.py index a3c9ff97..06c9e79e 100644 --- a/SimPEG/Examples/Forward_BasicDirectCurrent.py +++ b/SimPEG/Examples/Forward_BasicDirectCurrent.py @@ -73,5 +73,4 @@ def run(plotIt=True): if __name__ == '__main__': - Utils._makeExample(__file__) run() diff --git a/SimPEG/Examples/Inversion_Linear.py b/SimPEG/Examples/Inversion_Linear.py index 2446bd58..a8a0eddc 100644 --- a/SimPEG/Examples/Inversion_Linear.py +++ b/SimPEG/Examples/Inversion_Linear.py @@ -1,29 +1,37 @@ from SimPEG import * -class LinearSurvey(Survey.BaseSurvey): - def projectFields(self, u): - return u - -class LinearProblem(Problem.BaseProblem): - """docstring for LinearProblem""" - - surveyPair = LinearSurvey - - def __init__(self, mesh, G, **kwargs): - Problem.BaseProblem.__init__(self, mesh, **kwargs) - self.G = G - - def fields(self, m, u=None): - return self.G.dot(m) - - def Jvec(self, m, v, u=None): - return self.G.dot(v) - - def Jtvec(self, m, v, u=None): - return self.G.T.dot(v) - def run(N=100, plotIt=True): + """ + Inversion: Linear Problem + ========================= + + Here we go over the basics of creating a linear problem and inversion. + + """ + + class LinearSurvey(Survey.BaseSurvey): + def projectFields(self, u): + return u + + class LinearProblem(Problem.BaseProblem): + + surveyPair = LinearSurvey + + def __init__(self, mesh, G, **kwargs): + Problem.BaseProblem.__init__(self, mesh, **kwargs) + self.G = G + + def fields(self, m, u=None): + return self.G.dot(m) + + def Jvec(self, m, v, u=None): + return self.G.dot(v) + + def Jtvec(self, m, v, u=None): + return self.G.T.dot(v) + + np.random.seed(1) mesh = Mesh.TensorMesh([N]) @@ -79,5 +87,4 @@ def run(N=100, plotIt=True): return prob, survey, mesh, mrec if __name__ == '__main__': - Utils._makeExample(__file__) run() diff --git a/SimPEG/Examples/Mesh_Basic_PlotImage.py b/SimPEG/Examples/Mesh_Basic_PlotImage.py new file mode 100644 index 00000000..3154f281 --- /dev/null +++ b/SimPEG/Examples/Mesh_Basic_PlotImage.py @@ -0,0 +1,46 @@ +from SimPEG import * + +def run(plotIt=True): + """ + Mesh: Basic: PlotImage + ====================== + + You can use M.PlotImage to plot images on all of the Meshes. + + + """ + M = Mesh.TensorMesh([32,32]) + v = Utils.ModelBuilder.randomModel(M.vnC, seed=789) + v = Utils.mkvc(v) + + O = Mesh.TreeMesh([32,32]) + O.refine(1) + def function(cell): + if (cell.center[0] < 0.75 and cell.center[0] > 0.25 and + cell.center[1] < 0.75 and cell.center[1] > 0.25):return 5 + if (cell.center[0] < 0.9 and cell.center[0] > 0.1 and + cell.center[1] < 0.9 and cell.center[1] > 0.1):return 4 + return 3 + O.refine(function) + + P = M.getInterpolationMat(O.gridCC, 'CC') + + ov = P * v + + if plotIt: + import matplotlib.pyplot as plt + + fig, axes = plt.subplots(1,2,figsize=(10,5)) + + out = M.plotImage(v, grid=True, ax=axes[0]) + cb = plt.colorbar(out[0], ax=axes[0]); cb.set_label("Random Field") + axes[0].set_title('TensorMesh') + + out = O.plotImage(ov, grid=True, ax=axes[1], clim=[0,1]) + cb = plt.colorbar(out[0], ax=axes[1]); cb.set_label("Random Field") + axes[1].set_title('TreeMesh') + + plt.show() + +if __name__ == '__main__': + run() diff --git a/SimPEG/Examples/Mesh_ThreeMeshes.py b/SimPEG/Examples/Mesh_Basic_Types.py similarity index 82% rename from SimPEG/Examples/Mesh_ThreeMeshes.py rename to SimPEG/Examples/Mesh_Basic_Types.py index 2fc7fa8d..430fe698 100644 --- a/SimPEG/Examples/Mesh_ThreeMeshes.py +++ b/SimPEG/Examples/Mesh_Basic_Types.py @@ -1,6 +1,13 @@ from SimPEG import * def run(plotIt=True): + """ + Mesh: Basic: Types + ================== + + Here we show SimPEG used to create three different types of meshes. + + """ sz = [16,16] tM = Mesh.TensorMesh(sz) qM = Mesh.TreeMesh(sz) @@ -20,5 +27,4 @@ def run(plotIt=True): plt.show() if __name__ == '__main__': - Utils._makeExample(__file__) run() diff --git a/SimPEG/Examples/Mesh_QuadTree_Create.py b/SimPEG/Examples/Mesh_QuadTree_Creation.py similarity index 50% rename from SimPEG/Examples/Mesh_QuadTree_Create.py rename to SimPEG/Examples/Mesh_QuadTree_Creation.py index f4417855..ede69a63 100644 --- a/SimPEG/Examples/Mesh_QuadTree_Create.py +++ b/SimPEG/Examples/Mesh_QuadTree_Creation.py @@ -1,7 +1,18 @@ from SimPEG import * def run(plotIt=True): - from SimPEG import Mesh, np + """ + Mesh: QuadTree: Creation + ======================== + + You can give the refine method a function, which is evaluated on every cell + of the TreeMesh. + + Occasionally it is useful to initially refine to a constant level + (e.g. 3 in this 32x32 mesh). This means the function is first evaluated + on an 8x8 mesh (2^3). + + """ M = Mesh.TreeMesh([32,32]) M.refine(3) def function(cell): @@ -14,5 +25,4 @@ def run(plotIt=True): if plotIt: M.plotGrid(showIt=True) if __name__ == '__main__': - Utils._makeExample(__file__) run() diff --git a/SimPEG/Examples/Mesh_QuadTree_HangingNodes.py b/SimPEG/Examples/Mesh_QuadTree_HangingNodes.py new file mode 100644 index 00000000..11726995 --- /dev/null +++ b/SimPEG/Examples/Mesh_QuadTree_HangingNodes.py @@ -0,0 +1,32 @@ +from SimPEG import * + +def run(plotIt=True): + """ + Mesh: QuadTree: Hanging Nodes + ============================= + + You can give the refine method a function, which is evaluated on every cell + of the TreeMesh. + + Occasionally it is useful to initially refine to a constant level + (e.g. 3 in this 32x32 mesh). This means the function is first evaluated + on an 8x8 mesh (2^3). + + """ + M = Mesh.TreeMesh([8,8]) + def function(cell): + xyz = cell.center + dist = ((xyz - [0.25,0.25])**2).sum()**0.5 + if dist < 0.25: + return 3 + return 2 + M.refine(function); + M.number() + if plotIt: + import matplotlib.pyplot as plt + M.plotGrid(nodes=True, cells=True, facesX=True) + plt.legend(('Grid', 'Cell Centers', 'Nodes', 'Hanging Nodes', 'X faces', 'Hanging X faces')) + plt.show() + +if __name__ == '__main__': + run() diff --git a/SimPEG/Examples/Mesh_Tensor_Creation.py b/SimPEG/Examples/Mesh_Tensor_Creation.py new file mode 100644 index 00000000..31ad3d69 --- /dev/null +++ b/SimPEG/Examples/Mesh_Tensor_Creation.py @@ -0,0 +1,35 @@ +from SimPEG import * + +def run(plotIt=True): + """ + + Mesh: Tensor: Creation + ====================== + + For tensor meshes, there are some functions that can come + in handy. For example, creating mesh tensors can be a bit time + consuming, these can be created speedily by just giving numbers + and sizes of padding. See the example below, that follows this + notation:: + + h1 = ( + (cellSize, numPad, [, increaseFactor]), + (cellSize, numCore), + (cellSize, numPad, [, increaseFactor]) + ) + + .. note:: + + You can center your mesh by passing a 'C' for the x0[i] position. + A 'N' will make the entire mesh negative, and a '0' (or a 0) will + make the mesh start at zero. + + """ + h1 = [(10, 5, -1.3), (5, 20), (10, 3, 1.3)] + M = Mesh.TensorMesh([h1, h1], x0='CN') + if plotIt: + M.plotGrid(showIt=True) + +if __name__ == '__main__': + run() + diff --git a/SimPEG/Examples/__init__.py b/SimPEG/Examples/__init__.py index 53b1f9a8..ac219a7f 100644 --- a/SimPEG/Examples/__init__.py +++ b/SimPEG/Examples/__init__.py @@ -5,5 +5,63 @@ __all__ = [] for x in glob(p.join(p.dirname(__file__), '*.py')): if not p.basename(x).startswith('__'): __import__(p.basename(x)[:-3], globals(), locals()) - __all__ += [p.basename(x)] + __all__ += [p.basename(x)[:-3]] del glob, p, x + +if __name__ == '__main__': + """ + + Run the following to create the examples documentation. + + """ + + import shutil, os + from SimPEG import Examples + + def _makeExample(filePath, runFunction): + filePath = os.path.realpath(filePath) + name = filePath.split(os.path.sep)[-1].rstrip('.pyc').rstrip('.py') + + docstr = runFunction.__doc__ + if docstr is None: + doc = '%s\n%s'%(name.replace('_',' '),'='*len(name)) + else: + doc = '\n'.join([_[8:].rstrip() for _ in docstr.split('\n')]) + + out = """.. _examples_%s: + +.. --------------------------------- .. +.. .. +.. THIS FILE IS AUTO GENEREATED .. +.. .. +.. SimPEG/Examples/__init__.py .. +.. .. +.. --------------------------------- .. + +%s + +.. plot:: + + from SimPEG import Examples + Examples.%s.run() + +.. literalinclude:: ../../SimPEG/Examples/%s.py + :language: python + :linenos: +"""%(name,doc,name,name) + + rst = os.path.sep.join((filePath.split(os.path.sep)[:-3] + ['docs', 'examples', name + '.rst'])) + + f = open(rst, 'w') + f.write(out) + f.close() + + + docExamplesDir = os.path.sep.join(os.path.realpath(__file__).split(os.path.sep)[:-3] + ['docs', 'examples']) + shutil.rmtree(docExamplesDir) + os.makedirs(docExamplesDir) + + for ex in dir(Examples): + if ex.startswith('_'): continue + E = getattr(Examples,ex) + _makeExample(E.__file__, E.run) diff --git a/SimPEG/FLOW/Examples/__init__.py b/SimPEG/FLOW/Examples/__init__.py deleted file mode 100644 index 7a894e11..00000000 --- a/SimPEG/FLOW/Examples/__init__.py +++ /dev/null @@ -1 +0,0 @@ -import Celia1990 diff --git a/SimPEG/Mesh/TreeMesh.py b/SimPEG/Mesh/TreeMesh.py index 69a44fdc..263909e4 100644 --- a/SimPEG/Mesh/TreeMesh.py +++ b/SimPEG/Mesh/TreeMesh.py @@ -1975,24 +1975,28 @@ class TreeMesh(BaseTensorMesh, InnerProducts): fig = ax.figure if grid: + X, Y, Z = [], [], [] for ind in self._sortedCells: p = self._asPointer(ind) n = self._cellN(p) h = self._cellH(p) - x = [n[0] , n[0] + h[0], n[0] + h[0], n[0] , n[0]] - y = [n[1] , n[1] , n[1] + h[1], n[1] + h[1], n[1]] if self.dim == 2: - ax.plot(x,y, 'b-') + X += [n[0] , n[0] + h[0], n[0] + h[0], n[0] , n[0], np.nan] + Y += [n[1] , n[1] , n[1] + h[1], n[1] + h[1], n[1], np.nan] elif self.dim == 3: - ax.plot(x,y, 'b-', zs=[n[2]]*5) - z = [n[2] + h[2], n[2] + h[2], n[2] + h[2], n[2] + h[2], n[2] + h[2]] - ax.plot(x,y, 'b-', zs=z) + X += [n[0] , n[0] + h[0], n[0] + h[0], n[0] , n[0], np.nan]*2 + Y += [n[1] , n[1] , n[1] + h[1], n[1] + h[1], n[1], np.nan]*2 + Z += [n[2]]*5+[np.nan] + Z += [n[2] + h[2], n[2] + h[2], n[2] + h[2], n[2] + h[2], n[2] + h[2], np.nan] sides = [0,0], [h[0],0], [0,h[1]], [h[0],h[1]] for s in sides: - x = [n[0] + s[0], n[0] + s[0]] - y = [n[1] + s[1], n[1] + s[1]] - z = [n[2] , n[2] + h[2]] - ax.plot(x,y, 'b-', zs=z) + X += [n[0] + s[0], n[0] + s[0]] + Y += [n[1] + s[1], n[1] + s[1]] + Z += [n[2] , n[2] + h[2]] + if self.dim == 2: + ax.plot(X,Y, 'b-') + elif self.dim == 3: + ax.plot(X,Y, 'b-', zs=Z) if self.dim == 2: if cells: @@ -2004,11 +2008,11 @@ class TreeMesh(BaseTensorMesh, InnerProducts): ax.plot(self._gridN[:,0], self._gridN[:,1], 'ms') ax.plot(self._gridN[self._hangingN.keys(),0], self._gridN[self._hangingN.keys(),1], 'ms', ms=10, mfc='none', mec='m') if facesX: - ax.plot(self._gridFx[self._hangingFx.keys(),0], self._gridFx[self._hangingFx.keys(),1], 'gs', ms=10, mfc='none', mec='g') ax.plot(self._gridFx[:,0], self._gridFx[:,1], 'g>') + ax.plot(self._gridFx[self._hangingFx.keys(),0], self._gridFx[self._hangingFx.keys(),1], 'gs', ms=10, mfc='none', mec='g') if facesY: - ax.plot(self._gridFy[self._hangingFy.keys(),0], self._gridFy[self._hangingFy.keys(),1], 'gs', ms=10, mfc='none', mec='g') ax.plot(self._gridFy[:,0], self._gridFy[:,1], 'g^') + ax.plot(self._gridFy[self._hangingFy.keys(),0], self._gridFy[self._hangingFy.keys(),1], 'gs', ms=10, mfc='none', mec='g') ax.set_xlabel('x1') ax.set_ylabel('x2') elif self.dim == 3: @@ -2079,12 +2083,15 @@ class TreeMesh(BaseTensorMesh, InnerProducts): ax.grid(True) if showIt:plt.show() - def plotImage(self, I, ax=None, showIt=True, grid=False): + def plotImage(self, I, ax=None, showIt=False, grid=False, clim=None): if self.dim == 3: raise Exception('Use plot slice?') if ax is None: ax = plt.subplot(111) jet = cm = plt.get_cmap('jet') - cNorm = colors.Normalize(vmin=I.min(), vmax=I.max()) + cNorm = colors.Normalize( + vmin=I.min() if clim is None else clim[0], + vmax=I.max() if clim is None else clim[1]) + scalarMap = cmx.ScalarMappable(norm=cNorm, cmap=jet) ax.set_xlim((self.x0[0], self.h[0].sum())) ax.set_ylim((self.x0[1], self.h[1].sum())) @@ -2093,8 +2100,10 @@ class TreeMesh(BaseTensorMesh, InnerProducts): ax.add_patch(plt.Rectangle((x0[0], x0[1]), sz[0], sz[1], facecolor=scalarMap.to_rgba(I[ii]), edgecolor='k' if grid else 'none')) # if text: ax.text(self.center[0],self.center[1],self.num) scalarMap._A = [] # http://stackoverflow.com/questions/8342549/matplotlib-add-colorbar-to-a-sequence-of-line-plots - plt.colorbar(scalarMap) + ax.set_xlabel('x') + ax.set_ylabel('y') if showIt: plt.show() + return [scalarMap] def plotSlice(self, v, vType='CC', normal='Z', ind=None, grid=True, view='real', diff --git a/SimPEG/Utils/__init__.py b/SimPEG/Utils/__init__.py index 67e1b117..5280ae79 100644 --- a/SimPEG/Utils/__init__.py +++ b/SimPEG/Utils/__init__.py @@ -1,6 +1,5 @@ from matutils import * from codeutils import * -from codeutils import _makeExample from meshutils import exampleLrmGrid, meshTensor, closestPoints, readUBCTensorMesh, writeUBCTensorMesh, writeUBCTensorModel, readVTRFile, writeVTRFile from curvutils import volTetra, faceInfo, indexCube from interputils import interpmat diff --git a/SimPEG/Utils/codeutils.py b/SimPEG/Utils/codeutils.py index fef34282..4a9a28a7 100644 --- a/SimPEG/Utils/codeutils.py +++ b/SimPEG/Utils/codeutils.py @@ -227,35 +227,3 @@ def requires(var): return requiresVarWrapper return requiresVar - -def _makeExample(filePath): - - import os - name = filePath.split(os.path.sep)[-1][:-3] - out = """.. _examples_%s: - -.. ------------------------------ .. -.. THIS FILE IS AUTO GENEREATED .. -.. ------------------------------ .. - -%s -%s - -.. plot:: - - from SimPEG import Examples - Examples.%s.run() - -.. literalinclude:: ../../SimPEG/Examples/%s.py - :language: python - :linenos: -"""%(name,name.replace('_',' '),'='*len(name),name,name) - - rst = os.path.sep.join((filePath.split(os.path.sep)[:-3] + ['docs', 'examples', name + '.rst'])) - - f = open(rst, 'w') - f.write(out) - f.close() - - - diff --git a/docs/examples/EM_FDEM_1D_Inversion.rst b/docs/examples/EM_FDEM_1D_Inversion.rst new file mode 100644 index 00000000..acbc8cdc --- /dev/null +++ b/docs/examples/EM_FDEM_1D_Inversion.rst @@ -0,0 +1,26 @@ +.. _examples_EM_FDEM_1D_Inversion: + +.. --------------------------------- .. +.. .. +.. THIS FILE IS AUTO GENEREATED .. +.. .. +.. SimPEG/Examples/__init__.py .. +.. .. +.. --------------------------------- .. + + +EM: FDEM: 1D: Inversion +======================= + +Here we will create and run a FDEM 1D inversion. + + + +.. plot:: + + from SimPEG import Examples + Examples.EM_FDEM_1D_Inversion.run() + +.. literalinclude:: ../../SimPEG/Examples/EM_FDEM_1D_Inversion.py + :language: python + :linenos: diff --git a/docs/examples/FLOW_Richards_1D_Celia1990.rst b/docs/examples/FLOW_Richards_1D_Celia1990.rst new file mode 100644 index 00000000..d2e01c13 --- /dev/null +++ b/docs/examples/FLOW_Richards_1D_Celia1990.rst @@ -0,0 +1,52 @@ +.. _examples_FLOW_Richards_1D_Celia1990: + +.. --------------------------------- .. +.. .. +.. THIS FILE IS AUTO GENEREATED .. +.. .. +.. SimPEG/Examples/__init__.py .. +.. .. +.. --------------------------------- .. + + +FLOW: Richards: 1D: Celia1990 +============================= + +There are two different forms of Richards equation that differ +on how they deal with the non-linearity in the time-stepping term. + +The most fundamental form, referred to as the +'mixed'-form of Richards Equation Celia1990_ + +.. math:: + + \frac{\partial \theta(\psi)}{\partial t} - \nabla \cdot k(\psi) \nabla \psi - \frac{\partial k(\psi)}{\partial z} = 0 + \quad \psi \in \Omega + +where \\(\\theta\\) is water content, and \\(\\psi\\) is pressure head. +This formulation of Richards equation is called the +'mixed'-form because the equation is parameterized in \\(\\psi\\) +but the time-stepping is in terms of \\(\\theta\\). + +As noted in Celia1990_ the 'head'-based form of Richards +equation can be written in the continuous form as: + +.. math:: + + \frac{\partial \theta}{\partial \psi}\frac{\partial \psi}{\partial t} - \nabla \cdot k(\psi) \nabla \psi - \frac{\partial k(\psi)}{\partial z} = 0 \quad \psi \in \Omega + +However, it can be shown that this does not conserve mass in the discrete formulation. + +Here we reproduce the results from Celia1990_ demonstrating the head-based formulation and the mixed-formulation. + +.. _Celia1990: http://www.webpages.uidaho.edu/ch/papers/Celia.pdf + + +.. plot:: + + from SimPEG import Examples + Examples.FLOW_Richards_1D_Celia1990.run() + +.. literalinclude:: ../../SimPEG/Examples/FLOW_Richards_1D_Celia1990.py + :language: python + :linenos: diff --git a/docs/examples/Forward_BasicDirectCurrent.rst b/docs/examples/Forward_BasicDirectCurrent.rst index 6fc816e7..20b39eb8 100644 --- a/docs/examples/Forward_BasicDirectCurrent.rst +++ b/docs/examples/Forward_BasicDirectCurrent.rst @@ -1,8 +1,12 @@ .. _examples_Forward_BasicDirectCurrent: -.. ------------------------------ .. -.. THIS FILE IS AUTO GENEREATED .. -.. ------------------------------ .. +.. --------------------------------- .. +.. .. +.. THIS FILE IS AUTO GENEREATED .. +.. .. +.. SimPEG/Examples/__init__.py .. +.. .. +.. --------------------------------- .. Forward BasicDirectCurrent ========================== diff --git a/docs/examples/Inversion_Linear.rst b/docs/examples/Inversion_Linear.rst index 637315ae..d635d8e1 100644 --- a/docs/examples/Inversion_Linear.rst +++ b/docs/examples/Inversion_Linear.rst @@ -1,11 +1,20 @@ .. _examples_Inversion_Linear: -.. ------------------------------ .. -.. THIS FILE IS AUTO GENEREATED .. -.. ------------------------------ .. +.. --------------------------------- .. +.. .. +.. THIS FILE IS AUTO GENEREATED .. +.. .. +.. SimPEG/Examples/__init__.py .. +.. .. +.. --------------------------------- .. + + +Inversion: Linear Problem +========================= + +Here we go over the basics of creating a linear problem and inversion. + -Inversion Linear -================ .. plot:: diff --git a/docs/examples/Mesh_Basic_PlotImage.rst b/docs/examples/Mesh_Basic_PlotImage.rst new file mode 100644 index 00000000..a730f303 --- /dev/null +++ b/docs/examples/Mesh_Basic_PlotImage.rst @@ -0,0 +1,27 @@ +.. _examples_Mesh_Basic_PlotImage: + +.. --------------------------------- .. +.. .. +.. THIS FILE IS AUTO GENEREATED .. +.. .. +.. SimPEG/Examples/__init__.py .. +.. .. +.. --------------------------------- .. + + +Mesh: Basic: PlotImage +====================== + +You can use M.PlotImage to plot images on all of the Meshes. + + + + +.. plot:: + + from SimPEG import Examples + Examples.Mesh_Basic_PlotImage.run() + +.. literalinclude:: ../../SimPEG/Examples/Mesh_Basic_PlotImage.py + :language: python + :linenos: diff --git a/docs/examples/Mesh_Basic_Types.rst b/docs/examples/Mesh_Basic_Types.rst new file mode 100644 index 00000000..9bbce0e8 --- /dev/null +++ b/docs/examples/Mesh_Basic_Types.rst @@ -0,0 +1,26 @@ +.. _examples_Mesh_Basic_Types: + +.. --------------------------------- .. +.. .. +.. THIS FILE IS AUTO GENEREATED .. +.. .. +.. SimPEG/Examples/__init__.py .. +.. .. +.. --------------------------------- .. + + +Mesh: Basic: Types +================== + +Here we show SimPEG used to create three different types of meshes. + + + +.. plot:: + + from SimPEG import Examples + Examples.Mesh_Basic_Types.run() + +.. literalinclude:: ../../SimPEG/Examples/Mesh_Basic_Types.py + :language: python + :linenos: diff --git a/docs/examples/Mesh_QuadTree_Create.rst b/docs/examples/Mesh_QuadTree_Create.rst deleted file mode 100644 index 644593d5..00000000 --- a/docs/examples/Mesh_QuadTree_Create.rst +++ /dev/null @@ -1,17 +0,0 @@ -.. _examples_Mesh_QuadTree_Create: - -.. ------------------------------ .. -.. THIS FILE IS AUTO GENEREATED .. -.. ------------------------------ .. - -Mesh QuadTree Create -==================== - -.. plot:: - - from SimPEG import Examples - Examples.Mesh_QuadTree_Create.run() - -.. literalinclude:: ../../SimPEG/Examples/Mesh_QuadTree_Create.py - :language: python - :linenos: diff --git a/docs/examples/Mesh_QuadTree_Creation.rst b/docs/examples/Mesh_QuadTree_Creation.rst new file mode 100644 index 00000000..5db5a982 --- /dev/null +++ b/docs/examples/Mesh_QuadTree_Creation.rst @@ -0,0 +1,31 @@ +.. _examples_Mesh_QuadTree_Creation: + +.. --------------------------------- .. +.. .. +.. THIS FILE IS AUTO GENEREATED .. +.. .. +.. SimPEG/Examples/__init__.py .. +.. .. +.. --------------------------------- .. + + +Mesh: QuadTree: Creation +======================== + +You can give the refine method a function, which is evaluated on every cell +of the TreeMesh. + +Occasionally it is useful to initially refine to a constant level +(e.g. 3 in this 32x32 mesh). This means the function is first evaluated +on an 8x8 mesh (2^3). + + + +.. plot:: + + from SimPEG import Examples + Examples.Mesh_QuadTree_Creation.run() + +.. literalinclude:: ../../SimPEG/Examples/Mesh_QuadTree_Creation.py + :language: python + :linenos: diff --git a/docs/examples/Mesh_QuadTree_HangingNodes.rst b/docs/examples/Mesh_QuadTree_HangingNodes.rst new file mode 100644 index 00000000..93875478 --- /dev/null +++ b/docs/examples/Mesh_QuadTree_HangingNodes.rst @@ -0,0 +1,31 @@ +.. _examples_Mesh_QuadTree_HangingNodes: + +.. --------------------------------- .. +.. .. +.. THIS FILE IS AUTO GENEREATED .. +.. .. +.. SimPEG/Examples/__init__.py .. +.. .. +.. --------------------------------- .. + + +Mesh: QuadTree: Hanging Nodes +============================= + +You can give the refine method a function, which is evaluated on every cell +of the TreeMesh. + +Occasionally it is useful to initially refine to a constant level +(e.g. 3 in this 32x32 mesh). This means the function is first evaluated +on an 8x8 mesh (2^3). + + + +.. plot:: + + from SimPEG import Examples + Examples.Mesh_QuadTree_HangingNodes.run() + +.. literalinclude:: ../../SimPEG/Examples/Mesh_QuadTree_HangingNodes.py + :language: python + :linenos: diff --git a/docs/examples/Mesh_Tensor_Creation.rst b/docs/examples/Mesh_Tensor_Creation.rst new file mode 100644 index 00000000..e6cc5b67 --- /dev/null +++ b/docs/examples/Mesh_Tensor_Creation.rst @@ -0,0 +1,43 @@ +.. _examples_Mesh_Tensor_Creation: + +.. --------------------------------- .. +.. .. +.. THIS FILE IS AUTO GENEREATED .. +.. .. +.. SimPEG/Examples/__init__.py .. +.. .. +.. --------------------------------- .. + + + +Mesh: Tensor: Creation +====================== + +For tensor meshes, there are some functions that can come +in handy. For example, creating mesh tensors can be a bit time +consuming, these can be created speedily by just giving numbers +and sizes of padding. See the example below, that follows this +notation:: + + h1 = ( + (cellSize, numPad, [, increaseFactor]), + (cellSize, numCore), + (cellSize, numPad, [, increaseFactor]) + ) + +.. note:: + + You can center your mesh by passing a 'C' for the x0[i] position. + A 'N' will make the entire mesh negative, and a '0' (or a 0) will + make the mesh start at zero. + + + +.. plot:: + + from SimPEG import Examples + Examples.Mesh_Tensor_Creation.run() + +.. literalinclude:: ../../SimPEG/Examples/Mesh_Tensor_Creation.py + :language: python + :linenos: diff --git a/docs/examples/Mesh_ThreeMeshes.rst b/docs/examples/Mesh_ThreeMeshes.rst deleted file mode 100644 index 948398a7..00000000 --- a/docs/examples/Mesh_ThreeMeshes.rst +++ /dev/null @@ -1,17 +0,0 @@ -.. _examples_Mesh_ThreeMeshes: - -.. ------------------------------ .. -.. THIS FILE IS AUTO GENEREATED .. -.. ------------------------------ .. - -Mesh ThreeMeshes -================ - -.. plot:: - - from SimPEG import Examples - Examples.Mesh_ThreeMeshes.run() - -.. literalinclude:: ../../SimPEG/Examples/Mesh_ThreeMeshes.py - :language: python - :linenos: diff --git a/tests/em/examples/__init__.py b/tests/em/examples/__init__.py deleted file mode 100644 index 38d84328..00000000 --- a/tests/em/examples/__init__.py +++ /dev/null @@ -1,11 +0,0 @@ -if __name__ == '__main__': - import os - import glob - import unittest - test_file_strings = glob.glob('test_*.py') - module_strings = [str[0:len(str)-3] for str in test_file_strings] - suites = [unittest.defaultTestLoader.loadTestsFromName(str) for str - in module_strings] - testSuite = unittest.TestSuite(suites) - - unittest.TextTestRunner(verbosity=2).run(testSuite) diff --git a/tests/em/examples/test_Examples.py b/tests/em/examples/test_Examples.py deleted file mode 100644 index 5a601d3b..00000000 --- a/tests/em/examples/test_Examples.py +++ /dev/null @@ -1,10 +0,0 @@ -import unittest, os -from SimPEG.EM import Examples - -class EM_ExamplesRunning(unittest.TestCase): - - def test_CylInversion(self): - Examples.CylInversion.run(plotIt=False) - -if __name__ == '__main__': - unittest.main() diff --git a/tests/flow/test_examples.py b/tests/flow/test_examples.py deleted file mode 100644 index bc519e6d..00000000 --- a/tests/flow/test_examples.py +++ /dev/null @@ -1,12 +0,0 @@ -import unittest -import sys -from SimPEG.FLOW.Examples import Celia1990 -import numpy as np - -class TestCelia1990(unittest.TestCase): - def test_running(self): - Celia1990.run(plotIt=False) - self.assertTrue(True) - -if __name__ == '__main__': - unittest.main() From b2aab4163b655b93c1465012b9e5e2024e308e05 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Thu, 26 Nov 2015 13:56:04 -0700 Subject: [PATCH 19/28] New examples for Mesh --- .../Examples/Mesh_Operators_CahnHilliard.py | 105 ++++++++++++++++++ SimPEG/Examples/Mesh_QuadTree_FaceDiv.py | 49 ++++++++ docs/examples/Mesh_Operators_CahnHilliard.rst | 57 ++++++++++ docs/examples/Mesh_QuadTree_FaceDiv.rst | 26 +++++ 4 files changed, 237 insertions(+) create mode 100644 SimPEG/Examples/Mesh_Operators_CahnHilliard.py create mode 100644 SimPEG/Examples/Mesh_QuadTree_FaceDiv.py create mode 100644 docs/examples/Mesh_Operators_CahnHilliard.rst create mode 100644 docs/examples/Mesh_QuadTree_FaceDiv.rst diff --git a/SimPEG/Examples/Mesh_Operators_CahnHilliard.py b/SimPEG/Examples/Mesh_Operators_CahnHilliard.py new file mode 100644 index 00000000..8e42e617 --- /dev/null +++ b/SimPEG/Examples/Mesh_Operators_CahnHilliard.py @@ -0,0 +1,105 @@ +from SimPEG import * + +def run(plotIt=True, n=60): + """ + Mesh: Operators: Cahn Hilliard + ============================== + + This example is based on the example in the FiPy_ library. + Please see their documentation for more information about the Cahn-Hilliard equation. + + The "Cahn-Hilliard" equation separates a field \\\\( \\\\phi \\\\) into 0 and 1 with smooth transitions. + + .. math:: + + \\frac{\partial \phi}{\partial t} = \\nabla \cdot D \\nabla \left( \\frac{\partial f}{\partial \phi} - \epsilon^2 \\nabla^2 \phi \\right) + + Where \\\\( f \\\\) is the energy function \\\\( f = ( a^2 / 2 )\\\\phi^2(1 - \\\\phi)^2 \\\\) + which drives \\\\( \\\\phi \\\\) towards either 0 or 1, this competes with the term + \\\\(\\\\epsilon^2 \\\\nabla^2 \\\\phi \\\\) which is a diffusion term that creates smooth changes in \\\\( \\\\phi \\\\). + The equation can be factored: + + .. math:: + + \\frac{\partial \phi}{\partial t} = \\nabla \cdot D \\nabla \psi \\\\ + \psi = \\frac{\partial^2 f}{\partial \phi^2} (\phi - \phi^{\\text{old}}) + \\frac{\partial f}{\partial \phi} - \epsilon^2 \\nabla^2 \phi + + Here we will need the derivatives of \\\\( f \\\\): + + .. math:: + + \\frac{\partial f}{\partial \phi} = (a^2/2)2\phi(1-\phi)(1-2\phi) + \\frac{\partial^2 f}{\partial \phi^2} = (a^2/2)2[1-6\phi(1-\phi)] + + The implementation below uses backwards Euler in time with an exponentially increasing time step. + The initial \\\\( \\\\phi \\\\) is a normally distributed field with a standard deviation of 0.1 and mean of 0.5. + The grid is 60x60 and takes a few seconds to solve ~130 times. The results are seen below, and you can see the + field separating as the time increases. + + .. _FiPy: http://www.ctcms.nist.gov/fipy/examples/cahnHilliard/generated/examples.cahnHilliard.mesh2DCoupled.html + + """ + + np.random.seed(5) + + # Here we are going to rearrange the equations: + + # (phi_ - phi)/dt = A*(d2fdphi2*(phi_ - phi) + dfdphi - L*phi_) + # (phi_ - phi)/dt = A*(d2fdphi2*phi_ - d2fdphi2*phi + dfdphi - L*phi_) + # (phi_ - phi)/dt = A*d2fdphi2*phi_ + A*( - d2fdphi2*phi + dfdphi - L*phi_) + # phi_ - phi = dt*A*d2fdphi2*phi_ + dt*A*(- d2fdphi2*phi + dfdphi - L*phi_) + # phi_ - dt*A*d2fdphi2 * phi_ = dt*A*(- d2fdphi2*phi + dfdphi - L*phi_) + phi + # (I - dt*A*d2fdphi2) * phi_ = dt*A*(- d2fdphi2*phi + dfdphi - L*phi_) + phi + # (I - dt*A*d2fdphi2) * phi_ = dt*A*dfdphi - dt*A*d2fdphi2*phi - dt*A*L*phi_ + phi + # (dt*A*d2fdphi2 - I) * phi_ = dt*A*d2fdphi2*phi + dt*A*L*phi_ - phi - dt*A*dfdphi + # (dt*A*d2fdphi2 - I - dt*A*L) * phi_ = (dt*A*d2fdphi2 - I)*phi - dt*A*dfdphi + + h = [(0.25,n)] + M = Mesh.TensorMesh([h,h]) + + # Constants + D = a = epsilon = 1. + I = Utils.speye(M.nC) + + # Operators + A = D * M.faceDiv * M.cellGrad + L = epsilon**2 * M.faceDiv * M.cellGrad + + duration = 75 + elapsed = 0. + dexp = -5 + phi = np.random.normal(loc=0.5,scale=0.01,size=M.nC) + ii, jj = 0, 0 + PHIS = [] + capture = np.logspace(-1,np.log10(duration),8) + while elapsed < duration: + dt = min(100, np.exp(dexp)) + elapsed += dt + dexp += 0.05 + + dfdphi = a**2 * 2 * phi * (1 - phi) * (1 - 2 * phi) + d2fdphi2 = Utils.sdiag(a**2 * 2 * (1 - 6 * phi * (1 - phi))) + + MAT = (dt*A*d2fdphi2 - I - dt*A*L) + rhs = (dt*A*d2fdphi2 - I)*phi - dt*A*dfdphi + phi = Solver(MAT)*rhs + + if elapsed > capture[jj]: + PHIS += [(elapsed, phi.copy())] + jj += 1 + if ii % 10 == 0: print ii, elapsed + ii += 1 + + if plotIt: + import matplotlib.pyplot as plt + fig, axes = plt.subplots(2,4,figsize=(14,6)) + axes = np.array(axes).flatten().tolist() + for ii, ax in zip(np.linspace(0,len(PHIS)-1,len(axes)),axes): + ii = int(ii) + out = M.plotImage(PHIS[ii][1],ax=ax) + ax.axis('off') + ax.set_title('Elapsed Time: %4.1f'%PHIS[ii][0]) + plt.show() + +if __name__ == '__main__': + run() diff --git a/SimPEG/Examples/Mesh_QuadTree_FaceDiv.py b/SimPEG/Examples/Mesh_QuadTree_FaceDiv.py new file mode 100644 index 00000000..5bd67929 --- /dev/null +++ b/SimPEG/Examples/Mesh_QuadTree_FaceDiv.py @@ -0,0 +1,49 @@ +from SimPEG import * + +def run(plotIt=True, n=60): + """ + Mesh: QuadTree: FaceDiv + ======================= + + + + """ + + + M = Mesh.TreeMesh([[(1,16)],[(1,16)]], levels=4) + M._refineCell([0,0,0]) + M._refineCell([0,0,1]) + M._refineCell([4,4,2]) + M.__dirty__ = True + M.number() + + + if plotIt: + import matplotlib.pyplot as plt + fig, axes = plt.subplots(2,1,figsize=(10,10)) + + M.plotGrid(cells=True, nodes=False, ax=axes[0]) + axes[0].axis('off') + axes[0].set_title('Simple QuadTree Mesh') + axes[0].set_xlim([-1,17]) + axes[0].set_ylim([-1,17]) + + for ii, loc in zip(range(M.nC),M.gridCC): + axes[0].text(loc[0]+0.2,loc[1],'%d'%ii, color='r') + + axes[0].plot(M.gridFx[:,0],M.gridFx[:,1], 'g>') + for ii, loc in zip(range(M.nFx),M.gridFx): + axes[0].text(loc[0]+0.2,loc[1],'%d'%ii, color='g') + + axes[0].plot(M.gridFy[:,0],M.gridFy[:,1], 'm^') + for ii, loc in zip(range(M.nFy),M.gridFy): + axes[0].text(loc[0]+0.2,loc[1]+0.2,'%d'%(ii+M.nFx), color='m') + + axes[1].spy(M.faceDiv) + axes[1].set_title('Face Divergence') + axes[1].set_ylabel('Cell Number') + axes[1].set_xlabel('Face Number') + plt.show() + +if __name__ == '__main__': + run() diff --git a/docs/examples/Mesh_Operators_CahnHilliard.rst b/docs/examples/Mesh_Operators_CahnHilliard.rst new file mode 100644 index 00000000..9786e911 --- /dev/null +++ b/docs/examples/Mesh_Operators_CahnHilliard.rst @@ -0,0 +1,57 @@ +.. _examples_Mesh_Operators_CahnHilliard: + +.. --------------------------------- .. +.. .. +.. THIS FILE IS AUTO GENEREATED .. +.. .. +.. SimPEG/Examples/__init__.py .. +.. .. +.. --------------------------------- .. + + +Mesh: Operators: Cahn Hilliard +============================== + +This example is based on the example in the FiPy_ library. +Please see their documentation for more information about the Cahn-Hilliard equation. + +The "Cahn-Hilliard" equation separates a field \\( \\phi \\) into 0 and 1 with smooth transitions. + +.. math:: + + \frac{\partial \phi}{\partial t} = \nabla \cdot D \nabla \left( \frac{\partial f}{\partial \phi} - \epsilon^2 \nabla^2 \phi \right) + +Where \\( f \\) is the energy function \\( f = ( a^2 / 2 )\\phi^2(1 - \\phi)^2 \\) +which drives \\( \\phi \\) towards either 0 or 1, this competes with the term +\\(\\epsilon^2 \\nabla^2 \\phi \\) which is a diffusion term that creates smooth changes in \\( \\phi \\). +The equation can be factored: + +.. math:: + + \frac{\partial \phi}{\partial t} = \nabla \cdot D \nabla \psi \\ + \psi = \frac{\partial^2 f}{\partial \phi^2} (\phi - \phi^{\text{old}}) + \frac{\partial f}{\partial \phi} - \epsilon^2 \nabla^2 \phi + +Here we will need the derivatives of \\( f \\): + +.. math:: + + \frac{\partial f}{\partial \phi} = (a^2/2)2\phi(1-\phi)(1-2\phi) + \frac{\partial^2 f}{\partial \phi^2} = (a^2/2)2[1-6\phi(1-\phi)] + +The implementation below uses backwards Euler in time with an exponentially increasing time step. +The initial \\( \\phi \\) is a normally distributed field with a standard deviation of 0.1 and mean of 0.5. +The grid is 60x60 and takes a few seconds to solve ~130 times. The results are seen below, and you can see the +field separating as the time increases. + +.. _FiPy: http://www.ctcms.nist.gov/fipy/examples/cahnHilliard/generated/examples.cahnHilliard.mesh2DCoupled.html + + + +.. plot:: + + from SimPEG import Examples + Examples.Mesh_Operators_CahnHilliard.run() + +.. literalinclude:: ../../SimPEG/Examples/Mesh_Operators_CahnHilliard.py + :language: python + :linenos: diff --git a/docs/examples/Mesh_QuadTree_FaceDiv.rst b/docs/examples/Mesh_QuadTree_FaceDiv.rst new file mode 100644 index 00000000..6bfdd47f --- /dev/null +++ b/docs/examples/Mesh_QuadTree_FaceDiv.rst @@ -0,0 +1,26 @@ +.. _examples_Mesh_QuadTree_FaceDiv: + +.. --------------------------------- .. +.. .. +.. THIS FILE IS AUTO GENEREATED .. +.. .. +.. SimPEG/Examples/__init__.py .. +.. .. +.. --------------------------------- .. + + +Mesh: QuadTree: FaceDiv +======================= + + + + + +.. plot:: + + from SimPEG import Examples + Examples.Mesh_QuadTree_FaceDiv.run() + +.. literalinclude:: ../../SimPEG/Examples/Mesh_QuadTree_FaceDiv.py + :language: python + :linenos: From f0f3f6c06a0640bc5ded5359b6f5424404e1421d Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Sat, 2 Jan 2016 18:06:49 -0800 Subject: [PATCH 20/28] Update the examples init.py to explicitly reference things. You can run the main function to update the docs and imports. --- SimPEG/Examples/__init__.py | 62 ++++++++++++++++++++++++++++--------- docs/api_Tests.rst | 2 +- 2 files changed, 48 insertions(+), 16 deletions(-) diff --git a/SimPEG/Examples/__init__.py b/SimPEG/Examples/__init__.py index ac219a7f..191664bb 100644 --- a/SimPEG/Examples/__init__.py +++ b/SimPEG/Examples/__init__.py @@ -1,24 +1,60 @@ -# This will import everything in the directory into this file -from os import path as p -from glob import glob -__all__ = [] -for x in glob(p.join(p.dirname(__file__), '*.py')): - if not p.basename(x).startswith('__'): - __import__(p.basename(x)[:-3], globals(), locals()) - __all__ += [p.basename(x)[:-3]] -del glob, p, x +# Run this file to add imports. + +##### AUTOIMPORTS ##### +import EM_FDEM_1D_Inversion +import FLOW_Richards_1D_Celia1990 +import Forward_BasicDirectCurrent +import Inversion_Linear +import Mesh_Basic_PlotImage +import Mesh_Basic_Types +import Mesh_Operators_CahnHilliard +import Mesh_QuadTree_Creation +import Mesh_QuadTree_FaceDiv +import Mesh_QuadTree_HangingNodes +import Mesh_Tensor_Creation +##### AUTOIMPORTS ##### if __name__ == '__main__': """ - Run the following to create the examples documentation. + Run the following to create the examples documentation and add to the imports at the top. """ import shutil, os from SimPEG import Examples + # Create the examples dir in the docs folder. + docExamplesDir = os.path.sep.join(os.path.realpath(__file__).split(os.path.sep)[:-3] + ['docs', 'examples']) + shutil.rmtree(docExamplesDir) + os.makedirs(docExamplesDir) + + # Get all the python examples in this folder + thispath = os.path.sep.join(__file__.split(os.path.sep)[:-1]) + onlyfiles = [f[:-3] for f in os.listdir(thispath) if os.path.isfile(os.path.join(thispath, f)) and f.endswith('.py') and not f.startswith('_')] + + # Add the imports to the top in the AUTOIMPORTS section + f = file(__file__, 'r') + inimports = False + out = '' + for line in f: + if not inimports: + out += line + + if line == "##### AUTOIMPORTS #####\n": + inimports = not inimports + if inimports: + out += '\n'.join(["import %s"%_ for _ in onlyfiles]) + out += '\n##### AUTOIMPORTS #####\n' + f.close() + + f = file(__file__, 'w') + f.write(out) + f.close() + + def _makeExample(filePath, runFunction): + """Makes the example given a path of the file and the run function.""" filePath = os.path.realpath(filePath) name = filePath.split(os.path.sep)[-1].rstrip('.pyc').rstrip('.py') @@ -52,15 +88,11 @@ if __name__ == '__main__': rst = os.path.sep.join((filePath.split(os.path.sep)[:-3] + ['docs', 'examples', name + '.rst'])) + print 'Creating: %s.rst'%name f = open(rst, 'w') f.write(out) f.close() - - docExamplesDir = os.path.sep.join(os.path.realpath(__file__).split(os.path.sep)[:-3] + ['docs', 'examples']) - shutil.rmtree(docExamplesDir) - os.makedirs(docExamplesDir) - for ex in dir(Examples): if ex.startswith('_'): continue E = getattr(Examples,ex) diff --git a/docs/api_Tests.rst b/docs/api_Tests.rst index 614d8747..b0b2fbd0 100644 --- a/docs/api_Tests.rst +++ b/docs/api_Tests.rst @@ -3,6 +3,6 @@ Testing SimPEG ============== -.. automodule:: SimPEG.Tests.TestUtils +.. automodule:: SimPEG.Tests :members: :undoc-members: From a18d48348d83b800115ff0a00cd2e8559a1b128a Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Sat, 2 Jan 2016 18:20:58 -0800 Subject: [PATCH 21/28] Merge fast tests on travis. --- .travis.yml | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/.travis.yml b/.travis.yml index b9826bab..8c8637af 100644 --- a/.travis.yml +++ b/.travis.yml @@ -5,15 +5,13 @@ python: sudo: false env: + - TEST_DIR="tests/mesh tests/base tests/utils" + - TEST_DIR=tests/examples - TEST_DIR=tests/em/fdem/forward - TEST_DIR=tests/em/fdem/inverse/derivs - TEST_DIR=tests/em/fdem/inverse/adjoint - TEST_DIR=tests/em/tdem - - TEST_DIR=tests/mesh - TEST_DIR=tests/flow - - TEST_DIR=tests/utils - - TEST_DIR=tests/base - - TEST_DIR=tests/examples # Setup anaconda before_install: From d64fd4ae3558ae1771aa342ed1a8460184b2cd00 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Sat, 2 Jan 2016 18:29:04 -0800 Subject: [PATCH 22/28] List examples in the init file. --- SimPEG/Examples/__init__.py | 10 +++++++--- tests/examples/test_examples.py | 4 +--- 2 files changed, 8 insertions(+), 6 deletions(-) diff --git a/SimPEG/Examples/__init__.py b/SimPEG/Examples/__init__.py index 191664bb..8431e4ba 100644 --- a/SimPEG/Examples/__init__.py +++ b/SimPEG/Examples/__init__.py @@ -12,6 +12,9 @@ import Mesh_QuadTree_Creation import Mesh_QuadTree_FaceDiv import Mesh_QuadTree_HangingNodes import Mesh_Tensor_Creation + +__examples__ = ["EM_FDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Forward_BasicDirectCurrent", "Inversion_Linear", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation"] + ##### AUTOIMPORTS ##### if __name__ == '__main__': @@ -30,8 +33,8 @@ if __name__ == '__main__': os.makedirs(docExamplesDir) # Get all the python examples in this folder - thispath = os.path.sep.join(__file__.split(os.path.sep)[:-1]) - onlyfiles = [f[:-3] for f in os.listdir(thispath) if os.path.isfile(os.path.join(thispath, f)) and f.endswith('.py') and not f.startswith('_')] + thispath = os.path.sep.join(__file__.split(os.path.sep)[:-1]) + exfiles = [f[:-3] for f in os.listdir(thispath) if os.path.isfile(os.path.join(thispath, f)) and f.endswith('.py') and not f.startswith('_')] # Add the imports to the top in the AUTOIMPORTS section f = file(__file__, 'r') @@ -44,7 +47,8 @@ if __name__ == '__main__': if line == "##### AUTOIMPORTS #####\n": inimports = not inimports if inimports: - out += '\n'.join(["import %s"%_ for _ in onlyfiles]) + out += '\n'.join(["import %s"%_ for _ in exfiles]) + out += '\n\n__examples__ = ["' + '", "'.join(exfiles)+ '"]\n' out += '\n##### AUTOIMPORTS #####\n' f.close() diff --git a/tests/examples/test_examples.py b/tests/examples/test_examples.py index 48cbaa65..cffa2094 100644 --- a/tests/examples/test_examples.py +++ b/tests/examples/test_examples.py @@ -10,10 +10,8 @@ def get(test): self.assertTrue(True) return func attrs = dict() -tests = [_ for _ in dir(Examples) if not _.startswith('_')] -for test in tests: +for test in Examples.__examples__: attrs['test_'+test] = get(test) -del get, tests, _ TestExamples = type('TestExamples', (unittest.TestCase,), attrs) From dd5f5df69eefb8ed106376c73a37b0466fa3adfd Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Sun, 3 Jan 2016 12:48:01 -0800 Subject: [PATCH 23/28] my_function.__name__ must be set to 'test_' for nosetests --- tests/examples/test_examples.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/examples/test_examples.py b/tests/examples/test_examples.py index cffa2094..2e4803b1 100644 --- a/tests/examples/test_examples.py +++ b/tests/examples/test_examples.py @@ -4,11 +4,11 @@ from SimPEG import Examples import numpy as np def get(test): - def func(self): + def test_func(self): print '\nTesting %s.run(plotIt=False)\n'%test getattr(Examples, test).run(plotIt=False) self.assertTrue(True) - return func + return test_func attrs = dict() for test in Examples.__examples__: attrs['test_'+test] = get(test) From 9ebd0fcedc0aa830dcbad483d67d252ed543bf2d Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Sun, 3 Jan 2016 19:02:06 -0800 Subject: [PATCH 24/28] Use pymatsolver in TravisCI --- .travis.yml | 17 +++++++++++++++-- 1 file changed, 15 insertions(+), 2 deletions(-) diff --git a/.travis.yml b/.travis.yml index 8c8637af..db79b031 100644 --- a/.travis.yml +++ b/.travis.yml @@ -4,6 +4,16 @@ python: sudo: false +addons: + apt: + packages: + - gcc + - gfortran + - libopenmpi-dev + - libmumps-seq-dev + - libblas-dev + - liblapack-dev + env: - TEST_DIR="tests/mesh tests/base tests/utils" - TEST_DIR=tests/examples @@ -23,9 +33,12 @@ before_install: # Install packages install: - - conda install --yes pip python=$TRAVIS_PYTHON_VERSION numpy scipy matplotlib cython ipython + - conda install --yes pip python=$TRAVIS_PYTHON_VERSION numpy scipy matplotlib cython ipython nose - pip install nose-cov python-coveralls - # - pip install -r requirements.txt + + - git clone https://github.com/rowanc1/pymatsolver.git + - cd pymatsolver; python setup.py install; cd .. + - python setup.py install - python setup.py build_ext --inplace From 00a29d27aa6f1bb0815342a5e3d4cf862db34d36 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Sun, 3 Jan 2016 19:11:45 -0800 Subject: [PATCH 25/28] Decrease the number of iterations on the EM derivatives. This is to increase the speed at which travis runs tests. --- tests/em/fdem/inverse/derivs/test_FDEM_derivs.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py b/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py index 52108c4e..d3bcb218 100644 --- a/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py +++ b/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py @@ -38,7 +38,7 @@ def derivTest(fdemType, comp): 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=FLR) + return Tests.checkDerivative(fun, x0, num=2, plotIt=False, eps=FLR) class FDEM_DerivTests(unittest.TestCase): From d93a23306ce222eee6c3caab58ba08979e5b4154 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Sun, 10 Jan 2016 12:47:43 -0800 Subject: [PATCH 26/28] =?UTF-8?q?Bump=20version:=200.1.3=20=E2=86=92=200.1?= =?UTF-8?q?.4?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .bumpversion.cfg | 2 +- SimPEG/__init__.py | 2 +- docs/conf.py | 4 ++-- setup.py | 2 +- 4 files changed, 5 insertions(+), 5 deletions(-) diff --git a/.bumpversion.cfg b/.bumpversion.cfg index f739d080..f5efbf0d 100644 --- a/.bumpversion.cfg +++ b/.bumpversion.cfg @@ -1,4 +1,4 @@ [bumpversion] -current_version = 0.1.3 +current_version = 0.1.4 files = setup.py SimPEG/__init__.py docs/conf.py diff --git a/SimPEG/__init__.py b/SimPEG/__init__.py index 369134e7..94bf255f 100644 --- a/SimPEG/__init__.py +++ b/SimPEG/__init__.py @@ -15,7 +15,7 @@ import Directives import Inversion import Tests -__version__ = '0.1.3' +__version__ = '0.1.4' __author__ = 'Rowan Cockett' __license__ = 'MIT' __copyright__ = 'Copyright 2014 Rowan Cockett' diff --git a/docs/conf.py b/docs/conf.py index 3bb33621..4f3f4462 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -51,9 +51,9 @@ copyright = u'2013, SimPEG Developers' # built documents. # # The short X.Y version. -version = '0.1.3' +version = '0.1.4' # The full version, including alpha/beta/rc tags. -release = '0.1.3' +release = '0.1.4' # The language for content autogenerated by Sphinx. Refer to documentation # for a list of supported languages. diff --git a/setup.py b/setup.py index 531dc539..a90d6119 100644 --- a/setup.py +++ b/setup.py @@ -76,7 +76,7 @@ with open("README.rst") as f: setup( name = "SimPEG", - version = "0.1.3", + version = "0.1.4", packages = find_packages(), install_requires = ['numpy>=1.7', 'scipy>=0.13', From 8da717521cf95484c21812a94982c18d8c5356c8 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Sun, 10 Jan 2016 13:01:00 -0800 Subject: [PATCH 27/28] Update scripts for pip --- setup.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/setup.py b/setup.py index a90d6119..34c29703 100644 --- a/setup.py +++ b/setup.py @@ -64,6 +64,7 @@ cython_files = [ "SimPEG/Mesh/TreeUtils" ] extensions = [Extension(f, [f+ext]) for f in cython_files] +scripts = [f+'.pyx' for f in cython_files] if USE_CYTHON and "cleanall" not in args: from Cython.Build import cythonize @@ -95,5 +96,6 @@ setup( use_2to3 = False, include_dirs=[np.get_include()], ext_modules = extensions, + scripts=scripts, **cythonKwargs ) From 985d5b6469f39d56aea9d5d94ccb1e59fc10dee0 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Sun, 10 Jan 2016 13:01:05 -0800 Subject: [PATCH 28/28] =?UTF-8?q?Bump=20version:=200.1.7=20=E2=86=92=200.1?= =?UTF-8?q?.8=20(+7=20squashed=20commits)=20Squashed=20commits:=20[ac5bb36?= =?UTF-8?q?]=20Bump=20version:=200.1.8=20=E2=86=92=200.1.9=20[8acd6b2]=20I?= =?UTF-8?q?mportException=20-->=20ImportError=20[ac410a8]=20matplotlib.pyp?= =?UTF-8?q?lot=20has=20errors=20on=20import,=20put=20these=20in=20the=20fu?= =?UTF-8?q?nctions=20that=20rely=20on=20them=20directly.=20[f128a20]=20Bum?= =?UTF-8?q?p=20version:=200.1.6=20=E2=86=92=200.1.7=20[5866bea]=20Remove?= =?UTF-8?q?=20IPython=20utils.?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit These are out of date, and have problems on Linux (without a proper visual backend). [a519e56] Bump version: 0.1.5 → 0.1.6 [f45aa83] Bump version: 0.1.4 → 0.1.5 --- .bumpversion.cfg | 2 +- SimPEG/EM/Analytics/FDEM.py | 1 - SimPEG/Examples/EM_FDEM_1D_Inversion.py | 5 +- SimPEG/Examples/FLOW_Richards_1D_Celia1990.py | 2 +- SimPEG/Examples/Forward_BasicDirectCurrent.py | 7 +- SimPEG/Mesh/TreeMesh.py | 33 +++++++-- SimPEG/Mesh/View.py | 68 ++----------------- SimPEG/Tests.py | 2 +- SimPEG/Utils/__init__.py | 1 - SimPEG/Utils/ipythonutils.py | 28 -------- SimPEG/__init__.py | 2 +- docs/conf.py | 4 +- setup.py | 2 +- 13 files changed, 50 insertions(+), 107 deletions(-) delete mode 100644 SimPEG/Utils/ipythonutils.py diff --git a/.bumpversion.cfg b/.bumpversion.cfg index f5efbf0d..ea37cced 100644 --- a/.bumpversion.cfg +++ b/.bumpversion.cfg @@ -1,4 +1,4 @@ [bumpversion] -current_version = 0.1.4 +current_version = 0.1.9 files = setup.py SimPEG/__init__.py docs/conf.py diff --git a/SimPEG/EM/Analytics/FDEM.py b/SimPEG/EM/Analytics/FDEM.py index ce7e623d..9e776fdf 100644 --- a/SimPEG/EM/Analytics/FDEM.py +++ b/SimPEG/EM/Analytics/FDEM.py @@ -2,7 +2,6 @@ from __future__ import division import numpy as np from scipy.constants import mu_0, pi from scipy.special import erf -import matplotlib.pyplot as plt from SimPEG import Utils diff --git a/SimPEG/Examples/EM_FDEM_1D_Inversion.py b/SimPEG/Examples/EM_FDEM_1D_Inversion.py index fdf2750e..aba70f4b 100644 --- a/SimPEG/Examples/EM_FDEM_1D_Inversion.py +++ b/SimPEG/Examples/EM_FDEM_1D_Inversion.py @@ -1,7 +1,7 @@ from SimPEG import * import SimPEG.EM as EM from scipy.constants import mu_0 -import matplotlib.pyplot as plt + def run(plotIt=True): """ @@ -31,6 +31,7 @@ def run(plotIt=True): if plotIt: + import matplotlib.pyplot as plt fig, ax = plt.subplots(1,1, figsize = (3, 6)) plt.semilogx(sigma[active], mesh.vectorCCz[active]) ax.set_ylim(-600, 0) @@ -60,6 +61,7 @@ def run(plotIt=True): survey.Wd = 1/(abs(survey.dobs)*std) if plotIt: + import matplotlib.pyplot as plt fig, ax = plt.subplots(1,1, figsize = (10, 6)) ax.loglog(rx.times, dtrue, 'b.-') ax.loglog(rx.times, survey.dobs, 'r.-') @@ -88,6 +90,7 @@ def run(plotIt=True): mopt = inv.run(m0) if plotIt: + import matplotlib.pyplot as plt fig, ax = plt.subplots(1,1, figsize = (3, 6)) plt.semilogx(sigma[active], mesh.vectorCCz[active]) plt.semilogx(np.exp(mopt), mesh.vectorCCz[active]) diff --git a/SimPEG/Examples/FLOW_Richards_1D_Celia1990.py b/SimPEG/Examples/FLOW_Richards_1D_Celia1990.py index 7c663e6b..4c90ea7b 100644 --- a/SimPEG/Examples/FLOW_Richards_1D_Celia1990.py +++ b/SimPEG/Examples/FLOW_Richards_1D_Celia1990.py @@ -1,6 +1,5 @@ from SimPEG import * from SimPEG.FLOW import Richards -import matplotlib.pyplot as plt def run(plotIt=True): """ @@ -61,6 +60,7 @@ def run(plotIt=True): Hs_H120= getFields(120.,'head') if not plotIt:return + import matplotlib.pyplot as plt plt.figure(figsize=(13,5)) plt.subplot(121) plt.plot(40-M.gridCC, Hs_M10[-1],'b-') diff --git a/SimPEG/Examples/Forward_BasicDirectCurrent.py b/SimPEG/Examples/Forward_BasicDirectCurrent.py index 06c9e79e..efbb287c 100644 --- a/SimPEG/Examples/Forward_BasicDirectCurrent.py +++ b/SimPEG/Examples/Forward_BasicDirectCurrent.py @@ -1,7 +1,4 @@ from SimPEG import Mesh, Utils, np, SolverLU -import matplotlib.pyplot as plt -import matplotlib -from matplotlib.mlab import griddata ## 2D DC forward modeling example with Tensor and Curvilinear Meshes @@ -39,6 +36,10 @@ def run(plotIt=True): if not plotIt: return + import matplotlib.pyplot as plt + import matplotlib + from matplotlib.mlab import griddata + #Step4: Making Figure fig, axes = plt.subplots(1,2,figsize=(12*1.2,4*1.2)) label = ["(a)", "(b)"] diff --git a/SimPEG/Mesh/TreeMesh.py b/SimPEG/Mesh/TreeMesh.py index 997f7be9..1d5c8c8c 100644 --- a/SimPEG/Mesh/TreeMesh.py +++ b/SimPEG/Mesh/TreeMesh.py @@ -90,11 +90,6 @@ # from SimPEG import np, sp, Utils, Solver -import matplotlib.pyplot as plt -import matplotlib -from mpl_toolkits.mplot3d import Axes3D -import matplotlib.colors as colors -import matplotlib.cm as cmx try: import TreeUtils @@ -1973,6 +1968,13 @@ class TreeMesh(BaseTensorMesh, InnerProducts): facesX=False, facesY=False, facesZ=False, edgesX=False, edgesY=False, edgesZ=False): + + import matplotlib.pyplot as plt + import matplotlib + from mpl_toolkits.mplot3d import Axes3D + import matplotlib.colors as colors + import matplotlib.cm as cmx + # self.number() axOpts = {'projection':'3d'} if self.dim == 3 else {} @@ -2094,6 +2096,13 @@ class TreeMesh(BaseTensorMesh, InnerProducts): def plotImage(self, I, ax=None, showIt=False, grid=False, clim=None): if self.dim == 3: raise Exception('Use plot slice?') + + import matplotlib.pyplot as plt + import matplotlib + from mpl_toolkits.mplot3d import Axes3D + import matplotlib.colors as colors + import matplotlib.cm as cmx + if ax is None: ax = plt.subplot(111) jet = cm = plt.get_cmap('jet') cNorm = colors.Normalize( @@ -2123,6 +2132,13 @@ class TreeMesh(BaseTensorMesh, InnerProducts): assert vType in ['CC','F','E'] assert self.dim == 3 + + import matplotlib.pyplot as plt + import matplotlib + from mpl_toolkits.mplot3d import Axes3D + import matplotlib.colors as colors + import matplotlib.cm as cmx + szSliceDim = len(getattr(self, 'h'+normal.lower())) #: Size of the sliced dimension if ind is None: ind = int(szSliceDim/2) assert type(ind) in [int, long], 'ind must be an integer' @@ -2269,6 +2285,13 @@ class CellLookUpException(TreeException): if __name__ == '__main__': + + import matplotlib.pyplot as plt + import matplotlib + from mpl_toolkits.mplot3d import Axes3D + import matplotlib.colors as colors + import matplotlib.cm as cmx + def topo(x): return np.sin(x*(2.*np.pi))*0.3 + 0.5 diff --git a/SimPEG/Mesh/View.py b/SimPEG/Mesh/View.py index b078eb66..089d7d9a 100644 --- a/SimPEG/Mesh/View.py +++ b/SimPEG/Mesh/View.py @@ -1,8 +1,11 @@ import numpy as np -import matplotlib.pyplot as plt -import matplotlib -from mpl_toolkits.mplot3d import Axes3D -from SimPEG.Utils import mkvc, animate +from SimPEG.Utils import mkvc +try: + import matplotlib.pyplot as plt + import matplotlib + from mpl_toolkits.mplot3d import Axes3D +except ImportError, e: + print 'Trouble importing matplotlib.' class TensorView(object): @@ -479,63 +482,6 @@ class TensorView(object): ax.grid(True) if showIt: plt.show() - def slicer(mesh, var, imageType='CC', normal='z', index=0, ax=None, clim=None): - assert normal in 'xyz', 'normal must be x, y, or z' - if ax is None: ax = plt.subplot(111) - I = mesh.r(var,'CC','CC','M') - axes = [p for p in 'xyz' if p not in normal.lower()] - if normal is 'x': I = I[index,:,:] - if normal is 'y': I = I[:,index,:] - if normal is 'z': I = I[:,:,index] - if clim is None: clim = [I.min(),I.max()] - p = ax.pcolormesh(getattr(mesh,'vectorN'+axes[0]),getattr(mesh,'vectorN'+axes[1]),I.T,vmin=clim[0],vmax=clim[1]) - ax.axis('tight') - ax.set_xlabel(axes[0]) - ax.set_ylabel(axes[1]) - return p - - def videoSlicer(mesh,var,imageType='CC',normal='z',figsize=(10,8)): - assert mesh.dim > 2, 'This is for 3D meshes only.' - # First set up the figure, the axis, and the plot element we want to animate - fig = plt.figure(figsize=figsize) - ax = plt.axes() - clim = [var.min(),var.max()] - plt.colorbar(mesh.slicer(var, imageType=imageType, normal=normal, index=0, ax=ax, clim=clim)) - tlt = plt.title(normal) - - def animateFrame(i): - mesh.slicer(var, imageType=imageType, normal=normal, index=i, ax=ax, clim=clim) - tlt.set_text(normal.upper()+('-Slice: %d, %4.4f' % (i,getattr(mesh,'vectorCC'+normal)[i]))) - - return animate(fig, animateFrame, frames=mesh.vnC['xyz'.index(normal)]) - - def video(mesh, var, function, figsize=(10, 8), colorbar=True, skip=1): - """ - Call a function for a list of models to create a video. - - :: - - def function(var, ax, clim, tlt, i): - tlt.set_text('%d'%i) - return mesh.plotImage(var, imageType='CC', ax=ax, clim=clim) - - mesh.video([model1, model2, ..., modeln],function) - """ - # First set up the figure, the axis, and the plot element we want to animate - fig = plt.figure(figsize=figsize) - ax = plt.axes() - VAR = np.concatenate(var) - clim = [VAR.min(),VAR.max()] - tlt = plt.title('') - if colorbar: - plt.colorbar(function(var[0],ax,clim,tlt,0)) - - frames = np.arange(0,len(var),skip) - def animateFrame(j): - i = frames[j] - function(var[i],ax,clim,tlt,i) - - return animate(fig, animateFrame, frames=len(frames)) class CylView(object): diff --git a/SimPEG/Tests.py b/SimPEG/Tests.py index 6d414441..804ca6b8 100644 --- a/SimPEG/Tests.py +++ b/SimPEG/Tests.py @@ -1,5 +1,4 @@ import numpy as np -import matplotlib.pyplot as plt from numpy.linalg import norm from SimPEG.Utils import mkvc, sdiag, diagEst from SimPEG import Utils @@ -311,6 +310,7 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None, expectedOrder=2, tole if plotIt: + import matplotlib.pyplot as plt ax = ax or plt.subplot(111) ax.loglog(h, E0, 'b') ax.loglog(h, E1, 'g--') diff --git a/SimPEG/Utils/__init__.py b/SimPEG/Utils/__init__.py index 5280ae79..250146c1 100644 --- a/SimPEG/Utils/__init__.py +++ b/SimPEG/Utils/__init__.py @@ -3,7 +3,6 @@ from codeutils import * from meshutils import exampleLrmGrid, meshTensor, closestPoints, readUBCTensorMesh, writeUBCTensorMesh, writeUBCTensorModel, readVTRFile, writeVTRFile from curvutils import volTetra, faceInfo, indexCube from interputils import interpmat -from ipythonutils import easyAnimate as animate from CounterUtils import * import ModelBuilder import SolverUtils diff --git a/SimPEG/Utils/ipythonutils.py b/SimPEG/Utils/ipythonutils.py deleted file mode 100644 index 7b66e131..00000000 --- a/SimPEG/Utils/ipythonutils.py +++ /dev/null @@ -1,28 +0,0 @@ -from tempfile import NamedTemporaryFile -import matplotlib.pyplot as plt -from matplotlib import animation - -# http://jakevdp.github.io/blog/2013/05/12/embedding-matplotlib-animations/ -# http://www.renevolution.com/how-to-install-ffmpeg-on-mac-os-x/ - -VIDEO_TAG = """""" - -def anim_to_html(anim): - if not hasattr(anim, '_encoded_video'): - with NamedTemporaryFile(suffix='.mp4') as f: - anim.save(f.name, fps=20, extra_args=['-vcodec', 'libx264', '-pix_fmt', 'yuv420p']) - video = open(f.name, "rb").read() - anim._encoded_video = video.encode("base64") - - return VIDEO_TAG.format(anim._encoded_video) - -def display_animation(anim): - plt.close(anim._fig) - return anim_to_html(anim) - -animation.Animation._repr_html_ = display_animation - -easyAnimate = animation.FuncAnimation diff --git a/SimPEG/__init__.py b/SimPEG/__init__.py index 94bf255f..cc51fd1f 100644 --- a/SimPEG/__init__.py +++ b/SimPEG/__init__.py @@ -15,7 +15,7 @@ import Directives import Inversion import Tests -__version__ = '0.1.4' +__version__ = '0.1.9' __author__ = 'Rowan Cockett' __license__ = 'MIT' __copyright__ = 'Copyright 2014 Rowan Cockett' diff --git a/docs/conf.py b/docs/conf.py index 4f3f4462..fee262de 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -51,9 +51,9 @@ copyright = u'2013, SimPEG Developers' # built documents. # # The short X.Y version. -version = '0.1.4' +version = '0.1.9' # The full version, including alpha/beta/rc tags. -release = '0.1.4' +release = '0.1.9' # The language for content autogenerated by Sphinx. Refer to documentation # for a list of supported languages. diff --git a/setup.py b/setup.py index 34c29703..bcb5b8e3 100644 --- a/setup.py +++ b/setup.py @@ -77,7 +77,7 @@ with open("README.rst") as f: setup( name = "SimPEG", - version = "0.1.4", + version = "0.1.9", packages = find_packages(), install_requires = ['numpy>=1.7', 'scipy>=0.13',