deal.II version GIT relicensing-1721-g8100761196 2024-08-31 12:30:00+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
step-22.h
Go to the documentation of this file.
1) const
1340 *   {
1341 *   return Tensor<1, dim>();
1342 *   }
1343 *  
1344 *  
1345 *   template <int dim>
1346 *   void RightHandSide<dim>::value_list(const std::vector<Point<dim>> &vp,
1347 *   std::vector<Tensor<1, dim>> &values) const
1348 *   {
1349 *   for (unsigned int c = 0; c < vp.size(); ++c)
1350 *   {
1351 *   values[c] = RightHandSide<dim>::value(vp[c]);
1352 *   }
1353 *   }
1354 *  
1355 *  
1356 * @endcode
1357 *
1358 *
1359 * <a name="step_22-Linearsolversandpreconditioners"></a>
1360 * <h3>Linear solvers and preconditioners</h3>
1361 *
1362
1363 *
1364 * The linear solvers and preconditioners are discussed extensively in the
1365 * introduction. Here, we create the respective objects that will be used.
1366 *
1367
1368 *
1369 *
1370 * <a name="step_22-ThecodeInverseMatrixcodeclasstemplate"></a>
1371 * <h4>The <code>InverseMatrix</code> class template</h4>
1372 * The <code>InverseMatrix</code> class represents the data structure for an
1373 * inverse matrix. Unlike @ref step_20 "step-20", we implement this with a class instead of
1374 * the helper function inverse_linear_operator() we will apply this class to
1375 * different kinds of matrices that will require different preconditioners
1376 * (in @ref step_20 "step-20" we only used a non-identity preconditioner for the mass
1377 * matrix). The types of matrix and preconditioner are passed to this class
1378 * via template parameters, and matrix and preconditioner objects of these
1379 * types will then be passed to the constructor when an
1380 * <code>InverseMatrix</code> object is created. The member function
1381 * <code>vmult</code> is obtained by solving a linear system:
1382 *
1383 * @code
1384 *   template <class MatrixType, class PreconditionerType>
1385 *   class InverseMatrix : public Subscriptor
1386 *   {
1387 *   public:
1388 *   InverseMatrix(const MatrixType &m,
1389 *   const PreconditionerType &preconditioner);
1390 *  
1391 *   void vmult(Vector<double> &dst, const Vector<double> &src) const;
1392 *  
1393 *   private:
1395 *   const SmartPointer<const PreconditionerType> preconditioner;
1396 *   };
1397 *  
1398 *  
1399 *   template <class MatrixType, class PreconditionerType>
1400 *   InverseMatrix<MatrixType, PreconditionerType>::InverseMatrix(
1401 *   const MatrixType &m,
1402 *   const PreconditionerType &preconditioner)
1403 *   : matrix(&m)
1404 *   , preconditioner(&preconditioner)
1405 *   {}
1406 *  
1407 *  
1408 * @endcode
1409 *
1410 * This is the implementation of the <code>vmult</code> function.
1411 *
1412
1413 *
1414 * In this class we use a rather large tolerance for the solver control. The
1415 * reason for this is that the function is used very frequently, and hence,
1416 * any additional effort to make the residual in the CG solve smaller makes
1417 * the solution more expensive. Note that we do not only use this class as a
1418 * preconditioner for the Schur complement, but also when forming the
1419 * inverse of the Laplace matrix &ndash; which is hence directly responsible
1420 * for the accuracy of the solution itself, so we can't choose a too large
1421 * tolerance, either.
1422 *
1423 * @code
1424 *   template <class MatrixType, class PreconditionerType>
1425 *   void InverseMatrix<MatrixType, PreconditionerType>::vmult(
1426 *   Vector<double> &dst,
1427 *   const Vector<double> &src) const
1428 *   {
1429 *   SolverControl solver_control(src.size(), 1e-6 * src.l2_norm());
1430 *   SolverCG<Vector<double>> cg(solver_control);
1431 *  
1432 *   dst = 0;
1433 *  
1434 *   cg.solve(*matrix, dst, src, *preconditioner);
1435 *   }
1436 *  
1437 *  
1438 * @endcode
1439 *
1440 *
1441 * <a name="step_22-ThecodeSchurComplementcodeclasstemplate"></a>
1442 * <h4>The <code>SchurComplement</code> class template</h4>
1443 *
1444
1445 *
1446 * This class implements the Schur complement discussed in the introduction.
1447 * It is in analogy to @ref step_20 "step-20". Though, we now call it with a template
1448 * parameter <code>PreconditionerType</code> in order to access that when
1449 * specifying the respective type of the inverse matrix class. As a
1450 * consequence of the definition above, the declaration
1451 * <code>InverseMatrix</code> now contains the second template parameter for
1452 * a preconditioner class as above, which affects the
1453 * <code>SmartPointer</code> object <code>m_inverse</code> as well.
1454 *
1455 * @code
1456 *   template <class PreconditionerType>
1457 *   class SchurComplement : public Subscriptor
1458 *   {
1459 *   public:
1460 *   SchurComplement(
1461 *   const BlockSparseMatrix<double> &system_matrix,
1462 *   const InverseMatrix<SparseMatrix<double>, PreconditionerType> &A_inverse);
1463 *  
1464 *   void vmult(Vector<double> &dst, const Vector<double> &src) const;
1465 *  
1466 *   private:
1467 *   const SmartPointer<const BlockSparseMatrix<double>> system_matrix;
1468 *   const SmartPointer<
1469 *   const InverseMatrix<SparseMatrix<double>, PreconditionerType>>
1470 *   A_inverse;
1471 *  
1472 *   mutable Vector<double> tmp1, tmp2;
1473 *   };
1474 *  
1475 *  
1476 *  
1477 *   template <class PreconditionerType>
1478 *   SchurComplement<PreconditionerType>::SchurComplement(
1479 *   const BlockSparseMatrix<double> &system_matrix,
1480 *   const InverseMatrix<SparseMatrix<double>, PreconditionerType> &A_inverse)
1481 *   : system_matrix(&system_matrix)
1482 *   , A_inverse(&A_inverse)
1483 *   , tmp1(system_matrix.block(0, 0).m())
1484 *   , tmp2(system_matrix.block(0, 0).m())
1485 *   {}
1486 *  
1487 *  
1488 *   template <class PreconditionerType>
1489 *   void
1490 *   SchurComplement<PreconditionerType>::vmult(Vector<double> &dst,
1491 *   const Vector<double> &src) const
1492 *   {
1493 *   system_matrix->block(0, 1).vmult(tmp1, src);
1494 *   A_inverse->vmult(tmp2, tmp1);
1495 *   system_matrix->block(1, 0).vmult(dst, tmp2);
1496 *   }
1497 *  
1498 *  
1499 * @endcode
1500 *
1501 *
1502 * <a name="step_22-StokesProblemclassimplementation"></a>
1503 * <h3>StokesProblem class implementation</h3>
1504 *
1505
1506 *
1507 *
1508 * <a name="step_22-StokesProblemStokesProblem"></a>
1509 * <h4>StokesProblem::StokesProblem</h4>
1510 *
1511
1512 *
1513 * The constructor of this class looks very similar to the one of
1514 * @ref step_20 "step-20". The constructor initializes the variables for the polynomial
1515 * degree, triangulation, finite element system and the dof handler. The
1516 * underlying polynomial functions are of order <code>degree+1</code> for
1517 * the vector-valued velocity components and of order <code>degree</code>
1518 * for the pressure. This gives the LBB-stable element pair
1519 * @f$Q_{degree+1}^d\times Q_{degree}@f$, often referred to as the Taylor-Hood
1520 * element.
1521 *
1522
1523 *
1524 * Note that we initialize the triangulation with a MeshSmoothing argument,
1525 * which ensures that the refinement of cells is done in a way that the
1526 * approximation of the PDE solution remains well-behaved (problems arise if
1527 * grids are too unstructured), see the documentation of
1528 * <code>Triangulation::MeshSmoothing</code> for details.
1529 *
1530 * @code
1531 *   template <int dim>
1532 *   StokesProblem<dim>::StokesProblem(const unsigned int degree)
1533 *   : degree(degree)
1534 *   , triangulation(Triangulation<dim>::maximum_smoothing)
1535 *   , fe(FE_Q<dim>(degree + 1) ^ dim, FE_Q<dim>(degree))
1536 *   , dof_handler(triangulation)
1537 *   {}
1538 *  
1539 *  
1540 * @endcode
1541 *
1542 *
1543 * <a name="step_22-StokesProblemsetup_dofs"></a>
1544 * <h4>StokesProblem::setup_dofs</h4>
1545 *
1546
1547 *
1548 * Given a mesh, this function associates the degrees of freedom with it and
1549 * creates the corresponding matrices and vectors. At the beginning it also
1550 * releases the pointer to the preconditioner object (if the shared pointer
1551 * pointed at anything at all at this point) since it will definitely not be
1552 * needed any more after this point and will have to be re-computed after
1553 * assembling the matrix, and unties the sparse matrices from their sparsity
1554 * pattern objects.
1555 *
1556
1557 *
1558 * We then proceed with distributing degrees of freedom and renumbering
1559 * them: In order to make the ILU preconditioner (in 3d) work efficiently,
1560 * it is important to enumerate the degrees of freedom in such a way that it
1561 * reduces the bandwidth of the matrix, or maybe more importantly: in such a
1562 * way that the ILU is as close as possible to a real LU decomposition. On
1563 * the other hand, we need to preserve the block structure of velocity and
1564 * pressure already seen in @ref step_20 "step-20" and @ref step_21 "step-21". This is done in two
1565 * steps: First, all dofs are renumbered to improve the ILU and then we
1566 * renumber once again by components. Since
1567 * <code>DoFRenumbering::component_wise</code> does not touch the
1568 * renumbering within the individual blocks, the basic renumbering from the
1569 * first step remains. As for how the renumber degrees of freedom to improve
1570 * the ILU: deal.II has a number of algorithms that attempt to find
1571 * orderings to improve ILUs, or reduce the bandwidth of matrices, or
1572 * optimize some other aspect. The DoFRenumbering namespace shows a
1573 * comparison of the results we obtain with several of these algorithms
1574 * based on the testcase discussed here in this tutorial program. Here, we
1575 * will use the traditional Cuthill-McKee algorithm already used in some of
1576 * the previous tutorial programs. In the <a href="#improved-ilu">section
1577 * on improved ILU</a> we're going to discuss this issue in more detail.
1578 *
1579
1580 *
1581 * There is one more change compared to previous tutorial programs: There is
1582 * no reason in sorting the <code>dim</code> velocity components
1583 * individually. In fact, rather than first enumerating all @f$x@f$-velocities,
1584 * then all @f$y@f$-velocities, etc, we would like to keep all velocities at the
1585 * same location together and only separate between velocities (all
1586 * components) and pressures. By default, this is not what the
1587 * DoFRenumbering::component_wise function does: it treats each vector
1588 * component separately; what we have to do is group several components into
1589 * "blocks" and pass this block structure to that function. Consequently, we
1590 * allocate a vector <code>block_component</code> with as many elements as
1591 * there are components and describe all velocity components to correspond
1592 * to block 0, while the pressure component will form block 1:
1593 *
1594 * @code
1595 *   template <int dim>
1596 *   void StokesProblem<dim>::setup_dofs()
1597 *   {
1598 *   A_preconditioner.reset();
1599 *   system_matrix.clear();
1600 *   preconditioner_matrix.clear();
1601 *  
1602 *   dof_handler.distribute_dofs(fe);
1603 *   DoFRenumbering::Cuthill_McKee(dof_handler);
1604 *  
1605 *   std::vector<unsigned int> block_component(dim + 1, 0);
1606 *   block_component[dim] = 1;
1607 *   DoFRenumbering::component_wise(dof_handler, block_component);
1608 *  
1609 * @endcode
1610 *
1611 * Now comes the implementation of Dirichlet boundary conditions, which
1612 * should be evident after the discussion in the introduction. All that
1613 * changed is that the function already appears in the setup functions,
1614 * whereas we were used to see it in some assembly routine. Further down
1615 * below where we set up the mesh, we will associate the top boundary
1616 * where we impose Dirichlet boundary conditions with boundary indicator
1617 * 1. We will have to pass this boundary indicator as second argument to
1618 * the function below interpolating boundary values. There is one more
1619 * thing, though. The function describing the Dirichlet conditions was
1620 * defined for all components, both velocity and pressure. However, the
1621 * Dirichlet conditions are to be set for the velocity only. To this end,
1622 * we use a ComponentMask that only selects the velocity components. The
1623 * component mask is obtained from the finite element by specifying the
1624 * particular components we want. Since we use adaptively refined grids,
1625 * the affine constraints object needs to be first filled with hanging node
1626 * constraints generated from the DoF handler. Note the order of the two
1627 * functions &mdash; we first compute the hanging node constraints, and
1628 * then insert the boundary values into the constraints object. This makes
1629 * sure that we respect H<sup>1</sup> conformity on boundaries with
1630 * hanging nodes (in three space dimensions), where the hanging node needs
1631 * to dominate the Dirichlet boundary values.
1632 *
1633 * @code
1634 *   {
1635 *   constraints.clear();
1636 *  
1637 *   const FEValuesExtractors::Vector velocities(0);
1638 *   DoFTools::make_hanging_node_constraints(dof_handler, constraints);
1640 *   1,
1641 *   BoundaryValues<dim>(),
1642 *   constraints,
1643 *   fe.component_mask(velocities));
1644 *   }
1645 *  
1646 *   constraints.close();
1647 *  
1648 * @endcode
1649 *
1650 * In analogy to @ref step_20 "step-20", we count the dofs in the individual components.
1651 * We could do this in the same way as there, but we want to operate on
1652 * the block structure we used already for the renumbering: The function
1653 * <code>DoFTools::count_dofs_per_fe_block</code> does the same as
1654 * <code>DoFTools::count_dofs_per_fe_component</code>, but now grouped as
1655 * velocity and pressure block via <code>block_component</code>.
1656 *
1657 * @code
1658 *   const std::vector<types::global_dof_index> dofs_per_block =
1659 *   DoFTools::count_dofs_per_fe_block(dof_handler, block_component);
1660 *   const types::global_dof_index n_u = dofs_per_block[0];
1661 *   const types::global_dof_index n_p = dofs_per_block[1];
1662 *  
1663 *   std::cout << " Number of active cells: " << triangulation.n_active_cells()
1664 *   << std::endl
1665 *   << " Number of degrees of freedom: " << dof_handler.n_dofs()
1666 *   << " (" << n_u << '+' << n_p << ')' << std::endl;
1667 *  
1668 * @endcode
1669 *
1670 * The next task is to allocate a sparsity pattern for the system matrix we
1671 * will create and one for the preconditioner matrix. We could do this in
1672 * the same way as in @ref step_20 "step-20", i.e. directly build an object of type
1673 * SparsityPattern through DoFTools::make_sparsity_pattern. However, there
1674 * is a major reason not to do so:
1675 * In 3d, the function DoFTools::max_couplings_between_dofs yields a
1676 * conservative but rather large number for the coupling between the
1677 * individual dofs, so that the memory initially provided for the creation
1678 * of the sparsity pattern of the matrix is far too much -- so much actually
1679 * that the initial sparsity pattern won't even fit into the physical memory
1680 * of most systems already for moderately-sized 3d problems, see also the
1681 * discussion in @ref step_18 "step-18". Instead, we first build temporary objects that use
1682 * a different data structure that doesn't require allocating more memory
1683 * than necessary but isn't suitable for use as a basis of SparseMatrix or
1684 * BlockSparseMatrix objects; in a second step we then copy these objects
1685 * into objects of type BlockSparsityPattern. This is entirely analogous to
1686 * what we already did in @ref step_11 "step-11" and @ref step_18 "step-18". In particular, we make use of
1687 * the fact that we will never write into the @f$(1,1)@f$ block of the system
1688 * matrix and that this is the only block to be filled for the
1689 * preconditioner matrix.
1690 *
1691
1692 *
1693 * All this is done inside new scopes, which means that the memory of
1694 * <code>dsp</code> will be released once the information has been copied to
1695 * <code>sparsity_pattern</code>.
1696 *
1697 * @code
1698 *   {
1699 *   BlockDynamicSparsityPattern dsp(dofs_per_block, dofs_per_block);
1700 *  
1701 *   Table<2, DoFTools::Coupling> coupling(dim + 1, dim + 1);
1702 *   for (unsigned int c = 0; c < dim + 1; ++c)
1703 *   for (unsigned int d = 0; d < dim + 1; ++d)
1704 *   if (!((c == dim) && (d == dim)))
1705 *   coupling[c][d] = DoFTools::always;
1706 *   else
1707 *   coupling[c][d] = DoFTools::none;
1708 *  
1709 *   DoFTools::make_sparsity_pattern(
1710 *   dof_handler, coupling, dsp, constraints, false);
1711 *  
1712 *   sparsity_pattern.copy_from(dsp);
1713 *   }
1714 *  
1715 *   {
1716 *   BlockDynamicSparsityPattern preconditioner_dsp(dofs_per_block,
1717 *   dofs_per_block);
1718 *  
1719 *   Table<2, DoFTools::Coupling> preconditioner_coupling(dim + 1, dim + 1);
1720 *   for (unsigned int c = 0; c < dim + 1; ++c)
1721 *   for (unsigned int d = 0; d < dim + 1; ++d)
1722 *   if (((c == dim) && (d == dim)))
1723 *   preconditioner_coupling[c][d] = DoFTools::always;
1724 *   else
1725 *   preconditioner_coupling[c][d] = DoFTools::none;
1726 *  
1727 *   DoFTools::make_sparsity_pattern(dof_handler,
1728 *   preconditioner_coupling,
1729 *   preconditioner_dsp,
1730 *   constraints,
1731 *   false);
1732 *  
1733 *   preconditioner_sparsity_pattern.copy_from(preconditioner_dsp);
1734 *   }
1735 *  
1736 * @endcode
1737 *
1738 * Finally, the system matrix, the preconsitioner matrix, the solution and
1739 * the right hand side vector are created from the block structure similar
1740 * to the approach in @ref step_20 "step-20":
1741 *
1742 * @code
1743 *   system_matrix.reinit(sparsity_pattern);
1744 *   preconditioner_matrix.reinit(preconditioner_sparsity_pattern);
1745 *  
1746 *   solution.reinit(dofs_per_block);
1747 *   system_rhs.reinit(dofs_per_block);
1748 *   }
1749 *  
1750 *  
1751 * @endcode
1752 *
1753 *
1754 * <a name="step_22-StokesProblemassemble_system"></a>
1755 * <h4>StokesProblem::assemble_system</h4>
1756 *
1757
1758 *
1759 * The assembly process follows the discussion in @ref step_20 "step-20" and in the
1760 * introduction. We use the well-known abbreviations for the data structures
1761 * that hold the local matrices, right hand side, and global numbering of the
1762 * degrees of freedom for the present cell.
1763 *
1764 * @code
1765 *   template <int dim>
1766 *   void StokesProblem<dim>::assemble_system()
1767 *   {
1768 *   system_matrix = 0;
1769 *   system_rhs = 0;
1770 *   preconditioner_matrix = 0;
1771 *  
1772 *   const QGauss<dim> quadrature_formula(degree + 2);
1773 *  
1774 *   FEValues<dim> fe_values(fe,
1775 *   quadrature_formula,
1776 *   update_values | update_quadrature_points |
1777 *   update_JxW_values | update_gradients);
1778 *  
1779 *   const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
1780 *  
1781 *   const unsigned int n_q_points = quadrature_formula.size();
1782 *  
1783 *   FullMatrix<double> local_matrix(dofs_per_cell, dofs_per_cell);
1784 *   FullMatrix<double> local_preconditioner_matrix(dofs_per_cell,
1785 *   dofs_per_cell);
1786 *   Vector<double> local_rhs(dofs_per_cell);
1787 *  
1788 *   std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1789 *  
1790 *   const RightHandSide<dim> right_hand_side;
1791 *   std::vector<Tensor<1, dim>> rhs_values(n_q_points, Tensor<1, dim>());
1792 *  
1793 * @endcode
1794 *
1795 * Next, we need two objects that work as extractors for the FEValues
1796 * object. Their use is explained in detail in the report on @ref
1797 * vector_valued :
1798 *
1799 * @code
1800 *   const FEValuesExtractors::Vector velocities(0);
1801 *   const FEValuesExtractors::Scalar pressure(dim);
1802 *  
1803 * @endcode
1804 *
1805 * As an extension over @ref step_20 "step-20" and @ref step_21 "step-21", we include a few optimizations
1806 * that make assembly much faster for this particular problem. The
1807 * improvements are based on the observation that we do a few calculations
1808 * too many times when we do as in @ref step_20 "step-20": The symmetric gradient actually
1809 * has <code>dofs_per_cell</code> different values per quadrature point, but
1810 * we extract it <code>dofs_per_cell*dofs_per_cell</code> times from the
1811 * FEValues object - for both the loop over <code>i</code> and the inner
1812 * loop over <code>j</code>. In 3d, that means evaluating it @f$89^2=7921@f$
1813 * instead of @f$89@f$ times, a not insignificant difference.
1814 *
1815
1816 *
1817 * So what we're going to do here is to avoid such repeated calculations
1818 * by getting a vector of rank-2 tensors (and similarly for the divergence
1819 * and the basis function value on pressure) at the quadrature point prior
1820 * to starting the loop over the dofs on the cell. First, we create the
1821 * respective objects that will hold these values. Then, we start the loop
1822 * over all cells and the loop over the quadrature points, where we first
1823 * extract these values. There is one more optimization we implement here:
1824 * the local matrix (as well as the global one) is going to be symmetric,
1825 * since all the operations involved are symmetric with respect to @f$i@f$ and
1826 * @f$j@f$. This is implemented by simply running the inner loop not to
1827 * <code>dofs_per_cell</code>, but only up to <code>i</code>, the index of
1828 * the outer loop.
1829 *
1830 * @code
1831 *   std::vector<SymmetricTensor<2, dim>> symgrad_phi_u(dofs_per_cell);
1832 *   std::vector<double> div_phi_u(dofs_per_cell);
1833 *   std::vector<Tensor<1, dim>> phi_u(dofs_per_cell);
1834 *   std::vector<double> phi_p(dofs_per_cell);
1835 *  
1836 *   for (const auto &cell : dof_handler.active_cell_iterators())
1837 *   {
1838 *   fe_values.reinit(cell);
1839 *  
1840 *   local_matrix = 0;
1841 *   local_preconditioner_matrix = 0;
1842 *   local_rhs = 0;
1843 *  
1844 *   right_hand_side.value_list(fe_values.get_quadrature_points(),
1845 *   rhs_values);
1846 *  
1847 *   for (unsigned int q = 0; q < n_q_points; ++q)
1848 *   {
1849 *   for (unsigned int k = 0; k < dofs_per_cell; ++k)
1850 *   {
1851 *   symgrad_phi_u[k] =
1852 *   fe_values[velocities].symmetric_gradient(k, q);
1853 *   div_phi_u[k] = fe_values[velocities].divergence(k, q);
1854 *   phi_u[k] = fe_values[velocities].value(k, q);
1855 *   phi_p[k] = fe_values[pressure].value(k, q);
1856 *   }
1857 *  
1858 * @endcode
1859 *
1860 * Now finally for the bilinear forms of both the system matrix and
1861 * the matrix we use for the preconditioner. Recall that the
1862 * formulas for these two are
1863 * @f{align*}{
1864 * A_{ij} &= a(\varphi_i,\varphi_j)
1865 * \\ &= \underbrace{2(\varepsilon(\varphi_{i,\textbf{u}}),
1866 * \varepsilon(\varphi_{j,\textbf{u}}))_{\Omega}}
1867 * _{(1)}
1868 * \;
1869 * \underbrace{- (\textrm{div}\; \varphi_{i,\textbf{u}},
1870 * \varphi_{j,p})_{\Omega}}
1871 * _{(2)}
1872 * \;
1873 * \underbrace{- (\varphi_{i,p},
1874 * \textrm{div}\;
1875 * \varphi_{j,\textbf{u}})_{\Omega}}
1876 * _{(3)}
1877 * @f}
1878 * and
1879 * @f{align*}{
1880 * M_{ij} &= \underbrace{(\varphi_{i,p},
1881 * \varphi_{j,p})_{\Omega}}
1882 * _{(4)},
1883 * @f}
1884 * respectively, where @f$\varphi_{i,\textbf{u}}@f$ and @f$\varphi_{i,p}@f$
1885 * are the velocity and pressure components of the @f$i@f$th shape
1886 * function. The various terms above are then easily recognized in
1887 * the following implementation:
1888 *
1889 * @code
1890 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
1891 *   {
1892 *   for (unsigned int j = 0; j <= i; ++j)
1893 *   {
1894 *   local_matrix(i, j) +=
1895 *   (2 * (symgrad_phi_u[i] * symgrad_phi_u[j]) // (1)
1896 *   - div_phi_u[i] * phi_p[j] // (2)
1897 *   - phi_p[i] * div_phi_u[j]) // (3)
1898 *   * fe_values.JxW(q); // * dx
1899 *  
1900 *   local_preconditioner_matrix(i, j) +=
1901 *   (phi_p[i] * phi_p[j]) // (4)
1902 *   * fe_values.JxW(q); // * dx
1903 *   }
1904 * @endcode
1905 *
1906 * Note that in the implementation of (1) above, `operator*`
1907 * is overloaded for symmetric tensors, yielding the scalar
1908 * product between the two tensors.
1909 *
1910
1911 *
1912 * For the right-hand side, we need to multiply the (vector of)
1913 * velocity shape functions with the vector of body force
1914 * right-hand side components, both evaluated at the current
1915 * quadrature point. We have implemented the body forces as a
1916 * `TensorFunction<1,dim>`, so its values at quadrature points
1917 * are already tensors for which the application of `operator*`
1918 * against the velocity components of the shape function results
1919 * in the dot product, as intended.
1920 *
1921 * @code
1922 *   local_rhs(i) += phi_u[i] // phi_u_i(x_q)
1923 *   * rhs_values[q] // * f(x_q)
1924 *   * fe_values.JxW(q); // * dx
1925 *   }
1926 *   }
1927 *  
1928 * @endcode
1929 *
1930 * Before we can write the local data into the global matrix (and
1931 * simultaneously use the AffineConstraints object to apply
1932 * Dirichlet boundary conditions and eliminate hanging node constraints,
1933 * as we discussed in the introduction), we have to be careful about one
1934 * thing, though. We have only built half of the local matrices
1935 * because of symmetry, but we're going to save the full matrices
1936 * in order to use the standard functions for solving. This is done
1937 * by flipping the indices in case we are pointing into the empty part
1938 * of the local matrices.
1939 *
1940 * @code
1941 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
1942 *   for (unsigned int j = i + 1; j < dofs_per_cell; ++j)
1943 *   {
1944 *   local_matrix(i, j) = local_matrix(j, i);
1945 *   local_preconditioner_matrix(i, j) =
1946 *   local_preconditioner_matrix(j, i);
1947 *   }
1948 *  
1949 *   cell->get_dof_indices(local_dof_indices);
1950 *   constraints.distribute_local_to_global(local_matrix,
1951 *   local_rhs,
1952 *   local_dof_indices,
1953 *   system_matrix,
1954 *   system_rhs);
1955 *   constraints.distribute_local_to_global(local_preconditioner_matrix,
1956 *   local_dof_indices,
1957 *   preconditioner_matrix);
1958 *   }
1959 *  
1960 * @endcode
1961 *
1962 * Before we're going to solve this linear system, we generate a
1963 * preconditioner for the velocity-velocity matrix, i.e.,
1964 * <code>block(0,0)</code> in the system matrix. As mentioned above, this
1965 * depends on the spatial dimension. Since the two classes described by
1966 * the <code>InnerPreconditioner::type</code> alias have the same
1967 * interface, we do not have to do anything different whether we want to
1968 * use a sparse direct solver or an ILU:
1969 *
1970 * @code
1971 *   std::cout << " Computing preconditioner..." << std::endl << std::flush;
1972 *  
1973 *   A_preconditioner =
1974 *   std::make_shared<typename InnerPreconditioner<dim>::type>();
1975 *   A_preconditioner->initialize(
1976 *   system_matrix.block(0, 0),
1977 *   typename InnerPreconditioner<dim>::type::AdditionalData());
1978 *   }
1979 *  
1980 *  
1981 *  
1982 * @endcode
1983 *
1984 *
1985 * <a name="step_22-StokesProblemsolve"></a>
1986 * <h4>StokesProblem::solve</h4>
1987 *
1988
1989 *
1990 * After the discussion in the introduction and the definition of the
1991 * respective classes above, the implementation of the <code>solve</code>
1992 * function is rather straight-forward and done in a similar way as in
1993 * @ref step_20 "step-20". To start with, we need an object of the
1994 * <code>InverseMatrix</code> class that represents the inverse of the
1995 * matrix A. As described in the introduction, the inverse is generated with
1996 * the help of an inner preconditioner of type
1997 * <code>InnerPreconditioner::type</code>.
1998 *
1999 * @code
2000 *   template <int dim>
2001 *   void StokesProblem<dim>::solve()
2002 *   {
2003 *   const InverseMatrix<SparseMatrix<double>,
2004 *   typename InnerPreconditioner<dim>::type>
2005 *   A_inverse(system_matrix.block(0, 0), *A_preconditioner);
2006 *   Vector<double> tmp(solution.block(0).size());
2007 *  
2008 * @endcode
2009 *
2010 * This is as in @ref step_20 "step-20". We generate the right hand side @f$B A^{-1} F - G@f$
2011 * for the Schur complement and an object that represents the respective
2012 * linear operation @f$B A^{-1} B^T@f$, now with a template parameter
2013 * indicating the preconditioner - in accordance with the definition of
2014 * the class.
2015 *
2016 * @code
2017 *   {
2018 *   Vector<double> schur_rhs(solution.block(1).size());
2019 *   A_inverse.vmult(tmp, system_rhs.block(0));
2020 *   system_matrix.block(1, 0).vmult(schur_rhs, tmp);
2021 *   schur_rhs -= system_rhs.block(1);
2022 *  
2023 *   SchurComplement<typename InnerPreconditioner<dim>::type> schur_complement(
2024 *   system_matrix, A_inverse);
2025 *  
2026 * @endcode
2027 *
2028 * The usual control structures for the solver call are created...
2029 *
2030 * @code
2031 *   SolverControl solver_control(solution.block(1).size(),
2032 *   1e-6 * schur_rhs.l2_norm());
2033 *   SolverCG<Vector<double>> cg(solver_control);
2034 *  
2035 * @endcode
2036 *
2037 * Now to the preconditioner to the Schur complement. As explained in
2038 * the introduction, the preconditioning is done by a @ref GlossMassMatrix "mass matrix" in the
2039 * pressure variable.
2040 *
2041
2042 *
2043 * Actually, the solver needs to have the preconditioner in the form
2044 * @f$P^{-1}@f$, so we need to create an inverse operation. Once again, we
2045 * use an object of the class <code>InverseMatrix</code>, which
2046 * implements the <code>vmult</code> operation that is needed by the
2047 * solver. In this case, we have to invert the pressure mass matrix. As
2048 * it already turned out in earlier tutorial programs, the inversion of
2049 * a mass matrix is a rather cheap and straight-forward operation
2050 * (compared to, e.g., a Laplace matrix). The CG method with ILU
2051 * preconditioning converges in 5-10 steps, independently on the mesh
2052 * size. This is precisely what we do here: We choose another ILU
2053 * preconditioner and take it along to the InverseMatrix object via the
2054 * corresponding template parameter. A CG solver is then called within
2055 * the vmult operation of the inverse matrix.
2056 *
2057
2058 *
2059 * An alternative that is cheaper to build, but needs more iterations
2060 * afterwards, would be to choose a SSOR preconditioner with factor
2061 * 1.2. It needs about twice the number of iterations, but the costs for
2062 * its generation are almost negligible.
2063 *
2064 * @code
2065 *   SparseILU<double> preconditioner;
2066 *   preconditioner.initialize(preconditioner_matrix.block(1, 1),
2067 *   SparseILU<double>::AdditionalData());
2068 *  
2069 *   InverseMatrix<SparseMatrix<double>, SparseILU<double>> m_inverse(
2070 *   preconditioner_matrix.block(1, 1), preconditioner);
2071 *  
2072 * @endcode
2073 *
2074 * With the Schur complement and an efficient preconditioner at hand, we
2075 * can solve the respective equation for the pressure (i.e. block 0 in
2076 * the solution vector) in the usual way:
2077 *
2078 * @code
2079 *   cg.solve(schur_complement, solution.block(1), schur_rhs, m_inverse);
2080 *  
2081 * @endcode
2082 *
2083 * After this first solution step, the hanging node constraints have to
2084 * be distributed to the solution in order to achieve a consistent
2085 * pressure field.
2086 *
2087 * @code
2088 *   constraints.distribute(solution);
2089 *  
2090 *   std::cout << " " << solver_control.last_step()
2091 *   << " outer CG Schur complement iterations for pressure"
2092 *   << std::endl;
2093 *   }
2094 *  
2095 * @endcode
2096 *
2097 * As in @ref step_20 "step-20", we finally need to solve for the velocity equation where
2098 * we plug in the solution to the pressure equation. This involves only
2099 * objects we already know - so we simply multiply @f$p@f$ by @f$B^T@f$, subtract
2100 * the right hand side and multiply by the inverse of @f$A@f$. At the end, we
2101 * need to distribute the constraints from hanging nodes in order to
2102 * obtain a consistent flow field:
2103 *
2104 * @code
2105 *   {
2106 *   system_matrix.block(0, 1).vmult(tmp, solution.block(1));
2107 *   tmp *= -1;
2108 *   tmp += system_rhs.block(0);
2109 *  
2110 *   A_inverse.vmult(solution.block(0), tmp);
2111 *  
2112 *   constraints.distribute(solution);
2113 *   }
2114 *   }
2115 *  
2116 *  
2117 * @endcode
2118 *
2119 *
2120 * <a name="step_22-StokesProblemoutput_results"></a>
2121 * <h4>StokesProblem::output_results</h4>
2122 *
2123
2124 *
2125 * The next function generates graphical output. In this example, we are
2126 * going to use the VTK file format. We attach names to the individual
2127 * variables in the problem: <code>velocity</code> to the <code>dim</code>
2128 * components of velocity and <code>pressure</code> to the pressure.
2129 *
2130
2131 *
2132 * Not all visualization programs have the ability to group individual
2133 * vector components into a vector to provide vector plots; in particular,
2134 * this holds for some VTK-based visualization programs. In this case, the
2135 * logical grouping of components into vectors should already be described
2136 * in the file containing the data. In other words, what we need to do is
2137 * provide our output writers with a way to know which of the components of
2138 * the finite element logically form a vector (with @f$d@f$ components in @f$d@f$
2139 * space dimensions) rather than letting them assume that we simply have a
2140 * bunch of scalar fields. This is achieved using the members of the
2141 * <code>DataComponentInterpretation</code> namespace: as with the filename,
2142 * we create a vector in which the first <code>dim</code> components refer
2143 * to the velocities and are given the tag
2144 * DataComponentInterpretation::component_is_part_of_vector; we
2145 * finally push one tag
2146 * DataComponentInterpretation::component_is_scalar to describe
2147 * the grouping of the pressure variable.
2148 *
2149
2150 *
2151 * The rest of the function is then the same as in @ref step_20 "step-20".
2152 *
2153 * @code
2154 *   template <int dim>
2155 *   void
2156 *   StokesProblem<dim>::output_results(const unsigned int refinement_cycle) const
2157 *   {
2158 *   std::vector<std::string> solution_names(dim, "velocity");
2159 *   solution_names.emplace_back("pressure");
2160 *  
2161 *   std::vector<DataComponentInterpretation::DataComponentInterpretation>
2162 *   data_component_interpretation(
2163 *   dim, DataComponentInterpretation::component_is_part_of_vector);
2164 *   data_component_interpretation.push_back(
2165 *   DataComponentInterpretation::component_is_scalar);
2166 *  
2167 *   DataOut<dim> data_out;
2168 *   data_out.attach_dof_handler(dof_handler);
2169 *   data_out.add_data_vector(solution,
2170 *   solution_names,
2171 *   DataOut<dim>::type_dof_data,
2172 *   data_component_interpretation);
2173 *   data_out.build_patches();
2174 *  
2175 *   std::ofstream output(
2176 *   "solution-" + Utilities::int_to_string(refinement_cycle, 2) + ".vtk");
2177 *   data_out.write_vtk(output);
2178 *   }
2179 *  
2180 *  
2181 * @endcode
2182 *
2183 *
2184 * <a name="step_22-StokesProblemrefine_mesh"></a>
2185 * <h4>StokesProblem::refine_mesh</h4>
2186 *
2187
2188 *
2189 * This is the last interesting function of the <code>StokesProblem</code>
2190 * class. As indicated by its name, it takes the solution to the problem
2191 * and refines the mesh where this is needed. The procedure is the same as
2192 * in the respective step in @ref step_6 "step-6", with the exception that we base the
2193 * refinement only on the change in pressure, i.e., we call the Kelly error
2194 * estimator with a mask object of type ComponentMask that selects the
2195 * single scalar component for the pressure that we are interested in (we
2196 * get such a mask from the finite element class by specifying the component
2197 * we want). Additionally, we do not coarsen the grid again:
2198 *
2199 * @code
2200 *   template <int dim>
2201 *   void StokesProblem<dim>::refine_mesh()
2202 *   {
2203 *   Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
2204 *  
2205 *   const FEValuesExtractors::Scalar pressure(dim);
2206 *   KellyErrorEstimator<dim>::estimate(
2207 *   dof_handler,
2208 *   QGauss<dim - 1>(degree + 1),
2209 *   std::map<types::boundary_id, const Function<dim> *>(),
2210 *   solution,
2211 *   estimated_error_per_cell,
2212 *   fe.component_mask(pressure));
2213 *  
2214 *   GridRefinement::refine_and_coarsen_fixed_number(triangulation,
2215 *   estimated_error_per_cell,
2216 *   0.3,
2217 *   0.0);
2218 *   triangulation.execute_coarsening_and_refinement();
2219 *   }
2220 *  
2221 *  
2222 * @endcode
2223 *
2224 *
2225 * <a name="step_22-StokesProblemrun"></a>
2226 * <h4>StokesProblem::run</h4>
2227 *
2228
2229 *
2230 * The last step in the Stokes class is, as usual, the function that
2231 * generates the initial grid and calls the other functions in the
2232 * respective order.
2233 *
2234
2235 *
2236 * We start off with a rectangle of size @f$4 \times 1@f$ (in 2d) or @f$4 \times 1
2237 * \times 1@f$ (in 3d), placed in @f$R^2/R^3@f$ as @f$(-2,2)\times(-1,0)@f$ or
2238 * @f$(-2,2)\times(0,1)\times(-1,0)@f$, respectively. It is natural to start
2239 * with equal mesh size in each direction, so we subdivide the initial
2240 * rectangle four times in the first coordinate direction. To limit the
2241 * scope of the variables involved in the creation of the mesh to the range
2242 * where we actually need them, we put the entire block between a pair of
2243 * braces:
2244 *
2245 * @code
2246 *   template <int dim>
2247 *   void StokesProblem<dim>::run()
2248 *   {
2249 *   {
2250 *   std::vector<unsigned int> subdivisions(dim, 1);
2251 *   subdivisions[0] = 4;
2252 *  
2253 *   const Point<dim> bottom_left = (dim == 2 ?
2254 *   Point<dim>(-2, -1) : // 2d case
2255 *   Point<dim>(-2, 0, -1)); // 3d case
2256 *  
2257 *   const Point<dim> top_right = (dim == 2 ?
2258 *   Point<dim>(2, 0) : // 2d case
2259 *   Point<dim>(2, 1, 0)); // 3d case
2260 *  
2261 *   GridGenerator::subdivided_hyper_rectangle(triangulation,
2262 *   subdivisions,
2263 *   bottom_left,
2264 *   top_right);
2265 *   }
2266 *  
2267 * @endcode
2268 *
2269 * A boundary indicator of 1 is set to all boundaries that are subject to
2270 * Dirichlet boundary conditions, i.e. to faces that are located at 0 in
2271 * the last coordinate direction. See the example description above for
2272 * details.
2273 *
2274 * @code
2275 *   for (const auto &cell : triangulation.active_cell_iterators())
2276 *   for (const auto &face : cell->face_iterators())
2277 *   if (face->center()[dim - 1] == 0)
2278 *   face->set_all_boundary_ids(1);
2279 *  
2280 *  
2281 * @endcode
2282 *
2283 * We then apply an initial refinement before solving for the first
2284 * time. In 3d, there are going to be more degrees of freedom, so we
2285 * refine less there:
2286 *
2287 * @code
2288 *   triangulation.refine_global(4 - dim);
2289 *  
2290 * @endcode
2291 *
2292 * As first seen in @ref step_6 "step-6", we cycle over the different refinement levels
2293 * and refine (except for the first cycle), setup the degrees of freedom
2294 * and matrices, assemble, solve and create output:
2295 *
2296 * @code
2297 *   for (unsigned int refinement_cycle = 0; refinement_cycle < 6;
2298 *   ++refinement_cycle)
2299 *   {
2300 *   std::cout << "Refinement cycle " << refinement_cycle << std::endl;
2301 *  
2302 *   if (refinement_cycle > 0)
2303 *   refine_mesh();
2304 *  
2305 *   setup_dofs();
2306 *  
2307 *   std::cout << " Assembling..." << std::endl << std::flush;
2308 *   assemble_system();
2309 *  
2310 *   std::cout << " Solving..." << std::flush;
2311 *   solve();
2312 *  
2313 *   output_results(refinement_cycle);
2314 *  
2315 *   std::cout << std::endl;
2316 *   }
2317 *   }
2318 *   } // namespace Step22
2319 *  
2320 *  
2321 * @endcode
2322 *
2323 *
2324 * <a name="step_22-Thecodemaincodefunction"></a>
2325 * <h3>The <code>main</code> function</h3>
2326 *
2327
2328 *
2329 * The main function is the same as in @ref step_20 "step-20". We pass the element degree as
2330 * a parameter and choose the space dimension at the well-known template slot.
2331 *
2332 * @code
2333 *   int main()
2334 *   {
2335 *   try
2336 *   {
2337 *   using namespace Step22;
2338 *  
2339 *   StokesProblem<2> flow_problem(1);
2340 *   flow_problem.run();
2341 *   }
2342 *   catch (std::exception &exc)
2343 *   {
2344 *   std::cerr << std::endl
2345 *   << std::endl
2346 *   << "----------------------------------------------------"
2347 *   << std::endl;
2348 *   std::cerr << "Exception on processing: " << std::endl
2349 *   << exc.what() << std::endl
2350 *   << "Aborting!" << std::endl
2351 *   << "----------------------------------------------------"
2352 *   << std::endl;
2353 *  
2354 *   return 1;
2355 *   }
2356 *   catch (...)
2357 *   {
2358 *   std::cerr << std::endl
2359 *   << std::endl
2360 *   << "----------------------------------------------------"
2361 *   << std::endl;
2362 *   std::cerr << "Unknown exception!" << std::endl
2363 *   << "Aborting!" << std::endl
2364 *   << "----------------------------------------------------"
2365 *   << std::endl;
2366 *   return 1;
2367 *   }
2368 *  
2369 *   return 0;
2370 *   }
2371 * @endcode
2372<a name="step_22-Results"></a><h1>Results</h1>
2373
2374
2375<a name="step_22-Outputoftheprogramandgraphicalvisualization"></a><h3>Output of the program and graphical visualization</h3>
2376
2377
2378<a name="step_22-2Dcalculations"></a><h4>2D calculations</h4>
2379
2380
2381Running the program with the space dimension set to 2 in the <code>main</code>
2382function yields the following output (in "release mode",
2383See also <a href="https://www.math.colostate.edu/~bangerth/videos.676.18.html">video lecture 18</a>.):
2384@code
2385examples/step-22> make run
2386Refinement cycle 0
2387 Number of active cells: 64
2388 Number of degrees of freedom: 679 (594+85)
2389 Assembling...
2390 Computing preconditioner...
2391 Solving... 11 outer CG Schur complement iterations for pressure
2392
2393Refinement cycle 1
2394 Number of active cells: 160
2395 Number of degrees of freedom: 1683 (1482+201)
2396 Assembling...
2397 Computing preconditioner...
2398 Solving... 11 outer CG Schur complement iterations for pressure
2399
2400Refinement cycle 2
2401 Number of active cells: 376
2402 Number of degrees of freedom: 3813 (3370+443)
2403 Assembling...
2404 Computing preconditioner...
2405 Solving... 11 outer CG Schur complement iterations for pressure
2406
2407Refinement cycle 3
2408 Number of active cells: 880
2409 Number of degrees of freedom: 8723 (7722+1001)
2410 Assembling...
2411 Computing preconditioner...
2412 Solving... 11 outer CG Schur complement iterations for pressure
2413
2414Refinement cycle 4
2415 Number of active cells: 2008
2416 Number of degrees of freedom: 19383 (17186+2197)
2417 Assembling...
2418 Computing preconditioner...
2419 Solving... 11 outer CG Schur complement iterations for pressure
2420
2421Refinement cycle 5
2422 Number of active cells: 4288
2423 Number of degrees of freedom: 40855 (36250+4605)
2424 Assembling...
2425 Computing preconditioner...
2426 Solving... 11 outer CG Schur complement iterations for pressure
2427@endcode
2428
2429The entire computation above takes about 2 seconds on a reasonably
2430quick (for 2015 standards) machine.
2431
2432What we see immediately from this is that the number of (outer)
2433iterations does not increase as we refine the mesh. This confirms the
2434statement in the introduction that preconditioning the Schur
2435complement with the mass matrix indeed yields a matrix spectrally
2436equivalent to the identity matrix (i.e. with eigenvalues bounded above
2437and below independently of the mesh size or the relative sizes of
2438cells). In other words, the mass matrix and the Schur complement are
2439spectrally equivalent.
2440
2441In the images below, we show the grids for the first six refinement
2442steps in the program. Observe how the grid is refined in regions
2443where the solution rapidly changes: On the upper boundary, we have
2444Dirichlet boundary conditions that are -1 in the left half of the line
2445and 1 in the right one, so there is an abrupt change at @f$x=0@f$. Likewise,
2446there are changes from Dirichlet to Neumann data in the two upper
2447corners, so there is need for refinement there as well:
2448
2449<table width="60%" align="center">
2450 <tr>
2451 <td align="center">
2452 <img src="https://www.dealii.org/images/steps/developer/step-22.2d.mesh-0.png" alt="">
2453 </td>
2454 <td align="center">
2455 <img src="https://www.dealii.org/images/steps/developer/step-22.2d.mesh-1.png" alt="">
2456 </td>
2457 </tr>
2458 <tr>
2459 <td align="center">
2460 <img src="https://www.dealii.org/images/steps/developer/step-22.2d.mesh-2.png" alt="">
2461 </td>
2462 <td align="center">
2463 <img src="https://www.dealii.org/images/steps/developer/step-22.2d.mesh-3.png" alt="">
2464 </td>
2465 </tr>
2466 <tr>
2467 <td align="center">
2468 <img src="https://www.dealii.org/images/steps/developer/step-22.2d.mesh-4.png" alt="">
2469 </td>
2470 <td align="center">
2471 <img src="https://www.dealii.org/images/steps/developer/step-22.2d.mesh-5.png" alt="">
2472 </td>
2473 </tr>
2474</table>
2475
2476Finally, following is a plot of the flow field. It shows fluid
2477transported along with the moving upper boundary and being replaced by
2478material coming from below:
2479
2480<img src="https://www.dealii.org/images/steps/developer/step-22.2d.solution.png" alt="">
2481
2482This plot uses the capability of VTK-based visualization programs (in
2483this case of VisIt) to show vector data; this is the result of us
2484declaring the velocity components of the finite element in use to be a
2485set of vector components, rather than independent scalar components in
2486the <code>StokesProblem@<dim@>::%output_results</code> function of this
2487tutorial program.
2488
2489
2490
2491<a name="step_22-3Dcalculations"></a><h4>3D calculations</h4>
2492
2493
2494In 3d, the screen output of the program looks like this:
2495
2496@code
2497Refinement cycle 0
2498 Number of active cells: 32
2499 Number of degrees of freedom: 1356 (1275+81)
2500 Assembling...
2501 Computing preconditioner...
2502 Solving... 13 outer CG Schur complement iterations for pressure.
2503
2504Refinement cycle 1
2505 Number of active cells: 144
2506 Number of degrees of freedom: 5088 (4827+261)
2507 Assembling...
2508 Computing preconditioner...
2509 Solving... 14 outer CG Schur complement iterations for pressure.
2510
2511Refinement cycle 2
2512 Number of active cells: 704
2513 Number of degrees of freedom: 22406 (21351+1055)
2514 Assembling...
2515 Computing preconditioner...
2516 Solving... 14 outer CG Schur complement iterations for pressure.
2517
2518Refinement cycle 3
2519 Number of active cells: 3168
2520 Number of degrees of freedom: 93176 (89043+4133)
2521 Assembling...
2522 Computing preconditioner...
2523 Solving... 15 outer CG Schur complement iterations for pressure.
2524
2525Refinement cycle 4
2526 Number of active cells: 11456
2527 Number of degrees of freedom: 327808 (313659+14149)
2528 Assembling...
2529 Computing preconditioner...
2530 Solving... 15 outer CG Schur complement iterations for pressure.
2531
2532Refinement cycle 5
2533 Number of active cells: 45056
2534 Number of degrees of freedom: 1254464 (1201371+53093)
2535 Assembling...
2536 Computing preconditioner...
2537 Solving... 14 outer CG Schur complement iterations for pressure.
2538@endcode
2539
2540Again, we see that the number of outer iterations does not increase as
2541we refine the mesh. Nevertheless, the compute time increases
2542significantly: for each of the iterations above separately, it takes about
25430.14 seconds, 0.63 seconds, 4.8 seconds, 35 seconds, 2 minutes and 33 seconds,
2544and 13 minutes and 12 seconds. This overall superlinear (in the number of
2545unknowns) increase in runtime is due to the fact that our inner solver is not
2546@f${\cal O}(N)@f$: a simple experiment shows that as we keep refining the mesh, the
2547average number of ILU-preconditioned CG iterations to invert the
2548velocity-velocity block @f$A@f$ increases.
2549
2550We will address the question of how possibly to improve our solver <a
2551href="#improved-solver">below</a>.
2552
2553As for the graphical output, the grids generated during the solution
2554look as follow:
2555
2556<table width="60%" align="center">
2557 <tr>
2558 <td align="center">
2559 <img src="https://www.dealii.org/images/steps/developer/step-22.3d.mesh-0.png" alt="">
2560 </td>
2561 <td align="center">
2562 <img src="https://www.dealii.org/images/steps/developer/step-22.3d.mesh-1.png" alt="">
2563 </td>
2564 </tr>
2565 <tr>
2566 <td align="center">
2567 <img src="https://www.dealii.org/images/steps/developer/step-22.3d.mesh-2.png" alt="">
2568 </td>
2569 <td align="center">
2570 <img src="https://www.dealii.org/images/steps/developer/step-22.3d.mesh-3.png" alt="">
2571 </td>
2572 </tr>
2573 <tr>
2574 <td align="center">
2575 <img src="https://www.dealii.org/images/steps/developer/step-22.3d.mesh-4.png" alt="">
2576 </td>
2577 <td align="center">
2578 <img src="https://www.dealii.org/images/steps/developer/step-22.3d.mesh-5.png" alt="">
2579 </td>
2580 </tr>
2581</table>
2582
2583Again, they show essentially the location of singularities introduced
2584by boundary conditions. The vector field computed makes for an
2585interesting graph:
2586
2587<img src="https://www.dealii.org/images/steps/developer/step-22.3d.solution.png" alt="">
2588
2589The isocontours shown here as well are those of the pressure
2590variable, showing the singularity at the point of discontinuous
2591velocity boundary conditions.
2592
2593
2594
2595<a name="step_22-Sparsitypattern"></a><h3>Sparsity pattern</h3>
2596
2597
2598As explained during the generation of the sparsity pattern, it is
2599important to have the numbering of degrees of freedom in mind when
2600using preconditioners like incomplete LU decompositions. This is most
2601conveniently visualized using the distribution of nonzero elements in
2602the @ref GlossStiffnessMatrix "stiffness matrix".
2603
2604If we don't do anything special to renumber degrees of freedom (i.e.,
2605without using DoFRenumbering::Cuthill_McKee, but with using
2606DoFRenumbering::component_wise to ensure that degrees of freedom are
2607appropriately sorted into their corresponding blocks of the matrix and
2608vector), then we get the following image after the first adaptive
2609refinement in two dimensions:
2610
2611<img src="https://www.dealii.org/images/steps/developer/step-22.2d.sparsity-nor.png" alt="">
2612
2613In order to generate such a graph, you have to insert a piece of
2614code like the following to the end of the setup step.
2615@code
2616 {
2617 std::ofstream out ("sparsity_pattern.gpl");
2618 sparsity_pattern.print_gnuplot(out);
2619 }
2620@endcode
2621
2622It is clearly visible that the nonzero entries are spread over almost the
2623whole matrix. This makes preconditioning by ILU inefficient: ILU generates a
2624Gaussian elimination (LU decomposition) without fill-in elements, which means
2625that more tentative fill-ins left out will result in a worse approximation of
2626the complete decomposition.
2627
2628In this program, we have thus chosen a more advanced renumbering of
2629components. The renumbering with DoFRenumbering::Cuthill_McKee and grouping
2630the components into velocity and pressure yields the following output:
2631
2632<img src="https://www.dealii.org/images/steps/developer/step-22.2d.sparsity-ren.png" alt="">
2633
2634It is apparent that the situation has improved a lot. Most of the elements are
2635now concentrated around the diagonal in the (0,0) block in the matrix. Similar
2636effects are also visible for the other blocks. In this case, the ILU
2637decomposition will be much closer to the full LU decomposition, which improves
2638the quality of the preconditioner. (It may be interesting to note that the
2639sparse direct solver UMFPACK does some %internal renumbering of the equations
2640before actually generating a sparse LU decomposition; that procedure leads to
2641a very similar pattern to the one we got from the Cuthill-McKee algorithm.)
2642
2643Finally, we want to have a closer
2644look at a sparsity pattern in 3D. We show only the (0,0) block of the
2645matrix, again after one adaptive refinement. Apart from the fact that the matrix
2646size has increased, it is also visible that there are many more entries
2647in the matrix. Moreover, even for the optimized renumbering, there will be a
2648considerable amount of tentative fill-in elements. This illustrates why UMFPACK
2649is not a good choice in 3D - a full decomposition needs many new entries that
2650 eventually won't fit into the physical memory (RAM):
2651
2652<img src="https://www.dealii.org/images/steps/developer/step-22.3d.sparsity_uu-ren.png" alt="">
2653
2654
2655
2656<a name="step_22-Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
2657
2658
2659<a name="step-22-improved-solver">
2660<a name="step_22-Improvedlinearsolverin3D"></a><h4>Improved linear solver in 3D</h4>
2661
2662</a>
2663
2664We have seen in the section of computational results that the number of outer
2665iterations does not depend on the mesh size, which is optimal in a sense of
2666scalability. This does, however, not apply to the solver as a whole, as
2667mentioned above:
2668We did not look at the number of inner iterations when generating the inverse of
2669the matrix @f$A@f$ and the mass matrix @f$M_p@f$. Of course, this is unproblematic in
2670the 2D case where we precondition @f$A@f$ with a direct solver and the
2671<code>vmult</code> operation of the inverse matrix structure will converge in
2672one single CG step, but this changes in 3D where we only use an ILU
2673preconditioner. There, the number of required preconditioned CG steps to
2674invert @f$A@f$ increases as the mesh is refined, and each <code>vmult</code>
2675operation involves on average approximately 14, 23, 36, 59, 75 and 101 inner
2676CG iterations in the refinement steps shown above. (On the other hand,
2677the number of iterations for applying the inverse pressure mass matrix is
2678always around five, both in two and three dimensions.) To summarize, most work
2679is spent on solving linear systems with the same matrix @f$A@f$ over and over again.
2680What makes this look even worse is the fact that we
2681actually invert a matrix that is about 95 percent the size of the total system
2682matrix and stands for 85 percent of the non-zero entries in the sparsity
2683pattern. Hence, the natural question is whether it is reasonable to solve a
2684linear system with matrix @f$A@f$ for about 15 times when calculating the solution
2685to the block system.
2686
2687The answer is, of course, that we can do that in a few other (most of the time
2688better) ways.
2689Nevertheless, it has to be remarked that an indefinite system as the one
2690at hand puts indeed much higher
2691demands on the linear algebra than standard elliptic problems as we have seen
2692in the early tutorial programs. The improvements are still rather
2693unsatisfactory, if one compares with an elliptic problem of similar
2694size. Either way, we will introduce below a number of improvements to the
2695linear solver, a discussion that we will re-consider again with additional
2696options in the @ref step_31 "step-31" program.
2697
2698<a name="step-22-improved-ilu">
2699<a name="step_22-BetterILUdecompositionbysmartreordering"></a><h5>Better ILU decomposition by smart reordering</h5>
2700
2701</a>
2702A first attempt to improve the speed of the linear solution process is to choose
2703a dof reordering that makes the ILU being closer to a full LU decomposition, as
2704already mentioned in the in-code comments. The DoFRenumbering namespace compares
2705several choices for the renumbering of dofs for the Stokes equations. The best
2706result regarding the computing time was found for the King ordering, which is
2707accessed through the call DoFRenumbering::boost::king_ordering. With that
2708program, the inner solver needs considerably less operations, e.g. about 62
2709inner CG iterations for the inversion of @f$A@f$ at cycle 4 compared to about 75
2710iterations with the standard Cuthill-McKee-algorithm. Also, the computing time
2711at cycle 4 decreased from about 17 to 11 minutes for the <code>solve()</code>
2712call. However, the King ordering (and the orderings provided by the
2713DoFRenumbering::boost namespace in general) has a serious drawback - it uses
2714much more memory than the in-build deal versions, since it acts on abstract
2715graphs rather than the geometry provided by the triangulation. In the present
2716case, the renumbering takes about 5 times as much memory, which yields an
2717infeasible algorithm for the last cycle in 3D with 1.2 million
2718unknowns.
2719
2720<a name="step_22-BetterpreconditionerfortheinnerCGsolver"></a><h5>Better preconditioner for the inner CG solver</h5>
2721
2722Another idea to improve the situation even more would be to choose a
2723preconditioner that makes CG for the (0,0) matrix @f$A@f$ converge in a
2724mesh-independent number of iterations, say 10 to 30. We have seen such a
2725candidate in @ref step_16 "step-16": multigrid.
2726
2727<a name="step_22-BlockSchurcomplementpreconditioner"></a><h5>Block Schur complement preconditioner</h5>
2728
2729<a name="step-22-block-schur"></a>
2730Even with a good preconditioner for @f$A@f$, we still
2731need to solve of the same linear system repeatedly (with different
2732right hand sides, though) in order to make the Schur complement solve
2733converge. The approach we are going to discuss here is how inner iteration
2734and outer iteration can be combined. If we persist in calculating the Schur
2735complement, there is no other possibility.
2736
2737The alternative is to attack the block system at once and use an approximate
2738Schur complement as efficient preconditioner. The idea is as
2739follows: If we find a block preconditioner @f$P@f$ such that the matrix
2740@f{eqnarray*}{
2741 P^{-1}\left(\begin{array}{cc}
2742 A & B^T \\ B & 0
2743 \end{array}\right)
2744@f}
2745is simple, then an iterative solver with that preconditioner will converge in a
2746few iterations. Using the Schur complement @f$S = B A^{-1} B^T@f$, one finds that
2747@f{eqnarray*}{
2748 P^{-1}
2749 =
2750 \left(\begin{array}{cc}
2751 A^{-1} & 0 \\ S^{-1} B A^{-1} & -S^{-1}
2752 \end{array}\right)
2753@f}
2754would appear to be a good choice since
2755@f{eqnarray*}{
2756 P^{-1}\left(\begin{array}{cc}
2757 A & B^T \\ B & 0
2758 \end{array}\right)
2759 =
2760 \left(\begin{array}{cc}
2761 A^{-1} & 0 \\ S^{-1} B A^{-1} & -S^{-1}
2762 \end{array}\right)\cdot \left(\begin{array}{cc}
2763 A & B^T \\ B & 0
2764 \end{array}\right)
2765 =
2766 \left(\begin{array}{cc}
2767 I & A^{-1} B^T \\ 0 & I
2768 \end{array}\right).
2769@f}
2770This is the approach taken by the paper by Silvester and Wathen referenced
2771to in the introduction (with the exception that Silvester and Wathen use
2772right preconditioning). In this case, a Krylov-based iterative method would
2773converge in one step only if exact inverses of @f$A@f$ and @f$S@f$ were applied,
2774since all the eigenvalues are one (and the number of iterations in such a
2775method is bounded by the number of distinct eigenvalues). Below, we will
2776discuss the choice of an adequate solver for this problem. First, we are
2777going to have a closer look at the implementation of the preconditioner.
2778
2779Since @f$P@f$ is aimed to be a preconditioner only, we shall use approximations to
2780the inverse of the Schur complement @f$S@f$ and the matrix @f$A@f$. Hence, the Schur
2781complement will be approximated by the pressure mass matrix @f$M_p@f$, and we use
2782a preconditioner to @f$A@f$ (without an InverseMatrix class around it) for
2783approximating @f$A^{-1}@f$.
2784
2785Here comes the class that implements the block Schur
2786complement preconditioner. The <code>vmult</code> operation for block vectors
2787according to the derivation above can be specified by three successive
2788operations:
2789@code
2790template <class PreconditionerA, class PreconditionerMp>
2791class BlockSchurPreconditioner : public Subscriptor
2792{
2793 public:
2794 BlockSchurPreconditioner (const BlockSparseMatrix<double> &S,
2795 const InverseMatrix<SparseMatrix<double>,PreconditionerMp> &Mpinv,
2796 const PreconditionerA &Apreconditioner);
2797
2798 void vmult (BlockVector<double> &dst,
2799 const BlockVector<double> &src) const;
2800
2801 private:
2804 PreconditionerMp > > m_inverse;
2805 const PreconditionerA &a_preconditioner;
2806
2807 mutable Vector<double> tmp;
2808
2809};
2810
2811template <class PreconditionerA, class PreconditionerMp>
2812BlockSchurPreconditioner<PreconditionerA, PreconditionerMp>::BlockSchurPreconditioner(
2814 const InverseMatrix<SparseMatrix<double>,PreconditionerMp> &Mpinv,
2815 const PreconditionerA &Apreconditioner
2816 )
2817 :
2818 system_matrix (&S),
2819 m_inverse (&Mpinv),
2820 a_preconditioner (Apreconditioner),
2821 tmp (S.block(1,1).m())
2822{}
2823
2824 // Now the interesting function, the multiplication of
2825 // the preconditioner with a BlockVector.
2826template <class PreconditionerA, class PreconditionerMp>
2827void BlockSchurPreconditioner<PreconditionerA, PreconditionerMp>::vmult (
2829 const BlockVector<double> &src) const
2830{
2831 // Form u_new = A^{-1} u
2832 a_preconditioner.vmult (dst.block(0), src.block(0));
2833 // Form tmp = - B u_new + p
2834 // (<code>SparseMatrix::residual</code>
2835 // does precisely this)
2836 system_matrix->block(1,0).residual(tmp, dst.block(0), src.block(1));
2837 // Change sign in tmp
2838 tmp *= -1;
2839 // Multiply by approximate Schur complement
2840 // (i.e. a pressure mass matrix)
2841 m_inverse->vmult (dst.block(1), tmp);
2842}
2843@endcode
2844
2845Since we act on the whole block system now, we have to live with one
2846disadvantage: we need to perform the solver iterations on
2847the full block system instead of the smaller pressure space.
2848
2849Now we turn to the question which solver we should use for the block
2850system. The first observation is that the resulting preconditioned matrix cannot
2851be solved with CG since it is neither positive definite nor symmetric.
2852
2853The deal.II libraries implement several solvers that are appropriate for the
2854problem at hand. One choice is the solver @ref SolverBicgstab "BiCGStab", which
2855was used for the solution of the unsymmetric advection problem in @ref step_9 "step-9". The
2856second option, the one we are going to choose, is @ref SolverGMRES "GMRES"
2857(generalized minimum residual). Both methods have their pros and cons - there
2858are problems where one of the two candidates clearly outperforms the other, and
2859vice versa.
2860<a href="http://en.wikipedia.org/wiki/GMRES#Comparison_with_other_solvers">Wikipedia</a>'s
2861article on the GMRES method gives a comparative presentation.
2862A more comprehensive and well-founded comparison can be read e.g. in the book by
2863J.W. Demmel (Applied Numerical Linear Algebra, SIAM, 1997, section 6.6.6).
2864
2865For our specific problem with the ILU preconditioner for @f$A@f$, we certainly need
2866to perform hundreds of iterations on the block system for large problem sizes
2867(we won't beat CG!). Actually, this disfavors GMRES: During the GMRES
2868iterations, a basis of Krylov vectors is successively built up and some
2869operations are performed on these vectors. The more vectors are in this basis,
2870the more operations and memory will be needed. The number of operations scales
2871as @f${\cal O}(n + k^2)@f$ and memory as @f${\cal O}(kn)@f$, where @f$k@f$ is the number of
2872vectors in the Krylov basis and @f$n@f$ the size of the (block) matrix.
2873To not let these demands grow excessively, deal.II limits the size @f$k@f$ of the
2874basis to 30 vectors by default.
2875Then, the basis is rebuilt. This implementation of the GMRES method is called
2876GMRES(k), with default @f$k=30@f$. What we have gained by this restriction,
2877namely a bound on operations and memory requirements, will be compensated by
2878the fact that we use an incomplete basis - this will increase the number of
2879required iterations.
2880
2881BiCGStab, on the other hand, won't get slower when many iterations are needed
2882(one iteration uses only results from one preceding step and
2883not all the steps as GMRES). Besides the fact the BiCGStab is more expensive per
2884step since two matrix-vector products are needed (compared to one for
2885CG or GMRES), there is one main reason which makes BiCGStab not appropriate for
2886this problem: The preconditioner applies the inverse of the pressure
2887mass matrix by using the InverseMatrix class. Since the application of the
2888inverse matrix to a vector is done only in approximative way (an exact inverse
2889is too expensive), this will also affect the solver. In the case of BiCGStab,
2890the Krylov vectors will not be orthogonal due to that perturbation. While
2891this is uncritical for a small number of steps (up to about 50), it ruins the
2892performance of the solver when these perturbations have grown to a significant
2893magnitude in the coarse of iterations.
2894
2895We did some experiments with BiCGStab and found it to
2896be faster than GMRES up to refinement cycle 3 (in 3D), but it became very slow
2897for cycles 4 and 5 (even slower than the original Schur complement), so the
2898solver is useless in this situation. Choosing a sharper tolerance for the
2899inverse matrix class (<code>1e-10*src.l2_norm()</code> instead of
2900<code>1e-6*src.l2_norm()</code>) made BiCGStab perform well also for cycle 4,
2901but did not change the failure on the very large problems.
2902
2903GMRES is of course also effected by the approximate inverses, but it is not as
2904sensitive to orthogonality and retains a relatively good performance also for
2905large sizes, see the results below.
2906
2907With this said, we turn to the realization of the solver call with GMRES with
2908@f$k=100@f$ temporary vectors:
2909
2910@code
2911 const SparseMatrix<double> &pressure_mass_matrix
2912 = preconditioner_matrix.block(1,1);
2913 SparseILU<double> pmass_preconditioner;
2914 pmass_preconditioner.initialize (pressure_mass_matrix,
2915 SparseILU<double>::AdditionalData());
2916
2917 InverseMatrix<SparseMatrix<double>,SparseILU<double> >
2918 m_inverse (pressure_mass_matrix, pmass_preconditioner);
2919
2920 BlockSchurPreconditioner<typename InnerPreconditioner<dim>::type,
2921 SparseILU<double> >
2922 preconditioner (system_matrix, m_inverse, *A_preconditioner);
2923
2924 SolverControl solver_control (system_matrix.m(),
2925 1e-6*system_rhs.l2_norm());
2926 GrowingVectorMemory<BlockVector<double> > vector_memory;
2927 SolverGMRES<BlockVector<double> >::AdditionalData gmres_data;
2928 gmres_data.max_basis_size = 100;
2929
2930 SolverGMRES<BlockVector<double> > gmres(solver_control, vector_memory,
2931 gmres_data);
2932
2933 gmres.solve(system_matrix, solution, system_rhs,
2934 preconditioner);
2935
2936 constraints.distribute (solution);
2937
2938 std::cout << " "
2939 << solver_control.last_step()
2940 << " block GMRES iterations";
2941@endcode
2942
2943Obviously, one needs to add the include file @ref SolverGMRES
2944"<lac/solver_gmres.h>" in order to make this run.
2945We call the solver with a BlockVector template in order to enable
2946GMRES to operate on block vectors and matrices.
2947Note also that we need to set the (1,1) block in the system
2948matrix to zero (we saved the pressure mass matrix there which is not part of the
2949problem) after we copied the information to another matrix.
2950
2951Using the Timer class, we collect some statistics that compare the runtime
2952of the block solver with the one from the problem implementation above.
2953Besides the solution with the two options we also check if the solutions
2954of the two variants are close to each other (i.e. this solver gives indeed the
2955same solution as we had before) and calculate the infinity
2956norm of the vector difference.
2957
2958Let's first see the results in 2D:
2959@code
2960Refinement cycle 0
2961 Number of active cells: 64
2962 Number of degrees of freedom: 679 (594+85) [0.00162792 s]
2963 Assembling... [0.00108981 s]
2964 Computing preconditioner... [0.0025959 s]
2965 Solving...
2966 Schur complement: 11 outer CG iterations for p [0.00479603s ]
2967 Block Schur preconditioner: 12 GMRES iterations [0.00441718 s]
2968 l_infinity difference between solution vectors: 5.38258e-07
2969
2970Refinement cycle 1
2971 Number of active cells: 160
2972 Number of degrees of freedom: 1683 (1482+201) [0.00345707 s]
2973 Assembling... [0.00237417 s]
2974 Computing preconditioner... [0.00605702 s]
2975 Solving...
2976 Schur complement: 11 outer CG iterations for p [0.0123992s ]
2977 Block Schur preconditioner: 12 GMRES iterations [0.011909 s]
2978 l_infinity difference between solution vectors: 1.74658e-05
2979
2980Refinement cycle 2
2981 Number of active cells: 376
2982 Number of degrees of freedom: 3813 (3370+443) [0.00729299 s]
2983 Assembling... [0.00529909 s]
2984 Computing preconditioner... [0.0167508 s]
2985 Solving...
2986 Schur complement: 11 outer CG iterations for p [0.031672s ]
2987 Block Schur preconditioner: 12 GMRES iterations [0.029232 s]
2988 l_infinity difference between solution vectors: 7.81569e-06
2989
2990Refinement cycle 3
2991 Number of active cells: 880
2992 Number of degrees of freedom: 8723 (7722+1001) [0.017709 s]
2993 Assembling... [0.0126002 s]
2994 Computing preconditioner... [0.0435679 s]
2995 Solving...
2996 Schur complement: 11 outer CG iterations for p [0.0971651s ]
2997 Block Schur preconditioner: 12 GMRES iterations [0.0992041 s]
2998 l_infinity difference between solution vectors: 1.87249e-05
2999
3000Refinement cycle 4
3001 Number of active cells: 2008
3002 Number of degrees of freedom: 19383 (17186+2197) [0.039988 s]
3003 Assembling... [0.028281 s]
3004 Computing preconditioner... [0.118314 s]
3005 Solving...
3006 Schur complement: 11 outer CG iterations for p [0.252133s ]
3007 Block Schur preconditioner: 13 GMRES iterations [0.269125 s]
3008 l_infinity difference between solution vectors: 6.38657e-05
3009
3010Refinement cycle 5
3011 Number of active cells: 4288
3012 Number of degrees of freedom: 40855 (36250+4605) [0.0880702 s]
3013 Assembling... [0.0603511 s]
3014 Computing preconditioner... [0.278339 s]
3015 Solving...
3016 Schur complement: 11 outer CG iterations for p [0.53846s ]
3017 Block Schur preconditioner: 13 GMRES iterations [0.578667 s]
3018 l_infinity difference between solution vectors: 0.000173363
3019@endcode
3020
3021We see that there is no huge difference in the solution time between the
3022block Schur complement preconditioner solver and the Schur complement
3023itself. The reason is simple: we used a direct solve as preconditioner for
3024@f$A@f$ - so we cannot expect any gain by avoiding the inner iterations. We see
3025that the number of iterations has slightly increased for GMRES, but all in
3026all the two choices are fairly similar.
3027
3028The picture of course changes in 3D:
3029
3030@code
3031Refinement cycle 0
3032 Number of active cells: 32
3033 Number of degrees of freedom: 1356 (1275+81) [0.00845218 s]
3034 Assembling... [0.019372 s]
3035 Computing preconditioner... [0.00712395 s]
3036 Solving...
3037 Schur complement: 13 outer CG iterations for p [0.0320101s ]
3038 Block Schur preconditioner: 22 GMRES iterations [0.0048759 s]
3039 l_infinity difference between solution vectors: 2.15942e-05
3040
3041Refinement cycle 1
3042 Number of active cells: 144
3043 Number of degrees of freedom: 5088 (4827+261) [0.0346942 s]
3044 Assembling... [0.0857739 s]
3045 Computing preconditioner... [0.0465031 s]
3046 Solving...
3047 Schur complement: 14 outer CG iterations for p [0.349258s ]
3048 Block Schur preconditioner: 35 GMRES iterations [0.048759 s]
3049 l_infinity difference between solution vectors: 1.77657e-05
3050
3051Refinement cycle 2
3052 Number of active cells: 704
3053 Number of degrees of freedom: 22406 (21351+1055) [0.175669 s]
3054 Assembling... [0.437447 s]
3055 Computing preconditioner... [0.286435 s]
3056 Solving...
3057 Schur complement: 14 outer CG iterations for p [3.65519s ]
3058 Block Schur preconditioner: 63 GMRES iterations [0.497787 s]
3059 l_infinity difference between solution vectors: 5.08078e-05
3060
3061Refinement cycle 3
3062 Number of active cells: 3168
3063 Number of degrees of freedom: 93176 (89043+4133) [0.790985 s]
3064 Assembling... [1.97598 s]
3065 Computing preconditioner... [1.4325 s]
3066 Solving...
3067 Schur complement: 15 outer CG iterations for p [29.9666s ]
3068 Block Schur preconditioner: 128 GMRES iterations [5.02645 s]
3069 l_infinity difference between solution vectors: 0.000119671
3070
3071Refinement cycle 4
3072 Number of active cells: 11456
3073 Number of degrees of freedom: 327808 (313659+14149) [3.44995 s]
3074 Assembling... [7.54772 s]
3075 Computing preconditioner... [5.46306 s]
3076 Solving...
3077 Schur complement: 15 outer CG iterations for p [139.987s ]
3078 Block Schur preconditioner: 255 GMRES iterations [38.0946 s]
3079 l_infinity difference between solution vectors: 0.00020793
3080
3081Refinement cycle 5
3082 Number of active cells: 45056
3083 Number of degrees of freedom: 1254464 (1201371+53093) [19.6795 s]
3084 Assembling... [28.6586 s]
3085 Computing preconditioner... [22.401 s]
3086 Solving...
3087 Schur complement: 14 outer CG iterations for p [796.767s ]
3088 Block Schur preconditioner: 524 GMRES iterations [355.597 s]
3089 l_infinity difference between solution vectors: 0.000501219
3090@endcode
3091
3092Here, the block preconditioned solver is clearly superior to the Schur
3093complement, but the advantage gets less for more mesh points. This is
3094because GMRES(k) scales worse with the problem size than CG, as we discussed
3095above. Nonetheless, the improvement by a factor of 3-6 for moderate problem
3096sizes is quite impressive.
3097
3098
3099<a name="step_22-Combiningtheblockpreconditionerandmultigrid"></a><h5>Combining the block preconditioner and multigrid</h5>
3100
3101An ultimate linear solver for this problem could be imagined as a
3102combination of an optimal
3103preconditioner for @f$A@f$ (e.g. multigrid) and the block preconditioner
3104described above, which is the approach taken in the @ref step_31 "step-31"
3105and @ref step_32 "step-32" tutorial programs (where we use an algebraic multigrid
3106method) and @ref step_56 "step-56" (where we use a geometric multigrid method).
3107
3108
3109<a name="step_22-Noblockmatricesandvectors"></a><h5>No block matrices and vectors</h5>
3110
3111Another possibility that can be taken into account is to not set up a block
3112system, but rather solve the system of velocity and pressure all at once. The
3113options are direct solve with UMFPACK (2D) or GMRES with ILU
3114preconditioning (3D). It should be straightforward to try that.
3115
3116
3117
3118<a name="step_22-Moreinterestingtestcases"></a><h4>More interesting testcases</h4>
3119
3120
3121The program can of course also serve as a basis to compute the flow in more
3122interesting cases. The original motivation to write this program was for it to
3123be a starting point for some geophysical flow problems, such as the
3124movement of magma under places where continental plates drift apart (for
3125example mid-ocean ridges). Of course, in such places, the geometry is more
3126complicated than the examples shown above, but it is not hard to accommodate
3127for that.
3128
3129For example, by using the following modification of the boundary values
3130function
3131@code
3132template <int dim>
3133double
3134BoundaryValues<dim>::value (const Point<dim> &p,
3135 const unsigned int component) const
3136{
3137 Assert (component < this->n_components,
3138 ExcIndexRange (component, 0, this->n_components));
3139
3140 const double x_offset = std::atan(p[1]*4)/3;
3141
3142 if (component == 0)
3143 return (p[0] < x_offset ? -1 : (p[0] > x_offset ? 1 : 0));
3144 return 0;
3145}
3146@endcode
3147and the following way to generate the mesh as the domain
3148@f$[-2,2]\times[-2,2]\times[-1,0]@f$
3149@code
3150 std::vector<unsigned int> subdivisions (dim, 1);
3151 subdivisions[0] = 4;
3152 if (dim>2)
3153 subdivisions[1] = 4;
3154
3155 const Point<dim> bottom_left = (dim == 2 ?
3156 Point<dim>(-2,-1) :
3157 Point<dim>(-2,-2,-1));
3158 const Point<dim> top_right = (dim == 2 ?
3159 Point<dim>(2,0) :
3160 Point<dim>(2,2,0));
3161
3163 subdivisions,
3164 bottom_left,
3165 top_right);
3166@endcode
3167then we get images where the fault line is curved:
3168<table width="60%" align="center">
3169 <tr>
3170 <td align="center">
3171 <img src="https://www.dealii.org/images/steps/developer/step-22.3d-extension.png" alt="">
3172 </td>
3173 <td align="center">
3174 <img src="https://www.dealii.org/images/steps/developer/step-22.3d-grid-extension.png" alt="">
3175 </td>
3176 </tr>
3177</table>
3178 *
3179 *
3180<a name="step_22-PlainProg"></a>
3181<h1> The plain program</h1>
3182@include "step-22.cc"
3183*/
BlockType & block(const unsigned int i)
Definition point.h:111
Point< 2 > second
Definition grid_out.cc:4624
Point< 2 > first
Definition grid_out.cc:4623
#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())
Definition loop.h:564
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)
const Event initial
Definition event.cc:64
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 >())
std::vector< types::global_dof_index > count_dofs_per_fe_block(const DoFHandler< dim, spacedim > &dof, const std::vector< unsigned int > &target_block=std::vector< unsigned int >())
std::vector< types::global_dof_index > count_dofs_per_fe_component(const DoFHandler< dim, spacedim > &dof_handler, const bool vector_valued_once=false, const std::vector< unsigned int > &target_component={})
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.
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:191
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
constexpr ReturnType< rank, T >::value_type & extract(T &t, const ArrayType &indices)
VectorType::value_type * end(VectorType &V)
void interpolate_boundary_values(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const std::map< types::boundary_id, const Function< spacedim, number > * > &function_map, std::map< types::global_dof_index, number > &boundary_values, const ComponentMask &component_mask={})
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)
Definition types.h:32
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)