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

Merge pull request #194 from simpeg/dev

Browse Source 
Example, Pip, Zero, Identity, Faster Builds
This commit is contained in:
Rowan Cockett
2016-01-10 14:02:01 -08:00
parent 8b53fd5965 985d5b6469
commit 1ba96ba8d0
52 changed files with 1118 additions and 454 deletions
Download Patch File Download Diff File Expand all files Collapse all files
.bumpversion.cfg
+1 -1
View File
@@ -1,4 +1,4 @@
[bumpversion]
current_version = 0.1.3
current_version = 0.1.9
files = setup.py SimPEG/__init__.py docs/conf.py
.travis.yml
+17 -7
View File
@@ -4,17 +4,24 @@ python:
sudo: false
addons:
apt:
packages:
- gcc
- gfortran
- libopenmpi-dev
- libmumps-seq-dev
- libblas-dev
- liblapack-dev
env:
- TEST_DIR=tests/em/examples
- 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:
@@ -26,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
SimPEG/EM/Analytics/FDEM.py
-1
View File
@@ -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
SimPEG/EM/Examples/__init__.py
-1
View File
@@ -1 +0,0 @@
import CylInversion
SimPEG/EM/FDEM/FDEM.py
+48 -139
View File
@@ -66,34 +66,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(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)
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)
@@ -126,30 +120,23 @@ 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(src, dA_duIT, adjoint=True)
if dRHS_dmT is None:
du_dmT = - dA_dmT
else:
du_dmT = -dA_dmT + dRHS_dmT
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
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')
@@ -173,10 +160,8 @@ class BaseFDEMProblem(BaseEMProblem):
for i, src in enumerate(Srcs):
smi, sei = src.eval(self)
if smi is not None:
S_m[:,i] = Utils.mkvc(smi)
if sei is not None:
S_e[:,i] = Utils.mkvc(sei)
S_m[:,i] = S_m[:,i] + smi
S_e[:,i] = S_e[:,i] + sei
return S_m, S_e
@@ -249,39 +234,21 @@ class Problem_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)
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 Problem_b(BaseFDEMProblem):
@@ -362,7 +329,6 @@ class Problem_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 )
@@ -372,51 +338,28 @@ class Problem_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
# 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
MeSigmaIDeriv = self.MeSigmaIDeriv(S_e)
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:
RHSderiv += SrcDeriv
elif SrcDeriv is not None:
RHSderiv = SrcDeriv
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)
SrcDeriv = S_mDeriv(v) + self.MeSigmaI.T * (C.T * S_eDeriv(v))
if RHSderiv is not None:
if self._makeASymmetric is True and not adjoint:
return MfMui.T * RHSderiv
if self._makeASymmetric is True and not adjoint:
return MfMui.T * (SrcDeriv + RHSderiv)
return RHSderiv
return RHSderiv + SrcDeriv
@@ -519,7 +462,7 @@ class Problem_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)
@@ -528,27 +471,10 @@ class Problem_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
@@ -627,35 +553,18 @@ class Problem_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
RHSDeriv = None
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)
MfRhoDeriv = self.MfRhoDeriv(S_e)
if not adjoint:
RHSDeriv = C.T * (MfRhoDeriv * v)
elif adjoint:
RHSDeriv = MfRhoDeriv.T * (C * v)
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)
return RHSDeriv
return RHSDeriv + S_mDeriv(v) + C.T * (MfRho * S_eDeriv(v))
SimPEG/EM/FDEM/FieldsFDEM.py
+28 -47
View File
@@ -3,6 +3,7 @@ import scipy.sparse as sp
import SimPEG
from SimPEG import Utils
from SimPEG.EM.Utils import omega
from SimPEG.Utils import Zero, Identity
class Fields(SimPEG.Problem.Fields):
@@ -32,8 +33,7 @@ class Fields_e(Fields):
ePrimary = np.zeros_like(eSolution)
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, eSolution, srcList):
@@ -43,18 +43,17 @@ class Fields_e(Fields):
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)
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, eSolution, srcList):
@@ -63,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):
@@ -76,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
return 1./(1j * omega(src.freq)) * S_mDeriv
def _b(self, eSolution, srcList):
return self._bPrimary(eSolution, srcList) + self._bSecondary(eSolution, srcList)
@@ -118,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):
@@ -129,26 +124,24 @@ class Fields_b(Fields):
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)
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)
if S_e is not None:
e[:,i] += -self._MeSigmaI * S_e
e[:,i] = e[:,i]+ -self._MeSigmaI * S_e
return e
def _eSecondaryDeriv_u(self, src, v, adjoint=False):
@@ -166,8 +159,7 @@ 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)
w = w - Utils.mkvc(Me * S_e,2)
if not adjoint:
de_dm = self._MeSigmaIDeriv(w) * v
@@ -177,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 += -self._MeSigmaI * Se_Deriv
de_dm = de_dm - self._MeSigmaI * Se_Deriv
return de_dm
@@ -219,8 +210,7 @@ class Fields_j(Fields):
jPrimary = np.zeros_like(jSolution,dtype = complex)
for i, src in enumerate(srcList):
jp = src.jPrimary(self.prob)
if jp is not None:
jPrimary[:,i] += jp
jPrimary[:,i] = jPrimary[:,i] + jp
return jPrimary
def _jSecondary(self, jSolution, srcList):
@@ -230,18 +220,17 @@ class Fields_j(Fields):
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)
for i, src in enumerate(srcList):
hp = src.hPrimary(self.prob)
if hp is not None:
hPrimary[:,i] = hp
hPrimary[:,i] = hPrimary[:,i] + hp
return hPrimary
def _hSecondary(self, jSolution, srcList):
@@ -249,8 +238,7 @@ class Fields_j(Fields):
for i, src in enumerate(srcList):
h[:,i] *= -1./(1j*omega(src.freq))
S_m,_ = src.eval(self.prob)
if S_m is not None:
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):
@@ -276,12 +264,10 @@ class Fields_j(Fields):
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)
hDeriv_m = 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
hDeriv_m = hDeriv_m + 1./(1j*omega(src.freq)) * S_mDeriv
return hDeriv_m
@@ -320,9 +306,8 @@ class Fields_h(Fields):
hPrimary = np.zeros_like(hSolution,dtype = complex)
for i, src in enumerate(srcList):
hp = src.hPrimary(self.prob)
if hp is not None:
hPrimary[:,i] += hp
return hPrimary
hPrimary[:,i] = hPrimary[:,i] + hp
return hPrimary
def _hSecondary(self, hSolution, srcList):
return hSolution
@@ -331,26 +316,24 @@ class Fields_h(Fields):
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)
for i, src in enumerate(srcList):
jp = src.jPrimary(self.prob)
if jp is not None:
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)
if S_e is not None:
j[:,i] += -S_e
j[:,i] = j[:,i]+ -S_e
return j
def _jSecondaryDeriv_u(self, src, v, adjoint=False):
@@ -362,9 +345,7 @@ class Fields_h(Fields):
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
return -S_eDeriv
def _j(self, hSolution, srcList):
return self._jPrimary(hSolution, srcList) + self._jSecondary(hSolution, srcList)
SimPEG/EM/FDEM/SrcFDEM.py
+14 -45
View File
@@ -1,6 +1,7 @@
from SimPEG import Survey, Problem, Utils, np, sp
from scipy.constants import mu_0
from SimPEG.EM.Utils import *
from SimPEG.Utils import Zero
# from SurveyFDEM import Rx
@@ -18,28 +19,28 @@ class BaseSrc(Survey.BaseSrc):
return lambda v: self.S_mDeriv(prob,v,adjoint), lambda v: self.S_eDeriv(prob,v,adjoint)
def bPrimary(self, prob):
return None
return Zero()
def hPrimary(self, prob):
return None
return Zero()
def ePrimary(self, prob):
return None
return Zero()
def jPrimary(self, prob):
return None
return Zero()
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 None
return Zero()
def S_eDeriv(self, prob, v, adjoint = False):
return None
return Zero()
class RawVec_e(BaseSrc):
@@ -51,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):
"""
@@ -85,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):
"""
@@ -192,7 +161,7 @@ class MagDipole(BaseSrc):
def S_e(self, prob):
if all(np.r_[self.mu] == np.r_[prob.curModel.mu]):
return None
return Zero()
else:
eqLocs = prob._eqLocs
@@ -261,7 +230,7 @@ class MagDipole_Bfield(BaseSrc):
def S_e(self, prob):
if all(np.r_[self.mu] == np.r_[prob.curModel.mu]):
return None
return Zero()
else:
eqLocs = prob._eqLocs
@@ -329,7 +298,7 @@ class CircularLoop(BaseSrc):
def S_e(self, prob):
if all(np.r_[self.mu] == np.r_[prob.curModel.mu]):
return None
return Zero()
else:
eqLocs = prob._eqLocs
SimPEG/EM/FDEM/SurveyFDEM.py
+2
View File
@@ -1,8 +1,10 @@
import SimPEG
from SimPEG.EM.Utils import *
from scipy.constants import mu_0
from SimPEG.Utils import Zero, Identity
import SrcFDEM as Src
####################################################
# Receivers
####################################################
SimPEG/EM/Examples/CylInversion.py → SimPEG/Examples/EM_FDEM_1D_Inversion.py
+11 -1
View File
@@ -1,9 +1,16 @@
from SimPEG import *
import SimPEG.EM as EM
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)]
@@ -24,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)
@@ -53,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.-')
@@ -81,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])
SimPEG/FLOW/Examples/Celia1990.py → SimPEG/Examples/FLOW_Richards_1D_Celia1990.py
+35 -1
View File
@@ -1,8 +1,40 @@
from SimPEG import *
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
@@ -28,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-')
@@ -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()
SimPEG/Examples/DCfwd.py → SimPEG/Examples/Forward_BasicDirectCurrent.py
+6 -4
View File
@@ -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
@@ -12,7 +9,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 +35,11 @@ def run(plotIt=True):
phirM = AinvrM*rhsrM
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)"]
@@ -69,6 +70,7 @@ def run(plotIt=True):
else:
axes[i].set_ylabel(" ")
axes[i].set_xlabel("x")
plt.show()
if __name__ == '__main__':
SimPEG/Examples/Linear.py → SimPEG/Examples/Inversion_Linear.py
+36 -24
View File
@@ -1,29 +1,39 @@
from SimPEG import *
class LinearSurvey(Survey.BaseSurvey):
def projectFields(self, u):
return u
class LinearProblem(Problem.BaseProblem):
"""docstring for LinearProblem"""
def run(N=100, plotIt=True):
"""
Inversion: Linear Problem
=========================
surveyPair = LinearSurvey
Here we go over the basics of creating a linear problem and inversion.
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)
class LinearSurvey(Survey.BaseSurvey):
def projectFields(self, u):
return u
def Jvec(self, m, v, u=None):
return self.G.dot(v)
class LinearProblem(Problem.BaseProblem):
def Jtvec(self, m, v, u=None):
return self.G.T.dot(v)
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, plotIt=True):
np.random.seed(1)
mesh = Mesh.TensorMesh([N])
nk = 20
@@ -52,7 +62,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 +73,18 @@ 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)
run()
SimPEG/Examples/Mesh_Basic_PlotImage.py
+46
View File
@@ -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()
SimPEG/Examples/Mesh_Basic_Types.py
+30
View File
@@ -0,0 +1,30 @@
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)
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__':
run()
SimPEG/Examples/Mesh_Operators_CahnHilliard.py
+105
View File
@@ -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()
SimPEG/Examples/Mesh_QuadTree_Creation.py
+28
View File
@@ -0,0 +1,28 @@
from SimPEG import *
def run(plotIt=True):
"""
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):
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__':
run()
SimPEG/Examples/Mesh_QuadTree_FaceDiv.py
+49
View File
@@ -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()
SimPEG/Examples/Mesh_QuadTree_HangingNodes.py
+32
View File
@@ -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()
SimPEG/Examples/Mesh_Tensor_Creation.py
+35
View File
@@ -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()
SimPEG/Examples/__init__.py
+103 -1
View File
@@ -1 +1,103 @@
import Linear, DCfwd
# 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
__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__':
"""
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])
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')
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 exfiles])
out += '\n\n__examples__ = ["' + '", "'.join(exfiles)+ '"]\n'
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')
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']))
print 'Creating: %s.rst'%name
f = open(rst, 'w')
f.write(out)
f.close()
for ex in dir(Examples):
if ex.startswith('_'): continue
E = getattr(Examples,ex)
_makeExample(E.__file__, E.run)
SimPEG/FLOW/Examples/__init__.py
-1
View File
@@ -1 +0,0 @@
import Celia1990
SimPEG/Mesh/TreeMesh.py
+69 -25
View File
@@ -90,13 +90,14 @@
#
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
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 +109,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."
@@ -1960,11 +1963,18 @@ 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):
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 {}
@@ -1975,24 +1985,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 +2018,13 @@ 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:
if cells:
ax.plot(self.gridCC[:,0], self.gridCC[:,1], 'r.', zs=self.gridCC[:,2])
@@ -2056,7 +2072,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,15 +2087,28 @@ 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):
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(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()))
@@ -2089,8 +2117,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',
@@ -2102,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'
@@ -2191,7 +2228,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)
@@ -2248,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
SimPEG/Mesh/View.py
+7 -61
View File
@@ -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):
SimPEG/Tests.py
+1 -1
View File
@@ -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--')
SimPEG/Utils/SolverUtils.py
+8
View File
@@ -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])
SimPEG/Utils/__init__.py
-1
View File
@@ -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
SimPEG/Utils/ipythonutils.py
-28
View File
@@ -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 = """<video controls loop>
<source src="data:video/x-m4v;base64,{0}" type="video/mp4">
Your browser does not support the video tag.
</video>"""
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
SimPEG/Utils/matutils.py
+2
View File
@@ -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
SimPEG/__init__.py
+1 -1
View File
@@ -15,7 +15,7 @@ import Directives
import Inversion
import Tests
__version__ = '0.1.3'
__version__ = '0.1.9'
__author__ = 'Rowan Cockett'
__license__ = 'MIT'
__copyright__ = 'Copyright 2014 Rowan Cockett'
docs/api_Examples.rst
+9 -2
View File
@@ -3,8 +3,15 @@
Examples
********
Forward problem
===============
.. toctree::
:maxdepth: 1
:glob:
examples/*
External Notebooks
==================
* `Example 1: Direct Current <http://www.seogi.me/s/notebooks/DCEx.html>`_
* `Example 2: Seismic-Acoustic <http://www.seogi.me/s/notebooks/SeismicEx.html>`_
docs/api_Mesh.rst
+2 -16
View File
@@ -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
docs/api_Tests.rst
+1 -1
View File
@@ -3,6 +3,6 @@
Testing SimPEG
==============
.. automodule:: SimPEG.Tests.TestUtils
.. automodule:: SimPEG.Tests
:members:
:undoc-members:
docs/conf.py
+2 -2
View File
@@ -51,9 +51,9 @@ copyright = u'2013, SimPEG Developers'
# built documents.
#
# The short X.Y version.
version = '0.1.3'
version = '0.1.9'
# The full version, including alpha/beta/rc tags.
release = '0.1.3'
release = '0.1.9'
# The language for content autogenerated by Sphinx. Refer to documentation
# for a list of supported languages.
docs/examples/EM_FDEM_1D_Inversion.rst
+26
View File
@@ -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:
docs/examples/FLOW_Richards_1D_Celia1990.rst
+52
View File
@@ -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:
docs/examples/Forward_BasicDirectCurrent.rst
+21
View File
@@ -0,0 +1,21 @@
.. _examples_Forward_BasicDirectCurrent:
.. --------------------------------- ..
.. ..
.. THIS FILE IS AUTO GENEREATED ..
.. ..
.. SimPEG/Examples/__init__.py ..
.. ..
.. --------------------------------- ..
Forward BasicDirectCurrent
==========================
.. plot::
from SimPEG import Examples
Examples.Forward_BasicDirectCurrent.run()
.. literalinclude:: ../../SimPEG/Examples/Forward_BasicDirectCurrent.py
:language: python
:linenos:
docs/examples/Inversion_Linear.rst
+26
View File
@@ -0,0 +1,26 @@
.. _examples_Inversion_Linear:
.. --------------------------------- ..
.. ..
.. 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.
.. plot::
from SimPEG import Examples
Examples.Inversion_Linear.run()
.. literalinclude:: ../../SimPEG/Examples/Inversion_Linear.py
:language: python
:linenos:
docs/examples/Mesh_Basic_PlotImage.rst
+27
View File
@@ -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:
docs/examples/Mesh_Basic_Types.rst
+26
View File
@@ -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:
docs/examples/Mesh_Operators_CahnHilliard.rst
+57
View File
@@ -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:
docs/examples/Mesh_QuadTree_Creation.rst
+31
View File
@@ -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:
docs/examples/Mesh_QuadTree_FaceDiv.rst
+26
View File
@@ -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:
docs/examples/Mesh_QuadTree_HangingNodes.rst
+31
View File
@@ -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:
docs/examples/Mesh_Tensor_Creation.rst
+43
View File
@@ -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:
setup.py
+3 -1
View File
@@ -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
@@ -76,7 +77,7 @@ with open("README.rst") as f:
setup(
name = "SimPEG",
version = "0.1.3",
version = "0.1.9",
packages = find_packages(),
install_requires = ['numpy>=1.7',
'scipy>=0.13',
@@ -95,5 +96,6 @@ setup(
use_2to3 = False,
include_dirs=[np.get_include()],
ext_modules = extensions,
scripts=scripts,
**cythonKwargs
)
tests/em/examples/__init__.py
-11
View File
@@ -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)
tests/em/examples/test_Examples.py
-10
View File
@@ -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()
tests/em/fdem/inverse/derivs/test_FDEM_derivs.py
+1 -1
View File
@@ -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):
tests/examples/test_examples.py
+11 -8
View File
@@ -1,17 +1,20 @@
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):
def test_func(self):
print '\nTesting %s.run(plotIt=False)\n'%test
getattr(Examples, test).run(plotIt=False)
self.assertTrue(True)
return test_func
attrs = dict()
for test in Examples.__examples__:
attrs['test_'+test] = get(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()
tests/flow/test_Richards.py
+2
View File
@@ -12,6 +12,8 @@ except Exception, e:
TOL = 1E-8
np.random.seed(0)
class TestModels(unittest.TestCase):
def test_BaseHaverkamp_Theta(self):
tests/flow/test_examples.py
-12
View File
@@ -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()
tests/utils/test_Zero.py
+7
View File
@@ -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):