This notebook uses Python 2.7 and the open source package SimPEG. SimPEG can be installed using the python package manager PyPi and running:
pip install SimPEGAlternatively, these notebooks can be run on the web using binders
This tutorial consists of 3 parts, here, we introduce the problem, in divergence.ipynb we build the discrete divergence operator and in weakformulation.ipynb, we discretize and solve the DC equations using weak formulation.
Contents
DC resistivity surveys obtain information about subsurface electrical conductivity, $\sigma$. This physical property is often diagnostic in mineral exploration, geotechnical, environmental and hydrogeologic problems, where the target of interest has a significant electrical conductivity contrast from the background. In a DC resistivity survey, steady state currents are set up in the subsurface by injecting current through a positive electrode and completing the circuit with a return electrode.
Conservation of charge (which can be derived by taking the divergence of Ampere’s law at steady state) connects the divergence of the current density everywhere in space to the source term which consists of two point sources, one positive and one negative. The flow of current sets up electric fields according to Ohm’s law, which relates current density to electric fields through the electrical conductivity. From Faraday’s law for steady state fields, we can describe the electric field in terms of a scalar potential, $\phi$, which we sample at potential electrodes to obtain data in the form of potential differences.
Where are we going??
Here, we are going to do a run through of how to setup and solve the DC resistivity equations for a 2D problem using SimPEG. This is meant to give you a once-over of the whole picture. We will break down the steps to get here in the series of notebooks that follow...
In [1]:
# Import numpy, python's n-dimensional array package,
# the mesh class with differential operators from SimPEG
# matplotlib, the basic python plotting package
import numpy as np
from SimPEG import Mesh, Utils
import matplotlib.pyplot as plt
%matplotlib inline
plt.set_cmap(plt.get_cmap('viridis')) # use a nice colormap!
Where we solve things! See mesh.ipynb a discussion of how we construct a mesh and the associated properties we need.
In [2]:
# Define a unit-cell mesh
mesh = Mesh.TensorMesh([100, 80]) # setup a mesh on which to solve
print("The mesh has {nC} cells.".format(nC=mesh.nC))
mesh.plotGrid()
plt.axis('tight');
In [ ]:
In [3]:
# model parameters
sigma_background = 1. # Conductivity of the background, S/m
sigma_block = 10. # Conductivity of the block, S/m
# add a block to our model
x_block = np.r_[0.4, 0.6]
y_block = np.r_[0.4, 0.6]
# assign them on the mesh
sigma = sigma_background * np.ones(mesh.nC) # create a physical property model
block_indices = ((mesh.gridCC[:,0] >= x_block[0]) & # left boundary
(mesh.gridCC[:,0] <= x_block[1]) & # right boundary
(mesh.gridCC[:,1] >= y_block[0]) & # bottom boundary
(mesh.gridCC[:,1] <= y_block[1])) # top boundary
# add the block to the physical property model
sigma[block_indices] = sigma_block
# plot it!
plt.colorbar(mesh.plotImage(sigma)[0])
plt.title('electrical conductivity, $\sigma$')
Out[3]:
In [4]:
# Define a source
a_loc, b_loc = np.r_[0.2, 0.5], np.r_[0.8, 0.5]
source_locs = [a_loc, b_loc]
# locate it on the mesh
source_loc_inds = Utils.closestPoints(mesh, source_locs)
a_loc_mesh = mesh.gridCC[source_loc_inds[0],:]
b_loc_mesh = mesh.gridCC[source_loc_inds[1],:]
# plot it
plt.colorbar(mesh.plotImage(sigma)[0])
plt.plot(a_loc_mesh[0], a_loc_mesh[1],'wv', markersize=8) # a-electrode
plt.plot(b_loc_mesh[0], b_loc_mesh[1],'w^', markersize=8) # b-electrode
plt.title('electrical conductivity, $\sigma$')
Out[4]:
How we construct the divergence operator is discussed in divergence.ipynb, and the inner product matrix in weakformulation.ipynb. The final system is assembled and discussed in play.ipynb (with widgets!).
In [5]:
mesh.faceDiv??
In [6]:
# Assemble and solve the DC resistivity problem
Div = mesh.faceDiv
Sigma = mesh.getFaceInnerProduct(sigma, invProp=True, invMat=True)
Vol = Utils.sdiag(mesh.vol)
# assemble the system matrix
A = Vol * Div * Sigma * Div.T * Vol
# right hand side
q = np.zeros(mesh.nC)
q[source_loc_inds] = np.r_[+1, -1]
In [7]:
from SimPEG import Solver # import the default solver (LU)
In [8]:
# solve the DC resistivity problem
Ainv = Solver(A) # create a matrix that behaves like A inverse
phi = Ainv * q
In [10]:
# look at the results!
plt.colorbar(mesh.plotImage(phi)[0])
plt.title('Electric Potential, $\phi$');
In the notebooks that follow, we will