1346 *
template <
int dim>
1350 *
for (
unsigned int c = 0; c <
vp.size(); ++c)
1360 * <a name=
"step_22-Linearsolversandpreconditioners"></a>
1371 * <a name=
"step_22-ThecodeInverseMatrixcodeclasstemplate"></a>
1390 *
const PreconditionerType &preconditioner);
1400 *
template <
class MatrixType,
class PreconditionerType>
1402 *
const MatrixType &m,
1403 *
const PreconditionerType &preconditioner)
1405 *
, preconditioner(&preconditioner)
1422 * tolerance, either.
1425 * template <class MatrixType, class PreconditionerType>
1426 * void InverseMatrix<MatrixType, PreconditionerType>::vmult(
1427 * Vector<double> &dst,
1428 * const Vector<double> &src) const
1430 * SolverControl solver_control(src.size(), 1e-6 * src.l2_norm());
1431 * SolverCG<Vector<double>> cg(solver_control);
1435 * cg.solve(*matrix, dst, src, *preconditioner);
1442 * <a name="step_22-ThecodeSchurComplementcodeclasstemplate"></a>
1443 * <h4>The <code>SchurComplement</code> class template</h4>
1447 * This class implements the Schur complement discussed in the introduction.
1448 * It is in analogy to @ref step_20 "step-20". Though, we now call it with a template
1449 * parameter <code>PreconditionerType</code> in order to access that when
1450 * specifying the respective type of the inverse matrix class. As a
1451 * consequence of the definition above, the declaration
1452 * <code>InverseMatrix</code> now contains the second template parameter for
1453 * a preconditioner class as above, which affects the
1454 * <code>ObserverPointer</code> object <code>m_inverse</code> as well.
1457 * template <class PreconditionerType>
1458 * class SchurComplement : public EnableObserverPointer
1462 * const BlockSparseMatrix<double> &system_matrix,
1463 * const InverseMatrix<SparseMatrix<double>, PreconditionerType> &A_inverse);
1465 * void vmult(Vector<double> &dst, const Vector<double> &src) const;
1468 * const ObserverPointer<const BlockSparseMatrix<double>> system_matrix;
1469 * const ObserverPointer<
1470 * const InverseMatrix<SparseMatrix<double>, PreconditionerType>>
1473 * mutable Vector<double> tmp1, tmp2;
1478 * template <class PreconditionerType>
1479 * SchurComplement<PreconditionerType>::SchurComplement(
1480 * const BlockSparseMatrix<double> &system_matrix,
1481 * const InverseMatrix<SparseMatrix<double>, PreconditionerType> &A_inverse)
1482 * : system_matrix(&system_matrix)
1483 * , A_inverse(&A_inverse)
1484 * , tmp1(system_matrix.block(0, 0).m())
1485 * , tmp2(system_matrix.block(0, 0).m())
1489 * template <class PreconditionerType>
1491 * SchurComplement<PreconditionerType>::vmult(Vector<double> &dst,
1492 * const Vector<double> &src) const
1494 * system_matrix->block(0, 1).vmult(tmp1, src);
1495 * A_inverse->vmult(tmp2, tmp1);
1496 * system_matrix->block(1, 0).vmult(dst, tmp2);
1503 * <a name="step_22-StokesProblemclassimplementation"></a>
1504 * <h3>StokesProblem class implementation</h3>
1509 * <a name="step_22-StokesProblemStokesProblem"></a>
1510 * <h4>StokesProblem::StokesProblem</h4>
1514 * The constructor of this class looks very similar to the one of
1515 * @ref step_20 "step-20". The constructor initializes the variables for the polynomial
1516 * degree, triangulation, finite element system and the dof handler. The
1517 * underlying polynomial functions are of order <code>degree+1</code> for
1518 * the vector-valued velocity components and of order <code>degree</code>
1519 * for the pressure. This gives the LBB-stable element pair
1520 * @f$Q_{degree+1}^d\times Q_{degree}@f$, often referred to as the Taylor-Hood
1525 * Note that we initialize the triangulation with a MeshSmoothing argument,
1526 * which ensures that the refinement of cells is done in a way that the
1527 * approximation of the PDE solution remains well-behaved (problems arise if
1528 * grids are too unstructured), see the documentation of
1529 * <code>Triangulation::MeshSmoothing</code> for details.
1532 * template <int dim>
1533 * StokesProblem<dim>::StokesProblem(const unsigned int degree)
1535 * , triangulation(Triangulation<dim>::maximum_smoothing)
1536 * , fe(FE_Q<dim>(degree + 1) ^ dim, FE_Q<dim>(degree))
1537 * , dof_handler(triangulation)
1544 * <a name="step_22-StokesProblemsetup_dofs"></a>
1545 * <h4>StokesProblem::setup_dofs</h4>
1549 * Given a mesh, this function associates the degrees of freedom with it and
1550 * creates the corresponding matrices and vectors. At the beginning it also
1551 * releases the pointer to the preconditioner object (if the shared pointer
1552 * pointed at anything at all at this point) since it will definitely not be
1553 * needed any more after this point and will have to be re-computed after
1554 * assembling the matrix, and unties the sparse matrices from their sparsity
1559 * We then proceed with distributing degrees of freedom and renumbering
1560 * them: In order to make the ILU preconditioner (in 3d) work efficiently,
1561 * it is important to enumerate the degrees of freedom in such a way that it
1562 * reduces the bandwidth of the matrix, or maybe more importantly: in such a
1563 * way that the ILU is as close as possible to a real LU decomposition. On
1564 * the other hand, we need to preserve the block structure of velocity and
1565 * pressure already seen in @ref step_20 "step-20" and @ref step_21 "step-21". This is done in two
1566 * steps: First, all dofs are renumbered to improve the ILU and then we
1567 * renumber once again by components. Since
1568 * <code>DoFRenumbering::component_wise</code> does not touch the
1569 * renumbering within the individual blocks, the basic renumbering from the
1570 * first step remains. As for how the renumber degrees of freedom to improve
1571 * the ILU: deal.II has a number of algorithms that attempt to find
1572 * orderings to improve ILUs, or reduce the bandwidth of matrices, or
1573 * optimize some other aspect. The DoFRenumbering namespace shows a
1574 * comparison of the results we obtain with several of these algorithms
1575 * based on the testcase discussed here in this tutorial program. Here, we
1576 * will use the traditional Cuthill-McKee algorithm already used in some of
1577 * the previous tutorial programs. In the <a href="#improved-ilu">section
1596 *
template <
int dim>
1600 *
system_matrix.
clear();
1603 *
dof_handler.distribute_dofs(fe);
1636 *
constraints.
clear();
1647 *
constraints.close();
1664 *
std::cout <<
" Number of active cells: " <<
triangulation.n_active_cells()
1666 *
<<
" Number of degrees of freedom: " << dof_handler.n_dofs()
1667 *
<<
" (" <<
n_u <<
'+' <<
n_p <<
')' << std::endl;
1676 *
In 3
d,
the function DoFTools::max_couplings_between_dofs
yields a
1681 * of most systems already for moderately-sized 3d problems, see also the
1682 * discussion in @ref step_18 "step-18". Instead, we first build temporary objects that use
1685 * BlockSparseMatrix objects; in a second step we then copy these objects
1686 * into objects of type BlockSparsityPattern. This is entirely analogous to
1687 * what we already did in @ref step_11 "step-11" and @ref step_18 "step-18". In particular, we make use of
1688 * the fact that we will never write into the @f$(1,1)@f$ block of the system
1689 * matrix and that this is the only block to be filled for the
1690 * preconditioner matrix.
1694 * All this is done inside new scopes, which means that the memory of
1695 * <code>dsp</code> will be released once the information has been copied to
1696 * <code>sparsity_pattern</code>.
1700 * BlockDynamicSparsityPattern dsp(dofs_per_block, dofs_per_block);
1702 * Table<2, DoFTools::Coupling> coupling(dim + 1, dim + 1);
1703 * for (unsigned int c = 0; c < dim + 1; ++c)
1704 * for (unsigned int d = 0; d < dim + 1; ++d)
1705 * if (!((c == dim) && (d == dim)))
1706 * coupling[c][d] = DoFTools::always;
1708 * coupling[c][d] = DoFTools::none;
1710 * DoFTools::make_sparsity_pattern(
1711 * dof_handler, coupling, dsp, constraints, false);
1713 * sparsity_pattern.copy_from(dsp);
1717 * BlockDynamicSparsityPattern preconditioner_dsp(dofs_per_block,
1720 * Table<2, DoFTools::Coupling> preconditioner_coupling(dim + 1, dim + 1);
1721 * for (unsigned int c = 0; c < dim + 1; ++c)
1722 * for (unsigned int d = 0; d < dim + 1; ++d)
1723 * if (((c == dim) && (d == dim)))
1724 * preconditioner_coupling[c][d] = DoFTools::always;
1726 * preconditioner_coupling[c][d] = DoFTools::none;
1728 * DoFTools::make_sparsity_pattern(dof_handler,
1729 * preconditioner_coupling,
1730 * preconditioner_dsp,
1734 * preconditioner_sparsity_pattern.copy_from(preconditioner_dsp);
1739 * Finally, the system matrix, the preconsitioner matrix, the solution and
1740 * the right hand side vector are created from the block structure similar
1741 * to the approach in @ref step_20 "step-20":
1744 * system_matrix.reinit(sparsity_pattern);
1745 * preconditioner_matrix.reinit(preconditioner_sparsity_pattern);
1747 * solution.reinit(dofs_per_block);
1748 * system_rhs.reinit(dofs_per_block);
1755 * <a name="step_22-StokesProblemassemble_system"></a>
1756 * <h4>StokesProblem::assemble_system</h4>
1760 * The assembly process follows the discussion in @ref step_20 "step-20" and in the
1761 * introduction. We use the well-known abbreviations for the data structures
1762 * that hold the local matrices, right hand side, and global numbering of the
1763 * degrees of freedom for the present cell.
1766 * template <int dim>
1767 * void StokesProblem<dim>::assemble_system()
1769 * system_matrix = 0;
1771 * preconditioner_matrix = 0;
1773 * const QGauss<dim> quadrature_formula(degree + 2);
1775 * FEValues<dim> fe_values(fe,
1776 * quadrature_formula,
1777 * update_values | update_quadrature_points |
1778 * update_JxW_values | update_gradients);
1780 * const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
1782 * const unsigned int n_q_points = quadrature_formula.size();
1784 * FullMatrix<double> local_matrix(dofs_per_cell, dofs_per_cell);
1785 * FullMatrix<double> local_preconditioner_matrix(dofs_per_cell,
1787 * Vector<double> local_rhs(dofs_per_cell);
1789 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1791 * const RightHandSide<dim> right_hand_side;
1792 * std::vector<Tensor<1, dim>> rhs_values(n_q_points, Tensor<1, dim>());
1796 * Next, we need two objects that work as extractors for the FEValues
1797 * object. Their use is explained in detail in the report on @ref
1801 * const FEValuesExtractors::Vector velocities(0);
1802 * const FEValuesExtractors::Scalar pressure(dim);
1806 * As an extension over @ref step_20 "step-20" and @ref step_21 "step-21", we include a few optimizations
1807 * that make assembly much faster for this particular problem. The
1808 * improvements are based on the observation that we do a few calculations
1809 * too many times when we do as in @ref step_20 "step-20": The symmetric gradient actually
1810 * has <code>dofs_per_cell</code> different values per quadrature point, but
1811 * we extract it <code>dofs_per_cell*dofs_per_cell</code> times from the
1812 * FEValues object - for both the loop over <code>i</code> and the inner
1813 * loop over <code>j</code>. In 3d, that means evaluating it @f$89^2=7921@f$
1814 * instead of @f$89@f$ times, a not insignificant difference.
1832 *
std::vector<SymmetricTensor<2, dim>>
symgrad_phi_u(dofs_per_cell);
1833 *
std::vector<double>
div_phi_u(dofs_per_cell);
1834 *
std::vector<Tensor<1, dim>>
phi_u(dofs_per_cell);
1835 *
std::vector<double>
phi_p(dofs_per_cell);
1837 *
for (
const auto &cell : dof_handler.active_cell_iterators())
1839 *
fe_values.
reinit(cell);
1848 *
for (
unsigned int q = 0;
q < n_q_points; ++
q)
1850 *
for (
unsigned int k = 0;
k < dofs_per_cell; ++
k)
1891 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1893 *
for (
unsigned int j = 0;
j <= i; ++
j)
1899 *
* fe_values.JxW(
q);
1903 *
* fe_values.JxW(
q);
1925 *
* fe_values.JxW(
q);
1939 *
of the local matrices.
1942 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1943 *
for (
unsigned int j = i + 1;
j < dofs_per_cell; ++
j)
1950 *
cell->get_dof_indices(local_dof_indices);
1953 *
local_dof_indices,
1957 *
local_dof_indices,
1963 *
Before we're going to solve this linear system, we generate a
1964 * preconditioner for the velocity-velocity matrix, i.e.,
1965 * <code>block(0,0)</code> in the system matrix. As mentioned above, this
1966 * depends on the spatial dimension. Since the two classes described by
1967 * the <code>InnerPreconditioner::type</code> alias have the same
1968 * interface, we do not have to do anything different whether we want to
1969 * use a sparse direct solver or an ILU:
1972 * std::cout << " Computing preconditioner..." << std::endl << std::flush;
1974 * A_preconditioner =
1975 * std::make_shared<typename InnerPreconditioner<dim>::type>();
1976 * A_preconditioner->initialize(
1977 * system_matrix.block(0, 0),
1978 * typename InnerPreconditioner<dim>::type::AdditionalData());
1986 * <a name="step_22-StokesProblemsolve"></a>
1987 * <h4>StokesProblem::solve</h4>
1991 * After the discussion in the introduction and the definition of the
1992 * respective classes above, the implementation of the <code>solve</code>
1993 * function is rather straight-forward and done in a similar way as in
1994 * @ref step_20 "step-20". To start with, we need an object of the
1995 * <code>InverseMatrix</code> class that represents the inverse of the
1996 * matrix A. As described in the introduction, the inverse is generated with
1997 * the help of an inner preconditioner of type
1998 * <code>InnerPreconditioner::type</code>.
2001 * template <int dim>
2002 * void StokesProblem<dim>::solve()
2004 * const InverseMatrix<SparseMatrix<double>,
2005 * typename InnerPreconditioner<dim>::type>
2006 * A_inverse(system_matrix.block(0, 0), *A_preconditioner);
2007 * Vector<double> tmp(solution.block(0).size());
2011 * This is as in @ref step_20 "step-20". We generate the right hand side @f$B A^{-1} F - G@f$
2012 * for the Schur complement and an object that represents the respective
2013 * linear operation @f$B A^{-1} B^T@f$, now with a template parameter
2014 * indicating the preconditioner - in accordance with the definition of
2019 * Vector<double> schur_rhs(solution.block(1).size());
2020 * A_inverse.vmult(tmp, system_rhs.block(0));
2021 * system_matrix.block(1, 0).vmult(schur_rhs, tmp);
2022 * schur_rhs -= system_rhs.block(1);
2024 * SchurComplement<typename InnerPreconditioner<dim>::type> schur_complement(
2025 * system_matrix, A_inverse);
2029 * The usual control structures for the solver call are created...
2032 * SolverControl solver_control(solution.block(1).size(),
2033 * 1e-6 * schur_rhs.l2_norm());
2034 * SolverCG<Vector<double>> cg(solver_control);
2038 * Now to the preconditioner to the Schur complement. As explained in
2039 * the introduction, the preconditioning is done by a @ref GlossMassMatrix "mass matrix" in the
2040 * pressure variable.
2044 * Actually, the solver needs to have the preconditioner in the form
2045 * @f$P^{-1}@f$, so we need to create an inverse operation. Once again, we
2046 * use an object of the class <code>InverseMatrix</code>, which
2047 * implements the <code>vmult</code> operation that is needed by the
2048 * solver. In this case, we have to invert the pressure mass matrix. As
2049 * it already turned out in earlier tutorial programs, the inversion of
2050 * a mass matrix is a rather cheap and straight-forward operation
2051 * (compared to, e.g., a Laplace matrix). The CG method with ILU
2052 * preconditioning converges in 5-10 steps, independently on the mesh
2053 * size. This is precisely what we do here: We choose another ILU
2054 * preconditioner and take it along to the InverseMatrix object via the
2055 * corresponding template parameter. A CG solver is then called within
2056 * the vmult operation of the inverse matrix.
2060 * An alternative that is cheaper to build, but needs more iterations
2061 * afterwards, would be to choose a SSOR preconditioner with factor
2062 * 1.2. It needs about twice the number of iterations, but the costs for
2063 * its generation are almost negligible.
2066 * SparseILU<double> preconditioner;
2067 * preconditioner.initialize(preconditioner_matrix.block(1, 1),
2068 * SparseILU<double>::AdditionalData());
2070 * InverseMatrix<SparseMatrix<double>, SparseILU<double>> m_inverse(
2071 * preconditioner_matrix.block(1, 1), preconditioner);
2075 * With the Schur complement and an efficient preconditioner at hand, we
2076 * can solve the respective equation for the pressure (i.e. block 0 in
2077 * the solution vector) in the usual way:
2080 * cg.solve(schur_complement, solution.block(1), schur_rhs, m_inverse);
2084 * After this first solution step, the hanging node constraints have to
2085 * be distributed to the solution in order to achieve a consistent
2089 * constraints.distribute(solution);
2091 * std::cout << " " << solver_control.last_step()
2092 * << " outer CG Schur complement iterations for pressure"
2098 * As in @ref step_20 "step-20", we finally need to solve for the velocity equation where
2099 * we plug in the solution to the pressure equation. This involves only
2100 * objects we already know - so we simply multiply @f$p@f$ by @f$B^T@f$, subtract
2101 * the right hand side and multiply by the inverse of @f$A@f$. At the end, we
2102 * need to distribute the constraints from hanging nodes in order to
2103 * obtain a consistent flow field:
2107 * system_matrix.block(0, 1).vmult(tmp, solution.block(1));
2109 * tmp += system_rhs.block(0);
2111 * A_inverse.vmult(solution.block(0), tmp);
2113 * constraints.distribute(solution);
2121 * <a name="step_22-StokesProblemoutput_results"></a>
2122 * <h4>StokesProblem::output_results</h4>
2126 * The next function generates graphical output. In this example, we are
2127 * going to use the VTK file format. We attach names to the individual
2128 * variables in the problem: <code>velocity</code> to the <code>dim</code>
2129 * components of velocity and <code>pressure</code> to the pressure.
2133 * Not all visualization programs have the ability to group individual
2134 * vector components into a vector to provide vector plots; in particular,
2135 * this holds for some VTK-based visualization programs. In this case, the
2136 * logical grouping of components into vectors should already be described
2137 * in the file containing the data. In other words, what we need to do is
2138 * provide our output writers with a way to know which of the components of
2139 * the finite element logically form a vector (with @f$d@f$ components in @f$d@f$
2140 * space dimensions) rather than letting them assume that we simply have a
2141 * bunch of scalar fields. This is achieved using the members of the
2142 * <code>DataComponentInterpretation</code> namespace: as with the filename,
2143 * we create a vector in which the first <code>dim</code> components refer
2144 * to the velocities and are given the tag
2145 * DataComponentInterpretation::component_is_part_of_vector; we
2146 * finally push one tag
2147 * DataComponentInterpretation::component_is_scalar to describe
2148 * the grouping of the pressure variable.
2152 * The rest of the function is then the same as in @ref step_20 "step-20".
2155 * template <int dim>
2157 * StokesProblem<dim>::output_results(const unsigned int refinement_cycle) const
2159 * std::vector<std::string> solution_names(dim, "velocity");
2160 * solution_names.emplace_back("pressure");
2162 * std::vector<DataComponentInterpretation::DataComponentInterpretation>
2163 * data_component_interpretation(
2164 * dim, DataComponentInterpretation::component_is_part_of_vector);
2165 * data_component_interpretation.push_back(
2166 * DataComponentInterpretation::component_is_scalar);
2168 * DataOut<dim> data_out;
2169 * data_out.attach_dof_handler(dof_handler);
2170 * data_out.add_data_vector(solution,
2172 * DataOut<dim>::type_dof_data,
2173 * data_component_interpretation);
2174 * data_out.build_patches();
2176 * std::ofstream output(
2177 * "solution-" + Utilities::int_to_string(refinement_cycle, 2) + ".vtk");
2178 * data_out.write_vtk(output);
2185 * <a name="step_22-StokesProblemrefine_mesh"></a>
2186 * <h4>StokesProblem::refine_mesh</h4>
2190 * This is the last interesting function of the <code>StokesProblem</code>
2191 * class. As indicated by its name, it takes the solution to the problem
2192 * and refines the mesh where this is needed. The procedure is the same as
2193 * in the respective step in @ref step_6 "step-6", with the exception that we base the
2194 * refinement only on the change in pressure, i.e., we call the Kelly error
2195 * estimator with a mask object of type ComponentMask that selects the
2196 * single scalar component for the pressure that we are interested in (we
2197 * get such a mask from the finite element class by specifying the component
2198 * we want). Additionally, we do not coarsen the grid again:
2201 * template <int dim>
2202 * void StokesProblem<dim>::refine_mesh()
2204 * Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
2206 * const FEValuesExtractors::Scalar pressure(dim);
2207 * KellyErrorEstimator<dim>::estimate(
2209 * QGauss<dim - 1>(degree + 1),
2210 * std::map<types::boundary_id, const Function<dim> *>(),
2212 * estimated_error_per_cell,
2213 * fe.component_mask(pressure));
2215 * GridRefinement::refine_and_coarsen_fixed_number(triangulation,
2216 * estimated_error_per_cell,
2219 * triangulation.execute_coarsening_and_refinement();
2226 * <a name="step_22-StokesProblemrun"></a>
2227 * <h4>StokesProblem::run</h4>
2231 * The last step in the Stokes class is, as usual, the function that
2232 * generates the initial grid and calls the other functions in the
2237 * We start off with a rectangle of size @f$4 \times 1@f$ (in 2d) or @f$4 \times 1
2238 * \times 1@f$ (in 3d), placed in @f$R^2/R^3@f$ as @f$(-2,2)\times(-1,0)@f$ or
2239 * @f$(-2,2)\times(0,1)\times(-1,0)@f$, respectively. It is natural to start
2240 * with equal mesh size in each direction, so we subdivide the initial
2241 * rectangle four times in the first coordinate direction. To limit the
2242 * scope of the variables involved in the creation of the mesh to the range
2243 * where we actually need them, we put the entire block between a pair of
2247 * template <int dim>
2248 * void StokesProblem<dim>::run()
2251 * std::vector<unsigned int> subdivisions(dim, 1);
2252 * subdivisions[0] = 4;
2254 * const Point<dim> bottom_left = (dim == 2 ?
2255 * Point<dim>(-2, -1) : // 2d case
2256 * Point<dim>(-2, 0, -1)); // 3d case
2258 * const Point<dim> top_right = (dim == 2 ?
2259 * Point<dim>(2, 0) : // 2d case
2260 * Point<dim>(2, 1, 0)); // 3d case
2262 * GridGenerator::subdivided_hyper_rectangle(triangulation,
2270 * A boundary indicator of 1 is set to all boundaries that are subject to
2271 * Dirichlet boundary conditions, i.e. to faces that are located at 0 in
2272 * the last coordinate direction. See the example description above for
2276 * for (const auto &cell : triangulation.active_cell_iterators())
2277 * for (const auto &face : cell->face_iterators())
2278 * if (face->center()[dim - 1] == 0)
2279 * face->set_all_boundary_ids(1);
2284 * We then apply an initial refinement before solving for the first
2285 * time. In 3d, there are going to be more degrees of freedom, so we
2286 * refine less there:
2289 * triangulation.refine_global(4 - dim);
2293 * As first seen in @ref step_6 "step-6", we cycle over the different refinement levels
2294 * and refine (except for the first cycle), setup the degrees of freedom
2295 * and matrices, assemble, solve and create output:
2298 * for (unsigned int refinement_cycle = 0; refinement_cycle < 6;
2299 * ++refinement_cycle)
2301 * std::cout << "Refinement cycle " << refinement_cycle << std::endl;
2303 * if (refinement_cycle > 0)
2308 * std::cout << " Assembling..." << std::endl << std::flush;
2309 * assemble_system();
2311 * std::cout << " Solving..." << std::flush;
2314 * output_results(refinement_cycle);
2316 * std::cout << std::endl;
2319 * } // namespace Step22
2325 * <a name="step_22-Thecodemaincodefunction"></a>
2326 * <h3>The <code>main</code> function</h3>
2330 * The main function is the same as in @ref step_20 "step-20". We pass the element degree as
2331 * a parameter and choose the space dimension at the well-known template slot.
2338 * using namespace Step22;
2340 * StokesProblem<2> flow_problem(1);
2341 * flow_problem.run();
2343 * catch (std::exception &exc)
2345 * std::cerr << std::endl
2347 * << "----------------------------------------------------"
2349 * std::cerr << "Exception on processing: " << std::endl
2350 * << exc.what() << std::endl
2351 * << "Aborting!" << std::endl
2352 * << "----------------------------------------------------"
2359 * std::cerr << std::endl
2361 * << "----------------------------------------------------"
2363 * std::cerr << "Unknown exception!" << std::endl
2364 * << "Aborting!" << std::endl
2365 * << "----------------------------------------------------"
2373<a name="step_22-Results"></a><h1>Results</h1>
2376<a name="step_22-Outputoftheprogramandgraphicalvisualization"></a><h3>Output of the program and graphical visualization</h3>
2379<a name="step_22-2Dcalculations"></a><h4>2D calculations</h4>
2382Running the program with the space dimension set to 2 in the <code>main</code>
2383function yields the following output (in "release mode",
2384See also <a href="https://www.math.colostate.edu/~bangerth/videos.676.18.html">video lecture 18</a>.):
2386examples/step-22> make run
2388 Number of active cells: 64
2389 Number of degrees of freedom: 679 (594+85)
2391 Computing preconditioner...
2392 Solving... 11 outer CG Schur complement iterations for pressure
2395 Number of active cells: 160
2396 Number of degrees of freedom: 1683 (1482+201)
2398 Computing preconditioner...
2399 Solving... 11 outer CG Schur complement iterations for pressure
2402 Number of active cells: 376
2403 Number of degrees of freedom: 3813 (3370+443)
2405 Computing preconditioner...
2406 Solving... 11 outer CG Schur complement iterations for pressure
2409 Number of active cells: 880
2410 Number of degrees of freedom: 8723 (7722+1001)
2412 Computing preconditioner...
2413 Solving... 11 outer CG Schur complement iterations for pressure
2416 Number of active cells: 2008
2417 Number of degrees of freedom: 19383 (17186+2197)
2419 Computing preconditioner...
2420 Solving... 11 outer CG Schur complement iterations for pressure
2423 Number of active cells: 4288
2424 Number of degrees of freedom: 40855 (36250+4605)
2426 Computing preconditioner...
2427 Solving... 11 outer CG Schur complement iterations for pressure
2430The entire computation above takes about 2 seconds on a reasonably
2431quick (for 2015 standards) machine.
2433What we see immediately from this is that the number of (outer)
2434iterations does not increase as we refine the mesh. This confirms the
2435statement in the introduction that preconditioning the Schur
2436complement with the mass matrix indeed yields a matrix spectrally
2437equivalent to the identity matrix (i.e. with eigenvalues bounded above
2438and below independently of the mesh size or the relative sizes of
2439cells). In other words, the mass matrix and the Schur complement are
2440spectrally equivalent.
2442In the images below, we show the grids for the first six refinement
2443steps in the program. Observe how the grid is refined in regions
2444where the solution rapidly changes: On the upper boundary, we have
2445Dirichlet boundary conditions that are -1 in the left half of the line
2446and 1 in the right one, so there is an abrupt change at @f$x=0@f$. Likewise,
2447there are changes from Dirichlet to Neumann data in the two upper
2448corners, so there is need for refinement there as well:
2450<table width="60%" align="center">
2453 <img src="https://www.dealii.org/images/steps/developer/step-22.2d.mesh-0.png" alt="">
2456 <img src="https://www.dealii.org/images/steps/developer/step-22.2d.mesh-1.png" alt="">
2461 <img src="https://www.dealii.org/images/steps/developer/step-22.2d.mesh-2.png" alt="">
2464 <img src="https://www.dealii.org/images/steps/developer/step-22.2d.mesh-3.png" alt="">
2469 <img src="https://www.dealii.org/images/steps/developer/step-22.2d.mesh-4.png" alt="">
2472 <img src="https://www.dealii.org/images/steps/developer/step-22.2d.mesh-5.png" alt="">
2477Finally, following is a plot of the flow field. It shows fluid
2478transported along with the moving upper boundary and being replaced by
2479material coming from below:
2481<img src="https://www.dealii.org/images/steps/developer/step-22.2d.solution.png" alt="">
2483This plot uses the capability of VTK-based visualization programs (in
2484this case of VisIt) to show vector data; this is the result of us
2485declaring the velocity components of the finite element in use to be a
2486set of vector components, rather than independent scalar components in
2487the <code>StokesProblem@<dim@>::%output_results</code> function of this
2492<a name="step_22-3Dcalculations"></a><h4>3D calculations</h4>
2495In 3d, the screen output of the program looks like this:
2499 Number of active cells: 32
2500 Number of degrees of freedom: 1356 (1275+81)
2502 Computing preconditioner...
2503 Solving... 13 outer CG Schur complement iterations for pressure.
2506 Number of active cells: 144
2507 Number of degrees of freedom: 5088 (4827+261)
2509 Computing preconditioner...
2510 Solving... 14 outer CG Schur complement iterations for pressure.
2513 Number of active cells: 704
2514 Number of degrees of freedom: 22406 (21351+1055)
2516 Computing preconditioner...
2517 Solving... 14 outer CG Schur complement iterations for pressure.
2520 Number of active cells: 3168
2521 Number of degrees of freedom: 93176 (89043+4133)
2523 Computing preconditioner...
2524 Solving... 15 outer CG Schur complement iterations for pressure.
2527 Number of active cells: 11456
2528 Number of degrees of freedom: 327808 (313659+14149)
2530 Computing preconditioner...
2531 Solving... 15 outer CG Schur complement iterations for pressure.
2534 Number of active cells: 45056
2535 Number of degrees of freedom: 1254464 (1201371+53093)
2537 Computing preconditioner...
2538 Solving... 14 outer CG Schur complement iterations for pressure.
2541Again, we see that the number of outer iterations does not increase as
2542we refine the mesh. Nevertheless, the compute time increases
2543significantly: for each of the iterations above separately, it takes about
25440.14 seconds, 0.63 seconds, 4.8 seconds, 35 seconds, 2 minutes and 33 seconds,
2545and 13 minutes and 12 seconds. This overall superlinear (in the number of
2546unknowns) increase in runtime is due to the fact that our inner solver is not
2547@f${\cal O}(N)@f$: a simple experiment shows that as we keep refining the mesh, the
2548average number of ILU-preconditioned CG iterations to invert the
2549velocity-velocity block @f$A@f$ increases.
2551We will address the question of how possibly to improve our solver <a
2552href="#improved-solver">below</a>.
2554As for the graphical output, the grids generated during the solution
2557<table width="60%" align="center">
2560 <img src="https://www.dealii.org/images/steps/developer/step-22.3d.mesh-0.png" alt="">
2563 <img src="https://www.dealii.org/images/steps/developer/step-22.3d.mesh-1.png" alt="">
2568 <img src="https://www.dealii.org/images/steps/developer/step-22.3d.mesh-2.png" alt="">
2571 <img src="https://www.dealii.org/images/steps/developer/step-22.3d.mesh-3.png" alt="">
2576 <img src="https://www.dealii.org/images/steps/developer/step-22.3d.mesh-4.png" alt="">
2579 <img src="https://www.dealii.org/images/steps/developer/step-22.3d.mesh-5.png" alt="">
2584Again, they show essentially the location of singularities introduced
2585by boundary conditions. The vector field computed makes for an
2588<img src="https://www.dealii.org/images/steps/developer/step-22.3d.solution.png" alt="">
2590The isocontours shown here as well are those of the pressure
2591variable, showing the singularity at the point of discontinuous
2592velocity boundary conditions.
2596<a name="step_22-Sparsitypattern"></a><h3>Sparsity pattern</h3>
2599As explained during the generation of the sparsity pattern, it is
2600important to have the numbering of degrees of freedom in mind when
2601using preconditioners like incomplete LU decompositions. This is most
2602conveniently visualized using the distribution of nonzero elements in
2603the @ref GlossStiffnessMatrix "stiffness matrix".
2612<
img src=
"https://www.dealii.org/images/steps/developer/step-22.2d.sparsity-nor.png" alt=
"">
2618 std::ofstream out (
"sparsity_pattern.gpl");
2619 sparsity_pattern.print_gnuplot(out);
2633<
img src=
"https://www.dealii.org/images/steps/developer/step-22.2d.sparsity-ren.png" alt=
"">
2730<a name="step-22-block-
schur"></a>
2752 A^{-1} & 0 \\ S^{-1}
B A^{-1} & -S^{-1}
2761 \left(\begin{array}{
cc}
2762 A^{-1} & 0 \\ S^{-1}
B A^{-1} & -S^{-1}
2763 \
end{array}\right)\
cdot \left(\begin{array}{
cc}
2767 \left(\begin{array}{
cc}
2768 I &
A^{-1}
B^
T \\ 0 &
I
2791template <
class PreconditionerA,
class PreconditionerMp>
2812template <
class PreconditionerA,
class PreconditionerMp>
2822 tmp (S.block(1,1).m())
2827template <
class PreconditionerA,
class PreconditionerMp>
2837 system_matrix->block(1,0).residual(tmp, dst.block(0), src.block(1));
2861<a
href=
"http://en.wikipedia.org/wiki/GMRES#Comparison_with_other_solvers">
Wikipedia</a>'s
2862article on the GMRES method gives a comparative presentation.
2863A more comprehensive and well-founded comparison can be read e.g. in the book by
2864J.W. Demmel (Applied Numerical Linear Algebra, SIAM, 1997, section 6.6.6).
2866For our specific problem with the ILU preconditioner for @f$A@f$, we certainly need
2867to perform hundreds of iterations on the block system for large problem sizes
2875basis
to 30 vectors
by default.
2883(one iteration uses only results from one preceding step and
2884not all the steps as GMRES). Besides the fact the BiCGStab is more expensive per
2885step since two matrix-vector products are needed (compared to one for
2886CG or GMRES), there is one main reason which makes BiCGStab not appropriate for
2887this problem: The preconditioner applies the inverse of the pressure
2888mass matrix by using the InverseMatrix class. Since the application of the
2889inverse matrix to a vector is done only in approximative way (an exact inverse
2890is too expensive), this will also affect the solver. In the case of BiCGStab,
2891the Krylov vectors will not be orthogonal due to that perturbation. While
2892this is uncritical for a small number of steps (up to about 50), it ruins the
2893performance of the solver when these perturbations have grown to a significant
2894magnitude in the coarse of iterations.
2896We did some experiments with BiCGStab and found it to
2897be faster than GMRES up to refinement cycle 3 (in 3D), but it became very slow
2898for cycles 4 and 5 (even slower than the original Schur complement), so the
2899solver is useless in this situation. Choosing a sharper tolerance for the
2900inverse matrix class (<code>1e-10*src.l2_norm()</code> instead of
2901<code>1e-6*src.l2_norm()</code>) made BiCGStab perform well also for cycle 4,
2902but did not change the failure on the very large problems.
2904GMRES is of course also effected by the approximate inverses, but it is not as
2905sensitive to orthogonality and retains a relatively good performance also for
2906large sizes, see the results below.
2908With this said, we turn to the realization of the solver call with GMRES with
2909@f$k=100@f$ temporary vectors:
2912 const SparseMatrix<double> &pressure_mass_matrix
2913 = preconditioner_matrix.block(1,1);
2914 SparseILU<double> pmass_preconditioner;
2915 pmass_preconditioner.initialize (pressure_mass_matrix,
2916 SparseILU<double>::AdditionalData());
2918 InverseMatrix<SparseMatrix<double>,SparseILU<double> >
2919 m_inverse (pressure_mass_matrix, pmass_preconditioner);
2921 BlockSchurPreconditioner<typename InnerPreconditioner<dim>::type,
2923 preconditioner (system_matrix, m_inverse, *A_preconditioner);
2925 SolverControl solver_control (system_matrix.m(),
2926 1e-6*system_rhs.l2_norm());
2927 GrowingVectorMemory<BlockVector<double> > vector_memory;
2928 SolverGMRES<BlockVector<double> >::AdditionalData gmres_data;
2929 gmres_data.max_basis_size = 100;
2931 SolverGMRES<BlockVector<double> > gmres(solver_control, vector_memory,
2934 gmres.solve(system_matrix, solution, system_rhs,
2937 constraints.distribute (solution);
2940 << solver_control.last_step()
2941 << " block GMRES iterations";
2944Obviously, one needs to add the include file @ref SolverGMRES
2945"<lac/solver_gmres.h>" in order to make this run.
2946We call the solver with a BlockVector template in order to enable
2947GMRES to operate on block vectors and matrices.
2948Note also that we need to set the (1,1) block in the system
2949matrix to zero (we saved the pressure mass matrix there which is not part of the
2950problem) after we copied the information to another matrix.
2952Using the Timer class, we collect some statistics that compare the runtime
2953of the block solver with the one from the problem implementation above.
2954Besides the solution with the two options we also check if the solutions
2955of the two variants are close to each other (i.e. this solver gives indeed the
2956same solution as we had before) and calculate the infinity
2957norm of the vector difference.
2963 Number
of degrees
of freedom: 679 (594+85) [0.00162792 s]
2965 Computing preconditioner... [0.0025959 s]
2973 Number
of degrees
of freedom: 1683 (1482+201) [0.00345707 s]
2975 Computing preconditioner... [0.00605702 s]
2983 Number
of degrees
of freedom: 3813 (3370+443) [0.00729299 s]
2985 Computing preconditioner... [0.0167508 s]
2993 Number
of degrees
of freedom: 8723 (7722+1001) [0.017709 s]
2995 Computing preconditioner... [0.0435679 s]
3003 Number
of degrees
of freedom: 19383 (17186+2197) [0.039988 s]
3005 Computing preconditioner... [0.118314 s]
3013 Number
of degrees
of freedom: 40855 (36250+4605) [0.0880702 s]
3015 Computing preconditioner... [0.278339 s]
3034 Number
of degrees
of freedom: 1356 (1275+81) [0.00845218 s]
3036 Computing preconditioner... [0.00712395 s]
3044 Number
of degrees
of freedom: 5088 (4827+261) [0.0346942 s]
3046 Computing preconditioner... [0.0465031 s]
3054 Number
of degrees
of freedom: 22406 (21351+1055) [0.175669 s]
3056 Computing preconditioner... [0.286435 s]
3064 Number
of degrees
of freedom: 93176 (89043+4133) [0.790985 s]
3074 Number
of degrees
of freedom: 327808 (313659+14149) [3.44995 s]
3084 Number
of degrees
of freedom: 1254464 (1201371+53093) [19.6795 s]
3100<a name=
"step_22-Combiningtheblockpreconditionerandmultigrid"></a><
h5>
Combining the block preconditioner
and multigrid</
h5>
3110<a name=
"step_22-Noblockmatricesandvectors"></a><
h5>
No block matrices
and vectors</
h5>
3139 ExcIndexRange (component, 0, this->n_components));
3158 Point<dim>(-2,-2,-1));
3169<table width=
"60%" align=
"center">
3172 <
img src=
"https://www.dealii.org/images/steps/developer/step-22.3d-extension.png" alt=
"">
3175 <
img src=
"https://www.dealii.org/images/steps/developer/step-22.3d-grid-extension.png" alt=
"">
3181<a name=
"step_22-PlainProg"></a>
static constexpr unsigned int rank
std::conditional_t< rank_==1, std::array< Number, dim >, std::array< Tensor< rank_ - 1, dim, Number >, dim > > values
#define Assert(cond, exc)
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())
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)
std::vector< index_type > data
void component_wise(DoFHandler< dim, spacedim > &dof_handler, const std::vector< unsigned int > &target_component=std::vector< unsigned int >())
void Cuthill_McKee(DoFHandler< dim, spacedim > &dof_handler, const bool reversed_numbering=false, const bool use_constraints=false, const std::vector< types::global_dof_index > &starting_indices=std::vector< types::global_dof_index >())
void subdivided_hyper_rectangle(Triangulation< dim, spacedim > &tria, const std::vector< unsigned int > &repetitions, const Point< dim > &p1, const Point< dim > &p2, const bool colorize=false)
@ matrix
Contents is actually a matrix.
@ symmetric
Matrix is symmetric.
@ diagonal
Matrix is diagonal.
constexpr types::blas_int one
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
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)
constexpr ReturnType< rank, T >::value_type & extract(T &t, const ArrayType &indices)
VectorType::value_type * end(VectorType &V)
int(&) functions(const void *v1, const void *v2)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
inline ::VectorizedArray< Number, width > atan(const ::VectorizedArray< Number, width > &x)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)