521 *
prm.declare_entry(
"Epsilon",
524 *
"Diffusion parameter");
526 *
prm.declare_entry(
"Fe degree",
529 *
"Finite Element degree");
530 *
prm.declare_entry(
"Smoother type",
533 *
"Select smoother: SOR|Jacobi|block SOR|block Jacobi");
534 *
prm.declare_entry(
"Smoothing steps",
537 *
"Number of smoothing steps");
542 *
"Select DoF renumbering: none|downstream|upstream|random");
543 *
prm.declare_entry(
"With streamline diffusion",
546 *
"Enable streamline diffusion stabilization: true|false");
547 *
prm.declare_entry(
"Output",
550 *
"Generate graphical output: true|false");
557 *
false, ExcMessage(
"Please pass a .prm file as the first argument!"));
562 *
epsilon = prm.get_double(
"Epsilon");
563 *
fe_degree = prm.get_integer(
"Fe degree");
564 *
smoother_type = prm.get(
"Smoother type");
567 *
const std::string
renumbering = prm.get(
"DoF renumbering");
569 *
dof_renumbering = DoFRenumberingStrategy::none;
571 *
dof_renumbering = DoFRenumberingStrategy::downstream;
575 *
dof_renumbering = DoFRenumberingStrategy::random;
578 *
ExcMessage(
"The <DoF renumbering> parameter has "
579 *
"an invalid value."));
582 *
output = prm.get_bool(
"Output");
589 * <a name=
"step_63-Cellpermutations"></a>
618 *
std::vector<unsigned int>
621 *
const unsigned int level)
623 *
std::vector<typename DoFHandler<dim>::level_cell_iterator>
ordered_cells;
625 *
for (
const auto &cell : dof_handler.cell_iterators_on_level(
level))
645 *
std::vector<unsigned int>
649 *
std::vector<typename DoFHandler<dim>::active_cell_iterator>
ordered_cells;
650 *
ordered_cells.reserve(dof_handler.get_triangulation().n_active_cells());
651 *
for (
const auto &cell : dof_handler.active_cell_iterators())
660 *
ordered_indices.reserve(dof_handler.get_triangulation().n_active_cells());
682 *
std::vector<unsigned int>
684 *
const unsigned int level)
688 *
for (
const auto &cell : dof_handler.cell_iterators_on_level(
level))
702 *
std::vector<unsigned int>
706 *
ordered_cells.reserve(dof_handler.get_triangulation().n_active_cells());
707 *
for (
const auto &cell : dof_handler.active_cell_iterators())
722 * <a name=
"step_63-Righthandsideandboundaryvalues"></a>
729 *
href=
"https://global.oup.com/academic/product/finite-elements-and-fast-iterative-solvers-9780199678808">
749 *
virtual void value_list(
const std::vector<
Point<dim>> &points,
750 *
std::vector<double> &values,
751 *
const unsigned int component = 0)
const override;
758 *
const unsigned int component)
const
760 *
Assert(component == 0, ExcIndexRange(component, 0, 1));
770 *
std::vector<double> &values,
771 *
const unsigned int component)
const
775 *
for (
unsigned int i = 0; i < points.size(); ++i)
792 *
const unsigned int component = 0)
const override;
794 *
virtual void value_list(
const std::vector<
Point<dim>> &points,
795 *
std::vector<double> &
values,
796 *
const unsigned int component = 0)
const override;
803 *
const unsigned int component)
const
805 *
Assert(component == 0, ExcIndexRange(component, 0, 1));
813 *
if (std::fabs(p[0] - 1) < 1e-8 ||
814 *
(std::fabs(p[1] + 1) < 1e-8 && p[0] >= 0.5))
828 *
std::vector<double> &values,
829 *
const unsigned int component)
const
833 *
for (
unsigned int i = 0; i < points.size(); ++i)
842 * <a name=
"step_63-Streamlinediffusionimplementation"></a>
850 *
href=
"https://link.springer.com/chapter/10.1007/978-3-540-34288-5_27">
On
862 *
const double coth =
865 *
return hk / (2.0 * dir.norm() *
pk) * (coth - 1.0 /
Peclet);
872 * <a name=
"step_63-codeAdvectionProblemcodeclass"></a>
900 *
template <
class IteratorType>
903 *
CopyData ©_data);
909 *
void refine_grid();
937 *
std::unique_ptr<MGSmoother<Vector<double>>>
mg_smoother;
957 *
, fe(settings.fe_degree)
958 *
, mapping(settings.fe_degree)
959 *
, settings(settings)
972 * <a name=
"step_63-codeAdvectionProblemsetup_systemcode"></a>
993 *
dof_handler.distribute_dofs(fe);
995 *
solution.reinit(dof_handler.n_dofs());
998 *
constraints.
clear();
1005 *
constraints.close();
1013 *
sparsity_pattern.copy_from(
dsp);
1014 *
system_matrix.reinit(sparsity_pattern);
1016 *
dof_handler.distribute_mg_dofs();
1031 *
if (settings.smoother_type ==
"SOR" || settings.smoother_type ==
"Jacobi")
1033 *
if (settings.dof_renumbering ==
1034 *
Settings::DoFRenumberingStrategy::downstream ||
1035 *
settings.dof_renumbering ==
1036 *
Settings::DoFRenumberingStrategy::upstream)
1039 *
(settings.dof_renumbering ==
1040 *
Settings::DoFRenumberingStrategy::upstream ?
1051 *
else if (settings.dof_renumbering ==
1052 *
Settings::DoFRenumberingStrategy::random)
1069 *
mg_constrained_dofs.
clear();
1070 *
mg_constrained_dofs.initialize(dof_handler);
1072 *
mg_constrained_dofs.make_zero_boundary_constraints(dof_handler, {0, 1});
1087 *
dof_handler.n_dofs(
level));
1094 *
dof_handler.n_dofs(
level));
1096 *
mg_constrained_dofs,
1111 * <a name=
"step_63-codeAdvectionProblemassemble_cellcode"></a>
1129 *
ScratchData<dim> &scratch_data,
1130 *
CopyData ©_data)
1132 *
copy_data.level = cell->level();
1134 *
const unsigned int dofs_per_cell =
1135 *
scratch_data.fe_values.get_fe().n_dofs_per_cell();
1136 *
copy_data.dofs_per_cell = dofs_per_cell;
1137 *
copy_data.cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
1139 *
const unsigned int n_q_points =
1140 *
scratch_data.fe_values.get_quadrature().size();
1142 *
if (cell->is_level_cell() ==
false)
1143 *
copy_data.cell_rhs.reinit(dofs_per_cell);
1145 *
copy_data.local_dof_indices.resize(dofs_per_cell);
1146 *
cell->get_active_or_mg_dof_indices(copy_data.local_dof_indices);
1148 *
scratch_data.fe_values.reinit(cell);
1151 *
std::vector<double>
rhs_values(n_q_points);
1153 *
right_hand_side.value_list(scratch_data.fe_values.get_quadrature_points(),
1165 *
const double delta = (settings.with_streamline_diffusion ?
1169 *
settings.fe_degree) :
1173 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1175 *
for (
unsigned int j = 0;
j < dofs_per_cell; ++
j)
1183 *
copy_data.cell_matrix(i,
j) +=
1184 *
(settings.epsilon *
1185 *
scratch_data.fe_values.shape_grad(i,
q_point) *
1186 *
scratch_data.fe_values.shape_grad(
j,
q_point) *
1187 *
scratch_data.fe_values.JxW(
q_point)) +
1188 *
(scratch_data.fe_values.shape_value(i,
q_point) *
1190 *
scratch_data.fe_values.shape_grad(
j,
q_point)) *
1191 *
scratch_data.fe_values.JxW(
q_point))
1199 *
scratch_data.fe_values.shape_grad(
j,
q_point)) *
1201 *
scratch_data.fe_values.shape_grad(i,
q_point)) *
1202 *
scratch_data.fe_values.JxW(
q_point) -
1203 *
delta * settings.epsilon *
1206 *
scratch_data.fe_values.shape_grad(i,
q_point)) *
1207 *
scratch_data.fe_values.JxW(
q_point);
1209 *
if (cell->is_level_cell() ==
false)
1217 *
copy_data.cell_rhs(i) +=
1218 *
scratch_data.fe_values.shape_value(i,
q_point) *
1226 *
scratch_data.fe_values.shape_grad(i,
q_point) *
1227 *
scratch_data.fe_values.JxW(
q_point);
1236 * <a name=
"step_63-codeAdvectionProblemassemble_system_and_multigridcode"></a>
1252 *
[&](
const decltype(dof_handler.begin_active()) &cell,
1254 *
CopyData ©_data) {
1258 *
const auto copier_active = [&](
const CopyData ©_data) {
1259 *
constraints.distribute_local_to_global(copy_data.cell_matrix,
1260 *
copy_data.cell_rhs,
1261 *
copy_data.local_dof_indices,
1268 *
dof_handler.end(),
1288 *
dof_handler.locally_owned_mg_dofs(
level),
1292 *
mg_constrained_dofs.get_refinement_edge_indices(
level))
1295 *
mg_constrained_dofs.get_boundary_indices(
level))
1301 *
[&](
const decltype(dof_handler.begin_mg()) &cell,
1303 *
CopyData ©_data) {
1307 *
const auto copier_mg = [&](
const CopyData ©_data) {
1309 *
copy_data.cell_matrix,
1310 *
copy_data.local_dof_indices,
1316 * `interface_in` dof
pair.
Note:
For `interface_in`,
we load
1322 *
we must store both `interface_in`
and `interface_out`
1326 *
for (
unsigned int i = 0; i < copy_data.dofs_per_cell; ++i)
1327 *
for (
unsigned int j = 0;
j < copy_data.dofs_per_cell; ++
j)
1328 *
if (mg_constrained_dofs.is_interface_matrix_entry(
1330 *
copy_data.local_dof_indices[i],
1331 *
copy_data.local_dof_indices[
j]))
1334 *
copy_data.local_dof_indices[i],
1335 *
copy_data.local_dof_indices[
j],
1336 *
copy_data.cell_matrix(i,
j));
1338 *
copy_data.local_dof_indices[i],
1339 *
copy_data.local_dof_indices[
j],
1340 *
copy_data.cell_matrix(
j, i));
1345 *
dof_handler.end_mg(),
1357 * <a name=
"step_63-codeAdvectionProblemsetup_smoothercode"></a>
1365 * relaxation parameter.
1387 * carries new information to all its DoFs.
1391 * Finally, as mentioned above, the point smoothers only operate on DoFs, and
1392 * the block smoothers on cells, so only the block smoothers need to be given
1393 * information regarding cell orderings. DoF ordering for point smoothers has
1394 * already been taken care of in `setup_system()`.
1400 * template <int dim>
1401 * void AdvectionProblem<dim>::setup_smoother()
1403 * if (settings.smoother_type == "SOR")
1405 * using Smoother = PreconditionSOR<SparseMatrix<double>>;
1408 * std::make_unique<MGSmootherPrecondition<SparseMatrix<double>,
1410 * Vector<double>>>();
1411 * smoother->initialize(mg_matrices,
1412 * Smoother::AdditionalData(fe.degree == 1 ? 1.0 :
1414 * smoother->set_steps(settings.smoothing_steps);
1415 * mg_smoother = std::move(smoother);
1417 * else if (settings.smoother_type == "Jacobi")
1419 * using Smoother = PreconditionJacobi<SparseMatrix<double>>;
1421 * std::make_unique<MGSmootherPrecondition<SparseMatrix<double>,
1423 * Vector<double>>>();
1424 * smoother->initialize(mg_matrices,
1425 * Smoother::AdditionalData(fe.degree == 1 ? 0.6667 :
1427 * smoother->set_steps(settings.smoothing_steps);
1428 * mg_smoother = std::move(smoother);
1430 * else if (settings.smoother_type == "block SOR" ||
1431 * settings.smoother_type == "block Jacobi")
1433 * smoother_data.resize(0, triangulation.n_levels() - 1);
1435 * for (unsigned int level = 0; level < triangulation.n_levels(); ++level)
1437 * DoFTools::make_cell_patches(smoother_data[level].block_list,
1441 * smoother_data[level].relaxation =
1442 * (settings.smoother_type == "block SOR" ? 1.0 : 0.25);
1443 * smoother_data[level].inversion = PreconditionBlockBase<double>::svd;
1445 * std::vector<unsigned int> ordered_indices;
1446 * switch (settings.dof_renumbering)
1448 * case Settings::DoFRenumberingStrategy::downstream:
1450 * create_downstream_cell_ordering(dof_handler,
1451 * advection_direction,
1455 * case Settings::DoFRenumberingStrategy::upstream:
1457 * create_downstream_cell_ordering(dof_handler,
1458 * -1.0 * advection_direction,
1462 * case Settings::DoFRenumberingStrategy::random:
1464 * create_random_cell_ordering(dof_handler, level);
1467 * case Settings::DoFRenumberingStrategy::none:
1471 * AssertThrow(false, ExcNotImplemented());
1475 * smoother_data[level].order =
1476 * std::vector<std::vector<unsigned int>>(1, ordered_indices);
1479 * if (settings.smoother_type == "block SOR")
1481 * auto smoother = std::make_unique<MGSmootherPrecondition<
1482 * SparseMatrix<double>,
1483 * RelaxationBlockSOR<SparseMatrix<double>, double, Vector<double>>,
1484 * Vector<double>>>();
1485 * smoother->initialize(mg_matrices, smoother_data);
1486 * smoother->set_steps(settings.smoothing_steps);
1487 * mg_smoother = std::move(smoother);
1489 * else if (settings.smoother_type == "block Jacobi")
1491 * auto smoother = std::make_unique<
1492 * MGSmootherPrecondition<SparseMatrix<double>,
1493 * RelaxationBlockJacobi<SparseMatrix<double>,
1496 * Vector<double>>>();
1497 * smoother->initialize(mg_matrices, smoother_data);
1498 * smoother->set_steps(settings.smoothing_steps);
1499 * mg_smoother = std::move(smoother);
1503 * AssertThrow(false, ExcNotImplemented());
1510 * <a name="step_63-codeAdvectionProblemsolvecode"></a>
1511 * <h4><code>AdvectionProblem::solve()</code></h4>
1515 * Before we can solve the system, we must first set up the multigrid
1516 * preconditioner. This requires the setup of the transfer between levels,
1517 * the coarse matrix solver, and the smoother. This setup follows almost
1518 * identically to @ref step_16 "step-16", the main difference being the various smoothers
1519 * defined above and the fact that we need different interface edge matrices
1520 * for in and out since our problem is non-symmetric. (In reality, for this
1521 * tutorial these interface matrices are empty since we are only using global
1522 * refinement, and thus have no refinement edges. However, we have still
1523 * included both here since if one made the simple switch to an adaptively
1524 * refined method, the program would still run correctly.)
1528 * The last thing to note is that since our problem is non-symmetric, we must
1529 * use an appropriate Krylov subspace method. We choose here to
1530 * use GMRES since it offers the guarantee of residual reduction in each
1531 * iteration. The major disadvantage of GMRES is that, for each iteration,
1532 * the number of stored temporary vectors increases by one, and one also needs
1533 * to compute a scalar product with all previously stored vectors. This is
1534 * rather expensive. This requirement is relaxed by using the restarted GMRES
1535 * method which puts a cap on the number of vectors we are required to store
1536 * at any one time (here we restart after 50 temporary vectors, or 48
1537 * iterations). This then has the disadvantage that we lose information we
1538 * have gathered throughout the iteration and therefore we could see slower
1539 * convergence. As a consequence, where to restart is a question of balancing
1540 * memory consumption, CPU effort, and convergence speed.
1541 * However, the goal of this tutorial is to have very low
1542 * iteration counts by using a powerful GMG preconditioner, so we have picked
1543 * the restart length such that all of the results shown below converge prior
1544 * to restart happening, and thus we have a standard GMRES method. If the user
1545 * is interested, another suitable method offered in deal.II would be
1552 * template <int dim>
1553 * void AdvectionProblem<dim>::solve()
1555 * const unsigned int max_iters = 200;
1556 * const double solve_tolerance = 1e-8 * system_rhs.l2_norm();
1557 * SolverControl solver_control(max_iters, solve_tolerance, true, true);
1558 * solver_control.enable_history_data();
1560 * using Transfer = MGTransferPrebuilt<Vector<double>>;
1561 * Transfer mg_transfer(mg_constrained_dofs);
1562 * mg_transfer.build(dof_handler);
1564 * FullMatrix<double> coarse_matrix;
1565 * coarse_matrix.copy_from(mg_matrices[0]);
1566 * MGCoarseGridHouseholder<double, Vector<double>> coarse_grid_solver;
1567 * coarse_grid_solver.initialize(coarse_matrix);
1571 * mg_matrix.initialize(mg_matrices);
1572 * mg_interface_matrix_in.initialize(mg_interface_in);
1573 * mg_interface_matrix_out.initialize(mg_interface_out);
1575 * Multigrid<Vector<double>> mg(
1576 * mg_matrix, coarse_grid_solver, mg_transfer, *mg_smoother, *mg_smoother);
1577 * mg.set_edge_matrices(mg_interface_matrix_out, mg_interface_matrix_in);
1579 * PreconditionMG<dim, Vector<double>, Transfer> preconditioner(dof_handler,
1583 * std::cout << " Solving with GMRES to tol " << solve_tolerance << "..."
1585 * SolverGMRES<Vector<double>> solver(
1586 * solver_control, SolverGMRES<Vector<double>>::AdditionalData(50, true));
1590 * solver.solve(system_matrix, solution, system_rhs, preconditioner);
1593 * std::cout << " converged in " << solver_control.last_step()
1595 * << " in " << time.last_wall_time() << " seconds " << std::endl;
1597 * constraints.distribute(solution);
1599 * mg_smoother.release();
1606 * <a name="step_63-codeAdvectionProblemoutput_resultscode"></a>
1607 * <h4><code>AdvectionProblem::output_results()</code></h4>
1611 * The final function of interest generates graphical output.
1612 * Here we output the solution and cell ordering in a .vtu format.
1616 * At the top of the function, we generate an index for each cell to
1617 * visualize the ordering used by the smoothers. Note that we do
1618 * this only for the active cells instead of the levels, where the
1619 * smoothers are actually used. For the point smoothers we renumber
1620 * DoFs instead of cells, so this is only an approximation of what
1621 * happens in reality. Finally, the random ordering is not the
1622 * random ordering we actually use (see `create_smoother()` for that).
1626 * The (integer) ordering of cells is then copied into a (floating
1627 * point) vector for graphical output.
1630 * template <int dim>
1631 * void AdvectionProblem<dim>::output_results(const unsigned int cycle) const
1633 * const unsigned int n_active_cells = triangulation.n_active_cells();
1634 * Vector<double> cell_indices(n_active_cells);
1636 * std::vector<unsigned int> ordered_indices;
1637 * switch (settings.dof_renumbering)
1639 * case Settings::DoFRenumberingStrategy::downstream:
1641 * create_downstream_cell_ordering(dof_handler, advection_direction);
1644 * case Settings::DoFRenumberingStrategy::upstream:
1646 * create_downstream_cell_ordering(dof_handler,
1647 * -1.0 * advection_direction);
1650 * case Settings::DoFRenumberingStrategy::random:
1651 * ordered_indices = create_random_cell_ordering(dof_handler);
1654 * case Settings::DoFRenumberingStrategy::none:
1655 * ordered_indices.resize(n_active_cells);
1656 * for (unsigned int i = 0; i < n_active_cells; ++i)
1657 * ordered_indices[i] = i;
1661 * AssertThrow(false, ExcNotImplemented());
1665 * for (unsigned int i = 0; i < n_active_cells; ++i)
1666 * cell_indices(ordered_indices[i]) = static_cast<double>(i);
1671 * The remainder of the function is then straightforward, given
1672 * previous tutorial programs:
1675 * DataOut<dim> data_out;
1676 * data_out.attach_dof_handler(dof_handler);
1677 * data_out.add_data_vector(solution, "solution");
1678 * data_out.add_data_vector(cell_indices, "cell_index");
1679 * data_out.build_patches();
1681 * const std::string filename =
1682 * "solution-" + Utilities::int_to_string(cycle) + ".vtu";
1683 * std::ofstream output(filename);
1684 * data_out.write_vtu(output);
1691 * <a name="step_63-codeAdvectionProblemruncode"></a>
1692 * <h4><code>AdvectionProblem::run()</code></h4>
1696 * As in most tutorials, this function creates/refines the mesh and calls
1697 * the various functions defined above to set up, assemble, solve, and output
1702 * In cycle zero, we generate the mesh for the on the square
1703 * <code>[-1,1]^dim</code> with a hole of radius 3/10 units centered
1704 * at the origin. For objects with `manifold_id` equal to one
1705 * (namely, the faces adjacent to the hole), we assign a spherical
1712 * template <int dim>
1713 * void AdvectionProblem<dim>::run()
1715 * for (unsigned int cycle = 0; cycle < (settings.fe_degree == 1 ? 7u : 5u);
1718 * std::cout << " Cycle " << cycle << ':
' << std::endl;
1722 * GridGenerator::hyper_cube_with_cylindrical_hole(triangulation,
1726 * const SphericalManifold<dim> manifold_description(Point<dim>(0, 0));
1727 * triangulation.set_manifold(1, manifold_description);
1730 * triangulation.refine_global();
1734 * std::cout << " Number of active cells: "
1735 * << triangulation.n_active_cells() << " ("
1736 * << triangulation.n_levels() << " levels)" << std::endl;
1737 * std::cout << " Number of degrees of freedom: "
1738 * << dof_handler.n_dofs() << std::endl;
1740 * assemble_system_and_multigrid();
1744 * if (settings.output)
1745 * output_results(cycle);
1747 * std::cout << std::endl;
1750 * } // namespace Step63
1756 * <a name="step_63-Thecodemaincodefunction"></a>
1757 * <h3>The <code>main</code> function</h3>
1761 * Finally, the main function is like most tutorials. The only
1762 * interesting bit is that we require the user to pass a `.prm` file
1763 * as a sole command line argument. If no parameter file is given, the
1764 * program will output the contents of a sample parameter file with
1765 * all default values to the screen that the user can then copy and
1766 * paste into their own `.prm` file.
1772 * int main(int argc, char *argv[])
1776 * Step63::Settings settings;
1777 * settings.get_parameters((argc > 1) ? (argv[1]) : "");
1779 * Step63::AdvectionProblem<2> advection_problem_2d(settings);
1780 * advection_problem_2d.run();
1782 * catch (std::exception &exc)
1784 * std::cerr << std::endl
1786 * << "----------------------------------------------------"
1788 * std::cerr << "Exception on processing: " << std::endl
1789 * << exc.what() << std::endl
1790 * << "Aborting!" << std::endl
1791 * << "----------------------------------------------------"
1797 * std::cerr << std::endl
1799 * << "----------------------------------------------------"
1801 * std::cerr << "Unknown exception!" << std::endl
1802 * << "Aborting!" << std::endl
1803 * << "----------------------------------------------------"
1811<a name="step_63-Results"></a><h1>Results</h1>
1814<a name="step_63-GMRESIterationNumbers"></a><h3> GMRES Iteration Numbers </h3>
1817The major advantage for GMG is that it is an @f$\mathcal{O}(n)@f$ method,
1818that is, the complexity of the problem increases linearly with the
1819problem size. To show then that the linear solver presented in this
1820tutorial is in fact @f$\mathcal{O}(n)@f$, all one needs to do is show that
1821the iteration counts for the GMRES solve stay roughly constant as we
1824Each of the following tables gives the GMRES iteration counts to reduce the
1825initial residual by a factor of @f$10^8@f$. We selected a sufficient number of smoothing steps
1826(based on the method) to get iteration numbers independent of mesh size. As
1827can be seen from the tables below, the method is indeed @f$\mathcal{O}(n)@f$.
1829<a name="step_63-DoFCellRenumbering"></a><h4> DoF/Cell Renumbering </h4>
1832The point-wise smoothers ("Jacobi" and "SOR") get applied in the order the
1833DoFs are numbered on each level. We can influence this using the
1834DoFRenumbering namespace. The block smoothers are applied based on the
1835ordering we set in `setup_smoother()`. We can visualize this numbering. The
1836following pictures show the cell numbering of the active cells in downstream,
1837random, and upstream numbering (left to right):
1839<img src="https://www.dealii.org/images/steps/developer/step-63-cell-order.png" alt="">
1841Let us start with the additive smoothers. The following table shows
1842the number of iterations necessary to obtain convergence from GMRES:
1844<table align="center" class="doxtable">
1848 <th colspan="1">@f$Q_1@f$</th>
1849 <th colspan="7">Smoother (smoothing steps)</th>
1855 <th colspan="3">Jacobi (6)</th>
1857 <th colspan="3">Block Jacobi (3)</th>
1863 <th colspan="3">Renumbering Strategy</th>
1865 <th colspan="3">Renumbering Strategy</th>
1965We see that renumbering the
1966DoFs/cells has no effect on convergence speed. This is because these
1967smoothers compute operations on each DoF (point-smoother) or cell
1968(block-smoother) independently and add up the results. Since we can
1969define these smoothers as an application of a sum of matrices, and
1970matrix addition is commutative, the order at which we sum the
1971different components will not affect the end result.
1973On the other hand, the situation is different for multiplicative smoothers:
1975<table align="center" class="doxtable">
1979 <th colspan="1">@f$Q_1@f$</th>
1980 <th colspan="7">Smoother (smoothing steps)</th>
1986 <th colspan="3">SOR (3)</th>
1988 <th colspan="3">Block SOR (1)</th>
1994 <th colspan="3">Renumbering Strategy</th>
1996 <th colspan="3">Renumbering Strategy</th>
2096Here, we can speed up
2097convergence by renumbering the DoFs/cells in the advection direction,
2098and similarly, we can slow down convergence if we do the renumbering
2099in the opposite direction. This is because advection-dominated
2100problems have a directional flow of information (in the advection
2101direction) which, given the right renumbering of DoFs/cells,
2102multiplicative methods are able to capture.
2104This feature of multiplicative methods is, however, dependent on the
2105value of @f$\varepsilon@f$. As we increase @f$\varepsilon@f$ and the problem
2106becomes more diffusion-dominated, we have a more uniform propagation
2107of information over the mesh and there is a diminished advantage for
2108renumbering in the advection direction. On the opposite end, in the
2109extreme case of @f$\varepsilon=0@f$ (advection-only), we have a 1st-order
2110PDE and multiplicative methods with the right renumbering become
2111effective solvers: A correct downstream numbering may lead to methods
2112that require only a single iteration because information can be
2113propagated from the inflow boundary downstream, with no information
2114transport in the opposite direction. (Note, however, that in the case
2115of @f$\varepsilon=0@f$, special care must be taken for the boundary
2116conditions in this case).
2119<a name="step_63-Pointvsblocksmoothers"></a><h4> %Point vs. block smoothers </h4>
2122We will limit the results to runs using the downstream
2123renumbering. Here is a cross comparison of all four smoothers for both
2124@f$Q_1@f$ and @f$Q_3@f$ elements:
2126<table align="center" class="doxtable">
2130 <th colspan="1">@f$Q_1@f$</th>
2131 <th colspan="4">Smoother (smoothing steps)</th>
2133 <th colspan="1">@f$Q_3@f$</th>
2134 <th colspan="4">Smoother (smoothing steps)</th>
2137 <th colspan="1">Cells</th>
2139 <th colspan="1">DoFs</th>
2140 <th colspan="1">Jacobi (6)</th>
2141 <th colspan="1">Block Jacobi (3)</th>
2142 <th colspan="1">SOR (3)</th>
2143 <th colspan="1">Block SOR (1)</th>
2145 <th colspan="1">DoFs</th>
2146 <th colspan="1">Jacobi (6)</th>
2147 <th colspan="1">Block Jacobi (3)</th>
2148 <th colspan="1">SOR (3)</th>
2149 <th colspan="1">Block SOR (1)</th>
2248We see that for @f$Q_1@f$, both multiplicative smoothers require a smaller
2249combination of smoothing steps and iteration counts than either
2250additive smoother. However, when we increase the degree to a @f$Q_3@f$
2251element, there is a clear advantage for the block smoothers in terms
2252of the number of smoothing steps and iterations required to
2253solve. Specifically, the block SOR smoother gives constant iteration
2254counts over the degree, and the block Jacobi smoother only sees about
2255a 38% increase in iterations compared to 75% and 183% for Jacobi and
2258<a name="step_63-Cost"></a><h3> Cost </h3>
2261Iteration counts do not tell the full story in the optimality of a one
2262smoother over another. Obviously we must examine the cost of an
2263iteration. Block smoothers here are at a disadvantage as they are
2264having to construct and invert a cell matrix for each cell. Here is a
2265comparison of solve times for a @f$Q_3@f$ element with 74,496 DoFs:
2267<table align="center" class="doxtable">
2269 <th colspan="1">@f$Q_3@f$</th>
2270 <th colspan="4">Smoother (smoothing steps)</th>
2273 <th colspan="1">DoFs</th>
2274 <th colspan="1">Jacobi (6)</th>
2275 <th colspan="1">Block Jacobi (3)</th>
2276 <th colspan="1">SOR (3)</th>
2277 <th colspan="1">Block SOR (1)</th>
2288The smoother that requires the most iterations (Jacobi) actually takes
2289the shortest time (roughly 2/3 the time of the next fastest
2290method). This is because all that is required to apply a Jacobi
2291smoothing step is multiplication by a diagonal matrix which is very
2292cheap. On the other hand, while SOR requires over 3x more iterations
2293(each with 3x more smoothing steps) than block SOR, the times are
2294roughly equivalent, implying that a smoothing step of block SOR is
2295roughly 9x slower than a smoothing step of SOR. Lastly, block Jacobi
2296is almost 6x more expensive than block SOR, which intuitively makes
2297sense from the fact that 1 step of each method has the same cost
2298(inverting the cell matrices and either adding or multiply them
2299together), and block Jacobi has 3 times the number of smoothing steps per
2300iteration with 2 times the iterations.
2303<a name="step_63-Additionalpoints"></a><h3> Additional points </h3>
2306There are a few more important points to mention:
2309<li> For a mesh distributed in parallel, multiplicative methods cannot
2310be executed over the entire domain. This is because they operate one
2311cell at a time, and downstream cells can only be handled once upstream
2312cells have already been done. This is fine on a single processor: The
2313processor just goes through the list of cells one after the
2314other. However, in parallel, it would imply that some processors are
2315idle because upstream processors have not finished doing the work on
2316cells upstream from the ones owned by the current processor. Once the
2317upstream processors are done, the downstream ones can start, but by
2318that time the upstream processors have no work left. In other words,
2319most of the time during these smoother steps, most processors are in
2320fact idle. This is not how one obtains good parallel scalability!
2322One can use a hybrid method where
2323a multiplicative smoother is applied on each subdomain, but as you
2324increase the number of subdomains, the method approaches the behavior
2325of an additive method. This is a major disadvantage to these methods.
2328<li> Current research into block smoothers suggest that soon we will be
2329able to compute the inverse of the cell matrices much cheaper than
2330what is currently being done inside deal.II. This research is based on
2331the fast diagonalization method (dating back to the 1960s) and has
2332been used in the spectral community for around 20 years (see, e.g., <a
2333href="https://doi.org/10.1007/s10915-004-4787-3"> Hybrid
2334Multigrid/Schwarz Algorithms for the Spectral Element Method by Lottes
2335and Fischer</a>). There are currently efforts to generalize these
2336methods to DG and make them more robust. Also, it seems that one
2337should be able to take advantage of matrix-free implementations and
2338the fact that, in the interior of the domain, cell matrices tend to
2339look very similar, allowing fewer matrix inverse computations.
2343Combining 1. and 2. gives a good reason for expecting that a method
2344like block Jacobi could become very powerful in the future, even
2345though currently for these examples it is quite slow.
2348<a name="step_63-Possibilitiesforextensions"></a><h3> Possibilities for extensions </h3>
2351<a name="step_63-ConstantiterationsforQsub5sub"></a><h4> Constant iterations for Q<sub>5</sub> </h4>
2354Change the number of smoothing steps and the smoother relaxation
2355parameter (set in <code>Smoother::AdditionalData()</code> inside
2356<code>create_smoother()</code>, only necessary for point smoothers) so
2357that we maintain a constant number of iterations for a @f$Q_5@f$ element.
2359<a name="step_63-Effectivenessofrenumberingforchangingepsilon"></a><h4> Effectiveness of renumbering for changing epsilon </h4>
2362Increase/decrease the parameter "Epsilon" in the `.prm` files of the
2363multiplicative methods and observe for which values renumbering no
2364longer influences convergence speed.
2366<a name="step_63-Meshadaptivity"></a><h4> Mesh adaptivity </h4>
2369The code is set up to work correctly with an adaptively refined mesh (the
2370interface matrices are created and set). Devise a suitable refinement
2371criterium or try the KellyErrorEstimator class.
2374<a name="step_63-PlainProg"></a>
2375<h1> The plain program</h1>
2376@include "step-63.cc"
numbers::NumberTraits< Number >::real_type norm() const
std::conditional_t< rank_==1, std::array< Number, dim >, std::array< Tensor< rank_ - 1, dim, Number >, dim > > values
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertThrow(cond, exc)
void mesh_loop(const CellIteratorType &begin, const CellIteratorType &end, const CellWorkerFunctionType &cell_worker, const CopierType &copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const AssembleFlags flags=assemble_own_cells, const BoundaryWorkerFunctionType &boundary_worker=BoundaryWorkerFunctionType(), const FaceWorkerFunctionType &face_worker=FaceWorkerFunctionType(), const unsigned int queue_length=2 *MultithreadInfo::n_threads(), const unsigned int chunk_size=8)
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
#define DEAL_II_NOT_IMPLEMENTED()
std::vector< index_type > data
Expression coth(const Expression &x)
void downstream(DoFHandler< dim, spacedim > &dof_handler, const Tensor< 1, spacedim > &direction, const bool dof_wise_renumbering=false)
void random(DoFHandler< dim, spacedim > &dof_handler)
void refine(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold, const unsigned int max_to_mark=numbers::invalid_unsigned_int)
@ matrix
Contents is actually a matrix.
@ symmetric
Matrix is symmetric.
constexpr types::blas_int zero
constexpr types::blas_int one
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
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)
SymmetricTensor< 2, dim, Number > epsilon(const Tensor< 2, dim, Number > &Grad_u)
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
int(&) functions(const void *v1, const void *v2)
::VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
DEAL_II_HOST constexpr Number trace(const SymmetricTensor< 2, dim2, Number > &)