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

Merge pull request #192 from simpeg/feat/examples

Browse Source 
SimPEG examples
This commit is contained in:
Rowan Cockett
2016-01-03 14:02:38 -08:00
parent fe94bd3edf dd5f5df69e
commit 6747882e52
34 changed files with 930 additions and 112 deletions
Download Patch File Download Diff File Expand all files Collapse all files
.travis.yml
+2 -5
View File
@@ -5,16 +5,13 @@ python:
sudo: false
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:
SimPEG/EM/Examples/__init__.py
-1
View File
@@ -1 +0,0 @@
import CylInversion
SimPEG/EM/Examples/CylInversion.py → SimPEG/Examples/EM_FDEM_1D_Inversion.py
+7
View File
@@ -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)]
SimPEG/FLOW/Examples/Celia1990.py → SimPEG/Examples/FLOW_Richards_1D_Celia1990.py
+34
View File
@@ -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()
SimPEG/Examples/DCfwd.py → SimPEG/Examples/Forward_BasicDirectCurrent.py
+2 -1
View File
@@ -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,6 +69,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
+32 -19
View File
@@ -1968,7 +1968,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):
@@ -1983,24 +1983,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:
@@ -2012,11 +2016,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])
@@ -2064,7 +2070,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])
@@ -2080,15 +2085,21 @@ 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?')
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()))
@@ -2097,8 +2108,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',
@@ -2199,7 +2212,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)
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/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:
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/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_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()