<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Building our operators
Step1:  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.
Step2:  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
Step3:  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
Step4:  Working in Sparse
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.
Step5:  In SimPEG, this is stored as the faceDiv property on the mesh
Step6:  Moving to 2D
To move up in dimensionality, we build a 2D mesh which has both x and y faces
Step7:  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!
Step8:  But first... what does the matrix look like?
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", 
- 2 that are due to the x-contribution
- 2 that are due to the y-contribution
Here, we define a small 2D mesh so that it is easier to see the matrix structure.
Step9:  Now, lets look at the matrix structure
Step10:  How did we construct the matrix? - Kronecker products. 
There is a handy identity that relates the vectorized face function to its matrix form (<a href = "https
Step11:  Now that we have a discrete divergence, lets check out the divergence of  the face function we defined earlier.
Step12:  Are we right??
Since we chose a simple function,
$$
\vec{j} = - \sin(2\pi x) \hat{x} - \sin(2\pi y) \hat{y} 
$$
we know the continuous divergence...
$$
\nabla \cdot \vec{j} = -2\pi (\cos(2\pi x) + \cos(2\pi y))
$$
So lets plot it and take a look
Step13:  Those look similar 
<ASSISTANT_TASK:>
Python Code:
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!
# 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')
# 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');
Explanation: Building our operators: the Face Divergence
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.
<img src="./images/Divergence.png" width=80% align="center">
<h4 align="center">Figure 4. Geometrical definition of the divergence and the discretization.</h4>
Implementation
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.
End of explanation
# 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))
Explanation: 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.
End of explanation
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))
Explanation: 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!)
End of explanation
"There are {nnz} zeros (too many!) that we are storing".format(nnz = np.sum(faceDiv == 0))
Explanation: 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
End of explanation
import scipy.sparse as sp
from SimPEG.Utils import sdiag # we are often building sparse diagonal matrices, so we made a functio in SimPEG!
# 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))
Explanation: Working in Sparse
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.
End of explanation
print(mesh1D.faceDiv * face_vec) # and still gives us the same answer!
Explanation: In SimPEG, this is stored as the faceDiv property on the mesh
End of explanation
mesh2D = Mesh.TensorMesh([100,80])
mesh2D.plotGrid()
plt.axis('tight');
Explanation: Moving to 2D
To move up in dimensionality, we build a 2D mesh which has both x and y faces
End of explanation
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])
Explanation: 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!
End of explanation
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,:]))
Explanation: But first... what does the matrix look like?
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", 
- 2 that are due to the x-contribution
- 2 that are due to the y-contribution
Here, we define a small 2D mesh so that it is easier to see the matrix structure.
End of explanation
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));
Explanation: Now, lets look at the matrix structure
End of explanation
# small_mesh2D.faceDiv?? # check out the code!
Explanation: How did we construct the matrix? - Kronecker products. 
There is a handy identity that relates the vectorized face function to its matrix form (<a href = "https://en.wikipedia.org/wiki/Vectorization_(mathematics)#Compatibility_with_Kronecker_products">wikipedia link!</a>)
$$
\text{vec}(AUB^\top) = (B \otimes A) \text{vec}(U)
$$
For the x-contribution:
- A is our 1D differential operator ([-1, +1] on the diagonals)
- U is $j_x$ (the x-face function as a matrix) 
- B is just an identity
so 
$$
\text{Div}x \text{vec}(j_x) = (I \otimes Div{1D}) \text{vec}(j_x)
$$
For the y-contribution: 
- A is just an identity!
- U is $j_y$ (the y-face function as a matrix) 
- B is our 1D differential operator ([-1, +1] on the diagonals)
so
$$
\text{Div}y \text{vec}(j_y) = (\text{Div}{1D} \otimes I) \text{vec}(j_y)
$$
$$
\text{Div} \cdot j  = \text{Div}_x \cdot j_x + \text{Div}_y \cdot j_y = [\text{Div}_x, \text{Div}_y] \cdot [j_x; j_y]
$$
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??
End of explanation
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()
Explanation: Now that we have a discrete divergence, lets check out the divergence of  the face function we defined earlier.
End of explanation
# 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()
Explanation: Are we right??
Since we chose a simple function,
$$
\vec{j} = - \sin(2\pi x) \hat{x} - \sin(2\pi y) \hat{y} 
$$
we know the continuous divergence...
$$
\nabla \cdot \vec{j} = -2\pi (\cos(2\pi x) + \cos(2\pi y))
$$
So lets plot it and take a look
End of explanation
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 );
Explanation: Those look similar :)
Order Test
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!
End of explanation
<END_TASK>