Matrix-free mixed DG for the Poisson equation

This simple example discusses the most basic steps behind a matrix-free mixed Discontinuous Galerkin solver for the Poisson equation.

The standard "benchmark" for this solver (once again solving the Poisson equation on the unit square) can be found in tests/multigrid/matrixfree/poisson/mixed.py. Our discussion follows this script, though we omit many technical details such as the parsing of the command line) and focus on key numerical ingredients from a user perspective. Therefore, it is reasonable to study the present text in combination with the actual script. This page is more on the how to use the class to construct a working solver.

Project setup

Creating a project is a one-liner:

project = mghype.matrixfree.api.Project( project_name = "DGPointwise", 
                              namespace = [ "tests", "multigrid", "matrixfree", "poisson" ]
                            )

Create the solver

The solver construction also is a "one-liner" using the following observations:

• We solve a scalar equation system plus a vector-valued PDE in one rush. Therefore, \( K_c=d+1 \) (unknowns_per_cell_node).
• Within each cell, we have to project the solution onto the face plus one of the \( p \) values. Therefore, projections_per_face_node=2.
• We solve the numerical flux on the face to obtain both a \( \hat u \) value and a \( \widehat{(n,p)} \) value. We set solutions_per_face_node=2.

solver = mghype.matrixfree.solvers.api.DiscontinuousGalerkinDiscretisationPointWiseRiemannSolver(
  name = "DGPoisson",
  dimensions = args.dimensions,
  poly_degree = args.poly_degree,
  unknowns_per_cell_node = 1,
  solutions_per_face_node = 2,
  projections_per_face_node = 2,
  min_h = args.meshsize,
  max_h = args.meshsize,
  assembly_matrix= assembly_matrix, 
  assembly_matrix_scaling = assembly_matrix_scaling, 
  mass_matrix = mass_matrix, 
  mass_matrix_scaling = mass_matrix_scaling, 
  face_from_cell_projection = face_from_cell_projection, 
  face_from_cell_projection_scaling = face_from_cell_projection_scaling,
  cell_from_face_projection = cell_from_face_projection, 
  cell_from_face_projection_scaling = cell_from_face_projection_scaling,
  riemann_matrix = matrices.get_face_face_riemann_problem_matrix(),
  boundary_matrix = matrices.get_boundary_matrix(),
  cell_relaxation = args.omega_c,
  face_relaxation = args.omega_f,  # default here is 1, i.e. we solve exactly. See class docu
  approximate_system_matrix = approximate_system_matrix,
  approximate_system_matrix_scaling = approximate_system_matrix_scaling,
  solver_tolerance = 0.01,
)
 

Matrices

We have \( (p+1)^d \) Terminology and degree of freedom layout nodes in the cell. Consequently, the matrix-matrix part of the global assembly is from \( \mathbb{R}^{ (d+1)(p+1)^d \times (d+1)(p+1)^d } \). Conceptionally, this matrix resembles

\( \begin{bmatrix} X & X & X \\ 0 & X & 0 \\ 0 & 0 & X \end{bmatrix} \begin{bmatrix} u \\ p_x \\ p_y \end{bmatrix} = \left( h^d \begin{bmatrix} 0 & X & X \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} + h^{d-1} \begin{bmatrix} X & 0 & 0 \\ 0 & X & 0 \\ 0 & 0 & X \end{bmatrix} \right) \begin{bmatrix} u \\ p_x \\ p_y \end{bmatrix} = \begin{bmatrix} 0 \\ f_x \\ f_y \end{bmatrix} \)

but we store the entries as AoS according to our storage format. The sparsity pattern thus is subject to some permutations. The structure discussion above emphasises that the first line of the equation system holds entries with different scaling. The \( p \) entries enter with a \( h^d \) scaling, while the \( u \) terms have a scaling of \( h^{d-1} \). The PDE over the secondary variables is structurally simpler: All the entries on the left-hand side have a scaling of \( h^{d-1} \).

 
matrices = mghype.api.matrixgenerators.GLMatrixFree(
  args.dimensions,
  args.poly_degree,
  1,  # Unknowns per cell dof. Scalar PDE here
  2,  # We use the penalty formulation (see docu in tutorials)
  2
)
 
assembly_matrix, assembly_matrix_scaling = matrices.get_cell_system_matrix_for_laplacian()
mass_matrix, mass_matrix_scaling         = matrices.get_cell_mass_matrix()
face_from_cell_projection, \
face_from_cell_projection_scaling        = matrices.get_face_from_cell_matrix()
cell_from_face_projection, \
cell_from_face_projection_scaling        = matrices.get_cell_from_face_matrix()
approximate_system_matrix, \
approximate_system_matrix_scaling        = matrices.get_A_tilde()

Todo

Tobias and Alex should discuss.

The Riemann problem (numerical flux)

Todo

Tobias, please write.