Reference documentation for deal.II version GIT relicensing-422-gb369f187d8 2024-04-17 18:10:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
step-59.h
Go to the documentation of this file.
1true,
816 *   }
817 *  
818 *  
819 *  
820 * @endcode
821 *
822 * Since the Laplacian is symmetric, the `Tvmult()` (needed by the multigrid
823 * smoother interfaces) operation is simply forwarded to the `vmult()` case.
824 *
825
826 *
827 *
828 * @code
829 *   template <int dim, int fe_degree, typename number>
830 *   void LaplaceOperator<dim, fe_degree, number>::Tvmult(
832 *   const LinearAlgebra::distributed::Vector<number> &src) const
833 *   {
834 *   vmult(dst, src);
835 *   }
836 *  
837 *  
838 *  
839 * @endcode
840 *
841 * The cell operation is very similar to @ref step_37 "step-37". We do not use a
842 * coefficient here, though. The second difference is that we replaced the
843 * two steps of FEEvaluation::read_dof_values() followed by
844 * FEEvaluation::evaluate() by a single function call
845 * FEEvaluation::gather_evaluate() which internally calls the sequence of
846 * the two individual methods. Likewise, FEEvaluation::integrate_scatter()
847 * implements the sequence of FEEvaluation::integrate() followed by
848 * FEEvaluation::distribute_local_to_global(). In this case, these new
849 * functions merely save two lines of code. However, we use them for the
850 * analogy with FEFaceEvaluation where they are more important as
851 * explained below.
852 *
853
854 *
855 *
856 * @code
857 *   template <int dim, int fe_degree, typename number>
858 *   void LaplaceOperator<dim, fe_degree, number>::apply_cell(
859 *   const MatrixFree<dim, number> &data,
860 *   LinearAlgebra::distributed::Vector<number> &dst,
861 *   const LinearAlgebra::distributed::Vector<number> &src,
862 *   const std::pair<unsigned int, unsigned int> &cell_range) const
863 *   {
865 *   for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
866 *   {
867 *   phi.reinit(cell);
868 *   phi.gather_evaluate(src, EvaluationFlags::gradients);
869 *   for (const unsigned int q : phi.quadrature_point_indices())
870 *   phi.submit_gradient(phi.get_gradient(q), q);
871 *   phi.integrate_scatter(EvaluationFlags::gradients, dst);
872 *   }
873 *   }
874 *  
875 *  
876 *  
877 * @endcode
878 *
879 * The face operation implements the terms of the interior penalty method in
880 * analogy to @ref step_39 "step-39", as explained in the introduction. We need two
881 * evaluator objects for this task, one for handling the solution that comes
882 * from the cell on one of the two sides of an interior face, and one for
883 * handling the solution from the other side. The evaluators for face
884 * integrals are called FEFaceEvaluation and take a boolean argument in the
885 * second slot of the constructor to indicate which of the two sides the
886 * evaluator should belong two. In FEFaceEvaluation and MatrixFree, we call
887 * one of the two sides the `interior` one and the other the `exterior`
888 * one. The name `exterior` refers to the fact that the evaluator from both
889 * sides will return the same normal vector. For the `interior` side, the
890 * normal vector points outwards, whereas it points inwards on the other
891 * side, and is opposed to the outer normal vector of that cell. Apart from
892 * the new class name, we again get a range of items to work with in
893 * analogy to what was discussed in @ref step_37 "step-37", but for the interior faces in
894 * this case. Note that the data structure of MatrixFree forms batches of
895 * faces that are analogous to the batches of cells for the cell
896 * integrals. All faces within a batch involve different cell numbers but
897 * have the face number within the reference cell, have the same refinement
898 * configuration (no refinement or the same subface), and the same
899 * orientation, to keep SIMD operations simple and efficient.
900 *
901
902 *
903 * Note that there is no implied meaning in interior versus exterior except
904 * the logic decision of the orientation of the normal, which is pretty
905 * random internally. One can in no way rely on a certain pattern of
906 * assigning interior versus exterior flags, as the decision is made for the
907 * sake of access regularity and uniformity in the MatrixFree setup
908 * routines. Since most sane DG methods are conservative, i.e., fluxes look
909 * the same from both sides of an interface, the mathematics are unaltered
910 * if the interior/exterior flags are switched and normal vectors get the
911 * opposite sign.
912 *
913
914 *
915 *
916 * @code
917 *   template <int dim, int fe_degree, typename number>
918 *   void LaplaceOperator<dim, fe_degree, number>::apply_face(
919 *   const MatrixFree<dim, number> &data,
922 *   const std::pair<unsigned int, unsigned int> &face_range) const
923 *   {
925 *   true);
927 *   false);
928 *   for (unsigned int face = face_range.first; face < face_range.second; ++face)
929 *   {
930 * @endcode
931 *
932 * On a given batch of faces, we first update the pointers to the
933 * current face and then access the vector. As mentioned above, we
934 * combine the vector access with the evaluation. In the case of face
935 * integrals, the data access into the vector can be reduced for the
936 * special case of an FE_DGQHermite basis as explained for the data
937 * exchange above: Since only @f$2(k+1)^{d-1}@f$ out of the @f$(k+1)^d@f$ cell
938 * degrees of freedom get multiplied by a non-zero value or derivative
939 * of a shape function, this structure can be utilized for the
940 * evaluation, significantly reducing the data access. The reduction
941 * of the data access is not only beneficial because it reduces the
942 * data in flight and thus helps caching, but also because the data
943 * access to faces is often more irregular than for cell integrals when
944 * gathering values from cells that are farther apart in the index
945 * list of cells.
946 *
947 * @code
948 *   phi_inner.reinit(face);
949 *   phi_inner.gather_evaluate(src,
952 *   phi_outer.reinit(face);
953 *   phi_outer.gather_evaluate(src,
956 *  
957 * @endcode
958 *
959 * The next two statements compute the penalty parameter for the
960 * interior penalty method. As explained in the introduction, we would
961 * like to have a scaling like @f$\frac{1}{h_\text{i}}@f$ of the length
962 * @f$h_\text{i}@f$ normal to the face. For a general non-Cartesian mesh,
963 * this length must be computed by the product of the inverse Jacobian
964 * times the normal vector in real coordinates. From this vector of
965 * `dim` components, we must finally pick the component that is
966 * oriented normal to the reference cell. In the geometry data stored
967 * in MatrixFree, a permutation of the components in the Jacobian is
968 * applied such that this latter direction is always the last
969 * component `dim-1` (this is beneficial because reference-cell
970 * derivative sorting can be made agnostic of the direction of the
971 * face). This means that we can simply access the last entry `dim-1`
972 * and must not look up the local face number in
973 * `data.get_face_info(face).interior_face_no` and
974 * `data.get_face_info(face).exterior_face_no`. Finally, we must also
975 * take the absolute value of these factors as the normal could point
976 * into either positive or negative direction.
977 *
978 * @code
979 *   const VectorizedArray<number> inverse_length_normal_to_face =
980 *   0.5 * (std::abs((phi_inner.normal_vector(0) *
981 *   phi_inner.inverse_jacobian(0))[dim - 1]) +
982 *   std::abs((phi_outer.normal_vector(0) *
983 *   phi_outer.inverse_jacobian(0))[dim - 1]));
984 *   const VectorizedArray<number> sigma =
985 *   inverse_length_normal_to_face * get_penalty_factor();
986 *  
987 * @endcode
988 *
989 * In the loop over the quadrature points, we eventually compute all
990 * contributions to the interior penalty scheme. According to the
991 * formulas in the introduction, the value of the test function gets
992 * multiplied by the difference of the jump in the solution times the
993 * penalty parameter and the average of the normal derivative in real
994 * space. Since the two evaluators for interior and exterior sides get
995 * different signs due to the jump, we pass the result with a
996 * different sign here. The normal derivative of the test function
997 * gets multiplied by the negative jump in the solution between the
998 * interior and exterior side. This term, coined adjoint consistency
999 * term, must also include the factor of @f$\frac{1}{2}@f$ in the code in
1000 * accordance with its relation to the primal consistency term that
1001 * gets the factor of one half due to the average in the test function
1002 * slot.
1003 *
1004 * @code
1005 *   for (const unsigned int q : phi_inner.quadrature_point_indices())
1006 *   {
1007 *   const VectorizedArray<number> solution_jump =
1008 *   (phi_inner.get_value(q) - phi_outer.get_value(q));
1009 *   const VectorizedArray<number> average_normal_derivative =
1010 *   (phi_inner.get_normal_derivative(q) +
1011 *   phi_outer.get_normal_derivative(q)) *
1012 *   number(0.5);
1013 *   const VectorizedArray<number> test_by_value =
1014 *   solution_jump * sigma - average_normal_derivative;
1015 *  
1016 *   phi_inner.submit_value(test_by_value, q);
1017 *   phi_outer.submit_value(-test_by_value, q);
1018 *  
1019 *   phi_inner.submit_normal_derivative(-solution_jump * number(0.5), q);
1020 *   phi_outer.submit_normal_derivative(-solution_jump * number(0.5), q);
1021 *   }
1022 *  
1023 * @endcode
1024 *
1025 * Once we are done with the loop over quadrature points, we can do
1026 * the sum factorization operations for the integration loops on faces
1027 * and sum the results into the result vector, using the
1028 * `integrate_scatter` function. The name `scatter` reflects the
1029 * distribution of the vector data into scattered positions in the
1030 * vector using the same pattern as in `gather_evaluate`. Like before,
1031 * the combined integrate + write operation allows us to reduce the
1032 * data access.
1033 *
1034 * @code
1035 *   phi_inner.integrate_scatter(EvaluationFlags::values |
1037 *   dst);
1038 *   phi_outer.integrate_scatter(EvaluationFlags::values |
1040 *   dst);
1041 *   }
1042 *   }
1043 *  
1044 *  
1045 *  
1046 * @endcode
1047 *
1048 * The boundary face function follows by and large the interior face
1049 * function. The only difference is the fact that we do not have a separate
1050 * FEFaceEvaluation object that provides us with exterior values @f$u^+@f$, but
1051 * we must define them from the boundary conditions and interior values
1052 * @f$u^-@f$. As explained in the introduction, we use @f$u^+ = -u^- + 2
1053 * g_\text{D}@f$ and @f$\mathbf{n}^-\cdot \nabla u^+ = \mathbf{n}^-\cdot \nabla
1054 * u^-@f$ on Dirichlet boundaries and @f$u^+=u^-@f$ and @f$\mathbf{n}^-\cdot \nabla
1055 * u^+ = -\mathbf{n}^-\cdot \nabla u^- + 2 g_\text{N}@f$ on Neumann
1056 * boundaries. Since this operation implements the homogeneous part, i.e.,
1057 * the matrix-vector product, we must neglect the boundary functions
1058 * @f$g_\text{D}@f$ and @f$g_\text{N}@f$ here, and added them to the right hand side
1059 * in `LaplaceProblem::compute_rhs()`. Note that due to extension of the
1060 * solution @f$u^-@f$ to the exterior via @f$u^+@f$, we can keep all factors @f$0.5@f$
1061 * the same as in the inner face function, see also the discussion in
1062 * @ref step_39 "step-39".
1063 *
1064
1065 *
1066 * There is one catch at this point: The implementation below uses a boolean
1067 * variable `is_dirichlet` to switch between the Dirichlet and the Neumann
1068 * cases. However, we solve a problem where we also want to impose periodic
1069 * boundary conditions on some boundaries, namely along those in the @f$x@f$
1070 * direction. One might wonder how those conditions should be handled
1071 * here. The answer is that MatrixFree automatically treats periodic
1072 * boundaries as what they are technically, namely an inner face where the
1073 * solution values of two adjacent cells meet and must be treated by proper
1074 * numerical fluxes. Thus, all the faces on the periodic boundaries will
1075 * appear in the `apply_face()` function and not in this one.
1076 *
1077
1078 *
1079 *
1080 * @code
1081 *   template <int dim, int fe_degree, typename number>
1082 *   void LaplaceOperator<dim, fe_degree, number>::apply_boundary(
1083 *   const MatrixFree<dim, number> &data,
1086 *   const std::pair<unsigned int, unsigned int> &face_range) const
1087 *   {
1089 *   true);
1090 *   for (unsigned int face = face_range.first; face < face_range.second; ++face)
1091 *   {
1092 *   phi_inner.reinit(face);
1093 *   phi_inner.gather_evaluate(src,
1096 *  
1097 *   const VectorizedArray<number> inverse_length_normal_to_face = std::abs((
1098 *   phi_inner.normal_vector(0) * phi_inner.inverse_jacobian(0))[dim - 1]);
1099 *   const VectorizedArray<number> sigma =
1100 *   inverse_length_normal_to_face * get_penalty_factor();
1101 *  
1102 *   const bool is_dirichlet = (data.get_boundary_id(face) == 0);
1103 *  
1104 *   for (const unsigned int q : phi_inner.quadrature_point_indices())
1105 *   {
1106 *   const VectorizedArray<number> u_inner = phi_inner.get_value(q);
1107 *   const VectorizedArray<number> u_outer =
1108 *   is_dirichlet ? -u_inner : u_inner;
1109 *   const VectorizedArray<number> normal_derivative_inner =
1110 *   phi_inner.get_normal_derivative(q);
1111 *   const VectorizedArray<number> normal_derivative_outer =
1112 *   is_dirichlet ? normal_derivative_inner : -normal_derivative_inner;
1113 *   const VectorizedArray<number> solution_jump = (u_inner - u_outer);
1114 *   const VectorizedArray<number> average_normal_derivative =
1115 *   (normal_derivative_inner + normal_derivative_outer) * number(0.5);
1116 *   const VectorizedArray<number> test_by_value =
1117 *   solution_jump * sigma - average_normal_derivative;
1118 *   phi_inner.submit_normal_derivative(-solution_jump * number(0.5), q);
1119 *   phi_inner.submit_value(test_by_value, q);
1120 *   }
1121 *   phi_inner.integrate_scatter(EvaluationFlags::values |
1123 *   dst);
1124 *   }
1125 *   }
1126 *  
1127 *  
1128 *  
1129 * @endcode
1130 *
1131 * Next we turn to the preconditioner initialization. As explained in the
1132 * introduction, we want to construct an (approximate) inverse of the cell
1133 * matrices from a product of 1d mass and Laplace matrices. Our first task
1134 * is to compute the 1d matrices, which we do by first creating a 1d finite
1135 * element. Instead of anticipating FE_DGQHermite<1> here, we get the finite
1136 * element's name from DoFHandler, replace the @p dim argument (2 or 3) by 1
1137 * to create a 1d name, and construct the 1d element by using FETools.
1138 *
1139
1140 *
1141 *
1142 * @code
1143 *   template <int dim, int fe_degree, typename number>
1144 *   void PreconditionBlockJacobi<dim, fe_degree, number>::initialize(
1145 *   const LaplaceOperator<dim, fe_degree, number> &op)
1146 *   {
1147 *   data = op.get_matrix_free();
1148 *  
1149 *   std::string name = data->get_dof_handler().get_fe().get_name();
1150 *   name.replace(name.find('<') + 1, 1, "1");
1151 *   std::unique_ptr<FiniteElement<1>> fe_1d = FETools::get_fe_by_name<1>(name);
1152 *  
1153 * @endcode
1154 *
1155 * As for computing the 1d matrices on the unit element, we simply write
1156 * down what a typical assembly procedure over rows and columns of the
1157 * matrix as well as the quadrature points would do. We select the same
1158 * Laplace matrices once and for all using the coefficients 0.5 for
1159 * interior faces (but possibly scaled differently in different directions
1160 * as a result of the mesh). Thus, we make a slight mistake at the
1161 * Dirichlet boundary (where the correct factor would be 1 for the
1162 * derivative terms and 2 for the penalty term, see @ref step_39 "step-39") or at the
1163 * Neumann boundary where the factor should be zero. Since we only use
1164 * this class as a smoother inside a multigrid scheme, this error is not
1165 * going to have any significant effect and merely affects smoothing
1166 * quality.
1167 *
1168 * @code
1169 *   const unsigned int N = fe_degree + 1;
1170 *   FullMatrix<double> laplace_unscaled(N, N);
1171 *   std::array<Table<2, VectorizedArray<number>>, dim> mass_matrices;
1172 *   std::array<Table<2, VectorizedArray<number>>, dim> laplace_matrices;
1173 *   for (unsigned int d = 0; d < dim; ++d)
1174 *   {
1175 *   mass_matrices[d].reinit(N, N);
1176 *   laplace_matrices[d].reinit(N, N);
1177 *   }
1178 *  
1179 *   QGauss<1> quadrature(N);
1180 *   for (unsigned int i = 0; i < N; ++i)
1181 *   for (unsigned int j = 0; j < N; ++j)
1182 *   {
1183 *   double sum_mass = 0, sum_laplace = 0;
1184 *   for (unsigned int q = 0; q < quadrature.size(); ++q)
1185 *   {
1186 *   sum_mass += (fe_1d->shape_value(i, quadrature.point(q)) *
1187 *   fe_1d->shape_value(j, quadrature.point(q))) *
1188 *   quadrature.weight(q);
1189 *   sum_laplace += (fe_1d->shape_grad(i, quadrature.point(q))[0] *
1190 *   fe_1d->shape_grad(j, quadrature.point(q))[0]) *
1191 *   quadrature.weight(q);
1192 *   }
1193 *   for (unsigned int d = 0; d < dim; ++d)
1194 *   mass_matrices[d](i, j) = sum_mass;
1195 *  
1196 * @endcode
1197 *
1198 * The left and right boundary terms assembled by the next two
1199 * statements appear to have somewhat arbitrary signs, but those are
1200 * correct as can be verified by looking at @ref step_39 "step-39" and inserting
1201 * the value -1 and 1 for the normal vector in the 1d case.
1202 *
1203 * @code
1204 *   sum_laplace +=
1205 *   (1. * fe_1d->shape_value(i, Point<1>()) *
1206 *   fe_1d->shape_value(j, Point<1>()) * op.get_penalty_factor() +
1207 *   0.5 * fe_1d->shape_grad(i, Point<1>())[0] *
1208 *   fe_1d->shape_value(j, Point<1>()) +
1209 *   0.5 * fe_1d->shape_grad(j, Point<1>())[0] *
1210 *   fe_1d->shape_value(i, Point<1>()));
1211 *  
1212 *   sum_laplace +=
1213 *   (1. * fe_1d->shape_value(i, Point<1>(1.0)) *
1214 *   fe_1d->shape_value(j, Point<1>(1.0)) * op.get_penalty_factor() -
1215 *   0.5 * fe_1d->shape_grad(i, Point<1>(1.0))[0] *
1216 *   fe_1d->shape_value(j, Point<1>(1.0)) -
1217 *   0.5 * fe_1d->shape_grad(j, Point<1>(1.0))[0] *
1218 *   fe_1d->shape_value(i, Point<1>(1.0)));
1219 *  
1220 *   laplace_unscaled(i, j) = sum_laplace;
1221 *   }
1222 *  
1223 * @endcode
1224 *
1225 * Next, we go through the cells and pass the scaled matrices to
1226 * TensorProductMatrixSymmetricSum to actually compute the generalized
1227 * eigenvalue problem for representing the inverse: Since the matrix
1228 * approximation is constructed as @f$A\otimes M + M\otimes A@f$ and the
1229 * weights are constant for each element, we can apply all weights on the
1230 * Laplace matrix and simply keep the mass matrices unscaled. In the loop
1231 * over cells, we want to make use of the geometry compression provided by
1232 * the MatrixFree class and check if the current geometry is the same as
1233 * on the last cell batch, in which case there is nothing to do. This
1234 * compression can be accessed by
1235 * FEEvaluation::get_mapping_data_index_offset() once `reinit()` has been
1236 * called.
1237 *
1238
1239 *
1240 * Once we have accessed the inverse Jacobian through the FEEvaluation
1241 * access function (we take the one for the zeroth quadrature point as
1242 * they should be the same on all quadrature points for a Cartesian cell),
1243 * we check that it is diagonal and then extract the determinant of the
1244 * original Jacobian, i.e., the inverse of the determinant of the inverse
1245 * Jacobian, and set the weight as @f$\text{det}(J) / h_d^2@f$ according to
1246 * the 1d Laplacian times @f$d-1@f$ copies of the @ref GlossMassMatrix "mass matrix".
1247 *
1248 * @code
1249 *   cell_matrices.clear();
1250 *   FEEvaluation<dim, fe_degree, fe_degree + 1, 1, number> phi(*data);
1251 *   unsigned int old_mapping_data_index = numbers::invalid_unsigned_int;
1252 *   for (unsigned int cell = 0; cell < data->n_cell_batches(); ++cell)
1253 *   {
1254 *   phi.reinit(cell);
1255 *  
1256 *   if (phi.get_mapping_data_index_offset() == old_mapping_data_index)
1257 *   continue;
1258 *  
1259 *   Tensor<2, dim, VectorizedArray<number>> inverse_jacobian =
1260 *   phi.inverse_jacobian(0);
1261 *  
1262 *   for (unsigned int d = 0; d < dim; ++d)
1263 *   for (unsigned int e = 0; e < dim; ++e)
1264 *   if (d != e)
1265 *   for (unsigned int v = 0; v < VectorizedArray<number>::size(); ++v)
1266 *   AssertThrow(inverse_jacobian[d][e][v] == 0.,
1267 *   ExcNotImplemented());
1268 *  
1269 *   VectorizedArray<number> jacobian_determinant = inverse_jacobian[0][0];
1270 *   for (unsigned int e = 1; e < dim; ++e)
1271 *   jacobian_determinant *= inverse_jacobian[e][e];
1272 *   jacobian_determinant = 1. / jacobian_determinant;
1273 *  
1274 *   for (unsigned int d = 0; d < dim; ++d)
1275 *   {
1276 *   const VectorizedArray<number> scaling_factor =
1277 *   inverse_jacobian[d][d] * inverse_jacobian[d][d] *
1278 *   jacobian_determinant;
1279 *  
1280 * @endcode
1281 *
1282 * Once we know the factor by which we should scale the Laplace
1283 * matrix, we apply this weight to the unscaled DG Laplace matrix
1284 * and send the array to the class TensorProductMatrixSymmetricSum
1285 * for computing the generalized eigenvalue problem mentioned in
1286 * the introduction.
1287 *
1288
1289 *
1290 *
1291 * @code
1292 *   for (unsigned int i = 0; i < N; ++i)
1293 *   for (unsigned int j = 0; j < N; ++j)
1294 *   laplace_matrices[d](i, j) =
1295 *   scaling_factor * laplace_unscaled(i, j);
1296 *   }
1297 *   if (cell_matrices.size() <= phi.get_mapping_data_index_offset())
1298 *   cell_matrices.resize(phi.get_mapping_data_index_offset() + 1);
1299 *   cell_matrices[phi.get_mapping_data_index_offset()].reinit(
1300 *   mass_matrices, laplace_matrices);
1301 *   }
1302 *   }
1303 *  
1304 *  
1305 *  
1306 * @endcode
1307 *
1308 * The vmult function for the approximate block-Jacobi preconditioner is
1309 * very simple in the DG context: We simply need to read the values of the
1310 * current cell batch, apply the inverse for the given entry in the array of
1311 * tensor product matrix, and write the result back. In this loop, we
1312 * overwrite the content in `dst` rather than first setting the entries to
1313 * zero. This is legitimate for a DG method because every cell has
1314 * independent degrees of freedom. Furthermore, we manually write out the
1315 * loop over all cell batches, rather than going through
1316 * MatrixFree::cell_loop(). We do this because we know that we are not going
1317 * to need data exchange over the MPI network here as all computations are
1318 * done on the cells held locally on each processor.
1319 *
1320
1321 *
1322 *
1323 * @code
1324 *   template <int dim, int fe_degree, typename number>
1325 *   void PreconditionBlockJacobi<dim, fe_degree, number>::vmult(
1326 *   LinearAlgebra::distributed::Vector<number> &dst,
1327 *   const LinearAlgebra::distributed::Vector<number> &src) const
1328 *   {
1329 *   adjust_ghost_range_if_necessary(*data, dst);
1330 *   adjust_ghost_range_if_necessary(*data, src);
1331 *  
1332 *   FEEvaluation<dim, fe_degree, fe_degree + 1, 1, number> phi(*data);
1333 *   for (unsigned int cell = 0; cell < data->n_cell_batches(); ++cell)
1334 *   {
1335 *   phi.reinit(cell);
1336 *   phi.read_dof_values(src);
1337 *   cell_matrices[phi.get_mapping_data_index_offset()].apply_inverse(
1338 *   ArrayView<VectorizedArray<number>>(phi.begin_dof_values(),
1339 *   phi.dofs_per_cell),
1340 *   ArrayView<const VectorizedArray<number>>(phi.begin_dof_values(),
1341 *   phi.dofs_per_cell));
1342 *   phi.set_dof_values(dst);
1343 *   }
1344 *   }
1345 *  
1346 *  
1347 *  
1348 * @endcode
1349 *
1350 * The definition of the LaplaceProblem class is very similar to
1351 * @ref step_37 "step-37". One difference is the fact that we add the element degree as a
1352 * template argument to the class, which would allow us to more easily
1353 * include more than one degree in the same program by creating different
1354 * instances in the `main()` function. The second difference is the
1355 * selection of the element, FE_DGQHermite, which is specialized for this
1356 * kind of equations.
1357 *
1358
1359 *
1360 *
1361 * @code
1362 *   template <int dim, int fe_degree>
1363 *   class LaplaceProblem
1364 *   {
1365 *   public:
1366 *   LaplaceProblem();
1367 *   void run();
1368 *  
1369 *   private:
1370 *   void setup_system();
1371 *   void compute_rhs();
1372 *   void solve();
1373 *   void analyze_results() const;
1374 *  
1375 *   #ifdef DEAL_II_WITH_P4EST
1376 *   parallel::distributed::Triangulation<dim> triangulation;
1377 *   #else
1378 *   Triangulation<dim> triangulation;
1379 *   #endif
1380 *  
1381 *   FE_DGQHermite<dim> fe;
1382 *   DoFHandler<dim> dof_handler;
1383 *  
1384 *   MappingQ1<dim> mapping;
1385 *  
1386 *   using SystemMatrixType = LaplaceOperator<dim, fe_degree, double>;
1387 *   SystemMatrixType system_matrix;
1388 *  
1389 *   using LevelMatrixType = LaplaceOperator<dim, fe_degree, float>;
1390 *   MGLevelObject<LevelMatrixType> mg_matrices;
1391 *  
1392 *   LinearAlgebra::distributed::Vector<double> solution;
1393 *   LinearAlgebra::distributed::Vector<double> system_rhs;
1394 *  
1395 *   double setup_time;
1396 *   ConditionalOStream pcout;
1397 *   ConditionalOStream time_details;
1398 *   };
1399 *  
1400 *  
1401 *  
1402 *   template <int dim, int fe_degree>
1403 *   LaplaceProblem<dim, fe_degree>::LaplaceProblem()
1404 *   #ifdef DEAL_II_WITH_P4EST
1405 *   : triangulation(MPI_COMM_WORLD,
1406 *   Triangulation<dim>::limit_level_difference_at_vertices,
1407 *   parallel::distributed::Triangulation<
1408 *   dim>::construct_multigrid_hierarchy)
1409 *   #else
1410 *   : triangulation(Triangulation<dim>::limit_level_difference_at_vertices)
1411 *   #endif
1412 *   , fe(fe_degree)
1413 *   , dof_handler(triangulation)
1414 *   , setup_time(0.)
1415 *   , pcout(std::cout, Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
1416 *   , time_details(std::cout,
1417 *   false &&
1418 *   Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
1419 *   {}
1420 *  
1421 *  
1422 *  
1423 * @endcode
1424 *
1425 * The setup function differs in two aspects from @ref step_37 "step-37". The first is that
1426 * we do not need to interpolate any constraints for the discontinuous
1427 * ansatz space, and simply pass a dummy AffineConstraints object into
1428 * Matrixfree::reinit(). The second change arises because we need to tell
1429 * MatrixFree to also initialize the data structures for faces. We do this
1430 * by setting update flags for the inner and boundary faces,
1431 * respectively. On the boundary faces, we need both the function values,
1432 * their gradients, JxW values (for integration), the normal vectors, and
1433 * quadrature points (for the evaluation of the boundary conditions),
1434 * whereas we only need shape function values, gradients, JxW values, and
1435 * normal vectors for interior faces. The face data structures in MatrixFree
1436 * are always built as soon as one of `mapping_update_flags_inner_faces` or
1437 * `mapping_update_flags_boundary_faces` are different from the default
1438 * value `update_default` of UpdateFlags.
1439 *
1440
1441 *
1442 *
1443 * @code
1444 *   template <int dim, int fe_degree>
1445 *   void LaplaceProblem<dim, fe_degree>::setup_system()
1446 *   {
1447 *   Timer time;
1448 *   setup_time = 0;
1449 *  
1450 *   system_matrix.clear();
1451 *   mg_matrices.clear_elements();
1452 *  
1453 *   dof_handler.distribute_dofs(fe);
1454 *   dof_handler.distribute_mg_dofs();
1455 *  
1456 *   pcout << "Number of degrees of freedom: " << dof_handler.n_dofs()
1457 *   << std::endl;
1458 *  
1459 *   setup_time += time.wall_time();
1460 *   time_details << "Distribute DoFs " << time.wall_time() << " s"
1461 *   << std::endl;
1462 *   time.restart();
1463 *  
1464 *   AffineConstraints<double> dummy;
1465 *   dummy.close();
1466 *  
1467 *   {
1468 *   typename MatrixFree<dim, double>::AdditionalData additional_data;
1469 *   additional_data.tasks_parallel_scheme =
1470 *   MatrixFree<dim, double>::AdditionalData::none;
1471 *   additional_data.mapping_update_flags =
1472 *   (update_gradients | update_JxW_values | update_quadrature_points);
1473 *   additional_data.mapping_update_flags_inner_faces =
1474 *   (update_gradients | update_JxW_values | update_normal_vectors);
1475 *   additional_data.mapping_update_flags_boundary_faces =
1476 *   (update_gradients | update_JxW_values | update_normal_vectors |
1477 *   update_quadrature_points);
1478 *   const auto system_mf_storage =
1479 *   std::make_shared<MatrixFree<dim, double>>();
1480 *   system_mf_storage->reinit(
1481 *   mapping, dof_handler, dummy, QGauss<1>(fe.degree + 1), additional_data);
1482 *   system_matrix.initialize(system_mf_storage);
1483 *   }
1484 *  
1485 *   system_matrix.initialize_dof_vector(solution);
1486 *   system_matrix.initialize_dof_vector(system_rhs);
1487 *  
1488 *   setup_time += time.wall_time();
1489 *   time_details << "Setup matrix-free system " << time.wall_time() << " s"
1490 *   << std::endl;
1491 *   time.restart();
1492 *  
1493 *   const unsigned int nlevels = triangulation.n_global_levels();
1494 *   mg_matrices.resize(0, nlevels - 1);
1495 *  
1496 *   for (unsigned int level = 0; level < nlevels; ++level)
1497 *   {
1498 *   typename MatrixFree<dim, float>::AdditionalData additional_data;
1499 *   additional_data.tasks_parallel_scheme =
1500 *   MatrixFree<dim, float>::AdditionalData::none;
1501 *   additional_data.mapping_update_flags =
1502 *   (update_gradients | update_JxW_values);
1503 *   additional_data.mapping_update_flags_inner_faces =
1504 *   (update_gradients | update_JxW_values);
1505 *   additional_data.mapping_update_flags_boundary_faces =
1506 *   (update_gradients | update_JxW_values);
1507 *   additional_data.mg_level = level;
1508 *   const auto mg_mf_storage_level =
1509 *   std::make_shared<MatrixFree<dim, float>>();
1510 *   mg_mf_storage_level->reinit(mapping,
1511 *   dof_handler,
1512 *   dummy,
1513 *   QGauss<1>(fe.degree + 1),
1514 *   additional_data);
1515 *  
1516 *   mg_matrices[level].initialize(mg_mf_storage_level);
1517 *   }
1518 *   setup_time += time.wall_time();
1519 *   time_details << "Setup matrix-free levels " << time.wall_time() << " s"
1520 *   << std::endl;
1521 *   }
1522 *  
1523 *  
1524 *  
1525 * @endcode
1526 *
1527 * The computation of the right hand side is a bit more complicated than in
1528 * @ref step_37 "step-37". The cell term now consists of the negative Laplacian of the
1529 * analytical solution, `RightHandSide`, for which we need to first split up
1530 * the Point of VectorizedArray fields, i.e., a batch of points, into a
1531 * single point by evaluating all lanes in the VectorizedArray
1532 * separately. Remember that the number of lanes depends on the hardware; it
1533 * could be 1 for systems that do not offer vectorization (or where deal.II
1534 * does not have intrinsics), but it could also be 8 or 16 on AVX-512 of
1535 * recent Intel architectures.
1536 *
1537 * @code
1538 *   template <int dim, int fe_degree>
1539 *   void LaplaceProblem<dim, fe_degree>::compute_rhs()
1540 *   {
1541 *   Timer time;
1542 *   system_rhs = 0;
1543 *   const MatrixFree<dim, double> &data = *system_matrix.get_matrix_free();
1544 *   FEEvaluation<dim, fe_degree> phi(data);
1545 *   RightHandSide<dim> rhs_func;
1546 *   Solution<dim> exact_solution;
1547 *   for (unsigned int cell = 0; cell < data.n_cell_batches(); ++cell)
1548 *   {
1549 *   phi.reinit(cell);
1550 *   for (const unsigned int q : phi.quadrature_point_indices())
1551 *   {
1552 *   VectorizedArray<double> rhs_val = VectorizedArray<double>();
1553 *   Point<dim, VectorizedArray<double>> point_batch =
1554 *   phi.quadrature_point(q);
1555 *   for (unsigned int v = 0; v < VectorizedArray<double>::size(); ++v)
1556 *   {
1557 *   Point<dim> single_point;
1558 *   for (unsigned int d = 0; d < dim; ++d)
1559 *   single_point[d] = point_batch[d][v];
1560 *   rhs_val[v] = rhs_func.value(single_point);
1561 *   }
1562 *   phi.submit_value(rhs_val, q);
1563 *   }
1564 *   phi.integrate_scatter(EvaluationFlags::values, system_rhs);
1565 *   }
1566 *  
1567 * @endcode
1568 *
1569 * Secondly, we also need to apply the Dirichlet and Neumann boundary
1570 * conditions. This function is the missing part of to the function
1571 * `LaplaceOperator::apply_boundary()` function once the exterior solution
1572 * values @f$u^+ = -u^- + 2 g_\text{D}@f$ and @f$\mathbf{n}^-\cdot \nabla u^+ =
1573 * \mathbf{n}^-\cdot \nabla u^-@f$ on Dirichlet boundaries and @f$u^+=u^-@f$ and
1574 * @f$\mathbf{n}^-\cdot \nabla u^+ = -\mathbf{n}^-\cdot \nabla u^- + 2
1575 * g_\text{N}@f$ on Neumann boundaries are inserted and expanded in terms of
1576 * the boundary functions @f$g_\text{D}@f$ and @f$g_\text{N}@f$. One thing to
1577 * remember is that we move the boundary conditions to the right hand
1578 * side, so the sign is the opposite from what we imposed on the solution
1579 * part.
1580 *
1581
1582 *
1583 * We could have issued both the cell and the boundary part through a
1584 * MatrixFree::loop part, but we choose to manually write the full loop
1585 * over all faces to learn how the index layout of face indices is set up
1586 * in MatrixFree: Both the inner faces and the boundary faces share the
1587 * index range, and all batches of inner faces have lower numbers than the
1588 * batches of boundary cells. A single index for both variants allows us
1589 * to easily use the same data structure FEFaceEvaluation for both cases
1590 * that attaches to the same data field, just at different positions. The
1591 * number of inner face batches (where a batch is due to the combination
1592 * of several faces into one for vectorization) is given by
1593 * MatrixFree::n_inner_face_batches(), whereas the number of boundary face
1594 * batches is given by MatrixFree::n_boundary_face_batches().
1595 *
1596 * @code
1597 *   FEFaceEvaluation<dim, fe_degree> phi_face(data, true);
1598 *   for (unsigned int face = data.n_inner_face_batches();
1599 *   face < data.n_inner_face_batches() + data.n_boundary_face_batches();
1600 *   ++face)
1601 *   {
1602 *   phi_face.reinit(face);
1603 *  
1604 *   const VectorizedArray<double> inverse_length_normal_to_face = std::abs(
1605 *   (phi_face.normal_vector(0) * phi_face.inverse_jacobian(0))[dim - 1]);
1606 *   const VectorizedArray<double> sigma =
1607 *   inverse_length_normal_to_face * system_matrix.get_penalty_factor();
1608 *  
1609 *   for (const unsigned int q : phi_face.quadrature_point_indices())
1610 *   {
1611 *   VectorizedArray<double> test_value = VectorizedArray<double>(),
1612 *   test_normal_derivative =
1613 *   VectorizedArray<double>();
1614 *   Point<dim, VectorizedArray<double>> point_batch =
1615 *   phi_face.quadrature_point(q);
1616 *  
1617 *   for (unsigned int v = 0; v < VectorizedArray<double>::size(); ++v)
1618 *   {
1619 *   Point<dim> single_point;
1620 *   for (unsigned int d = 0; d < dim; ++d)
1621 *   single_point[d] = point_batch[d][v];
1622 *  
1623 * @endcode
1624 *
1625 * The MatrixFree class lets us query the boundary_id of the
1626 * current face batch. Remember that MatrixFree sets up the
1627 * batches for vectorization such that all faces within a
1628 * batch have the same properties, which includes their
1629 * `boundary_id`. Thus, we can query that id here for the
1630 * current face index `face` and either impose the Dirichlet
1631 * case (where we add something to the function value) or the
1632 * Neumann case (where we add something to the normal
1633 * derivative).
1634 *
1635 * @code
1636 *   if (data.get_boundary_id(face) == 0)
1637 *   test_value[v] = 2.0 * exact_solution.value(single_point);
1638 *   else
1639 *   {
1640 *   Tensor<1, dim> normal;
1641 *   for (unsigned int d = 0; d < dim; ++d)
1642 *   normal[d] = phi_face.normal_vector(q)[d][v];
1643 *   test_normal_derivative[v] =
1644 *   -normal * exact_solution.gradient(single_point);
1645 *   }
1646 *   }
1647 *   phi_face.submit_value(test_value * sigma - test_normal_derivative,
1648 *   q);
1649 *   phi_face.submit_normal_derivative(-0.5 * test_value, q);
1650 *   }
1651 *   phi_face.integrate_scatter(EvaluationFlags::values |
1652 *   EvaluationFlags::gradients,
1653 *   system_rhs);
1654 *   }
1655 *  
1656 * @endcode
1657 *
1658 * Since we have manually run the loop over cells rather than using
1659 * MatrixFree::loop(), we must not forget to perform the data exchange
1660 * with MPI - or actually, we would not need that for DG elements here
1661 * because each cell carries its own degrees of freedom and cell and
1662 * boundary integrals only evaluate quantities on the locally owned
1663 * cells. The coupling to neighboring subdomain only comes in by the inner
1664 * face integrals, which we have not done here. That said, it does not
1665 * hurt to call this function here, so we do it as a reminder of what
1666 * happens inside MatrixFree::loop().
1667 *
1668 * @code
1669 *   system_rhs.compress(VectorOperation::add);
1670 *   setup_time += time.wall_time();
1671 *   time_details << "Compute right hand side " << time.wall_time()
1672 *   << " s\n";
1673 *   }
1674 *  
1675 *  
1676 *  
1677 * @endcode
1678 *
1679 * The `solve()` function is copied almost verbatim from @ref step_37 "step-37". We set up
1680 * the same multigrid ingredients, namely the level transfer, a smoother,
1681 * and a coarse grid solver. The only difference is the fact that we do not
1682 * use the diagonal of the Laplacian for the preconditioner of the Chebyshev
1683 * iteration used for smoothing, but instead our newly resolved class
1684 * `%PreconditionBlockJacobi`. The mechanisms are the same, though.
1685 *
1686 * @code
1687 *   template <int dim, int fe_degree>
1688 *   void LaplaceProblem<dim, fe_degree>::solve()
1689 *   {
1690 *   Timer time;
1691 *   MGTransferMatrixFree<dim, float> mg_transfer;
1692 *   mg_transfer.build(dof_handler);
1693 *   setup_time += time.wall_time();
1694 *   time_details << "MG build transfer time " << time.wall_time()
1695 *   << " s\n";
1696 *   time.restart();
1697 *  
1698 *   using SmootherType =
1699 *   PreconditionChebyshev<LevelMatrixType,
1700 *   LinearAlgebra::distributed::Vector<float>,
1701 *   PreconditionBlockJacobi<dim, fe_degree, float>>;
1702 *   mg::SmootherRelaxation<SmootherType,
1703 *   LinearAlgebra::distributed::Vector<float>>
1704 *   mg_smoother;
1705 *   MGLevelObject<typename SmootherType::AdditionalData> smoother_data;
1706 *   smoother_data.resize(0, triangulation.n_global_levels() - 1);
1707 *   for (unsigned int level = 0; level < triangulation.n_global_levels();
1708 *   ++level)
1709 *   {
1710 *   if (level > 0)
1711 *   {
1712 *   smoother_data[level].smoothing_range = 15.;
1713 *   smoother_data[level].degree = 3;
1714 *   smoother_data[level].eig_cg_n_iterations = 10;
1715 *   }
1716 *   else
1717 *   {
1718 *   smoother_data[0].smoothing_range = 2e-2;
1719 *   smoother_data[0].degree = numbers::invalid_unsigned_int;
1720 *   smoother_data[0].eig_cg_n_iterations = mg_matrices[0].m();
1721 *   }
1722 *   smoother_data[level].preconditioner =
1723 *   std::make_shared<PreconditionBlockJacobi<dim, fe_degree, float>>();
1724 *   smoother_data[level].preconditioner->initialize(mg_matrices[level]);
1725 *   }
1726 *   mg_smoother.initialize(mg_matrices, smoother_data);
1727 *  
1728 *   MGCoarseGridApplySmoother<LinearAlgebra::distributed::Vector<float>>
1729 *   mg_coarse;
1730 *   mg_coarse.initialize(mg_smoother);
1731 *  
1732 *   mg::Matrix<LinearAlgebra::distributed::Vector<float>> mg_matrix(
1733 *   mg_matrices);
1734 *  
1735 *   Multigrid<LinearAlgebra::distributed::Vector<float>> mg(
1736 *   mg_matrix, mg_coarse, mg_transfer, mg_smoother, mg_smoother);
1737 *  
1738 *   PreconditionMG<dim,
1739 *   LinearAlgebra::distributed::Vector<float>,
1740 *   MGTransferMatrixFree<dim, float>>
1741 *   preconditioner(dof_handler, mg, mg_transfer);
1742 *  
1743 *   SolverControl solver_control(10000, 1e-12 * system_rhs.l2_norm());
1744 *   SolverCG<LinearAlgebra::distributed::Vector<double>> cg(solver_control);
1745 *   setup_time += time.wall_time();
1746 *   time_details << "MG build smoother time " << time.wall_time()
1747 *   << "s\n";
1748 *   pcout << "Total setup time " << setup_time << " s\n";
1749 *  
1750 *   time.reset();
1751 *   time.start();
1752 *   cg.solve(system_matrix, solution, system_rhs, preconditioner);
1753 *  
1754 *   pcout << "Time solve (" << solver_control.last_step() << " iterations) "
1755 *   << time.wall_time() << " s" << std::endl;
1756 *   }
1757 *  
1758 *  
1759 *  
1760 * @endcode
1761 *
1762 * Since we have solved a problem with analytic solution, we want to verify
1763 * the correctness of our implementation by computing the L2 error of the
1764 * numerical result against the analytic solution.
1765 *
1766
1767 *
1768 *
1769 * @code
1770 *   template <int dim, int fe_degree>
1771 *   void LaplaceProblem<dim, fe_degree>::analyze_results() const
1772 *   {
1773 *   Vector<float> error_per_cell(triangulation.n_active_cells());
1774 *   VectorTools::integrate_difference(mapping,
1775 *   dof_handler,
1776 *   solution,
1777 *   Solution<dim>(),
1778 *   error_per_cell,
1779 *   QGauss<dim>(fe.degree + 2),
1780 *   VectorTools::L2_norm);
1781 *   pcout << "Verification via L2 error: "
1782 *   << std::sqrt(
1783 *   Utilities::MPI::sum(error_per_cell.norm_sqr(), MPI_COMM_WORLD))
1784 *   << std::endl;
1785 *   }
1786 *  
1787 *  
1788 *  
1789 * @endcode
1790 *
1791 * The `run()` function sets up the initial grid and then runs the multigrid
1792 * program in the usual way. As a domain, we choose a rectangle with
1793 * periodic boundary conditions in the @f$x@f$-direction, a Dirichlet condition
1794 * on the front face in @f$y@f$ direction (i.e., the face with index number 2,
1795 * with boundary id equal to 0), and Neumann conditions on the back face as
1796 * well as the two faces in @f$z@f$ direction for the 3d case (with boundary id
1797 * equal to 1). The extent of the domain is a bit different in the @f$x@f$
1798 * direction (where we want to achieve a periodic solution given the
1799 * definition of `Solution`) as compared to the @f$y@f$ and @f$z@f$ directions.
1800 *
1801
1802 *
1803 *
1804 * @code
1805 *   template <int dim, int fe_degree>
1806 *   void LaplaceProblem<dim, fe_degree>::run()
1807 *   {
1808 *   const unsigned int n_ranks =
1809 *   Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD);
1810 *   pcout << "Running with " << n_ranks << " MPI process"
1811 *   << (n_ranks > 1 ? "es" : "") << ", element " << fe.get_name()
1812 *   << std::endl
1813 *   << std::endl;
1814 *   for (unsigned int cycle = 0; cycle < 9 - dim; ++cycle)
1815 *   {
1816 *   pcout << "Cycle " << cycle << std::endl;
1817 *  
1818 *   if (cycle == 0)
1819 *   {
1820 *   Point<dim> upper_right;
1821 *   upper_right[0] = 2.5;
1822 *   for (unsigned int d = 1; d < dim; ++d)
1823 *   upper_right[d] = 2.8;
1824 *   GridGenerator::hyper_rectangle(triangulation,
1825 *   Point<dim>(),
1826 *   upper_right);
1827 *   triangulation.begin_active()->face(0)->set_boundary_id(10);
1828 *   triangulation.begin_active()->face(1)->set_boundary_id(11);
1829 *   triangulation.begin_active()->face(2)->set_boundary_id(0);
1830 *   for (unsigned int f = 3;
1831 *   f < triangulation.begin_active()->n_faces();
1832 *   ++f)
1833 *   triangulation.begin_active()->face(f)->set_boundary_id(1);
1834 *  
1835 *   std::vector<GridTools::PeriodicFacePair<
1836 *   typename Triangulation<dim>::cell_iterator>>
1837 *   periodic_faces;
1838 *   GridTools::collect_periodic_faces(
1839 *   triangulation, 10, 11, 0, periodic_faces);
1840 *   triangulation.add_periodicity(periodic_faces);
1841 *  
1842 *   triangulation.refine_global(6 - 2 * dim);
1843 *   }
1844 *   triangulation.refine_global(1);
1845 *   setup_system();
1846 *   compute_rhs();
1847 *   solve();
1848 *   analyze_results();
1849 *   pcout << std::endl;
1850 *   };
1851 *   }
1852 *   } // namespace Step59
1853 *  
1854 *  
1855 *  
1856 * @endcode
1857 *
1858 * There is nothing unexpected in the `main()` function. We call `MPI_Init()`
1859 * through the `MPI_InitFinalize` class, pass on the two parameters on the
1860 * dimension and the degree set at the top of the file, and run the Laplace
1861 * problem.
1862 *
1863
1864 *
1865 *
1866 * @code
1867 *   int main(int argc, char *argv[])
1868 *   {
1869 *   try
1870 *   {
1871 *   using namespace Step59;
1872 *  
1873 *   Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1);
1874 *  
1875 *   LaplaceProblem<dimension, degree_finite_element> laplace_problem;
1876 *   laplace_problem.run();
1877 *   }
1878 *   catch (std::exception &exc)
1879 *   {
1880 *   std::cerr << std::endl
1881 *   << std::endl
1882 *   << "----------------------------------------------------"
1883 *   << std::endl;
1884 *   std::cerr << "Exception on processing: " << std::endl
1885 *   << exc.what() << std::endl
1886 *   << "Aborting!" << std::endl
1887 *   << "----------------------------------------------------"
1888 *   << std::endl;
1889 *   return 1;
1890 *   }
1891 *   catch (...)
1892 *   {
1893 *   std::cerr << std::endl
1894 *   << std::endl
1895 *   << "----------------------------------------------------"
1896 *   << std::endl;
1897 *   std::cerr << "Unknown exception!" << std::endl
1898 *   << "Aborting!" << std::endl
1899 *   << "----------------------------------------------------"
1900 *   << std::endl;
1901 *   return 1;
1902 *   }
1903 *  
1904 *   return 0;
1905 *   }
1906 * @endcode
1907<a name="step_59-Results"></a><h1>Results</h1>
1908
1909
1910<a name="step_59-Programoutput"></a><h3>Program output</h3>
1911
1912
1913Like in @ref step_37 "step-37", we evaluate the multigrid solver in terms of run time. In
1914two space dimensions with elements of degree 8, a possible output could look
1915as follows:
1916@code
1917Running with 12 MPI processes, element FE_DGQHermite<2>(8)
1918
1919Cycle 0
1920Number of degrees of freedom: 5184
1921Total setup time 0.0282445 s
1922Time solve (14 iterations) 0.0110712 s
1923Verification via L2 error: 1.66232e-07
1924
1925Cycle 1
1926Number of degrees of freedom: 20736
1927Total setup time 0.0126282 s
1928Time solve (14 iterations) 0.0157021 s
1929Verification via L2 error: 2.91505e-10
1930
1931Cycle 2
1932Number of degrees of freedom: 82944
1933Total setup time 0.0227573 s
1934Time solve (14 iterations) 0.026568 s
1935Verification via L2 error: 6.64514e-13
1936
1937Cycle 3
1938Number of degrees of freedom: 331776
1939Total setup time 0.0604685 s
1940Time solve (14 iterations) 0.0628356 s
1941Verification via L2 error: 5.57513e-13
1942
1943Cycle 4
1944Number of degrees of freedom: 1327104
1945Total setup time 0.154359 s
1946Time solve (13 iterations) 0.219555 s
1947Verification via L2 error: 3.08139e-12
1948
1949Cycle 5
1950Number of degrees of freedom: 5308416
1951Total setup time 0.467764 s
1952Time solve (13 iterations) 1.1821 s
1953Verification via L2 error: 3.90334e-12
1954
1955Cycle 6
1956Number of degrees of freedom: 21233664
1957Total setup time 1.73263 s
1958Time solve (13 iterations) 5.21054 s
1959Verification via L2 error: 4.94543e-12
1960@endcode
1961
1962Like in @ref step_37 "step-37", the number of CG iterations remains constant with increasing
1963problem size. The iteration counts are a bit higher, which is because we use a
1964lower degree of the Chebyshev polynomial (2 vs 5 in @ref step_37 "step-37") and because the
1965interior penalty discretization has a somewhat larger spread in
1966eigenvalues. Nonetheless, 13 iterations to reduce the residual by 12 orders of
1967magnitude, or almost a factor of 9 per iteration, indicates an overall very
1968efficient method. In particular, we can solve a system with 21 million degrees
1969of freedom in 5 seconds when using 12 cores, which is a very good
1970efficiency. Of course, in 2D we are well inside the regime of roundoff for a
1971polynomial degree of 8 &ndash; as a matter of fact, around 83k DoFs or 0.025s
1972would have been enough to fully converge this (simple) analytic solution
1973here.
1974
1975Not much changes if we run the program in three spatial dimensions, except for
1976the fact that we now use do something more useful with the higher polynomial
1977degree and increasing mesh sizes, as the roundoff errors are only obtained at
1978the finest mesh. Still, it is remarkable that we can solve a 3D Laplace
1979problem with a wave of three periods to roundoff accuracy on a twelve-core
1980machine pretty easily - using about 3.5 GB of memory in total for the second
1981to largest case with 24m DoFs, taking not more than eight seconds. The largest
1982case uses 30GB of memory with 191m DoFs.
1983
1984@code
1985Running with 12 MPI processes, element FE_DGQHermite<3>(8)
1986
1987Cycle 0
1988Number of degrees of freedom: 5832
1989Total setup time 0.0210681 s
1990Time solve (15 iterations) 0.0956945 s
1991Verification via L2 error: 0.0297194
1992
1993Cycle 1
1994Number of degrees of freedom: 46656
1995Total setup time 0.0452428 s
1996Time solve (15 iterations) 0.113827 s
1997Verification via L2 error: 9.55733e-05
1998
1999Cycle 2
2000Number of degrees of freedom: 373248
2001Total setup time 0.190423 s
2002Time solve (15 iterations) 0.218309 s
2003Verification via L2 error: 2.6868e-07
2004
2005Cycle 3
2006Number of degrees of freedom: 2985984
2007Total setup time 0.627914 s
2008Time solve (15 iterations) 1.0595 s
2009Verification via L2 error: 4.6918e-10
2010
2011Cycle 4
2012Number of degrees of freedom: 23887872
2013Total setup time 2.85215 s
2014Time solve (15 iterations) 8.30576 s
2015Verification via L2 error: 9.38583e-13
2016
2017Cycle 5
2018Number of degrees of freedom: 191102976
2019Total setup time 16.1324 s
2020Time solve (15 iterations) 65.57 s
2021Verification via L2 error: 3.17875e-13
2022@endcode
2023
2024<a name="step_59-Comparisonofefficiencyatdifferentpolynomialdegrees"></a><h3>Comparison of efficiency at different polynomial degrees</h3>
2025
2026
2027In the introduction and in-code comments, it was mentioned several times that
2028high orders are treated very efficiently with the FEEvaluation and
2029FEFaceEvaluation evaluators. Now, we want to substantiate these claims by
2030looking at the throughput of the 3D multigrid solver for various polynomial
2031degrees. We collect the times as follows: We first run a solver at problem
2032size close to ten million, indicated in the first four table rows, and record
2033the timings. Then, we normalize the throughput by recording the number of
2034million degrees of freedom solved per second (MDoFs/s) to be able to compare
2035the efficiency of the different degrees, which is computed by dividing the
2036number of degrees of freedom by the solver time.
2037
2038<table align="center" class="doxtable">
2039 <tr>
2040 <th>degree</th>
2041 <th>1</th>
2042 <th>2</th>
2043 <th>3</th>
2044 <th>4</th>
2045 <th>5</th>
2046 <th>6</th>
2047 <th>7</th>
2048 <th>8</th>
2049 <th>9</th>
2050 <th>10</th>
2051 <th>11</th>
2052 <th>12</th>
2053 </tr>
2054 <tr>
2055 <th>Number of DoFs</th>
2056 <td>2097152</td>
2057 <td>7077888</td>
2058 <td>16777216</td>
2059 <td>32768000</td>
2060 <td>7077888</td>
2061 <td>11239424</td>
2062 <td>16777216</td>
2063 <td>23887872</td>
2064 <td>32768000</td>
2065 <td>43614208</td>
2066 <td>7077888</td>
2067 <td>8998912</td>
2068 </tr>
2069 <tr>
2070 <th>Number of iterations</th>
2071 <td>13</td>
2072 <td>12</td>
2073 <td>12</td>
2074 <td>12</td>
2075 <td>13</td>
2076 <td>13</td>
2077 <td>15</td>
2078 <td>15</td>
2079 <td>17</td>
2080 <td>19</td>
2081 <td>18</td>
2082 <td>18</td>
2083 </tr>
2084 <tr>
2085 <th>Solver time [s]</th>
2086 <td>0.713</td>
2087 <td>2.150</td>
2088 <td>4.638</td>
2089 <td>8.803</td>
2090 <td>2.041</td>
2091 <td>3.295</td>
2092 <td>5.723</td>
2093 <td>8.306</td>
2094 <td>12.75</td>
2095 <td>19.25</td>
2096 <td>3.530</td>
2097 <td>4.814</td>
2098 </tr>
2099 <tr>
2100 <th>MDoFs/s</th>
2101 <td>2.94</td>
2102 <td>3.29</td>
2103 <td>3.62</td>
2104 <td>3.72</td>
2105 <td>3.47</td>
2106 <td>3.41</td>
2107 <td>2.93</td>
2108 <td>2.88</td>
2109 <td>2.57</td>
2110 <td>2.27</td>
2111 <td>2.01</td>
2112 <td>1.87</td>
2113 </tr>
2114</table>
2115
2116We clearly see how the efficiency per DoF initially improves until it reaches
2117a maximum for the polynomial degree @f$k=4@f$. This effect is surprising, not only
2118because higher polynomial degrees often yield a vastly better solution, but
2119especially also when having matrix-based schemes in mind where the denser
2120coupling at higher degree leads to a monotonously decreasing throughput (and a
2121drastic one in 3D, with @f$k=4@f$ being more than ten times slower than
2122@f$k=1@f$!). For higher degrees, the throughput decreases a bit, which is both due
2123to an increase in the number of iterations (going from 12 at @f$k=2,3,4@f$ to 19
2124at @f$k=10@f$) and due to the @f$\mathcal O(k)@f$ complexity of operator
2125evaluation. Nonetheless, efficiency as the time to solution would be still
2126better for higher polynomial degrees because they have better convergence rates (at least
2127for problems as simple as this one): For @f$k=12@f$, we reach roundoff accuracy
2128already with 1 million DoFs (solver time less than a second), whereas for @f$k=8@f$
2129we need 24 million DoFs and 8 seconds. For @f$k=5@f$, the error is around
2130@f$10^{-9}@f$ with 57m DoFs and thus still far away from roundoff, despite taking 16
2131seconds.
2132
2133Note that the above numbers are a bit pessimistic because they include the
2134time it takes the Chebyshev smoother to compute an eigenvalue estimate, which
2135is around 10 percent of the solver time. If the system is solved several times
2136(as e.g. common in fluid dynamics), this eigenvalue cost is only paid once and
2137faster times become available.
2138
2139<a name="step_59-Evaluationofefficiencyofingredients"></a><h3>Evaluation of efficiency of ingredients</h3>
2140
2141
2142Finally, we take a look at some of the special ingredients presented in this
2143tutorial program, namely the FE_DGQHermite basis in particular and the
2144specification of MatrixFree::DataAccessOnFaces. In the following table, the
2145third row shows the optimized solver above, the fourth row shows the timings
2146with only the MatrixFree::DataAccessOnFaces set to `unspecified` rather than
2147the optimal `gradients`, and the last one with replacing FE_DGQHermite by the
2148basic FE_DGQ elements where both the MPI exchange are more expensive and the
2149operations done by FEFaceEvaluation::gather_evaluate() and
2150FEFaceEvaluation::integrate_scatter().
2151
2152<table align="center" class="doxtable">
2153 <tr>
2154 <th>degree</th>
2155 <th>1</th>
2156 <th>2</th>
2157 <th>3</th>
2158 <th>4</th>
2159 <th>5</th>
2160 <th>6</th>
2161 <th>7</th>
2162 <th>8</th>
2163 <th>9</th>
2164 <th>10</th>
2165 <th>11</th>
2166 <th>12</th>
2167 </tr>
2168 <tr>
2169 <th>Number of DoFs</th>
2170 <td>2097152</td>
2171 <td>7077888</td>
2172 <td>16777216</td>
2173 <td>32768000</td>
2174 <td>7077888</td>
2175 <td>11239424</td>
2176 <td>16777216</td>
2177 <td>23887872</td>
2178 <td>32768000</td>
2179 <td>43614208</td>
2180 <td>7077888</td>
2181 <td>8998912</td>
2182 </tr>
2183 <tr>
2184 <th>Solver time optimized as in tutorial [s]</th>
2185 <td>0.713</td>
2186 <td>2.150</td>
2187 <td>4.638</td>
2188 <td>8.803</td>
2189 <td>2.041</td>
2190 <td>3.295</td>
2191 <td>5.723</td>
2192 <td>8.306</td>
2193 <td>12.75</td>
2194 <td>19.25</td>
2195 <td>3.530</td>
2196 <td>4.814</td>
2197 </tr>
2198 <tr>
2199 <th>Solver time MatrixFree::DataAccessOnFaces::unspecified [s]</th>
2200 <td>0.711</td>
2201 <td>2.151</td>
2202 <td>4.675</td>
2203 <td>8.968</td>
2204 <td>2.243</td>
2205 <td>3.655</td>
2206 <td>6.277</td>
2207 <td>9.082</td>
2208 <td>13.50</td>
2209 <td>20.05</td>
2210 <td>3.817</td>
2211 <td>5.178</td>
2212 </tr>
2213 <tr>
2214 <th>Solver time FE_DGQ [s]</th>
2215 <td>0.712</td>
2216 <td>2.041</td>
2217 <td>5.066</td>
2218 <td>9.335</td>
2219 <td>2.379</td>
2220 <td>3.802</td>
2221 <td>6.564</td>
2222 <td>9.714</td>
2223 <td>14.54</td>
2224 <td>22.76</td>
2225 <td>4.148</td>
2226 <td>5.857</td>
2227 </tr>
2228</table>
2229
2230The data in the table shows that not using MatrixFree::DataAccessOnFaces
2231increases costs by around 10% for higher polynomial degrees. For lower
2232degrees, the difference is obviously less pronounced because the
2233volume-to-surface ratio is more beneficial and less data needs to be
2234exchanged. The difference is larger when looking at the matrix-vector product
2235only, rather than the full multigrid solver shown here, with around 20% worse
2236timings just because of the MPI communication.
2237
2238For @f$k=1@f$ and @f$k=2@f$, the Hermite-like basis functions do obviously not really
2239pay off (indeed, for @f$k=1@f$ the polynomials are exactly the same as for FE_DGQ)
2240and the results are similar as with the FE_DGQ basis. However, for degrees
2241starting at three, we see an increasing advantage for FE_DGQHermite, showing
2242the effectiveness of these basis functions.
2243
2244<a name="step_59-Possibilitiesforextension"></a><h3>Possibilities for extension</h3>
2245
2246
2247As mentioned in the introduction, the fast diagonalization method as realized
2248here is tied to a Cartesian mesh with constant coefficients. When dealing with
2249meshes that contain deformed cells or with variable coefficients, it is common
2250to determine a nearby Cartesian mesh cell as an approximation. This can be
2251done with the class TensorProductMatrixSymmetricSumCollection. Here, one can
2252insert cell matrices similarly to the PreconditionBlockJacobi::initialize()
2253function of this tutorial program. The benefit of the collection class is that
2254cells on which the coefficient of the PDE has the same value can re-use the
2255same Laplacian matrix, which reduces the memory consumption for the inverse
2256matrices. As compared to the algorithm implemented in this tutorial program,
2257one would define the length scales as the distances between opposing
2258faces. For continuous elements, the code project <a
2259href=https://github.com/peterrum/dealii-dd-and-schwarz">Cache-optimized and
2260low-overhead implementations of multigrid smoothers for high-order FEM
2261computations</a> presents the computation for continuous elements. There is
2262currently no infrastructure in deal.II to automatically generate the 1D
2263matrices for discontinuous elements with SIP-DG discretization, as opposed to
2264continuous elements, where we provide
2265TensorProductMatrixCreator::create_laplace_tensor_product_matrix().
2266
2267Another way of extending the program would be to include support for adaptive
2268meshes. While the classical approach of defining interface operations at edges
2269of different refinement level, as discussed in @ref step_39 "step-39", is one possibility,
2270for Poisson-type problems another option is typically more beneficial. Using
2271the class MGTransferGlobalCoarsening, which is explained in the @ref step_75 "step-75"
2272tutorial program, one can deal with meshes of hanging nodes on all levels. An
2273algorithmic improvement can be obtained by combining the discontinuous
2274function space with the auxiliary continuous finite element space of the same
2275polynomial degree. This idea, introduced by Antonietti et al.
2276@cite antonietti2016uniform in 2016, allows making the multigrid convergence
2277independent of the penalty parameter. As demonstrated by Fehn et al.
2278@cite fehn2020hybrid, this also gives considerably lower iteration counts than
2279a multigrid solver directly working on levels with discontinuous function
2280spaces. The latter work also proposes p-multigrid techniques and combination
2281with algebraic multigrid coarse spaces as a means to efficiently solve Poisson
2282problems with high-order discontinuous Galerkin discretizations on complicated
2283geometries, representing the current state-of-the-art for simple Poisson-type
2284problems. The class MGTransferGlobalCoarsening provides features for each of
2285these three coarsening variants, the discontinuous-continuous auxiliary
2286function concept, p-multigrid, and traditional h-multigrid. The main
2287ingredient is to define an appropriate MGTwoLevelTransfer object and call
2288MGTwoLevelTransfer::reinit_geometric_transfer() or
2289MGTwoLevelTransfer::reinit_polynomial_transfer(), respectively.
2290 *
2291 *
2292<a name="step_59-PlainProg"></a>
2293<h1> The plain program</h1>
2294@include "step-59.cc"
2295*/
gradient_type get_gradient(const unsigned int q_point) const
void submit_gradient(const gradient_type grad_in, const unsigned int q_point)
void read_dof_values(const VectorType &src, const unsigned int first_index=0, const std::bitset< n_lanes > &mask=std::bitset< n_lanes >().flip())
std_cxx20::ranges::iota_view< unsigned int, unsigned int > quadrature_point_indices() const
Point< 2 > second
Definition grid_out.cc:4624
Point< 2 > first
Definition grid_out.cc:4623
__global__ void reduction(Number *result, const Number *v, const size_type N)
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(DOFINFO &, DOFINFO &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
Definition loop.h:442
Expression sign(const Expression &x)
void random(DoFHandler< dim, spacedim > &dof_handler)
spacedim & mesh
Definition grid_tools.h:989
@ matrix
Contents is actually a matrix.
@ symmetric
Matrix is symmetric.
@ general
No special properties.
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:191
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
T sum(const T &t, const MPI_Comm mpi_communicator)
T reduce(const T &local_value, const MPI_Comm comm, const std::function< T(const T &, const T &)> &combiner, const unsigned int root_process=0)
T scatter(const MPI_Comm comm, const std::vector< T > &objects_to_send, const unsigned int root_process=0)
int(&) functions(const void *v1, const void *v2)
const InputIterator OutputIterator out
Definition parallel.h:167
const InputIterator OutputIterator const Function & function
Definition parallel.h:168
STL namespace.
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)