 |
Reference documentation for deal.II version 9.2.0
|
|
Reference documentation for deal.II version 9.2.0
|
This tutorial depends on step-16, step-7, step-39.
| Table of contents | |
|---|---|
This is a variant of step-16 with the only change that we are using the MeshWorker framework with the pre-made LocalIntegrator helper classes instead of assembling the face terms using FEInterfaceValues.
The details of this framework on how it is used in practice will be explained as part of this tutorial program.
The problem we solve here is the same as the one in step-12.
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.
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.
We are going to use gradients as refinement indicator.
Here come the new include files for using the MeshWorker framework. The first contains the class MeshWorker::DoFInfo, which provides local integrators with a mapping between local and global degrees of freedom. It stores the results of local integrals as well in its base class MeshWorker::LocalResults. In the second of these files, we find an object of type MeshWorker::IntegrationInfo, which is mostly a wrapper around a group of FEValues objects. The file meshworker/simple.h contains classes assembling locally integrated data into a global system containing only a single matrix. Finally, we will need the file that runs the loop over all mesh cells and faces.
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 \([0,1]^2\) 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 \(\beta=\beta(\mathbf x)\). As explained in the introduction, we will use a rotational field around the origin in 2d. In 3d, we simply leave the \(z\)-component unset (i.e., at zero), whereas the function can not be used in 1d in its current implementation:
After this preparations, we proceed with the main class of this program, called AdvectionProblem. It is basically the main class of step-6. We do not have an AffineConstraints object, because there are no hanging node constraints in DG discretizations.
Major differences will only come up in the implementation of the assemble functions, since here, we not only need to cover the flux integrals over faces, we also use the MeshWorker interface to simplify the loops involved.
Furthermore we want to use DG elements of degree 1 (but this is only specified in the constructor). If you want to use a DG method of a different degree the whole program stays the same, only replace 1 in the constructor by the desired polynomial degree.
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.
Finally, we have to provide functions that assemble the cell, boundary, and inner face terms. Within the MeshWorker framework, the loop over all cells and much of the setup of operations will be done outside this class, so all we have to provide are these three operations. They will then work on intermediate objects for which first, we here define alias to the info objects handed to the local integration functions in order to make our life easier below.
The following three functions are then the ones that get called inside the generic loop over all cells and faces. They are the ones doing the actual integration.
In our code below, these functions do not access member variables of the current class, so we can mark them as static and simply pass pointers to these functions to the MeshWorker framework. If, however, these functions would want to access member variables (or needed additional arguments beyond the ones specified below), we could use the facilities of lambda functions to provide the MeshWorker framework with objects that act as if they had the required number and types of arguments, but have in fact other arguments already bound.
We start with the constructor. The 1 in the constructor call of fe is the polynomial degree.
In the function that sets up the usual finite element data structures, we first need to distribute the DoFs.
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.
To build the sparsity pattern for DG discretizations, we can call the function analogue to DoFTools::make_sparsity_pattern, which is called DoFTools::make_flux_sparsity_pattern:
Finally, we set up the structure of all components of the linear system.
Here we see the major difference to assembling by hand. Instead of writing loops over cells and faces, we leave all this to the MeshWorker framework. In order to do so, we just have to define local integration functions and use one of the classes in namespace MeshWorker::Assembler to build the global system.
This is the magic object, which knows everything about the data structures and local integration. This is the object doing the work in the function MeshWorker::loop(), which is implicitly called by MeshWorker::integration_loop() below. After the functions to which we provide pointers did the local integration, the MeshWorker::Assembler::SystemSimple object distributes these into the global sparse matrix and the right hand side vector.
First, we initialize the quadrature formulae and the update flags in the worker base class. For quadrature, we play safe and use a QGauss formula with number of points one higher than the polynomial degree used. Since the quadratures for cells, boundary and interior faces can be selected independently, we have to hand over this value three times.
These are the types of values we need for integrating our system. They are added to the flags used on cells, boundary and interior faces, as well as interior neighbor faces, which is forced by the four true values.
After preparing all data in info_box, we initialize the FEValues objects in there.
The object created so far helps us do the local integration on each cell and face. Now, we need an object which receives the integrated (local) data and forwards them to the assembler.
Now, we have to create the assembler object and tell it, where to put the local data. These will be our system matrix and the right hand side.
Finally, the integration loop over all active cells (determined by the first argument, which is an active iterator).
As noted in the discussion when declaring the local integration functions in the class declaration, the arguments expected by the assembling integrator class are not actually function pointers. Rather, they are objects that can be called like functions with a certain number of arguments. Consequently, we could also pass objects with appropriate operator() implementations here, or lambda functions if the local integrators were, for example, non-static member functions.
These are the functions given to the MeshWorker::integration_loop() called just above. They compute the local contributions to the system matrix and right hand side on cells and faces.
First, let us retrieve some of the objects used here from info. Note that these objects can handle much more complex structures, thus the access here looks more complicated than might seem necessary.
With these objects, we continue local integration like always. First, we loop over the quadrature points and compute the advection vector in the current point.
We solve a homogeneous equation, thus no right hand side shows up in the cell term. What's left is integrating the matrix entries.
Now the same for the boundary terms. Note that now we use FEValuesBase, the base class for both FEFaceValues and FESubfaceValues, in order to get access to normal vectors.
Finally, the interior face terms. The difference here is that we receive two info objects, one for each cell adjacent to the face and we assemble four matrices, one for each cell and two for coupling back and forth.
For quadrature points, weights, etc., we use the FEValuesBase object of the first argument.
For additional shape functions, we have to ask the neighbors FEValuesBase.
Then we get references to the four local matrices. The letters u and v refer to trial and test functions, respectively. The numbers indicate the cells provided by info1 and info2. By convention, the two matrices in each info object refer to the test functions on the respective cell. The first matrix contains the interior couplings of that cell, while the second contains the couplings between cells.
Here, following the previous functions, we would have the local right hand side vectors. Fortunately, the interface terms only involve the solution and the right hand side does not receive any contributions.
This term we've already seen:
We additionally assemble the term \((\beta\cdot n u,\hat v)_{\partial \kappa_+}\),
This one we've already seen, too:
And this is another new one: \((\beta\cdot n \hat u,\hat v)_{\partial \kappa_-}\):
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.
Here we create the preconditioner,
then assign the matrix to it and set the right block size:
After these preparations we are ready to start the linear solver.
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 \(h^{1+d/2}|\nabla_h u_h|\). Furthermore we note that we do not consider approximate second derivatives because solutions to the linear advection equation are in general not in \(H^2\) but only in \(H^1\) (or, to be more precise: in \(H^1_\beta\), i.e., the space of functions whose derivatives in direction \(\beta\) are square integrable).
The DerivativeApproximation class computes the gradients to float precision. This is sufficient as they are approximate and serve as refinement indicators only.
Now the approximate gradients are computed
and they are cell-wise scaled by the factor \(h^{1+d/2}\)
Finally they serve as refinement indicator.
The output of this program consists of eps-files of the adaptively refined grids and the numerical solutions given in gnuplot format.
First write the grid in eps format.
Then output the solution in gnuplot format.
The following run function is similar to previous examples.
The following main function is similar to previous examples as well, and need not be commented on.
The output of this program is very similar to step-16 and we are not repeating the output here.
We show the solutions on the initial mesh, the mesh after two and after five adaptive refinement steps.
Then we show the final grid (after 5 refinement steps) and the solution again, this time with a nicer 3d rendering (obtained using the DataOutBase::write_vtk function and the VTK-based VisIt visualization program) that better shows the sharpness of the jump on the refined mesh and the over- and undershoots of the solution along the interface:
And finally we show a plot of a 3d computation.
For ideas for further extensions, please see see step-12.
1.8.17