DG Method From Scratch (in Python) - An Introduction
Understanding the overall program architecture, math notation, and data structures

last update: Mar 2, 2025

Article Motivation

  • Outline the structure of this series of articles.
  • Understand the discontinuous Galerkin method.
  • List the mathematical notation used.
  • Show examples of data structures used in the algorithm.

Introduction

Hyperbolic PDEs

  • 1D hyperbolic system
\begin{equation*} \end{equation*} \begin{equation*} \frac{\partial \mathbf{u}}{\partial t} + \mathcal{A} \frac{\partial \mathbf{u}}{\partial x} = 0 \end{equation*}

Where \(\mathcal{A}\) is diagonalizable.

  • 2D hyperbolic system
\begin{equation*} \frac{\partial \mathbf{u}}{\partial t} + \mathcal{A}_x \frac{\partial \mathbf{u}}{\partial x} + \mathcal{A}_y \frac{\partial \mathbf{u}}{\partial y} = 0 \end{equation*}
  • scalar conservation law in 1D
\begin{equation*} \frac{\partial u}{\partial t} + \frac{\partial f(u)}{\partial x} = 0 \end{equation*}
  • 3D conservation law
\begin{equation*} \frac{\partial \mathbf{u}}{\partial t} = \nabla \cdot F(\mathbf{u}) = \frac{\partial F(\mathbf{u})}{\partial x} + \frac{\partial F(\mathbf{u})}{\partial y} + \frac{\partial F(\mathbf{u})}{\partial z} \end{equation*}

Examples of data structures on a simple mesh

For reference, below are some examples of the matrix data structures used in the algorithm developed in this series of articles.

Some fields, like p, the pressure, are defined over all nodes in the simulation, including on the faces of elements, and in the interior of elements.

Other fields, like dp, the jump in pressure across a face, are defined only over face nodes.

All the matrices below are based on the simple five element tetrahedral mesh of a cubic domain.

That means number_of_cells = 5, nodes_per_cell = 10, nodes_per_face = \(\frac{(n+1)(n+2)}{2} =\) 6, faces_per_cell = 4.

cell_to_vertices(number_of_cells, faces_per_cell)

  • cell_to_vertices[i,j] - the global vertex number of the ith vertex of the jth cell
array([[0, 1, 3, 4],
       [1, 2, 3, 6],
       [1, 6, 4, 5],
       [3, 4, 6, 7],
       [1, 6, 3, 4]])

cell_to_cells(number_of_cells, faces_per_cell)

  • cell_to_cell[i,j] - the global cell number of the cell that touches the jth face of cell i.
  • If I have a cell and a face of that cell, then I can find the adjacent cell like this, adjacent_cell = cell_to_cell[cell, face]
array([[0, 0, 4, 0],
       [1, 1, 1, 4],
       [4, 2, 2, 2],
       [4, 3, 3, 3],
       [1, 2, 3, 0]])

cell_to_faces(number_of_cells, faces_per_cell)

  • cell_to_face[i,j] - the local face number of the cell that touches cell i on it's jth face.
  • If I have a cell and a face of that cell, then I can find the corresponding face of the adjacent cell like this, adjacent_face = cell_to_face[cell, face]
  • This can be used together with cell_to_cell to identify which faces of which cells are adjacent.
array([[0, 1, 3, 3],
       [0, 1, 2, 0],
       [1, 1, 2, 3],
       [2, 1, 2, 3],
       [3, 0, 0, 2]])

node_ids(nodes_per_cell, number_of_cells)

  • node_ids[i,j] is the global node id number of the ith node on the jth mesh cell.
  • provides a global reference for the id numbers of every cell on the mesh.
np.array([[ 0, 10, 20, 30, 40],
          [ 1, 11, 21, 31, 41],
          [ 2, 12, 22, 32, 42],
          [ 3, 13, 23, 33, 43],
          [ 4, 14, 24, 34, 44],
          [ 5, 15, 25, 35, 45],
          [ 6, 16, 26, 36, 46],
          [ 7, 17, 27, 37, 47],
          [ 8, 18, 28, 38, 48],
          [ 9, 19, 29, 39, 49]])

interior_face_node_indices(nodes_per_face * number_of_cells)

  • Called vmapM by Hesthaven and warburton [1].
  • Provides the indices for the nodes on the interior side of a face of a particular cell. Can be used to index all interior face nodes to create fields defined only over face nodes.
array([ 0,  3,  5,  6,  8,  9,  0,  1,  2,  6,  7,  9,  2,  4,  5,  7,  8,
        9,  0,  1,  2,  3,  4,  5, 10, 13, 15, 16, 18, 19, 10, 11, 12, 16,
       17, 19, 12, 14, 15, 17, 18, 19, 10, 11, 12, 13, 14, 15, 20, 23, 25,
       26, 28, 29, 20, 21, 22, 26, 27, 29, 22, 24, 25, 27, 28, 29, 20, 21,
       22, 23, 24, 25, 30, 33, 35, 36, 38, 39, 30, 31, 32, 36, 37, 39, 32,
       34, 35, 37, 38, 39, 30, 31, 32, 33, 34, 35, 40, 43, 45, 46, 48, 49,
       40, 41, 42, 46, 47, 49, 42, 44, 45, 47, 48, 49, 40, 41, 42, 43, 44,
       45])

