Note
Go to the end to download the full example code.
Differential Operators#
When solving PDEs using the finite volume approach, inner products may
contain differential operators. Where
In this section, we demonstrate how to go from the inner product to the discrete approximation for each case. In doing so, we must construct discrete differential operators, inner product matricies and consider boundary conditions.
Import Packages#
Here we import the packages required for this tutorial
from discretize.utils import sdiag
from discretize import TensorMesh
import numpy as np
import matplotlib.pyplot as plt
Gradient#
Where
Inner Product at edges:
In the case that
Inner Product at faces:
In the case that
# Make basic mesh
h = np.ones(10)
mesh = TensorMesh([h, h, h])
# Items required to perform u.T*(Me*Gn*phi)
Me = mesh.get_edge_inner_product() # Basic inner product matrix (edges)
Gn = mesh.nodal_gradient # Nodes to edges gradient
# Items required to perform u.T*(Mf*Gc*phi)
Mf = mesh.get_face_inner_product() # Basic inner product matrix (faces)
mesh.set_cell_gradient_BC(
["neumann", "dirichlet", "neumann"]
) # Set boundary conditions
Gc = mesh.cell_gradient # Cells to faces gradient
# Plot Sparse Representation
fig = plt.figure(figsize=(5, 6))
ax1 = fig.add_subplot(121)
ax1.spy(Me * Gn, markersize=0.5)
ax1.set_title("Me*Gn")
ax2 = fig.add_subplot(122)
ax2.spy(Mf * Gc, markersize=0.5)
ax2.set_title("Mf*Gc")

Text(0.5, 1.0513468013468013, 'Mf*Gc')
Divergence#
Where
The divergence defines a measure of the flux leaving/entering a volume. As a
result, it is natural for
# Make basic mesh
h = np.ones(10)
mesh = TensorMesh([h, h, h])
# Items required to perform psi.T*(Mc*D*v)
Mc = sdiag(mesh.cell_volumes) # Basic inner product matrix (centers)
D = mesh.face_divergence # Faces to centers divergence
# Plot sparse representation
fig = plt.figure(figsize=(8, 5))
ax1 = fig.add_subplot(111)
ax1.spy(Mc * D, markersize=0.5)
ax1.set_title("Mc*D", pad=20)

Text(0.5, 1.0, 'Mc*D')
Curl#
Where
Inner Product at Faces:
Let
Inner Product at Edges:
Now let
# Make basic mesh
h = np.ones(10)
mesh = TensorMesh([h, h, h])
# Items required to perform u.T*(Mf*Ce*v)
Mf = mesh.get_face_inner_product() # Basic inner product matrix (faces)
Ce = mesh.edge_curl # Edges to faces curl
# Items required to perform u.T*(Me*Cf*v)
Me = mesh.get_edge_inner_product() # Basic inner product matrix (edges)
Cf = mesh.edge_curl.T # Faces to edges curl (assumes Dirichlet)
# Plot Sparse Representation
fig = plt.figure(figsize=(9, 5))
ax1 = fig.add_subplot(121)
ax1.spy(Mf * Ce, markersize=0.5)
ax1.set_title("Mf*Ce", pad=10)
ax2 = fig.add_subplot(122)
ax2.spy(Me * Cf, markersize=0.5)
ax2.set_title("Me*Cf", pad=10)

Text(0.5, 1.068020708880924, 'Me*Cf')
Scalar Laplacian#
Where
Using
In this case, the surface integral can be eliminated if we can assume a
Neumann condition of
Inner Prodcut at Edges:
Let
Inner Product at Faces:
Let
# Make basic mesh
h = np.ones(10)
mesh = TensorMesh([h, h, h])
# Items required to perform psi.T*(Gn.T*Me*Gn*phi)
Me = mesh.get_edge_inner_product() # Basic inner product matrix (edges)
Gn = mesh.nodal_gradient # Nodes to edges gradient
# Items required to perform psi.T*(Gc.T*Mf*Gc*phi)
Mf = mesh.get_face_inner_product() # Basic inner product matrix (faces)
mesh.set_cell_gradient_BC(["dirichlet", "dirichlet", "dirichlet"])
Gc = mesh.cell_gradient # Centers to faces gradient
# Plot Sparse Representation
fig = plt.figure(figsize=(9, 4))
ax1 = fig.add_subplot(121)
ax1.spy(Gn.T * Me * Gn, markersize=0.5)
ax1.set_title("Gn.T*Me*Gn", pad=5)
ax2 = fig.add_subplot(122)
ax2.spy(Gc.T * Mf * Gc, markersize=0.5)
ax2.set_title("Gc.T*Mf*Gc", pad=5)

Text(0.5, 1.0770202020202022, 'Gc.T*Mf*Gc')
Total running time of the script: (0 minutes 0.448 seconds)