:uid: weakform-dc :title: Discretizing DC Resistivity :description: Using finite volume :tooltip: None :tag: :group: :license: CC-BY-4.0 :source: https://api.github.com/gists/fe7b1d8a3b905bab158d
We will start with the formulation of the Direct Current (DC) resistivity problem in geophysics.
\begin{equation} \label{dc2} \frac{1}{\sigma}\vec{j} = \nabla \phi \end{equation}
\begin{equation} \label{dc1} \nabla\cdot \vec{j} = q \end{equation}
In the following discretization,
from SimPEG import *
M = Mesh.TensorMesh([10,10])
For equation \ref{dc2} the integration with a general face function
\begin{equation} \label{dc2wf} \left(\frac{1}{\sigma}\vec{j},\vec{f}\right) = \left(\nabla\phi,\vec{f}\right) \end{equation}
Where the inner product is defined as:
where
We will first look at the left-side of this equation. The matrial property
$$ \newcommand{\MsigI}{\mathbf{M}\frac{1}{\sigma}} \newcommand{\Pbcin}{\mathbf{P{in}}} \newcommand{\Pbcout}{\mathbf{P_{out}}} $$
Here
\text{diag}\left( {\mathbf{A}_v^{f2cc}}^\top \left( \mathbf{v} \circ \frac{1}{\boldsymbol{\sigma}} \right) \right) $$
Or in code:
Msig = sdiag(M.aveF2CC.T*(M.vol*(1/sigma))
Here we will use integration by parts to simplify the right-side of equation \ref{dc2wf}:
Integrating and rearranging gives:
$$\int \nabla \phi \cdot \vec{f} dv =
- \int \phi \nabla \cdot \vec{f} dv
- \int \nabla \cdot(\phi \vec{f}) dv $$
We will use the Divergence theorem,
$$\left(\nabla\phi,\vec{f}\right) = \int \nabla \phi \cdot \vec{f} dv =
- \int \phi \nabla \cdot \vec{f} dv
- \oint_{\partial \Omega} \phi \vec{f}\cdot\vec{n} ds $$
Where the discretization is:
\begin{equation} -\boldsymbol{\phi}^\top\text{diag}(\mathbf{v})\mathbf{D}\mathbf{f} + \boldsymbol{\phi}{bc}^\top\text{diag}(\mathbf{a}{bc}) \mathbf{B}\mathbf{f} \end{equation}
Here
{:note: In 1D area doesn't really make sense and can be substituted for a vector of ones.}
here we will define $\mathbf{P}{bc} = \mathbf{B}^\top \text{diag}(\mathbf{a}{bc})$ and take the transpose (which is the same because it is a scalar!):
\begin{equation} -\mathbf{f}^\top \mathbf{D}^\top\text{diag}(\mathbf{v})\boldsymbol{\phi} + \mathbf{f}^\top \mathbf{P}{bc}\boldsymbol{\phi}{bc} \end{equation}
The two sides of the equations can be combined and generalized by removing the face-function,
\begin{equation} \label{discDC2} \MsigI \mathbf{j} = - \mathbf{D}^\top\text{diag}(\mathbf{v})\boldsymbol{\phi} + \mathbf{P}{bc}\boldsymbol{\phi}{bc} \end{equation}
For equation \ref{dc1} the integration with a general cell function
\begin{equation} \label{dc1wf} \left(\nabla\cdot \vec{j}, w\right) = \left(q, w\right) \end{equation}
Or equivalently:
Discretizing the divergence operator as the matrix
or the transpose (which is the same because it is a scalar!):
This can be further simplified by eliminating the
The discretized matrix equation is very similar to the original equation \ref{dc1}:
However, we may want to enforce Neumann boundary conditions on
D = M.faceDiv
Boundary conditions are often important to consider, in a DC resistivity survey we may create a domain that is sufficiently padded to assume that at the boundaries the flux into or out of the domain is zero. If this is the case, we can reduce the number of unknowns in equation \ref{dc1disc}. To illustrate the separation of the fluxes on the boundary $\mathbf{j}{bc}$ and the fluxes inside the domain $\mathbf{j}{in}$, we will use a simple 1D example.
The discretization for the divergence for this example is:
Where
This can be written similarly using projection matrices (reduced identities):
Where
Care needs to be taken when putting these equations together to ensure that the boundary conditions match. The natural boundary conditions that occur in this system using a cell-centered discretization are homogeneous Dirichlet (
In equation \ref{discDC2} we need to make sure that the divergence that we apply on the general face function
$$\int \nabla \phi \cdot \vec{f} dv =
- \int \phi \nabla \cdot \vec{f} dv
- \oint_{\partial \Omega} \phi \vec{f}\cdot\vec{n} ds $$
In this case we note that
$$\mathbf{j} =
- \mathbf{M}_{\frac{1}{\sigma}}^{-1}\Pbcin^\top\Pbcin\mathbf{D}^\top\text{diag}(\mathbf{v})\boldsymbol{\phi}$$
We can substitute this into equation \ref{discDC2} by defining
$$ -\mathbf{D}{in}\mathbf{M}{\frac{1}{\sigma}}^{-1}\mathbf{D}{in}^\top\text{diag}(\mathbf{v})\boldsymbol{\phi} + \mathbf{D}\Pbcout^\top\mathbf{j}{bc} = \mathbf{q} $$
Here we know
{:note: The matrix has a null space of a constant vector! This is expected with homogeneous Neumann boundary conditions, but we need to not forget about this when we are solving the matrix system.}
from SimPEG import *
M = Mesh.TensorMesh([60,60,10])
Pbc, Pin, Pout = M.getBCProjWF('neumann')
V = Utils.sdiag(M.vol)
Msig = M.getFaceInnerProduct()
MsigI = Utils.sdInv(Msig)
Din = M.faceDiv*Pin.T*Pin
A = V*Din*MsigI*Din.T*V
q = np.zeros(M.vnC)
q[[30,30],[20,40],5] = [-1,1]
q = -V*Utils.mkvc(q)
phi = Solver(A).solve(q)
# Plot it!
figsize(10,6)
M.plotSlice(-Din.T*phi,'F',normal='Z',view='vec')
xlim([0,1]); ylim([0,1])