mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-15 11:26:09 +08:00
New Examples
- Put all examples in same directory - Make a single test - use __init__.py to create the docs automatically
This commit is contained in:
@@ -5,7 +5,6 @@ python:
|
||||
sudo: false
|
||||
|
||||
env:
|
||||
- TEST_DIR=tests/em/examples
|
||||
- TEST_DIR=tests/em/fdem/forward
|
||||
- TEST_DIR=tests/em/fdem/inverse/derivs
|
||||
- TEST_DIR=tests/em/fdem/inverse/adjoint
|
||||
|
||||
@@ -1 +0,0 @@
|
||||
import CylInversion
|
||||
@@ -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)]
|
||||
@@ -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()
|
||||
@@ -73,5 +73,4 @@ def run(plotIt=True):
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
Utils._makeExample(__file__)
|
||||
run()
|
||||
|
||||
@@ -1,29 +1,37 @@
|
||||
from SimPEG import *
|
||||
|
||||
class LinearSurvey(Survey.BaseSurvey):
|
||||
def projectFields(self, u):
|
||||
return u
|
||||
|
||||
class LinearProblem(Problem.BaseProblem):
|
||||
"""docstring for LinearProblem"""
|
||||
|
||||
surveyPair = LinearSurvey
|
||||
|
||||
def __init__(self, mesh, G, **kwargs):
|
||||
Problem.BaseProblem.__init__(self, mesh, **kwargs)
|
||||
self.G = G
|
||||
|
||||
def fields(self, m, u=None):
|
||||
return self.G.dot(m)
|
||||
|
||||
def Jvec(self, m, v, u=None):
|
||||
return self.G.dot(v)
|
||||
|
||||
def Jtvec(self, m, v, u=None):
|
||||
return self.G.T.dot(v)
|
||||
|
||||
|
||||
def run(N=100, plotIt=True):
|
||||
"""
|
||||
Inversion: Linear Problem
|
||||
=========================
|
||||
|
||||
Here we go over the basics of creating a linear problem and inversion.
|
||||
|
||||
"""
|
||||
|
||||
class LinearSurvey(Survey.BaseSurvey):
|
||||
def projectFields(self, u):
|
||||
return u
|
||||
|
||||
class LinearProblem(Problem.BaseProblem):
|
||||
|
||||
surveyPair = LinearSurvey
|
||||
|
||||
def __init__(self, mesh, G, **kwargs):
|
||||
Problem.BaseProblem.__init__(self, mesh, **kwargs)
|
||||
self.G = G
|
||||
|
||||
def fields(self, m, u=None):
|
||||
return self.G.dot(m)
|
||||
|
||||
def Jvec(self, m, v, u=None):
|
||||
return self.G.dot(v)
|
||||
|
||||
def Jtvec(self, m, v, u=None):
|
||||
return self.G.T.dot(v)
|
||||
|
||||
|
||||
np.random.seed(1)
|
||||
|
||||
mesh = Mesh.TensorMesh([N])
|
||||
@@ -79,5 +87,4 @@ def run(N=100, plotIt=True):
|
||||
return prob, survey, mesh, mrec
|
||||
|
||||
if __name__ == '__main__':
|
||||
Utils._makeExample(__file__)
|
||||
run()
|
||||
|
||||
@@ -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()
|
||||
@@ -1,6 +1,13 @@
|
||||
from SimPEG import *
|
||||
|
||||
def run(plotIt=True):
|
||||
"""
|
||||
Mesh: Basic: Types
|
||||
==================
|
||||
|
||||
Here we show SimPEG used to create three different types of meshes.
|
||||
|
||||
"""
|
||||
sz = [16,16]
|
||||
tM = Mesh.TensorMesh(sz)
|
||||
qM = Mesh.TreeMesh(sz)
|
||||
@@ -20,5 +27,4 @@ def run(plotIt=True):
|
||||
plt.show()
|
||||
|
||||
if __name__ == '__main__':
|
||||
Utils._makeExample(__file__)
|
||||
run()
|
||||
+12
-2
@@ -1,7 +1,18 @@
|
||||
from SimPEG import *
|
||||
|
||||
def run(plotIt=True):
|
||||
from SimPEG import Mesh, np
|
||||
"""
|
||||
Mesh: QuadTree: Creation
|
||||
========================
|
||||
|
||||
You can give the refine method a function, which is evaluated on every cell
|
||||
of the TreeMesh.
|
||||
|
||||
Occasionally it is useful to initially refine to a constant level
|
||||
(e.g. 3 in this 32x32 mesh). This means the function is first evaluated
|
||||
on an 8x8 mesh (2^3).
|
||||
|
||||
"""
|
||||
M = Mesh.TreeMesh([32,32])
|
||||
M.refine(3)
|
||||
def function(cell):
|
||||
@@ -14,5 +25,4 @@ def run(plotIt=True):
|
||||
if plotIt: M.plotGrid(showIt=True)
|
||||
|
||||
if __name__ == '__main__':
|
||||
Utils._makeExample(__file__)
|
||||
run()
|
||||
@@ -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()
|
||||
@@ -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()
|
||||
|
||||
@@ -5,5 +5,63 @@ __all__ = []
|
||||
for x in glob(p.join(p.dirname(__file__), '*.py')):
|
||||
if not p.basename(x).startswith('__'):
|
||||
__import__(p.basename(x)[:-3], globals(), locals())
|
||||
__all__ += [p.basename(x)]
|
||||
__all__ += [p.basename(x)[:-3]]
|
||||
del glob, p, x
|
||||
|
||||
if __name__ == '__main__':
|
||||
"""
|
||||
|
||||
Run the following to create the examples documentation.
|
||||
|
||||
"""
|
||||
|
||||
import shutil, os
|
||||
from SimPEG import Examples
|
||||
|
||||
def _makeExample(filePath, runFunction):
|
||||
filePath = os.path.realpath(filePath)
|
||||
name = filePath.split(os.path.sep)[-1].rstrip('.pyc').rstrip('.py')
|
||||
|
||||
docstr = runFunction.__doc__
|
||||
if docstr is None:
|
||||
doc = '%s\n%s'%(name.replace('_',' '),'='*len(name))
|
||||
else:
|
||||
doc = '\n'.join([_[8:].rstrip() for _ in docstr.split('\n')])
|
||||
|
||||
out = """.. _examples_%s:
|
||||
|
||||
.. --------------------------------- ..
|
||||
.. ..
|
||||
.. THIS FILE IS AUTO GENEREATED ..
|
||||
.. ..
|
||||
.. SimPEG/Examples/__init__.py ..
|
||||
.. ..
|
||||
.. --------------------------------- ..
|
||||
|
||||
%s
|
||||
|
||||
.. plot::
|
||||
|
||||
from SimPEG import Examples
|
||||
Examples.%s.run()
|
||||
|
||||
.. literalinclude:: ../../SimPEG/Examples/%s.py
|
||||
:language: python
|
||||
:linenos:
|
||||
"""%(name,doc,name,name)
|
||||
|
||||
rst = os.path.sep.join((filePath.split(os.path.sep)[:-3] + ['docs', 'examples', name + '.rst']))
|
||||
|
||||
f = open(rst, 'w')
|
||||
f.write(out)
|
||||
f.close()
|
||||
|
||||
|
||||
docExamplesDir = os.path.sep.join(os.path.realpath(__file__).split(os.path.sep)[:-3] + ['docs', 'examples'])
|
||||
shutil.rmtree(docExamplesDir)
|
||||
os.makedirs(docExamplesDir)
|
||||
|
||||
for ex in dir(Examples):
|
||||
if ex.startswith('_'): continue
|
||||
E = getattr(Examples,ex)
|
||||
_makeExample(E.__file__, E.run)
|
||||
|
||||
@@ -1 +0,0 @@
|
||||
import Celia1990
|
||||
+24
-15
@@ -1975,24 +1975,28 @@ class TreeMesh(BaseTensorMesh, InnerProducts):
|
||||
fig = ax.figure
|
||||
|
||||
if grid:
|
||||
X, Y, Z = [], [], []
|
||||
for ind in self._sortedCells:
|
||||
p = self._asPointer(ind)
|
||||
n = self._cellN(p)
|
||||
h = self._cellH(p)
|
||||
x = [n[0] , n[0] + h[0], n[0] + h[0], n[0] , n[0]]
|
||||
y = [n[1] , n[1] , n[1] + h[1], n[1] + h[1], n[1]]
|
||||
if self.dim == 2:
|
||||
ax.plot(x,y, 'b-')
|
||||
X += [n[0] , n[0] + h[0], n[0] + h[0], n[0] , n[0], np.nan]
|
||||
Y += [n[1] , n[1] , n[1] + h[1], n[1] + h[1], n[1], np.nan]
|
||||
elif self.dim == 3:
|
||||
ax.plot(x,y, 'b-', zs=[n[2]]*5)
|
||||
z = [n[2] + h[2], n[2] + h[2], n[2] + h[2], n[2] + h[2], n[2] + h[2]]
|
||||
ax.plot(x,y, 'b-', zs=z)
|
||||
X += [n[0] , n[0] + h[0], n[0] + h[0], n[0] , n[0], np.nan]*2
|
||||
Y += [n[1] , n[1] , n[1] + h[1], n[1] + h[1], n[1], np.nan]*2
|
||||
Z += [n[2]]*5+[np.nan]
|
||||
Z += [n[2] + h[2], n[2] + h[2], n[2] + h[2], n[2] + h[2], n[2] + h[2], np.nan]
|
||||
sides = [0,0], [h[0],0], [0,h[1]], [h[0],h[1]]
|
||||
for s in sides:
|
||||
x = [n[0] + s[0], n[0] + s[0]]
|
||||
y = [n[1] + s[1], n[1] + s[1]]
|
||||
z = [n[2] , n[2] + h[2]]
|
||||
ax.plot(x,y, 'b-', zs=z)
|
||||
X += [n[0] + s[0], n[0] + s[0]]
|
||||
Y += [n[1] + s[1], n[1] + s[1]]
|
||||
Z += [n[2] , n[2] + h[2]]
|
||||
if self.dim == 2:
|
||||
ax.plot(X,Y, 'b-')
|
||||
elif self.dim == 3:
|
||||
ax.plot(X,Y, 'b-', zs=Z)
|
||||
|
||||
if self.dim == 2:
|
||||
if cells:
|
||||
@@ -2004,11 +2008,11 @@ class TreeMesh(BaseTensorMesh, InnerProducts):
|
||||
ax.plot(self._gridN[:,0], self._gridN[:,1], 'ms')
|
||||
ax.plot(self._gridN[self._hangingN.keys(),0], self._gridN[self._hangingN.keys(),1], 'ms', ms=10, mfc='none', mec='m')
|
||||
if facesX:
|
||||
ax.plot(self._gridFx[self._hangingFx.keys(),0], self._gridFx[self._hangingFx.keys(),1], 'gs', ms=10, mfc='none', mec='g')
|
||||
ax.plot(self._gridFx[:,0], self._gridFx[:,1], 'g>')
|
||||
ax.plot(self._gridFx[self._hangingFx.keys(),0], self._gridFx[self._hangingFx.keys(),1], 'gs', ms=10, mfc='none', mec='g')
|
||||
if facesY:
|
||||
ax.plot(self._gridFy[self._hangingFy.keys(),0], self._gridFy[self._hangingFy.keys(),1], 'gs', ms=10, mfc='none', mec='g')
|
||||
ax.plot(self._gridFy[:,0], self._gridFy[:,1], 'g^')
|
||||
ax.plot(self._gridFy[self._hangingFy.keys(),0], self._gridFy[self._hangingFy.keys(),1], 'gs', ms=10, mfc='none', mec='g')
|
||||
ax.set_xlabel('x1')
|
||||
ax.set_ylabel('x2')
|
||||
elif self.dim == 3:
|
||||
@@ -2079,12 +2083,15 @@ class TreeMesh(BaseTensorMesh, InnerProducts):
|
||||
ax.grid(True)
|
||||
if showIt:plt.show()
|
||||
|
||||
def plotImage(self, I, ax=None, showIt=True, grid=False):
|
||||
def plotImage(self, I, ax=None, showIt=False, grid=False, clim=None):
|
||||
if self.dim == 3: raise Exception('Use plot slice?')
|
||||
|
||||
if ax is None: ax = plt.subplot(111)
|
||||
jet = cm = plt.get_cmap('jet')
|
||||
cNorm = colors.Normalize(vmin=I.min(), vmax=I.max())
|
||||
cNorm = colors.Normalize(
|
||||
vmin=I.min() if clim is None else clim[0],
|
||||
vmax=I.max() if clim is None else clim[1])
|
||||
|
||||
scalarMap = cmx.ScalarMappable(norm=cNorm, cmap=jet)
|
||||
ax.set_xlim((self.x0[0], self.h[0].sum()))
|
||||
ax.set_ylim((self.x0[1], self.h[1].sum()))
|
||||
@@ -2093,8 +2100,10 @@ class TreeMesh(BaseTensorMesh, InnerProducts):
|
||||
ax.add_patch(plt.Rectangle((x0[0], x0[1]), sz[0], sz[1], facecolor=scalarMap.to_rgba(I[ii]), edgecolor='k' if grid else 'none'))
|
||||
# if text: ax.text(self.center[0],self.center[1],self.num)
|
||||
scalarMap._A = [] # http://stackoverflow.com/questions/8342549/matplotlib-add-colorbar-to-a-sequence-of-line-plots
|
||||
plt.colorbar(scalarMap)
|
||||
ax.set_xlabel('x')
|
||||
ax.set_ylabel('y')
|
||||
if showIt: plt.show()
|
||||
return [scalarMap]
|
||||
|
||||
def plotSlice(self, v, vType='CC',
|
||||
normal='Z', ind=None, grid=True, view='real',
|
||||
|
||||
@@ -1,6 +1,5 @@
|
||||
from matutils import *
|
||||
from codeutils import *
|
||||
from codeutils import _makeExample
|
||||
from meshutils import exampleLrmGrid, meshTensor, closestPoints, readUBCTensorMesh, writeUBCTensorMesh, writeUBCTensorModel, readVTRFile, writeVTRFile
|
||||
from curvutils import volTetra, faceInfo, indexCube
|
||||
from interputils import interpmat
|
||||
|
||||
@@ -227,35 +227,3 @@ def requires(var):
|
||||
|
||||
return requiresVarWrapper
|
||||
return requiresVar
|
||||
|
||||
def _makeExample(filePath):
|
||||
|
||||
import os
|
||||
name = filePath.split(os.path.sep)[-1][:-3]
|
||||
out = """.. _examples_%s:
|
||||
|
||||
.. ------------------------------ ..
|
||||
.. THIS FILE IS AUTO GENEREATED ..
|
||||
.. ------------------------------ ..
|
||||
|
||||
%s
|
||||
%s
|
||||
|
||||
.. plot::
|
||||
|
||||
from SimPEG import Examples
|
||||
Examples.%s.run()
|
||||
|
||||
.. literalinclude:: ../../SimPEG/Examples/%s.py
|
||||
:language: python
|
||||
:linenos:
|
||||
"""%(name,name.replace('_',' '),'='*len(name),name,name)
|
||||
|
||||
rst = os.path.sep.join((filePath.split(os.path.sep)[:-3] + ['docs', 'examples', name + '.rst']))
|
||||
|
||||
f = open(rst, 'w')
|
||||
f.write(out)
|
||||
f.close()
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -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:
|
||||
@@ -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:
|
||||
@@ -1,8 +1,12 @@
|
||||
.. _examples_Forward_BasicDirectCurrent:
|
||||
|
||||
.. ------------------------------ ..
|
||||
.. THIS FILE IS AUTO GENEREATED ..
|
||||
.. ------------------------------ ..
|
||||
.. --------------------------------- ..
|
||||
.. ..
|
||||
.. THIS FILE IS AUTO GENEREATED ..
|
||||
.. ..
|
||||
.. SimPEG/Examples/__init__.py ..
|
||||
.. ..
|
||||
.. --------------------------------- ..
|
||||
|
||||
Forward BasicDirectCurrent
|
||||
==========================
|
||||
|
||||
@@ -1,11 +1,20 @@
|
||||
.. _examples_Inversion_Linear:
|
||||
|
||||
.. ------------------------------ ..
|
||||
.. THIS FILE IS AUTO GENEREATED ..
|
||||
.. ------------------------------ ..
|
||||
.. --------------------------------- ..
|
||||
.. ..
|
||||
.. THIS FILE IS AUTO GENEREATED ..
|
||||
.. ..
|
||||
.. SimPEG/Examples/__init__.py ..
|
||||
.. ..
|
||||
.. --------------------------------- ..
|
||||
|
||||
|
||||
Inversion: Linear Problem
|
||||
=========================
|
||||
|
||||
Here we go over the basics of creating a linear problem and inversion.
|
||||
|
||||
|
||||
Inversion Linear
|
||||
================
|
||||
|
||||
.. plot::
|
||||
|
||||
|
||||
@@ -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:
|
||||
@@ -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:
|
||||
@@ -1,17 +0,0 @@
|
||||
.. _examples_Mesh_QuadTree_Create:
|
||||
|
||||
.. ------------------------------ ..
|
||||
.. THIS FILE IS AUTO GENEREATED ..
|
||||
.. ------------------------------ ..
|
||||
|
||||
Mesh QuadTree Create
|
||||
====================
|
||||
|
||||
.. plot::
|
||||
|
||||
from SimPEG import Examples
|
||||
Examples.Mesh_QuadTree_Create.run()
|
||||
|
||||
.. literalinclude:: ../../SimPEG/Examples/Mesh_QuadTree_Create.py
|
||||
:language: python
|
||||
:linenos:
|
||||
@@ -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:
|
||||
@@ -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:
|
||||
@@ -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:
|
||||
@@ -1,17 +0,0 @@
|
||||
.. _examples_Mesh_ThreeMeshes:
|
||||
|
||||
.. ------------------------------ ..
|
||||
.. THIS FILE IS AUTO GENEREATED ..
|
||||
.. ------------------------------ ..
|
||||
|
||||
Mesh ThreeMeshes
|
||||
================
|
||||
|
||||
.. plot::
|
||||
|
||||
from SimPEG import Examples
|
||||
Examples.Mesh_ThreeMeshes.run()
|
||||
|
||||
.. literalinclude:: ../../SimPEG/Examples/Mesh_ThreeMeshes.py
|
||||
:language: python
|
||||
:linenos:
|
||||
@@ -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)
|
||||
@@ -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()
|
||||
@@ -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()
|
||||
Reference in New Issue
Block a user