Reference documentation for deal.II version 9.6.0
\(\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-9.h
Go to the documentation of this file.
1
1067false);
1068 *   sparsity_pattern.copy_from(dsp);
1069 *  
1070 *   system_matrix.reinit(sparsity_pattern);
1071 *  
1072 *   solution.reinit(dof_handler.n_dofs());
1073 *   system_rhs.reinit(dof_handler.n_dofs());
1074 *   }
1075 *  
1076 *  
1077 *  
1078 * @endcode
1079 *
1080 * In the following function, the matrix and right hand side are
1081 * assembled. As stated in the documentation of the main class above, it
1082 * does not do this itself, but rather delegates to the function following
1083 * next, utilizing the WorkStream concept discussed in @ref threads .
1084 *
1085
1086 *
1087 * If you have looked through the @ref threads topic, you will have
1088 * seen that assembling in parallel does not take an incredible
1089 * amount of extra code as long as you diligently describe what the
1090 * scratch and copy data objects are, and if you define suitable
1091 * functions for the local assembly and the copy operation from local
1092 * contributions to global objects. This done, the following will do
1093 * all the heavy lifting to get these operations done on multiple
1094 * threads on as many cores as you have in your system:
1095 *
1096 * @code
1097 *   template <int dim>
1098 *   void AdvectionProblem<dim>::assemble_system()
1099 *   {
1100 *   WorkStream::run(dof_handler.begin_active(),
1101 *   dof_handler.end(),
1102 *   *this,
1103 *   &AdvectionProblem::local_assemble_system,
1104 *   &AdvectionProblem::copy_local_to_global,
1105 *   AssemblyScratchData(fe),
1106 *   AssemblyCopyData());
1107 *   }
1108 *  
1109 *  
1110 *  
1111 * @endcode
1112 *
1113 * As already mentioned above, we need to have scratch objects for
1114 * the parallel computation of local contributions. These objects
1115 * contain FEValues and FEFaceValues objects (as well as some arrays), and so
1116 * we will need to have constructors and copy constructors that allow us to
1117 * create them. For the cell terms we need the values
1118 * and gradients of the shape functions, the quadrature points in
1119 * order to determine the source density and the advection field at
1120 * a given point, and the weights of the quadrature points times the
1121 * determinant of the Jacobian at these points. In contrast, for the
1122 * boundary integrals, we don't need the gradients, but rather the
1123 * normal vectors to the cells. This determines which update flags
1124 * we will have to pass to the constructors of the members of the
1125 * class:
1126 *
1127 * @code
1128 *   template <int dim>
1129 *   AdvectionProblem<dim>::AssemblyScratchData::AssemblyScratchData(
1130 *   const FiniteElement<dim> &fe)
1131 *   : fe_values(fe,
1132 *   QGauss<dim>(fe.degree + 1),
1134 *   update_JxW_values)
1135 *   , fe_face_values(fe,
1136 *   QGauss<dim - 1>(fe.degree + 1),
1139 *   , rhs_values(fe_values.get_quadrature().size())
1140 *   , advection_directions(fe_values.get_quadrature().size())
1141 *   , face_boundary_values(fe_face_values.get_quadrature().size())
1142 *   , face_advection_directions(fe_face_values.get_quadrature().size())
1143 *   {}
1144 *  
1145 *  
1146 *  
1147 *   template <int dim>
1148 *   AdvectionProblem<dim>::AssemblyScratchData::AssemblyScratchData(
1149 *   const AssemblyScratchData &scratch_data)
1150 *   : fe_values(scratch_data.fe_values.get_fe(),
1151 *   scratch_data.fe_values.get_quadrature(),
1153 *   update_JxW_values)
1154 *   , fe_face_values(scratch_data.fe_face_values.get_fe(),
1155 *   scratch_data.fe_face_values.get_quadrature(),
1158 *   , rhs_values(scratch_data.rhs_values.size())
1159 *   , advection_directions(scratch_data.advection_directions.size())
1160 *   , face_boundary_values(scratch_data.face_boundary_values.size())
1161 *   , face_advection_directions(scratch_data.face_advection_directions.size())
1162 *   {}
1163 *  
1164 *  
1165 *  
1166 * @endcode
1167 *
1168 * Now, this is the function that does the actual work. It is not very
1169 * different from the <code>assemble_system</code> functions of previous
1170 * example programs, so we will again only comment on the differences. The
1171 * mathematical stuff closely follows what we have said in the introduction.
1172 *
1173
1174 *
1175 * There are a number of points worth mentioning here, though. The
1176 * first one is that we have moved the FEValues and FEFaceValues
1177 * objects into the ScratchData object. We have done so because the
1178 * alternative would have been to simply create one every time we
1179 * get into this function -- i.e., on every cell. It now turns out
1180 * that the FEValues classes were written with the explicit goal of
1181 * moving everything that remains the same from cell to cell into
1182 * the construction of the object, and only do as little work as
1183 * possible in FEValues::reinit() whenever we move to a new
1184 * cell. What this means is that it would be very expensive to
1185 * create a new object of this kind in this function as we would
1186 * have to do it for every cell -- exactly the thing we wanted to
1187 * avoid with the FEValues class. Instead, what we do is create it
1188 * only once (or a small number of times) in the scratch objects and
1189 * then re-use it as often as we can.
1190 *
1191
1192 *
1193 * This begs the question of whether there are other objects we
1194 * create in this function whose creation is expensive compared to
1195 * its use. Indeed, at the top of the function, we declare all sorts
1196 * of objects. The <code>AdvectionField</code>,
1197 * <code>RightHandSide</code> and <code>BoundaryValues</code> do not
1198 * cost much to create, so there is no harm here. However,
1199 * allocating memory in creating the <code>rhs_values</code> and
1200 * similar variables below typically costs a significant amount of
1201 * time, compared to just accessing the (temporary) values we store
1202 * in them. Consequently, these would be candidates for moving into
1203 * the <code>AssemblyScratchData</code> class. We will leave this as
1204 * an exercise.
1205 *
1206 * @code
1207 *   template <int dim>
1208 *   void AdvectionProblem<dim>::local_assemble_system(
1209 *   const typename DoFHandler<dim>::active_cell_iterator &cell,
1210 *   AssemblyScratchData &scratch_data,
1211 *   AssemblyCopyData &copy_data)
1212 *   {
1213 * @endcode
1214 *
1215 * We define some abbreviations to avoid unnecessarily long lines:
1216 *
1217 * @code
1218 *   const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
1219 *   const unsigned int n_q_points =
1220 *   scratch_data.fe_values.get_quadrature().size();
1221 *   const unsigned int n_face_q_points =
1222 *   scratch_data.fe_face_values.get_quadrature().size();
1223 *  
1224 * @endcode
1225 *
1226 * We declare cell matrix and cell right hand side...
1227 *
1228 * @code
1229 *   copy_data.cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
1230 *   copy_data.cell_rhs.reinit(dofs_per_cell);
1231 *  
1232 * @endcode
1233 *
1234 * ... an array to hold the global indices of the degrees of freedom of
1235 * the cell on which we are presently working...
1236 *
1237 * @code
1238 *   copy_data.local_dof_indices.resize(dofs_per_cell);
1239 *  
1240 * @endcode
1241 *
1242 * ... then initialize the <code>FEValues</code> object...
1243 *
1244 * @code
1245 *   scratch_data.fe_values.reinit(cell);
1246 *  
1247 * @endcode
1248 *
1249 * ... obtain the values of right hand side and advection directions
1250 * at the quadrature points...
1251 *
1252 * @code
1253 *   scratch_data.advection_field.value_list(
1254 *   scratch_data.fe_values.get_quadrature_points(),
1255 *   scratch_data.advection_directions);
1256 *   scratch_data.right_hand_side.value_list(
1257 *   scratch_data.fe_values.get_quadrature_points(), scratch_data.rhs_values);
1258 *  
1259 * @endcode
1260 *
1261 * ... set the value of the streamline diffusion parameter as
1262 * described in the introduction...
1263 *
1264 * @code
1265 *   const double delta = 0.1 * cell->diameter();
1266 *  
1267 * @endcode
1268 *
1269 * ... and assemble the local contributions to the system matrix and
1270 * right hand side as also discussed above:
1271 *
1272 * @code
1273 *   for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
1274 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
1275 *   {
1276 * @endcode
1277 *
1278 * Alias the AssemblyScratchData object to keep the lines from
1279 * getting too long:
1280 *
1281 * @code
1282 *   const auto &sd = scratch_data;
1283 *   for (unsigned int j = 0; j < dofs_per_cell; ++j)
1284 *   copy_data.cell_matrix(i, j) +=
1285 *   ((sd.fe_values.shape_value(i, q_point) + // (phi_i +
1286 *   delta * (sd.advection_directions[q_point] * // delta beta
1287 *   sd.fe_values.shape_grad(i, q_point))) * // grad phi_i)
1288 *   sd.advection_directions[q_point] * // beta
1289 *   sd.fe_values.shape_grad(j, q_point)) * // grad phi_j
1290 *   sd.fe_values.JxW(q_point); // dx
1291 *  
1292 *   copy_data.cell_rhs(i) +=
1293 *   (sd.fe_values.shape_value(i, q_point) + // (phi_i +
1294 *   delta * (sd.advection_directions[q_point] * // delta beta
1295 *   sd.fe_values.shape_grad(i, q_point))) * // grad phi_i)
1296 *   sd.rhs_values[q_point] * // f
1297 *   sd.fe_values.JxW(q_point); // dx
1298 *   }
1299 *  
1300 * @endcode
1301 *
1302 * Besides the cell terms which we have built up now, the bilinear
1303 * form of the present problem also contains terms on the boundary of
1304 * the domain. Therefore, we have to check whether any of the faces of
1305 * this cell are on the boundary of the domain, and if so assemble the
1306 * contributions of this face as well. Of course, the bilinear form
1307 * only contains contributions from the <code>inflow</code> part of
1308 * the boundary, but to find out whether a certain part of a face of
1309 * the present cell is part of the inflow boundary, we have to have
1310 * information on the exact location of the quadrature points and on
1311 * the direction of flow at this point; we obtain this information
1312 * using the FEFaceValues object and only decide within the main loop
1313 * whether a quadrature point is on the inflow boundary.
1314 *
1315 * @code
1316 *   for (const auto &face : cell->face_iterators())
1317 *   if (face->at_boundary())
1318 *   {
1319 * @endcode
1320 *
1321 * Ok, this face of the present cell is on the boundary of the
1322 * domain. Just as for the usual FEValues object which we have
1323 * used in previous examples and also above, we have to
1324 * reinitialize the FEFaceValues object for the present face:
1325 *
1326 * @code
1327 *   scratch_data.fe_face_values.reinit(cell, face);
1328 *  
1329 * @endcode
1330 *
1331 * For the quadrature points at hand, we ask for the values of
1332 * the inflow function and for the direction of flow:
1333 *
1334 * @code
1335 *   scratch_data.boundary_values.value_list(
1336 *   scratch_data.fe_face_values.get_quadrature_points(),
1337 *   scratch_data.face_boundary_values);
1338 *   scratch_data.advection_field.value_list(
1339 *   scratch_data.fe_face_values.get_quadrature_points(),
1340 *   scratch_data.face_advection_directions);
1341 *  
1342 * @endcode
1343 *
1344 * Now loop over all quadrature points and see whether this face is on
1345 * the inflow or outflow part of the boundary. The normal
1346 * vector points out of the cell: since the face is at
1347 * the boundary, the normal vector points out of the domain,
1348 * so if the advection direction points into the domain, its
1349 * scalar product with the normal vector must be negative (to see why
1350 * this is true, consider the scalar product definition that uses a
1351 * cosine):
1352 *
1353 * @code
1354 *   for (unsigned int q_point = 0; q_point < n_face_q_points; ++q_point)
1355 *   if (scratch_data.fe_face_values.normal_vector(q_point) *
1356 *   scratch_data.face_advection_directions[q_point] <
1357 *   0.)
1358 * @endcode
1359 *
1360 * If the face is part of the inflow boundary, then compute the
1361 * contributions of this face to the global matrix and right
1362 * hand side, using the values obtained from the
1363 * FEFaceValues object and the formulae discussed in the
1364 * introduction:
1365 *
1366 * @code
1367 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
1368 *   {
1369 *   for (unsigned int j = 0; j < dofs_per_cell; ++j)
1370 *   copy_data.cell_matrix(i, j) -=
1371 *   (scratch_data.face_advection_directions[q_point] *
1372 *   scratch_data.fe_face_values.normal_vector(q_point) *
1373 *   scratch_data.fe_face_values.shape_value(i, q_point) *
1374 *   scratch_data.fe_face_values.shape_value(j, q_point) *
1375 *   scratch_data.fe_face_values.JxW(q_point));
1376 *  
1377 *   copy_data.cell_rhs(i) -=
1378 *   (scratch_data.face_advection_directions[q_point] *
1379 *   scratch_data.fe_face_values.normal_vector(q_point) *
1380 *   scratch_data.face_boundary_values[q_point] *
1381 *   scratch_data.fe_face_values.shape_value(i, q_point) *
1382 *   scratch_data.fe_face_values.JxW(q_point));
1383 *   }
1384 *   }
1385 *  
1386 * @endcode
1387 *
1388 * The final piece of information the copy routine needs is the global
1389 * indices of the degrees of freedom on this cell, so we end by writing
1390 * them to the local array:
1391 *
1392 * @code
1393 *   cell->get_dof_indices(copy_data.local_dof_indices);
1394 *   }
1395 *  
1396 *  
1397 *  
1398 * @endcode
1399 *
1400 * The second function we needed to write was the one that copies
1401 * the local contributions the previous function computed (and
1402 * put into the AssemblyCopyData object) into the global matrix and right
1403 * hand side vector objects. This is essentially what we always had
1404 * as the last block of code when assembling something on every
1405 * cell. The following should therefore be pretty obvious:
1406 *
1407 * @code
1408 *   template <int dim>
1409 *   void
1410 *   AdvectionProblem<dim>::copy_local_to_global(const AssemblyCopyData &copy_data)
1411 *   {
1412 *   hanging_node_constraints.distribute_local_to_global(
1413 *   copy_data.cell_matrix,
1414 *   copy_data.cell_rhs,
1415 *   copy_data.local_dof_indices,
1416 *   system_matrix,
1417 *   system_rhs);
1418 *   }
1419 *  
1420 * @endcode
1421 *
1422 * The next function is the linear solver routine. As the system is no longer
1423 * symmetric positive definite as in all the previous examples, we cannot
1424 * use the Conjugate Gradient method any more. Rather, we use a solver that
1425 * is more general and does not rely on any special properties of the
1426 * matrix: the GMRES method. GMRES, like the conjugate gradient method,
1427 * requires a decent preconditioner: we use a Jacobi preconditioner here,
1428 * which works well enough for this problem.
1429 *
1430 * @code
1431 *   template <int dim>
1432 *   void AdvectionProblem<dim>::solve()
1433 *   {
1434 *   SolverControl solver_control(std::max<std::size_t>(1000,
1435 *   system_rhs.size() / 10),
1436 *   1e-10 * system_rhs.l2_norm());
1437 *   SolverGMRES<Vector<double>> solver(solver_control);
1438 *   PreconditionJacobi<SparseMatrix<double>> preconditioner;
1439 *   preconditioner.initialize(system_matrix, 1.0);
1440 *   solver.solve(system_matrix, solution, system_rhs, preconditioner);
1441 *  
1442 *   hanging_node_constraints.distribute(solution);
1443 *  
1444 *   std::cout << " Iterations required for convergence: "
1445 *   << solver_control.last_step() << '\n'
1446 *   << " Norm of residual at convergence: "
1447 *   << solver_control.last_value() << '\n';
1448 *   }
1449 *  
1450 * @endcode
1451 *
1452 * The following function refines the grid according to the quantity
1453 * described in the introduction. The respective computations are made in
1454 * the class <code>GradientEstimation</code>.
1455 *
1456 * @code
1457 *   template <int dim>
1458 *   void AdvectionProblem<dim>::refine_grid()
1459 *   {
1460 *   Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
1461 *  
1462 *   GradientEstimation::estimate(dof_handler,
1463 *   solution,
1464 *   estimated_error_per_cell);
1465 *  
1467 *   estimated_error_per_cell,
1468 *   0.3,
1469 *   0.03);
1470 *  
1472 *   }
1473 *  
1474 * @endcode
1475 *
1476 * This function is similar to the one in @ref step_6 "step-6", but since we use a higher
1477 * degree finite element we save the solution in a different
1478 * way. Visualization programs like VisIt and Paraview typically only
1479 * understand data that is associated with nodes: they cannot plot
1480 * fifth-degree basis functions, which results in a very inaccurate picture
1481 * of the solution we computed. To get around this we save multiple
1482 * <em>patches</em> per cell: in 2d we save @f$8\times 8=64@f$ bilinear
1483 * `sub-cells' to the VTU file for each cell, and in 3d we save
1484 * @f$8\times 8\times 8 = 512@f$. The end result is that the
1485 * visualization program will use a piecewise linear interpolation of the
1486 * cubic basis functions on a 3 times refined mesh:
1487 * This captures the solution detail and, with most
1488 * screen resolutions, looks smooth. We save the grid in a separate step
1489 * with no extra patches so that we have a visual representation of the cell
1490 * faces.
1491 *
1492
1493 *
1494 * @note Version 9.1 of deal.II gained the ability to write higher degree
1495 * polynomials (i.e., write piecewise bicubic visualization data for our
1496 * piecewise bicubic solution) VTK and VTU output: however, not all recent
1497 * versions of ParaView and VisIt (as of 2018) can read this format, so we
1498 * use the older, more general (but less efficient) approach here.
1499 *
1500 * @code
1501 *   template <int dim>
1502 *   void AdvectionProblem<dim>::output_results(const unsigned int cycle) const
1503 *   {
1504 *   {
1505 *   GridOut grid_out;
1506 *   const std::string filename = "grid-" + std::to_string(cycle) + ".vtu";
1507 *   std::ofstream output(filename);
1508 *   grid_out.write_vtu(triangulation, output);
1509 *   std::cout << " Grid written to " << filename << std::endl;
1510 *   }
1511 *  
1512 *   {
1513 *   DataOut<dim> data_out;
1514 *   data_out.attach_dof_handler(dof_handler);
1515 *   data_out.add_data_vector(solution, "solution");
1516 *   data_out.build_patches(8);
1517 *  
1518 * @endcode
1519 *
1520 * VTU output can be expensive, both to compute and to write to
1521 * disk. Here we ask ZLib, a compression library, to compress the data
1522 * in a way that maximizes throughput.
1523 *
1524 * @code
1525 *   DataOutBase::VtkFlags vtk_flags;
1526 *   vtk_flags.compression_level = DataOutBase::CompressionLevel::best_speed;
1527 *   data_out.set_flags(vtk_flags);
1528 *  
1529 *   const std::string filename = "solution-" + std::to_string(cycle) + ".vtu";
1530 *   std::ofstream output(filename);
1531 *   data_out.write_vtu(output);
1532 *   std::cout << " Solution written to " << filename << std::endl;
1533 *   }
1534 *   }
1535 *  
1536 *  
1537 * @endcode
1538 *
1539 * ... as is the main loop (setup -- solve -- refine), aside from the number
1540 * of cycles and the initial grid:
1541 *
1542 * @code
1543 *   template <int dim>
1544 *   void AdvectionProblem<dim>::run()
1545 *   {
1546 *   for (unsigned int cycle = 0; cycle < 10; ++cycle)
1547 *   {
1548 *   std::cout << "Cycle " << cycle << ':' << std::endl;
1549 *  
1550 *   if (cycle == 0)
1551 *   {
1552 *   GridGenerator::hyper_cube(triangulation, -1, 1);
1553 *   triangulation.refine_global(3);
1554 *   }
1555 *   else
1556 *   {
1557 *   refine_grid();
1558 *   }
1559 *  
1560 *  
1561 *   std::cout << " Number of active cells: "
1562 *   << triangulation.n_active_cells() << std::endl;
1563 *  
1564 *   setup_system();
1565 *  
1566 *   std::cout << " Number of degrees of freedom: "
1567 *   << dof_handler.n_dofs() << std::endl;
1568 *  
1569 *   assemble_system();
1570 *   solve();
1571 *   output_results(cycle);
1572 *   }
1573 *   }
1574 *  
1575 *  
1576 *  
1577 * @endcode
1578 *
1579 *
1580 * <a name="step_9-GradientEstimationclassimplementation"></a>
1581 * <h3>GradientEstimation class implementation</h3>
1582 *
1583
1584 *
1585 * Now for the implementation of the <code>GradientEstimation</code> class.
1586 * Let us start by defining constructors for the
1587 * <code>EstimateScratchData</code> class used by the
1588 * <code>estimate_cell()</code> function:
1589 *
1590 * @code
1591 *   template <int dim>
1592 *   GradientEstimation::EstimateScratchData<dim>::EstimateScratchData(
1593 *   const FiniteElement<dim> &fe,
1594 *   const Vector<double> &solution,
1595 *   Vector<float> &error_per_cell)
1596 *   : fe_midpoint_value(fe,
1597 *   QMidpoint<dim>(),
1598 *   update_values | update_quadrature_points)
1599 *   , solution(solution)
1600 *   , error_per_cell(error_per_cell)
1601 *   , cell_midpoint_value(1)
1602 *   , neighbor_midpoint_value(1)
1603 *   {
1604 * @endcode
1605 *
1606 * We allocate a vector to hold iterators to all active neighbors of
1607 * a cell. We reserve the maximal number of active neighbors in order to
1608 * avoid later reallocations. Note how this maximal number of active
1609 * neighbors is computed here.
1610 *
1611 * @code
1612 *   active_neighbors.reserve(GeometryInfo<dim>::faces_per_cell *
1613 *   GeometryInfo<dim>::max_children_per_face);
1614 *   }
1615 *  
1616 *  
1617 *   template <int dim>
1618 *   GradientEstimation::EstimateScratchData<dim>::EstimateScratchData(
1619 *   const EstimateScratchData &scratch_data)
1620 *   : fe_midpoint_value(scratch_data.fe_midpoint_value.get_fe(),
1621 *   scratch_data.fe_midpoint_value.get_quadrature(),
1622 *   update_values | update_quadrature_points)
1623 *   , solution(scratch_data.solution)
1624 *   , error_per_cell(scratch_data.error_per_cell)
1625 *   , cell_midpoint_value(1)
1626 *   , neighbor_midpoint_value(1)
1627 *   {}
1628 *  
1629 *  
1630 * @endcode
1631 *
1632 * Next comes the implementation of the <code>GradientEstimation</code>
1633 * class. The first function does not much except for delegating work to the
1634 * other function, but there is a bit of setup at the top.
1635 *
1636
1637 *
1638 * Before starting with the work, we check that the vector into
1639 * which the results are written has the right size, using the `Assert`
1640 * macro and the exception class we declared above. Programming
1641 * mistakes in which one forgets to size arguments correctly at the
1642 * calling site are quite common. Because the resulting damage from
1643 * not catching such errors is often subtle (e.g., corruption of
1644 * data somewhere in memory, or non-reproducible results), it is
1645 * well worth the effort to check for such things.
1646 *
1647 * @code
1648 *   template <int dim>
1649 *   void GradientEstimation::estimate(const DoFHandler<dim> &dof_handler,
1650 *   const Vector<double> &solution,
1651 *   Vector<float> &error_per_cell)
1652 *   {
1653 *   Assert(
1654 *   error_per_cell.size() == dof_handler.get_triangulation().n_active_cells(),
1655 *   ExcInvalidVectorLength(error_per_cell.size(),
1656 *   dof_handler.get_triangulation().n_active_cells()));
1657 *  
1658 *   WorkStream::run(dof_handler.begin_active(),
1659 *   dof_handler.end(),
1660 *   &GradientEstimation::template estimate_cell<dim>,
1661 *   std::function<void(const EstimateCopyData &)>(),
1662 *   EstimateScratchData<dim>(dof_handler.get_fe(),
1663 *   solution,
1664 *   error_per_cell),
1665 *   EstimateCopyData());
1666 *   }
1667 *  
1668 *  
1669 * @endcode
1670 *
1671 * Here comes the function that estimates the local error by computing the
1672 * finite difference approximation of the gradient. The function first
1673 * computes the list of active neighbors of the present cell and then
1674 * computes the quantities described in the introduction for each of
1675 * the neighbors. The reason for this order is that it is not a one-liner
1676 * to find a given neighbor with locally refined meshes. In principle, an
1677 * optimized implementation would find neighbors and the quantities
1678 * depending on them in one step, rather than first building a list of
1679 * neighbors and in a second step their contributions but we will gladly
1680 * leave this as an exercise. As discussed before, the worker function
1681 * passed to WorkStream::run works on "scratch" objects that keep all
1682 * temporary objects. This way, we do not need to create and initialize
1683 * objects that are expensive to initialize within the function that does
1684 * the work every time it is called for a given cell. Such an argument is
1685 * passed as the second argument. The third argument would be a "copy-data"
1686 * object (see @ref threads for more information) but we do not actually use
1687 * any of these here. Since WorkStream::run() insists on passing three
1688 * arguments, we declare this function with three arguments, but simply
1689 * ignore the last one.
1690 *
1691
1692 *
1693 * (This is unsatisfactory from an aesthetic perspective. It can be avoided
1694 * by using an anonymous (lambda) function. If you allow, let us here show
1695 * how. First, assume that we had declared this function to only take two
1696 * arguments by omitting the unused last one. Now, WorkStream::run still
1697 * wants to call this function with three arguments, so we need to find a
1698 * way to "forget" the third argument in the call. Simply passing
1699 * WorkStream::run the pointer to the function as we do above will not do
1700 * this -- the compiler will complain that a function declared to have two
1701 * arguments is called with three arguments. However, we can do this by
1702 * passing the following as the third argument to WorkStream::run():
1703 * <div class=CodeFragmentInTutorialComment>
1704 * @code
1705 * [](const typename DoFHandler<dim>::active_cell_iterator &cell,
1706 * EstimateScratchData<dim> & scratch_data,
1707 * EstimateCopyData &)
1708 * {
1709 * GradientEstimation::estimate_cell<dim>(cell, scratch_data);
1710 * }
1711 * @endcode
1712 * </div>
1713 * This is not much better than the solution implemented below: either the
1714 * routine itself must take three arguments or it must be wrapped by
1715 * something that takes three arguments. We don't use this since adding the
1716 * unused argument at the beginning is simpler.
1717 *
1718
1719 *
1720 * Now for the details:
1721 *
1722 * @code
1723 *   template <int dim>
1724 *   void GradientEstimation::estimate_cell(
1725 *   const typename DoFHandler<dim>::active_cell_iterator &cell,
1726 *   EstimateScratchData<dim> &scratch_data,
1727 *   const EstimateCopyData &)
1728 *   {
1729 * @endcode
1730 *
1731 * We need space for the tensor <code>Y</code>, which is the sum of
1732 * outer products of the y-vectors.
1733 *
1734 * @code
1735 *   Tensor<2, dim> Y;
1736 *  
1737 * @endcode
1738 *
1739 * First initialize the <code>FEValues</code> object, as well as the
1740 * <code>Y</code> tensor:
1741 *
1742 * @code
1743 *   scratch_data.fe_midpoint_value.reinit(cell);
1744 *  
1745 * @endcode
1746 *
1747 * Now, before we go on, we first compute a list of all active neighbors
1748 * of the present cell. We do so by first looping over all faces and see
1749 * whether the neighbor there is active, which would be the case if it
1750 * is on the same level as the present cell or one level coarser (note
1751 * that a neighbor can only be once coarser than the present cell, as
1752 * we only allow a maximal difference of one refinement over a face in
1753 * deal.II). Alternatively, the neighbor could be on the same level
1754 * and be further refined; then we have to find which of its children
1755 * are next to the present cell and select these (note that if a child
1756 * of a neighbor of an active cell that is next to this active cell,
1757 * needs necessarily be active itself, due to the one-refinement rule
1758 * cited above).
1759 *
1760
1761 *
1762 * Things are slightly different in one space dimension, as there the
1763 * one-refinement rule does not exist: neighboring active cells may
1764 * differ in as many refinement levels as they like. In this case, the
1765 * computation becomes a little more difficult, but we will explain
1766 * this below.
1767 *
1768
1769 *
1770 * Before starting the loop over all neighbors of the present cell, we
1771 * have to clear the array storing the iterators to the active
1772 * neighbors, of course.
1773 *
1774 * @code
1775 *   scratch_data.active_neighbors.clear();
1776 *   for (const auto face_n : cell->face_indices())
1777 *   if (!cell->at_boundary(face_n))
1778 *   {
1779 * @endcode
1780 *
1781 * First define an abbreviation for the iterator to the face and
1782 * the neighbor
1783 *
1784 * @code
1785 *   const auto face = cell->face(face_n);
1786 *   const auto neighbor = cell->neighbor(face_n);
1787 *  
1788 * @endcode
1789 *
1790 * Then check whether the neighbor is active. If it is, then it
1791 * is on the same level or one level coarser (if we are not in
1792 * 1d), and we are interested in it in any case.
1793 *
1794 * @code
1795 *   if (neighbor->is_active())
1796 *   scratch_data.active_neighbors.push_back(neighbor);
1797 *   else
1798 *   {
1799 * @endcode
1800 *
1801 * If the neighbor is not active, then check its children.
1802 *
1803 * @code
1804 *   if (dim == 1)
1805 *   {
1806 * @endcode
1807 *
1808 * To find the child of the neighbor which bounds to the
1809 * present cell, successively go to its right child if
1810 * we are left of the present cell (n==0), or go to the
1811 * left child if we are on the right (n==1), until we
1812 * find an active cell.
1813 *
1814 * @code
1815 *   auto neighbor_child = neighbor;
1816 *   while (neighbor_child->has_children())
1817 *   neighbor_child = neighbor_child->child(face_n == 0 ? 1 : 0);
1818 *  
1819 * @endcode
1820 *
1821 * As this used some non-trivial geometrical intuition,
1822 * we might want to check whether we did it right,
1823 * i.e., check whether the neighbor of the cell we found
1824 * is indeed the cell we are presently working
1825 * on. Checks like this are often useful and have
1826 * frequently uncovered errors both in algorithms like
1827 * the line above (where it is simple to involuntarily
1828 * exchange <code>n==1</code> for <code>n==0</code> or
1829 * the like) and in the library (the assumptions
1830 * underlying the algorithm above could either be wrong,
1831 * wrongly documented, or are violated due to an error
1832 * in the library). One could in principle remove such
1833 * checks after the program works for some time, but it
1834 * might be a good things to leave it in anyway to check
1835 * for changes in the library or in the algorithm above.
1836 *
1837
1838 *
1839 * Note that if this check fails, then this is certainly
1840 * an error that is irrecoverable and probably qualifies
1841 * as an internal error. We therefore use a predefined
1842 * exception class to throw here.
1843 *
1844 * @code
1845 *   Assert(neighbor_child->neighbor(face_n == 0 ? 1 : 0) == cell,
1846 *   ExcInternalError());
1847 *  
1848 * @endcode
1849 *
1850 * If the check succeeded, we push the active neighbor
1851 * we just found to the stack we keep:
1852 *
1853 * @code
1854 *   scratch_data.active_neighbors.push_back(neighbor_child);
1855 *   }
1856 *   else
1857 * @endcode
1858 *
1859 * If we are not in 1d, we collect all neighbor children
1860 * `behind' the subfaces of the current face and move on:
1861 *
1862 * @code
1863 *   for (unsigned int subface_n = 0; subface_n < face->n_children();
1864 *   ++subface_n)
1865 *   scratch_data.active_neighbors.push_back(
1866 *   cell->neighbor_child_on_subface(face_n, subface_n));
1867 *   }
1868 *   }
1869 *  
1870 * @endcode
1871 *
1872 * OK, now that we have all the neighbors, lets start the computation
1873 * on each of them. First we do some preliminaries: find out about the
1874 * center of the present cell and the solution at this point. The
1875 * latter is obtained as a vector of function values at the quadrature
1876 * points, of which there are only one, of course. Likewise, the
1877 * position of the center is the position of the first (and only)
1878 * quadrature point in real space.
1879 *
1880 * @code
1881 *   const Point<dim> this_center =
1882 *   scratch_data.fe_midpoint_value.quadrature_point(0);
1883 *  
1884 *   scratch_data.fe_midpoint_value.get_function_values(
1885 *   scratch_data.solution, scratch_data.cell_midpoint_value);
1886 *  
1887 * @endcode
1888 *
1889 * Now loop over all active neighbors and collect the data we
1890 * need.
1891 *
1892 * @code
1893 *   Tensor<1, dim> projected_gradient;
1894 *   for (const auto &neighbor : scratch_data.active_neighbors)
1895 *   {
1896 * @endcode
1897 *
1898 * Then get the center of the neighbor cell and the value of the
1899 * finite element function at that point. Note that for this
1900 * information we have to reinitialize the <code>FEValues</code>
1901 * object for the neighbor cell.
1902 *
1903 * @code
1904 *   scratch_data.fe_midpoint_value.reinit(neighbor);
1905 *   const Point<dim> neighbor_center =
1906 *   scratch_data.fe_midpoint_value.quadrature_point(0);
1907 *  
1908 *   scratch_data.fe_midpoint_value.get_function_values(
1909 *   scratch_data.solution, scratch_data.neighbor_midpoint_value);
1910 *  
1911 * @endcode
1912 *
1913 * Compute the vector <code>y</code> connecting the centers of the
1914 * two cells. Note that as opposed to the introduction, we denote
1915 * by <code>y</code> the normalized difference vector, as this is
1916 * the quantity used everywhere in the computations.
1917 *
1918 * @code
1919 *   Tensor<1, dim> y = neighbor_center - this_center;
1920 *   const double distance = y.norm();
1921 *   y /= distance;
1922 *  
1923 * @endcode
1924 *
1925 * Then add up the contribution of this cell to the Y matrix...
1926 *
1927 * @code
1928 *   for (unsigned int i = 0; i < dim; ++i)
1929 *   for (unsigned int j = 0; j < dim; ++j)
1930 *   Y[i][j] += y[i] * y[j];
1931 *  
1932 * @endcode
1933 *
1934 * ... and update the sum of difference quotients:
1935 *
1936 * @code
1937 *   projected_gradient += (scratch_data.neighbor_midpoint_value[0] -
1938 *   scratch_data.cell_midpoint_value[0]) /
1939 *   distance * y;
1940 *   }
1941 *  
1942 * @endcode
1943 *
1944 * If now, after collecting all the information from the neighbors, we
1945 * can determine an approximation of the gradient for the present
1946 * cell, then we need to have passed over vectors <code>y</code> which
1947 * span the whole space, otherwise we would not have all components of
1948 * the gradient. This is indicated by the invertibility of the matrix.
1949 *
1950
1951 *
1952 * If the matrix is not invertible, then the present
1953 * cell had an insufficient number of active neighbors. In contrast to
1954 * all previous cases (where we raised exceptions) this is, however,
1955 * not a programming error: it is a runtime error that can happen in
1956 * optimized mode even if it ran well in debug mode, so it is
1957 * reasonable to try to catch this error also in optimized mode. For
1958 * this case, there is the <code>AssertThrow</code> macro: it checks
1959 * the condition like the <code>Assert</code> macro, but not only in
1960 * debug mode; it then outputs an error message, but instead of
1961 * aborting the program as in the case of the <code>Assert</code>
1962 * macro, the exception is thrown using the <code>throw</code> command
1963 * of C++. This way, one has the possibility to catch this error and
1964 * take reasonable counter actions. One such measure would be to
1965 * refine the grid globally, as the case of insufficient directions
1966 * can not occur if every cell of the initial grid has been refined at
1967 * least once.
1968 *
1969 * @code
1970 *   AssertThrow(determinant(Y) != 0, ExcInsufficientDirections());
1971 *  
1972 * @endcode
1973 *
1974 * If, on the other hand, the matrix is invertible, then invert it,
1975 * multiply the other quantity with it, and compute the estimated error
1976 * using this quantity and the correct powers of the mesh width:
1977 *
1978 * @code
1979 *   const Tensor<2, dim> Y_inverse = invert(Y);
1980 *  
1981 *   const Tensor<1, dim> gradient = Y_inverse * projected_gradient;
1982 *  
1983 * @endcode
1984 *
1985 * The last part of this function is the one where we write into
1986 * the element of the output vector what we have just
1987 * computed. The address of this vector has been stored in the
1988 * scratch data object, and all we have to do is know how to get
1989 * at the correct element inside this vector -- but we can ask the
1990 * cell we're on the how-manyth active cell it is for this:
1991 *
1992 * @code
1993 *   scratch_data.error_per_cell(cell->active_cell_index()) =
1994 *   (std::pow(cell->diameter(), 1 + 1.0 * dim / 2) * gradient.norm());
1995 *   }
1996 *   } // namespace Step9
1997 *  
1998 *  
1999 * @endcode
2000 *
2001 *
2002 * <a name="step_9-Mainfunction"></a>
2003 * <h3>Main function</h3>
2004 *
2005
2006 *
2007 * The <code>main</code> function is similar to the previous examples. The
2008 * primary difference is that we use MultithreadInfo to set the maximum
2009 * number of threads (see the documentation topic @ref threads
2010 * "Parallel computing with multiple processors accessing shared memory"
2011 * for more information). The number of threads used is the minimum of the
2012 * environment variable DEAL_II_NUM_THREADS and the parameter of
2013 * <code>set_thread_limit</code>. If no value is given to
2014 * <code>set_thread_limit</code>, the default value from the Intel Threading
2015 * Building Blocks (TBB) library is used. If the call to
2016 * <code>set_thread_limit</code> is omitted, the number of threads will be
2017 * chosen by TBB independently of DEAL_II_NUM_THREADS.
2018 *
2019 * @code
2020 *   int main()
2021 *   {
2022 *   using namespace dealii;
2023 *   try
2024 *   {
2026 *  
2027 *   Step9::AdvectionProblem<2> advection_problem_2d;
2028 *   advection_problem_2d.run();
2029 *   }
2030 *   catch (std::exception &exc)
2031 *   {
2032 *   std::cerr << std::endl
2033 *   << std::endl
2034 *   << "----------------------------------------------------"
2035 *   << std::endl;
2036 *   std::cerr << "Exception on processing: " << std::endl
2037 *   << exc.what() << std::endl
2038 *   << "Aborting!" << std::endl
2039 *   << "----------------------------------------------------"
2040 *   << std::endl;
2041 *   return 1;
2042 *   }
2043 *   catch (...)
2044 *   {
2045 *   std::cerr << std::endl
2046 *   << std::endl
2047 *   << "----------------------------------------------------"
2048 *   << std::endl;
2049 *   std::cerr << "Unknown exception!" << std::endl
2050 *   << "Aborting!" << std::endl
2051 *   << "----------------------------------------------------"
2052 *   << std::endl;
2053 *   return 1;
2054 *   }
2055 *  
2056 *   return 0;
2057 *   }
2058 * @endcode
2059<a name="step_9-Results"></a><h1>Results</h1>
2060
2061
2062
2063The results of this program are not particularly spectacular. They
2064consist of the console output, some grid files, and the solution on
2065each of these grids. First for the console output:
2066@code
2067Cycle 0:
2068 Number of active cells: 64
2069 Number of degrees of freedom: 1681
2070 Iterations required for convergence: 338
2071 Norm of residual at convergence: 2.9177e-11
2072 Grid written to grid-0.vtu
2073 Solution written to solution-0.vtu
2074Cycle 1:
2075 Number of active cells: 121
2076 Number of degrees of freedom: 3436
2077 Iterations required for convergence: 448
2078 Norm of residual at convergence: 2.40051e-11
2079 Grid written to grid-1.vtu
2080 Solution written to solution-1.vtu
2081Cycle 2:
2082 Number of active cells: 238
2083 Number of degrees of freedom: 6487
2084 Iterations required for convergence: 535
2085 Norm of residual at convergence: 2.06897e-11
2086 Grid written to grid-2.vtu
2087 Solution written to solution-2.vtu
2088Cycle 3:
2089 Number of active cells: 481
2090 Number of degrees of freedom: 13510
2091 Iterations required for convergence: 669
2092 Norm of residual at convergence: 1.41502e-11
2093 Grid written to grid-3.vtu
2094 Solution written to solution-3.vtu
2095Cycle 4:
2096 Number of active cells: 958
2097 Number of degrees of freedom: 26137
2098 Iterations required for convergence: 1039
2099 Norm of residual at convergence: 1.27042e-11
2100 Grid written to grid-4.vtu
2101 Solution written to solution-4.vtu
2102Cycle 5:
2103 Number of active cells: 1906
2104 Number of degrees of freedom: 52832
2105 Iterations required for convergence: 1345
2106 Norm of residual at convergence: 1.00694e-11
2107 Grid written to grid-5.vtu
2108 Solution written to solution-5.vtu
2109Cycle 6:
2110 Number of active cells: 3829
2111 Number of degrees of freedom: 104339
2112 Iterations required for convergence: 1976
2113 Norm of residual at convergence: 7.43452e-12
2114 Grid written to grid-6.vtu
2115 Solution written to solution-6.vtu
2116Cycle 7:
2117 Number of active cells: 7414
2118 Number of degrees of freedom: 201946
2119 Iterations required for convergence: 2256
2120 Norm of residual at convergence: 7.88403e-12
2121 Grid written to grid-7.vtu
2122 Solution written to solution-7.vtu
2123Cycle 8:
2124 Number of active cells: 14413
2125 Number of degrees of freedom: 389558
2126 Iterations required for convergence: 2968
2127 Norm of residual at convergence: 6.72725e-12
2128 Grid written to grid-8.vtu
2129 Solution written to solution-8.vtu
2130Cycle 9:
2131 Number of active cells: 28141
2132 Number of degrees of freedom: 750187
2133 Iterations required for convergence: 3885
2134 Norm of residual at convergence: 5.86246e-12
2135 Grid written to grid-9.vtu
2136 Solution written to solution-9.vtu
2137@endcode
2138
2139Quite a number of cells are used on the finest level to resolve the features of
2140the solution. Here are the fourth and tenth grids:
2141<div class="twocolumn" style="width: 80%">
2142 <div>
2143 <img src="https://www.dealii.org/images/steps/developer/step-9-grid-3.png"
2144 alt="Fourth grid in the refinement cycle, showing some adaptivity to features."
2145 width="400" height="400">
2146 </div>
2147 <div>
2148 <img src="https://www.dealii.org/images/steps/developer/step-9-grid-9.png"
2149 alt="Tenth grid in the refinement cycle, showing that the waves are fully captured."
2150 width="400" height="400">
2151 </div>
2152</div>
2153and the fourth and tenth solutions:
2154<div class="twocolumn" style="width: 80%">
2155 <div>
2156 <img src="https://www.dealii.org/images/steps/developer/step-9-solution-3.png"
2157 alt="Fourth solution, showing that we resolve most features but some
2158 are sill unresolved and appear blury."
2159 width="400" height="400">
2160 </div>
2161 <div>
2162 <img src="https://www.dealii.org/images/steps/developer/step-9-solution-9.png"
2163 alt="Tenth solution, showing a fully resolved flow."
2164 width="400" height="400">
2165 </div>
2166</div>
2167and both the grid and solution zoomed in:
2168<div class="twocolumn" style="width: 80%">
2169 <div>
2170 <img src="https://www.dealii.org/images/steps/developer/step-9-solution-3-zoom.png"
2171 alt="Detail of the fourth solution, showing that we resolve most
2172 features but some are sill unresolved and appear blury. In particular,
2173 the larger cells need to be refined."
2174 width="400" height="400">
2175 </div>
2176 <div>
2177 <img src="https://www.dealii.org/images/steps/developer/step-9-solution-9-zoom.png"
2178 alt="Detail of the tenth solution, showing that we needed a lot more
2179 cells than were present in the fourth solution."
2180 width="400" height="400">
2181 </div>
2182</div>
2183
2184The solution is created by that part that is transported along the wiggly
2185advection field from the left and lower boundaries to the top right, and the
2186part that is created by the source in the lower left corner, and the results of
2187which are also transported along. The grid shown above is well-adapted to
2188resolve these features. The comparison between plots shows that, even though we
2189are using a high-order approximation, we still need adaptive mesh refinement to
2190fully resolve the wiggles.
2191 *
2192 *
2193<a name="step_9-PlainProg"></a>
2194<h1> The plain program</h1>
2195@include "step-9.cc"
2196*/
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, level_dof_access > > &cell)
static void set_thread_limit(const unsigned int max_threads=numbers::invalid_unsigned_int)
void initialize(const MatrixType &A, const AdditionalData &parameters=AdditionalData())
unsigned int n_active_cells() const
virtual void execute_coarsening_and_refinement() override
Definition tria.cc:3320
Point< 2 > second
Definition grid_out.cc:4624
Point< 2 > first
Definition grid_out.cc:4623
unsigned int level
Definition grid_out.cc:4626
__global__ void set(Number *val, const Number s, const size_type N)
#define Assert(cond, exc)
typename ActiveSelector::active_cell_iterator active_cell_iterator
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
Definition loop.h:564
@ update_values
Shape function values.
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
void refine_and_coarsen_fixed_number(Triangulation< dim, spacedim > &triangulation, const Vector< Number > &criteria, const double top_fraction_of_cells, const double bottom_fraction_of_cells, const unsigned int max_n_cells=std::numeric_limits< unsigned int >::max())
@ 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 > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
VectorType::value_type * end(VectorType &V)
T sum(const T &t, const MPI_Comm mpi_communicator)
void run(const std::vector< std::vector< Iterator > > &colored_iterators, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length=2 *MultithreadInfo::n_threads(), const unsigned int chunk_size=8)
bool check(const ConstraintKinds kind_in, const unsigned int dim)
void copy(const T *begin, const T *end, U *dest)
int(& functions)(const void *v1, const void *v2)
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
Definition loop.h:70
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
DEAL_II_HOST constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)