 |
Peano
|
|
Peano
|
We solve the Poisson equation
\( -\Delta u = 4 \pi^2 \prod_{d} \sin (2 \pi x_d) \)
over the unit square or cube \( \Omega = [0,1]^d \), subject to Dirichlet boundary conditions \( u |_{\partial \Omega} = 0 \). The exact solution of this problem is \( u_{\mathrm{exact}} = \prod_d \sin (2 \pi x_d) \)
Things that you should be familiar with before you read this tutorial:
This simple example discusses the steps and code snippets from the file benchmarks/multigrid/petsc/poisson/dg.py. We omit some technical details (such as the parsing of the command line) and focus on key ingredients. Therefore, it is reasonable to study this description in combination with the actual script. It is best to read all high-level documentation in combination with the documentation of the class multigrid.petsc.solvers.DiscontinuousGalerkinDiscretisation which provides a lot of details on the actual implementation and mathematical ingredients. This page is more on the how to use the class to construct a working solver.
First, we have to create a multigrid project.
This is kind of an empty hull until we fill it with some global information (such as the domain size) and also tell the API to scan the global Makefile.
We will insert the solvers into this project in a minute. Once this is done, we need to configure the project
A DG solver is created through
We discuss the details of the parameters below, as they depend - to some degree - upon the chosen Riemann solver. For the time being, it is sufficient to recognise that we create a solver object, instruct it about the matrices/operators that are to be used, and the add it to the PETSc project.
Finally, we can ask the PETSc project to build a Peano project which we can translate into C++. So we have a multilevel architecture: A high-level, PETSc-specific view is used to construct the solver. Peano's Python API breaks it down into a Peano representation which is a conceptual/logical representation of the final C++ code. This Peano project then is used to actually write all code snippets and glue code to end up with a genuine C++ application. You can then either compile the code via a plain Makefile through the Python API or you switch to the command line. All of these steps follow the generic architecture of any Peano project.
Once the script has terminated, a lot of subdirectories will be there in your project repository. They contain glue code and usually don't have to be inspected. There are four files which represent your solver: An abstract solver class AbstractDGPoisson (with a header and an implementation) and its implementation DGPoisson (header plus implementation).
You should never edit the abstract class, as it is overwritten every time you invoke the Python script. However, it is worth having a look: This abstract superclass contains all the matrices, constants, ... That is, it picks up the details from the Python script.
The actual implementation is something you can edit and alter. Unless you prescribe a right-hand side (not done in the example discussed), this is the place where you introduce your problem-specific right-hand side. It is also the place where you inject Dirichlet Boundary conditions for example.
Here we demonstrate how a simple weak formulation resulting from the classical Poisson equation can be implemented within the solver interface.
Throughout the tutorial, we follow the notation as used for multigrid.petsc.solvers.DiscontinuousGalerkinDiscretisation and discuss 2D Poisson equation with \( \mathcal{L}(u) = -\Delta u = - \nabla ^2 u = - \partial_x^2 u - \partial_y^2 u \). The extension to higher dimensions is straightforward. Multiplying by a test-function \( v \) and integrating by parts, our weak formulation becomes
\begin{eqnarray*} \int _\Omega \mathcal{L}(u) \cdot v \ d\boldsymbol{x} &=& \int _\Omega ( -\Delta u) \cdot v \ d\boldsymbol{x} \\ &=& \sum_{\text{cells } K} \int _{K} (- \Delta u) \cdot v \ d\boldsymbol{x} \\ &=& \sum_{\text{cells } K} \left( \int _{K} ( \nabla u, \nabla v) \ d\boldsymbol{x} - \int _{\partial K} (\widehat{\nabla u},n) \cdot v \ ds \right) = \sum_{\text{cells } K} \int _{K} f \cdot v \ d\boldsymbol{x} \end{eqnarray*}
The Dirichlet boundary condition \( u |_{\partial \Omega} = 0 \) is enforced weakly. However, we will defer a discussion of how this is implemented to the next section where we introduce the interior penalty method. In the facet integrals \( \int_{\partial K} \) of the above expression, \( n \) denotes the outward normal of the cell \( K \). The term \( ( \nabla u, n) \), which results from the integration by parts, needs to be evaluated on the cell boundary \( \partial K \). For interior faces this is not defined since the solution is allowed to be discontinuous across cells' boundaries in the DG approach. Hence, we need to replace \( \nabla u \) by some \( \widehat{\nabla u} \) as discussed in the following. First, we introduce the vector \( n_F\) on each facet. On boundary facets we set \( n_F = n \) such that \( n_F \) is the outward normal vector. For all interior facets there are two outward normals, one associated with each of the two neighbouring cells, and we need to make a choice. If the facet is oriented such that its normal points in coordinate direction \( j \) then \( n_F\) is the unit vector with 1 at position \( j \) and zeros everywhere else. For example, in two dimensions we have that \( n_F \in \{(1,0)^T, (0,1)^T\} \) on interior facets. Observe that in this case \( n_F \) and \( n \) only differ by a sign: \( n = \pm n_F \) or more explicitly \( n = n_F (n\cdot n_F) \).
In Peano, we need proceed as follows to deal with the term \( ( \nabla u, n) = ( \nabla u, n_F) (n\cdot n_F) \):
\begin{equation} q^\pm = ( \nabla u, n_F)|_{\partial K} \end{equation}
Mathematically, the projections can be constructed through L-2 projection or extrapolation of the cell values. From the perspective of each individual cell, the projections are distributed as shown in the picture in enumeration discussion: left and bottom facets are assigned the "+" projection, and top and right facets receive the "-" projections. As a result, each facet stores both \( q^+ \) and \( q^- \) coming from its two neighboring cells (the case of boundary facets will be discussed separately).\begin{eqnarray*} \sum_{\text{cells } K} \left( \int _{K} ( \nabla u, \nabla v) \ d\boldsymbol{x} - \int_{\partial K} q^f v(n\cdot n_F)\, ds \right) = \sum_{\text{cells } K} \int _{K} f \cdot v \ d\boldsymbol{x} \end{eqnarray*}
Now, we consider the discretization of the problem at the level of an individual cell and its facets, as it gives us the information required for the arguments in multigrid.petsc.solvers.DiscontinuousGalerkinDiscretisation(). We illustrate each components of the solver with an example of a 1st order discretization. As the problem is discretized by choosing appropriate test functions and calculating the integrals, the discussed above weak formulation for an individual cell results in a set of linear equations, which can be written in the matrix form:
\begin{equation} \left( \begin{array}{cc} A_{c \gets c} & A_{c \gets f} \\ \end{array} \right) \left( \begin{array}{c} \bar{u} \\ \bar{q}^f \end{array} \right) = M \bar{f} \end{equation}
where \( M \) is the mass matrix that arises from the right hand side of the weak formulation above. The projections of the normal gradient can then be represented as follows
\begin{equation} \bar{q}^\pm = P_{f \gets c} \bar{u} \end{equation}
The bar notation of \( \bar{u} \), \( \bar{q}^f \), \( \bar{q}^\pm \) here denotes vectors of nodal values of corresponding unknowns for a single representative cell and its faces as illustrated below.
Here, \( \bar{u} = \left( u_0, u_1, u_2, u_3 \right) \) (4 unknowns), \( \bar{q}^f = \left( q^f_0, ... , q^f_7 \right) \) (8 unknowns) and \( \bar{q}^\pm = \left( q^+_0, ..., q^+_3 , q^-_4, ..., q^-_7 \right) \) (8 unknowns). Note that the vector \( \bar{q}^\pm \) is composed only of those projections \( q^+_i \) and \( q^-_i \) to which the values from the given cell are projected. As such, we do not count her, for example, \( q^-_0 \) and \( q^+_4 \), as they will be assigned values when accessed from neighbouring cells. On the right-hand side, \( \bar{f} \) is the vector of nodal values of the given function \( f \).
The ingredients required for a DG solver are the following:
polynomial_degree . In our 2D 1st order example, \( p = 1 \), \( d = 2 \), so each cell has 4 nodes. cell_unknowns = 1If you have vector-valued PDEs, it is worth studying the documentation of the DG solver itself and notably the section "Data storage" therein.
face_unknowns = 1The number of \( q^f \) unknowns per facet is \( (p+1)^{d-1} \), which gives 2 for \( p = 1 \), \( d = 2 \). Additionally, on each facet we store the same number of \( q^- \) and \( q^+ \) projections, which gives extra \( 2(p+1)^{d-1} \) unknowns per facet.
The operator \( A_{c \gets c} \) takes the dofs within a cell and couples them with each other. It arises from the discretization of the term \( \int _{K} ( \nabla u, \nabla v) \ d\boldsymbol{x} \) and represents a \( (p+1)^{d} \times (p+1)^{d} \) matrix, which should be specified in cell_cell_lhs_matrix . In our example, this is a \( 4 \times 4 \) matrix.
Additionally, since the supplied matrix is supposed to be mesh-independent and computed for a reference unit square, one needs to indicate how it should be scaled when mapped onto a particular mesh cell with spacing \( h \). This should be done by specifying a power of \( h \) in cell_cell_lhs_matrix_scaling . It is easy to show that for the term \( \int _{K} ( \nabla u, \nabla v) \ d\boldsymbol{x} \) this scaling factor is \( h^{d-2} \), i.e.
cell_cell_lhs_matrix_scaling = args.dimensions-2
For two dimensions, it is 0.
The operator \( A_{c \gets f} \) takes the \( 2d(p+1)^{d-1} \) values of \( q^f \) on the faces and couples them with \( (p+1)^{d} \) cell values of \( u \). This is responsible for the face integral term \( \int_{\partial K} q^f v \, ds \) in our weak formulation. Thus, it represents a \( (p+1)^{d} \times 2d(p+1)^{d-1} \) matrix, which is passed on as face_to_cell_matrix with scaling power of \( h \) specified in face_to_cell_matrix_scaling .
In the 2D 1st order case, this is a \( 4 \times 8 \) matrix, where the first two columns correspond to the integral over the left facet, the next two columns correspond to the bottom facet, followed by the next two columns for the right facet, and the last two columns for the top. See this page for more details on the enumeration order.
cell_cell_rhs_matrix_scaling = args.dimensions
cell_to_face_matrix with scaling cell_to_face_matrix_scaling . As we compute derivatives, the matrix needs to be scaled by \( 1/h \), so cell_to_face_matrix_scaling = -1For \( p = 1 \) and \( d = 2 \), \( P_{f \leftarrow c} \) is an \( 8 \times 4 \) matrix with the first two lines corresponding to the normal gradient computed for the left facet, with the remaining lines following the same order as earlier discussed for \( A_{c \leftarrow f} \).
face_face_Riemann_problem_matrix that defines the linear relation \( q^f = \mathcal{R} \left( q^- , q^+ \right) = \frac{1}{2} \left( q^- + q^+ \right) \). Unlike the operators \( A_{c \leftarrow f} \) and \( P_{f \leftarrow c} \), which employ unknowns on all the facets of a given cell, this needs to be specified for one representative facet only, assuming that the relation is the same for all interior facets. As such, it requires a \( (p+1)^{d-1} \times 2(p+1)^{d-1} \) matrix. Let us illustrate this with our reference case. Coupling the unknowns on the left facet gives \begin{equation} \left( \begin{array}{c} q^f_0 \\ q^f_1 \end{array} \right) = \left( \begin{array}{cccc} \frac{1}{2} & \frac{1}{2} & 0 & 0 \\ 0 & 0 & \frac{1}{2} & \frac{1}{2} \end{array} \right) \left( \begin{array}{c} q_0^- \\ q_0^+ \\ q_1^- \\ q_1^+ \end{array} \right) \end{equation}
This \( 2 \times 4 \) matrix needs to be specified asface_face_Riemann_problem_matrix in our case. See Section "The simplest Riemann solver: averaging" for a general discussion.Denoting the global dof-vectors vectors \( \bar{u}^{\text{global}} \), \( \bar{q}^{f,\text{global}} \) and \( \bar{q}^{\pm,\text{global}} \), the weak formulation is equivalent to the matrix equation
\begin{equation} \begin{pmatrix} A_{c\leftarrow c}^{\text{global}} & A_{c\leftarrow f}^{\text{global}} & 0 \\ P_{f\leftarrow c}^{\text{global}} & 0 & -\text{Id} \\ 0 & -\text{Id} & R_{f\leftarrow f}^{\text{global}} \end{pmatrix} \begin{pmatrix} \bar{u}^{\text{global}} \\ \bar{q}^{f,\text{global}} \\ \bar{q}^{\pm,\text{global}} \\ \end{pmatrix} = \begin{pmatrix} M^{\text{global}}\bar{f}^{\text{global}} \\ 0 \\ 0 \end{pmatrix} \end{equation}
In this expression \( A_{c\leftarrow c}^{\text{global}} \), \( A_{c\leftarrow f}^{\text{global}} \), \( P_{f\leftarrow c}^{\text{global}} \) and \( M^{\text{global}} \) are the globally assembled versions of the cell-local matrices \( A_{c\leftarrow c} \), \( A_{c\leftarrow f} \), \( P_{f\leftarrow c} \) and \( M \) defined above. The matrix \( R_{f\leftarrow f}^{\text{global}} \) arises from the construction \( q^f = \mathcal{R}(q^-,q^+) \) on the facets.
It should be stressed that the construction of the global matrix (or, in a matrix-free approach, the machinery for applying it to a given vector) is handled by Peano, so the matrix equation is written down in abstract form here for future reference and to highlight its general structure. The user only needs to specify the local building block as explained above.
The above discretisation, which averages the normal gradient of two neighbouring cells to obtain a value on the facets has the serious drawback that it leads to a singular matrix. Mathematically, the reason for this is that the weak form does not have a unique solution: given a solution \( u \), we can always obtain another solution \( u' = u+\delta u \) by adding an arbitrary function \( \delta u \) which is piecewise constant on the interior cells. In order to make the DG discretization of the Poisson equation well-defined, we need to modify the weak formulation by introducing two additional terms, "penalizing" jumps of the solution across the faces. This method is known as the interior penalty method, and we will use the formulation summarized in Bastian, Blatt, Scheichl: Num. Lin. Alg. with Appl., 19(2), pp.367-388 (2012).
First, we need to introduce some notation. It can be shown that the sum of face integrals in the weak form discussed in the previous paragraph, i.e. \( - \sum_{\text{cells } K} \int _{\partial K} (\widehat{\nabla u},n) \cdot v \ ds \), is equivalent to the following, for the interior facets:
\begin{equation} - \sum_{\text{facets } F} \int _{F} \{ (\nabla u, n_F) \} [\![ v ]\!] \ ds \end{equation}
Unlike the summation over the cells \( K \), where each facet of the boundary \( \partial K \) is integrated over twice, in this case, we sum over the facets, so each of them is covered only once. This results in the jump and average terms under the integral:
\begin{equation} [\![ v ]\!] := v^- - v^+, \quad \{ (\nabla u, n_F) \} := \frac{1}{2} \left( (\nabla u, n_F)^- + (\nabla u, n_F)^+ \right), \end{equation}
where the "-"s and "+"s are the values from the two cells sharing the facet, evaluated on this facet. For the vertical facets, the "+" and "-" cells are the right and the left one, respectively, and for the horizontal ones, they correspond to the upper and the lower cells. In our framework, we will use the projection variables for the "+" and "-" values discussed earlier .
Note that the above equivalence of integral terms is consistent with how we defined the averaging flux \( \mathcal{R}\left( q^- , q^+ \right) = \left( q^- + q^+ \right) / 2 \). If we want to use another relation for \( \mathcal{R} \), we will need to redefine the \( \{ \cdot \} \).
In these terms, the full formulation for our problem reads as follows
\begin{multline} \sum_{\text{cells } K} \int _{K} ( \nabla u, \nabla v) \ d\boldsymbol{x} + \sum_{\text{interior facets } F_i} \int _{F_i} \biggl[ - \{ (\nabla u, n_F) \} [\![ v ]\!] + \theta [\![ u ]\!] \{ (\nabla v, n_F) \} + \gamma_F [\![ u ]\!] [\![ v ]\!] \biggr] \ ds + \\ \sum_{\text{boundary facets } F_b} \int _{F_b} \biggl[ - (\nabla u, n_F) \, v + \theta \, u \, (\nabla v, n_F) + \gamma_F u v \biggr] \ ds = \sum_{\text{cells } K} \int _{K} f v \ d\boldsymbol{x}, \end{multline}
where \( \theta \) and \( \gamma_F \) are user-defined constant penalty parameters. We will consider \( \theta = 1 \) and \( \gamma_F = p (p + d - 1) / h = 2/h \), where \( h \) is the grid spacing. See Bastian, Blatt, Scheichl: Num. Lin. Alg. with Appl., 19(2), pp.367-388 (2012) and reference therein for details on the choice of penalty parameters.
To implement this formulation in a solver, we will use the same notation and follow the same steps as for the "naive discretization" discussed earlier.
First, in addition to the gradient projections \( q^\pm \), we introduce the projections of the solution
\begin{equation} u^+ := - u |_{\text{left and bottom facet}}, \quad u^- := u |_{\text{right and top facet}}. \end{equation}
It is convenient to introduce the opposite signs for the "-" and "+" projections, as we did above, because it will provide the required jumps when summed up. We also introduce \( u^f \), which couples with the projections on each facet in the same way as \( q^f \), i.e. \( u^f = \mathcal{R} \left( u^-, u^+ \right) = \left( u^- + u^+ \right) / 2 \).
The discretized weak formulation takes the following form
\begin{equation} \left( \begin{array}{cc} A_{c \gets c} & A^{\text{penalty}}_{c \gets f} \\ \end{array} \right) \left( \begin{array}{c} \bar{u} \\ \bar{w}^f \end{array} \right) = M \bar{f} \end{equation}
We introduced \( \bar{w}^f \) to denote the collection of \( q^f \) and \( u^f \) values as they are intertwined due to the enumeration system used: we first count all the values on the left facet, followed by the bottom, and so on, i.e. \( \bar{w}^f = \left( q_0^f, q_1^f, u_0^f, u_1^f, q_2^f, q_3^f, u_2^f, u_3^f, ... \right)\), which makes 16 unknowns in our example.
The projection equations become
\begin{equation} \bar{w}^\pm = P_{f \gets c} \bar{u}, \end{equation}
where, similarly, the vector \( \bar{w}^\pm \) combines the projection components following the same order: \( \bar{w}^\pm = \left( q_0^+, q_1^+, u_0^+, u_1^+, \quad q_2^+, q_3^+, u_2^+, u_3^+, \quad q_4^-, q_5^-, ... \right) \). Again, here we only include the unknowns that are assigned values when accessed from the given cell: these are the "+" projections for the left and bottom facets, and "-" projections for the right and top. 16 in total.
Now, let us review the components required to build the solver:
polynomial_degree and set cell_unknowns = 1Number of cell unknowns stays the same.
face_unknowns = 2The number of \( q^f \) and \( u^f \) unknowns together per facet is \( 2(p+1)^{d-1} \), which gives 4 for \( p = 1 \), \( d = 2 \). Additionally, on each facet we store the same number of \( q^-, u^- \) and \( q^+, u^+ \) projections, which gives extra \( 4(p+1)^{d-1} \) unknowns per facet.
The operator \( A^{\text{penalty}}_{c \gets f} \) is an extended version of the operator \( A_{c \gets f} \) of the classical weak formulation, which incorporates the penalty terms. It takes the \( 4d(p+1)^{d-1} \) values of \( q^f \) and \( u^f \) on the faces and couples them with \( (p+1)^{d} \) cell values of \( u \). It represents a \( (p+1)^{d} \times 4d(p+1)^{d-1} \) matrix, which is passed on as face_to_cell_matrix .
This operator now has a part that couples \( u \) with \( q^f \) and another part that couples \( u \) with \( u^f \). The first, "classical", component under the boundary integral is responsible for coupling \( u \) and \( q^f \): When computed for a given cell boundary, it gives the terms \( \pm \int q^f v \, ds \).
Both penalty terms couple \( u \) and \( u^f \). They result from the integrals like \( \pm \theta \int u^f ( \nabla v, n) \, ds \pm 2 \gamma_F \int u^f v \, ds \). In the 2D 1st order case, this is a \( 4 \times 16 \) matrix, which is divided into segments as shown in the picture: blue segments result from the standard flux integral, and the green ones result from the summed-up penalty terms.
As one can see, the blue segments needs to be scaled by \( h \). A potential issue with the current approach to scaling in the interface becomes evident when we look at the green segments. Each green entity is a sum of two penalty terms, one of which should be scaled by \( h^0 \) (no scaling) and another one by \( h \). It means that, in the general case, we cannot separate this scaling without passing the information about \( h \) to the matrix generator, which we want to avoid.
In our case, this issue is resolved by choosing the penalty parameter \( \gamma_F \sim 1/h \), which is consistent with the literature and ensures that the entire green segment has the same scaling factor \( h^0 \).
However, we still need different scaling factors for the blue ( \( h^1 \)) and green ( \( h^0 \)) segments, so we cannot do with just one argument face_to_cell_matrix_scaling , and we split the matrix internally for the purpose of scaling.
The operator \( P_{f \leftarrow c} \) needs to be changed to accommodate the projection \( u^\pm \). It is now a \( 4d(p+1)^{d-1} \times (p+1)^{d} \) (i.e. \( 16 \times 4 \)) matrix with the lines arranged according to the order of unknowns in \( \bar{w}^\pm \). This is specified in cell_to_face_matrix .
Similar to \( A^{\text{penalty}}_{c \gets f} \), the projection matrix will be divided into the segments corresponding to \( q^\pm \) and segments corresponding to \( u^\pm \). The gradient-related \( q^\pm \) segments, as previously, should be scaled by \( h^{-1} \), while the \( u^\pm \) part does not need to be scaled ( \( h^0 \)).
face_face_Riemann_problem_matrix to include the relation \( u^f = \left( u^- + u^+ \right) / 2 \). This gives a \( (p+1)^{d-1} \times 4(p+1)^{d-1} \) matrix: \begin{equation} \left( \begin{array}{c} q^f_0 \\ u^f_0 \\ q^f_1 \\ u^f_1 \end{array} \right) = \left( \begin{array}{cccccccc} \frac{1}{2} & 0 & 0 & 0 & \frac{1}{2} & 0 & 0 & 0 \\ 0 & \frac{1}{2} & 0 & 0 & 0 & \frac{1}{2} & 0 & 0 \\ 0 & 0 & \frac{1}{2} & 0 & 0 & 0 & \frac{1}{2} & 0 \\ 0 & 0 & 0 & \frac{1}{2} & 0 & 0 & 0 & \frac{1}{2} \\ \end{array} \right) \left( \begin{array}{c} q_0^- \\ u_0^- \\ q_1^- \\ u_1^- \\ q_0^+ \\ u_0^+ \\ q_1^+ \\ u_1^+ \\ \end{array} \right) \end{equation}
See Section "The simplest Riemann solver: averaging" for a general discussion on the implementation of the face-face coupling.While on the interior faces we construct \( q^f=\mathcal{R}(q^-,q^+)\), \( u^f =\mathcal{R}(u^-,u^+)\) from pairs of values \( q^\pm \), \( u^\pm \) obtained from the two neighbouring cells, the boundary faces need to be treated differently. For this we set \( q^f = q^+ \), \( u^f = u^+ \) for the left and bottom facets and \( q^f = q^- \), \( u^f = u^- \) for the right and top facets.
Since Peano stores pairs of variables \( u^\pm \) and \( q^\pm \) on all faces, this means that the linear system no longer contains equations for the variables \( q^- \), \( u^- \) for the left and bottom facets and \( q^+ \), \( u^+ \) for the right and top facets. This is expressing the fact that these variables do not have any physical meaning and are decoupled from the solution. As a consequence, they can be set to arbitrary values. To prevent the system matrix from becoming singular, we simply insert 1s into the relevant positions of the global matrix. This implies that if the corresponding entries of the right hand side vector are set to zero, the redundant variables will also be zero. This treatment of boundary values is implemented in python/multigrid/petsc/actionsets/ImposeDirichletBoundaryConditionsWithInteriorPenaltyMethod.py
To see understand the treatment of boundary conditions it is instructive to consider the general case where \( u|_{\partial \Omega} = g \) for some arbitrary function \( g \). As described in Bastian, Blatt, Scheichl: Num. Lin. Alg. with Appl., 19(2), pp.367-388 (2012), the boundary condition is enforced weakly by adding the following expression to the right hand of the weak formulation:
\begin{equation} \sum_{\text{boundary facets } F_b} \int _{F_b} \biggl[ \theta \, g \, (\nabla v, n_F) + \gamma_F g v \biggr] \ ds \end{equation}
In our case \( g=0 \), so this term vanishes. It is also worth pointing out that Neumann boundary conditions can be enforced in a similar way, see Bastian, Blatt, Scheichl: Num. Lin. Alg. with Appl., 19(2), pp.367-388 (2012) for details.
Again, we can write down a the matrix equation for the global dof-vectors vectors \( \bar{u}^{\text{global}} \), \( \bar{w}^{f,\text{global}} \) and \( \bar{w}^{\pm,\text{global}} \) which is equivalent to the weak formulation. This matrix equation is
\begin{equation} \begin{pmatrix} A_{c\leftarrow c}^{\text{global}} & A_{c\leftarrow f}^{\text{penalty,global}} & 0 \\ P_{f\leftarrow c}^{\text{global}} & 0 & -\text{Id} \\ 0 & -\text{Id} & R_{f\leftarrow f}^{\text{global}} \end{pmatrix} \begin{pmatrix} \bar{u}^{\text{global}} \\ \bar{w}^{f,\text{global}} \\ \bar{w}^{\pm,\text{global}} \\ \end{pmatrix} = \begin{pmatrix} M^{\text{global}}\bar{f}^{\text{global}} \\ 0 \\ 0 \end{pmatrix} \end{equation}
where construction of the global matrices are is handled by Peano. The matrix \( R_{f\leftarrow f}^{\text{global}} \) takes into account the correct treatment of redundant variables on boundary faces as discussed above.
Tobias. Continue writing.
This section needs to be revisited or moved.
face_unknowns = 1These are not really unknowns in the linear algebra sense. They are helpers.
face_face_Riemann_problem_matrix .The treatment of the face terms differs from Riemann solver to Riemann solver.
This discussion follows the Discontinuous Galerking discussion. The only difference is that we have to extend our averaging such that both solution components, i.e. \( u \) and the normal projection of the \( \phi \) are taken into account:
\begin{eqnarray*} \int_{\partial \text{cell}_{\psi _\phi }} (n \cdot \phi) \psi _\phi dS(x) & \approx & \int_{\partial \text{cell}_{\psi _\phi }} [[ (\phi,n) ]] \psi _\phi dS(x) = \pm \frac{1}{2} \int_{\partial \text{cell}_{\psi _\phi }} \left( \phi ^+ + \phi ^- \right) \psi _\phi dS(x) \\ \int _{\partial \text{cell}_{\psi _u}} u \cdot (n, \psi _u) dS(x) & \approx & \int _{\partial \text{cell}_{\psi _u}} [[u]] \cdot (n, \psi _u) dS(x) = \frac{1}{2} \int _{\partial \text{cell}_{\psi _u}} (u^- + u^+) \cdot ( n, \psi _u ) dS(x) \end{eqnarray*}
The image below illustrates this:
Neither the solution nor the gradient exist at the cell boundaries. We have to average both of them.
There are two ways to implement this scheme:
The first variant is close-to-trivial and certainly a good sanity check. The face-to-face matrix here is empty/zero, but the user has to take the averaging into account within the projection matrix.
The second variant is more sophsticiated, yet allows users later to implement more sophisticated schemes aka Riemann solvers. We have already gone down this route by introducing dedicated "helper" dofs on the face.
\begin{eqnarray*} \begin{pmatrix} \frac{1}{2} & 0 & 0 & 0 & \frac{1}{2} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & \frac{1}{2} & 0 & 0 & 0 & \frac{1}{2} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{1}{2} & 0 & 0 & 0 & \frac{1}{2} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{1}{2} & 0 & 0 & 0 & \frac{1}{2} & 0 & 0 \end{pmatrix} \end{eqnarray*}
The first two rows refer to projected values from the left. They are sole input values (should be brought to the right by the solver in this context) and hence not subject to any additional equations. The third line averages the \( u \) values from left and right, the fourth line addresses the normal projections of \( \phi \).
In the example above, we use d-linear shape functions arranged along the Gauss-Lobatto integration points. So the two nodes carrying the shape function weights are spread out over the unit interval.
1.10.0