Using finite volume

weakform-dc.md

: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, $\sigma$ and $\phi$ will be discretized on the cell-centers and the flux, $j$, will be on the faces. We will use the weak formulation to discretize the DC resistivity equation. {: icare="<5"}

from SimPEG import *
M = Mesh.TensorMesh([10,10])

Discretizing Equation

For equation \ref{dc2} the integration with a general face function $\vec{f}$ results in:

\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:

$$\left(a,b\right) = \int_\Omega{a \cdot b}{\partial v}$$

where $a$ and $b$ are either scalars or vectors.

---

Discretizing Equation \ref{dc2}: Left-Side

We will first look at the left-side of this equation. The matrial property $\boldsymbol{\sigma}$ and the flux $\mathbf{j}$ live in different locations on the mesh. Therefore, we have to do some linear interpolation for $\boldsymbol{\sigma}$.

$$ \newcommand{\MsigI}{\mathbf{M}_\frac{1}{\sigma}} \newcommand{\Pbcin}{\mathbf{P_{in}}} \newcommand{\Pbcout}{\mathbf{P_{out}}} $$

$$ \mathbf{j}^\top \text{diag}\left( {\mathbf{A}_v^{f2cc}}^\top \left( \mathbf{v} \circ \frac{1}{\boldsymbol{\sigma}} \right) \right) \mathbf{f} = \mathbf{f}^\top \MsigI \mathbf{j} $$

Here $\MsigI$ is the mass matrix that incorporates the conductivity structure and allows multiplication with $\mathbf{j}$ which lives on the faces.

$$ \MsigI

\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))

$\mathbf{A}_v^{f2cc}$ is an averaging matrix that takes you from faces to cell-centers, $\circ$ is a Hadamard product. A Hadamard product indicates point wise multiplication or equivalently:

$$\mathbf{a}\circ\mathbf{b} = \text{diag}(\mathbf{a})\mathbf{b} = \text{diag}(\mathbf{b})\mathbf{a}$$

---

Discretizing Equation \ref{dc1}: Right-Side

Here we will use integration by parts to simplify the right-side of equation \ref{dc2wf}:

$$\nabla\cdot(\phi\vec{f})=\nabla\phi\cdot\vec{f} + \phi\nabla\cdot\vec{f}$$

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, $\int_\Omega \nabla \cdot \mathbf{F} , dv = \oint_{\partial \Omega} \mathbf{F} \cdot \mathbf{n} , ds $, in the last integral:

$$\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 $\mathbf{B}$ must take into account the direction of the normal (which points out of the face), so we need some $\pm1$ in the right places. For a 1D mesh $\mathbf{B}$ would have the form:

$$ \mathbf{B} = \left[ \begin{array}{cc} -1 & 0\\ 0 & 0\\ \vdots & \vdots\\ 0 & 0\\ 0 & 1\\ \end{array} \right] $$

{: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}

---

Discretizing Equation 2: Putting it Together

The two sides of the equations can be combined and generalized by removing the face-function, $\mathbf{f}^\top$, from all terms.

\begin{equation} \label{discDC2} \MsigI \mathbf{j} = - \mathbf{D}^\top\text{diag}(\mathbf{v})\boldsymbol{\phi} + \mathbf{P}{bc}\boldsymbol{\phi}{bc} \end{equation}

---

Discretizing Equation \ref{dc1}

For equation \ref{dc1} the integration with a general cell function $w$ results in:

\begin{equation} \label{dc1wf} \left(\nabla\cdot \vec{j}, w\right) = \left(q, w\right) \end{equation}

Or equivalently:

$$\int \nabla\cdot \vec{j} w dv = \int q w dv$$

Discretizing the divergence operator as the matrix $\mathbf{D}$ this results in the following:

$$\mathbf{j}^\top\mathbf{D}^\top\text{diag}(\mathbf{v})\mathbf{w} = \mathbf{q}^\top\text{diag}(\mathbf{v})\mathbf{w}$$

or the transpose (which is the same because it is a scalar!):

$$\mathbf{w}^\top\text{diag}(\mathbf{v})\mathbf{D}\mathbf{j} = \mathbf{w}^\top\text{diag}(\mathbf{v})\mathbf{q}$$

