Mixed Galerkin for the Poisson equation with PETSc

Table of Contents

• Remarks on DG vs CG
• Mapping onto Peano's PETSc multigrid extension

A brief, informal discussion how to implement mixed Galerkin methods with PETSc within Peano.

This is a brief discussion of key ideas behind Mixed Galerkin (DG) methods. While we discuss some of the mathematics behind mixed formulations, the text is informal (from a maths point of view) but puts special emphasis on how numerical concepts are mapped onto Peano's data structures and mesh traversal. We focus notably on the combination of mixed methods with Discontinous Galerkin methods. One can use mixed continuous Galerkin methods, but they are not particular widespread due to multiple numerical issues. The script code resulting from this tutorial can be found in benchmarks/multigrid/petsc/mixed.py.

Remarks on DG vs CG

If we use a mixed formulation within globally continuous shape functions, the boundary integrals over \( u \) disappear as they form a telescope sum. The fact that \( u \) is continuous however does not imply that the gradients are continuous, too. Unless we enforce this, we have to use a DG method for the "mixed" part of the equations even though we kicked off with a CG model for the primary unknown.

In general, having a continuous space for both types of quantities implies that the spaces are not compatible with each other: A continuous \( \phi \) space is too restrictive to represent the derivatives. If we try to "approximate" them therein, we will end up with some inconsistent representations.

Therefore DG is a natural choice for mixed formulations. As we need DG for the helper PDEs anyway, it is reasonably to start with DG from the first place. It makes no sense to change your working horse.

Mapping onto Peano's PETSc multigrid extension

Once we commit to the formulation

\begin{eqnarray*} - \text{div} \ \phi & = & f, \\ \nabla u & = & \phi. \end{eqnarray*}

and omit, for the time being, the right-hand side, as see that this system of PDEs translates 1:1 onto two Discontinuous Galerkin solvers:

primary_solver = multigrid.petsc.solvers.DiscontinuousGalerkinDiscretisation(name                 = "Poisson_PrimaryVariable",
                                                                     polynomial_degree    = args.degree,
                                                                     cell_unknowns        = 1,                      # unknowns_per_cell_dof,
                                                                     face_unknowns        = 1,                      # unknowns_per_face_dof,
                                                                     dimensions           = args.dimensions,
                                                                     cell_cell_rhs_matrix   = [...],     
                                                                     cell_cell_lhs_matrix   = [...],
                                                                     cell_cell_lhs_matrix_scaling = args.dimensions-1,
                                                                     cell_cell_rhs_matrix_scaling = args.dimensions,
                                                                     [...]
                                                                     )
gradient_solver = multigrid.petsc.solvers.DiscontinuousGalerkinDiscretisation(name                 = "Poisson_Gradients",
                                                                     polynomial_degree    = args.degree,
                                                                     cell_unknowns        = args.dimensions,        # unknowns_per_cell_dof,
                                                                     face_unknowns        = 1,                      # unknowns_per_face_dof,
                                                                     dimensions           = args.dimensions,
                                                                     cell_cell_rhs_matrix   = [...],     
                                                                     cell_cell_lhs_matrix   = [...],
                                                                     cell_cell_lhs_matrix_scaling = args.dimensions-1,
                                                                     cell_cell_rhs_matrix_scaling = [...],
                                                                     [...]
                                                                     )

The primary solver is a scalar PDE, where we have a derivative on the left-hand side. As we integrate over the domain left and right, we get a scaling of \( h^{d-1} \) on the left-hand and a scaling of \( h^{d} \) on the right-hand side. These scalings are encoded within cell_cell_lhs_matrix_scaling and cell_cell_rhs_matrix_scaling.

The gradient solver represents a system of PDEs with \( d \) unknowns/equations. Again we have one derivative on the left-hand side and therefore a scaling of \( h^{d-1} \). Different to the primary solver, we do not project the "whole" approximation onto the face (which is only a scalar anyway, so the term "whole" makes limited sense"), but we only the normal component. Therefore, this solver uses @icode face_unknowns = 1 @endicode two even though the cells' hold vector unknowns. The other ingredients, i.e. the Riemann solver and the cell-to-face projections all follow 1:1 the discussion for the @ref tutorials_documentation_multigrid_petsc_discontinuous_galerkin_poisson "Discontinuous Galerkin" solver. We have to set them per PDE system, but they are totally analogous.

The coupling between the two solvers is realised through the right-hand side of the gradient_solver. This part requires some work (it is actually the only part that requires additional work compared to the discontinuous Galerkin), as the solvers, so far, know nothing of each other:

gradient_solver.right_hand_size_assembly = xxxxx