The particular concern of this program are the loops of DG methods. These turn out to be especially complex, primarily because for the face terms, we have to distinguish the cases of boundary, regular interior faces and interior faces with hanging nodes, respectively. The MeshWorker::mesh_loop() handles the complexity on iterating over cells and faces and allows specifying "workers" for the different cell and face terms. The integration of face terms itself, including on adaptively refined faces, is done using the FEInterfaceValues class.
The equation
The model problem solved in this example is the linear advection equation
subject to the boundary conditions
on the inflow part of the boundary of the domain. Here, denotes a vector field, the (scalar) solution function, a boundary value function,
the inflow part of the boundary of the domain and denotes the unit outward normal to the boundary . This equation is the conservative version of the advection equation already considered in step-9 of this tutorial.
On each cell , we multiply by a test function from the left and integrate by parts to get:
When summing this expression over all cells , the boundary integral is done over all internal and external faces and as such there are three cases:
outer boundary on the inflow (we replace by given ):
outer boundary on the outflow:
inner faces (integral from two sides turns into jump, we use the upwind velocity):
Here, the jump is defined as , where the superscripts refer to the left ('+') and right ('-') values at the face. The upwind value is defined to be if and otherwise.
As a result, the mesh-dependent weak form reads:
Here, is the set of all active cells of the triangulation and is the set of all active interior faces. This formulation is known as the upwind discontinuous Galerkin method.
In order to implement this bilinear form, we need to compute the cell terms (first sum) using the usual way to achieve integration on a cell, the interface terms (second sum) using FEInterfaceValues, and the boundary terms (the other two terms). The summation of all those is done by MeshWorker::mesh_loop().
The test problem
We solve the advection equation on with representing a circular counterclockwise flow field, and on and on .
We solve on a sequence of meshes by refining the mesh adaptively by estimating the norm of the gradient on each cell. After solving on each mesh, we output the solution in vtk format and compute the norm of the solution. As the exact solution is either 0 or 1, we can measure the magnitude of the overshoot of the numerical solution with this.
The commented program
The first few files have already been covered in previous examples and will thus not be further commented on:
Here the discontinuous finite elements are defined. They are used in the same way as all other finite elements, though – as you have seen in previous tutorial programs – there isn't much user interaction with finite element classes at all: they are passed to DoFHandler and FEValues objects, and that is about it.
#include <deal.II/fe/fe_dgq.h>
This header is needed for FEInterfaceValues to compute integrals on interfaces:
#include <deal.II/fe/fe_interface_values.h>
We are going to use the simplest possible solver, called Richardson iteration, that represents a simple defect correction. This, in combination with a block SSOR preconditioner (defined in precondition_block.h), that uses the special block matrix structure of system matrices arising from DG discretizations.
#include <deal.II/lac/solver_richardson.h>
#include <deal.II/lac/precondition_block.h>
We are going to use gradients as refinement indicator.
Finally, the new include file for using the mesh_loop from the MeshWorker framework
#include <deal.II/meshworker/mesh_loop.h>
Like in all programs, we finish this section by including the needed C++ headers and declaring we want to use objects in the dealii namespace without prefix.
First, we define a class describing the inhomogeneous boundary data. Since only its values are used, we implement value_list(), but leave all other functions of Function undefined.
Given the flow direction, the inflow boundary of the unit square are the right and the lower boundaries. We prescribe discontinuous boundary values 1 and 0 on the x-axis and value 0 on the right boundary. The values of this function on the outflow boundaries will not be used within the DG scheme.
Finally, a function that computes and returns the wind field . As explained in the introduction, we will use a rotational field around the origin in 2d. In 3d, we simply leave the -component unset (i.e., at zero), whereas the function can not be used in 1d in its current implementation:
The following objects are the scratch and copy objects we use in the call to MeshWorker::mesh_loop(). The new object is the FEInterfaceValues object, that works similar to FEValues or FEFaceValues, except that it acts on an interface between two cells and allows us to assemble the interface terms in our weak form.
The next four members represent the linear system to be solved. system_matrix and right_hand_side are generated by assemble_system(), the solution is computed in solve(). The sparsity_pattern is used to determine the location of nonzero elements in system_matrix.
In the function that sets up the usual finite element data structures, we first need to distribute the DoFs.
dof_handler.distribute_dofs(fe);
We start by generating the sparsity pattern. To this end, we first fill an intermediate object of type DynamicSparsityPattern with the couplings appearing in the system. After building the pattern, this object is copied to sparsity_pattern and can be discarded.
Finally, we set up the structure of all components of the linear system.
system_matrix.reinit(sparsity_pattern);
solution.reinit(dof_handler.n_dofs());
right_hand_side.reinit(dof_handler.n_dofs());
}
The assemble_system function
Here we see the major difference to assembling by hand. Instead of writing loops over cells and faces, the logic is contained in the call to MeshWorker::mesh_loop() and we only need to specify what should happen on each cell, each boundary face, and each interior face. These three tasks are handled by the lambda functions inside the function below.
We solve a homogeneous equation, thus no right hand side shows up in the cell term. What's left is integrating the matrix entries.
for (unsignedint point = 0; point < fe_v.n_quadrature_points; ++point)
{
auto beta_q = beta(q_points[point]);
for (unsignedint i = 0; i < n_dofs; ++i)
for (unsignedint j = 0; j < n_dofs; ++j)
{
copy_data.cell_matrix(i, j) +=
-beta_q // -\beta
* fe_v.shape_grad(i, point) // \nabla \phi_i
* fe_v.shape_value(j, point) // \phi_j
* JxW[point]; // dx
}
}
};
This is the function called for boundary faces and consists of a normal integration using FEFaceValues. New is the logic to decide if the term goes into the system matrix (outflow) or the right-hand side (inflow).
This is the function called on interior faces. The arguments specify cells, face and subface indices (for adaptive refinement). We just pass them along to the reinit() function of FEInterfaceValues.
void reinit(const CellIteratorType &cell, const unsigned int face_no, const unsigned int sub_face_no, const CellNeighborIteratorType &cell_neighbor, const unsigned int face_no_neighbor, const unsigned int sub_face_no_neighbor, const unsigned int q_index=numbers::invalid_unsigned_int, const unsigned int mapping_index=numbers::invalid_unsigned_int, const unsigned int fe_index=numbers::invalid_unsigned_int)
The following lambda function will handle copying the data from the cell and face assembly into the global matrix and right-hand side.
While we would not need an AffineConstraints object, because there are no hanging node constraints in DG discretizations, we use an empty object here as this allows us to use its copy_local_to_global functionality.
Here, we finally handle the assembly. We pass in ScratchData and CopyData objects, the lambda functions from above, an specify that we want to assemble interior faces once.
For this simple problem we use the simplest possible solver, called Richardson iteration, that represents a simple defect correction. This, in combination with a block SSOR preconditioner, that uses the special block matrix structure of system matrices arising from DG discretizations. The size of these blocks are the number of DoFs per cell. Here, we use a SSOR preconditioning as we have not renumbered the DoFs according to the flow field. If the DoFs are renumbered in the downstream direction of the flow, then a block Gauss-Seidel preconditioner (see the PreconditionBlockSOR class with relaxation=1) does a much better job.
std::cout << " Solver converged in " << solver_control.last_step()
<< " iterations." << std::endl;
}
We refine the grid according to a very simple refinement criterion, namely an approximation to the gradient of the solution. As here we consider the DG(1) method (i.e. we use piecewise bilinear shape functions) we could simply compute the gradients on each cell. But we do not want to base our refinement indicator on the gradients on each cell only, but want to base them also on jumps of the discontinuous solution function over faces between neighboring cells. The simplest way of doing that is to compute approximative gradients by difference quotients including the cell under consideration and its neighbors. This is done by the DerivativeApproximation class that computes the approximate gradients in a way similar to the GradientEstimation described in step-9 of this tutorial. In fact, the DerivativeApproximation class was developed following the GradientEstimation class of step-9. Relating to the discussion in step-9, here we consider . Furthermore we note that we do not consider approximate second derivatives because solutions to the linear advection equation are in general not in but only in (or, to be more precise: in , i.e., the space of functions whose derivatives in direction are square integrable).
template <int dim>
void AdvectionProblem<dim>::refine_grid()
{
The DerivativeApproximation class computes the gradients to float precision. This is sufficient as they are approximate and serve as refinement indicators only.
The output of this program consists of a vtk file of the adaptively refined grids and the numerical solutions. Finally, we also compute the L-infinity norm of the solution using VectorTools::integrate_difference().
The output of this program consist of the console output and solutions in vtk format:
Cycle 0
Number of active cells: 64
Number of degrees of freedom: 256
Solver converged in 4 iterations.
Writing solution to <solution-0.vtk>
L-infinity norm: 1.09057
Cycle 1
Number of active cells: 112
Number of degrees of freedom: 448
Solver converged in 9 iterations.
Writing solution to <solution-1.vtk>
L-infinity norm: 1.10402
Cycle 2
Number of active cells: 214
Number of degrees of freedom: 856
Solver converged in 16 iterations.
Writing solution to <solution-2.vtk>
L-infinity norm: 1.09813
Cycle 3
Number of active cells: 415
Number of degrees of freedom: 1660
Solver converged in 26 iterations.
Writing solution to <solution-3.vtk>
L-infinity norm: 1.09579
Cycle 4
Number of active cells: 796
Number of degrees of freedom: 3184
Solver converged in 44 iterations.
Writing solution to <solution-4.vtk>
L-infinity norm: 1.09612
Cycle 5
Number of active cells: 1561
Number of degrees of freedom: 6244
Solver converged in 81 iterations.
Writing solution to <solution-5.vtk>
We show the solutions on the initial mesh, the mesh after two and after five adaptive refinement steps.
And finally we show a plot of a 3d computation.
Why use discontinuous elements
In this program we have used discontinuous elements. It is a legitimate question to ask why not simply use the normal, continuous ones. Of course, to everyone with a background in numerical methods, the answer is obvious: the continuous Galerkin (cG) method is not stable for the transport equation, unless one specifically adds stabilization terms. The DG method, however, is stable. Illustrating this with the current program is not very difficult; in fact, only the following minor modifications are necessary:
Add handling of hanging node constraints in exactly the same way as step-6.
We need a different solver; the direct solver in step-29 is a convenient choice. An experienced deal.II user will be able to do this in less than 10 minutes.
While the 2d solution has been shown above, containing a number of small spikes at the interface that are, however, stable in height under mesh refinement, results look much different when using a continuous element:
0
1
2
3
4
5
In refinement iteration 5, the image can't be plotted in a reasonable way any more as a 3d plot. We thus show a color plot with a range of (the solution values of the exact solution lie in , of course). In any case, it is clear that the continuous Galerkin solution exhibits oscillatory behavior that gets worse and worse as the mesh is refined more and more.
There are a number of strategies to stabilize the cG method, if one wants to use continuous elements for some reason. Discussing these methods is beyond the scope of this tutorial program; an interested reader could, for example, take a look at step-31.
Possibilities for extensions
Given that the exact solution is known in this case, one interesting avenue for further extensions would be to confirm the order of convergence for this program. In the current case, the solution is non-smooth, and so we can not expect to get a particularly high order of convergence, even if we used higher order elements. But even if the solution is smooth, the equation is not elliptic and so it is not immediately clear that we should obtain a convergence order that equals that of the optimal interpolation estimates (i.e. for example that we would get convergence in the norm by using quadratic elements).
In fact, for hyperbolic equations, theoretical predictions often indicate that the best one can hope for is an order one half below the interpolation estimate. For example, for the streamline diffusion method (an alternative method to the DG method used here to stabilize the solution of the transport equation), one can prove that for elements of degree , the order of convergence is on arbitrary meshes. While the observed order is frequently on uniformly refined meshes, one can construct so-called Peterson meshes on which the worse theoretical bound is actually attained. This should be relatively simple to verify, for example using the VectorTools::integrate_difference function.
A different direction is to observe that the solution of transport problems often has discontinuities and that therefore a mesh in which we bisect every cell in every coordinate direction may not be optimal. Rather, a better strategy would be to only cut cells in the direction parallel to the discontinuity. This is called anisotropic mesh refinement and is the subject of step-30.