This can be further simplified by eliminating the $\mathbf{w}^\top\text{diag}(\mathbf{v})$ from both sides of the equations.

The discretized matrix equation is very similar to the original equation \ref{dc1}:

$$\mathbf{D}\mathbf{j} = \mathbf{q} \label{dc1disc}$$

However, we may want to enforce Neumann boundary conditions on $\mathbf{j}$, which requires a few tweaks in our matrix equations. We can do this using projection matrices!

D = M.faceDiv

---

Boundary conditions on $j$

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.

Figure 1: Small one dimensional mesh, showing location of $\mathbf{j}$, $\boldsymbol{\sigma}$, and $\boldsymbol{\phi}$

The discretization for the divergence for this example is:

$$ \mathbf{D}\mathbf{j} = \frac{1}{h}\left[ \begin{array}{c|cc|c} -1 & 1 & & \\ & -1 & 1 & \\ & & -1 & 1 \end{array} \right] \left[ \begin{array}{c} j_0\j_1\j_2\j_3 \end{array} \right] $$

Where $h$ is the cell width, which in this case is constant and can be separated from the stencil matrix of $\pm 1$. Notice here that we can separate out the fluxes on the boundaries.

$$ \mathbf{D}\mathbf{j} = \frac{1}{h}\left[ \begin{array}{cc} 1 & \\ -1 & 1\\ & -1 \end{array} \right] \left[ \begin{array}{c} j_1\j_2 \end{array} \right] + \frac{1}{h}\left[ \begin{array}{cc} -1 & \\ & \\ & 1 \end{array} \right] \left[ \begin{array}{c} j_0\j_3 \end{array} \right] $$

This can be written similarly using projection matrices (reduced identities):

$$ \mathbf{Dj} = \mathbf{D}\Pbcin^\top\Pbcin\mathbf{j} + \mathbf{D}\Pbcout^\top\mathbf{j}_{bc} $$

Where $\mathbf{j}_{bc} = \Pbcout\mathbf{j}$ and the projection matrices for the 1D example are:

$$ \Pbcin = \left[ \begin{array}{cccc} 0 & 1 & 0 & 0\\ 0 & 0 & 1 & 0\\ \end{array} \right], \qquad \Pbcout = \left[ \begin{array}{cccc} 1 & 0 & 0 & 0\\ 0 & 0 & 0 & 1\\ \end{array} \right] $$

$$\label{discDC1} \mathbf{D}\Pbcin^\top\Pbcin\mathbf{j} + \mathbf{D}\Pbcout^\top\mathbf{j}_{bc} = \mathbf{q}$$

---

Putting the Equations Together

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 ($\phi = 0 \in \partial \Omega$). In the DC resistivity problem, however, we are often interested in homogeneous Neumann boundary conditions ($j = 0 \in \partial \Omega$).

In equation \ref{discDC2} we need to make sure that the divergence that we apply on the general face function $\mathbf{f}^\top$ match those of the DC resistivity equations. Here we will rewrite the integration and make a few notes:

$$\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 $\vec{f}\cdot\vec{n} = 0 \in \partial\Omega$ and the $\nabla\cdot\vec{f}$ must enforce the correct boundary conditions. In the discretized form this translates to using the projected divergence matrix and dropping the boundary condition term for $\boldsymbol{\phi}_{bc}$. We can plug the projected divergence matrix into equation \ref{discDC2} and rearrange to solve for $\mathbf{j}$:

$$\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{D}\Pbcin^\top\Pbcin$:

$$ -\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 $\mathbf{j}_{bc}$ is zero on all of the boundaries, so we will drop that term, and multiply on the left by a $\text{diag}(\mathbf{v})$ to make the matrix symmetric.

$$\text{diag}(\mathbf{v})\mathbf{D}_{in}\mathbf{M}_{\frac{1}{\sigma}}^{-1}\mathbf{D}_{in}^\top\text{diag}(\mathbf{v})\boldsymbol{\phi} = -\text{diag}(\mathbf{v})\mathbf{q}$$

{: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])