523 *
"Diffusion parameter");
525 * prm.declare_entry(
"Fe degree",
528 *
"Finite Element degree");
529 * prm.declare_entry(
"Smoother type",
532 *
"Select smoother: SOR|Jacobi|block SOR|block Jacobi");
533 * prm.declare_entry(
"Smoothing steps",
536 *
"Number of smoothing steps");
541 *
"Select DoF renumbering: none|downstream|upstream|random");
542 * prm.declare_entry(
"With streamline diffusion",
545 *
"Enable streamline diffusion stabilization: true|false");
546 * prm.declare_entry(
"Output",
549 *
"Generate graphical output: true|false");
552 *
if (prm_filename.empty())
556 *
false, ExcMessage(
"Please pass a .prm file as the first argument!"));
559 * prm.parse_input(prm_filename);
561 *
epsilon = prm.get_double(
"Epsilon");
562 * fe_degree = prm.get_integer(
"Fe degree");
563 * smoother_type = prm.get(
"Smoother type");
564 * smoothing_steps = prm.get_integer(
"Smoothing steps");
566 *
const std::string renumbering = prm.get(
"DoF renumbering");
567 *
if (renumbering ==
"none")
568 * dof_renumbering = DoFRenumberingStrategy::none;
569 *
else if (renumbering ==
"downstream")
570 * dof_renumbering = DoFRenumberingStrategy::downstream;
571 *
else if (renumbering ==
"upstream")
572 * dof_renumbering = DoFRenumberingStrategy::upstream;
573 *
else if (renumbering ==
"random")
574 * dof_renumbering = DoFRenumberingStrategy::random;
577 * ExcMessage(
"The <DoF renumbering> parameter has "
578 *
"an invalid value."));
580 * with_streamline_diffusion = prm.get_bool(
"With streamline diffusion");
581 * output = prm.get_bool(
"Output");
588 * <a name=
"step_63-Cellpermutations"></a>
589 * <h3>Cell permutations</h3>
593 * The ordering in which cells and degrees of freedom are traversed
594 * will play a role in the speed of convergence
for multiplicative
595 * methods. Here we define
functions which
return a specific ordering
596 * of cells to be used by the block smoothers.
600 * For each type of cell ordering, we define a function
for the
601 * active mesh and one
for a
level mesh (i.e.,
for the cells at one
602 *
level of a multigrid hierarchy). While the only reordering
603 * necessary
for solving the system will be on the
level meshes, we
604 * include the active reordering
for visualization purposes in
610 * array with all of the relevant cells that we then sort in
611 *
downstream direction
using a
"comparator" object. The output of
612 * the
functions is then simply an array of the indices of the cells
613 * in the just computed order.
617 * std::vector<unsigned int>
620 *
const unsigned int level)
622 * std::vector<typename DoFHandler<dim>::level_cell_iterator> ordered_cells;
623 * ordered_cells.reserve(dof_handler.get_triangulation().n_cells(
level));
624 *
for (
const auto &cell : dof_handler.cell_iterators_on_level(
level))
625 * ordered_cells.push_back(cell);
628 * CompareDownstream<typename DoFHandler<dim>::level_cell_iterator, dim>
629 * comparator(direction);
630 * std::sort(ordered_cells.begin(), ordered_cells.end(), comparator);
632 * std::vector<unsigned> ordered_indices;
633 * ordered_indices.reserve(dof_handler.get_triangulation().n_cells(
level));
635 *
for (
const auto &cell : ordered_cells)
636 * ordered_indices.push_back(cell->
index());
638 *
return ordered_indices;
644 * std::vector<unsigned int>
648 * std::vector<typename DoFHandler<dim>::active_cell_iterator> ordered_cells;
649 * ordered_cells.reserve(dof_handler.get_triangulation().n_active_cells());
650 *
for (
const auto &cell : dof_handler.active_cell_iterators())
651 * ordered_cells.push_back(cell);
654 * CompareDownstream<typename DoFHandler<dim>::active_cell_iterator, dim>
655 * comparator(direction);
656 * std::sort(ordered_cells.begin(), ordered_cells.end(), comparator);
658 * std::vector<unsigned int> ordered_indices;
659 * ordered_indices.reserve(dof_handler.get_triangulation().n_active_cells());
661 *
for (
const auto &cell : ordered_cells)
662 * ordered_indices.push_back(cell->
index());
664 *
return ordered_indices;
671 * spirit in that they
first put information about all cells into an
672 * array. But then, instead of sorting them, they shuffle the
673 * elements randomly
using the facilities
C++ offers to generate
674 *
random numbers. The way
this is done is by iterating over all
675 * elements of the array, drawing a
random number
for another
676 * element before that, and then exchanging these elements. The
677 * result is a
random shuffle of the elements of the array.
681 * std::vector<unsigned int>
683 *
const unsigned int level)
685 * std::vector<unsigned int> ordered_cells;
686 * ordered_cells.reserve(dof_handler.get_triangulation().n_cells(
level));
687 *
for (
const auto &cell : dof_handler.cell_iterators_on_level(
level))
688 * ordered_cells.push_back(cell->
index());
690 * std::mt19937 random_number_generator;
691 * std::shuffle(ordered_cells.begin(),
692 * ordered_cells.end(),
693 * random_number_generator);
695 *
return ordered_cells;
701 * std::vector<unsigned int>
704 * std::vector<unsigned int> ordered_cells;
705 * ordered_cells.reserve(dof_handler.get_triangulation().n_active_cells());
706 *
for (
const auto &cell : dof_handler.active_cell_iterators())
707 * ordered_cells.push_back(cell->
index());
709 * std::mt19937 random_number_generator;
710 * std::shuffle(ordered_cells.begin(),
711 * ordered_cells.end(),
712 * random_number_generator);
714 *
return ordered_cells;
721 * <a name=
"step_63-Righthandsideandboundaryvalues"></a>
722 * <h3>Right-hand side and boundary
values</h3>
726 * The problem solved in
this tutorial is an adaptation of Ex. 3.1.3 found
728 * href=
"https://global.oup.com/academic/product/finite-elements-and-fast-iterative-solvers-9780199678808">
729 * Finite Elements and Fast Iterative Solvers: with Applications in
730 * Incompressible Fluid Dynamics by Elman, Silvester, and Wathen</a>. The
731 * main difference being that we add a hole in the
center of our domain with
732 * zero Dirichlet boundary conditions.
736 * For a complete description, we need classes that implement the
737 * zero right-hand side
first (we could of course have just used
742 * class RightHandSide : public
Function<dim>
746 * const unsigned
int component = 0) const override;
748 *
virtual void value_list(
const std::vector<
Point<dim>> &points,
749 * std::vector<double> &values,
750 *
const unsigned int component = 0)
const override;
756 *
double RightHandSide<dim>::value(
const Point<dim> &,
757 *
const unsigned int component)
const
759 *
Assert(component == 0, ExcIndexRange(component, 0, 1));
768 *
void RightHandSide<dim>::value_list(
const std::vector<
Point<dim>> &points,
769 * std::vector<double> &values,
770 *
const unsigned int component)
const
774 *
for (
unsigned int i = 0; i < points.size(); ++i)
775 * values[i] = RightHandSide<dim>::value(points[i], component);
781 * We also have Dirichlet boundary conditions. On a connected portion of the
782 * outer, square boundary we set the
value to 1, and we set the
value to 0
783 * everywhere
else (including the inner, circular boundary):
787 *
class BoundaryValues :
public Function<dim>
791 *
const unsigned int component = 0)
const override;
794 * std::vector<double> &values,
795 *
const unsigned int component = 0)
const override;
801 *
double BoundaryValues<dim>::value(
const Point<dim> &p,
802 *
const unsigned int component)
const
804 *
Assert(component == 0, ExcIndexRange(component, 0, 1));
809 * Set boundary to 1
if @f$x=1@f$, or
if @f$x>0.5@f$ and @f$y=-1@f$.
812 *
if (std::fabs(p[0] - 1) < 1e-8 ||
813 * (std::fabs(p[1] + 1) < 1e-8 && p[0] >= 0.5))
826 *
void BoundaryValues<dim>::value_list(
const std::vector<
Point<dim>> &points,
827 * std::vector<double> &values,
828 *
const unsigned int component)
const
832 *
for (
unsigned int i = 0; i < points.size(); ++i)
833 * values[i] = BoundaryValues<dim>::value(points[i], component);
841 * <a name=
"step_63-Streamlinediffusionimplementation"></a>
842 * <h3>Streamline diffusion implementation</h3>
846 * The streamline diffusion method has a stabilization
constant that
847 * we need to be able to compute. The choice of how
this parameter
848 * is computed is taken from <a
849 * href=
"https://link.springer.com/chapter/10.1007/978-3-540-34288-5_27">On
850 * Discontinuity-Capturing Methods
for Convection-Diffusion
851 * Equations by Volker John and Petr Knobloch</a>.
855 *
double compute_stabilization_delta(
const double hk,
860 *
const double Peclet = dir.norm() * hk / (2.0 *
eps * pk);
861 *
const double coth =
864 *
return hk / (2.0 * dir.norm() * pk) * (coth - 1.0 / Peclet);
871 * <a name=
"step_63-codeAdvectionProblemcodeclass"></a>
872 * <h3><code>AdvectionProblem</code>
class</h3>
876 * This is the main
class of the program, and should look very similar to
877 * @ref step_16
"step-16". The major difference is that, since we are defining our multigrid
878 * smoother at runtime, we choose to define a function `create_smoother()` and
879 * a
class object `mg_smoother` which is a `
std::unique_ptr` to a smoother
880 * that is derived from
MGSmoother. Note that for smoothers derived from
882 * This will contain information about the cell ordering and the method of
883 * inverting cell matrices.
890 * class AdvectionProblem
893 * AdvectionProblem(
const Settings &settings);
897 *
void setup_system();
899 *
template <
class IteratorType>
900 *
void assemble_cell(
const IteratorType &cell,
901 * ScratchData<dim> &scratch_data,
902 * CopyData ©_data);
903 *
void assemble_system_and_multigrid();
905 *
void setup_smoother();
908 *
void refine_grid();
909 *
void output_results(
const unsigned int cycle)
const;
936 * std::unique_ptr<MGSmoother<Vector<double>>> mg_smoother;
938 *
using SmootherType =
940 *
using SmootherAdditionalDataType = SmootherType::AdditionalData;
953 * AdvectionProblem<dim>::AdvectionProblem(
const Settings &settings)
956 * , fe(settings.fe_degree)
957 * , mapping(settings.fe_degree)
958 * , settings(settings)
971 * <a name=
"step_63-codeAdvectionProblemsetup_systemcode"></a>
972 * <h4><code>AdvectionProblem::setup_system()</code></h4>
981 * We could renumber the active DoFs with the
DoFRenumbering class,
982 * but the smoothers only act on multigrid levels and as such,
this
983 * would not matter
for the computations. Instead, we will renumber the
984 * DoFs on each multigrid
level below.
988 *
void AdvectionProblem<dim>::setup_system()
992 * dof_handler.distribute_dofs(fe);
994 * solution.reinit(dof_handler.n_dofs());
995 * system_rhs.reinit(dof_handler.n_dofs());
997 * constraints.clear();
1001 * mapping, dof_handler, 0, BoundaryValues<dim>(), constraints);
1003 * mapping, dof_handler, 1, BoundaryValues<dim>(), constraints);
1004 * constraints.close();
1012 * sparsity_pattern.copy_from(dsp);
1013 * system_matrix.reinit(sparsity_pattern);
1015 * dof_handler.distribute_mg_dofs();
1019 * Having enumerated the global degrees of freedom as well as (in
1020 * the last line above) the
level degrees of freedom, let us
1021 * renumber the
level degrees of freedom to get a better smoother
1022 * as explained in the introduction. The
first block below
1024 * direction
if needed. This is only necessary
for point smoothers
1025 * (SOR and Jacobi) as the block smoothers operate on cells (see
1026 * `create_smoother()`). The blocks below then also implement
1030 * if (settings.smoother_type ==
"SOR" || settings.smoother_type ==
"Jacobi")
1032 *
if (settings.dof_renumbering ==
1033 * Settings::DoFRenumberingStrategy::downstream ||
1034 * settings.dof_renumbering ==
1035 * Settings::DoFRenumberingStrategy::upstream)
1038 * (settings.dof_renumbering ==
1039 * Settings::DoFRenumberingStrategy::upstream ?
1042 * advection_direction;
1050 *
else if (settings.dof_renumbering ==
1051 * Settings::DoFRenumberingStrategy::random)
1062 * The rest of the function just sets up data structures. The last
1063 * lines of the code below is unlike the other GMG tutorials, as
1064 * it sets up both the
interface in and out matrices. We need this
1068 * mg_constrained_dofs.clear();
1069 * mg_constrained_dofs.initialize(dof_handler);
1071 * mg_constrained_dofs.make_zero_boundary_constraints(dof_handler, {0, 1});
1073 * mg_matrices.resize(0, n_levels - 1);
1074 * mg_matrices.clear_elements();
1075 * mg_interface_in.resize(0, n_levels - 1);
1076 * mg_interface_in.clear_elements();
1077 * mg_interface_out.resize(0, n_levels - 1);
1078 * mg_interface_out.clear_elements();
1079 * mg_sparsity_patterns.resize(0, n_levels - 1);
1080 * mg_interface_sparsity_patterns.resize(0, n_levels - 1);
1086 * dof_handler.n_dofs(
level));
1088 * mg_sparsity_patterns[
level].copy_from(dsp);
1089 * mg_matrices[
level].reinit(mg_sparsity_patterns[
level]);
1093 * dof_handler.n_dofs(
level));
1095 * mg_constrained_dofs,
1098 * mg_interface_sparsity_patterns[
level].copy_from(dsp);
1100 * mg_interface_in[
level].reinit(mg_interface_sparsity_patterns[
level]);
1101 * mg_interface_out[
level].reinit(mg_interface_sparsity_patterns[
level]);
1110 * <a name=
"step_63-codeAdvectionProblemassemble_cellcode"></a>
1111 * <h4><code>AdvectionProblem::assemble_cell()</code></h4>
1115 * Here we define the assembly of the linear system on each cell to
1116 * be used by the
mesh_loop() function below. This one function
1117 * assembles the cell matrix for either an active or a
level cell
1118 * (whatever it is passed as its
first argument), and only assembles
1119 * a right-hand side if called with an active cell.
1125 * template <
int dim>
1126 * template <class IteratorType>
1127 *
void AdvectionProblem<dim>::assemble_cell(const IteratorType &cell,
1128 * ScratchData<dim> &scratch_data,
1129 * CopyData ©_data)
1131 * copy_data.level = cell->level();
1133 *
const unsigned int dofs_per_cell =
1134 * scratch_data.fe_values.get_fe().n_dofs_per_cell();
1135 * copy_data.dofs_per_cell = dofs_per_cell;
1136 * copy_data.cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
1138 *
const unsigned int n_q_points =
1139 * scratch_data.fe_values.get_quadrature().size();
1141 *
if (cell->is_level_cell() ==
false)
1142 * copy_data.cell_rhs.reinit(dofs_per_cell);
1144 * copy_data.local_dof_indices.resize(dofs_per_cell);
1145 * cell->get_active_or_mg_dof_indices(copy_data.local_dof_indices);
1147 * scratch_data.fe_values.reinit(cell);
1149 * RightHandSide<dim> right_hand_side;
1150 * std::vector<double> rhs_values(n_q_points);
1152 * right_hand_side.value_list(scratch_data.fe_values.get_quadrature_points(),
1157 * If we are
using streamline diffusion we must add its contribution
1158 * to both the cell
matrix and the cell right-hand side. If we are not
1159 *
using streamline diffusion, setting @f$\delta=0@f$ negates
this contribution
1160 * below and we are left with the standard, Galerkin finite element
1164 *
const double delta = (settings.with_streamline_diffusion ?
1165 * compute_stabilization_delta(cell->diameter(),
1167 * advection_direction,
1168 * settings.fe_degree) :
1171 *
for (
unsigned int q_point = 0; q_point < n_q_points; ++q_point)
1172 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1174 *
for (
unsigned int j = 0; j < dofs_per_cell; ++j)
1178 * The assembly of the local
matrix has two parts. First
1179 * the Galerkin contribution:
1182 * copy_data.cell_matrix(i, j) +=
1183 * (settings.epsilon *
1184 * scratch_data.fe_values.shape_grad(i, q_point) *
1185 * scratch_data.fe_values.shape_grad(j, q_point) *
1186 * scratch_data.fe_values.JxW(q_point)) +
1187 * (scratch_data.fe_values.shape_value(i, q_point) *
1188 * (advection_direction *
1189 * scratch_data.fe_values.shape_grad(j, q_point)) *
1190 * scratch_data.fe_values.JxW(q_point))
1193 * and then the streamline diffusion contribution:
1197 * (advection_direction *
1198 * scratch_data.fe_values.shape_grad(j, q_point)) *
1199 * (advection_direction *
1200 * scratch_data.fe_values.shape_grad(i, q_point)) *
1201 * scratch_data.fe_values.JxW(q_point) -
1202 * delta * settings.epsilon *
1203 *
trace(scratch_data.fe_values.shape_hessian(j, q_point)) *
1204 * (advection_direction *
1205 * scratch_data.fe_values.shape_grad(i, q_point)) *
1206 * scratch_data.fe_values.JxW(q_point);
1208 *
if (cell->is_level_cell() ==
false)
1212 * The same applies to the right hand side. First the
1213 * Galerkin contribution:
1216 * copy_data.cell_rhs(i) +=
1217 * scratch_data.fe_values.shape_value(i, q_point) *
1218 * rhs_values[q_point] * scratch_data.fe_values.JxW(q_point)
1221 * and then the streamline diffusion contribution:
1224 * + delta * rhs_values[q_point] * advection_direction *
1225 * scratch_data.fe_values.shape_grad(i, q_point) *
1226 * scratch_data.fe_values.JxW(q_point);
1235 * <a name=
"step_63-codeAdvectionProblemassemble_system_and_multigridcode"></a>
1236 * <h4><code>AdvectionProblem::assemble_system_and_multigrid()</code></h4>
1241 * system_matrix, system_rhs, and all mg_matrices for us.
1247 * template <
int dim>
1248 *
void AdvectionProblem<dim>::assemble_system_and_multigrid()
1250 *
const auto cell_worker_active =
1251 * [&](
const decltype(dof_handler.begin_active()) &cell,
1252 * ScratchData<dim> &scratch_data,
1253 * CopyData ©_data) {
1254 * this->assemble_cell(cell, scratch_data, copy_data);
1257 *
const auto copier_active = [&](
const CopyData ©_data) {
1258 * constraints.distribute_local_to_global(copy_data.cell_matrix,
1259 * copy_data.cell_rhs,
1260 * copy_data.local_dof_indices,
1267 * dof_handler.end(),
1268 * cell_worker_active,
1270 * ScratchData<dim>(fe, fe.degree + 1),
1276 * Unlike the constraints
for the active
level, we choose to create
1277 * constraint objects
for each multigrid
level local to
this function
1278 * since they are never needed elsewhere in the program.
1281 * std::vector<AffineConstraints<double>> boundary_constraints(
1286 * boundary_constraints[
level].reinit(
1287 * dof_handler.locally_owned_mg_dofs(
level),
1291 * mg_constrained_dofs.get_refinement_edge_indices(
level))
1292 * boundary_constraints[
level].constrain_dof_to_zero(dof_index);
1294 * mg_constrained_dofs.get_boundary_indices(
level))
1295 * boundary_constraints[
level].constrain_dof_to_zero(dof_index);
1296 * boundary_constraints[
level].close();
1299 *
const auto cell_worker_mg =
1300 * [&](
const decltype(dof_handler.begin_mg()) &cell,
1301 * ScratchData<dim> &scratch_data,
1302 * CopyData ©_data) {
1303 * this->assemble_cell(cell, scratch_data, copy_data);
1306 *
const auto copier_mg = [&](
const CopyData ©_data) {
1307 * boundary_constraints[copy_data.level].distribute_local_to_global(
1308 * copy_data.cell_matrix,
1309 * copy_data.local_dof_indices,
1310 * mg_matrices[copy_data.level]);
1314 * If @f$(i,j)@f$ is an `interface_out` dof pair, then @f$(j,i)@f$ is an
1315 * `interface_in` dof pair. Note: For `interface_in`, we load
1316 * the
transpose of the
interface entries, i.
e., the entry for
1317 * dof pair @f$(j,i)@f$ is stored in `interface_in(i,j)`. This is an
1318 * optimization
for the
symmetric case which allows only one
1319 *
matrix to be used when setting the edge_matrices in
1320 * solve(). Here, however, since our problem is non-
symmetric,
1321 * we must store both `interface_in` and `interface_out`
1325 *
for (
unsigned int i = 0; i < copy_data.dofs_per_cell; ++i)
1326 *
for (
unsigned int j = 0; j < copy_data.dofs_per_cell; ++j)
1327 *
if (mg_constrained_dofs.is_interface_matrix_entry(
1329 * copy_data.local_dof_indices[i],
1330 * copy_data.local_dof_indices[j]))
1332 * mg_interface_out[copy_data.level].add(
1333 * copy_data.local_dof_indices[i],
1334 * copy_data.local_dof_indices[j],
1335 * copy_data.cell_matrix(i, j));
1336 * mg_interface_in[copy_data.level].add(
1337 * copy_data.local_dof_indices[i],
1338 * copy_data.local_dof_indices[j],
1339 * copy_data.cell_matrix(j, i));
1344 * dof_handler.end_mg(),
1347 * ScratchData<dim>(fe, fe.degree + 1),
1356 * <a name=
"step_63-codeAdvectionProblemsetup_smoothercode"></a>
1357 * <h4><code>AdvectionProblem::setup_smoother()</code></h4>
1361 * Next, we set up the smoother based on the settings in the `.prm` file. The
1362 * two options that are of significance is the number of pre- and
1363 * post-smoothing steps on each
level of the multigrid v-cycle and the
1364 * relaxation parameter.
1368 * Since multiplicative methods tend to be more powerful than additive method,
1369 * fewer smoothing steps are required to see convergence
independent of mesh
1370 * size. The same holds
for block smoothers over
point smoothers. This is
1371 * reflected in the choice
for the number of smoothing steps
for each type of
1376 * The relaxation parameter
for point smoothers is chosen based on trial and
1377 * error, and reflects
values necessary to keep the iteration counts in
1378 * the GMRES solve
constant (or as close as possible) as we
refine the mesh.
1379 * The two
values given
for both
"Jacobi" and
"SOR" in the `.prm` files are
1380 *
for degree 1 and degree 3 finite elements. If the user wants to change to
1381 * another degree, they may need to adjust these
numbers. For block smoothers,
1382 *
this parameter has a more straightforward interpretation, namely that
for
1383 * additive methods in 2
d, a DoF can have a repeated contribution from up to 4
1384 * cells, therefore we must relax these methods by 0.25 to compensate. This is
1385 * not an issue
for multiplicative methods as each cell
's inverse application
1386 * carries new information to all its DoFs.
1390 * Finally, as mentioned above, the point smoothers only operate on DoFs, and
1391 * the block smoothers on cells, so only the block smoothers need to be given
1392 * information regarding cell orderings. DoF ordering for point smoothers has
1393 * already been taken care of in `setup_system()`.
1399 * template <int dim>
1400 * void AdvectionProblem<dim>::setup_smoother()
1402 * if (settings.smoother_type == "SOR")
1404 * using Smoother = PreconditionSOR<SparseMatrix<double>>;
1407 * std::make_unique<MGSmootherPrecondition<SparseMatrix<double>,
1409 * Vector<double>>>();
1410 * smoother->initialize(mg_matrices,
1411 * Smoother::AdditionalData(fe.degree == 1 ? 1.0 :
1413 * smoother->set_steps(settings.smoothing_steps);
1414 * mg_smoother = std::move(smoother);
1416 * else if (settings.smoother_type == "Jacobi")
1418 * using Smoother = PreconditionJacobi<SparseMatrix<double>>;
1420 * std::make_unique<MGSmootherPrecondition<SparseMatrix<double>,
1422 * Vector<double>>>();
1423 * smoother->initialize(mg_matrices,
1424 * Smoother::AdditionalData(fe.degree == 1 ? 0.6667 :
1426 * smoother->set_steps(settings.smoothing_steps);
1427 * mg_smoother = std::move(smoother);
1429 * else if (settings.smoother_type == "block SOR" ||
1430 * settings.smoother_type == "block Jacobi")
1432 * smoother_data.resize(0, triangulation.n_levels() - 1);
1434 * for (unsigned int level = 0; level < triangulation.n_levels(); ++level)
1436 * DoFTools::make_cell_patches(smoother_data[level].block_list,
1440 * smoother_data[level].relaxation =
1441 * (settings.smoother_type == "block SOR" ? 1.0 : 0.25);
1442 * smoother_data[level].inversion = PreconditionBlockBase<double>::svd;
1444 * std::vector<unsigned int> ordered_indices;
1445 * switch (settings.dof_renumbering)
1447 * case Settings::DoFRenumberingStrategy::downstream:
1449 * create_downstream_cell_ordering(dof_handler,
1450 * advection_direction,
1454 * case Settings::DoFRenumberingStrategy::upstream:
1456 * create_downstream_cell_ordering(dof_handler,
1457 * -1.0 * advection_direction,
1461 * case Settings::DoFRenumberingStrategy::random:
1463 * create_random_cell_ordering(dof_handler, level);
1466 * case Settings::DoFRenumberingStrategy::none:
1470 * AssertThrow(false, ExcNotImplemented());
1474 * smoother_data[level].order =
1475 * std::vector<std::vector<unsigned int>>(1, ordered_indices);
1478 * if (settings.smoother_type == "block SOR")
1480 * auto smoother = std::make_unique<MGSmootherPrecondition<
1481 * SparseMatrix<double>,
1482 * RelaxationBlockSOR<SparseMatrix<double>, double, Vector<double>>,
1483 * Vector<double>>>();
1484 * smoother->initialize(mg_matrices, smoother_data);
1485 * smoother->set_steps(settings.smoothing_steps);
1486 * mg_smoother = std::move(smoother);
1488 * else if (settings.smoother_type == "block Jacobi")
1490 * auto smoother = std::make_unique<
1491 * MGSmootherPrecondition<SparseMatrix<double>,
1492 * RelaxationBlockJacobi<SparseMatrix<double>,
1495 * Vector<double>>>();
1496 * smoother->initialize(mg_matrices, smoother_data);
1497 * smoother->set_steps(settings.smoothing_steps);
1498 * mg_smoother = std::move(smoother);
1502 * AssertThrow(false, ExcNotImplemented());
1509 * <a name="step_63-codeAdvectionProblemsolvecode"></a>
1510 * <h4><code>AdvectionProblem::solve()</code></h4>
1514 * Before we can solve the system, we must first set up the multigrid
1515 * preconditioner. This requires the setup of the transfer between levels,
1516 * the coarse matrix solver, and the smoother. This setup follows almost
1517 * identically to @ref step_16 "step-16", the main difference being the various smoothers
1518 * defined above and the fact that we need different interface edge matrices
1519 * for in and out since our problem is non-symmetric. (In reality, for this
1520 * tutorial these interface matrices are empty since we are only using global
1521 * refinement, and thus have no refinement edges. However, we have still
1522 * included both here since if one made the simple switch to an adaptively
1523 * refined method, the program would still run correctly.)
1527 * The last thing to note is that since our problem is non-symmetric, we must
1528 * use an appropriate Krylov subspace method. We choose here to
1529 * use GMRES since it offers the guarantee of residual reduction in each
1530 * iteration. The major disadvantage of GMRES is that, for each iteration,
1531 * the number of stored temporary vectors increases by one, and one also needs
1532 * to compute a scalar product with all previously stored vectors. This is
1533 * rather expensive. This requirement is relaxed by using the restarted GMRES
1534 * method which puts a cap on the number of vectors we are required to store
1535 * at any one time (here we restart after 50 temporary vectors, or 48
1536 * iterations). This then has the disadvantage that we lose information we
1537 * have gathered throughout the iteration and therefore we could see slower
1538 * convergence. As a consequence, where to restart is a question of balancing
1539 * memory consumption, CPU effort, and convergence speed.
1540 * However, the goal of this tutorial is to have very low
1541 * iteration counts by using a powerful GMG preconditioner, so we have picked
1542 * the restart length such that all of the results shown below converge prior
1543 * to restart happening, and thus we have a standard GMRES method. If the user
1544 * is interested, another suitable method offered in deal.II would be
1551 * template <int dim>
1552 * void AdvectionProblem<dim>::solve()
1554 * const unsigned int max_iters = 200;
1555 * const double solve_tolerance = 1e-8 * system_rhs.l2_norm();
1556 * SolverControl solver_control(max_iters, solve_tolerance, true, true);
1557 * solver_control.enable_history_data();
1559 * using Transfer = MGTransferPrebuilt<Vector<double>>;
1560 * Transfer mg_transfer(mg_constrained_dofs);
1561 * mg_transfer.build(dof_handler);
1563 * FullMatrix<double> coarse_matrix;
1564 * coarse_matrix.copy_from(mg_matrices[0]);
1565 * MGCoarseGridHouseholder<double, Vector<double>> coarse_grid_solver;
1566 * coarse_grid_solver.initialize(coarse_matrix);
1570 * mg_matrix.initialize(mg_matrices);
1571 * mg_interface_matrix_in.initialize(mg_interface_in);
1572 * mg_interface_matrix_out.initialize(mg_interface_out);
1574 * Multigrid<Vector<double>> mg(
1575 * mg_matrix, coarse_grid_solver, mg_transfer, *mg_smoother, *mg_smoother);
1576 * mg.set_edge_matrices(mg_interface_matrix_out, mg_interface_matrix_in);
1578 * PreconditionMG<dim, Vector<double>, Transfer> preconditioner(dof_handler,
1582 * std::cout << " Solving with GMRES to tol " << solve_tolerance << "..."
1584 * SolverGMRES<Vector<double>> solver(
1585 * solver_control, SolverGMRES<Vector<double>>::AdditionalData(50, true));
1589 * solver.solve(system_matrix, solution, system_rhs, preconditioner);
1592 * std::cout << " converged in " << solver_control.last_step()
1594 * << " in " << time.last_wall_time() << " seconds " << std::endl;
1596 * constraints.distribute(solution);
1598 * mg_smoother.release();
1605 * <a name="step_63-codeAdvectionProblemoutput_resultscode"></a>
1606 * <h4><code>AdvectionProblem::output_results()</code></h4>
1610 * The final function of interest generates graphical output.
1611 * Here we output the solution and cell ordering in a .vtu format.
1615 * At the top of the function, we generate an index for each cell to
1616 * visualize the ordering used by the smoothers. Note that we do
1617 * this only for the active cells instead of the levels, where the
1618 * smoothers are actually used. For the point smoothers we renumber
1619 * DoFs instead of cells, so this is only an approximation of what
1620 * happens in reality. Finally, the random ordering is not the
1621 * random ordering we actually use (see `create_smoother()` for that).
1625 * The (integer) ordering of cells is then copied into a (floating
1626 * point) vector for graphical output.
1629 * template <int dim>
1630 * void AdvectionProblem<dim>::output_results(const unsigned int cycle) const
1632 * const unsigned int n_active_cells = triangulation.n_active_cells();
1633 * Vector<double> cell_indices(n_active_cells);
1635 * std::vector<unsigned int> ordered_indices;
1636 * switch (settings.dof_renumbering)
1638 * case Settings::DoFRenumberingStrategy::downstream:
1640 * create_downstream_cell_ordering(dof_handler, advection_direction);
1643 * case Settings::DoFRenumberingStrategy::upstream:
1645 * create_downstream_cell_ordering(dof_handler,
1646 * -1.0 * advection_direction);
1649 * case Settings::DoFRenumberingStrategy::random:
1650 * ordered_indices = create_random_cell_ordering(dof_handler);
1653 * case Settings::DoFRenumberingStrategy::none:
1654 * ordered_indices.resize(n_active_cells);
1655 * for (unsigned int i = 0; i < n_active_cells; ++i)
1656 * ordered_indices[i] = i;
1660 * AssertThrow(false, ExcNotImplemented());
1664 * for (unsigned int i = 0; i < n_active_cells; ++i)
1665 * cell_indices(ordered_indices[i]) = static_cast<double>(i);
1670 * The remainder of the function is then straightforward, given
1671 * previous tutorial programs:
1674 * DataOut<dim> data_out;
1675 * data_out.attach_dof_handler(dof_handler);
1676 * data_out.add_data_vector(solution, "solution");
1677 * data_out.add_data_vector(cell_indices, "cell_index");
1678 * data_out.build_patches();
1680 * const std::string filename =
1681 * "solution-" + Utilities::int_to_string(cycle) + ".vtu";
1682 * std::ofstream output(filename);
1683 * data_out.write_vtu(output);
1690 * <a name="step_63-codeAdvectionProblemruncode"></a>
1691 * <h4><code>AdvectionProblem::run()</code></h4>
1695 * As in most tutorials, this function creates/refines the mesh and calls
1696 * the various functions defined above to set up, assemble, solve, and output
1701 * In cycle zero, we generate the mesh for the on the square
1702 * <code>[-1,1]^dim</code> with a hole of radius 3/10 units centered
1703 * at the origin. For objects with `manifold_id` equal to one
1704 * (namely, the faces adjacent to the hole), we assign a spherical
1711 * template <int dim>
1712 * void AdvectionProblem<dim>::run()
1714 * for (unsigned int cycle = 0; cycle < (settings.fe_degree == 1 ? 7u : 5u);
1717 * std::cout << " Cycle " << cycle << ':
' << std::endl;
1721 * GridGenerator::hyper_cube_with_cylindrical_hole(triangulation,
1725 * const SphericalManifold<dim> manifold_description(Point<dim>(0, 0));
1726 * triangulation.set_manifold(1, manifold_description);
1729 * triangulation.refine_global();
1733 * std::cout << " Number of active cells: "
1734 * << triangulation.n_active_cells() << " ("
1735 * << triangulation.n_levels() << " levels)" << std::endl;
1736 * std::cout << " Number of degrees of freedom: "
1737 * << dof_handler.n_dofs() << std::endl;
1739 * assemble_system_and_multigrid();
1743 * if (settings.output)
1744 * output_results(cycle);
1746 * std::cout << std::endl;
1749 * } // namespace Step63
1755 * <a name="step_63-Thecodemaincodefunction"></a>
1756 * <h3>The <code>main</code> function</h3>
1760 * Finally, the main function is like most tutorials. The only
1761 * interesting bit is that we require the user to pass a `.prm` file
1762 * as a sole command line argument. If no parameter file is given, the
1763 * program will output the contents of a sample parameter file with
1764 * all default values to the screen that the user can then copy and
1765 * paste into their own `.prm` file.
1771 * int main(int argc, char *argv[])
1775 * Step63::Settings settings;
1776 * settings.get_parameters((argc > 1) ? (argv[1]) : "");
1778 * Step63::AdvectionProblem<2> advection_problem_2d(settings);
1779 * advection_problem_2d.run();
1781 * catch (std::exception &exc)
1783 * std::cerr << std::endl
1785 * << "----------------------------------------------------"
1787 * std::cerr << "Exception on processing: " << std::endl
1788 * << exc.what() << std::endl
1789 * << "Aborting!" << std::endl
1790 * << "----------------------------------------------------"
1796 * std::cerr << std::endl
1798 * << "----------------------------------------------------"
1800 * std::cerr << "Unknown exception!" << std::endl
1801 * << "Aborting!" << std::endl
1802 * << "----------------------------------------------------"
1810<a name="step_63-Results"></a><h1>Results</h1>
1813<a name="step_63-GMRESIterationNumbers"></a><h3> GMRES Iteration Numbers </h3>
1816The major advantage for GMG is that it is an @f$\mathcal{O}(n)@f$ method,
1817that is, the complexity of the problem increases linearly with the
1818problem size. To show then that the linear solver presented in this
1819tutorial is in fact @f$\mathcal{O}(n)@f$, all one needs to do is show that
1820the iteration counts for the GMRES solve stay roughly constant as we
1823Each of the following tables gives the GMRES iteration counts to reduce the
1824initial residual by a factor of @f$10^8@f$. We selected a sufficient number of smoothing steps
1825(based on the method) to get iteration numbers independent of mesh size. As
1826can be seen from the tables below, the method is indeed @f$\mathcal{O}(n)@f$.
1828<a name="step_63-DoFCellRenumbering"></a><h4> DoF/Cell Renumbering </h4>
1831The point-wise smoothers ("Jacobi" and "SOR") get applied in the order the
1832DoFs are numbered on each level. We can influence this using the
1833DoFRenumbering namespace. The block smoothers are applied based on the
1834ordering we set in `setup_smoother()`. We can visualize this numbering. The
1835following pictures show the cell numbering of the active cells in downstream,
1836random, and upstream numbering (left to right):
1838<img src="https://www.dealii.org/images/steps/developer/step-63-cell-order.png" alt="">
1840Let us start with the additive smoothers. The following table shows
1841the number of iterations necessary to obtain convergence from GMRES:
1843<table align="center" class="doxtable">
1847 <th colspan="1">@f$Q_1@f$</th>
1848 <th colspan="7">Smoother (smoothing steps)</th>
1854 <th colspan="3">Jacobi (6)</th>
1856 <th colspan="3">Block Jacobi (3)</th>
1862 <th colspan="3">Renumbering Strategy</th>
1864 <th colspan="3">Renumbering Strategy</th>
1964We see that renumbering the
1965DoFs/cells has no effect on convergence speed. This is because these
1966smoothers compute operations on each DoF (point-smoother) or cell
1967(block-smoother) independently and add up the results. Since we can
1968define these smoothers as an application of a sum of matrices, and
1969matrix addition is commutative, the order at which we sum the
1970different components will not affect the end result.
1972On the other hand, the situation is different for multiplicative smoothers:
1974<table align="center" class="doxtable">
1978 <th colspan="1">@f$Q_1@f$</th>
1979 <th colspan="7">Smoother (smoothing steps)</th>
1985 <th colspan="3">SOR (3)</th>
1987 <th colspan="3">Block SOR (1)</th>
1993 <th colspan="3">Renumbering Strategy</th>
1995 <th colspan="3">Renumbering Strategy</th>
2095Here, we can speed up
2096convergence by renumbering the DoFs/cells in the advection direction,
2097and similarly, we can slow down convergence if we do the renumbering
2098in the opposite direction. This is because advection-dominated
2099problems have a directional flow of information (in the advection
2100direction) which, given the right renumbering of DoFs/cells,
2101multiplicative methods are able to capture.
2103This feature of multiplicative methods is, however, dependent on the
2104value of @f$\varepsilon@f$. As we increase @f$\varepsilon@f$ and the problem
2105becomes more diffusion-dominated, we have a more uniform propagation
2106of information over the mesh and there is a diminished advantage for
2107renumbering in the advection direction. On the opposite end, in the
2108extreme case of @f$\varepsilon=0@f$ (advection-only), we have a 1st-order
2109PDE and multiplicative methods with the right renumbering become
2110effective solvers: A correct downstream numbering may lead to methods
2111that require only a single iteration because information can be
2112propagated from the inflow boundary downstream, with no information
2113transport in the opposite direction. (Note, however, that in the case
2114of @f$\varepsilon=0@f$, special care must be taken for the boundary
2115conditions in this case).
2118<a name="step_63-Pointvsblocksmoothers"></a><h4> %Point vs. block smoothers </h4>
2121We will limit the results to runs using the downstream
2122renumbering. Here is a cross comparison of all four smoothers for both
2123@f$Q_1@f$ and @f$Q_3@f$ elements:
2125<table align="center" class="doxtable">
2129 <th colspan="1">@f$Q_1@f$</th>
2130 <th colspan="4">Smoother (smoothing steps)</th>
2132 <th colspan="1">@f$Q_3@f$</th>
2133 <th colspan="4">Smoother (smoothing steps)</th>
2136 <th colspan="1">Cells</th>
2138 <th colspan="1">DoFs</th>
2139 <th colspan="1">Jacobi (6)</th>
2140 <th colspan="1">Block Jacobi (3)</th>
2141 <th colspan="1">SOR (3)</th>
2142 <th colspan="1">Block SOR (1)</th>
2144 <th colspan="1">DoFs</th>
2145 <th colspan="1">Jacobi (6)</th>
2146 <th colspan="1">Block Jacobi (3)</th>
2147 <th colspan="1">SOR (3)</th>
2148 <th colspan="1">Block SOR (1)</th>
2247We see that for @f$Q_1@f$, both multiplicative smoothers require a smaller
2248combination of smoothing steps and iteration counts than either
2249additive smoother. However, when we increase the degree to a @f$Q_3@f$
2250element, there is a clear advantage for the block smoothers in terms
2251of the number of smoothing steps and iterations required to
2252solve. Specifically, the block SOR smoother gives constant iteration
2253counts over the degree, and the block Jacobi smoother only sees about
2254a 38% increase in iterations compared to 75% and 183% for Jacobi and
2257<a name="step_63-Cost"></a><h3> Cost </h3>
2260Iteration counts do not tell the full story in the optimality of a one
2261smoother over another. Obviously we must examine the cost of an
2262iteration. Block smoothers here are at a disadvantage as they are
2263having to construct and invert a cell matrix for each cell. Here is a
2264comparison of solve times for a @f$Q_3@f$ element with 74,496 DoFs:
2266<table align="center" class="doxtable">
2268 <th colspan="1">@f$Q_3@f$</th>
2269 <th colspan="4">Smoother (smoothing steps)</th>
2272 <th colspan="1">DoFs</th>
2273 <th colspan="1">Jacobi (6)</th>
2274 <th colspan="1">Block Jacobi (3)</th>
2275 <th colspan="1">SOR (3)</th>
2276 <th colspan="1">Block SOR (1)</th>
2287The smoother that requires the most iterations (Jacobi) actually takes
2288the shortest time (roughly 2/3 the time of the next fastest
2289method). This is because all that is required to apply a Jacobi
2290smoothing step is multiplication by a diagonal matrix which is very
2291cheap. On the other hand, while SOR requires over 3x more iterations
2292(each with 3x more smoothing steps) than block SOR, the times are
2293roughly equivalent, implying that a smoothing step of block SOR is
2294roughly 9x slower than a smoothing step of SOR. Lastly, block Jacobi
2295is almost 6x more expensive than block SOR, which intuitively makes
2296sense from the fact that 1 step of each method has the same cost
2297(inverting the cell matrices and either adding or multiply them
2298together), and block Jacobi has 3 times the number of smoothing steps per
2299iteration with 2 times the iterations.
2302<a name="step_63-Additionalpoints"></a><h3> Additional points </h3>
2305There are a few more important points to mention:
2308<li> For a mesh distributed in parallel, multiplicative methods cannot
2309be executed over the entire domain. This is because they operate one
2310cell at a time, and downstream cells can only be handled once upstream
2311cells have already been done. This is fine on a single processor: The
2312processor just goes through the list of cells one after the
2313other. However, in parallel, it would imply that some processors are
2314idle because upstream processors have not finished doing the work on
2315cells upstream from the ones owned by the current processor. Once the
2316upstream processors are done, the downstream ones can start, but by
2317that time the upstream processors have no work left. In other words,
2318most of the time during these smoother steps, most processors are in
2319fact idle. This is not how one obtains good parallel scalability!
2321One can use a hybrid method where
2322a multiplicative smoother is applied on each subdomain, but as you
2323increase the number of subdomains, the method approaches the behavior
2324of an additive method. This is a major disadvantage to these methods.
2327<li> Current research into block smoothers suggest that soon we will be
2328able to compute the inverse of the cell matrices much cheaper than
2329what is currently being done inside deal.II. This research is based on
2330the fast diagonalization method (dating back to the 1960s) and has
2331been used in the spectral community for around 20 years (see, e.g., <a
2332href="https://doi.org/10.1007/s10915-004-4787-3"> Hybrid
2333Multigrid/Schwarz Algorithms for the Spectral Element Method by Lottes
2334and Fischer</a>). There are currently efforts to generalize these
2335methods to DG and make them more robust. Also, it seems that one
2336should be able to take advantage of matrix-free implementations and
2337the fact that, in the interior of the domain, cell matrices tend to
2338look very similar, allowing fewer matrix inverse computations.
2342Combining 1. and 2. gives a good reason for expecting that a method
2343like block Jacobi could become very powerful in the future, even
2344though currently for these examples it is quite slow.
2347<a name="step_63-Possibilitiesforextensions"></a><h3> Possibilities for extensions </h3>
2350<a name="step_63-ConstantiterationsforQsub5sub"></a><h4> Constant iterations for Q<sub>5</sub> </h4>
2353Change the number of smoothing steps and the smoother relaxation
2354parameter (set in <code>Smoother::AdditionalData()</code> inside
2355<code>create_smoother()</code>, only necessary for point smoothers) so
2356that we maintain a constant number of iterations for a @f$Q_5@f$ element.
2358<a name="step_63-Effectivenessofrenumberingforchangingepsilon"></a><h4> Effectiveness of renumbering for changing epsilon </h4>
2361Increase/decrease the parameter "Epsilon" in the `.prm` files of the
2362multiplicative methods and observe for which values renumbering no
2363longer influences convergence speed.
2365<a name="step_63-Meshadaptivity"></a><h4> Mesh adaptivity </h4>
2368The code is set up to work correctly with an adaptively refined mesh (the
2369interface matrices are created and set). Devise a suitable refinement
2370criterium or try the KellyErrorEstimator class.
2373<a name="step_63-PlainProg"></a>
2374<h1> The plain program</h1>
2375@include "step-63.cc"
virtual void value_list(const std::vector< Point< dim > > &points, std::vector< RangeNumberType > &values, const unsigned int component=0) const
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
void declare_entry(const std::string &entry, const std::string &default_value, const Patterns::PatternBase &pattern=Patterns::Anything(), const std::string &documentation="", const bool has_to_be_set=false)
#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()
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.
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)
static constexpr double PI
::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 > &)