mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-05 01:42:56 +08:00
organizing the docs - put the content in a content folder. put the SimPEG core api docs in core_api
This commit is contained in:
Lindsey Heagy
@@ -0,0 +1,52 @@
|
||||
.. _api_DataMisfit:
|
||||
|
||||
|
||||
Data Misfit
|
||||
***********
|
||||
|
||||
The data misfit using an l_2 norm is:
|
||||
|
||||
.. math::
|
||||
|
||||
\mu_\text{data} = {1\over 2}\left| \mathbf{W}_d (\mathbf{d}_\text{pred} - \mathbf{d}_\text{obs}) \right|_2^2
|
||||
|
||||
If the field, u, is provided, the calculation of the data is fast:
|
||||
|
||||
.. math::
|
||||
|
||||
\mathbf{d}_\text{pred} = \mathbf{Pu(m)}
|
||||
|
||||
\mathbf{R} = \mathbf{W}_d (\mathbf{d}_\text{pred} - \mathbf{d}_\text{obs})
|
||||
|
||||
Where P is a projection matrix that brings the field on the full domain to the data measurement locations;
|
||||
u is the field of interest; d_obs is the observed data; and \\\(\\mathbf{W}_d\\\) is the weighting matrix.
|
||||
|
||||
The derivative of this, with respect to the model, is:
|
||||
|
||||
.. math::
|
||||
|
||||
\frac{\partial \mu_\text{data}}{\partial \mathbf{m}} = \mathbf{J}^\top \mathbf{W}_d \mathbf{R}
|
||||
|
||||
The second derivative is:
|
||||
|
||||
.. math::
|
||||
|
||||
\frac{\partial^2 \mu_\text{data}}{\partial^2 \mathbf{m}} = \mathbf{J}^\top \mathbf{W}_d \mathbf{W}_d \mathbf{J}
|
||||
|
||||
|
||||
The API
|
||||
=======
|
||||
|
||||
.. autoclass:: SimPEG.DataMisfit.BaseDataMisfit
|
||||
:members:
|
||||
:undoc-members:
|
||||
|
||||
Common Data Misfits
|
||||
===================
|
||||
|
||||
l2 norm
|
||||
-------
|
||||
|
||||
.. autoclass:: SimPEG.DataMisfit.l2_DataMisfit
|
||||
:members:
|
||||
:undoc-members:
|
||||
@@ -0,0 +1,8 @@
|
||||
.. _api_DiffOps:
|
||||
|
||||
Differential Operators
|
||||
======================
|
||||
|
||||
.. automodule:: SimPEG.Mesh.DiffOperators
|
||||
:members:
|
||||
:undoc-members:
|
||||
@@ -0,0 +1,18 @@
|
||||
.. _api_Examples:
|
||||
|
||||
Examples
|
||||
********
|
||||
|
||||
.. 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>`_
|
||||
|
||||
@@ -0,0 +1,27 @@
|
||||
.. _api_FiniteVolume:
|
||||
|
||||
Finite Volume
|
||||
*************
|
||||
|
||||
Any numerical implementation requires the discretization of continuous
|
||||
functions into discrete approximations. These approximations are typically
|
||||
organized in a mesh, which defines boundaries, locations, and connectivity. Of
|
||||
specific interest to geophysical simulations, we require that averaging,
|
||||
interpolation and differential operators be defined for any mesh. In SimPEG,
|
||||
we have implemented a staggered mimetic finite volume approach (`Hyman and
|
||||
Shashkov, 1999 <http://math.lanl.gov/~mac/papers/numerics/HS99B.pdf>`_). This
|
||||
approach requires the definitions of variables at either cell-centers, nodes,
|
||||
faces, or edges as seen in the figure below.
|
||||
|
||||
.. image:: ../../images/finitevolrealestate.png
|
||||
:width: 400 px
|
||||
:alt: FiniteVolume
|
||||
:align: center
|
||||
|
||||
|
||||
.. toctree::
|
||||
:maxdepth: 2
|
||||
|
||||
api_Mesh
|
||||
api_DiffOps
|
||||
api_InnerProducts
|
||||
@@ -0,0 +1,113 @@
|
||||
.. _api_ForwardProblem:
|
||||
|
||||
|
||||
Forward Problem
|
||||
***************
|
||||
|
||||
Problem Class
|
||||
=============
|
||||
|
||||
The problem is a partial differential equation of the form:
|
||||
|
||||
.. math::
|
||||
|
||||
c(m, u) = 0
|
||||
|
||||
Here, \\(m\\) is the model and u is the field (or fields).
|
||||
Given the model, \\(m\\), we can calculate the fields \\(u(m)\\),
|
||||
however, the data we collect is a subset of the fields,
|
||||
and can be defined by a linear projection, \\(P\\).
|
||||
|
||||
.. math::
|
||||
|
||||
d_\text{pred} = P u(m)
|
||||
|
||||
For the inverse problem, we are interested in how changing the model transforms the data,
|
||||
as such we can take write the Taylor expansion:
|
||||
|
||||
.. math::
|
||||
|
||||
Pu(m + hv) = Pu(m) + hP\frac{\partial u(m)}{\partial m} v + \mathcal{O}(h^2 \left\| v \right\| )
|
||||
|
||||
We can linearize and define the sensitivity matrix as:
|
||||
|
||||
.. math::
|
||||
|
||||
J = P\frac{\partial u}{\partial m}
|
||||
|
||||
The sensitivity matrix, and it's transpose will be used in the inverse problem
|
||||
to (locally) find how model parameters change the data, and optimize!
|
||||
|
||||
|
||||
Working with the general PDE, \\(c(m, u) = 0\\), where m is the model and u is the field,
|
||||
the sensitivity is defined as:
|
||||
|
||||
.. math::
|
||||
|
||||
J = P\frac{\partial u}{\partial m}
|
||||
|
||||
We can take the derivative of the PDE:
|
||||
|
||||
.. math::
|
||||
|
||||
\nabla_m c(m, u) \partial m + \nabla_u c(m, u) \partial u = 0
|
||||
|
||||
If the forward problem is invertible, then we can rearrange for
|
||||
\\(\\frac{\\partial u}{\\partial m}\\):
|
||||
|
||||
.. math::
|
||||
|
||||
J = - P \left( \nabla_u c(m, u) \right)^{-1} \nabla_m c(m, u)
|
||||
|
||||
This can often be computed given a vector (i.e. \\(J(v)\\)) rather than
|
||||
stored, as \\(J\\) is a large dense matrix.
|
||||
|
||||
|
||||
|
||||
The API
|
||||
=======
|
||||
|
||||
Problem
|
||||
-------
|
||||
|
||||
.. autoclass:: SimPEG.Problem.BaseProblem
|
||||
:members:
|
||||
:undoc-members:
|
||||
|
||||
.. autoclass:: SimPEG.Problem.BaseTimeProblem
|
||||
:members:
|
||||
:undoc-members:
|
||||
|
||||
Fields
|
||||
------
|
||||
|
||||
.. autoclass:: SimPEG.Fields.Fields
|
||||
:members:
|
||||
:undoc-members:
|
||||
|
||||
.. autoclass:: SimPEG.Fields.TimeFields
|
||||
:members:
|
||||
:undoc-members:
|
||||
|
||||
Survey
|
||||
------
|
||||
|
||||
.. autoclass:: SimPEG.Survey.BaseSurvey
|
||||
:members:
|
||||
:undoc-members:
|
||||
|
||||
.. autoclass:: SimPEG.Survey.BaseSrc
|
||||
:members:
|
||||
:undoc-members:
|
||||
|
||||
.. autoclass:: SimPEG.Survey.BaseRx
|
||||
:members:
|
||||
:undoc-members:
|
||||
|
||||
.. autoclass:: SimPEG.Survey.BaseTimeRx
|
||||
:members:
|
||||
:undoc-members:
|
||||
|
||||
.. autoclass:: SimPEG.Survey.Data
|
||||
:members:
|
||||
:undoc-members:
|
||||
@@ -0,0 +1,301 @@
|
||||
.. _api_InnerProducts:
|
||||
|
||||
|
||||
Inner Products
|
||||
**************
|
||||
|
||||
By using the weak formulation of many of the PDEs in geophysical applications,
|
||||
we can rapidly develop discretizations. Much of this work, however, needs a
|
||||
good understanding of how to approximate inner products on our discretized
|
||||
meshes. We will define the inner product as:
|
||||
|
||||
.. math::
|
||||
|
||||
\left(a,b\right) = \int_\Omega{a \cdot b}{\partial v}
|
||||
|
||||
where a and b are either scalars or vectors.
|
||||
|
||||
.. note::
|
||||
|
||||
The InnerProducts class is a base class providing inner product matrices
|
||||
for meshes and cannot run on its own.
|
||||
|
||||
|
||||
Example problem for DC resistivity
|
||||
----------------------------------
|
||||
|
||||
We will start with the formulation of the Direct Current (DC) resistivity
|
||||
problem in geophysics.
|
||||
|
||||
|
||||
.. math::
|
||||
|
||||
\frac{1}{\sigma}\vec{j} = \nabla \phi \\
|
||||
|
||||
\nabla\cdot \vec{j} = q
|
||||
|
||||
In the following discretization, :math:`\sigma` and :math:`\phi`
|
||||
will be discretized on the cell-centers and the flux, :math:`\vec{j}`,
|
||||
will be on the faces. We will use the weak formulation to discretize
|
||||
the DC resistivity equation.
|
||||
|
||||
We can define in weak form by integrating with a general face function
|
||||
:math:`\vec{f}`:
|
||||
|
||||
.. math::
|
||||
|
||||
\int_{\Omega}{\sigma^{-1}\vec{j} \cdot \vec{f}} = \int_{\Omega}{\nabla \phi \cdot \vec{f}}
|
||||
|
||||
Here we can integrate the right side by parts,
|
||||
|
||||
.. math::
|
||||
|
||||
\nabla\cdot(\phi\vec{f})=\nabla\phi\cdot\vec{f} + \phi\nabla\cdot\vec{f}
|
||||
|
||||
and rearrange it, and apply the Divergence theorem.
|
||||
|
||||
.. math::
|
||||
|
||||
\int_{\Omega}{\sigma^{-1}\vec{j} \cdot \vec{f}} =
|
||||
- \int_{\Omega}{(\phi \nabla \cdot \vec{f})}
|
||||
+ \int_{\partial \Omega}{ \phi \vec{f} \cdot \mathbf{n}}
|
||||
|
||||
We can then discretize for every cell:
|
||||
|
||||
.. math::
|
||||
|
||||
v_{\text{cell}} \sigma^{-1} (\mathbf{J}_x \mathbf{F}_x +\mathbf{J}_y \mathbf{F}_y + \mathbf{J}_z \mathbf{F}_z ) = -\phi^{\top} v_{\text{cell}} \mathbf{D}_{\text{cell}} \mathbf{F} + \text{BC}
|
||||
|
||||
.. note::
|
||||
|
||||
We have discretized the dot product above, but remember that we do not
|
||||
really have a single vector :math:`\mathbf{J}`, but approximations of
|
||||
:math:`\vec{j}` on each face of our cell. In 2D that means 2
|
||||
approximations of :math:`\mathbf{J}_x` and 2 approximations of
|
||||
:math:`\mathbf{J}_y`. In 3D we also have 2 approximations of
|
||||
:math:`\mathbf{J}_z`.
|
||||
|
||||
Regardless of how we choose to approximate this dot product, we can represent
|
||||
this in vector form (again this is for every cell), and will generalize for
|
||||
the case of anisotropic (tensor) sigma.
|
||||
|
||||
.. math::
|
||||
|
||||
\mathbf{F}_c^{\top} (\sqrt{v_{\text{cell}}} \Sigma^{-1} \sqrt{v_{\text{cell}}}) \mathbf{J}_c =
|
||||
-\phi^{\top} v_{\text{cell}} \mathbf{D}_{\text{cell}} \mathbf{F})
|
||||
+ \text{BC}
|
||||
|
||||
We multiply by square-root of volume on each side of the tensor conductivity
|
||||
to keep symmetry in the system. Here :math:`\mathbf{J}_c` is the Cartesian
|
||||
:math:`\mathbf{J}` (on the faces that we choose to use in our approximation)
|
||||
and must be calculated differently depending on the mesh:
|
||||
|
||||
.. math::
|
||||
\mathbf{J}_c = \mathbf{Q}_{(i)}\mathbf{J}_\text{TENSOR} \\
|
||||
\mathbf{J}_c = \mathbf{N}_{(i)}^{-1}\mathbf{Q}_{(i)}\mathbf{J}_\text{Curv}
|
||||
|
||||
Here the :math:`i` index refers to where we choose to approximate this integral, as discussed in the note above.
|
||||
We will approximate this integral by taking the fluxes clustered around every node of the cell, there are 8 combinations in 3D, and 4 in 2D. We will use a projection matrix :math:`\mathbf{Q}_{(i)}` to pick the appropriate fluxes. So, now that we have 8 approximations of this integral, we will just take the average. For the TensorMesh, this looks like:
|
||||
|
||||
.. math::
|
||||
|
||||
\mathbf{F}^{\top}
|
||||
{1\over 8}
|
||||
\left(\sum_{i=1}^8
|
||||
\mathbf{Q}_{(i)}^{\top} \sqrt{v_{\text{cell}}} \Sigma^{-1} \sqrt{v_{\text{cell}}} \mathbf{Q}_{(i)}
|
||||
\right)
|
||||
\mathbf{J}
|
||||
=
|
||||
-\mathbf{F}^{\top} \mathbf{D}_{\text{cell}}^{\top} v_{\text{cell}} \phi + \text{BC}
|
||||
|
||||
Or, when generalizing to the entire mesh and dropping our general face function:
|
||||
|
||||
.. math::
|
||||
|
||||
\mathbf{M}^f_{\Sigma^{-1}} \mathbf{J}
|
||||
=
|
||||
- \mathbf{D}^{\top} \text{diag}(\mathbf{v}) \phi + \text{BC}
|
||||
|
||||
By defining the faceInnerProduct (8 combinations of fluxes in 3D, 4 in 2D, 2 in 1D) to be:
|
||||
|
||||
.. math::
|
||||
|
||||
\mathbf{M}^f_{\Sigma^{-1}} =
|
||||
\sum_{i=1}^{2^d}
|
||||
\mathbf{P}_{(i)}^{\top} \Sigma^{-1} \mathbf{P}_{(i)}
|
||||
|
||||
Where :math:`d` is the dimension of the mesh.
|
||||
The :math:`\mathbf{M}^f` is returned when given the input of :math:`\Sigma^{-1}`.
|
||||
|
||||
Here each :math:`\mathbf{P} ~ \in ~ \mathbb{R}^{(d*nC, nF)}` is a combination
|
||||
of the projection, volume, and any normalization to Cartesian coordinates
|
||||
(where the dot product is well defined):
|
||||
|
||||
.. math::
|
||||
|
||||
\mathbf{P}_{(i)} = \sqrt{ \frac{1}{2^d} \mathbf{I}^d \otimes \text{diag}(\mathbf{v})} \overbrace{\mathbf{N}_{(i)}^{-1}}^{\text{Curv only}} \mathbf{Q}_{(i)}
|
||||
|
||||
.. note::
|
||||
|
||||
This is actually completed for each cell in the mesh at the same time, and the full matrices are returned.
|
||||
|
||||
If ``returnP=True`` is requested in any of these methods the projection matrices are returned as a list ordered by nodes around which the fluxes were approximated::
|
||||
|
||||
# In 3D
|
||||
P = [P000, P100, P010, P110, P001, P101, P011, P111]
|
||||
# In 2D
|
||||
P = [P00, P10, P01, P11]
|
||||
# In 1D
|
||||
P = [P0, P1]
|
||||
|
||||
The derivation for ``edgeInnerProducts`` is exactly the same, however, when we
|
||||
approximate the integral using the fields around each node, the projection
|
||||
matrices look a bit different because we have 12 edges in 3D instead of just 6
|
||||
faces. The interface to the code is exactly the same.
|
||||
|
||||
|
||||
Defining Tensor Properties
|
||||
--------------------------
|
||||
|
||||
**For 3D:**
|
||||
|
||||
Depending on the number of columns (either 1, 3, or 6) of mu, the material
|
||||
property is interpreted as follows:
|
||||
|
||||
.. math::
|
||||
|
||||
\vec{\mu} = \left[\begin{matrix} \mu_{1} & 0 & 0 \\ 0 & \mu_{1} & 0 \\ 0 & 0 & \mu_{1} \end{matrix}\right]
|
||||
|
||||
\vec{\mu} = \left[\begin{matrix} \mu_{1} & 0 & 0 \\ 0 & \mu_{2} & 0 \\ 0 & 0 & \mu_{3} \end{matrix}\right]
|
||||
|
||||
\vec{\mu} = \left[\begin{matrix} \mu_{1} & \mu_{4} & \mu_{5} \\ \mu_{4} & \mu_{2} & \mu_{6} \\ \mu_{5} & \mu_{6} & \mu_{3} \end{matrix}\right]
|
||||
|
||||
**For 2D:**
|
||||
|
||||
Depending on the number of columns (either 1, 2, or 3) of mu, the material property is interpreted as follows:
|
||||
|
||||
.. math::
|
||||
\vec{\mu} = \left[\begin{matrix} \mu_{1} & 0 \\ 0 & \mu_{1} \end{matrix}\right]
|
||||
|
||||
\vec{\mu} = \left[\begin{matrix} \mu_{1} & 0 \\ 0 & \mu_{2} \end{matrix}\right]
|
||||
|
||||
\vec{\mu} = \left[\begin{matrix} \mu_{1} & \mu_{3} \\ \mu_{3} & \mu_{2} \end{matrix}\right]
|
||||
|
||||
|
||||
Structure of Matrices
|
||||
---------------------
|
||||
|
||||
Both the isotropic, and anisotropic material properties result in a diagonal mass matrix.
|
||||
Which is nice and easy to invert if necessary, however, in the fully anisotropic case which is not aligned with the grid, the matrix is not diagonal. This can be seen for a 3D mesh in the figure below.
|
||||
|
||||
.. plot::
|
||||
|
||||
from SimPEG import *
|
||||
mesh = Mesh.TensorMesh([10,50,3])
|
||||
m1 = np.random.rand(mesh.nC)
|
||||
m2 = np.random.rand(mesh.nC,3)
|
||||
m3 = np.random.rand(mesh.nC,6)
|
||||
M = range(3)
|
||||
M[0] = mesh.getFaceInnerProduct(m1)
|
||||
M[1] = mesh.getFaceInnerProduct(m2)
|
||||
M[2] = mesh.getFaceInnerProduct(m3)
|
||||
plt.figure(figsize=(13,5))
|
||||
for i, lab in enumerate(['Isotropic','Anisotropic','Tensor']):
|
||||
plt.subplot(131 + i)
|
||||
plt.spy(M[i],ms=0.5,color='k')
|
||||
plt.tick_params(axis='both',which='both',labeltop='off',labelleft='off')
|
||||
plt.title(lab + ' Material Property')
|
||||
plt.show()
|
||||
|
||||
|
||||
Taking Derivatives
|
||||
------------------
|
||||
|
||||
We will take the derivative of the fully anisotropic tensor for a 3D mesh, the
|
||||
other cases are easier and will not be discussed here. Let us start with one
|
||||
part of the sum which makes up :math:`\mathbf{M}^f_\Sigma` and take the
|
||||
derivative when this is multiplied by some vector :math:`\mathbf{v}`:
|
||||
|
||||
.. math::
|
||||
|
||||
\mathbf{P}^\top \boldsymbol{\Sigma} \mathbf{Pv}
|
||||
|
||||
Here we will let :math:`\mathbf{Pv} = \mathbf{y}` and :math:`\mathbf{y}` will have the form:
|
||||
|
||||
.. math::
|
||||
|
||||
\mathbf{y} = \mathbf{Pv} =
|
||||
\left[
|
||||
\begin{matrix}
|
||||
\mathbf{y}_1\\
|
||||
\mathbf{y}_2\\
|
||||
\mathbf{y}_3\\
|
||||
\end{matrix}
|
||||
\right]
|
||||
|
||||
.. math::
|
||||
|
||||
\mathbf{P}^\top\Sigma\mathbf{y} =
|
||||
\mathbf{P}^\top
|
||||
\left[\begin{matrix}
|
||||
\boldsymbol{\sigma}_{1} & \boldsymbol{\sigma}_{4} & \boldsymbol{\sigma}_{5} \\
|
||||
\boldsymbol{\sigma}_{4} & \boldsymbol{\sigma}_{2} & \boldsymbol{\sigma}_{6} \\
|
||||
\boldsymbol{\sigma}_{5} & \boldsymbol{\sigma}_{6} & \boldsymbol{\sigma}_{3}
|
||||
\end{matrix}\right]
|
||||
\left[
|
||||
\begin{matrix}
|
||||
\mathbf{y}_1\\
|
||||
\mathbf{y}_2\\
|
||||
\mathbf{y}_3\\
|
||||
\end{matrix}
|
||||
\right]
|
||||
=
|
||||
\mathbf{P}^\top
|
||||
\left[
|
||||
\begin{matrix}
|
||||
\boldsymbol{\sigma}_{1}\circ \mathbf{y}_1 + \boldsymbol{\sigma}_{4}\circ \mathbf{y}_2 + \boldsymbol{\sigma}_{5}\circ \mathbf{y}_3\\
|
||||
\boldsymbol{\sigma}_{4}\circ \mathbf{y}_1 + \boldsymbol{\sigma}_{2}\circ \mathbf{y}_2 + \boldsymbol{\sigma}_{6}\circ \mathbf{y}_3\\
|
||||
\boldsymbol{\sigma}_{5}\circ \mathbf{y}_1 + \boldsymbol{\sigma}_{6}\circ \mathbf{y}_2 + \boldsymbol{\sigma}_{3}\circ \mathbf{y}_3\\
|
||||
\end{matrix}
|
||||
\right]
|
||||
|
||||
Now it is easy to take the derivative with respect to any one of the
|
||||
parameters, for example,
|
||||
:math:`\frac{\partial}{\partial\boldsymbol{\sigma}_1}`
|
||||
|
||||
.. math::
|
||||
\frac{\partial}{\partial \boldsymbol{\sigma}_1}\left(\mathbf{P}^\top\Sigma\mathbf{y}\right)
|
||||
=
|
||||
\mathbf{P}^\top
|
||||
\left[
|
||||
\begin{matrix}
|
||||
\text{diag}(\mathbf{y}_1)\\
|
||||
0\\
|
||||
0\\
|
||||
\end{matrix}
|
||||
\right]
|
||||
|
||||
Whereas :math:`\frac{\partial}{\partial\boldsymbol{\sigma}_4}`, for
|
||||
example, is:
|
||||
|
||||
.. math::
|
||||
\frac{\partial}{\partial \boldsymbol{\sigma}_4}\left(\mathbf{P}^\top\Sigma\mathbf{y}\right)
|
||||
=
|
||||
\mathbf{P}^\top
|
||||
\left[
|
||||
\begin{matrix}
|
||||
\text{diag}(\mathbf{y}_2)\\
|
||||
\text{diag}(\mathbf{y}_1)\\
|
||||
0\\
|
||||
\end{matrix}
|
||||
\right]
|
||||
|
||||
These are computed for each of the 8 projections, horizontally concatenated,
|
||||
and returned.
|
||||
|
||||
The API
|
||||
-------
|
||||
|
||||
.. autoclass:: SimPEG.Mesh.InnerProducts.InnerProducts
|
||||
:members:
|
||||
:undoc-members:
|
||||
@@ -0,0 +1,27 @@
|
||||
|
||||
|
||||
InvProblem
|
||||
**********
|
||||
|
||||
.. autoclass:: SimPEG.InvProblem.BaseInvProblem
|
||||
:show-inheritance:
|
||||
:members:
|
||||
:undoc-members:
|
||||
|
||||
|
||||
Inversion
|
||||
*********
|
||||
|
||||
.. autoclass:: SimPEG.Inversion.BaseInversion
|
||||
:show-inheritance:
|
||||
:members:
|
||||
:undoc-members:
|
||||
|
||||
Directives
|
||||
**********
|
||||
|
||||
.. automodule:: SimPEG.Directives
|
||||
:show-inheritance:
|
||||
:members:
|
||||
:undoc-members:
|
||||
|
||||
@@ -0,0 +1,11 @@
|
||||
Inversion Components
|
||||
********************
|
||||
|
||||
.. toctree::
|
||||
:maxdepth: 3
|
||||
|
||||
api_DataMisfit
|
||||
api_Regularization
|
||||
api_Optimization
|
||||
api_Inversion
|
||||
|
||||
@@ -0,0 +1,216 @@
|
||||
.. _api_Maps:
|
||||
|
||||
|
||||
SimPEG Maps
|
||||
***********
|
||||
|
||||
That's not a map...?!
|
||||
=====================
|
||||
|
||||
A SimPEG Map operates on a vector and transforms it to another space.
|
||||
We will use an example commonly applied in electromagnetics (EM) of the
|
||||
log-conductivity model.
|
||||
|
||||
.. math::
|
||||
|
||||
m = \log(\sigma)
|
||||
|
||||
Here we require a *mapping* to get from \\\(m\\\) to \\\(\\sigma\\\),
|
||||
we will call this map \\\(\\mathcal{M}\\\).
|
||||
|
||||
.. math::
|
||||
|
||||
\sigma = \mathcal{M}(m) = \exp(m)
|
||||
|
||||
In SimPEG, we use a (:class:`SimPEG.Maps.ExpMap`) to describe how to map
|
||||
back to conductivity. This is a relatively trivial example (we are just taking
|
||||
the exponential!) but by defining maps we can start to combine and manipulate
|
||||
exactly what we think about as our model, \\\(m\\\). In code, this looks like
|
||||
|
||||
.. code-block:: python
|
||||
:linenos:
|
||||
|
||||
M = Mesh.TensorMesh([100]) # Create a mesh
|
||||
expMap = Maps.ExpMap(M) # Create a mapping
|
||||
m = np.zeros(M.nC) # Create a model vector
|
||||
m[M.vectorCCx>0.5] = 1.0 # Set half of it to 1.0
|
||||
sig = expMap * m # Apply the mapping using *
|
||||
print m
|
||||
# [ 0. 0. 0. 1. 1. 1. ]
|
||||
print sig
|
||||
# [ 1. 1. 1. 2.718 2.718 2.718]
|
||||
|
||||
Combining Maps
|
||||
==============
|
||||
|
||||
We will use an example where we want a 1D layered earth as
|
||||
our model, but we want to map this to a 2D discretization to do our forward
|
||||
modeling. We will also assume that we are working in log conductivity still,
|
||||
so after the transformation we want to map to conductivity space.
|
||||
To do this we will introduce the vertical 1D map (:class:`SimPEG.Maps.SurjectVertical1D`),
|
||||
which does the first part of what we just described. The second part will be
|
||||
done by the :class:`SimPEG.Maps.ExpMap` described above.
|
||||
|
||||
.. code-block:: python
|
||||
:linenos:
|
||||
|
||||
M = Mesh.TensorMesh([7,5])
|
||||
v1dMap = Maps.SurjectVertical1D(M)
|
||||
expMap = Maps.ExpMap(M)
|
||||
myMap = expMap * v1dMap
|
||||
m = np.r_[0.2,1,0.1,2,2.9] # only 5 model parameters!
|
||||
sig = myMap * m
|
||||
|
||||
.. plot::
|
||||
|
||||
from SimPEG import *
|
||||
import matplotlib.pyplot as plt
|
||||
M = Mesh.TensorMesh([7,5])
|
||||
v1dMap = Maps.SurjectVertical1D(M)
|
||||
expMap = Maps.ExpMap(M)
|
||||
myMap = expMap * v1dMap
|
||||
m = np.r_[0.2,1,0.1,2,2.9] # only 5 model parameters!
|
||||
sig = myMap * m
|
||||
figs, axs = plt.subplots(1,2)
|
||||
axs[0].plot(m, M.vectorCCy, 'b-o')
|
||||
axs[0].set_title('Model')
|
||||
axs[0].set_ylabel('Depth, y')
|
||||
axs[0].set_xlabel('Value, $m_i$')
|
||||
axs[0].set_xlim(0,3)
|
||||
axs[0].set_ylim(0,1)
|
||||
clbar = plt.colorbar(M.plotImage(sig,ax=axs[1],grid=True,gridOpts=dict(color='grey'))[0])
|
||||
axs[1].set_title('Physical Property')
|
||||
axs[1].set_ylabel('Depth, y')
|
||||
clbar.set_label('$\sigma = \exp(\mathbf{P}m)$')
|
||||
plt.tight_layout()
|
||||
|
||||
If you noticed, it was pretty easy to combine maps. What is even cooler is
|
||||
that the derivatives also are made for you (if everything goes right).
|
||||
Just to be sure that the derivative is correct, you should always run the test
|
||||
on the mapping that you create.
|
||||
|
||||
Taking Derivatives
|
||||
==================
|
||||
|
||||
Now that we have wrapped up the mapping, we can ensure that it is easy to take
|
||||
derivatives (or at least have access to them!). In the :class:`SimPEG.Maps.ExpMap`
|
||||
there are no dependencies between model parameters, so it will be a diagonal matrix:
|
||||
|
||||
.. math::
|
||||
|
||||
\left(\frac{\partial \mathcal{M}(m)}{\partial m}\right)_{ii} = \frac{\partial \exp(m_i)}{\partial m} = \exp(m_i)
|
||||
|
||||
Or equivalently:
|
||||
|
||||
.. math::
|
||||
|
||||
\frac{\partial \mathcal{M}(m)}{\partial m} = \text{diag}( \exp(m) )
|
||||
|
||||
The mapping API makes this really easy to test that you have got the derivative correct.
|
||||
When these are used in the inverse problem, this is extremely important!!
|
||||
|
||||
.. plot::
|
||||
:include-source:
|
||||
|
||||
from SimPEG import *
|
||||
import matplotlib.pyplot as plt
|
||||
M = Mesh.TensorMesh([100])
|
||||
expMap = Maps.ExpMap(M)
|
||||
m = np.zeros(M.nC)
|
||||
m[M.vectorCCx>0.5] = 1.0
|
||||
expMap.test(m, plotIt=True)
|
||||
|
||||
|
||||
The API
|
||||
=======
|
||||
|
||||
The :code:`IdentityMap` is the base class for all mappings, and it does absolutely nothing.
|
||||
|
||||
.. autoclass:: SimPEG.Maps.IdentityMap
|
||||
:members:
|
||||
:undoc-members:
|
||||
|
||||
|
||||
Common Maps
|
||||
===========
|
||||
|
||||
Exponential Map
|
||||
---------------
|
||||
|
||||
Electrical conductivity varies over many orders of magnitude, so it is a common
|
||||
technique when solving the inverse problem to parameterize and optimize in terms
|
||||
of log conductivity. This makes sense not only because it ensures all conductivities
|
||||
will be positive, but because this is fundamentally the space where conductivity
|
||||
lives (i.e. it varies logarithmically).
|
||||
|
||||
.. autoclass:: SimPEG.Maps.ExpMap
|
||||
:members:
|
||||
:undoc-members:
|
||||
|
||||
|
||||
Vertical 1D Map
|
||||
---------------
|
||||
|
||||
.. autoclass:: SimPEG.Maps.SurjectVertical1D
|
||||
:members:
|
||||
:undoc-members:
|
||||
|
||||
|
||||
Map 2D Cross-Section to 3D Model
|
||||
--------------------------------
|
||||
|
||||
.. autoclass:: SimPEG.Maps.Map2Dto3D
|
||||
:members:
|
||||
:undoc-members:
|
||||
|
||||
|
||||
Mesh to Mesh Map
|
||||
----------------
|
||||
|
||||
|
||||
.. plot::
|
||||
|
||||
from SimPEG import *
|
||||
import matplotlib.pyplot as plt
|
||||
M = Mesh.TensorMesh([100,100])
|
||||
h1 = Utils.meshTensor([(6,7,-1.5),(6,10),(6,7,1.5)])
|
||||
h1 = h1/h1.sum()
|
||||
M2 = Mesh.TensorMesh([h1,h1])
|
||||
V = Utils.ModelBuilder.randomModel(M.vnC, seed=79, its=50)
|
||||
v = Utils.mkvc(V)
|
||||
modh = Maps.Mesh2Mesh([M,M2])
|
||||
modH = Maps.Mesh2Mesh([M2,M])
|
||||
H = modH * v
|
||||
h = modh * H
|
||||
ax = plt.subplot(131)
|
||||
M.plotImage(v, ax=ax)
|
||||
ax.set_title('Fine Mesh (Original)')
|
||||
ax = plt.subplot(132)
|
||||
M2.plotImage(H,clim=[0,1],ax=ax)
|
||||
ax.set_title('Course Mesh')
|
||||
ax = plt.subplot(133)
|
||||
M.plotImage(h,clim=[0,1],ax=ax)
|
||||
ax.set_title('Fine Mesh (Interpolated)')
|
||||
plt.show()
|
||||
|
||||
|
||||
.. autoclass:: SimPEG.Maps.Mesh2Mesh
|
||||
:members:
|
||||
:undoc-members:
|
||||
|
||||
|
||||
Under the Hood
|
||||
==============
|
||||
|
||||
Combo Map
|
||||
---------
|
||||
|
||||
The ComboMap holds the information for multiplying and combining
|
||||
maps. It also uses the chain rule to create the derivative.
|
||||
Remember, any time that you make your own combination of mappings
|
||||
be sure to test that the derivative is correct.
|
||||
|
||||
.. autoclass:: SimPEG.Maps.ComboMap
|
||||
:members:
|
||||
:undoc-members:
|
||||
|
||||
@@ -0,0 +1,193 @@
|
||||
.. _api_Mesh:
|
||||
|
||||
SimPEG Meshes
|
||||
*************
|
||||
|
||||
The Mesh objects in SimPEG provide a numerical grid on which to solve
|
||||
differential equations. Each mesh type has a similar API to make switching
|
||||
between different meshes relatively simple.
|
||||
|
||||
Overview of Meshes Available
|
||||
============================
|
||||
|
||||
The following meshes are available for use:
|
||||
|
||||
.. toctree::
|
||||
:maxdepth: 2
|
||||
|
||||
api_MeshCode
|
||||
|
||||
Each mesh code follows the guiding principles that are present in this
|
||||
tutorial, but the details, advantages and disadvantages differ between
|
||||
the implementations.
|
||||
|
||||
.. plot::
|
||||
|
||||
from SimPEG import Examples
|
||||
Examples.Mesh_Basic_Types.run()
|
||||
|
||||
|
||||
Variable Locations and Terminology
|
||||
==================================
|
||||
|
||||
We will go over the basics of using a TensorMesh, but these skills are transferable
|
||||
to the other meshes available in SimPEG. All of the mesh generation code is located
|
||||
in the Mesh package in SimPEG (i.e. SimPEG.Mesh).
|
||||
|
||||
|
||||
To create a TensorMesh we need to create mesh tensors, the widths of
|
||||
each cell of the mesh in each dimension. We will call these tensors h,
|
||||
and these will be define the constant widths of cells in each dimension
|
||||
of the TensorMesh.
|
||||
|
||||
.. plot::
|
||||
:include-source:
|
||||
|
||||
from SimPEG import Mesh, np
|
||||
import matplotlib.pyplot as plt
|
||||
hx = np.r_[3,2,1,1,1,1,2,3]
|
||||
hy = np.r_[3,1,1,3]
|
||||
M = Mesh.TensorMesh([hx, hy])
|
||||
M.plotGrid(centers=True)
|
||||
plt.show()
|
||||
|
||||
|
||||
In this simple mesh, the hx vector defines the widths of the cell
|
||||
in the x dimension, and starts counting from the origin (0,0). The
|
||||
resulting mesh is divided into cells, and the cell-centers are
|
||||
plotted above as red circles. Other terminology for this mesh are:
|
||||
|
||||
- cell-centers
|
||||
- nodes
|
||||
- faces
|
||||
- edges
|
||||
|
||||
.. plot::
|
||||
:include-source:
|
||||
|
||||
from SimPEG import Mesh, np
|
||||
import matplotlib.pyplot as plt
|
||||
hx = np.r_[3,2,1,1,1,1,2,3]
|
||||
hy = np.r_[3,1,1,3]
|
||||
M = Mesh.TensorMesh([hx, hy])
|
||||
M.plotGrid(faces=True, nodes=True)
|
||||
plt.title('Cell faces in the x- and y-directions.')
|
||||
plt.legend(('Nodes', 'X-Faces', 'Y-Faces'))
|
||||
plt.show()
|
||||
|
||||
Generally, the faces are used to discretize fluxes, quantities that
|
||||
leave or enter the cells. As such, these fluxes have a direction to
|
||||
them, which is normal to the cell (i.e. directly out of the cell face).
|
||||
The plot above shows that x-faces point in the x-direction, and
|
||||
y-faces point in the y-direction. The nodes are shown in blue,
|
||||
and lie at the intersection of the grid lines. In a two-dimensional
|
||||
mesh, the edges actually live in the same location as the faces,
|
||||
however, they align (or are tangent to) the face. This is easier to
|
||||
see in 3D, when the edges do not live in the same location as the faces.
|
||||
In the 3D plot below, the edge variables are seen as black triangles,
|
||||
and live on the edges(!) of the cell.
|
||||
|
||||
.. plot::
|
||||
:include-source:
|
||||
|
||||
from SimPEG import Mesh
|
||||
Mesh.TensorMesh([1,1,1]).plotGrid(faces=True, edges=True, centers=True, showIt=True)
|
||||
|
||||
How many of each?
|
||||
-----------------
|
||||
|
||||
When making variables that live in each of these locations, it is
|
||||
important to know how many of each variable type you are dealing with.
|
||||
SimPEG makes this pretty easy:
|
||||
|
||||
::
|
||||
|
||||
In [1]: print M
|
||||
---- 2-D TensorMesh ----
|
||||
x0: 0.00
|
||||
y0: 0.00
|
||||
nCx: 8
|
||||
nCy: 4
|
||||
hx: 3.00, 2.00, 4*1.00, 2.00, 3.00
|
||||
hy: 3.00, 2*1.00, 3.00
|
||||
|
||||
In [2]: count = {'numCells': M.nC,
|
||||
....: 'numCells_xDir': M.nCx,
|
||||
....: 'numCells_yDir': M.nCy,
|
||||
....: 'numCells_vector': M.vnC}
|
||||
|
||||
In [3]: print 'This mesh has %(numCells)d cells, which is %(numCells_xDir)d*%(numCells_yDir)d!!' % count
|
||||
|
||||
This mesh has 32 cells, which is 8*4!!
|
||||
|
||||
In [4]: print count
|
||||
|
||||
{
|
||||
'numCells_vector': array([8, 4]),
|
||||
'numCells_yDir': 4,
|
||||
'numCells_xDir': 8,
|
||||
'numCells': 32
|
||||
}
|
||||
|
||||
SimPEG also counts the nodes, faces, and edges.
|
||||
|
||||
::
|
||||
|
||||
Nodes: M.nN, M.nNx, M.nNy, M.nNz, M.vnN
|
||||
Faces: M.nF, M.nFx, M.nFy, M.nFz, M.vnF, M.vnFx, M.vnFy, M.vnFz
|
||||
Edges: M.nE, M.nEx, M.nEy, M.nEz, M.vnE, M.vnEx, M.vnEy, M.vnEz
|
||||
|
||||
Face and edge variables have different counts depending on
|
||||
the dimension of the direction that you are interested in.
|
||||
In a 4x5 mesh, for example, there is a 5x5 grid of x-faces,
|
||||
and a 4x6 grid of y-faces. You can count them below!
|
||||
As such, the vnF(x,y,z) and vnE(x,y,z) properties give the
|
||||
vector grid size.
|
||||
|
||||
.. plot::
|
||||
:include-source:
|
||||
|
||||
from SimPEG import Mesh
|
||||
Mesh.TensorMesh([4,5]).plotGrid(faces=True, showIt=True)
|
||||
|
||||
|
||||
Making Tensors
|
||||
--------------
|
||||
|
||||
For tensor meshes, there are some additional 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])
|
||||
)
|
||||
|
||||
.. plot::
|
||||
:include-source:
|
||||
|
||||
from SimPEG import Mesh, Utils
|
||||
h1 = [(10, 5, -1.3), (5, 20), (10, 3, 1.3)]
|
||||
M = Mesh.TensorMesh([h1, h1], x0='CN')
|
||||
M.plotGrid(showIt=True)
|
||||
|
||||
.. 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.
|
||||
|
||||
Hopefully, you now know how to create TensorMesh objects in SimPEG,
|
||||
and by extension you are also familiar with how to create and use
|
||||
other types of meshes in this SimPEG framework.
|
||||
|
||||
|
||||
The API
|
||||
=======
|
||||
|
||||
.. autoclass:: SimPEG.Mesh.BaseMesh.BaseMesh
|
||||
:members:
|
||||
:undoc-members:
|
||||
@@ -0,0 +1,68 @@
|
||||
.. _api_MeshCode:
|
||||
|
||||
Tensor Mesh
|
||||
===========
|
||||
|
||||
.. autoclass:: SimPEG.Mesh.TensorMesh
|
||||
:members:
|
||||
:undoc-members:
|
||||
:show-inheritance:
|
||||
|
||||
Cylindrical Mesh
|
||||
================
|
||||
|
||||
.. autoclass:: SimPEG.Mesh.CylMesh
|
||||
:members:
|
||||
:undoc-members:
|
||||
:show-inheritance:
|
||||
|
||||
Tree Mesh
|
||||
=========
|
||||
|
||||
.. autoclass:: SimPEG.Mesh.TreeMesh
|
||||
:members:
|
||||
:undoc-members:
|
||||
:show-inheritance:
|
||||
|
||||
Curvilinear Mesh
|
||||
================
|
||||
|
||||
.. autoclass:: SimPEG.Mesh.CurvilinearMesh
|
||||
:members:
|
||||
:undoc-members:
|
||||
:show-inheritance:
|
||||
|
||||
|
||||
Base Rectangular Mesh
|
||||
=====================
|
||||
|
||||
.. autoclass:: SimPEG.Mesh.BaseMesh.BaseRectangularMesh
|
||||
:members:
|
||||
:undoc-members:
|
||||
:show-inheritance:
|
||||
|
||||
Base Tensor Mesh
|
||||
================
|
||||
|
||||
.. autoclass:: SimPEG.Mesh.TensorMesh.BaseTensorMesh
|
||||
:members:
|
||||
:undoc-members:
|
||||
:show-inheritance:
|
||||
|
||||
|
||||
Mesh IO
|
||||
=======
|
||||
|
||||
.. automodule:: SimPEG.Mesh.MeshIO
|
||||
:members:
|
||||
:undoc-members:
|
||||
:show-inheritance:
|
||||
|
||||
|
||||
Mesh Viewing
|
||||
============
|
||||
|
||||
.. automodule:: SimPEG.Mesh.View
|
||||
:members:
|
||||
:undoc-members:
|
||||
:show-inheritance:
|
||||
@@ -0,0 +1,9 @@
|
||||
|
||||
Optimize
|
||||
********
|
||||
|
||||
.. automodule:: SimPEG.Optimization
|
||||
:show-inheritance:
|
||||
:private-members:
|
||||
:members:
|
||||
:undoc-members:
|
||||
@@ -0,0 +1,29 @@
|
||||
SimPEG PropMaps
|
||||
***************
|
||||
|
||||
The API
|
||||
=======
|
||||
|
||||
Property
|
||||
--------
|
||||
|
||||
.. autoclass:: SimPEG.PropMaps.Property
|
||||
:members:
|
||||
:undoc-members:
|
||||
|
||||
|
||||
PropMap
|
||||
-------
|
||||
|
||||
.. autoclass:: SimPEG.PropMaps.PropMap
|
||||
:members:
|
||||
:undoc-members:
|
||||
|
||||
|
||||
PropModel
|
||||
---------
|
||||
|
||||
.. autoclass:: SimPEG.PropMaps.PropModel
|
||||
:members:
|
||||
:undoc-members:
|
||||
|
||||
@@ -0,0 +1,100 @@
|
||||
|
||||
Regularization
|
||||
**************
|
||||
|
||||
If there is one model that has a misfit that equals the desired tolerance, then there are infinitely many other models which can fit to the same degree. The challenge is to find that model which has the desired characteristics and is compatible with a priori information. A single model can be selected from an infinite ensemble by measuring the length, or norm, of each model. Then a smallest, or sometimes largest, member can be isolated. Our goal is to design a norm that embodies our prior knowledge and, when minimized, yields a realistic candidate for the solution of our problem. The norm can penalize variation from a reference model, spatial derivatives of the model, or some combination of these.
|
||||
|
||||
Tikhonov Regularization
|
||||
=======================
|
||||
|
||||
Here we will define regularization of a model, m, in general however, this should be thought of as (m-m_ref) but otherwise it is exactly the same:
|
||||
|
||||
.. math::
|
||||
|
||||
R(m) = \int_\Omega \frac{\alpha_x}{2}\left(\frac{\partial m}{\partial x}\right)^2 + \frac{\alpha_y}{2}\left(\frac{\partial m}{\partial y}\right)^2 \partial v
|
||||
|
||||
Our discrete gradient operator works on cell centers and gives the derivative on the cell faces, which is not where we want to be evaluating this integral. We need to average the values back to the cell-centers before we integrate. To avoid null spaces, we square first and then average. In 2D with ij notation it looks like this:
|
||||
|
||||
.. math::
|
||||
|
||||
R(m) \approx \sum_{ij} \left[\frac{\alpha_x}{2}\left[\left(\frac{m_{i+1,j} - m_{i,j}}{h}\right)^2 + \left(\frac{m_{i,j} - m_{i-1,j}}{h}\right)^2\right] \\
|
||||
+ \frac{\alpha_y}{2}\left[\left(\frac{m_{i,j+1} - m_{i,j}}{h}\right)^2 + \left(\frac{m_{i,j} - m_{i,j-1}}{h}\right)^2\right]
|
||||
\right]h^2
|
||||
|
||||
If we let D_1 be the derivative matrix in the x direction
|
||||
|
||||
.. math::
|
||||
|
||||
\mathbf{D}_1 = \mathbf{I}_2\otimes\mathbf{d}_1
|
||||
|
||||
.. math::
|
||||
|
||||
\mathbf{D}_2 = \mathbf{d}_2\otimes\mathbf{I}_1
|
||||
|
||||
Where d_1 is the one dimensional derivative:
|
||||
|
||||
.. math::
|
||||
|
||||
\mathbf{d}_1 = \frac{1}{h} \left[ \begin{array}{cccc}
|
||||
-1 & 1 & & \\
|
||||
& \ddots & \ddots&\\
|
||||
& & -1 & 1\end{array} \right]
|
||||
|
||||
.. math::
|
||||
|
||||
R(m) \approx \mathbf{v}^\top \left[\frac{\alpha_x}{2}\mathbf{A}_1 (\mathbf{D}_1 m) \odot (\mathbf{D}_1 m) + \frac{\alpha_y}{2}\mathbf{A}_2 (\mathbf{D}_2 m) \odot (\mathbf{D}_2 m) \right]
|
||||
|
||||
Recall that this is really a just point wise multiplication, or a diagonal matrix times a vector. When we multiply by something in a diagonal we can interchange and it gives the same results (i.e. it is point wise)
|
||||
|
||||
.. math::
|
||||
|
||||
\mathbf{a\odot b} = \text{diag}(\mathbf{a})\mathbf{b} = \text{diag}(\mathbf{b})\mathbf{a} = \mathbf{b\odot a}
|
||||
|
||||
and the transpose also is true (but the sizes have to make sense...):
|
||||
|
||||
.. math::
|
||||
|
||||
\mathbf{a}^\top\text{diag}(\mathbf{b}) = \mathbf{b}^\top\text{diag}(\mathbf{a})
|
||||
|
||||
So R(m) can simplify to:
|
||||
|
||||
.. math::
|
||||
|
||||
R(m) \approx \mathbf{m}^\top \left[\frac{\alpha_x}{2}\mathbf{D}_1^\top \text{diag}(\mathbf{A}_1^\top\mathbf{v}) \mathbf{D}_1 + \frac{\alpha_y}{2}\mathbf{D}_2^\top \text{diag}(\mathbf{A}_2^\top \mathbf{v}) \mathbf{D}_2 \right] \mathbf{m}
|
||||
|
||||
We will define W_x as:
|
||||
|
||||
.. math::
|
||||
|
||||
\mathbf{W}_x = \sqrt{\alpha_x}\text{diag}\left(\sqrt{\mathbf{A}_1^\top\mathbf{v}}\right) \mathbf{D}_1
|
||||
|
||||
|
||||
And then W as a tall matrix of all of the different regularization terms:
|
||||
|
||||
.. math::
|
||||
|
||||
\mathbf{W} = \left[ \begin{array}{c}
|
||||
\mathbf{W}_s\\
|
||||
\mathbf{W}_x\\
|
||||
\mathbf{W}_y\end{array} \right]
|
||||
|
||||
Then we can write
|
||||
|
||||
.. math::
|
||||
|
||||
R(m) \approx \frac{1}{2}\mathbf{m^\top W^\top W m}
|
||||
|
||||
The API
|
||||
-------
|
||||
|
||||
.. autoclass:: SimPEG.Regularization.BaseRegularization
|
||||
:members:
|
||||
:undoc-members:
|
||||
|
||||
|
||||
.. autoclass:: SimPEG.Regularization.Tikhonov
|
||||
:show-inheritance:
|
||||
:members:
|
||||
|
||||
|
||||
|
||||
@@ -0,0 +1,53 @@
|
||||
.. _api_Solver:
|
||||
|
||||
|
||||
Solver
|
||||
******
|
||||
|
||||
BYOS
|
||||
====
|
||||
|
||||
The numerical linear algebra solver that you use will ultimately be the
|
||||
bottleneck of your large scale inversion. To be the most flexible, SimPEG
|
||||
provides wrappers rather than a comprehensive set of solvers (i.e. BYOS).
|
||||
|
||||
The interface is as follows::
|
||||
|
||||
A # Where A is a sparse matrix (or linear operator)
|
||||
Ainv = Solver(A, **solverOpts) # Create a solver object with key word arguments
|
||||
x = Ainv * b # Where b is a numpy array of shape (n,) or (n,*)
|
||||
Ainv.clean() # This cleans the memory footprint (if any)
|
||||
|
||||
.. note::
|
||||
|
||||
This is somewhat an abuse of notation for solvers as we never actually
|
||||
create A inverse. Instead we are creating an object that acts like A
|
||||
inverse, whether that be a Krylov subspace solver or an LU decomposition.
|
||||
|
||||
To wrap up solvers in scipy.sparse.linalg it takes one line of code::
|
||||
|
||||
|
||||
Solver = SolverWrapD(sp.linalg.spsolve, factorize=False)
|
||||
SolverLU = SolverWrapD(sp.linalg.splu, factorize=True)
|
||||
SolverCG = SolverWrapI(sp.linalg.cg)
|
||||
|
||||
.. note::
|
||||
|
||||
The above solvers are loaded into the base name space of SimPEG.
|
||||
|
||||
.. seealso::
|
||||
|
||||
- https://bitbucket.org/petsc/petsc4py
|
||||
- https://github.com/bfroehle/pymumps
|
||||
- https://github.com/rowanc1/pymatsolver
|
||||
|
||||
|
||||
The API
|
||||
=======
|
||||
|
||||
.. autofunction:: SimPEG.Utils.SolverUtils.SolverWrapD
|
||||
:noindex:
|
||||
|
||||
.. autofunction:: SimPEG.Utils.SolverUtils.SolverWrapI
|
||||
:noindex:
|
||||
|
||||
@@ -0,0 +1,8 @@
|
||||
.. _api_Tests:
|
||||
|
||||
Testing SimPEG
|
||||
==============
|
||||
|
||||
.. automodule:: SimPEG.Tests
|
||||
:members:
|
||||
:undoc-members:
|
||||
@@ -0,0 +1,11 @@
|
||||
Utilities
|
||||
*********
|
||||
|
||||
.. toctree::
|
||||
:maxdepth: 2
|
||||
|
||||
api_Solver
|
||||
api_Maps
|
||||
api_PropMaps
|
||||
api_Utils
|
||||
api_Tests
|
||||
@@ -0,0 +1,92 @@
|
||||
Utils
|
||||
*****
|
||||
|
||||
.. automodule:: SimPEG.Utils
|
||||
:members:
|
||||
:undoc-members:
|
||||
|
||||
Matrix Utilities
|
||||
================
|
||||
|
||||
.. automodule:: SimPEG.Utils.matutils
|
||||
:members:
|
||||
:undoc-members:
|
||||
|
||||
|
||||
Solver Utilities
|
||||
================
|
||||
|
||||
.. automodule:: SimPEG.Utils.SolverUtils
|
||||
:members:
|
||||
:undoc-members:
|
||||
|
||||
Curv Utilities
|
||||
==============
|
||||
|
||||
.. automodule:: SimPEG.Utils.curvutils
|
||||
:members:
|
||||
:undoc-members:
|
||||
|
||||
Mesh Utilities
|
||||
==============
|
||||
|
||||
.. automodule:: SimPEG.Utils.meshutils
|
||||
:members:
|
||||
:undoc-members:
|
||||
|
||||
Model Builder Utilities
|
||||
=======================
|
||||
|
||||
.. automodule:: SimPEG.Utils.ModelBuilder
|
||||
:members:
|
||||
:undoc-members:
|
||||
|
||||
Interpolation Utilities
|
||||
=======================
|
||||
|
||||
.. automodule:: SimPEG.Utils.interputils
|
||||
:members:
|
||||
:undoc-members:
|
||||
|
||||
Counter Utilities
|
||||
=================
|
||||
|
||||
.. code-block:: python
|
||||
:linenos:
|
||||
|
||||
class MyClass(object):
|
||||
def __init__(self, url):
|
||||
self.counter = Counter()
|
||||
|
||||
@count
|
||||
def MyMethod(self):
|
||||
pass
|
||||
|
||||
@timeIt
|
||||
def MySecondMethod(self):
|
||||
pass
|
||||
|
||||
c = MyClass('blah')
|
||||
for i in range(100): c.MyMethod()
|
||||
for i in range(300): c.MySecondMethod()
|
||||
c.counter.summary()
|
||||
|
||||
|
||||
.. code-block:: text
|
||||
:linenos:
|
||||
|
||||
Counters:
|
||||
MyClass.MyMethod : 100
|
||||
|
||||
Times: mean sum
|
||||
MyClass.MySecondMethod : 1.70e-06, 5.10e-04, 300x
|
||||
|
||||
|
||||
|
||||
The API
|
||||
-------
|
||||
|
||||
.. automodule:: SimPEG.Utils.CounterUtils
|
||||
:members:
|
||||
:undoc-members:
|
||||
|
||||
@@ -0,0 +1,69 @@
|
||||
Why SimPEG?
|
||||
===========
|
||||
|
||||
Our essential functions as researchers are the pursuit and dissemination of knowledge through research and education. As scientists we
|
||||
seek to find models that reproduce the observations that we make in the world. In geophysics, we use inverse theory to mathematically
|
||||
create models of the earth from measured data. It is a difficult problem with many moving pieces: physics, discretization, simulation,
|
||||
regularization, optimization, computer science, linear algebra, geology. Exploring each of these disciplines can take a career, if you
|
||||
are so inclined, but as geophysicists we care about the combination: how to pull these disciplines together to answer our questions.
|
||||
This is the first problem we hope to help solve: to create a toolbox for the geophysicist that allows you to work at a high level and
|
||||
keep your geophysical question in focus. However, a toolbox is not enough. The research questions that we are interested in surround
|
||||
the integration of information to make better decisions.
|
||||
|
||||
We believe that the feedback loops in the geosciences could use some serious work. For example, collect multiple data-sets from the
|
||||
same field area (geology, seismic, electromagnetics, hydrogeology), process the data separately, and then reconvene with your
|
||||
multidisciplinary team. You may be rather surprised (or not) that the everyone has a (completely!?) different model. Dissonant at best,
|
||||
but often conflicting in the details. Therein lies the second problem: how do we integrate these geoscience fields? Not by force or
|
||||
even by default, but at least to have the option of quantitative communication and built in feedback loops. What we require is an
|
||||
implementation that is inherently and unequivocally modular, with all pieces available to manipulation. Black-box software, where the
|
||||
implementations are hidden, obfuscated, or difficult to manipulate, do not promote experimentation and investigation. We are working on
|
||||
a framework that exposes the details of the implementation to the geophysicist in a manner that promotes productivity and question
|
||||
based interrogation. This framework can be easily extended to encompass many geophysical problems and is built with the inverse problem
|
||||
as the fundamental goal.
|
||||
|
||||
The future we see is a mix of tools that span our disciplines, and a framework that allows us to integrate many different types of
|
||||
geophysical data so that we can communicate effectively and experiment efficiently. A toolbox combined with a framework that allows you
|
||||
to solve your own problems, and creates opportunities for us to work together to better image and understand the subsurface. What we
|
||||
are building is called SimPEG, simulation and parameter estimation in geophysics. We are building it in the open. We are testing it.
|
||||
Breaking it. Building it. Fixing it. Using it. If you believe, like we do, that geophysics can be more innovative and informative in
|
||||
the open and that these tools are necessary and invaluable in education as well as research, then you should get in touch. There is a
|
||||
lot of work to do!
|
||||
|
||||
The Big Picture
|
||||
---------------
|
||||
|
||||
Defining a well-posed inverse problem and solving it is a complex task that requires many components that must interact. It is helpful
|
||||
to view this task as a workflow in which various elements are explicitly identified and integrated. The figure below outlines the inversion components that consists of inputs, implementation, and evaluation. The inputs are composed of the geophysical data, the equations which are a mathematical description of the governing physics, and prior knowledge or assumptions about the setting. The implementation consists of two broad categories: the forward simulation and the inversion. The **forward simulation** is the means by which we solve the governing equations given a model and the **inversion components** evaluate and update this model. We are considering a gradient based approach, which updates the model through an optimization routine. The output of this implementation is a model, which, prior to interpretation, must be evaluated. This requires considering, and often re-assessing, the choices and assumptions made in both the input and implementation stages.
|
||||
|
||||
.. image:: ../../images/InversionWorkflow-PreSimPEG.png
|
||||
:width: 400 px
|
||||
:alt: Components
|
||||
:align: center
|
||||
|
||||
|
||||
A Comprehensive Framework
|
||||
-------------------------
|
||||
|
||||
There are an overwhelming amount of choices to be made as one works through the forward modeling and inversion process (see figure above). As a result, software implementations of this workflow often become complex and highly interdependent, making it difficult to interact with and to ask other scientists to pick up and change. Our approach to handling this complexity is to propose a framework, (see below), that compartmentalizes the implementation of inversions into various units. We present it in this specific modular style, as each unit contains a targeted subset of choices crucial to the inversion process.
|
||||
|
||||
.. image:: ../../images/InversionWorkflow.png
|
||||
:width: 400 px
|
||||
:alt: Framework
|
||||
:align: center
|
||||
|
||||
The process of obtaining an acceptable model from an inversion generally requires the geophysicist to perform several iterations of the inversion workflow, rethinking and redesigning each piece of the framework to ensure it is appropriate in the current context. Inversions are experimental and empirical by nature and our software package is designed to facilitate this iterative process. To accomplish this, we have divided the inversion methodology into eight major components (See figure above). The :class:`SimPEG.Mesh.BaseMesh.BaseMesh` class handles the discretization of the earth and also provides numerical operators. The forward simulation is split into two classes, the :class:`SimPEG.Survey.BaseSurvey` and the :class:`SimPEG.Problem.BaseProblem`. The :class:`SimPEG.Survey.BaseSurvey` class handles the geometry of a geophysical problem as well as sources. The :class:`SimPEG.Problem.BaseProblem` class handles the simulation of the physics for the geophysical problem of interest. Although created independently, these two classes must be paired to form all of the components necessary for a geophysical forward simulation and calculation of the sensitivity. The :class:`SimPEG.Problem.BaseProblem` creates geophysical fields given a source from the :class:`SimPEG.Survey.BaseSurvey`. The :class:`SimPEG.Survey.BaseSurvey` interpolates these fields to the receiver locations and converts them to the appropriate data type, for example, by selecting only the measured components of the field. Each of these operations may have associated derivatives with respect to the model and the computed field; these are included in the calculation of the sensitivity. For the inversion, a :class:`SimPEG.DataMisfit.BaseDataMisfit` is chosen to capture the goodness of fit of the predicted data and a :class:`SimPEG.Regularization.BaseRegularization` is chosen to handle the non-uniqueness. These inversion elements and an Optimization routine are combined into an inverse problem class :class:`SimPEG.InvProblem.BaseInvProblem`. :class:`SimPEG.InvProblem.BaseInvProblem` is the mathematical statement that will be numerically solved by running an Inversion. The :class:`SimPEG.Inversion.BaseInversion` class handles organization and dispatch of directives between all of the various pieces of the framework.
|
||||
|
||||
The arrows in the figure above indicate what each class takes as a primary argument. For example, both the :class:`SimPEG.Problem.BaseProblem` and :class:`SimPEG.Regularization.BaseRegularization` classes take a :class:`SimPEG.Mesh.BaseMesh.BaseMesh` class as an argument. The diagram does not show class inheritance, as each of the base classes outlined have many subtypes that can be interchanged. The :class:`SimPEG.Mesh.BaseMesh.BaseMesh` class, for example, could be a regular Cartesian mesh :class:`SimPEG.Mesh.TensorMesh` or a cylindrical coordinate mesh :class:`SimPEG.Mesh.CylMesh`, which have many properties in common. These common features, such as both meshes being created from tensor products, can be exploited through inheritance of base classes, and differences can be expressed through subtype polymorphism. Please look at the documentation here for more in-depth information.
|
||||
|
||||
|
||||
.. include:: ../../../CITATION.rst
|
||||
|
||||
Authors
|
||||
-------
|
||||
|
||||
.. include:: ../../../AUTHORS.rst
|
||||
|
||||
License
|
||||
-------
|
||||
|
||||
.. include:: ../../../LICENSE
|
||||
@@ -0,0 +1,80 @@
|
||||
.. _api_installing:
|
||||
|
||||
Getting Started with SimPEG
|
||||
***************************
|
||||
|
||||
Dependencies
|
||||
============
|
||||
|
||||
- Python 2.7
|
||||
- NumPy 1.8 (or greater)
|
||||
- SciPy 0.13 (or greater)
|
||||
- matplotlib 1.3 (or greater)
|
||||
- Cython 0.20 (or greater)
|
||||
|
||||
Installing Python
|
||||
=================
|
||||
|
||||
Python is available on all major operating systems, but if you are getting started with python
|
||||
it is best to use a package manager such as
|
||||
`Continuum Anaconda <https://store.continuum.io/cshop/anaconda/>`_ or
|
||||
`Enthought Canopy <https://www.enthought.com/products/canopy/>`_.
|
||||
You can download the package manager and use it to install the dependencies above.
|
||||
|
||||
|
||||
Installing SimPEG
|
||||
=================
|
||||
|
||||
SimPEG is on pip!!::
|
||||
|
||||
pip install SimPEG
|
||||
|
||||
|
||||
Installing from Source
|
||||
----------------------
|
||||
|
||||
First (you need git)::
|
||||
|
||||
git clone https://github.com/simpeg/simpeg
|
||||
|
||||
Second (from the root of the simpeg repository)::
|
||||
|
||||
python setup.py install
|
||||
|
||||
|
||||
Useful Links
|
||||
============
|
||||
An enormous amount of information (including tutorials and examples) can be found on the official websites of the packages
|
||||
|
||||
* `Python Website <http://www.python.org/>`_
|
||||
* `Numpy Website <http://www.numpy.org/>`_
|
||||
* `SciPy Website <http://www.scipy.org/>`_
|
||||
* `Matplotlib <http://matplotlib.org/>`_
|
||||
|
||||
Python for scientific computing
|
||||
-------------------------------
|
||||
|
||||
* `Python for Scientists <https://sites.google.com/site/pythonforscientists/>`_ Links to commonly used packages, Matlab to Python comparison
|
||||
* `Python Wiki <http://wiki.python.org/moin/NumericAndScientific>`_ Lists packages and resources for scientific computing in Python
|
||||
|
||||
Numpy and Matlab
|
||||
----------------
|
||||
|
||||
* `NumPy for Matlab Users <http://wiki.scipy.org/NumPy_for_Matlab_Users>`_
|
||||
* `Python vs Matlab <https://sites.google.com/site/pythonforscientists/python-vs-matlab>`_
|
||||
|
||||
Lessons in Python
|
||||
-----------------
|
||||
|
||||
* `Software Carpentry <http://software-carpentry.org/v4/python/index.html>`_
|
||||
* `Introduction to NumPy and Matplotlib <http://www.youtube.com/watch?v=3Fp1zn5ao2M>`_
|
||||
|
||||
Editing Python
|
||||
--------------
|
||||
|
||||
There are numerous ways to edit and test Python (see `PythonWiki <http://wiki.python.org/moin/PythonEditors>`_ for an overview) and in our group at least the following options are being used:
|
||||
|
||||
* `Sublime <http://www.sublimetext.com/>`_
|
||||
* `iPython Notebook <http://ipython.org/notebook.html>`_
|
||||
* `iPython <http://ipython.org/>`_
|
||||
* `Enthought Canopy <https://www.enthought.com/products/canopy/>`_
|
||||
Reference in New Issue
Block a user