step-59.h

Go to the documentation of this file.

 
true,
 *                  MatrixFree<dim, number>::DataAccessOnFaces::gradients,
 *                  MatrixFree<dim, number>::DataAccessOnFaces::gradients);
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * Since the Laplacian is symmetric, the `Tvmult()` (needed by the multigrid
 * smoother interfaces) operation is simply forwarded to the `vmult()` case.
 * 
 
 * 
 * 
 * @code
 *     template <int dim, int fe_degree, typename number>
 *     void LaplaceOperator<dim, fe_degree, number>::Tvmult(
 *       LinearAlgebra::distributed::Vector<number>       &dst,
 *       const LinearAlgebra::distributed::Vector<number> &src) const
 *     {
 *       vmult(dst, src);
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * The cell operation is very similar to @ref step_37 "step-37". We do not use a
 * coefficient here, though. The second difference is that we replaced the
 * two steps of FEEvaluation::read_dof_values() followed by
 * FEEvaluation::evaluate() by a single function call
 * FEEvaluation::gather_evaluate() which internally calls the sequence of
 * the two individual methods. Likewise, FEEvaluation::integrate_scatter()
 * implements the sequence of FEEvaluation::integrate() followed by
 * FEEvaluation::distribute_local_to_global(). In this case, these new
 * functions merely save two lines of code. However, we use them for the
 * analogy with FEFaceEvaluation where they are more important as
 * explained below.
 * 
 
 * 
 * 
 * @code
 *     template <int dim, int fe_degree, typename number>
 *     void LaplaceOperator<dim, fe_degree, number>::apply_cell(
 *       const MatrixFree<dim, number>                    &data,
 *       LinearAlgebra::distributed::Vector<number>       &dst,
 *       const LinearAlgebra::distributed::Vector<number> &src,
 *       const std::pair<unsigned int, unsigned int>      &cell_range) const
 *     {
 *       FEEvaluation<dim, fe_degree, fe_degree + 1, 1, number> phi(data);
 *       for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
 *         {
 *           phi.reinit(cell);
 *           phi.gather_evaluate(src, EvaluationFlags::gradients);
 *           for (const unsigned int q : phi.quadrature_point_indices())
 *             phi.submit_gradient(phi.get_gradient(q), q);
 *           phi.integrate_scatter(EvaluationFlags::gradients, dst);
 *         }
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * The face operation implements the terms of the interior penalty method in
 * analogy to @ref step_39 "step-39", as explained in the introduction. We need two
 * evaluator objects for this task, one for handling the solution that comes
 * from the cell on one of the two sides of an interior face, and one for
 * handling the solution from the other side. The evaluators for face
 * integrals are called FEFaceEvaluation and take a boolean argument in the
 * second slot of the constructor to indicate which of the two sides the
 * evaluator should belong two. In FEFaceEvaluation and MatrixFree, we call
 * one of the two sides the `interior` one and the other the `exterior`
 * one. The name `exterior` refers to the fact that the evaluator from both
 * sides will return the same normal vector. For the `interior` side, the
 * normal vector points outwards, whereas it points inwards on the other
 * side, and is opposed to the outer normal vector of that cell. Apart from
 * the new class name, we again get a range of items to work with in
 * analogy to what was discussed in @ref step_37 "step-37", but for the interior faces in
 * this case. Note that the data structure of MatrixFree forms batches of
 * faces that are analogous to the batches of cells for the cell
 * integrals. All faces within a batch involve different cell numbers but
 * have the face number within the reference cell, have the same refinement
 * configuration (no refinement or the same subface), and the same
 * orientation, to keep SIMD operations simple and efficient.
 *   
 
 * 
 * Note that there is no implied meaning in interior versus exterior except
 * the logic decision of the orientation of the normal, which is pretty
 * random internally. One can in no way rely on a certain pattern of
 * assigning interior versus exterior flags, as the decision is made for the
 * sake of access regularity and uniformity in the MatrixFree setup
 * routines. Since most sane DG methods are conservative, i.e., fluxes look
 * the same from both sides of an interface, the mathematics are unaltered
 * if the interior/exterior flags are switched and normal vectors get the
 * opposite sign.
 * 
 
 * 
 * 
 * @code
 *     template <int dim, int fe_degree, typename number>
 *     void LaplaceOperator<dim, fe_degree, number>::apply_face(
 *       const MatrixFree<dim, number>                    &data,
 *       LinearAlgebra::distributed::Vector<number>       &dst,
 *       const LinearAlgebra::distributed::Vector<number> &src,
 *       const std::pair<unsigned int, unsigned int>      &face_range) const
 *     {
 *       FEFaceEvaluation<dim, fe_degree, fe_degree + 1, 1, number> phi_inner(data,
 *                                                                            true);
 *       FEFaceEvaluation<dim, fe_degree, fe_degree + 1, 1, number> phi_outer(data,
 *                                                                            false);
 *       for (unsigned int face = face_range.first; face < face_range.second; ++face)
 *         {
 * @endcode
 * 
 * On a given batch of faces, we first update the pointers to the
 * current face and then access the vector. As mentioned above, we
 * combine the vector access with the evaluation. In the case of face
 * integrals, the data access into the vector can be reduced for the
 * special case of an FE_DGQHermite basis as explained for the data
 * exchange above: Since only @f$2(k+1)^{d-1}@f$ out of the @f$(k+1)^d@f$ cell
 * degrees of freedom get multiplied by a non-zero value or derivative
 * of a shape function, this structure can be utilized for the
 * evaluation, significantly reducing the data access. The reduction
 * of the data access is not only beneficial because it reduces the
 * data in flight and thus helps caching, but also because the data
 * access to faces is often more irregular than for cell integrals when
 * gathering values from cells that are farther apart in the index
 * list of cells.
 * 
 * @code
 *           phi_inner.reinit(face);
 *           phi_inner.gather_evaluate(src,
 *                                     EvaluationFlags::values |
 *                                       EvaluationFlags::gradients);
 *           phi_outer.reinit(face);
 *           phi_outer.gather_evaluate(src,
 *                                     EvaluationFlags::values |
 *                                       EvaluationFlags::gradients);
 *   
 * @endcode
 * 
 * The next two statements compute the penalty parameter for the
 * interior penalty method. As explained in the introduction, we would
 * like to have a scaling like @f$\frac{1}{h_\text{i}}@f$ of the length
 * @f$h_\text{i}@f$ normal to the face. For a general non-Cartesian mesh,
 * this length must be computed by the product of the inverse Jacobian
 * times the normal vector in real coordinates. From this vector of
 * `dim` components, we must finally pick the component that is
 * oriented normal to the reference cell. In the geometry data stored
 * in MatrixFree, a permutation of the components in the Jacobian is
 * applied such that this latter direction is always the last
 * component `dim-1` (this is beneficial because reference-cell
 * derivative sorting can be made agnostic of the direction of the
 * face). This means that we can simply access the last entry `dim-1`
 * and must not look up the local face number in
 * `data.get_face_info(face).interior_face_no` and
 * `data.get_face_info(face).exterior_face_no`. Finally, we must also
 * take the absolute value of these factors as the normal could point
 * into either positive or negative direction.
 * 
 * @code
 *           const VectorizedArray<number> inverse_length_normal_to_face =
 *             0.5 * (std::abs((phi_inner.normal_vector(0) *
 *                              phi_inner.inverse_jacobian(0))[dim - 1]) +
 *                    std::abs((phi_outer.normal_vector(0) *
 *                              phi_outer.inverse_jacobian(0))[dim - 1]));
 *           const VectorizedArray<number> sigma =
 *             inverse_length_normal_to_face * get_penalty_factor();
 *   
 * @endcode
 * 
 * In the loop over the quadrature points, we eventually compute all
 * contributions to the interior penalty scheme. According to the
 * formulas in the introduction, the value of the test function gets
 * multiplied by the difference of the jump in the solution times the
 * penalty parameter and the average of the normal derivative in real
 * space. Since the two evaluators for interior and exterior sides get
 * different signs due to the jump, we pass the result with a
 * different sign here. The normal derivative of the test function
 * gets multiplied by the negative jump in the solution between the
 * interior and exterior side. This term, coined adjoint consistency
 * term, must also include the factor of @f$\frac{1}{2}@f$ in the code in
 * accordance with its relation to the primal consistency term that
 * gets the factor of one half due to the average in the test function
 * slot.
 * 
 * @code
 *           for (const unsigned int q : phi_inner.quadrature_point_indices())
 *             {
 *               const VectorizedArray<number> solution_jump =
 *                 (phi_inner.get_value(q) - phi_outer.get_value(q));
 *               const VectorizedArray<number> average_normal_derivative =
 *                 (phi_inner.get_normal_derivative(q) +
 *                  phi_outer.get_normal_derivative(q)) *
 *                 number(0.5);
 *               const VectorizedArray<number> test_by_value =
 *                 solution_jump * sigma - average_normal_derivative;
 *   
 *               phi_inner.submit_value(test_by_value, q);
 *               phi_outer.submit_value(-test_by_value, q);
 *   
 *               phi_inner.submit_normal_derivative(-solution_jump * number(0.5), q);
 *               phi_outer.submit_normal_derivative(-solution_jump * number(0.5), q);
 *             }
 *   
 * @endcode
 * 
 * Once we are done with the loop over quadrature points, we can do
 * the sum factorization operations for the integration loops on faces
 * and sum the results into the result vector, using the
 * `integrate_scatter` function. The name `scatter` reflects the
 * distribution of the vector data into scattered positions in the
 * vector using the same pattern as in `gather_evaluate`. Like before,
 * the combined integrate + write operation allows us to reduce the
 * data access.
 * 
 * @code
 *           phi_inner.integrate_scatter(EvaluationFlags::values |
 *                                         EvaluationFlags::gradients,
 *                                       dst);
 *           phi_outer.integrate_scatter(EvaluationFlags::values |
 *                                         EvaluationFlags::gradients,
 *                                       dst);
 *         }
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * The boundary face function follows by and large the interior face
 * function. The only difference is the fact that we do not have a separate
 * FEFaceEvaluation object that provides us with exterior values @f$u^+@f$, but
 * we must define them from the boundary conditions and interior values
 * @f$u^-@f$. As explained in the introduction, we use @f$u^+ = -u^- + 2
 * g_\text{D}@f$ and @f$\mathbf{n}^-\cdot \nabla u^+ = \mathbf{n}^-\cdot \nabla
 * u^-@f$ on Dirichlet boundaries and @f$u^+=u^-@f$ and @f$\mathbf{n}^-\cdot \nabla
 * u^+ = -\mathbf{n}^-\cdot \nabla u^- + 2 g_\text{N}@f$ on Neumann
 * boundaries. Since this operation implements the homogeneous part, i.e.,
 * the matrix-vector product, we must neglect the boundary functions
 * @f$g_\text{D}@f$ and @f$g_\text{N}@f$ here, and added them to the right hand side
 * in `LaplaceProblem::compute_rhs()`. Note that due to extension of the
 * solution @f$u^-@f$ to the exterior via @f$u^+@f$, we can keep all factors @f$0.5@f$
 * the same as in the inner face function, see also the discussion in
 * @ref step_39 "step-39".
 *   
 
 * 
 * There is one catch at this point: The implementation below uses a boolean
 * variable `is_dirichlet` to switch between the Dirichlet and the Neumann
 * cases. However, we solve a problem where we also want to impose periodic
 * boundary conditions on some boundaries, namely along those in the @f$x@f$
 * direction. One might wonder how those conditions should be handled
 * here. The answer is that MatrixFree automatically treats periodic
 * boundaries as what they are technically, namely an inner face where the
 * solution values of two adjacent cells meet and must be treated by proper
 * numerical fluxes. Thus, all the faces on the periodic boundaries will
 * appear in the `apply_face()` function and not in this one.
 * 
 
 * 
 * 
 * @code
 *     template <int dim, int fe_degree, typename number>
 *     void LaplaceOperator<dim, fe_degree, number>::apply_boundary(
 *       const MatrixFree<dim, number>                    &data,
 *       LinearAlgebra::distributed::Vector<number>       &dst,
 *       const LinearAlgebra::distributed::Vector<number> &src,
 *       const std::pair<unsigned int, unsigned int>      &face_range) const
 *     {
 *       FEFaceEvaluation<dim, fe_degree, fe_degree + 1, 1, number> phi_inner(data,
 *                                                                            true);
 *       for (unsigned int face = face_range.first; face < face_range.second; ++face)
 *         {
 *           phi_inner.reinit(face);
 *           phi_inner.gather_evaluate(src,
 *                                     EvaluationFlags::values |
 *                                       EvaluationFlags::gradients);
 *   
 *           const VectorizedArray<number> inverse_length_normal_to_face = std::abs((
 *             phi_inner.normal_vector(0) * phi_inner.inverse_jacobian(0))[dim - 1]);
 *           const VectorizedArray<number> sigma =
 *             inverse_length_normal_to_face * get_penalty_factor();
 *   
 *           const bool is_dirichlet = (data.get_boundary_id(face) == 0);
 *   
 *           for (const unsigned int q : phi_inner.quadrature_point_indices())
 *             {
 *               const VectorizedArray<number> u_inner = phi_inner.get_value(q);
 *               const VectorizedArray<number> u_outer =
 *                 is_dirichlet ? -u_inner : u_inner;
 *               const VectorizedArray<number> normal_derivative_inner =
 *                 phi_inner.get_normal_derivative(q);
 *               const VectorizedArray<number> normal_derivative_outer =
 *                 is_dirichlet ? normal_derivative_inner : -normal_derivative_inner;
 *               const VectorizedArray<number> solution_jump = (u_inner - u_outer);
 *               const VectorizedArray<number> average_normal_derivative =
 *                 (normal_derivative_inner + normal_derivative_outer) * number(0.5);
 *               const VectorizedArray<number> test_by_value =
 *                 solution_jump * sigma - average_normal_derivative;
 *               phi_inner.submit_normal_derivative(-solution_jump * number(0.5), q);
 *               phi_inner.submit_value(test_by_value, q);
 *             }
 *           phi_inner.integrate_scatter(EvaluationFlags::values |
 *                                         EvaluationFlags::gradients,
 *                                       dst);
 *         }
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * Next we turn to the preconditioner initialization. As explained in the
 * introduction, we want to construct an (approximate) inverse of the cell
 * matrices from a product of 1d mass and Laplace matrices. Our first task
 * is to compute the 1d matrices, which we do by first creating a 1d finite
 * element. Instead of anticipating FE_DGQHermite<1> here, we get the finite
 * element's name from DoFHandler, replace the @p dim argument (2 or 3) by 1
 * to create a 1d name, and construct the 1d element by using FETools.
 * 
 
 * 
 * 
 * @code
 *     template <int dim, int fe_degree, typename number>
 *     void PreconditionBlockJacobi<dim, fe_degree, number>::initialize(
 *       const LaplaceOperator<dim, fe_degree, number> &op)
 *     {
 *       data = op.get_matrix_free();
 *   
 *       std::string name = data->get_dof_handler().get_fe().get_name();
 *       name.replace(name.find('<') + 1, 1, "1");
 *       std::unique_ptr<FiniteElement<1>> fe_1d = FETools::get_fe_by_name<1>(name);
 *   
 * @endcode
 * 
 * As for computing the 1d matrices on the unit element, we simply write
 * down what a typical assembly procedure over rows and columns of the
 * matrix as well as the quadrature points would do. We select the same
 * Laplace matrices once and for all using the coefficients 0.5 for
 * interior faces (but possibly scaled differently in different directions
 * as a result of the mesh). Thus, we make a slight mistake at the
 * Dirichlet boundary (where the correct factor would be 1 for the
 * derivative terms and 2 for the penalty term, see @ref step_39 "step-39") or at the
 * Neumann boundary where the factor should be zero. Since we only use
 * this class as a smoother inside a multigrid scheme, this error is not
 * going to have any significant effect and merely affects smoothing
 * quality.
 * 
 * @code
 *       const unsigned int                                 N = fe_degree + 1;
 *       FullMatrix<double>                                 laplace_unscaled(N, N);
 *       std::array<Table<2, VectorizedArray<number>>, dim> mass_matrices;
 *       std::array<Table<2, VectorizedArray<number>>, dim> laplace_matrices;
 *       for (unsigned int d = 0; d < dim; ++d)
 *         {
 *           mass_matrices[d].reinit(N, N);
 *           laplace_matrices[d].reinit(N, N);
 *         }
 *   
 *       const QGauss<1> quadrature(N);
 *       for (unsigned int i = 0; i < N; ++i)
 *         for (unsigned int j = 0; j < N; ++j)
 *           {
 *             double sum_mass = 0, sum_laplace = 0;
 *             for (unsigned int q = 0; q < quadrature.size(); ++q)
 *               {
 *                 sum_mass += (fe_1d->shape_value(i, quadrature.point(q)) *
 *                              fe_1d->shape_value(j, quadrature.point(q))) *
 *                             quadrature.weight(q);
 *                 sum_laplace += (fe_1d->shape_grad(i, quadrature.point(q))[0] *
 *                                 fe_1d->shape_grad(j, quadrature.point(q))[0]) *
 *                                quadrature.weight(q);
 *               }
 *             for (unsigned int d = 0; d < dim; ++d)
 *               mass_matrices[d](i, j) = sum_mass;
 *   
 * @endcode
 * 
 * The left and right boundary terms assembled by the next two
 * statements appear to have somewhat arbitrary signs, but those are
 * correct as can be verified by looking at @ref step_39 "step-39" and inserting
 * the value -1 and 1 for the normal vector in the 1d case.
 * 
 * @code
 *             sum_laplace +=
 *               (1. * fe_1d->shape_value(i, Point<1>()) *
 *                  fe_1d->shape_value(j, Point<1>()) * op.get_penalty_factor() +
 *                0.5 * fe_1d->shape_grad(i, Point<1>())[0] *
 *                  fe_1d->shape_value(j, Point<1>()) +
 *                0.5 * fe_1d->shape_grad(j, Point<1>())[0] *
 *                  fe_1d->shape_value(i, Point<1>()));
 *   
 *             sum_laplace +=
 *               (1. * fe_1d->shape_value(i, Point<1>(1.0)) *
 *                  fe_1d->shape_value(j, Point<1>(1.0)) * op.get_penalty_factor() -
 *                0.5 * fe_1d->shape_grad(i, Point<1>(1.0))[0] *
 *                  fe_1d->shape_value(j, Point<1>(1.0)) -
 *                0.5 * fe_1d->shape_grad(j, Point<1>(1.0))[0] *
 *                  fe_1d->shape_value(i, Point<1>(1.0)));
 *   
 *             laplace_unscaled(i, j) = sum_laplace;
 *           }
 *   
 * @endcode
 * 
 * Next, we go through the cells and pass the scaled matrices to
 * TensorProductMatrixSymmetricSum to actually compute the generalized
 * eigenvalue problem for representing the inverse: Since the matrix
 * approximation is constructed as @f$A\otimes M + M\otimes A@f$ and the
 * weights are constant for each element, we can apply all weights on the
 * Laplace matrix and simply keep the mass matrices unscaled. In the loop
 * over cells, we want to make use of the geometry compression provided by
 * the MatrixFree class and check if the current geometry is the same as
 * on the last cell batch, in which case there is nothing to do. This
 * compression can be accessed by
 * FEEvaluation::get_mapping_data_index_offset() once `reinit()` has been
 * called.
 *     
 
 * 
 * Once we have accessed the inverse Jacobian through the FEEvaluation
 * access function (we take the one for the zeroth quadrature point as
 * they should be the same on all quadrature points for a Cartesian cell),
 * we check that it is diagonal and then extract the determinant of the
 * original Jacobian, i.e., the inverse of the determinant of the inverse
 * Jacobian, and set the weight as @f$\text{det}(J) / h_d^2@f$ according to
 * the 1d Laplacian times @f$d-1@f$ copies of the @ref GlossMassMatrix "mass matrix".
 * 
 * @code
 *       cell_matrices.clear();
 *       FEEvaluation<dim, fe_degree, fe_degree + 1, 1, number> phi(*data);
 *       unsigned int old_mapping_data_index = numbers::invalid_unsigned_int;
 *       for (unsigned int cell = 0; cell < data->n_cell_batches(); ++cell)
 *         {
 *           phi.reinit(cell);
 *   
 *           if (phi.get_mapping_data_index_offset() == old_mapping_data_index)
 *             continue;
 *   
 *           Tensor<2, dim, VectorizedArray<number>> inverse_jacobian =
 *             phi.inverse_jacobian(0);
 *   
 *           for (unsigned int d = 0; d < dim; ++d)
 *             for (unsigned int e = 0; e < dim; ++e)
 *               if (d != e)
 *                 for (unsigned int v = 0; v < VectorizedArray<number>::size(); ++v)
 *                   AssertThrow(inverse_jacobian[d][e][v] == 0.,
 *                               ExcNotImplemented());
 *   
 *           VectorizedArray<number> jacobian_determinant = inverse_jacobian[0][0];
 *           for (unsigned int e = 1; e < dim; ++e)
 *             jacobian_determinant *= inverse_jacobian[e][e];
 *           jacobian_determinant = 1. / jacobian_determinant;
 *   
 *           for (unsigned int d = 0; d < dim; ++d)
 *             {
 *               const VectorizedArray<number> scaling_factor =
 *                 inverse_jacobian[d][d] * inverse_jacobian[d][d] *
 *                 jacobian_determinant;
 *   
 * @endcode
 * 
 * Once we know the factor by which we should scale the Laplace
 * matrix, we apply this weight to the unscaled DG Laplace matrix
 * and send the array to the class TensorProductMatrixSymmetricSum
 * for computing the generalized eigenvalue problem mentioned in
 * the introduction.
 * 
 
 * 
 * 
 * @code
 *               for (unsigned int i = 0; i < N; ++i)
 *                 for (unsigned int j = 0; j < N; ++j)
 *                   laplace_matrices[d](i, j) =
 *                     scaling_factor * laplace_unscaled(i, j);
 *             }
 *           if (cell_matrices.size() <= phi.get_mapping_data_index_offset())
 *             cell_matrices.resize(phi.get_mapping_data_index_offset() + 1);
 *           cell_matrices[phi.get_mapping_data_index_offset()].reinit(
 *             mass_matrices, laplace_matrices);
 *         }
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * The vmult function for the approximate block-Jacobi preconditioner is
 * very simple in the DG context: We simply need to read the values of the
 * current cell batch, apply the inverse for the given entry in the array of
 * tensor product matrix, and write the result back. In this loop, we
 * overwrite the content in `dst` rather than first setting the entries to
 * zero. This is legitimate for a DG method because every cell has
 * independent degrees of freedom. Furthermore, we manually write out the
 * loop over all cell batches, rather than going through
 * MatrixFree::cell_loop(). We do this because we know that we are not going
 * to need data exchange over the MPI network here as all computations are
 * done on the cells held locally on each processor.
 * 
 
 * 
 * 
 * @code
 *     template <int dim, int fe_degree, typename number>
 *     void PreconditionBlockJacobi<dim, fe_degree, number>::vmult(
 *       LinearAlgebra::distributed::Vector<number>       &dst,
 *       const LinearAlgebra::distributed::Vector<number> &src) const
 *     {
 *       adjust_ghost_range_if_necessary(*data, dst);
 *       adjust_ghost_range_if_necessary(*data, src);
 *   
 *       FEEvaluation<dim, fe_degree, fe_degree + 1, 1, number> phi(*data);
 *       for (unsigned int cell = 0; cell < data->n_cell_batches(); ++cell)
 *         {
 *           phi.reinit(cell);
 *           phi.read_dof_values(src);
 *           cell_matrices[phi.get_mapping_data_index_offset()].apply_inverse(
 *             ArrayView<VectorizedArray<number>>(phi.begin_dof_values(),
 *                                                phi.dofs_per_cell),
 *             ArrayView<const VectorizedArray<number>>(phi.begin_dof_values(),
 *                                                      phi.dofs_per_cell));
 *           phi.set_dof_values(dst);
 *         }
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * The definition of the LaplaceProblem class is very similar to
 * @ref step_37 "step-37". One difference is the fact that we add the element degree as a
 * template argument to the class, which would allow us to more easily
 * include more than one degree in the same program by creating different
 * instances in the `main()` function. The second difference is the
 * selection of the element, FE_DGQHermite, which is specialized for this
 * kind of equations.
 * 
 
 * 
 * 
 * @code
 *     template <int dim, int fe_degree>
 *     class LaplaceProblem
 *     {
 *     public:
 *       LaplaceProblem();
 *       void run();
 *   
 *     private:
 *       void setup_system();
 *       void compute_rhs();
 *       void solve();
 *       void analyze_results() const;
 *   
 *   #ifdef DEAL_II_WITH_P4EST
 *       parallel::distributed::Triangulation<dim> triangulation;
 *   #else
 *       Triangulation<dim> triangulation;
 *   #endif
 *   
 *       const FE_DGQHermite<dim> fe;
 *       DoFHandler<dim>          dof_handler;
 *   
 *       MappingQ1<dim> mapping;
 *   
 *       using SystemMatrixType = LaplaceOperator<dim, fe_degree, double>;
 *       SystemMatrixType system_matrix;
 *   
 *       using LevelMatrixType = LaplaceOperator<dim, fe_degree, float>;
 *       MGLevelObject<LevelMatrixType> mg_matrices;
 *   
 *       LinearAlgebra::distributed::Vector<double> solution;
 *       LinearAlgebra::distributed::Vector<double> system_rhs;
 *   
 *       double             setup_time;
 *       ConditionalOStream pcout;
 *       ConditionalOStream time_details;
 *     };
 *   
 *   
 *   
 *     template <int dim, int fe_degree>
 *     LaplaceProblem<dim, fe_degree>::LaplaceProblem()
 *   #ifdef DEAL_II_WITH_P4EST
 *       : triangulation(MPI_COMM_WORLD,
 *                       Triangulation<dim>::limit_level_difference_at_vertices,
 *                       parallel::distributed::Triangulation<
 *                         dim>::construct_multigrid_hierarchy)
 *   #else
 *       : triangulation(Triangulation<dim>::limit_level_difference_at_vertices)
 *   #endif
 *       , fe(fe_degree)
 *       , dof_handler(triangulation)
 *       , setup_time(0.)
 *       , pcout(std::cout, Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
 *       , time_details(std::cout,
 *                      false &&
 *                        Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
 *     {}
 *   
 *   
 *   
 * @endcode
 * 
 * The setup function differs in two aspects from @ref step_37 "step-37". The first is that
 * we do not need to interpolate any constraints for the discontinuous
 * ansatz space, and simply pass a dummy AffineConstraints object into
 * Matrixfree::reinit(). The second change arises because we need to tell
 * MatrixFree to also initialize the data structures for faces. We do this
 * by setting update flags for the inner and boundary faces,
 * respectively. On the boundary faces, we need both the function values,
 * their gradients, JxW values (for integration), the normal vectors, and
 * quadrature points (for the evaluation of the boundary conditions),
 * whereas we only need shape function values, gradients, JxW values, and
 * normal vectors for interior faces. The face data structures in MatrixFree
 * are always built as soon as one of `mapping_update_flags_inner_faces` or
 * `mapping_update_flags_boundary_faces` are different from the default
 * value `update_default` of UpdateFlags.
 * 
 
 * 
 * 
 * @code
 *     template <int dim, int fe_degree>
 *     void LaplaceProblem<dim, fe_degree>::setup_system()
 *     {
 *       Timer time;
 *       setup_time = 0;
 *   
 *       system_matrix.clear();
 *       mg_matrices.clear_elements();
 *   
 *       dof_handler.distribute_dofs(fe);
 *       dof_handler.distribute_mg_dofs();
 *   
 *       pcout << "Number of degrees of freedom: " << dof_handler.n_dofs()
 *             << std::endl;
 *   
 *       setup_time += time.wall_time();
 *       time_details << "Distribute DoFs               " << time.wall_time() << " s"
 *                    << std::endl;
 *       time.restart();
 *   
 *       AffineConstraints<double> dummy;
 *       dummy.close();
 *   
 *       {
 *         typename MatrixFree<dim, double>::AdditionalData additional_data;
 *         additional_data.tasks_parallel_scheme =
 *           MatrixFree<dim, double>::AdditionalData::none;
 *         additional_data.mapping_update_flags =
 *           (update_gradients | update_JxW_values | update_quadrature_points);
 *         additional_data.mapping_update_flags_inner_faces =
 *           (update_gradients | update_JxW_values | update_normal_vectors);
 *         additional_data.mapping_update_flags_boundary_faces =
 *           (update_gradients | update_JxW_values | update_normal_vectors |
 *            update_quadrature_points);
 *         const auto system_mf_storage =
 *           std::make_shared<MatrixFree<dim, double>>();
 *         system_mf_storage->reinit(
 *           mapping, dof_handler, dummy, QGauss<1>(fe.degree + 1), additional_data);
 *         system_matrix.initialize(system_mf_storage);
 *       }
 *   
 *       system_matrix.initialize_dof_vector(solution);
 *       system_matrix.initialize_dof_vector(system_rhs);
 *   
 *       setup_time += time.wall_time();
 *       time_details << "Setup matrix-free system      " << time.wall_time() << " s"
 *                    << std::endl;
 *       time.restart();
 *   
 *       const unsigned int nlevels = triangulation.n_global_levels();
 *       mg_matrices.resize(0, nlevels - 1);
 *   
 *       for (unsigned int level = 0; level < nlevels; ++level)
 *         {
 *           typename MatrixFree<dim, float>::AdditionalData additional_data;
 *           additional_data.tasks_parallel_scheme =
 *             MatrixFree<dim, float>::AdditionalData::none;
 *           additional_data.mapping_update_flags =
 *             (update_gradients | update_JxW_values);
 *           additional_data.mapping_update_flags_inner_faces =
 *             (update_gradients | update_JxW_values);
 *           additional_data.mapping_update_flags_boundary_faces =
 *             (update_gradients | update_JxW_values);
 *           additional_data.mg_level = level;
 *           const auto mg_mf_storage_level =
 *             std::make_shared<MatrixFree<dim, float>>();
 *           mg_mf_storage_level->reinit(mapping,
 *                                       dof_handler,
 *                                       dummy,
 *                                       QGauss<1>(fe.degree + 1),
 *                                       additional_data);
 *   
 *           mg_matrices[level].initialize(mg_mf_storage_level);
 *         }
 *       setup_time += time.wall_time();
 *       time_details << "Setup matrix-free levels      " << time.wall_time() << " s"
 *                    << std::endl;
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * The computation of the right hand side is a bit more complicated than in
 * @ref step_37 "step-37". The cell term now consists of the negative Laplacian of the
 * analytical solution, `RightHandSide`, for which we need to first split up
 * the Point of VectorizedArray fields, i.e., a batch of points, into a
 * single point by evaluating all lanes in the VectorizedArray
 * separately. Remember that the number of lanes depends on the hardware; it
 * could be 1 for systems that do not offer vectorization (or where deal.II
 * does not have intrinsics), but it could also be 8 or 16 on AVX-512 of
 * recent Intel architectures.
 * 
 * @code
 *     template <int dim, int fe_degree>
 *     void LaplaceProblem<dim, fe_degree>::compute_rhs()
 *     {
 *       Timer time;
 *       system_rhs                          = 0;
 *       const MatrixFree<dim, double> &data = *system_matrix.get_matrix_free();
 *       FEEvaluation<dim, fe_degree>   phi(data);
 *       RightHandSide<dim>             rhs_func;
 *       Solution<dim>                  exact_solution;
 *       for (unsigned int cell = 0; cell < data.n_cell_batches(); ++cell)
 *         {
 *           phi.reinit(cell);
 *           for (const unsigned int q : phi.quadrature_point_indices())
 *             {
 *               VectorizedArray<double> rhs_val = VectorizedArray<double>();
 *               Point<dim, VectorizedArray<double>> point_batch =
 *                 phi.quadrature_point(q);
 *               for (unsigned int v = 0; v < VectorizedArray<double>::size(); ++v)
 *                 {
 *                   Point<dim> single_point;
 *                   for (unsigned int d = 0; d < dim; ++d)
 *                     single_point[d] = point_batch[d][v];
 *                   rhs_val[v] = rhs_func.value(single_point);
 *                 }
 *               phi.submit_value(rhs_val, q);
 *             }
 *           phi.integrate_scatter(EvaluationFlags::values, system_rhs);
 *         }
 *   
 * @endcode
 * 
 * Secondly, we also need to apply the Dirichlet and Neumann boundary
 * conditions. This function is the missing part of to the function
 * `LaplaceOperator::apply_boundary()` function once the exterior solution
 * values @f$u^+ = -u^- + 2 g_\text{D}@f$ and @f$\mathbf{n}^-\cdot \nabla u^+ =
 * \mathbf{n}^-\cdot \nabla u^-@f$ on Dirichlet boundaries and @f$u^+=u^-@f$ and
 * @f$\mathbf{n}^-\cdot \nabla u^+ = -\mathbf{n}^-\cdot \nabla u^- + 2
 * g_\text{N}@f$ on Neumann boundaries are inserted and expanded in terms of
 * the boundary functions @f$g_\text{D}@f$ and @f$g_\text{N}@f$. One thing to
 * remember is that we move the boundary conditions to the right hand
 * side, so the sign is the opposite from what we imposed on the solution
 * part.
 *     
 
 * 
 * We could have issued both the cell and the boundary part through a
 * MatrixFree::loop part, but we choose to manually write the full loop
 * over all faces to learn how the index layout of face indices is set up
 * in MatrixFree: Both the inner faces and the boundary faces share the
 * index range, and all batches of inner faces have lower numbers than the
 * batches of boundary cells. A single index for both variants allows us
 * to easily use the same data structure FEFaceEvaluation for both cases
 * that attaches to the same data field, just at different positions. The
 * number of inner face batches (where a batch is due to the combination
 * of several faces into one for vectorization) is given by
 * MatrixFree::n_inner_face_batches(), whereas the number of boundary face
 * batches is given by MatrixFree::n_boundary_face_batches().
 * 
 * @code
 *       FEFaceEvaluation<dim, fe_degree> phi_face(data, true);
 *       for (unsigned int face = data.n_inner_face_batches();
 *            face < data.n_inner_face_batches() + data.n_boundary_face_batches();
 *            ++face)
 *         {
 *           phi_face.reinit(face);
 *   
 *           const VectorizedArray<double> inverse_length_normal_to_face = std::abs(
 *             (phi_face.normal_vector(0) * phi_face.inverse_jacobian(0))[dim - 1]);
 *           const VectorizedArray<double> sigma =
 *             inverse_length_normal_to_face * system_matrix.get_penalty_factor();
 *   
 *           for (const unsigned int q : phi_face.quadrature_point_indices())
 *             {
 *               VectorizedArray<double> test_value = VectorizedArray<double>(),
 *                                       test_normal_derivative =
 *                                         VectorizedArray<double>();
 *               Point<dim, VectorizedArray<double>> point_batch =
 *                 phi_face.quadrature_point(q);
 *   
 *               for (unsigned int v = 0; v < VectorizedArray<double>::size(); ++v)
 *                 {
 *                   Point<dim> single_point;
 *                   for (unsigned int d = 0; d < dim; ++d)
 *                     single_point[d] = point_batch[d][v];
 *   
 * @endcode
 * 
 * The MatrixFree class lets us query the boundary_id of the
 * current face batch. Remember that MatrixFree sets up the
 * batches for vectorization such that all faces within a
 * batch have the same properties, which includes their
 * `boundary_id`. Thus, we can query that id here for the
 * current face index `face` and either impose the Dirichlet
 * case (where we add something to the function value) or the
 * Neumann case (where we add something to the normal
 * derivative).
 * 
 * @code
 *                   if (data.get_boundary_id(face) == 0)
 *                     test_value[v] = 2.0 * exact_solution.value(single_point);
 *                   else
 *                     {
 *                       Tensor<1, dim> normal;
 *                       for (unsigned int d = 0; d < dim; ++d)
 *                         normal[d] = phi_face.normal_vector(q)[d][v];
 *                       test_normal_derivative[v] =
 *                         -normal * exact_solution.gradient(single_point);
 *                     }
 *                 }
 *               phi_face.submit_value(test_value * sigma - test_normal_derivative,
 *                                     q);
 *               phi_face.submit_normal_derivative(-0.5 * test_value, q);
 *             }
 *           phi_face.integrate_scatter(EvaluationFlags::values |
 *                                        EvaluationFlags::gradients,
 *                                      system_rhs);
 *         }
 *   
 * @endcode
 * 
 * Since we have manually run the loop over cells rather than using
 * MatrixFree::loop(), we must not forget to perform the data exchange
 * with MPI - or actually, we would not need that for DG elements here
 * because each cell carries its own degrees of freedom and cell and
 * boundary integrals only evaluate quantities on the locally owned
 * cells. The coupling to neighboring subdomain only comes in by the inner
 * face integrals, which we have not done here. That said, it does not
 * hurt to call this function here, so we do it as a reminder of what
 * happens inside MatrixFree::loop().
 * 
 * @code
 *       system_rhs.compress(VectorOperation::add);
 *       setup_time += time.wall_time();
 *       time_details << "Compute right hand side       " << time.wall_time()
 *                    << " s\n";
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * The `solve()` function is copied almost verbatim from @ref step_37 "step-37". We set up
 * the same multigrid ingredients, namely the level transfer, a smoother,
 * and a coarse grid solver. The only difference is the fact that we do not
 * use the diagonal of the Laplacian for the preconditioner of the Chebyshev
 * iteration used for smoothing, but instead our newly resolved class
 * `%PreconditionBlockJacobi`. The mechanisms are the same, though.
 * 
 * @code
 *     template <int dim, int fe_degree>
 *     void LaplaceProblem<dim, fe_degree>::solve()
 *     {
 *       Timer                            time;
 *       MGTransferMatrixFree<dim, float> mg_transfer;
 *       mg_transfer.build(dof_handler);
 *       setup_time += time.wall_time();
 *       time_details << "MG build transfer time        " << time.wall_time()
 *                    << " s\n";
 *       time.restart();
 *   
 *       using SmootherType =
 *         PreconditionChebyshev<LevelMatrixType,
 *                               LinearAlgebra::distributed::Vector<float>,
 *                               PreconditionBlockJacobi<dim, fe_degree, float>>;
 *       mg::SmootherRelaxation<SmootherType,
 *                              LinearAlgebra::distributed::Vector<float>>
 *                                                            mg_smoother;
 *       MGLevelObject<typename SmootherType::AdditionalData> smoother_data;
 *       smoother_data.resize(0, triangulation.n_global_levels() - 1);
 *       for (unsigned int level = 0; level < triangulation.n_global_levels();
 *            ++level)
 *         {
 *           if (level > 0)
 *             {
 *               smoother_data[level].smoothing_range     = 15.;
 *               smoother_data[level].degree              = 3;
 *               smoother_data[level].eig_cg_n_iterations = 10;
 *             }
 *           else
 *             {
 *               smoother_data[0].smoothing_range = 2e-2;
 *               smoother_data[0].degree          = numbers::invalid_unsigned_int;
 *               smoother_data[0].eig_cg_n_iterations = mg_matrices[0].m();
 *             }
 *           smoother_data[level].preconditioner =
 *             std::make_shared<PreconditionBlockJacobi<dim, fe_degree, float>>();
 *           smoother_data[level].preconditioner->initialize(mg_matrices[level]);
 *         }
 *       mg_smoother.initialize(mg_matrices, smoother_data);
 *   
 *       MGCoarseGridApplySmoother<LinearAlgebra::distributed::Vector<float>>
 *         mg_coarse;
 *       mg_coarse.initialize(mg_smoother);
 *   
 *       mg::Matrix<LinearAlgebra::distributed::Vector<float>> mg_matrix(
 *         mg_matrices);
 *   
 *       Multigrid<LinearAlgebra::distributed::Vector<float>> mg(
 *         mg_matrix, mg_coarse, mg_transfer, mg_smoother, mg_smoother);
 *   
 *       PreconditionMG<dim,
 *                      LinearAlgebra::distributed::Vector<float>,
 *                      MGTransferMatrixFree<dim, float>>
 *         preconditioner(dof_handler, mg, mg_transfer);
 *   
 *       SolverControl solver_control(10000, 1e-12 * system_rhs.l2_norm());
 *       SolverCG<LinearAlgebra::distributed::Vector<double>> cg(solver_control);
 *       setup_time += time.wall_time();
 *       time_details << "MG build smoother time        " << time.wall_time()
 *                    << "s\n";
 *       pcout << "Total setup time              " << setup_time << " s\n";
 *   
 *       time.reset();
 *       time.start();
 *       cg.solve(system_matrix, solution, system_rhs, preconditioner);
 *   
 *       pcout << "Time solve (" << solver_control.last_step() << " iterations)    "
 *             << time.wall_time() << " s" << std::endl;
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * Since we have solved a problem with analytic solution, we want to verify
 * the correctness of our implementation by computing the L2 error of the
 * numerical result against the analytic solution.
 * 
 
 * 
 * 
 * @code
 *     template <int dim, int fe_degree>
 *     void LaplaceProblem<dim, fe_degree>::analyze_results() const
 *     {
 *       Vector<float> error_per_cell(triangulation.n_active_cells());
 *       VectorTools::integrate_difference(mapping,
 *                                         dof_handler,
 *                                         solution,
 *                                         Solution<dim>(),
 *                                         error_per_cell,
 *                                         QGauss<dim>(fe.degree + 2),
 *                                         VectorTools::L2_norm);
 *       pcout << "Verification via L2 error:    "
 *             << std::sqrt(
 *                  Utilities::MPI::sum(error_per_cell.norm_sqr(), MPI_COMM_WORLD))
 *             << std::endl;
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * The `run()` function sets up the initial grid and then runs the multigrid
 * program in the usual way. As a domain, we choose a rectangle with
 * periodic boundary conditions in the @f$x@f$-direction, a Dirichlet condition
 * on the front face in @f$y@f$ direction (i.e., the face with index number 2,
 * with boundary id equal to 0), and Neumann conditions on the back face as
 * well as the two faces in @f$z@f$ direction for the 3d case (with boundary id
 * equal to 1). The extent of the domain is a bit different in the @f$x@f$
 * direction (where we want to achieve a periodic solution given the
 * definition of `Solution`) as compared to the @f$y@f$ and @f$z@f$ directions.
 * 
 
 * 
 * 
 * @code
 *     template <int dim, int fe_degree>
 *     void LaplaceProblem<dim, fe_degree>::run()
 *     {
 *       const unsigned int n_ranks =
 *         Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD);
 *       pcout << "Running with " << n_ranks << " MPI process"
 *             << (n_ranks > 1 ? "es" : "") << ", element " << fe.get_name()
 *             << std::endl
 *             << std::endl;
 *       for (unsigned int cycle = 0; cycle < 9 - dim; ++cycle)
 *         {
 *           pcout << "Cycle " << cycle << std::endl;
 *   
 *           if (cycle == 0)
 *             {
 *               Point<dim> upper_right;
 *               upper_right[0] = 2.5;
 *               for (unsigned int d = 1; d < dim; ++d)
 *                 upper_right[d] = 2.8;
 *               GridGenerator::hyper_rectangle(triangulation,
 *                                              Point<dim>(),
 *                                              upper_right);
 *               triangulation.begin_active()->face(0)->set_boundary_id(10);
 *               triangulation.begin_active()->face(1)->set_boundary_id(11);
 *               triangulation.begin_active()->face(2)->set_boundary_id(0);
 *               for (unsigned int f = 3;
 *                    f < triangulation.begin_active()->n_faces();
 *                    ++f)
 *                 triangulation.begin_active()->face(f)->set_boundary_id(1);
 *   
 *               std::vector<GridTools::PeriodicFacePair<
 *                 typename Triangulation<dim>::cell_iterator>>
 *                 periodic_faces;
 *               GridTools::collect_periodic_faces(
 *                 triangulation, 10, 11, 0, periodic_faces);
 *               triangulation.add_periodicity(periodic_faces);
 *   
 *               triangulation.refine_global(6 - 2 * dim);
 *             }
 *           triangulation.refine_global(1);
 *           setup_system();
 *           compute_rhs();
 *           solve();
 *           analyze_results();
 *           pcout << std::endl;
 *         };
 *     }
 *   } // namespace Step59
 *   
 *   
 *   
 * @endcode
 * 
 * There is nothing unexpected in the `main()` function. We call `MPI_Init()`
 * through the `MPI_InitFinalize` class, pass on the two parameters on the
 * dimension and the degree set at the top of the file, and run the Laplace
 * problem.
 * 
 
 * 
 * 
 * @code
 *   int main(int argc, char *argv[])
 *   {
 *     try
 *       {
 *         using namespace Step59;
 *   
 *         Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1);
 *   
 *         LaplaceProblem<dimension, degree_finite_element> laplace_problem;
 *         laplace_problem.run();
 *       }
 *     catch (std::exception &exc)
 *       {
 *         std::cerr << std::endl
 *                   << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         std::cerr << "Exception on processing: " << std::endl
 *                   << exc.what() << std::endl
 *                   << "Aborting!" << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         return 1;
 *       }
 *     catch (...)
 *       {
 *         std::cerr << std::endl
 *                   << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         std::cerr << "Unknown exception!" << std::endl
 *                   << "Aborting!" << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         return 1;
 *       }
 *   
 *     return 0;
 *   }
 * @endcode
<a name="step_59-Results"></a><h1>Results</h1>
 
 
<a name="step_59-Programoutput"></a><h3>Program output</h3>
 
 
Like in @ref step_37 "step-37", we evaluate the multigrid solver in terms of run time.  In
two space dimensions with elements of degree 8, a possible output could look
as follows:
@code
Running with 12 MPI processes, element FE_DGQHermite<2>(8)
 
Cycle 0
Number of degrees of freedom: 5184
Total setup time              0.0282445 s
Time solve (14 iterations)    0.0110712 s
Verification via L2 error:    1.66232e-07
 
Cycle 1
Number of degrees of freedom: 20736
Total setup time              0.0126282 s
Time solve (14 iterations)    0.0157021 s
Verification via L2 error:    2.91505e-10
 
Cycle 2
Number of degrees of freedom: 82944
Total setup time              0.0227573 s
Time solve (14 iterations)    0.026568 s
Verification via L2 error:    6.64514e-13
 
Cycle 3
Number of degrees of freedom: 331776
Total setup time              0.0604685 s
Time solve (14 iterations)    0.0628356 s
Verification via L2 error:    5.57513e-13
 
Cycle 4
Number of degrees of freedom: 1327104
Total setup time              0.154359 s
Time solve (13 iterations)    0.219555 s
Verification via L2 error:    3.08139e-12
 
Cycle 5
Number of degrees of freedom: 5308416
Total setup time              0.467764 s
Time solve (13 iterations)    1.1821 s
Verification via L2 error:    3.90334e-12
 
Cycle 6
Number of degrees of freedom: 21233664
Total setup time              1.73263 s
Time solve (13 iterations)    5.21054 s
Verification via L2 error:    4.94543e-12
@endcode
 
Like in @ref step_37 "step-37", the number of CG iterations remains constant with increasing
problem size. The iteration counts are a bit higher, which is because we use a
lower degree of the Chebyshev polynomial (2 vs 5 in @ref step_37 "step-37") and because the
interior penalty discretization has a somewhat larger spread in
eigenvalues. Nonetheless, 13 iterations to reduce the residual by 12 orders of
magnitude, or almost a factor of 9 per iteration, indicates an overall very
efficient method. In particular, we can solve a system with 21 million degrees
of freedom in 5 seconds when using 12 cores, which is a very good
efficiency. Of course, in 2D we are well inside the regime of roundoff for a
polynomial degree of 8 &ndash; as a matter of fact, around 83k DoFs or 0.025s
would have been enough to fully converge this (simple) analytic solution
here.
 
Not much changes if we run the program in three spatial dimensions, except for
the fact that we now use do something more useful with the higher polynomial
degree and increasing mesh sizes, as the roundoff errors are only obtained at
the finest mesh. Still, it is remarkable that we can solve a 3D Laplace
problem with a wave of three periods to roundoff accuracy on a twelve-core
machine pretty easily - using about 3.5 GB of memory in total for the second
to largest case with 24m DoFs, taking not more than eight seconds. The largest
case uses 30GB of memory with 191m DoFs.
 
@code
Running with 12 MPI processes, element FE_DGQHermite<3>(8)
 
Cycle 0
Number of degrees of freedom: 5832
Total setup time              0.0210681 s
Time solve (15 iterations)    0.0956945 s
Verification via L2 error:    0.0297194
 
Cycle 1
Number of degrees of freedom: 46656
Total setup time              0.0452428 s
Time solve (15 iterations)    0.113827 s
Verification via L2 error:    9.55733e-05
 
Cycle 2
Number of degrees of freedom: 373248
Total setup time              0.190423 s
Time solve (15 iterations)    0.218309 s
Verification via L2 error:    2.6868e-07
 
Cycle 3
Number of degrees of freedom: 2985984
Total setup time              0.627914 s
Time solve (15 iterations)    1.0595 s
Verification via L2 error:    4.6918e-10
 
Cycle 4
Number of degrees of freedom: 23887872
Total setup time              2.85215 s
Time solve (15 iterations)    8.30576 s
Verification via L2 error:    9.38583e-13
 
Cycle 5
Number of degrees of freedom: 191102976
Total setup time              16.1324 s
Time solve (15 iterations)    65.57 s
Verification via L2 error:    3.17875e-13
@endcode
 
<a name="step_59-Comparisonofefficiencyatdifferentpolynomialdegrees"></a><h3>Comparison of efficiency at different polynomial degrees</h3>
 
 
In the introduction and in-code comments, it was mentioned several times that
high orders are treated very efficiently with the FEEvaluation and
FEFaceEvaluation evaluators. Now, we want to substantiate these claims by
looking at the throughput of the 3D multigrid solver for various polynomial
degrees. We collect the times as follows: We first run a solver at problem
size close to ten million, indicated in the first four table rows, and record
the timings. Then, we normalize the throughput by recording the number of
million degrees of freedom solved per second (MDoFs/s) to be able to compare
the efficiency of the different degrees, which is computed by dividing the
number of degrees of freedom by the solver time.
 
<table align="center" class="doxtable">
  <tr>
   <th>degree</th>
   <th>1</th>
   <th>2</th>
   <th>3</th>
   <th>4</th>
   <th>5</th>
   <th>6</th>
   <th>7</th>
   <th>8</th>
   <th>9</th>
   <th>10</th>
   <th>11</th>
   <th>12</th>
  </tr>
  <tr>
   <th>Number of DoFs</th>
   <td>2097152</td>
   <td>7077888</td>
   <td>16777216</td>
   <td>32768000</td>
   <td>7077888</td>
   <td>11239424</td>
   <td>16777216</td>
   <td>23887872</td>
   <td>32768000</td>
   <td>43614208</td>
   <td>7077888</td>
   <td>8998912</td>
  </tr>
  <tr>
   <th>Number of iterations</th>
   <td>13</td>
   <td>12</td>
   <td>12</td>
   <td>12</td>
   <td>13</td>
   <td>13</td>
   <td>15</td>
   <td>15</td>
   <td>17</td>
   <td>19</td>
   <td>18</td>
   <td>18</td>
  </tr>
  <tr>
   <th>Solver time [s]</th>
   <td>0.713</td>
   <td>2.150</td>
   <td>4.638</td>
   <td>8.803</td>
   <td>2.041</td>
   <td>3.295</td>
   <td>5.723</td>
   <td>8.306</td>
   <td>12.75</td>
   <td>19.25</td>
   <td>3.530</td>
   <td>4.814</td>
  </tr>
  <tr>
   <th>MDoFs/s</th>
   <td>2.94</td>
   <td>3.29</td>
   <td>3.62</td>
   <td>3.72</td>
   <td>3.47</td>
   <td>3.41</td>
   <td>2.93</td>
   <td>2.88</td>
   <td>2.57</td>
   <td>2.27</td>
   <td>2.01</td>
   <td>1.87</td>
  </tr>
</table>
 
We clearly see how the efficiency per DoF initially improves until it reaches
a maximum for the polynomial degree @f$k=4@f$. This effect is surprising, not only
because higher polynomial degrees often yield a vastly better solution, but
especially also when having matrix-based schemes in mind where the denser
coupling at higher degree leads to a monotonously decreasing throughput (and a
drastic one in 3D, with @f$k=4@f$ being more than ten times slower than
@f$k=1@f$!). For higher degrees, the throughput decreases a bit, which is both due
to an increase in the number of iterations (going from 12 at @f$k=2,3,4@f$ to 19
at @f$k=10@f$) and due to the @f$\mathcal O(k)@f$ complexity of operator
evaluation. Nonetheless, efficiency as the time to solution would be still
better for higher polynomial degrees because they have better convergence rates (at least
for problems as simple as this one): For @f$k=12@f$, we reach roundoff accuracy
already with 1 million DoFs (solver time less than a second), whereas for @f$k=8@f$
we need 24 million DoFs and 8 seconds. For @f$k=5@f$, the error is around
@f$10^{-9}@f$ with 57m DoFs and thus still far away from roundoff, despite taking 16
seconds.
 
Note that the above numbers are a bit pessimistic because they include the
time it takes the Chebyshev smoother to compute an eigenvalue estimate, which
is around 10 percent of the solver time. If the system is solved several times
(as e.g. common in fluid dynamics), this eigenvalue cost is only paid once and
faster times become available.
 
<a name="step_59-Evaluationofefficiencyofingredients"></a><h3>Evaluation of efficiency of ingredients</h3>
 
 
Finally, we take a look at some of the special ingredients presented in this
tutorial program, namely the FE_DGQHermite basis in particular and the
specification of MatrixFree::DataAccessOnFaces. In the following table, the
third row shows the optimized solver above, the fourth row shows the timings
with only the MatrixFree::DataAccessOnFaces set to `unspecified` rather than
the optimal `gradients`, and the last one with replacing FE_DGQHermite by the
basic FE_DGQ elements where both the MPI exchange are more expensive and the
operations done by FEFaceEvaluation::gather_evaluate() and
FEFaceEvaluation::integrate_scatter().
 
<table align="center" class="doxtable">
  <tr>
   <th>degree</th>
   <th>1</th>
   <th>2</th>
   <th>3</th>
   <th>4</th>
   <th>5</th>
   <th>6</th>
   <th>7</th>
   <th>8</th>
   <th>9</th>
   <th>10</th>
   <th>11</th>
   <th>12</th>
  </tr>
  <tr>
   <th>Number of DoFs</th>
   <td>2097152</td>
   <td>7077888</td>
   <td>16777216</td>
   <td>32768000</td>
   <td>7077888</td>
   <td>11239424</td>
   <td>16777216</td>
   <td>23887872</td>
   <td>32768000</td>
   <td>43614208</td>
   <td>7077888</td>
   <td>8998912</td>
  </tr>
  <tr>
   <th>Solver time optimized as in tutorial [s]</th>
   <td>0.713</td>
   <td>2.150</td>
   <td>4.638</td>
   <td>8.803</td>
   <td>2.041</td>
   <td>3.295</td>
   <td>5.723</td>
   <td>8.306</td>
   <td>12.75</td>
   <td>19.25</td>
   <td>3.530</td>
   <td>4.814</td>
  </tr>
  <tr>
   <th>Solver time MatrixFree::DataAccessOnFaces::unspecified [s]</th>
   <td>0.711</td>
   <td>2.151</td>
   <td>4.675</td>
   <td>8.968</td>
   <td>2.243</td>
   <td>3.655</td>
   <td>6.277</td>
   <td>9.082</td>
   <td>13.50</td>
   <td>20.05</td>
   <td>3.817</td>
   <td>5.178</td>
  </tr>
  <tr>
   <th>Solver time FE_DGQ [s]</th>
   <td>0.712</td>
   <td>2.041</td>
   <td>5.066</td>
   <td>9.335</td>
   <td>2.379</td>
   <td>3.802</td>
   <td>6.564</td>
   <td>9.714</td>
   <td>14.54</td>
   <td>22.76</td>
   <td>4.148</td>
   <td>5.857</td>
  </tr>
</table>
 
The data in the table shows that not using MatrixFree::DataAccessOnFaces
increases costs by around 10% for higher polynomial degrees. For lower
degrees, the difference is obviously less pronounced because the
volume-to-surface ratio is more beneficial and less data needs to be
exchanged. The difference is larger when looking at the matrix-vector product
only, rather than the full multigrid solver shown here, with around 20% worse
timings just because of the MPI communication.
 
For @f$k=1@f$ and @f$k=2@f$, the Hermite-like basis functions do obviously not really
pay off (indeed, for @f$k=1@f$ the polynomials are exactly the same as for FE_DGQ)
and the results are similar as with the FE_DGQ basis. However, for degrees
starting at three, we see an increasing advantage for FE_DGQHermite, showing
the effectiveness of these basis functions.
 
<a name="step_59-Possibilitiesforextension"></a><h3>Possibilities for extension</h3>
 
 
As mentioned in the introduction, the fast diagonalization method as realized
here is tied to a Cartesian mesh with constant coefficients. When dealing with
meshes that contain deformed cells or with variable coefficients, it is common
to determine a nearby Cartesian mesh cell as an approximation. This can be
done with the class TensorProductMatrixSymmetricSumCollection. Here, one can
insert cell matrices similarly to the PreconditionBlockJacobi::initialize()
function of this tutorial program. The benefit of the collection class is that
cells on which the coefficient of the PDE has the same value can re-use the
same Laplacian matrix, which reduces the memory consumption for the inverse
matrices. As compared to the algorithm implemented in this tutorial program,
one would define the length scales as the distances between opposing
faces. For continuous elements, the code project <a
href=https://github.com/peterrum/dealii-dd-and-schwarz">Cache-optimized and
low-overhead implementations of multigrid smoothers for high-order FEM
computations</a> presents the computation for continuous elements. There is
currently no infrastructure in deal.II to automatically generate the 1D
matrices for discontinuous elements with SIP-DG discretization, as opposed to
continuous elements, where we provide
TensorProductMatrixCreator::create_laplace_tensor_product_matrix().
 
Another way of extending the program would be to include support for adaptive
meshes. While the classical approach of defining interface operations at edges
of different refinement level, as discussed in @ref step_39 "step-39", is one possibility,
for Poisson-type problems another option is typically more beneficial. Using
the class MGTransferGlobalCoarsening, which is explained in the @ref step_75 "step-75"
tutorial program, one can deal with meshes of hanging nodes on all levels. An
algorithmic improvement can be obtained by combining the discontinuous
function space with the auxiliary continuous finite element space of the same
polynomial degree. This idea, introduced by Antonietti et al.
@cite antonietti2016uniform in 2016, allows making the multigrid convergence
independent of the penalty parameter. As demonstrated by Fehn et al.
@cite fehn2020hybrid, this also gives considerably lower iteration counts than
a multigrid solver directly working on levels with discontinuous function
spaces. The latter work also proposes p-multigrid techniques and combination
with algebraic multigrid coarse spaces as a means to efficiently solve Poisson
problems with high-order discontinuous Galerkin discretizations on complicated
geometries, representing the current state-of-the-art for simple Poisson-type
problems. The class MGTransferGlobalCoarsening provides features for each of
these three coarsening variants, the discontinuous-continuous auxiliary
function concept, p-multigrid, and traditional h-multigrid. The main
ingredient is to define an appropriate MGTwoLevelTransfer object and call
MGTwoLevelTransfer::reinit_geometric_transfer() or
MGTwoLevelTransfer::reinit_polynomial_transfer(), respectively.
 *
 *
<a name="step_59-PlainProg"></a>
<h1> The plain program</h1>
@include "step-59.cc"
*/