The divergence is the integral of a flux through a closed surface as that enclosed volume shrinks to a point. Since we have discretized and no longer have continuous functions, we cannot fully take the limit to a point; instead, we approximate it around some (finite!) volume: a cell. The flux out of the surface ($\vec{j} \cdot \vec{n}$) is actually how we discretized $\vec{j}$ onto our mesh (i.e. $\bf{j}$) except that the face normal points out of the cell (rather than in the axes direction). After fixing the direction of the face normal (multiplying by $\pm 1$), we only need to calculate the face areas and cell volume to create the discrete divergence matrix.
Although this is a really helpful way to think about conceptually what is happening, the implementation of that would be a huge for loop over each cell. In practice, this would be slow, so instead, we will take advantage of linear algebra. Let's start by looking at this in 1 dimension using the SimPEG Mesh class.
In [1]:
import numpy as np
from SimPEG import Mesh
import matplotlib.pyplot as plt
%matplotlib inline
plt.set_cmap(plt.get_cmap('viridis')) # use a nice colormap!
In [2]:
# define a 1D mesh
mesh1D = Mesh.TensorMesh([5]) # with 5 cells
fig, ax = plt.subplots(1,1, figsize=(12,2))
ax.plot(mesh1D.gridN, np.zeros(mesh1D.nN),'-k',marker='|',markeredgewidth=2, markersize=16)
ax.plot(mesh1D.gridCC,np.zeros(mesh1D.nC),'o')
ax.plot(mesh1D.gridFx,np.zeros(mesh1D.nFx),'>')
ax.set_title('1D Mesh')
Out[2]:
In [3]:
# and define a vector of fluxes that live on the faces of the 1D mesh
face_vec = np.r_[0., 1., 2., 2., 1., 0.] # vector of fluxes that live on the faces of the mesh
print("The flux on the faces is {}".format(face_vec))
plt.plot(mesh1D.gridFx, face_vec, '-o')
plt.ylim([face_vec.min()-0.5, face_vec.max()+0.5])
plt.grid(which='both')
plt.title('face_vec');
Over a single cell, the divergence is
$$ \nabla \cdot \vec{j}(p) = \lim_{v \to \{p\}} = \int \int_{S(v)} \frac{\vec{j}\cdot \vec{n}}{v} dS $$in 1D, this collapses to taking a single difference - how much is going out of the cell vs coming in?
$$ \nabla \cdot \vec{j} \approx \frac{1}{v}(-j_{\text{left}} + j_{\text{right}}) $$Since the normal of the x-face on the left side of the cell points in the positive x-direction, we multiply by -1 to get the flux going out of the cell. On the right, the normal defining the x-face is point out of the cell, so it is positive.
In [4]:
# We can take the divergence over the entire mesh by looping over each cell
div_face_vec = np.zeros(mesh1D.nC) # allocate for each cell
for i in range(mesh1D.nC): # loop over each cell and
div_face_vec[i] = 1.0/mesh1D.vol[i] * (-face_vec[i] + face_vec[i+1])
print("The face div of the 1D flux is {}".format(div_face_vec))
Doing it as a for loop is easy to program for the first time, but is difficult to see what is going on and could be slow! Instead, we can build a faceDiv matrix (note: this is a silly way to do this!)
In [5]:
faceDiv = np.zeros([mesh1D.nC, mesh1D.nF]) # allocate space for a face div matrix
for i in range(mesh1D.nC): # loop over each cell
faceDiv[i, [i, i+1]] = 1.0/mesh1D.vol[i] * np.r_[-1,+1]
print("The 1D face div matrix for this mesh is \n{}".format(faceDiv))
assert np.all( faceDiv.dot(face_vec) == div_face_vec ) # make sure we get the same result!
print("\nThe face div of the 1D flux is still {}!".format(div_face_vec))
the above is still a loop... (and python is not a fan of loops). Also, if the mesh gets big, we are storing a lot of unnecessary zeros
In [6]:
"There are {nnz} zeros (too many!) that we are storing".format(nnz = np.sum(faceDiv == 0))
Out[6]:
We will use instead sparse matrices instead. These are in scipy and act almost the same as numpy arrays (except they default to matrix multiplication), and they don't store all of those pesky zeros! We use scipy.sparse to build these matrices.
In [7]:
import scipy.sparse as sp
from SimPEG.Utils import sdiag # we are often building sparse diagonal matrices, so we made a functio in SimPEG!
In [8]:
# construct differencing matrix with diagonals -1, +1
sparse_diff = sp.spdiags((np.ones((mesh1D.nC+1, 1))*[-1, 1]).T, [0, 1], mesh1D.nC, mesh1D.nC+1, format="csr")
print("the sparse differencing matrix is \n{}".format(sparse_diff.todense()))
# account for the volume
faceDiv_sparse = sdiag(1./mesh1D.vol) * sparse_diff # account for volume
print("\n and the face divergence is \n{}".format(faceDiv_sparse.todense()))
print("\n but now we are only storing {nnz} nonzeros".format(nnz=faceDiv_sparse.nnz))
assert np.all(faceDiv_sparse.dot(face_vec) == div_face_vec)
print("\n and we get the same answer! {}".format(faceDiv_sparse * face_vec))
In SimPEG, this is stored as the faceDiv property on the mesh
In [9]:
print(mesh1D.faceDiv * face_vec) # and still gives us the same answer!
To move up in dimensionality, we build a 2D mesh which has both x and y faces
In [10]:
mesh2D = Mesh.TensorMesh([100,80])
mesh2D.plotGrid()
plt.axis('tight');
We define 2 face functions, one in the x-direction and one in the y-direction. Here, we choose to work with sine functions as the continuous divergence is easy to compute, meaning we can test it!
In [11]:
jx_fct = lambda x, y: -np.sin(2.*np.pi*x)
jy_fct = lambda x, y: -np.sin(2.*np.pi*y)
jx_vec = jx_fct(mesh2D.gridFx[:,0], mesh2D.gridFx[:,1])
jy_vec = jy_fct(mesh2D.gridFy[:,0], mesh2D.gridFy[:,1])
j_vec = np.r_[jx_vec, jy_vec]
print("There are {nFx} x-faces and {nFy} y-faces, so the length of the "
"face function, j, is {lenj}".format(
nFx=mesh2D.nFx,
nFy=mesh2D.nFy,
lenj=len(j_vec)
))
plt.colorbar(mesh2D.plotImage(j_vec, 'F', view='vec')[0])
Out[11]:
Now, we know that we do not want to loop over each of the cells and instead want to work with matrix-vector products. In this case, each row of the divergence matrix should pick out the two relevant faces in the x-direction and two in the y-direction (4 total).
When we unwrap our face function, we unwrap using column major ordering, so all of the x-faces are adjacent to one another, while the y-faces are separated by the number of cells in the x-direction (see mesh.ipynb for more details!).
When we plot the divergence matrix, there will be 4 "diagonals",
Here, we define a small 2D mesh so that it is easier to see the matrix structure.
In [12]:
small_mesh2D = Mesh.TensorMesh([3,4])
print("Each y-face is {} entries apart".format(small_mesh2D.nCx))
print("and the total number of x-faces is {}".format(small_mesh2D.nFx))
print("So in the first row of the faceDiv, we have non-zero entries at \n{}".format(
small_mesh2D.faceDiv[0,:]))
Now, lets look at the matrix structure
In [13]:
fig, ax = plt.subplots(1,2, figsize=(12,4))
# plot the non-zero entries in the faceDiv
ax[0].spy(small_mesh2D.faceDiv, ms=2)
ax[0].set_xlabel('2D faceDiv')
small_mesh2D.plotGrid(ax=ax[1])
# Number the faces and plot. (We should really add this to SimPEG... pull request anyone!?)
xys = zip(
small_mesh2D.gridFx[:,0],
small_mesh2D.gridFx[:,1],
range(small_mesh2D.nFx)
)
for x,y,ii in xys:
ax[1].plot(x, y, 'r>')
ax[1].text(x+0.01, y-0.02, ii, color='r')
xys = zip(
small_mesh2D.gridFy[:,0],
small_mesh2D.gridFy[:,1],
range(small_mesh2D.nFy)
)
for x,y,ii in xys:
ax[1].plot(x, y, 'g^')
ax[1].text(x-0.02, y+0.02, ii+small_mesh2D.nFx, color='g')
ax[1].set_xlim((-0.1,1.1));
ax[1].set_ylim((-0.1,1.1));
How did we construct the matrix? - Kronecker products. There is a handy identity that relates the vectorized face function to its matrix form (wikipedia link!) $$ \text{vec}(AUB^\top) = (B \otimes A) \text{vec}(U) $$
For the x-contribution:
For the y-contribution:
And $j$ is just $[j_x; j_y]$, so we can horizontally stack $\text{Div}_x$, $\text{Div}_y$
$$ \text{Div} = [\text{Div}_x, \text{Div}_y] $$You can check this out in the SimPEG docs by running small_mesh2D.faceDiv??
In [14]:
# small_mesh2D.faceDiv?? # check out the code!
Now that we have a discrete divergence, lets check out the divergence of the face function we defined earlier.
In [15]:
Div_j = mesh2D.faceDiv * j_vec
fig, ax = plt.subplots(1,2, figsize=(8,4))
plt.colorbar(mesh2D.plotImage(j_vec, 'F', view='vec', ax=ax[0])[0],ax=ax[0])
plt.colorbar(mesh2D.plotImage(Div_j, ax=ax[1])[0],ax=ax[1])
ax[0].set_title('j')
ax[1].set_title('Div j')
plt.tight_layout()
In [16]:
# from earlier
# jx_fct = lambda x, y: -np.sin(2*np.pi*x)
# jy_fct = lambda x, y: -np.sin(2*np.pi*y)
sol = lambda x, y: -2*np.pi*(np.cos(2*np.pi*x)+np.cos(2*np.pi*y))
cont_div_j = sol(mesh2D.gridCC[:,0], mesh2D.gridCC[:,1])
Div_j = mesh2D.faceDiv * j_vec
fig, ax = plt.subplots(1,2, figsize=(8,4))
plt.colorbar(mesh2D.plotImage(Div_j, ax=ax[0])[0],ax=ax[0])
plt.colorbar(mesh2D.plotImage(cont_div_j, ax=ax[1])[0],ax=ax[1])
ax[0].set_title('Discrete Div j')
ax[1].set_title('Continuous Div j')
plt.tight_layout()
Those look similar :)
We can do better than just an eye-ball comparison - since we are using a a staggered grid, with centered differences, the discretization should be second-order ($\mathcal{O}(h^2)$). That is, as we refine the mesh, our approximation of the divergence should improve by a factor of 2.
SimPEG has a number of testing functions for derivatives and order of convergence to make our lives easier!
In [17]:
import unittest
from SimPEG.Tests import OrderTest
jx = lambda x, y: -np.sin(2*np.pi*x)
jy = lambda x, y: -np.sin(2*np.pi*y)
sol = lambda x, y: -2*np.pi*(np.cos(2*np.pi*x)+np.cos(2*np.pi*y))
class Testify(OrderTest):
meshDimension = 2
def getError(self):
j = np.r_[jx(self.M.gridFx[:,0], self.M.gridFx[:,1]),
jy(self.M.gridFy[:,0], self.M.gridFy[:,1])]
num = self.M.faceDiv * j # numeric answer
ans = sol(self.M.gridCC[:,0], self.M.gridCC[:,1]) # note M is a 2D mesh
return np.linalg.norm((num - ans), np.inf) # look at the infinity norm
# (as we refine the mesh, the number of cells
# changes, so need to be careful if using a 2-norm)
def test_order(self):
self.orderTest()
# This just runs the unittest:
suite = unittest.TestLoader().loadTestsFromTestCase( Testify )
unittest.TextTestRunner().run( suite );
Looks good - Second order convergence!
In the next notebook, we will explore how to use the weak formulation to discretize the DC equations.