This program was contributed by Zhuoran Wang. Some more information about this program, as well as more numerical results, are presented in [Wang2019] .
Introduction
This tutorial program presents an implementation of the "weak Galerkin" finite element method for the Poisson equation. In some sense, the motivation for considering this method starts from the same point as in step-51: We would like to consider discontinuous shape functions, but then need to address the fact that the resulting problem has a much larger number of degrees of freedom compared to the usual continuous Galerkin method (because, for example, each vertex carries as many degrees of freedom as there are adjacent cells). We also have to address the fact that, unlike in the continuous Galerkin method, every degree of freedom on one cell couples with all of the degrees of freedom on each of its face neighbor cells. Consequently, the matrix one gets from the "traditional" discontinuous Galerkin methods are both large and relatively dense.
Both the hybridized discontinuous Galerkin method (HDG) in step-51 and the weak Galerkin (WG) method in this tutorial address the issue of coupling by introducing additional degrees of freedom whose shape functions only live on a face between cells (i.e., on the "skeleton" of the mesh), and which therefore "insulate" the degrees of freedom on the adjacent cells from each other: cell degrees of freedom only couple with other cell degrees of freedom on the same cell, as well as face degrees of freedom, but not with cell degrees of freedom on neighboring cells. Consequently, the coupling of shape functions for these cell degrees of freedom indeed couple on exactly one cell and the degrees of freedom defined on its faces.
For a given equation, say the second order Poisson equation, the difference between the HDG and the WG method is how precisely one formulates the problem that connects all of these different shape functions. (Indeed, for some WG and HDG formulation, it is possible to show that they are equivalent.) The HDG does things by reformulating second order problems in terms of a system of first order equations and then conceptually considers the face degrees of freedom to be "fluxes" of this first order system. In contrast, the WG method keeps things in second order form and considers the face degrees of freedom as of the same type as the primary solution variable, just restricted to the lower-dimensional faces. For the purposes of the equation, one then needs to somehow "extend" these shape functions into the interior of the cell when defining what it means to apply a differential operator to them. Compared to the HDG, the method has the advantage that it does not lead to a proliferation of unknowns due to rewriting the equation as a first-order system, but it is also not quite as easy to implement. However, as we will see in the following, this additional effort is not prohibitive.
Weak Galerkin finite element methods
Weak Galerkin Finite Element Methods (WGFEMs) use discrete weak functions to approximate scalar unknowns, and discrete weak gradients to approximate classical gradients. The method was originally introduced by Junping Wang and Xiu Ye in the paper A weak Galerkin finite element method for second order elliptic problems, J. Comput. Appl. Math., 103-115, 2013. Compared to the continuous Galerkin method, the weak Galerkin method satisfies important physical properties, namely local mass conservation and bulk normal flux continuity. It results in a SPD linear system, and optimal convergence rates can be obtained with mesh refinement.
The equation to solve
This program solves the Poisson equation using the weak Galerkin finite element method:
where is a bounded domain. In the context of the flow of a fluid through a porous medium, is the pressure, is a permeability tensor, is the source term, and represent Dirichlet and Neumann boundary conditions. We can introduce a flux, , that corresponds to the Darcy velocity (in the way we did in step-20) and this variable will be important in the considerations below.
In this program, we will consider a test case where the exact pressure is on the unit square domain, with homogeneous Dirichelet boundary conditions and the identity matrix. Then we will calculate errors of pressure, velocity, and flux.
Weak Galerkin scheme
The Poisson equation above has a solution that needs to satisfy the weak formulation of the problem,
for all test functions , where
and
Here, we have integrated by parts in the bilinear form, and we are evaluating the gradient of in the interior and the values of on the boundary of the domain. All of this is well defined because we assume that the solution is in for which taking the gradient and evaluating boundary values are valid operations.
The idea of the weak Galerkin method is now to approximate the exact solution with a discontinuous function. This function may only be discontinuous along interfaces between cells, and because we will want to evaluate this function also along interfaces, we have to prescribe not only what values it is supposed to have in the cell interiors but also its values along interfaces. We do this by saying that is actually a tuple, , though it's really just a single function that is either equal to or , depending on whether it is evaluated at a point that lies in the cell interior or on cell interfaces.
We would then like to simply stick this approximation into the bilinear form above. This works for the case where we have to evaluate the test function on the boundary (where we would simply take its interface part ) but we have to be careful with the gradient because that is only defined in cell interiors. Consequently, the weak Galerkin scheme for the Poisson equation is defined by
for all discrete test functions , where
and
The key point is that here, we have replaced the gradient by the discrete weak gradient operator that makes sense for our peculiarly defined approximation .
The question is then how that operator works. For this, let us first say how we think of the discrete approximation of the pressure. As mentioned above, the "function" actually consists of two parts: the values in the interior of cells, and on the interfaces. We have to define discrete (finite-dimensional) function spaces for both of these; in this program, we will use FE_DGQ for as the space in the interior of cells (defined on each cell, but in general discontinuous along interfaces), and FE_FaceQ for as the space on the interfaces.
Then let us consider just a single cell (because the integrals above are all defined cell-wise, and because the weak discrete gradient is defined cell-by-cell). The restriction of to cell , then consists of the pair . In essence, we can think of of some function defined on that approximates the gradient; in particular, if was the restriction of a differentiable function (to the interior and boundary of – which would make it continuous between the interior and boundary), then would simply be the exact gradient . But, since is not continuous between interior and boundary of , we need a more general definition; furthermore, we can not deal with arbitrary functions, and so require that is also in a finite element space (which, since the gradient is a vector, has to be vector-valued, and because the weak gradient is defined on each cell separately, will also be discontinuous between cells).
The way this is done is to define this weak gradient operator (where is the vector-valued Raviart-Thomas space of order on cell ) in the following way:
for all test functions . This is, in essence, simply an application of the integration-by-parts formula. In other words, for a given , we need to think of as that Raviart-Thomas function of degree for which the left hand side and right hand side are equal for all test functions.
A key point to make is then the following: While the usual gradient is a local operator that computes derivatives based simply on the value of a function at a point and its (infinitesimal) neighborhood, the weak discrete gradient does not have this property: It depends on the values of the function it is applied to on the entire cell, including the cell's boundary. Both are, however, linear operators as is clear from the definition of above, and that will allow us to represent via a matrix in the discussion below.
Note
It may be worth pointing out that while the weak discrete gradient is an element of the Raviart-Thomas space on each cell , it is discontinuous between cells. On the other hand, the Raviart-Thomas space defined on the entire mesh and implemented by the FE_RaviartThomas class represents functions that have continuous normal components at interfaces between cells. This means that globally, is not in , even though it is on every cell in . Rather, it is in a "broken" Raviart-Thomas space that below we will represent by the symbol . (The term "broken" here refers to the process of "breaking something apart", and not to the synonym to the expression "not functional".) One might therefore (rightfully) argue that the notation used in the weak Galerkin literature is a bit misleading, but as so often it all depends on the context in which a certain notation is used – in the current context, references to the Raviart-Thomas space or element are always understood to be to the "broken" spaces.
deal.II happens to have an implementation of this broken Raviart-Thomas space: The FE_DGRaviartThomas class. As a consequence, in this tutorial we will simply always use the FE_DGRaviartThomas class, even though in all of those places where we have to compute cell-local matrices and vectors, it makes no difference.
Representing the weak gradient
Since is an element of a finite element space, we can expand it in a basis as we always do, i.e., we can write
Here, since has two components (the interior and the interface components), the same must hold true for the basis functions , which we can write as . If you've followed the descriptions in step-8, step-20, and the documentation topic on vector-valued problems, it will be no surprise that for some values of , will be zero, whereas for other values of , will be zero – i.e., shape functions will be of either one or the other kind. That is not important, here, however. What is important is that we need to wonder how we can represent because that is clearly what will appear in the problem when we want to implement the bilinear form
The key point is that is known to be a member of the "broken" Raviart-Thomas space . What this means is that we can represent (on each cell separately)
where the functions , and where is a matrix of dimension
(That the weak discrete gradient can be represented as a matrix should not come as a surprise: It is a linear operator from one finite dimensional space to another finite dimensional space. If one chooses bases for both of these spaces, then every linear operator can of course be written as a matrix mapping the vector of expansion coefficients with regards to the basis of the domain space of the operator, to the vector of expansion coefficients with regards to the basis in the image space.)
Using this expansion, we can easily use the definition of the weak discrete gradient above to define what the matrix is going to be:
for all test functions .
This clearly leads to a linear system of the form
with
and consequently
(In this last step, we have assumed that the indices only range over those degrees of freedom active on cell , thereby ensuring that the mass matrix on the space is invertible.) Equivalently, using the symmetry of the matrix , we have that
Also worth pointing out is that the matrices and are of course not square but rectangular.
Assembling the linear system
Having explained how the weak discrete gradient is defined, we can now come back to the question of how the linear system for the equation in question should be assembled. Specifically, using the definition of the bilinear form shown above, we then need to compute the elements of the local contribution to the global matrix,
As explained above, we can expand in terms of the Raviart-Thomas basis on each cell, and similarly for :
By re-arranging sums, this yields the following expression:
So, if we have the matrix for each cell , then we can easily compute the contribution for cell to the matrix as follows:
Here,
which is really just the mass matrix on cell using the Raviart-Thomas basis and weighting by the permeability tensor . The derivation here then shows that the weak Galerkin method really just requires us to compute these and matrices on each cell , and then , which is easily computed. The code to be shown below does exactly this.
Having so computed the contribution of cell to the global matrix, all we have to do is to "distribute" these local contributions into the global matrix. How this is done is first shown in step-3 and step-4. In the current program, this will be facilitated by calling AffineConstraints::distribute_local_to_global().
A linear system of course also needs a right hand side. There is no difficulty associated with computing the right hand side here other than the fact that we only need to use the cell-interior part for each shape function .
Post-processing and L2-errors
The discussions in the previous sections have given us a linear system that we can solve for the numerical pressure . We can use this to compute an approximation to the variable that corresponds to the velocity with which the medium flows in a porous medium if this is the model we are trying to solve. This kind of step – computing a derived quantity from the solution of the discrete problem – is typically called "post-processing".
Here, instead of using the exact gradient of , let us instead use the discrete weak gradient of to calculate the velocity on each element. As discussed above, on each element the gradient of the numerical pressure can be approximated by discrete weak gradients :
On cell , the numerical velocity can be written as
where is the expansion matrix from above, and is the basis function of the space on a cell.
Unfortunately, may not be in the space (unless, of course, if is constant times the identity matrix). So, in order to represent it in a finite element program, we need to project it back into a finite dimensional space we can work with. Here, we will use the -projection to project it back to the (broken) space.
We define the projection as on each cell . For any , So, rather than the formula shown above, the numerical velocity on cell instead becomes
and we have the following system to solve for the coefficients :
In the implementation below, the matrix with elements is called cell_matrix_D, whereas the matrix with elements is called cell_matrix_E.
Then the elementwise velocity is
where is called cell_velocity in the code.
Using this velocity obtained by "postprocessing" the solution, we can define the -errors of pressure, velocity, and flux by the following formulas:
where is the area of the element, are faces of the element, are unit normal vectors of each face. The last of these norms measures the accuracy of the normal component of the velocity vectors over the interfaces between the cells of the mesh. The scaling factor is chosen so as to scale out the difference in the length (or area) of the collection of interfaces as the mesh size changes.
The first of these errors above is easily computed using VectorTools::integrate_difference. The others require a bit more work and are implemented in the code below.
The commented program
Include files
This program is based on step-7, step-20 and step-51, so most of the following header files are familiar. We need the following, of which only the one that imports the FE_DGRaviartThomas class (namely, deal.II/fe/fe_dg_vector.h) is really new; the FE_DGRaviartThomas implements the "broken" Raviart-Thomas space discussed in the introduction:
#include <deal.II/base/quadrature.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/tensor_function.h>
#include <deal.II/base/function.h>
#include <deal.II/base/point.h>
#include <deal.II/lac/block_vector.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/block_sparse_matrix.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/fe_raviart_thomas.h>
#include <deal.II/fe/fe_dg_vector.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/fe_face.h>
#include <deal.II/fe/component_mask.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/data_out_faces.h>
#include <fstream>
#include <iostream>
Our first step, as always, is to put everything related to this tutorial program into its own namespace:
This is the main class of this program. We will solve for the numerical pressure in the interior and on faces using the weak Galerkin (WG) method, and calculate the error of pressure. In the post-processing step, we will also calculate -errors of the velocity and flux.
The structure of the class is not fundamentally different from that of previous tutorial programs, so there is little need to comment on the details with one exception: The class has a member variable fe_dgrt that corresponds to the "broken" Raviart-Thomas space mentioned in the introduction. There is a matching dof_handler_dgrt that represents a global enumeration of a finite element field created from this element, and a vector darcy_velocity that holds nodal values for this field. We will use these three variables after solving for the pressure to compute a postprocessed velocity field for which we can then evaluate the error and which we can output for visualization.
Right hand side, boundary values, and exact solution
Next, we define the coefficient matrix (here, the identity matrix), Dirichlet boundary conditions, the right-hand side , and the exact solution that corresponds to these choices for and , namely .
The class that implements the exact pressure solution has an oddity in that we implement it as a vector-valued one with two components. (We say that it has two components in the constructor where we call the constructor of the base Function class.) In the value() function, we do not test for the value of the component argument, which implies that we return the same value for both components of the vector-valued function. We do this because we describe the finite element in use in this program as a vector-valued system that contains the interior and the interface pressures, and when we compute errors, we will want to use the same pressure solution to test both of these components.
In this constructor, we create a finite element space for vector valued functions, which will here include the ones used for the interior and interface pressures, and .
After we have created the mesh above, we distribute degrees of freedom and resize matrices and vectors. The only piece of interest in this function is how we interpolate the boundary values for the pressure. Since the pressure consists of interior and interface components, we need to make sure that we only interpolate onto that component of the vector-valued solution space that corresponds to the interface pressures (as these are the only ones that are defined on the boundary of the domain). We do this via a component mask object for only the interface pressures.
template <int dim>
void WGDarcyEquation<dim>::setup_system()
{
dof_handler.distribute_dofs(fe);
dof_handler_dgrt.distribute_dofs(fe_dgrt);
std::cout << " Number of pressure degrees of freedom: "
In the bilinear form, there is no integration term over faces between two neighboring cells, so we can just use DoFTools::make_sparsity_pattern to calculate the sparse matrix.
This function is more interesting. As detailed in the introduction, the assembly of the linear system requires us to evaluate the weak gradient of the shape functions, which is an element in the Raviart-Thomas space. As a consequence, we need to define a Raviart-Thomas finite element object, and have FEValues objects that evaluate it at quadrature points. We then need to compute the matrix on every cell , for which we need the matrices and mentioned in the introduction.
A point that may not be obvious is that in all previous tutorial programs, we have always called FEValues::reinit() with a cell iterator from a DoFHandler. This is so that one can call functions such as FEValuesBase::get_function_values() that extract the values of a finite element function (represented by a vector of DoF values) on the quadrature points of a cell. For this operation to work, one needs to know which vector elements correspond to the degrees of freedom on a given cell – i.e., exactly the kind of information and operation provided by the DoFHandler class.
We could create a DoFHandler object for the "broken" Raviart-Thomas space (using the FE_DGRaviartThomas class), but we really don't want to here: At least in the current function, we have no need for any globally defined degrees of freedom associated with this broken space, but really only need to reference the shape functions of such a space on the current cell. As a consequence, we use the fact that one can call FEValues::reinit() also with cell iterators into Triangulation objects (rather than DoFHandler objects). In this case, FEValues can of course only provide us with information that only references information about cells, rather than degrees of freedom enumerated on these cells. So we can't use FEValuesBase::get_function_values(), but we can use FEValues::shape_value() to obtain the values of shape functions at quadrature points on the current cell. It is this kind of functionality we will make use of below. The variable that will give us this information about the Raviart-Thomas functions below is then the fe_values_rt (and corresponding fe_face_values_rt) object.
Given this introduction, the following declarations should be pretty obvious:
This finally gets us in position to loop over all cells. On each cell, we will first calculate the various cell matrices used to construct the local matrix – as they depend on the cell in question, they need to be re-computed on each cell. We need shape functions for the Raviart-Thomas space as well, for which we need to create first an iterator to the cell of the triangulation, which we can obtain by assignment from the cell pointing into the DoFHandler.
for (constauto &cell : dof_handler.active_cell_iterators())
The first cell matrix we will compute is the mass matrix for the Raviart-Thomas space. Hence, we need to loop over all the quadrature points for the velocity FEValues object.
cell_matrix_M = 0;
for (unsignedint q = 0; q < n_q_points_dgrt; ++q)
for (unsignedint i = 0; i < dofs_per_cell_dgrt; ++i)
Next we take the inverse of this matrix by using FullMatrix::gauss_jordan(). It will be used to calculate the coefficient matrix later. It is worth recalling later that cell_matrix_M actually contains the inverse of after this call.
cell_matrix_M.gauss_jordan();
From the introduction, we know that the right hand side of the equation that defines is the difference between a face integral and a cell integral. Here, we approximate the negative of the contribution in the interior. Each component of this matrix is the integral of a product between a basis function of the polynomial space and the divergence of a basis function of the Raviart-Thomas space. These basis functions are defined in the interior.
cell_matrix_G = 0;
for (unsignedint q = 0; q < n_q_points; ++q)
for (unsignedint i = 0; i < dofs_per_cell_dgrt; ++i)
{
constdouble div_v_i =
fe_values_dgrt[velocities].divergence(i, q);
for (unsignedint j = 0; j < dofs_per_cell; ++j)
{
constdouble phi_j_interior =
fe_values[pressure_interior].value(j, q);
cell_matrix_G(i, j) -=
(div_v_i * phi_j_interior * fe_values.JxW(q));
}
}
Next, we approximate the integral on faces by quadrature. Each component is the integral of a product between a basis function of the polynomial space and the dot product of a basis function of the Raviart-Thomas space and the normal vector. So we loop over all the faces of the element and obtain the normal vector.
Finally we can compute the local matrix . Element is given by . We have calculated the coefficients in the previous step, and so obtain the following after suitably re-arranging the loops:
local_matrix = 0;
for (unsignedint q = 0; q < n_q_points_dgrt; ++q)
{
for (unsignedint k = 0; k < dofs_per_cell_dgrt; ++k)
The last step is to distribute components of the local matrix into the system matrix and transfer components of the cell right hand side into the system right hand side:
In this function, compute the velocity field from the pressure solution previously computed. The velocity is defined as , which requires us to compute many of the same terms as in the assembly of the system matrix. There are also the matrices we need to assemble (see the introduction) but they really just follow the same kind of pattern.
Computing the same matrices here as we have already done in the assemble_system() function is of course wasteful in terms of CPU time. Likewise, we copy some of the code from there to this function, and this is also generally a poor idea. A better implementation might provide for a function that encapsulates this duplicated code. One could also think of using the classic trade-off between computing efficiency and memory efficiency to only compute the matrices once per cell during the assembly, storing them somewhere on the side, and re-using them here. (This is what step-51 does, for example, where the assemble_system() function takes an argument that determines whether the local matrices are recomputed, and a similar approach – maybe with storing local matrices elsewhere – could be adapted for the current program.)
In the introduction, we explained how to calculate the numerical velocity on the cell. We need the pressure solution values on each cell, coefficients of the Gram matrix and coefficients of the projection. We have already calculated the global solution, so we will extract the cell solution from the global solution. The coefficients of the Gram matrix have been calculated when we assembled the system matrix for the pressures. We will do the same way here. For the coefficients of the projection, we do matrix multiplication, i.e., the inverse of the Gram matrix times the matrix with as components. Then, we multiply all these coefficients and call them beta. The numerical velocity is the product of beta and the basis functions of the Raviart-Thomas space.
Then we also need, again, to compute the matrix that is used to evaluate the weak discrete gradient. This is the exact same code as used in the assembly of the system matrix, so we just copy it from there:
cell_matrix_G = 0;
for (unsignedint q = 0; q < n_q_points; ++q)
for (unsignedint i = 0; i < dofs_per_cell_dgrt; ++i)
This part is to calculate the error of the pressure. We define a vector that holds the norm of the error on each cell. Next, we use VectorTool::integrate_difference() to compute the error in the norm on each cell. However, we really only care about the error in the interior component of the solution vector (we can't even evaluate the interface pressures at the quadrature points because these are all located in the interior of cells) and consequently have to use a weight function that ensures that the interface component of the solution variable is ignored. This is done by using the ComponentSelectFunction whose arguments indicate which component we want to select (component zero, i.e., the interior pressures) and how many components there are in total (two).
In this function, we evaluate errors for the velocity on each cell, and errors for the flux on faces. The function relies on the compute_postprocessed_velocity() function having previous computed, which computes the velocity field based on the pressure solution that has previously been computed.
We are going to evaluate velocities on each cell and calculate the difference between numerical and exact velocities.
Having previously computed the postprocessed velocity, we here only have to extract the corresponding values on each cell and face and compare it to the exact values.
for (constauto &cell_dgrt : dof_handler_dgrt.active_cell_iterators())
{
fe_values_dgrt.reinit(cell_dgrt);
First compute the error between the postprocessed velocity field and the exact one:
For reconstructing the flux we need the size of cells and faces. Since fluxes are calculated on faces, we have the loop over all four faces of each cell. To calculate the face velocity, we extract values at the quadrature points from the darcy_velocity which we have computed previously. Then, we calculate the squared velocity error in normal direction. Finally, we calculate the flux error on the cell by appropriately scaling with face and cell areas and add it to the global error.
constdouble cell_area = cell_dgrt->measure();
for (constauto &face_dgrt : cell_dgrt->face_iterators())
We have two sets of results to output: the interior solution and the skeleton solution. We use DataOut to visualize interior results. The graphical output for the skeleton results is done by using the DataOutFaces class.
In both of the output files, both the interior and the face variables are stored. For the interface output, the output file simply contains the interpolation of the interior pressures onto the faces, but because it is undefined which of the two interior pressure variables you get from the two adjacent cells, it is best to ignore the interior pressure in the interface output file. Conversely, for the cell interior output file, it is of course impossible to show any interface pressures , because these are only available on interfaces and not cell interiors. Consequently, you will see them shown as an invalid value (such as an infinity).
For the cell interior output, we also want to output the velocity variables. This is a bit tricky since it lives on the same mesh but uses a different DoFHandler object (the pressure variables live on the dof_handler object, the Darcy velocity on the dof_handler_dgrt object). Fortunately, there are variations of the DataOut::add_data_vector() function that allow specifying which DoFHandler a vector corresponds to, and consequently we can visualize the data from both DoFHandler objects within the same file.
We run the program with a right hand side that will produce the solution and with homogeneous Dirichlet boundary conditions in the domain . In addition, we choose the coefficient matrix in the differential operator as the identity matrix. We test this setup using , and element combinations, which one can select by using the appropriate constructor argument for the WGDarcyEquation object in main(). We will then visualize pressure values in interiors of cells and on faces. We want to see that the pressure maximum is around 1 and the minimum is around 0. With mesh refinement, the convergence rates of pressure, velocity and flux should then be around 1 for , 2 for , and 3 for .
Test results on WG(Q0,Q0;RT[0])
The following figures show interior pressures and face pressures using the element. The mesh is refined 2 times (top) and 4 times (bottom), respectively. (This number can be adjusted in the make_grid() function.) When the mesh is coarse, one can see the face pressures neatly between the values of the interior pressures on the two adjacent cells.
From the figures, we can see that with the mesh refinement, the maximum and minimum pressure values are approaching the values we expect. Since the mesh is a rectangular mesh and numbers of cells in each direction is even, we have symmetric solutions. From the 3d figures on the right, we can see that on , the pressure is a constant in the interior of the cell, as expected.
Convergence table for k=0
We run the code with differently refined meshes (chosen in the make_grid() function) and get the following convergence rates of pressure, velocity, and flux (as defined in the introduction).
number of refinements
2
1.587e-01
5.113e-01
7.062e-01
3
8.000e-02
2.529e-01
3.554e-01
4
4.006e-02
1.260e-01
1.780e-01
5
2.004e-02
6.297e-02
8.902e-02
Conv.rate
1.00
1.00
1.00
We can see that the convergence rates of are around 1. This, of course, matches our theoretical expectations.
Test results on WG(Q1,Q1;RT[1])
We can repeat the experiment from above using the next higher polynomial degree: The following figures are interior pressures and face pressures implemented using . The mesh is refined 4 times. Compared to the previous figures using , on each cell, the solution is no longer constant on each cell, as we now use bilinear polynomials to do the approximation. Consequently, there are 4 pressure values in one interior, 2 pressure values on each face.
Compared to the corresponding image for the combination, the solution is now substantially more accurate and, in particular so close to being continuous at the interfaces that we can no longer distinguish the interface pressures from the interior pressures on the adjacent cells.
Convergence table for k=1
The following are the convergence rates of pressure, velocity, and flux we obtain from using the element combination:
number of refinements
2
1.613e-02
5.093e-02
7.167e-02
3
4.056e-03
1.276e-02
1.802e-02
4
1.015e-03
3.191e-03
4.512e-03
5
2.540e-04
7.979e-04
1.128e-03
Conv.rate
2.00
2.00
2.00
The convergence rates of are around 2, as expected.
Test results on WG(Q2,Q2;RT[2])
Let us go one polynomial degree higher. The following are interior pressures and face pressures implemented using , with mesh size (i.e., 5 global mesh refinement steps). In the program, we use data_out_face.build_patches(fe.degree) when generating graphical output (see the documentation of DataOut::build_patches()), which here implies that we divide each 2d cell interior into 4 subcells in order to provide a better visualization of the quadratic polynomials.
Convergence table for k=2
As before, we can generate convergence data for the errors of pressure, velocity, and flux using the combination:
number of refinements
2
1.072e-03
3.375e-03
4.762e-03
3
1.347e-04
4.233e-04
5.982e-04
4
1.685e-05
5.295e-05
7.487e-05
5
2.107e-06
6.620e-06
9.362e-06
Conv.rate
3.00
3.00
3.00
Once more, the convergence rates of is as expected, with values around 3.