Weakform DC Resistivity Derivation

For equation (1) the integration with a general face function $\vec{f}$ results in:

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.

\left(\frac{1}{\sigma}\vec{j},\vec{f}\right) = \left(\nabla\phi,\vec{f}\right)

Discretizing Equation (1): Left-Side

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

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

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

\mathbf{M}_\frac{1}{\sigma} = \text{diag}\left( {\mathbf{A}_v^{f2cc}}^\top \left( \mathbf{v} \circ \frac{1}{\sigma} \right) \right) Msig = sdiag(mesh.aveF2CC.T*(mesh.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 (1): Right-Side

Here we will use integration by parts to simplify the right-side of the weakform equation:

\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 Divergence theorem: \int_\Omega \nabla \cdot \mathbf{F} \, dv = \oint_{\partial \Omega} \mathbf{F} \cdot \mathbf{n} \, ds

We will use the Divergence theorem 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:

-\phi^\top\text{diag}(\mathbf{v})\mathbf{D}\mathbf{f} + \phi_{bc}^\top\text{diag}(\mathbf{a}_{bc}) \mathbf{B}\mathbf{f}

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: Note: In 1D 'area' doesn't really make sense and can be substituted for a vector of ones. \mathbf{B} = \left[ \begin{array}{cc} -1 & 0\\ 0 & 0\\ \vdots & \vdots\\ 0 & 0\\ 0 & 1 \end{array} \right] 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!):

-\mathbf{f}^\top \mathbf{D}^\top\text{diag}(\mathbf{v})\phi + \mathbf{f}^\top \mathbf{P}_{bc}\phi_{bc}

Discretizing Equation (1): 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.

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. 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. Small one dimensional mesh, showing location of $\mathbf{j}$, $\sigma$, and $\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}\mathbf{P}_{in}^\top\mathbf{P}_{in}\mathbf{j} + \mathbf{D}\mathbf{P}_{out}^\top\mathbf{j}_{bc} Where $\mathbf{j}_{bc} = \mathbf{P}_{out}\mathbf{j}$ and the projection matrices for the 1D example are: \mathbf{P}_{in} = \left[ \begin{array}{cccc} 0 & 1 & 0 & 0\\ 0 & 0 & 1 & 0 \end{array} \right], \qquad \mathbf{P}_{out} = \left[ \begin{array}{cccc} 1 & 0 & 0 & 0\\ 0 & 0 & 0 & 1 \end{array} \right]