exterior_face_node_indices(nodes_per_face * number_of_cells)

  • Called vmapP by Hesthaven and warburton [1].
  • Provides the indices for the nodes on the exterior side of a face of a particular cell. Can be used to index all exterior face nodes to create fields defined only over exterior face nodes.
array([ 0,  3,  5,  6,  8,  9,  0,  1,  2,  6,  7,  9, 42, 44, 45, 41, 43,
       40,  0,  1,  2,  3,  4,  5, 10, 13, 15, 16, 18, 19, 10, 11, 12, 16,
       17, 19, 12, 14, 15, 17, 18, 19, 40, 46, 49, 43, 48, 45, 40, 41, 42,
       46, 47, 49, 20, 21, 22, 26, 27, 29, 22, 24, 25, 27, 28, 29, 20, 21,
       22, 23, 24, 25, 45, 48, 49, 44, 47, 42, 30, 31, 32, 36, 37, 39, 32,
       34, 35, 37, 38, 39, 30, 31, 32, 33, 34, 35, 10, 13, 15, 11, 14, 12,
       20, 23, 25, 26, 28, 29, 39, 36, 30, 38, 33, 35,  9,  7,  2,  8,  4,
        5])

boundary_face_node_indices(number_of_nodes_on_boundary)

  • Called vmapB by Hesthaven and Warburton [1].
  • Use to index boundary nodes on a field defined only on face nodes (like dp, the jump in pressure across face nodes)
array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 18, 19, 20, 21, 22,
       23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39,
       40, 41, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68,
       69, 70, 71, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91,
       92, 93, 94, 95])

boundary_node_indices(number_of_nodes_on_boundary)

  • Called mapB by Hesthaven and Warburton [1].
  • Use to index boundary nodes on a field defined on all nodes.
array([ 0,  3,  5,  6,  8,  9,  0,  1,  2,  6,  7,  9,  0,  1,  2,  3,  4,
        5, 10, 13, 15, 16, 18, 19, 10, 11, 12, 16, 17, 19, 12, 14, 15, 17,
       18, 19, 20, 21, 22, 26, 27, 29, 22, 24, 25, 27, 28, 29, 20, 21, 22,
       23, 24, 25, 30, 31, 32, 36, 37, 39, 32, 34, 35, 37, 38, 39, 30, 31,
       32, 33, 34, 35])

r_differentiation_marix(number_of_nodes, number_of_nodes)

  • Differentiate with respect to r on the reference tetrahedron.
  • For example, if u is defined over all node (i.e. with size (number_of_nodes, number_of_cells)) and you do matrix multiplication, r_differentiation_matrix @ u, you will get the derivative of u with respect to r.
  • In order to turn this into a derivative in the physical space you must multiply the result by the appropriate metric constant for the mapping from reference space to physical space, \(\frac{\partial r}{\partial x}\), named drdx.

lift_matrix(number_of_nodes, faces_per_cell * nodes_per_face)

  • Operator that mimics the extraction of a surface integral like this:
\begin{equation*} \mathcal{E} [\boldsymbol{g^1},\boldsymbol{g^2},\boldsymbol{g^3},\boldsymbol{g^4}]^T = \int_{\partial \mathsf{I}} \mathbf{\hat{n}} \cdot \boldsymbol{g}_h \ell_i(\boldsymbol{r}) d \boldsymbol{r} \end{equation*}
  • Here \([\boldsymbol{g^1},\boldsymbol{g^2},\boldsymbol{g^3},\boldsymbol{g^4}]^T\) is a matrix of size (points_per_cell, number_of_cells) where each row represents all the nodal and each column represents an element.
  • By multiplying the lift matrix, \(\mathcal{E}\) by a matrix defined by a field over all face nodes, we can compute the surface integral of that term.
  • In the code we'll define the lift matrix as \(\mathcal{M^{-1}}\mathcal{E}\). This just simplifies things because the inverse of the mass matrix always turns up in the surface integral calculations.
  • The lift matrix is formed by first initializing a matrix of the correct size with all zeros. Then finding the mass matrix for each face of the reference tetrahedron, and inserting those elements into the corresponding face nodes of the lift matrix.
  • On page 187 of Hesthaven and Warburton they explain that \(\boldsymbol{g}\) is a polynomial trace function, composed of the numerical flux or the jump in flux, depending on whether we use the weak or the strong form

Conclusion

References

1.
Hesthaven JS, Warburton T. Nodal discontinuous galerkin methods: Algorithms, analysis, and applications. 1st ed. Springer New York, NY; 2008. (Texts in applied mathematics; vol. 54).