Reference documentation for deal.II version GIT relicensing-214-g6e74dec06b 2024-03-27 18:10:01+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-31.h
Go to the documentation of this file.
1,
1270 *   const unsigned int /*component*/ = 0) const override
1271 *   {
1272 *   return 0;
1273 *   }
1274 *  
1275 *   virtual void vector_value(const Point<dim> &p,
1276 *   Vector<double> &value) const override
1277 *   {
1278 *   for (unsigned int c = 0; c < this->n_components; ++c)
1279 *   value(c) = TemperatureInitialValues<dim>::value(p, c);
1280 *   }
1281 *   };
1282 *  
1283 *  
1284 *  
1285 *   template <int dim>
1286 *   class TemperatureRightHandSide : public Function<dim>
1287 *   {
1288 *   public:
1289 *   TemperatureRightHandSide()
1290 *   : Function<dim>(1)
1291 *   {}
1292 *  
1293 *   virtual double value(const Point<dim> &p,
1294 *   const unsigned int component = 0) const override
1295 *   {
1296 *   (void)component;
1297 *   Assert(component == 0,
1298 *   ExcMessage("Invalid operation for a scalar function."));
1299 *  
1300 *   Assert((dim == 2) || (dim == 3), ExcNotImplemented());
1301 *  
1302 *   static const Point<dim> source_centers[3] = {
1303 *   (dim == 2 ? Point<dim>(.3, .1) : Point<dim>(.3, .5, .1)),
1304 *   (dim == 2 ? Point<dim>(.45, .1) : Point<dim>(.45, .5, .1)),
1305 *   (dim == 2 ? Point<dim>(.75, .1) : Point<dim>(.75, .5, .1))};
1306 *   static const double source_radius = (dim == 2 ? 1. / 32 : 1. / 8);
1307 *  
1308 *   return ((source_centers[0].distance(p) < source_radius) ||
1309 *   (source_centers[1].distance(p) < source_radius) ||
1310 *   (source_centers[2].distance(p) < source_radius) ?
1311 *   1 :
1312 *   0);
1313 *   }
1314 *  
1315 *   virtual void vector_value(const Point<dim> &p,
1316 *   Vector<double> &value) const override
1317 *   {
1318 *   for (unsigned int c = 0; c < this->n_components; ++c)
1319 *   value(c) = TemperatureRightHandSide<dim>::value(p, c);
1320 *   }
1321 *   };
1322 *   } // namespace EquationData
1323 *  
1324 *  
1325 *  
1326 * @endcode
1327 *
1328 *
1329 * <a name="step_31-Linearsolversandpreconditioners"></a>
1330 * <h3>Linear solvers and preconditioners</h3>
1331 *
1332
1333 *
1334 * This section introduces some objects that are used for the solution of
1335 * the linear equations of the Stokes system that we need to solve in each
1336 * time step. Many of the ideas used here are the same as in @ref step_20 "step-20", where
1337 * Schur complement based preconditioners and solvers have been introduced,
1338 * with the actual interface taken from @ref step_22 "step-22" (in particular the
1339 * discussion in the "Results" section of @ref step_22 "step-22", in which we introduce
1340 * alternatives to the direct Schur complement approach). Note, however,
1341 * that here we don't use the Schur complement to solve the Stokes
1342 * equations, though an approximate Schur complement (the mass matrix on the
1343 * pressure space) appears in the preconditioner.
1344 *
1345 * @code
1346 *   namespace LinearSolvers
1347 *   {
1348 * @endcode
1349 *
1350 *
1351 * <a name="step_31-ThecodeInverseMatrixcodeclasstemplate"></a>
1352 * <h4>The <code>InverseMatrix</code> class template</h4>
1353 *
1354
1355 *
1356 * This class is an interface to calculate the action of an "inverted"
1357 * matrix on a vector (using the <code>vmult</code> operation) in the same
1358 * way as the corresponding class in @ref step_22 "step-22": when the product of an
1359 * object of this class is requested, we solve a linear equation system
1360 * with that matrix using the CG method, accelerated by a preconditioner
1361 * of (templated) class <code>PreconditionerType</code>.
1362 *
1363
1364 *
1365 * In a minor deviation from the implementation of the same class in
1366 * @ref step_22 "step-22", we make the <code>vmult</code> function take any
1367 * kind of vector type (it will yield compiler errors, however, if the
1368 * matrix does not allow a matrix-vector product with this kind of
1369 * vector).
1370 *
1371
1372 *
1373 * Secondly, we catch any exceptions that the solver may have thrown. The
1374 * reason is as follows: When debugging a program like this one
1375 * occasionally makes a mistake of passing an indefinite or nonsymmetric
1376 * matrix or preconditioner to the current class. The solver will, in that
1377 * case, not converge and throw a run-time exception. If not caught here
1378 * it will propagate up the call stack and may end up in
1379 * <code>main()</code> where we output an error message that will say that
1380 * the CG solver failed. The question then becomes: Which CG solver? The
1381 * one that inverted the mass matrix? The one that inverted the top left
1382 * block with the Laplace operator? Or a CG solver in one of the several
1383 * other nested places where we use linear solvers in the current code? No
1384 * indication about this is present in a run-time exception because it
1385 * doesn't store the stack of calls through which we got to the place
1386 * where the exception was generated.
1387 *
1388
1389 *
1390 * So rather than letting the exception propagate freely up to
1391 * <code>main()</code>, we acknowledge that in the current context
1392 * there is little that an outer function can do if the inner
1393 * solver fails. As a consequence, we deal with the situation by
1394 * catching the exception and letting the program fail by
1395 * triggering an assertion with a `false` condition (which of
1396 * course always fails) and that uses the error message associated
1397 * with the caught exception (returned by calling `e.what()`) as
1398 * error text. In other words, instead of letting the error
1399 * message be produced by `main()`, we rather report it here where
1400 * we can abort the program preserving information about where the
1401 * problem happened.
1402 *
1403 * @code
1404 *   template <class MatrixType, class PreconditionerType>
1405 *   class InverseMatrix : public Subscriptor
1406 *   {
1407 *   public:
1408 *   InverseMatrix(const MatrixType &m,
1409 *   const PreconditionerType &preconditioner);
1410 *  
1411 *  
1412 *   template <typename VectorType>
1413 *   void vmult(VectorType &dst, const VectorType &src) const;
1414 *  
1415 *   private:
1417 *   const PreconditionerType &preconditioner;
1418 *   };
1419 *  
1420 *  
1421 *   template <class MatrixType, class PreconditionerType>
1422 *   InverseMatrix<MatrixType, PreconditionerType>::InverseMatrix(
1423 *   const MatrixType &m,
1424 *   const PreconditionerType &preconditioner)
1425 *   : matrix(&m)
1426 *   , preconditioner(preconditioner)
1427 *   {}
1428 *  
1429 *  
1430 *  
1431 *   template <class MatrixType, class PreconditionerType>
1432 *   template <typename VectorType>
1433 *   void InverseMatrix<MatrixType, PreconditionerType>::vmult(
1434 *   VectorType &dst,
1435 *   const VectorType &src) const
1436 *   {
1437 *   SolverControl solver_control(src.size(), 1e-7 * src.l2_norm());
1438 *   SolverCG<VectorType> cg(solver_control);
1439 *  
1440 *   dst = 0;
1441 *  
1442 *   try
1443 *   {
1444 *   cg.solve(*matrix, dst, src, preconditioner);
1445 *   }
1446 *   catch (std::exception &e)
1447 *   {
1448 *   Assert(false, ExcMessage(e.what()));
1449 *   }
1450 *   }
1451 *  
1452 * @endcode
1453 *
1454 *
1455 * <a name="step_31-Schurcomplementpreconditioner"></a>
1456 * <h4>Schur complement preconditioner</h4>
1457 *
1458
1459 *
1460 * This is the implementation of the Schur complement preconditioner as
1461 * described in detail in the introduction. As opposed to @ref step_20 "step-20" and
1462 * @ref step_22 "step-22", we solve the block system all-at-once using GMRES, and use the
1463 * Schur complement of the block structured matrix to build a good
1464 * preconditioner instead.
1465 *
1466
1467 *
1468 * Let's have a look at the ideal preconditioner matrix
1469 * @f$P=\left(\begin{array}{cc} A & 0 \\ B & -S \end{array}\right)@f$
1470 * described in the introduction. If we apply this matrix in the solution
1471 * of a linear system, convergence of an iterative GMRES solver will be
1472 * governed by the matrix @f{eqnarray*} P^{-1}\left(\begin{array}{cc} A &
1473 * B^T \\ B & 0 \end{array}\right) = \left(\begin{array}{cc} I & A^{-1}
1474 * B^T \\ 0 & I \end{array}\right), @f} which indeed is very simple. A
1475 * GMRES solver based on exact matrices would converge in one iteration,
1476 * since all eigenvalues are equal (any Krylov method takes at most as
1477 * many iterations as there are distinct eigenvalues). Such a
1478 * preconditioner for the blocked Stokes system has been proposed by
1479 * Silvester and Wathen ("Fast iterative solution of stabilised Stokes
1480 * systems part II. Using general block preconditioners", SIAM
1481 * J. Numer. Anal., 31 (1994), pp. 1352-1367).
1482 *
1483
1484 *
1485 * Replacing @f$P@f$ by @f$\tilde{P}@f$ keeps that spirit alive: the product
1486 * @f$P^{-1} A@f$ will still be close to a matrix with eigenvalues 1 with a
1487 * distribution that does not depend on the problem size. This lets us
1488 * hope to be able to get a number of GMRES iterations that is
1489 * problem-size independent.
1490 *
1491
1492 *
1493 * The deal.II users who have already gone through the @ref step_20 "step-20" and @ref step_22 "step-22"
1494 * tutorials can certainly imagine how we're going to implement this. We
1495 * replace the exact inverse matrices in @f$P^{-1}@f$ by some approximate
1496 * inverses built from the InverseMatrix class, and the inverse Schur
1497 * complement will be approximated by the pressure mass matrix @f$M_p@f$
1498 * (weighted by @f$\eta^{-1}@f$ as mentioned in the introduction). As pointed
1499 * out in the results section of @ref step_22 "step-22", we can replace the exact inverse
1500 * of @f$A@f$ by just the application of a preconditioner, in this case
1501 * on a vector Laplace matrix as was explained in the introduction. This
1502 * does increase the number of (outer) GMRES iterations, but is still
1503 * significantly cheaper than an exact inverse, which would require
1504 * between 20 and 35 CG iterations for <em>each</em> outer solver step
1505 * (using the AMG preconditioner).
1506 *
1507
1508 *
1509 * Having the above explanations in mind, we define a preconditioner class
1510 * with a <code>vmult</code> functionality, which is all we need for the
1511 * interaction with the usual solver functions further below in the
1512 * program code.
1513 *
1514
1515 *
1516 * First the declarations. These are similar to the definition of the
1517 * Schur complement in @ref step_20 "step-20", with the difference that we need some more
1518 * preconditioners in the constructor and that the matrices we use here
1519 * are built upon Trilinos:
1520 *
1521 * @code
1522 *   template <class PreconditionerTypeA, class PreconditionerTypeMp>
1523 *   class BlockSchurPreconditioner : public Subscriptor
1524 *   {
1525 *   public:
1526 *   BlockSchurPreconditioner(
1527 *   const TrilinosWrappers::BlockSparseMatrix &S,
1528 *   const InverseMatrix<TrilinosWrappers::SparseMatrix,
1529 *   PreconditionerTypeMp> &Mpinv,
1530 *   const PreconditionerTypeA &Apreconditioner);
1531 *  
1532 *   void vmult(TrilinosWrappers::MPI::BlockVector &dst,
1533 *   const TrilinosWrappers::MPI::BlockVector &src) const;
1534 *  
1535 *   private:
1537 *   stokes_matrix;
1538 *   const SmartPointer<const InverseMatrix<TrilinosWrappers::SparseMatrix,
1539 *   PreconditionerTypeMp>>
1540 *   m_inverse;
1541 *   const PreconditionerTypeA &a_preconditioner;
1542 *  
1543 *   mutable TrilinosWrappers::MPI::Vector tmp;
1544 *   };
1545 *  
1546 *  
1547 *  
1548 * @endcode
1549 *
1550 * When using a TrilinosWrappers::MPI::Vector or a
1551 * TrilinosWrappers::MPI::BlockVector, the Vector is initialized using an
1552 * IndexSet. IndexSet is used not only to resize the
1553 * TrilinosWrappers::MPI::Vector but it also associates an index in the
1554 * TrilinosWrappers::MPI::Vector with a degree of freedom (see @ref step_40 "step-40" for
1555 * a more detailed explanation). The function complete_index_set() creates
1556 * an IndexSet where every valid index is part of the set. Note that this
1557 * program can only be run sequentially and will throw an exception if used
1558 * in parallel.
1559 *
1560 * @code
1561 *   template <class PreconditionerTypeA, class PreconditionerTypeMp>
1562 *   BlockSchurPreconditioner<PreconditionerTypeA, PreconditionerTypeMp>::
1563 *   BlockSchurPreconditioner(
1564 *   const TrilinosWrappers::BlockSparseMatrix &S,
1565 *   const InverseMatrix<TrilinosWrappers::SparseMatrix,
1566 *   PreconditionerTypeMp> &Mpinv,
1567 *   const PreconditionerTypeA &Apreconditioner)
1568 *   : stokes_matrix(&S)
1569 *   , m_inverse(&Mpinv)
1570 *   , a_preconditioner(Apreconditioner)
1571 *   , tmp(complete_index_set(stokes_matrix->block(1, 1).m()))
1572 *   {}
1573 *  
1574 *  
1575 * @endcode
1576 *
1577 * Next is the <code>vmult</code> function. We implement the action of
1578 * @f$P^{-1}@f$ as described above in three successive steps. In formulas, we
1579 * want to compute @f$Y=P^{-1}X@f$ where @f$X,Y@f$ are both vectors with two block
1580 * components.
1581 *
1582
1583 *
1584 * The first step multiplies the velocity part of the vector by a
1585 * preconditioner of the matrix @f$A@f$, i.e., we compute @f$Y_0={\tilde
1586 * A}^{-1}X_0@f$. The resulting velocity vector is then multiplied by @f$B@f$
1587 * and subtracted from the pressure, i.e., we want to compute @f$X_1-BY_0@f$.
1588 * This second step only acts on the pressure vector and is accomplished
1589 * by the residual function of our matrix classes, except that the sign is
1590 * wrong. Consequently, we change the sign in the temporary pressure
1591 * vector and finally multiply by the inverse pressure mass matrix to get
1592 * the final pressure vector, completing our work on the Stokes
1593 * preconditioner:
1594 *
1595 * @code
1596 *   template <class PreconditionerTypeA, class PreconditionerTypeMp>
1597 *   void
1598 *   BlockSchurPreconditioner<PreconditionerTypeA, PreconditionerTypeMp>::vmult(
1600 *   const TrilinosWrappers::MPI::BlockVector &src) const
1601 *   {
1602 *   a_preconditioner.vmult(dst.block(0), src.block(0));
1603 *   stokes_matrix->block(1, 0).residual(tmp, dst.block(0), src.block(1));
1604 *   tmp *= -1;
1605 *   m_inverse->vmult(dst.block(1), tmp);
1606 *   }
1607 *   } // namespace LinearSolvers
1608 *  
1609 *  
1610 *  
1611 * @endcode
1612 *
1613 *
1614 * <a name="step_31-ThecodeBoussinesqFlowProblemcodeclasstemplate"></a>
1615 * <h3>The <code>BoussinesqFlowProblem</code> class template</h3>
1616 *
1617
1618 *
1619 * The definition of the class that defines the top-level logic of solving
1620 * the time-dependent Boussinesq problem is mainly based on the @ref step_22 "step-22"
1621 * tutorial program. The main differences are that now we also have to solve
1622 * for the temperature equation, which forces us to have a second DoFHandler
1623 * object for the temperature variable as well as matrices, right hand
1624 * sides, and solution vectors for the current and previous time steps. As
1625 * mentioned in the introduction, all linear algebra objects are going to
1626 * use wrappers of the corresponding Trilinos functionality.
1627 *
1628
1629 *
1630 * The member functions of this class are reminiscent of @ref step_21 "step-21", where we
1631 * also used a staggered scheme that first solve the flow equations (here
1632 * the Stokes equations, in @ref step_21 "step-21" Darcy flow) and then update the advected
1633 * quantity (here the temperature, there the saturation). The functions that
1634 * are new are mainly concerned with determining the time step, as well as
1635 * the proper size of the artificial viscosity stabilization.
1636 *
1637
1638 *
1639 * The last three variables indicate whether the various matrices or
1640 * preconditioners need to be rebuilt the next time the corresponding build
1641 * functions are called. This allows us to move the corresponding
1642 * <code>if</code> into the respective function and thereby keeping our main
1643 * <code>run()</code> function clean and easy to read.
1644 *
1645 * @code
1646 *   template <int dim>
1647 *   class BoussinesqFlowProblem
1648 *   {
1649 *   public:
1650 *   BoussinesqFlowProblem();
1651 *   void run();
1652 *  
1653 *   private:
1654 *   void setup_dofs();
1655 *   void assemble_stokes_preconditioner();
1656 *   void build_stokes_preconditioner();
1657 *   void assemble_stokes_system();
1658 *   void assemble_temperature_system(const double maximal_velocity);
1659 *   void assemble_temperature_matrix();
1660 *   double get_maximal_velocity() const;
1661 *   std::pair<double, double> get_extrapolated_temperature_range() const;
1662 *   void solve();
1663 *   void output_results() const;
1664 *   void refine_mesh(const unsigned int max_grid_level);
1665 *  
1666 *   double compute_viscosity(
1667 *   const std::vector<double> &old_temperature,
1668 *   const std::vector<double> &old_old_temperature,
1669 *   const std::vector<Tensor<1, dim>> &old_temperature_grads,
1670 *   const std::vector<Tensor<1, dim>> &old_old_temperature_grads,
1671 *   const std::vector<double> &old_temperature_laplacians,
1672 *   const std::vector<double> &old_old_temperature_laplacians,
1673 *   const std::vector<Tensor<1, dim>> &old_velocity_values,
1674 *   const std::vector<Tensor<1, dim>> &old_old_velocity_values,
1675 *   const std::vector<double> &gamma_values,
1676 *   const double global_u_infty,
1677 *   const double global_T_variation,
1678 *   const double cell_diameter) const;
1679 *  
1680 *  
1682 *   double global_Omega_diameter;
1683 *  
1684 *   const unsigned int stokes_degree;
1685 *   FESystem<dim> stokes_fe;
1686 *   DoFHandler<dim> stokes_dof_handler;
1687 *   AffineConstraints<double> stokes_constraints;
1688 *  
1689 *   std::vector<IndexSet> stokes_partitioning;
1690 *   TrilinosWrappers::BlockSparseMatrix stokes_matrix;
1691 *   TrilinosWrappers::BlockSparseMatrix stokes_preconditioner_matrix;
1692 *  
1693 *   TrilinosWrappers::MPI::BlockVector stokes_solution;
1694 *   TrilinosWrappers::MPI::BlockVector old_stokes_solution;
1695 *   TrilinosWrappers::MPI::BlockVector stokes_rhs;
1696 *  
1697 *  
1698 *   const unsigned int temperature_degree;
1699 *   FE_Q<dim> temperature_fe;
1700 *   DoFHandler<dim> temperature_dof_handler;
1701 *   AffineConstraints<double> temperature_constraints;
1702 *  
1703 *   TrilinosWrappers::SparseMatrix temperature_mass_matrix;
1704 *   TrilinosWrappers::SparseMatrix temperature_stiffness_matrix;
1705 *   TrilinosWrappers::SparseMatrix temperature_matrix;
1706 *  
1707 *   TrilinosWrappers::MPI::Vector temperature_solution;
1708 *   TrilinosWrappers::MPI::Vector old_temperature_solution;
1709 *   TrilinosWrappers::MPI::Vector old_old_temperature_solution;
1710 *   TrilinosWrappers::MPI::Vector temperature_rhs;
1711 *  
1712 *  
1713 *   double time_step;
1714 *   double old_time_step;
1715 *   unsigned int timestep_number;
1716 *  
1717 *   std::shared_ptr<TrilinosWrappers::PreconditionAMG> Amg_preconditioner;
1718 *   std::shared_ptr<TrilinosWrappers::PreconditionIC> Mp_preconditioner;
1719 *  
1720 *   bool rebuild_stokes_matrix;
1721 *   bool rebuild_temperature_matrices;
1722 *   bool rebuild_stokes_preconditioner;
1723 *   };
1724 *  
1725 *  
1726 * @endcode
1727 *
1728 *
1729 * <a name="step_31-BoussinesqFlowProblemclassimplementation"></a>
1730 * <h3>BoussinesqFlowProblem class implementation</h3>
1731 *
1732
1733 *
1734 *
1735 * <a name="step_31-BoussinesqFlowProblemBoussinesqFlowProblem"></a>
1736 * <h4>BoussinesqFlowProblem::BoussinesqFlowProblem</h4>
1737 *
1738
1739 *
1740 * The constructor of this class is an extension of the constructor in
1741 * @ref step_22 "step-22". We need to add the various variables that concern the
1742 * temperature. As discussed in the introduction, we are going to use
1743 * @f$Q_2^d\times Q_1@f$ (Taylor-Hood) elements again for the Stokes part, and
1744 * @f$Q_2@f$ elements for the temperature. However, by using variables that
1745 * store the polynomial degree of the Stokes and temperature finite
1746 * elements, it is easy to consistently modify the degree of the elements as
1747 * well as all quadrature formulas used on them downstream. Moreover, we
1748 * initialize the time stepping as well as the options for matrix assembly
1749 * and preconditioning:
1750 *
1751 * @code
1752 *   template <int dim>
1753 *   BoussinesqFlowProblem<dim>::BoussinesqFlowProblem()
1755 *   , global_Omega_diameter(std::numeric_limits<double>::quiet_NaN())
1756 *   , stokes_degree(1)
1757 *   , stokes_fe(FE_Q<dim>(stokes_degree + 1) ^ dim, FE_Q<dim>(stokes_degree))
1758 *   , stokes_dof_handler(triangulation)
1759 *   ,
1760 *  
1761 *   temperature_degree(2)
1762 *   , temperature_fe(temperature_degree)
1763 *   , temperature_dof_handler(triangulation)
1764 *   ,
1765 *  
1766 *   time_step(0)
1767 *   , old_time_step(0)
1768 *   , timestep_number(0)
1769 *   , rebuild_stokes_matrix(true)
1770 *   , rebuild_temperature_matrices(true)
1771 *   , rebuild_stokes_preconditioner(true)
1772 *   {}
1773 *  
1774 *  
1775 *  
1776 * @endcode
1777 *
1778 *
1779 * <a name="step_31-BoussinesqFlowProblemget_maximal_velocity"></a>
1780 * <h4>BoussinesqFlowProblem::get_maximal_velocity</h4>
1781 *
1782
1783 *
1784 * Starting the real functionality of this class is a helper function that
1785 * determines the maximum (@f$L_\infty@f$) velocity in the domain (at the
1786 * quadrature points, in fact). How it works should be relatively obvious to
1787 * all who have gotten to this point of the tutorial. Note that since we are
1788 * only interested in the velocity, rather than using
1789 * <code>stokes_fe_values.get_function_values</code> to get the values of
1790 * the entire Stokes solution (velocities and pressures) we use
1791 * <code>stokes_fe_values[velocities].get_function_values</code> to extract
1792 * only the velocities part. This has the additional benefit that we get it
1793 * as a Tensor<1,dim>, rather than some components in a Vector<double>,
1794 * allowing us to process it right away using the <code>norm()</code>
1795 * function to get the magnitude of the velocity.
1796 *
1797
1798 *
1799 * The only point worth thinking about a bit is how to choose the quadrature
1800 * points we use here. Since the goal of this function is to find the
1801 * maximal velocity over a domain by looking at quadrature points on each
1802 * cell. So we should ask how we should best choose these quadrature points
1803 * on each cell. To this end, recall that if we had a single @f$Q_1@f$ field
1804 * (rather than the vector-valued field of higher order) then the maximum
1805 * would be attained at a vertex of the mesh. In other words, we should use
1806 * the QTrapezoid class that has quadrature points only at the vertices of
1807 * cells.
1808 *
1809
1810 *
1811 * For higher order shape functions, the situation is more complicated: the
1812 * maxima and minima may be attained at points between the support points of
1813 * shape functions (for the usual @f$Q_p@f$ elements the support points are the
1814 * equidistant Lagrange interpolation points); furthermore, since we are
1815 * looking for the maximum magnitude of a vector-valued quantity, we can
1816 * even less say with certainty where the set of potential maximal points
1817 * are. Nevertheless, intuitively if not provably, the Lagrange
1818 * interpolation points appear to be a better choice than the Gauss points.
1819 *
1820
1821 *
1822 * There are now different methods to produce a quadrature formula with
1823 * quadrature points equal to the interpolation points of the finite
1824 * element. One option would be to use the
1825 * FiniteElement::get_unit_support_points() function, reduce the output to a
1826 * unique set of points to avoid duplicate function evaluations, and create
1827 * a Quadrature object using these points. Another option, chosen here, is
1828 * to use the QTrapezoid class and combine it with the QIterated class that
1829 * repeats the QTrapezoid formula on a number of sub-cells in each coordinate
1830 * direction. To cover all support points, we need to iterate it
1831 * <code>stokes_degree+1</code> times since this is the polynomial degree of
1832 * the Stokes element in use:
1833 *
1834 * @code
1835 *   template <int dim>
1836 *   double BoussinesqFlowProblem<dim>::get_maximal_velocity() const
1837 *   {
1838 *   const QIterated<dim> quadrature_formula(QTrapezoid<1>(), stokes_degree + 1);
1839 *   const unsigned int n_q_points = quadrature_formula.size();
1840 *  
1841 *   FEValues<dim> fe_values(stokes_fe, quadrature_formula, update_values);
1842 *   std::vector<Tensor<1, dim>> velocity_values(n_q_points);
1843 *   double max_velocity = 0;
1844 *  
1845 *   const FEValuesExtractors::Vector velocities(0);
1846 *  
1847 *   for (const auto &cell : stokes_dof_handler.active_cell_iterators())
1848 *   {
1849 *   fe_values.reinit(cell);
1850 *   fe_values[velocities].get_function_values(stokes_solution,
1851 *   velocity_values);
1852 *  
1853 *   for (unsigned int q = 0; q < n_q_points; ++q)
1854 *   max_velocity = std::max(max_velocity, velocity_values[q].norm());
1855 *   }
1856 *  
1857 *   return max_velocity;
1858 *   }
1859 *  
1860 *  
1861 *  
1862 * @endcode
1863 *
1864 *
1865 * <a name="step_31-BoussinesqFlowProblemget_extrapolated_temperature_range"></a>
1866 * <h4>BoussinesqFlowProblem::get_extrapolated_temperature_range</h4>
1867 *
1868
1869 *
1870 * Next a function that determines the minimum and maximum temperature at
1871 * quadrature points inside @f$\Omega@f$ when extrapolated from the two previous
1872 * time steps to the current one. We need this information in the
1873 * computation of the artificial viscosity parameter @f$\nu@f$ as discussed in
1874 * the introduction.
1875 *
1876
1877 *
1878 * The formula for the extrapolated temperature is
1879 * @f$\left(1+\frac{k_n}{k_{n-1}} \right)T^{n-1} + \frac{k_n}{k_{n-1}}
1880 * T^{n-2}@f$. The way to compute it is to loop over all quadrature points and
1881 * update the maximum and minimum value if the current value is
1882 * bigger/smaller than the previous one. We initialize the variables that
1883 * store the max and min before the loop over all quadrature points by the
1884 * smallest and the largest number representable as a double. Then we know
1885 * for a fact that it is larger/smaller than the minimum/maximum and that
1886 * the loop over all quadrature points is ultimately going to update the
1887 * initial value with the correct one.
1888 *
1889
1890 *
1891 * The only other complication worth mentioning here is that in the first
1892 * time step, @f$T^{k-2}@f$ is not yet available of course. In that case, we can
1893 * only use @f$T^{k-1}@f$ which we have from the initial temperature. As
1894 * quadrature points, we use the same choice as in the previous function
1895 * though with the difference that now the number of repetitions is
1896 * determined by the polynomial degree of the temperature field.
1897 *
1898 * @code
1899 *   template <int dim>
1900 *   std::pair<double, double>
1901 *   BoussinesqFlowProblem<dim>::get_extrapolated_temperature_range() const
1902 *   {
1903 *   const QIterated<dim> quadrature_formula(QTrapezoid<1>(),
1904 *   temperature_degree);
1905 *   const unsigned int n_q_points = quadrature_formula.size();
1906 *  
1907 *   FEValues<dim> fe_values(temperature_fe, quadrature_formula, update_values);
1908 *   std::vector<double> old_temperature_values(n_q_points);
1909 *   std::vector<double> old_old_temperature_values(n_q_points);
1910 *  
1911 *   if (timestep_number != 0)
1912 *   {
1913 *   double min_temperature = std::numeric_limits<double>::max(),
1914 *   max_temperature = -std::numeric_limits<double>::max();
1915 *  
1916 *   for (const auto &cell : temperature_dof_handler.active_cell_iterators())
1917 *   {
1918 *   fe_values.reinit(cell);
1919 *   fe_values.get_function_values(old_temperature_solution,
1920 *   old_temperature_values);
1921 *   fe_values.get_function_values(old_old_temperature_solution,
1922 *   old_old_temperature_values);
1923 *  
1924 *   for (unsigned int q = 0; q < n_q_points; ++q)
1925 *   {
1926 *   const double temperature =
1927 *   (1. + time_step / old_time_step) * old_temperature_values[q] -
1928 *   time_step / old_time_step * old_old_temperature_values[q];
1929 *  
1930 *   min_temperature = std::min(min_temperature, temperature);
1931 *   max_temperature = std::max(max_temperature, temperature);
1932 *   }
1933 *   }
1934 *  
1935 *   return std::make_pair(min_temperature, max_temperature);
1936 *   }
1937 *   else
1938 *   {
1939 *   double min_temperature = std::numeric_limits<double>::max(),
1940 *   max_temperature = -std::numeric_limits<double>::max();
1941 *  
1942 *   for (const auto &cell : temperature_dof_handler.active_cell_iterators())
1943 *   {
1944 *   fe_values.reinit(cell);
1945 *   fe_values.get_function_values(old_temperature_solution,
1946 *   old_temperature_values);
1947 *  
1948 *   for (unsigned int q = 0; q < n_q_points; ++q)
1949 *   {
1950 *   const double temperature = old_temperature_values[q];
1951 *  
1952 *   min_temperature = std::min(min_temperature, temperature);
1953 *   max_temperature = std::max(max_temperature, temperature);
1954 *   }
1955 *   }
1956 *  
1957 *   return std::make_pair(min_temperature, max_temperature);
1958 *   }
1959 *   }
1960 *  
1961 *  
1962 *  
1963 * @endcode
1964 *
1965 *
1966 * <a name="step_31-BoussinesqFlowProblemcompute_viscosity"></a>
1967 * <h4>BoussinesqFlowProblem::compute_viscosity</h4>
1968 *
1969
1970 *
1971 * The last of the tool functions computes the artificial viscosity
1972 * parameter @f$\nu|_K@f$ on a cell @f$K@f$ as a function of the extrapolated
1973 * temperature, its gradient and Hessian (second derivatives), the velocity,
1974 * the right hand side @f$\gamma@f$ all on the quadrature points of the current
1975 * cell, and various other parameters as described in detail in the
1976 * introduction.
1977 *
1978
1979 *
1980 * There are some universal constants worth mentioning here. First, we need
1981 * to fix @f$\beta@f$; we choose @f$\beta=0.017\cdot dim@f$, a choice discussed in
1982 * detail in the results section of this tutorial program. The second is the
1983 * exponent @f$\alpha@f$; @f$\alpha=1@f$ appears to work fine for the current
1984 * program, even though some additional benefit might be expected from
1985 * choosing @f$\alpha = 2@f$. Finally, there is one thing that requires special
1986 * casing: In the first time step, the velocity equals zero, and the formula
1987 * for @f$\nu|_K@f$ is not defined. In that case, we return @f$\nu|_K=5\cdot 10^3
1988 * \cdot h_K@f$, a choice admittedly more motivated by heuristics than
1989 * anything else (it is in the same order of magnitude, however, as the
1990 * value returned for most cells on the second time step).
1991 *
1992
1993 *
1994 * The rest of the function should be mostly obvious based on the material
1995 * discussed in the introduction:
1996 *
1997 * @code
1998 *   template <int dim>
1999 *   double BoussinesqFlowProblem<dim>::compute_viscosity(
2000 *   const std::vector<double> &old_temperature,
2001 *   const std::vector<double> &old_old_temperature,
2002 *   const std::vector<Tensor<1, dim>> &old_temperature_grads,
2003 *   const std::vector<Tensor<1, dim>> &old_old_temperature_grads,
2004 *   const std::vector<double> &old_temperature_laplacians,
2005 *   const std::vector<double> &old_old_temperature_laplacians,
2006 *   const std::vector<Tensor<1, dim>> &old_velocity_values,
2007 *   const std::vector<Tensor<1, dim>> &old_old_velocity_values,
2008 *   const std::vector<double> &gamma_values,
2009 *   const double global_u_infty,
2010 *   const double global_T_variation,
2011 *   const double cell_diameter) const
2012 *   {
2013 *   constexpr double beta = 0.017 * dim;
2014 *   constexpr double alpha = 1.0;
2015 *  
2016 *   if (global_u_infty == 0)
2017 *   return 5e-3 * cell_diameter;
2018 *  
2019 *   const unsigned int n_q_points = old_temperature.size();
2020 *  
2021 *   double max_residual = 0;
2022 *   double max_velocity = 0;
2023 *  
2024 *   for (unsigned int q = 0; q < n_q_points; ++q)
2025 *   {
2026 *   const Tensor<1, dim> u =
2027 *   (old_velocity_values[q] + old_old_velocity_values[q]) / 2;
2028 *  
2029 *   const double dT_dt =
2030 *   (old_temperature[q] - old_old_temperature[q]) / old_time_step;
2031 *   const double u_grad_T =
2032 *   u * (old_temperature_grads[q] + old_old_temperature_grads[q]) / 2;
2033 *  
2034 *   const double kappa_Delta_T =
2035 *   EquationData::kappa *
2036 *   (old_temperature_laplacians[q] + old_old_temperature_laplacians[q]) /
2037 *   2;
2038 *  
2039 *   const double residual =
2040 *   std::abs((dT_dt + u_grad_T - kappa_Delta_T - gamma_values[q]) *
2041 *   std::pow((old_temperature[q] + old_old_temperature[q]) / 2,
2042 *   alpha - 1.));
2043 *  
2044 *   max_residual = std::max(residual, max_residual);
2045 *   max_velocity = std::max(std::sqrt(u * u), max_velocity);
2046 *   }
2047 *  
2048 *   const double c_R = std::pow(2., (4. - 2 * alpha) / dim);
2049 *   const double global_scaling = c_R * global_u_infty * global_T_variation *
2050 *   std::pow(global_Omega_diameter, alpha - 2.);
2051 *  
2052 *   return (
2053 *   beta * max_velocity *
2054 *   std::min(cell_diameter,
2055 *   std::pow(cell_diameter, alpha) * max_residual / global_scaling));
2056 *   }
2057 *  
2058 *  
2059 *  
2060 * @endcode
2061 *
2062 *
2063 * <a name="step_31-BoussinesqFlowProblemsetup_dofs"></a>
2064 * <h4>BoussinesqFlowProblem::setup_dofs</h4>
2065 *
2066
2067 *
2068 * This is the function that sets up the DoFHandler objects we have here
2069 * (one for the Stokes part and one for the temperature part) as well as set
2070 * to the right sizes the various objects required for the linear algebra in
2071 * this program. Its basic operations are similar to what we do in @ref step_22 "step-22".
2072 *
2073
2074 *
2075 * The body of the function first enumerates all degrees of freedom for the
2076 * Stokes and temperature systems. For the Stokes part, degrees of freedom
2077 * are then sorted to ensure that velocities precede pressure DoFs so that
2078 * we can partition the Stokes matrix into a @f$2\times 2@f$ matrix. As a
2079 * difference to @ref step_22 "step-22", we do not perform any additional DoF
2080 * renumbering. In that program, it paid off since our solver was heavily
2081 * dependent on ILU's, whereas we use AMG here which is not sensitive to the
2082 * DoF numbering. The IC preconditioner for the inversion of the pressure
2083 * mass matrix would of course take advantage of a Cuthill-McKee like
2084 * renumbering, but its costs are low compared to the velocity portion, so
2085 * the additional work does not pay off.
2086 *
2087
2088 *
2089 * We then proceed with the generation of the hanging node constraints that
2090 * arise from adaptive grid refinement for both DoFHandler objects. For the
2091 * velocity, we impose no-flux boundary conditions @f$\mathbf{u}\cdot
2092 * \mathbf{n}=0@f$ by adding constraints to the object that already stores the
2093 * hanging node constraints matrix. The second parameter in the function
2094 * describes the first of the velocity components in the total dof vector,
2095 * which is zero here. The variable <code>no_normal_flux_boundaries</code>
2096 * denotes the boundary indicators for which to set the no flux boundary
2097 * conditions; here, this is boundary indicator zero.
2098 *
2099
2100 *
2101 * After having done so, we count the number of degrees of freedom in the
2102 * various blocks:
2103 *
2104 * @code
2105 *   template <int dim>
2106 *   void BoussinesqFlowProblem<dim>::setup_dofs()
2107 *   {
2108 *   std::vector<unsigned int> stokes_sub_blocks(dim + 1, 0);
2109 *   stokes_sub_blocks[dim] = 1;
2110 *  
2111 *   {
2112 *   stokes_dof_handler.distribute_dofs(stokes_fe);
2113 *   DoFRenumbering::component_wise(stokes_dof_handler, stokes_sub_blocks);
2114 *  
2115 *   stokes_constraints.clear();
2116 *   DoFTools::make_hanging_node_constraints(stokes_dof_handler,
2117 *   stokes_constraints);
2118 *   const std::set<types::boundary_id> no_normal_flux_boundaries = {0};
2119 *   VectorTools::compute_no_normal_flux_constraints(stokes_dof_handler,
2120 *   0,
2121 *   no_normal_flux_boundaries,
2122 *   stokes_constraints);
2123 *   stokes_constraints.close();
2124 *   }
2125 *   {
2126 *   temperature_dof_handler.distribute_dofs(temperature_fe);
2127 *  
2128 *   temperature_constraints.clear();
2129 *   DoFTools::make_hanging_node_constraints(temperature_dof_handler,
2130 *   temperature_constraints);
2131 *   temperature_constraints.close();
2132 *   }
2133 *  
2134 *   const std::vector<types::global_dof_index> stokes_dofs_per_block =
2135 *   DoFTools::count_dofs_per_fe_block(stokes_dof_handler, stokes_sub_blocks);
2136 *  
2137 *   const types::global_dof_index n_u = stokes_dofs_per_block[0],
2138 *   n_p = stokes_dofs_per_block[1],
2139 *   n_T = temperature_dof_handler.n_dofs();
2140 *  
2141 *   std::cout << "Number of active cells: " << triangulation.n_active_cells()
2142 *   << " (on " << triangulation.n_levels() << " levels)" << std::endl
2143 *   << "Number of degrees of freedom: " << n_u + n_p + n_T << " ("
2144 *   << n_u << '+' << n_p << '+' << n_T << ')' << std::endl
2145 *   << std::endl;
2146 *  
2147 * @endcode
2148 *
2149 * The next step is to create the sparsity pattern for the Stokes and
2150 * temperature system matrices as well as the preconditioner matrix from
2151 * which we build the Stokes preconditioner. As in @ref step_22 "step-22", we choose to
2152 * create the pattern by
2153 * using the blocked version of DynamicSparsityPattern.
2154 *
2155
2156 *
2157 * So, we first release the memory stored in the matrices, then set up an
2158 * object of type BlockDynamicSparsityPattern consisting of
2159 * @f$2\times 2@f$ blocks (for the Stokes system matrix and preconditioner) or
2160 * DynamicSparsityPattern (for the temperature part). We then
2161 * fill these objects with the nonzero pattern, taking into account that
2162 * for the Stokes system matrix, there are no entries in the
2163 * pressure-pressure block (but all velocity vector components couple with
2164 * each other and with the pressure). Similarly, in the Stokes
2165 * preconditioner matrix, only the diagonal blocks are nonzero, since we
2166 * use the vector Laplacian as discussed in the introduction. This
2167 * operator only couples each vector component of the Laplacian with
2168 * itself, but not with the other vector components. (Application of the
2169 * constraints resulting from the no-flux boundary conditions will couple
2170 * vector components at the boundary again, however.)
2171 *
2172
2173 *
2174 * When generating the sparsity pattern, we directly apply the constraints
2175 * from hanging nodes and no-flux boundary conditions. This approach was
2176 * already used in @ref step_27 "step-27", but is different from the one in early
2177 * tutorial programs where we first built the original sparsity pattern
2178 * and only then added the entries resulting from constraints. The reason
2179 * for doing so is that later during assembly we are going to distribute
2180 * the constraints immediately when transferring local to global
2181 * dofs. Consequently, there will be no data written at positions of
2182 * constrained degrees of freedom, so we can let the
2183 * DoFTools::make_sparsity_pattern function omit these entries by setting
2184 * the last Boolean flag to <code>false</code>. Once the sparsity pattern
2185 * is ready, we can use it to initialize the Trilinos matrices. Since the
2186 * Trilinos matrices store the sparsity pattern internally, there is no
2187 * need to keep the sparsity pattern around after the initialization of
2188 * the matrix.
2189 *
2190 * @code
2191 *   stokes_partitioning.resize(2);
2192 *   stokes_partitioning[0] = complete_index_set(n_u);
2193 *   stokes_partitioning[1] = complete_index_set(n_p);
2194 *   {
2195 *   stokes_matrix.clear();
2196 *  
2197 *   BlockDynamicSparsityPattern dsp(stokes_dofs_per_block,
2198 *   stokes_dofs_per_block);
2199 *  
2200 *   Table<2, DoFTools::Coupling> coupling(dim + 1, dim + 1);
2201 *  
2202 *   for (unsigned int c = 0; c < dim + 1; ++c)
2203 *   for (unsigned int d = 0; d < dim + 1; ++d)
2204 *   if (!((c == dim) && (d == dim)))
2205 *   coupling[c][d] = DoFTools::always;
2206 *   else
2207 *   coupling[c][d] = DoFTools::none;
2208 *  
2209 *   DoFTools::make_sparsity_pattern(
2210 *   stokes_dof_handler, coupling, dsp, stokes_constraints, false);
2211 *  
2212 *   stokes_matrix.reinit(dsp);
2213 *   }
2214 *  
2215 *   {
2216 *   Amg_preconditioner.reset();
2217 *   Mp_preconditioner.reset();
2218 *   stokes_preconditioner_matrix.clear();
2219 *  
2220 *   BlockDynamicSparsityPattern dsp(stokes_dofs_per_block,
2221 *   stokes_dofs_per_block);
2222 *  
2223 *   Table<2, DoFTools::Coupling> coupling(dim + 1, dim + 1);
2224 *   for (unsigned int c = 0; c < dim + 1; ++c)
2225 *   for (unsigned int d = 0; d < dim + 1; ++d)
2226 *   if (c == d)
2227 *   coupling[c][d] = DoFTools::always;
2228 *   else
2229 *   coupling[c][d] = DoFTools::none;
2230 *  
2231 *   DoFTools::make_sparsity_pattern(
2232 *   stokes_dof_handler, coupling, dsp, stokes_constraints, false);
2233 *  
2234 *   stokes_preconditioner_matrix.reinit(dsp);
2235 *   }
2236 *  
2237 * @endcode
2238 *
2239 * The creation of the temperature matrix (or, rather, matrices, since we
2240 * provide a temperature mass matrix and a temperature @ref GlossStiffnessMatrix "stiffness matrix",
2241 * that will be added together for time discretization) follows the
2242 * generation of the Stokes matrix &ndash; except that it is much easier
2243 * here since we do not need to take care of any blocks or coupling
2244 * between components. Note how we initialize the three temperature
2245 * matrices: We only use the sparsity pattern for reinitialization of the
2246 * first matrix, whereas we use the previously generated matrix for the
2247 * two remaining reinits. The reason for doing so is that reinitialization
2248 * from an already generated matrix allows Trilinos to reuse the sparsity
2249 * pattern instead of generating a new one for each copy. This saves both
2250 * some time and memory.
2251 *
2252 * @code
2253 *   {
2254 *   temperature_mass_matrix.clear();
2255 *   temperature_stiffness_matrix.clear();
2256 *   temperature_matrix.clear();
2257 *  
2258 *   DynamicSparsityPattern dsp(n_T, n_T);
2259 *   DoFTools::make_sparsity_pattern(temperature_dof_handler,
2260 *   dsp,
2261 *   temperature_constraints,
2262 *   false);
2263 *  
2264 *   temperature_matrix.reinit(dsp);
2265 *   temperature_mass_matrix.reinit(temperature_matrix);
2266 *   temperature_stiffness_matrix.reinit(temperature_matrix);
2267 *   }
2268 *  
2269 * @endcode
2270 *
2271 * Lastly, we set the vectors for the Stokes solutions @f$\mathbf u^{n-1}@f$
2272 * and @f$\mathbf u^{n-2}@f$, as well as for the temperatures @f$T^{n}@f$,
2273 * @f$T^{n-1}@f$ and @f$T^{n-2}@f$ (required for time stepping) and all the system
2274 * right hand sides to their correct sizes and block structure:
2275 *
2276 * @code
2277 *   IndexSet temperature_partitioning = complete_index_set(n_T);
2278 *   stokes_solution.reinit(stokes_partitioning, MPI_COMM_WORLD);
2279 *   old_stokes_solution.reinit(stokes_partitioning, MPI_COMM_WORLD);
2280 *   stokes_rhs.reinit(stokes_partitioning, MPI_COMM_WORLD);
2281 *  
2282 *   temperature_solution.reinit(temperature_partitioning, MPI_COMM_WORLD);
2283 *   old_temperature_solution.reinit(temperature_partitioning, MPI_COMM_WORLD);
2284 *   old_old_temperature_solution.reinit(temperature_partitioning,
2285 *   MPI_COMM_WORLD);
2286 *  
2287 *   temperature_rhs.reinit(temperature_partitioning, MPI_COMM_WORLD);
2288 *   }
2289 *  
2290 *  
2291 *  
2292 * @endcode
2293 *
2294 *
2295 * <a name="step_31-BoussinesqFlowProblemassemble_stokes_preconditioner"></a>
2296 * <h4>BoussinesqFlowProblem::assemble_stokes_preconditioner</h4>
2297 *
2298
2299 *
2300 * This function assembles the matrix we use for preconditioning the Stokes
2301 * system. What we need are a vector Laplace matrix on the velocity
2302 * components and a mass matrix weighted by @f$\eta^{-1}@f$ on the pressure
2303 * component. We start by generating a quadrature object of appropriate
2304 * order, the FEValues object that can give values and gradients at the
2305 * quadrature points (together with quadrature weights). Next we create data
2306 * structures for the cell matrix and the relation between local and global
2307 * DoFs. The vectors <code>grad_phi_u</code> and <code>phi_p</code> are
2308 * going to hold the values of the basis functions in order to faster build
2309 * up the local matrices, as was already done in @ref step_22 "step-22". Before we start
2310 * the loop over all active cells, we have to specify which components are
2311 * pressure and which are velocity.
2312 *
2313 * @code
2314 *   template <int dim>
2315 *   void BoussinesqFlowProblem<dim>::assemble_stokes_preconditioner()
2316 *   {
2317 *   stokes_preconditioner_matrix = 0;
2318 *  
2319 *   const QGauss<dim> quadrature_formula(stokes_degree + 2);
2320 *   FEValues<dim> stokes_fe_values(stokes_fe,
2321 *   quadrature_formula,
2322 *   update_JxW_values | update_values |
2323 *   update_gradients);
2324 *  
2325 *   const unsigned int dofs_per_cell = stokes_fe.n_dofs_per_cell();
2326 *   const unsigned int n_q_points = quadrature_formula.size();
2327 *  
2328 *   FullMatrix<double> local_matrix(dofs_per_cell, dofs_per_cell);
2329 *   std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
2330 *  
2331 *   std::vector<Tensor<2, dim>> grad_phi_u(dofs_per_cell);
2332 *   std::vector<double> phi_p(dofs_per_cell);
2333 *  
2334 *   const FEValuesExtractors::Vector velocities(0);
2335 *   const FEValuesExtractors::Scalar pressure(dim);
2336 *  
2337 *   for (const auto &cell : stokes_dof_handler.active_cell_iterators())
2338 *   {
2339 *   stokes_fe_values.reinit(cell);
2340 *   local_matrix = 0;
2341 *  
2342 * @endcode
2343 *
2344 * The creation of the local matrix is rather simple. There are only a
2345 * Laplace term (on the velocity) and a mass matrix weighted by
2346 * @f$\eta^{-1}@f$ to be generated, so the creation of the local matrix is
2347 * done in two lines. Once the local matrix is ready (loop over rows
2348 * and columns in the local matrix on each quadrature point), we get
2349 * the local DoF indices and write the local information into the
2350 * global matrix. We do this as in @ref step_27 "step-27", i.e., we directly apply the
2351 * constraints from hanging nodes locally. By doing so, we don't have
2352 * to do that afterwards, and we don't also write into entries of the
2353 * matrix that will actually be set to zero again later when
2354 * eliminating constraints.
2355 *
2356 * @code
2357 *   for (unsigned int q = 0; q < n_q_points; ++q)
2358 *   {
2359 *   for (unsigned int k = 0; k < dofs_per_cell; ++k)
2360 *   {
2361 *   grad_phi_u[k] = stokes_fe_values[velocities].gradient(k, q);
2362 *   phi_p[k] = stokes_fe_values[pressure].value(k, q);
2363 *   }
2364 *  
2365 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
2366 *   for (unsigned int j = 0; j < dofs_per_cell; ++j)
2367 *   local_matrix(i, j) +=
2368 *   (EquationData::eta *
2369 *   scalar_product(grad_phi_u[i], grad_phi_u[j]) +
2370 *   (1. / EquationData::eta) * phi_p[i] * phi_p[j]) *
2371 *   stokes_fe_values.JxW(q);
2372 *   }
2373 *  
2374 *   cell->get_dof_indices(local_dof_indices);
2375 *   stokes_constraints.distribute_local_to_global(
2376 *   local_matrix, local_dof_indices, stokes_preconditioner_matrix);
2377 *   }
2378 *   }
2379 *  
2380 *  
2381 *  
2382 * @endcode
2383 *
2384 *
2385 * <a name="step_31-BoussinesqFlowProblembuild_stokes_preconditioner"></a>
2386 * <h4>BoussinesqFlowProblem::build_stokes_preconditioner</h4>
2387 *
2388
2389 *
2390 * This function generates the inner preconditioners that are going to be
2391 * used for the Schur complement block preconditioner. Since the
2392 * preconditioners need only to be regenerated when the matrices change,
2393 * this function does not have to do anything in case the matrices have not
2394 * changed (i.e., the flag <code>rebuild_stokes_preconditioner</code> has
2395 * the value <code>false</code>). Otherwise its first task is to call
2396 * <code>assemble_stokes_preconditioner</code> to generate the
2397 * preconditioner matrices.
2398 *
2399
2400 *
2401 * Next, we set up the preconditioner for the velocity-velocity matrix
2402 * @f$A@f$. As explained in the introduction, we are going to use an AMG
2403 * preconditioner based on a vector Laplace matrix @f$\hat{A}@f$ (which is
2404 * spectrally close to the Stokes matrix @f$A@f$). Usually, the
2405 * TrilinosWrappers::PreconditionAMG class can be seen as a good black-box
2406 * preconditioner which does not need any special knowledge. In this case,
2407 * however, we have to be careful: since we build an AMG for a vector
2408 * problem, we have to tell the preconditioner setup which dofs belong to
2409 * which vector component. We do this using the function
2410 * DoFTools::extract_constant_modes, a function that generates a set of
2411 * <code>dim</code> vectors, where each one has ones in the respective
2412 * component of the vector problem and zeros elsewhere. Hence, these are the
2413 * constant modes on each component, which explains the name of the
2414 * variable.
2415 *
2416 * @code
2417 *   template <int dim>
2418 *   void BoussinesqFlowProblem<dim>::build_stokes_preconditioner()
2419 *   {
2420 *   if (rebuild_stokes_preconditioner == false)
2421 *   return;
2422 *  
2423 *   std::cout << " Rebuilding Stokes preconditioner..." << std::flush;
2424 *  
2425 *   assemble_stokes_preconditioner();
2426 *  
2427 *   Amg_preconditioner = std::make_shared<TrilinosWrappers::PreconditionAMG>();
2428 *  
2429 *   std::vector<std::vector<bool>> constant_modes;
2430 *   const FEValuesExtractors::Vector velocity_components(0);
2431 *   DoFTools::extract_constant_modes(stokes_dof_handler,
2432 *   stokes_fe.component_mask(
2433 *   velocity_components),
2434 *   constant_modes);
2435 *   TrilinosWrappers::PreconditionAMG::AdditionalData amg_data;
2436 *   amg_data.constant_modes = constant_modes;
2437 *  
2438 * @endcode
2439 *
2440 * Next, we set some more options of the AMG preconditioner. In
2441 * particular, we need to tell the AMG setup that we use quadratic basis
2442 * functions for the velocity matrix (this implies more nonzero elements
2443 * in the matrix, so that a more robust algorithm needs to be chosen
2444 * internally). Moreover, we want to be able to control how the coarsening
2445 * structure is build up. The way the Trilinos smoothed aggregation AMG
2446 * does this is to look which matrix entries are of similar size as the
2447 * diagonal entry in order to algebraically build a coarse-grid
2448 * structure. By setting the parameter <code>aggregation_threshold</code>
2449 * to 0.02, we specify that all entries that are more than two percent of
2450 * size of some diagonal pivots in that row should form one coarse grid
2451 * point. This parameter is rather ad hoc, and some fine-tuning of it can
2452 * influence the performance of the preconditioner. As a rule of thumb,
2453 * larger values of <code>aggregation_threshold</code> will decrease the
2454 * number of iterations, but increase the costs per iteration. A look at
2455 * the Trilinos documentation will provide more information on these
2456 * parameters. With this data set, we then initialize the preconditioner
2457 * with the matrix we want it to apply to.
2458 *
2459
2460 *
2461 * Finally, we also initialize the preconditioner for the inversion of the
2462 * pressure mass matrix. This matrix is symmetric and well-behaved, so we
2463 * can chose a simple preconditioner. We stick with an incomplete Cholesky
2464 * (IC) factorization preconditioner, which is designed for symmetric
2465 * matrices. We could have also chosen an SSOR preconditioner with
2466 * relaxation factor around 1.2, but IC is cheaper for our example. We
2467 * wrap the preconditioners into a <code>std::shared_ptr</code>
2468 * pointer, which makes it easier to recreate the preconditioner next time
2469 * around since we do not have to care about destroying the previously
2470 * used object.
2471 *
2472 * @code
2473 *   amg_data.elliptic = true;
2474 *   amg_data.higher_order_elements = true;
2475 *   amg_data.smoother_sweeps = 2;
2476 *   amg_data.aggregation_threshold = 0.02;
2477 *   Amg_preconditioner->initialize(stokes_preconditioner_matrix.block(0, 0),
2478 *   amg_data);
2479 *  
2480 *   Mp_preconditioner = std::make_shared<TrilinosWrappers::PreconditionIC>();
2481 *   Mp_preconditioner->initialize(stokes_preconditioner_matrix.block(1, 1));
2482 *  
2483 *   std::cout << std::endl;
2484 *  
2485 *   rebuild_stokes_preconditioner = false;
2486 *   }
2487 *  
2488 *  
2489 *  
2490 * @endcode
2491 *
2492 *
2493 * <a name="step_31-BoussinesqFlowProblemassemble_stokes_system"></a>
2494 * <h4>BoussinesqFlowProblem::assemble_stokes_system</h4>
2495 *
2496
2497 *
2498 * The time lag scheme we use for advancing the coupled Stokes-temperature
2499 * system forces us to split up the assembly (and the solution of linear
2500 * systems) into two step. The first one is to create the Stokes system
2501 * matrix and right hand side, and the second is to create matrix and right
2502 * hand sides for the temperature dofs, which depends on the result of the
2503 * linear system for the velocity.
2504 *
2505
2506 *
2507 * This function is called at the beginning of each time step. In the first
2508 * time step or if the mesh has changed, indicated by the
2509 * <code>rebuild_stokes_matrix</code>, we need to assemble the Stokes
2510 * matrix; on the other hand, if the mesh hasn't changed and the matrix is
2511 * already available, this is not necessary and all we need to do is
2512 * assemble the right hand side vector which changes in each time step.
2513 *
2514
2515 *
2516 * Regarding the technical details of implementation, not much has changed
2517 * from @ref step_22 "step-22". We reset matrix and vector, create a quadrature formula on
2518 * the cells, and then create the respective FEValues object. For the update
2519 * flags, we require basis function derivatives only in case of a full
2520 * assembly, since they are not needed for the right hand side; as always,
2521 * choosing the minimal set of flags depending on what is currently needed
2522 * makes the call to FEValues::reinit further down in the program more
2523 * efficient.
2524 *
2525
2526 *
2527 * There is one thing that needs to be commented &ndash; since we have a
2528 * separate finite element and DoFHandler for the temperature, we need to
2529 * generate a second FEValues object for the proper evaluation of the
2530 * temperature solution. This isn't too complicated to realize here: just
2531 * use the temperature structures and set an update flag for the basis
2532 * function values which we need for evaluation of the temperature
2533 * solution. The only important part to remember here is that the same
2534 * quadrature formula is used for both FEValues objects to ensure that we
2535 * get matching information when we loop over the quadrature points of the
2536 * two objects.
2537 *
2538
2539 *
2540 * The declarations proceed with some shortcuts for array sizes, the
2541 * creation of the local matrix and right hand side as well as the vector
2542 * for the indices of the local dofs compared to the global system.
2543 *
2544 * @code
2545 *   template <int dim>
2546 *   void BoussinesqFlowProblem<dim>::assemble_stokes_system()
2547 *   {
2548 *   std::cout << " Assembling..." << std::flush;
2549 *  
2550 *   if (rebuild_stokes_matrix == true)
2551 *   stokes_matrix = 0;
2552 *  
2553 *   stokes_rhs = 0;
2554 *  
2555 *   const QGauss<dim> quadrature_formula(stokes_degree + 2);
2556 *   FEValues<dim> stokes_fe_values(
2557 *   stokes_fe,
2558 *   quadrature_formula,
2559 *   update_values | update_quadrature_points | update_JxW_values |
2560 *   (rebuild_stokes_matrix == true ? update_gradients : UpdateFlags(0)));
2561 *  
2562 *   FEValues<dim> temperature_fe_values(temperature_fe,
2563 *   quadrature_formula,
2564 *   update_values);
2565 *  
2566 *   const unsigned int dofs_per_cell = stokes_fe.n_dofs_per_cell();
2567 *   const unsigned int n_q_points = quadrature_formula.size();
2568 *  
2569 *   FullMatrix<double> local_matrix(dofs_per_cell, dofs_per_cell);
2570 *   Vector<double> local_rhs(dofs_per_cell);
2571 *  
2572 *   std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
2573 *  
2574 * @endcode
2575 *
2576 * Next we need a vector that will contain the values of the temperature
2577 * solution at the previous time level at the quadrature points to
2578 * assemble the source term in the right hand side of the momentum
2579 * equation. Let's call this vector <code>old_solution_values</code>.
2580 *
2581
2582 *
2583 * The set of vectors we create next hold the evaluations of the basis
2584 * functions as well as their gradients and symmetrized gradients that
2585 * will be used for creating the matrices. Putting these into their own
2586 * arrays rather than asking the FEValues object for this information each
2587 * time it is needed is an optimization to accelerate the assembly
2588 * process, see @ref step_22 "step-22" for details.
2589 *
2590
2591 *
2592 * The last two declarations are used to extract the individual blocks
2593 * (velocity, pressure, temperature) from the total FE system.
2594 *
2595 * @code
2596 *   std::vector<double> old_temperature_values(n_q_points);
2597 *  
2598 *   std::vector<Tensor<1, dim>> phi_u(dofs_per_cell);
2599 *   std::vector<SymmetricTensor<2, dim>> grads_phi_u(dofs_per_cell);
2600 *   std::vector<double> div_phi_u(dofs_per_cell);
2601 *   std::vector<double> phi_p(dofs_per_cell);
2602 *  
2603 *   const FEValuesExtractors::Vector velocities(0);
2604 *   const FEValuesExtractors::Scalar pressure(dim);
2605 *  
2606 * @endcode
2607 *
2608 * Now start the loop over all cells in the problem. We are working on two
2609 * different DoFHandlers for this assembly routine, so we must have two
2610 * different cell iterators for the two objects in use. This might seem a
2611 * bit peculiar, since both the Stokes system and the temperature system
2612 * use the same grid, but that's the only way to keep degrees of freedom
2613 * in sync. The first statements within the loop are again all very
2614 * familiar, doing the update of the finite element data as specified by
2615 * the update flags, zeroing out the local arrays and getting the values
2616 * of the old solution at the quadrature points. Then we are ready to loop
2617 * over the quadrature points on the cell.
2618 *
2619 * @code
2620 *   auto cell = stokes_dof_handler.begin_active();
2621 *   const auto endc = stokes_dof_handler.end();
2622 *   auto temperature_cell = temperature_dof_handler.begin_active();
2623 *  
2624 *   for (; cell != endc; ++cell, ++temperature_cell)
2625 *   {
2626 *   stokes_fe_values.reinit(cell);
2627 *   temperature_fe_values.reinit(temperature_cell);
2628 *  
2629 *   local_matrix = 0;
2630 *   local_rhs = 0;
2631 *  
2632 *   temperature_fe_values.get_function_values(old_temperature_solution,
2633 *   old_temperature_values);
2634 *  
2635 *   for (unsigned int q = 0; q < n_q_points; ++q)
2636 *   {
2637 *   const double old_temperature = old_temperature_values[q];
2638 *  
2639 * @endcode
2640 *
2641 * Next we extract the values and gradients of basis functions
2642 * relevant to the terms in the inner products. As shown in
2643 * @ref step_22 "step-22" this helps accelerate assembly.
2644 *
2645
2646 *
2647 * Once this is done, we start the loop over the rows and columns
2648 * of the local matrix and feed the matrix with the relevant
2649 * products. The right hand side is filled with the forcing term
2650 * driven by temperature in direction of gravity (which is
2651 * vertical in our example). Note that the right hand side term
2652 * is always generated, whereas the matrix contributions are only
2653 * updated when it is requested by the
2654 * <code>rebuild_matrices</code> flag.
2655 *
2656 * @code
2657 *   for (unsigned int k = 0; k < dofs_per_cell; ++k)
2658 *   {
2659 *   phi_u[k] = stokes_fe_values[velocities].value(k, q);
2660 *   if (rebuild_stokes_matrix)
2661 *   {
2662 *   grads_phi_u[k] =
2663 *   stokes_fe_values[velocities].symmetric_gradient(k, q);
2664 *   div_phi_u[k] =
2665 *   stokes_fe_values[velocities].divergence(k, q);
2666 *   phi_p[k] = stokes_fe_values[pressure].value(k, q);
2667 *   }
2668 *   }
2669 *  
2670 *   if (rebuild_stokes_matrix)
2671 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
2672 *   for (unsigned int j = 0; j < dofs_per_cell; ++j)
2673 *   local_matrix(i, j) +=
2674 *   (EquationData::eta * 2 * (grads_phi_u[i] * grads_phi_u[j]) -
2675 *   div_phi_u[i] * phi_p[j] - phi_p[i] * div_phi_u[j]) *
2676 *   stokes_fe_values.JxW(q);
2677 *  
2678 *   const Point<dim> gravity =
2679 *   -((dim == 2) ? (Point<dim>(0, 1)) : (Point<dim>(0, 0, 1)));
2680 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
2681 *   local_rhs(i) += (-EquationData::density * EquationData::beta *
2682 *   gravity * phi_u[i] * old_temperature) *
2683 *   stokes_fe_values.JxW(q);
2684 *   }
2685 *  
2686 * @endcode
2687 *
2688 * The last step in the loop over all cells is to enter the local
2689 * contributions into the global matrix and vector structures to the
2690 * positions specified in <code>local_dof_indices</code>. Again, we
2691 * let the AffineConstraints class do the insertion of the cell
2692 * matrix elements to the global matrix, which already condenses the
2693 * hanging node constraints.
2694 *
2695 * @code
2696 *   cell->get_dof_indices(local_dof_indices);
2697 *  
2698 *   if (rebuild_stokes_matrix == true)
2699 *   stokes_constraints.distribute_local_to_global(local_matrix,
2700 *   local_rhs,
2701 *   local_dof_indices,
2702 *   stokes_matrix,
2703 *   stokes_rhs);
2704 *   else
2705 *   stokes_constraints.distribute_local_to_global(local_rhs,
2706 *   local_dof_indices,
2707 *   stokes_rhs);
2708 *   }
2709 *  
2710 *   rebuild_stokes_matrix = false;
2711 *  
2712 *   std::cout << std::endl;
2713 *   }
2714 *  
2715 *  
2716 *  
2717 * @endcode
2718 *
2719 *
2720 * <a name="step_31-BoussinesqFlowProblemassemble_temperature_matrix"></a>
2721 * <h4>BoussinesqFlowProblem::assemble_temperature_matrix</h4>
2722 *
2723
2724 *
2725 * This function assembles the matrix in the temperature equation. The
2726 * temperature matrix consists of two parts, a mass matrix and the time step
2727 * size times a stiffness matrix given by a Laplace term times the amount of
2728 * diffusion. Since the matrix depends on the time step size (which varies
2729 * from one step to another), the temperature matrix needs to be updated
2730 * every time step. We could simply regenerate the matrices in every time
2731 * step, but this is not really efficient since mass and Laplace matrix do
2732 * only change when we change the mesh. Hence, we do this more efficiently
2733 * by generating two separate matrices in this function, one for the mass
2734 * matrix and one for the stiffness (diffusion) matrix. We will then sum up
2735 * the matrix plus the stiffness matrix times the time step size once we
2736 * know the actual time step.
2737 *
2738
2739 *
2740 * So the details for this first step are very simple. In case we need to
2741 * rebuild the matrix (i.e., the mesh has changed), we zero the data
2742 * structures, get a quadrature formula and a FEValues object, and create
2743 * local matrices, local dof indices and evaluation structures for the basis
2744 * functions.
2745 *
2746 * @code
2747 *   template <int dim>
2748 *   void BoussinesqFlowProblem<dim>::assemble_temperature_matrix()
2749 *   {
2750 *   if (rebuild_temperature_matrices == false)
2751 *   return;
2752 *  
2753 *   temperature_mass_matrix = 0;
2754 *   temperature_stiffness_matrix = 0;
2755 *  
2756 *   QGauss<dim> quadrature_formula(temperature_degree + 2);
2757 *   FEValues<dim> temperature_fe_values(temperature_fe,
2758 *   quadrature_formula,
2759 *   update_values | update_gradients |
2760 *   update_JxW_values);
2761 *  
2762 *   const unsigned int dofs_per_cell = temperature_fe.n_dofs_per_cell();
2763 *   const unsigned int n_q_points = quadrature_formula.size();
2764 *  
2765 *   FullMatrix<double> local_mass_matrix(dofs_per_cell, dofs_per_cell);
2766 *   FullMatrix<double> local_stiffness_matrix(dofs_per_cell, dofs_per_cell);
2767 *  
2768 *   std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
2769 *  
2770 *   std::vector<double> phi_T(dofs_per_cell);
2771 *   std::vector<Tensor<1, dim>> grad_phi_T(dofs_per_cell);
2772 *  
2773 * @endcode
2774 *
2775 * Now, let's start the loop over all cells in the triangulation. We need
2776 * to zero out the local matrices, update the finite element evaluations,
2777 * and then loop over the rows and columns of the matrices on each
2778 * quadrature point, where we then create the mass matrix and the
2779 * stiffness matrix (Laplace terms times the diffusion
2780 * <code>EquationData::kappa</code>. Finally, we let the constraints
2781 * object insert these values into the global matrix, and directly
2782 * condense the constraints into the matrix.
2783 *
2784 * @code
2785 *   for (const auto &cell : temperature_dof_handler.active_cell_iterators())
2786 *   {
2787 *   local_mass_matrix = 0;
2788 *   local_stiffness_matrix = 0;
2789 *  
2790 *   temperature_fe_values.reinit(cell);
2791 *  
2792 *   for (unsigned int q = 0; q < n_q_points; ++q)
2793 *   {
2794 *   for (unsigned int k = 0; k < dofs_per_cell; ++k)
2795 *   {
2796 *   grad_phi_T[k] = temperature_fe_values.shape_grad(k, q);
2797 *   phi_T[k] = temperature_fe_values.shape_value(k, q);
2798 *   }
2799 *  
2800 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
2801 *   for (unsigned int j = 0; j < dofs_per_cell; ++j)
2802 *   {
2803 *   local_mass_matrix(i, j) +=
2804 *   (phi_T[i] * phi_T[j] * temperature_fe_values.JxW(q));
2805 *   local_stiffness_matrix(i, j) +=
2806 *   (EquationData::kappa * grad_phi_T[i] * grad_phi_T[j] *
2807 *   temperature_fe_values.JxW(q));
2808 *   }
2809 *   }
2810 *  
2811 *   cell->get_dof_indices(local_dof_indices);
2812 *  
2813 *   temperature_constraints.distribute_local_to_global(
2814 *   local_mass_matrix, local_dof_indices, temperature_mass_matrix);
2815 *   temperature_constraints.distribute_local_to_global(
2816 *   local_stiffness_matrix,
2817 *   local_dof_indices,
2818 *   temperature_stiffness_matrix);
2819 *   }
2820 *  
2821 *   rebuild_temperature_matrices = false;
2822 *   }
2823 *  
2824 *  
2825 *  
2826 * @endcode
2827 *
2828 *
2829 * <a name="step_31-BoussinesqFlowProblemassemble_temperature_system"></a>
2830 * <h4>BoussinesqFlowProblem::assemble_temperature_system</h4>
2831 *
2832
2833 *
2834 * This function does the second part of the assembly work on the
2835 * temperature matrix, the actual addition of pressure mass and stiffness
2836 * matrix (where the time step size comes into play), as well as the
2837 * creation of the velocity-dependent right hand side. The declarations for
2838 * the right hand side assembly in this function are pretty much the same as
2839 * the ones used in the other assembly routines, except that we restrict
2840 * ourselves to vectors this time. We are going to calculate residuals on
2841 * the temperature system, which means that we have to evaluate second
2842 * derivatives, specified by the update flag <code>update_hessians</code>.
2843 *
2844
2845 *
2846 * The temperature equation is coupled to the Stokes system by means of the
2847 * fluid velocity. These two parts of the solution are associated with
2848 * different DoFHandlers, so we again need to create a second FEValues
2849 * object for the evaluation of the velocity at the quadrature points.
2850 *
2851 * @code
2852 *   template <int dim>
2853 *   void BoussinesqFlowProblem<dim>::assemble_temperature_system(
2854 *   const double maximal_velocity)
2855 *   {
2856 *   const bool use_bdf2_scheme = (timestep_number != 0);
2857 *  
2858 *   if (use_bdf2_scheme == true)
2859 *   {
2860 *   temperature_matrix.copy_from(temperature_mass_matrix);
2861 *   temperature_matrix *=
2862 *   (2 * time_step + old_time_step) / (time_step + old_time_step);
2863 *   temperature_matrix.add(time_step, temperature_stiffness_matrix);
2864 *   }
2865 *   else
2866 *   {
2867 *   temperature_matrix.copy_from(temperature_mass_matrix);
2868 *   temperature_matrix.add(time_step, temperature_stiffness_matrix);
2869 *   }
2870 *  
2871 *   temperature_rhs = 0;
2872 *  
2873 *   const QGauss<dim> quadrature_formula(temperature_degree + 2);
2874 *   FEValues<dim> temperature_fe_values(temperature_fe,
2875 *   quadrature_formula,
2877 *   update_hessians |
2879 *   update_JxW_values);
2880 *   FEValues<dim> stokes_fe_values(stokes_fe,
2881 *   quadrature_formula,
2882 *   update_values);
2883 *  
2884 *   const unsigned int dofs_per_cell = temperature_fe.n_dofs_per_cell();
2885 *   const unsigned int n_q_points = quadrature_formula.size();
2886 *  
2887 *   Vector<double> local_rhs(dofs_per_cell);
2888 *  
2889 *   std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
2890 *  
2891 * @endcode
2892 *
2893 * Next comes the declaration of vectors to hold the old and older
2894 * solution values (as a notation for time levels @f$n-1@f$ and
2895 * @f$n-2@f$, respectively) and gradients at quadrature points of the
2896 * current cell. We also declare an object to hold the temperature right
2897 * hand side values (<code>gamma_values</code>), and we again use
2898 * shortcuts for the temperature basis functions. Eventually, we need to
2899 * find the temperature extrema and the diameter of the computational
2900 * domain which will be used for the definition of the stabilization
2901 * parameter (we got the maximal velocity as an input to this function).
2902 *
2903 * @code
2904 *   std::vector<Tensor<1, dim>> old_velocity_values(n_q_points);
2905 *   std::vector<Tensor<1, dim>> old_old_velocity_values(n_q_points);
2906 *   std::vector<double> old_temperature_values(n_q_points);
2907 *   std::vector<double> old_old_temperature_values(n_q_points);
2908 *   std::vector<Tensor<1, dim>> old_temperature_grads(n_q_points);
2909 *   std::vector<Tensor<1, dim>> old_old_temperature_grads(n_q_points);
2910 *   std::vector<double> old_temperature_laplacians(n_q_points);
2911 *   std::vector<double> old_old_temperature_laplacians(n_q_points);
2912 *  
2913 *   EquationData::TemperatureRightHandSide<dim> temperature_right_hand_side;
2914 *   std::vector<double> gamma_values(n_q_points);
2915 *  
2916 *   std::vector<double> phi_T(dofs_per_cell);
2917 *   std::vector<Tensor<1, dim>> grad_phi_T(dofs_per_cell);
2918 *  
2919 *   const std::pair<double, double> global_T_range =
2920 *   get_extrapolated_temperature_range();
2921 *  
2922 *   const FEValuesExtractors::Vector velocities(0);
2923 *  
2924 * @endcode
2925 *
2926 * Now, let's start the loop over all cells in the triangulation. Again,
2927 * we need two cell iterators that walk in parallel through the cells of
2928 * the two involved DoFHandler objects for the Stokes and temperature
2929 * part. Within the loop, we first set the local rhs to zero, and then get
2930 * the values and derivatives of the old solution functions at the
2931 * quadrature points, since they are going to be needed for the definition
2932 * of the stabilization parameters and as coefficients in the equation,
2933 * respectively. Note that since the temperature has its own DoFHandler
2934 * and FEValues object we get the entire solution at the quadrature point
2935 * (which is the scalar temperature field only anyway) whereas for the
2936 * Stokes part we restrict ourselves to extracting the velocity part (and
2937 * ignoring the pressure part) by using
2938 * <code>stokes_fe_values[velocities].get_function_values</code>.
2939 *
2940 * @code
2941 *   auto cell = temperature_dof_handler.begin_active();
2942 *   const auto endc = temperature_dof_handler.end();
2943 *   auto stokes_cell = stokes_dof_handler.begin_active();
2944 *  
2945 *   for (; cell != endc; ++cell, ++stokes_cell)
2946 *   {
2947 *   local_rhs = 0;
2948 *  
2949 *   temperature_fe_values.reinit(cell);
2950 *   stokes_fe_values.reinit(stokes_cell);
2951 *  
2952 *   temperature_fe_values.get_function_values(old_temperature_solution,
2953 *   old_temperature_values);
2954 *   temperature_fe_values.get_function_values(old_old_temperature_solution,
2955 *   old_old_temperature_values);
2956 *  
2957 *   temperature_fe_values.get_function_gradients(old_temperature_solution,
2958 *   old_temperature_grads);
2959 *   temperature_fe_values.get_function_gradients(
2960 *   old_old_temperature_solution, old_old_temperature_grads);
2961 *  
2962 *   temperature_fe_values.get_function_laplacians(
2963 *   old_temperature_solution, old_temperature_laplacians);
2964 *   temperature_fe_values.get_function_laplacians(
2965 *   old_old_temperature_solution, old_old_temperature_laplacians);
2966 *  
2967 *   temperature_right_hand_side.value_list(
2968 *   temperature_fe_values.get_quadrature_points(), gamma_values);
2969 *  
2970 *   stokes_fe_values[velocities].get_function_values(stokes_solution,
2971 *   old_velocity_values);
2972 *   stokes_fe_values[velocities].get_function_values(
2973 *   old_stokes_solution, old_old_velocity_values);
2974 *  
2975 * @endcode
2976 *
2977 * Next, we calculate the artificial viscosity for stabilization
2978 * according to the discussion in the introduction using the dedicated
2979 * function. With that at hand, we can get into the loop over
2980 * quadrature points and local rhs vector components. The terms here
2981 * are quite lengthy, but their definition follows the time-discrete
2982 * system developed in the introduction of this program. The BDF-2
2983 * scheme needs one more term from the old time step (and involves
2984 * more complicated factors) than the backward Euler scheme that is
2985 * used for the first time step. When all this is done, we distribute
2986 * the local vector into the global one (including hanging node
2987 * constraints).
2988 *
2989 * @code
2990 *   const double nu =
2991 *   compute_viscosity(old_temperature_values,
2992 *   old_old_temperature_values,
2993 *   old_temperature_grads,
2994 *   old_old_temperature_grads,
2995 *   old_temperature_laplacians,
2996 *   old_old_temperature_laplacians,
2997 *   old_velocity_values,
2998 *   old_old_velocity_values,
2999 *   gamma_values,
3000 *   maximal_velocity,
3001 *   global_T_range.second - global_T_range.first,
3002 *   cell->diameter());
3003 *  
3004 *   for (unsigned int q = 0; q < n_q_points; ++q)
3005 *   {
3006 *   for (unsigned int k = 0; k < dofs_per_cell; ++k)
3007 *   {
3008 *   grad_phi_T[k] = temperature_fe_values.shape_grad(k, q);
3009 *   phi_T[k] = temperature_fe_values.shape_value(k, q);
3010 *   }
3011 *  
3012 *   const double T_term_for_rhs =
3013 *   (use_bdf2_scheme ?
3014 *   (old_temperature_values[q] * (1 + time_step / old_time_step) -
3015 *   old_old_temperature_values[q] * (time_step * time_step) /
3016 *   (old_time_step * (time_step + old_time_step))) :
3017 *   old_temperature_values[q]);
3018 *  
3019 *   const Tensor<1, dim> ext_grad_T =
3020 *   (use_bdf2_scheme ?
3021 *   (old_temperature_grads[q] * (1 + time_step / old_time_step) -
3022 *   old_old_temperature_grads[q] * time_step / old_time_step) :
3023 *   old_temperature_grads[q]);
3024 *  
3025 *   const Tensor<1, dim> extrapolated_u =
3026 *   (use_bdf2_scheme ?
3027 *   (old_velocity_values[q] * (1 + time_step / old_time_step) -
3028 *   old_old_velocity_values[q] * time_step / old_time_step) :
3029 *   old_velocity_values[q]);
3030 *  
3031 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
3032 *   local_rhs(i) +=
3033 *   (T_term_for_rhs * phi_T[i] -
3034 *   time_step * extrapolated_u * ext_grad_T * phi_T[i] -
3035 *   time_step * nu * ext_grad_T * grad_phi_T[i] +
3036 *   time_step * gamma_values[q] * phi_T[i]) *
3037 *   temperature_fe_values.JxW(q);
3038 *   }
3039 *  
3040 *   cell->get_dof_indices(local_dof_indices);
3041 *   temperature_constraints.distribute_local_to_global(local_rhs,
3042 *   local_dof_indices,
3043 *   temperature_rhs);
3044 *   }
3045 *   }
3046 *  
3047 *  
3048 *  
3049 * @endcode
3050 *
3051 *
3052 * <a name="step_31-BoussinesqFlowProblemsolve"></a>
3053 * <h4>BoussinesqFlowProblem::solve</h4>
3054 *
3055
3056 *
3057 * This function solves the linear systems of equations. Following the
3058 * introduction, we start with the Stokes system, where we need to generate
3059 * our block Schur preconditioner. Since all the relevant actions are
3060 * implemented in the class <code>BlockSchurPreconditioner</code>, all we
3061 * have to do is to initialize the class appropriately. What we need to pass
3062 * down is an <code>InverseMatrix</code> object for the pressure mass
3063 * matrix, which we set up using the respective class together with the IC
3064 * preconditioner we already generated, and the AMG preconditioner for the
3065 * velocity-velocity matrix. Note that both <code>Mp_preconditioner</code>
3066 * and <code>Amg_preconditioner</code> are only pointers, so we use
3067 * <code>*</code> to pass down the actual preconditioner objects.
3068 *
3069
3070 *
3071 * Once the preconditioner is ready, we create a GMRES solver for the block
3072 * system. Since we are working with Trilinos data structures, we have to
3073 * set the respective template argument in the solver. GMRES needs to
3074 * internally store temporary vectors for each iteration (see the discussion
3075 * in the results section of @ref step_22 "step-22") &ndash; the more vectors it can use,
3076 * the better it will generally perform. To keep memory demands in check, we
3077 * set the number of vectors to 100. This means that up to 100 solver
3078 * iterations, every temporary vector can be stored. If the solver needs to
3079 * iterate more often to get the specified tolerance, it will work on a
3080 * reduced set of vectors by restarting at every 100 iterations.
3081 *
3082
3083 *
3084 * With this all set up, we solve the system and distribute the constraints
3085 * in the Stokes system, i.e., hanging nodes and no-flux boundary condition,
3086 * in order to have the appropriate solution values even at constrained
3087 * dofs. Finally, we write the number of iterations to the screen.
3088 *
3089 * @code
3090 *   template <int dim>
3091 *   void BoussinesqFlowProblem<dim>::solve()
3092 *   {
3093 *   std::cout << " Solving..." << std::endl;
3094 *  
3095 *   {
3096 *   const LinearSolvers::InverseMatrix<TrilinosWrappers::SparseMatrix,
3097 *   TrilinosWrappers::PreconditionIC>
3098 *   mp_inverse(stokes_preconditioner_matrix.block(1, 1),
3099 *   *Mp_preconditioner);
3100 *  
3101 *   const LinearSolvers::BlockSchurPreconditioner<
3102 *   TrilinosWrappers::PreconditionAMG,
3103 *   TrilinosWrappers::PreconditionIC>
3104 *   preconditioner(stokes_matrix, mp_inverse, *Amg_preconditioner);
3105 *  
3106 *   SolverControl solver_control(stokes_matrix.m(),
3107 *   1e-6 * stokes_rhs.l2_norm());
3108 *  
3109 *   SolverGMRES<TrilinosWrappers::MPI::BlockVector> gmres(
3110 *   solver_control,
3111 *   SolverGMRES<TrilinosWrappers::MPI::BlockVector>::AdditionalData(100));
3112 *  
3113 *   for (unsigned int i = 0; i < stokes_solution.size(); ++i)
3114 *   if (stokes_constraints.is_constrained(i))
3115 *   stokes_solution(i) = 0;
3116 *  
3117 *   gmres.solve(stokes_matrix, stokes_solution, stokes_rhs, preconditioner);
3118 *  
3119 *   stokes_constraints.distribute(stokes_solution);
3120 *  
3121 *   std::cout << " " << solver_control.last_step()
3122 *   << " GMRES iterations for Stokes subsystem." << std::endl;
3123 *   }
3124 *  
3125 * @endcode
3126 *
3127 * Once we know the Stokes solution, we can determine the new time step
3128 * from the maximal velocity. We have to do this to satisfy the CFL
3129 * condition since convection terms are treated explicitly in the
3130 * temperature equation, as discussed in the introduction. The exact form
3131 * of the formula used here for the time step is discussed in the results
3132 * section of this program.
3133 *
3134
3135 *
3136 * There is a snatch here. The formula contains a division by the maximum
3137 * value of the velocity. However, at the start of the computation, we
3138 * have a constant temperature field (we start with a constant
3139 * temperature, and it will be nonconstant only after the first time step
3140 * during which the source acts). Constant temperature means that no
3141 * buoyancy acts, and so the velocity is zero. Dividing by it will not
3142 * likely lead to anything good.
3143 *
3144
3145 *
3146 * To avoid the resulting infinite time step, we ask whether the maximal
3147 * velocity is very small (in particular smaller than the values we
3148 * encounter during any of the following time steps) and if so rather than
3149 * dividing by zero we just divide by a small value, resulting in a large
3150 * but finite time step.
3151 *
3152 * @code
3153 *   old_time_step = time_step;
3154 *   const double maximal_velocity = get_maximal_velocity();
3155 *  
3156 *   if (maximal_velocity >= 0.01)
3157 *   time_step = 1. / (1.7 * dim * std::sqrt(1. * dim)) / temperature_degree *
3158 *   GridTools::minimal_cell_diameter(triangulation) /
3159 *   maximal_velocity;
3160 *   else
3161 *   time_step = 1. / (1.7 * dim * std::sqrt(1. * dim)) / temperature_degree *
3162 *   GridTools::minimal_cell_diameter(triangulation) / .01;
3163 *  
3164 *   std::cout << " "
3165 *   << "Time step: " << time_step << std::endl;
3166 *  
3167 *   temperature_solution = old_temperature_solution;
3168 *  
3169 * @endcode
3170 *
3171 * Next we set up the temperature system and the right hand side using the
3172 * function <code>assemble_temperature_system()</code>. Knowing the
3173 * matrix and right hand side of the temperature equation, we set up a
3174 * preconditioner and a solver. The temperature matrix is a mass matrix
3175 * (with eigenvalues around one) plus a Laplace matrix (with eigenvalues
3176 * between zero and @f$ch^{-2}@f$) times a small number proportional to the
3177 * time step @f$k_n@f$. Hence, the resulting symmetric and positive definite
3178 * matrix has eigenvalues in the range @f$[1,1+k_nh^{-2}]@f$ (up to
3179 * constants). This matrix is only moderately ill conditioned even for
3180 * small mesh sizes and we get a reasonably good preconditioner by simple
3181 * means, for example with an incomplete Cholesky decomposition
3182 * preconditioner (IC) as we also use for preconditioning the pressure
3183 * mass matrix solver. As a solver, we choose the conjugate gradient
3184 * method CG. As before, we tell the solver to use Trilinos vectors via
3185 * the template argument <code>TrilinosWrappers::MPI::Vector</code>.
3186 * Finally, we solve, distribute the hanging node constraints and write out
3187 * the number of iterations.
3188 *
3189 * @code
3190 *   assemble_temperature_system(maximal_velocity);
3191 *   {
3192 *   SolverControl solver_control(temperature_matrix.m(),
3193 *   1e-8 * temperature_rhs.l2_norm());
3194 *   SolverCG<TrilinosWrappers::MPI::Vector> cg(solver_control);
3195 *  
3196 *   TrilinosWrappers::PreconditionIC preconditioner;
3197 *   preconditioner.initialize(temperature_matrix);
3198 *  
3199 *   cg.solve(temperature_matrix,
3200 *   temperature_solution,
3201 *   temperature_rhs,
3202 *   preconditioner);
3203 *  
3204 *   temperature_constraints.distribute(temperature_solution);
3205 *  
3206 *   std::cout << " " << solver_control.last_step()
3207 *   << " CG iterations for temperature." << std::endl;
3208 *  
3209 * @endcode
3210 *
3211 * At the end of this function, we step through the vector and read out
3212 * the maximum and minimum temperature value, which we also want to
3213 * output. This will come in handy when determining the correct constant
3214 * in the choice of time step as discuss in the results section of this
3215 * program.
3216 *
3217 * @code
3218 *   double min_temperature = temperature_solution(0),
3219 *   max_temperature = temperature_solution(0);
3220 *   for (unsigned int i = 0; i < temperature_solution.size(); ++i)
3221 *   {
3222 *   min_temperature =
3223 *   std::min<double>(min_temperature, temperature_solution(i));
3224 *   max_temperature =
3225 *   std::max<double>(max_temperature, temperature_solution(i));
3226 *   }
3227 *  
3228 *   std::cout << " Temperature range: " << min_temperature << ' '
3229 *   << max_temperature << std::endl;
3230 *   }
3231 *   }
3232 *  
3233 *  
3234 *  
3235 * @endcode
3236 *
3237 *
3238 * <a name="step_31-BoussinesqFlowProblemoutput_results"></a>
3239 * <h4>BoussinesqFlowProblem::output_results</h4>
3240 *
3241
3242 *
3243 * This function writes the solution to a VTK output file for visualization,
3244 * which is done every tenth time step. This is usually quite a simple task,
3245 * since the deal.II library provides functions that do almost all the job
3246 * for us. There is one new function compared to previous examples: We want
3247 * to visualize both the Stokes solution and the temperature as one data
3248 * set, but we have done all the calculations based on two different
3249 * DoFHandler objects. Luckily, the DataOut class is prepared to deal with
3250 * it. All we have to do is to not attach one single DoFHandler at the
3251 * beginning and then use that for all added vector, but specify the
3252 * DoFHandler to each vector separately. The rest is done as in @ref step_22 "step-22". We
3253 * create solution names (that are going to appear in the visualization
3254 * program for the individual components). The first <code>dim</code>
3255 * components are the vector velocity, and then we have pressure for the
3256 * Stokes part, whereas temperature is scalar. This information is read out
3257 * using the DataComponentInterpretation helper class. Next, we actually
3258 * attach the data vectors with their DoFHandler objects, build patches
3259 * according to the degree of freedom, which are (sub-) elements that
3260 * describe the data for visualization programs. Finally, we open a file
3261 * (that includes the time step number) and write the vtk data into it.
3262 *
3263 * @code
3264 *   template <int dim>
3265 *   void BoussinesqFlowProblem<dim>::output_results() const
3266 *   {
3267 *   if (timestep_number % 10 != 0)
3268 *   return;
3269 *  
3270 *   std::vector<std::string> stokes_names(dim, "velocity");
3271 *   stokes_names.emplace_back("p");
3272 *   std::vector<DataComponentInterpretation::DataComponentInterpretation>
3273 *   stokes_component_interpretation(
3274 *   dim + 1, DataComponentInterpretation::component_is_scalar);
3275 *   for (unsigned int i = 0; i < dim; ++i)
3276 *   stokes_component_interpretation[i] =
3277 *   DataComponentInterpretation::component_is_part_of_vector;
3278 *  
3279 *   DataOut<dim> data_out;
3280 *   data_out.add_data_vector(stokes_dof_handler,
3281 *   stokes_solution,
3282 *   stokes_names,
3283 *   stokes_component_interpretation);
3284 *   data_out.add_data_vector(temperature_dof_handler,
3285 *   temperature_solution,
3286 *   "T");
3287 *   data_out.build_patches(std::min(stokes_degree, temperature_degree));
3288 *  
3289 *   std::ofstream output("solution-" +
3290 *   Utilities::int_to_string(timestep_number, 4) + ".vtk");
3291 *   data_out.write_vtk(output);
3292 *   }
3293 *  
3294 *  
3295 *  
3296 * @endcode
3297 *
3298 *
3299 * <a name="step_31-BoussinesqFlowProblemrefine_mesh"></a>
3300 * <h4>BoussinesqFlowProblem::refine_mesh</h4>
3301 *
3302
3303 *
3304 * This function takes care of the adaptive mesh refinement. The three tasks
3305 * this function performs is to first find out which cells to
3306 * refine/coarsen, then to actually do the refinement and eventually
3307 * transfer the solution vectors between the two different grids. The first
3308 * task is simply achieved by using the well-established Kelly error
3309 * estimator on the temperature (it is the temperature we're mainly
3310 * interested in for this program, and we need to be accurate in regions of
3311 * high temperature gradients, also to not have too much numerical
3312 * diffusion). The second task is to actually do the remeshing. That
3313 * involves only basic functions as well, such as the
3314 * <code>refine_and_coarsen_fixed_fraction</code> that refines those cells
3315 * with the largest estimated error that together make up 80 per cent of the
3316 * error, and coarsens those cells with the smallest error that make up for
3317 * a combined 10 per cent of the error.
3318 *
3319
3320 *
3321 * If implemented like this, we would get a program that will not make much
3322 * progress: Remember that we expect temperature fields that are nearly
3323 * discontinuous (the diffusivity @f$\kappa@f$ is very small after all) and
3324 * consequently we can expect that a freely adapted mesh will refine further
3325 * and further into the areas of large gradients. This decrease in mesh size
3326 * will then be accompanied by a decrease in time step, requiring an
3327 * exceedingly large number of time steps to solve to a given final time. It
3328 * will also lead to meshes that are much better at resolving
3329 * discontinuities after several mesh refinement cycles than in the
3330 * beginning.
3331 *
3332
3333 *
3334 * In particular to prevent the decrease in time step size and the
3335 * correspondingly large number of time steps, we limit the maximal
3336 * refinement depth of the mesh. To this end, after the refinement indicator
3337 * has been applied to the cells, we simply loop over all cells on the
3338 * finest level and unselect them from refinement if they would result in
3339 * too high a mesh level.
3340 *
3341 * @code
3342 *   template <int dim>
3343 *   void
3344 *   BoussinesqFlowProblem<dim>::refine_mesh(const unsigned int max_grid_level)
3345 *   {
3346 *   Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
3347 *  
3348 *   KellyErrorEstimator<dim>::estimate(temperature_dof_handler,
3349 *   QGauss<dim - 1>(temperature_degree + 1),
3350 *   {},
3351 *   temperature_solution,
3352 *   estimated_error_per_cell);
3353 *  
3355 *   estimated_error_per_cell,
3356 *   0.8,
3357 *   0.1);
3358 *   if (triangulation.n_levels() > max_grid_level)
3359 *   for (auto &cell :
3360 *   triangulation.active_cell_iterators_on_level(max_grid_level))
3361 *   cell->clear_refine_flag();
3362 *  
3363 * @endcode
3364 *
3365 * As part of mesh refinement we need to transfer the solution vectors
3366 * from the old mesh to the new one. To this end we use the
3367 * SolutionTransfer class and we have to prepare the solution vectors that
3368 * should be transferred to the new grid (we will lose the old grid once
3369 * we have done the refinement so the transfer has to happen concurrently
3370 * with refinement). What we definitely need are the current and the old
3371 * temperature (BDF-2 time stepping requires two old solutions). Since the
3372 * SolutionTransfer objects only support to transfer one object per dof
3373 * handler, we need to collect the two temperature solutions in one data
3374 * structure. Moreover, we choose to transfer the Stokes solution, too,
3375 * since we need the velocity at two previous time steps, of which only
3376 * one is calculated on the fly.
3377 *
3378
3379 *
3380 * Consequently, we initialize two SolutionTransfer objects for the Stokes
3381 * and temperature DoFHandler objects, by attaching them to the old dof
3382 * handlers. With this at place, we can prepare the triangulation and the
3383 * data vectors for refinement (in this order).
3384 *
3385 * @code
3386 *   const std::vector<TrilinosWrappers::MPI::Vector> x_temperature = {
3387 *   temperature_solution, old_temperature_solution};
3388 *   TrilinosWrappers::MPI::BlockVector x_stokes = stokes_solution;
3389 *  
3391 *   temperature_dof_handler);
3393 *   stokes_dof_handler);
3394 *  
3395 *   triangulation.prepare_coarsening_and_refinement();
3396 *   temperature_trans.prepare_for_coarsening_and_refinement(x_temperature);
3397 *   stokes_trans.prepare_for_coarsening_and_refinement(x_stokes);
3398 *  
3399 * @endcode
3400 *
3401 * Now everything is ready, so do the refinement and recreate the dof
3402 * structure on the new grid, and initialize the matrix structures and the
3403 * new vectors in the <code>setup_dofs</code> function. Next, we actually
3404 * perform the interpolation of the solutions between the grids. We create
3405 * another copy of temporary vectors for temperature (now corresponding to
3406 * the new grid), and let the interpolate function do the job. Then, the
3407 * resulting array of vectors is written into the respective vector member
3408 * variables.
3409 *
3410
3411 *
3412 * Remember that the set of constraints will be updated for the new
3413 * triangulation in the setup_dofs() call.
3414 *
3415 * @code
3416 *   triangulation.execute_coarsening_and_refinement();
3417 *   setup_dofs();
3418 *  
3419 *   std::vector<TrilinosWrappers::MPI::Vector> tmp = {
3420 *   TrilinosWrappers::MPI::Vector(temperature_solution),
3421 *   TrilinosWrappers::MPI::Vector(temperature_solution)};
3422 *   temperature_trans.interpolate(x_temperature, tmp);
3423 *  
3424 *   temperature_solution = tmp[0];
3425 *   old_temperature_solution = tmp[1];
3426 *  
3427 * @endcode
3428 *
3429 * After the solution has been transferred we then enforce the constraints
3430 * on the transferred solution.
3431 *
3432 * @code
3433 *   temperature_constraints.distribute(temperature_solution);
3434 *   temperature_constraints.distribute(old_temperature_solution);
3435 *  
3436 * @endcode
3437 *
3438 * For the Stokes vector, everything is just the same &ndash; except that
3439 * we do not need another temporary vector since we just interpolate a
3440 * single vector. In the end, we have to tell the program that the matrices
3441 * and preconditioners need to be regenerated, since the mesh has changed.
3442 *
3443 * @code
3444 *   stokes_trans.interpolate(x_stokes, stokes_solution);
3445 *  
3446 *   stokes_constraints.distribute(stokes_solution);
3447 *  
3448 *   rebuild_stokes_matrix = true;
3449 *   rebuild_temperature_matrices = true;
3450 *   rebuild_stokes_preconditioner = true;
3451 *   }
3452 *  
3453 *  
3454 *  
3455 * @endcode
3456 *
3457 *
3458 * <a name="step_31-BoussinesqFlowProblemrun"></a>
3459 * <h4>BoussinesqFlowProblem::run</h4>
3460 *
3461
3462 *
3463 * This function performs all the essential steps in the Boussinesq
3464 * program. It starts by setting up a grid (depending on the spatial
3465 * dimension, we choose some different level of initial refinement and
3466 * additional adaptive refinement steps, and then create a cube in
3467 * <code>dim</code> dimensions and set up the dofs for the first time. Since
3468 * we want to start the time stepping already with an adaptively refined
3469 * grid, we perform some pre-refinement steps, consisting of all assembly,
3470 * solution and refinement, but without actually advancing in time. Rather,
3471 * we use the vilified <code>goto</code> statement to jump out of the time
3472 * loop right after mesh refinement to start all over again on the new mesh
3473 * beginning at the <code>start_time_iteration</code> label. (The use of the
3474 * <code>goto</code> is discussed in @ref step_26 "step-26".)
3475 *
3476
3477 *
3478 * Before we start, we project the initial values to the grid and obtain the
3479 * first data for the <code>old_temperature_solution</code> vector. Then, we
3480 * initialize time step number and time step and start the time loop.
3481 *
3482 * @code
3483 *   template <int dim>
3484 *   void BoussinesqFlowProblem<dim>::run()
3485 *   {
3486 *   const unsigned int initial_refinement = (dim == 2 ? 4 : 2);
3487 *   const unsigned int n_pre_refinement_steps = (dim == 2 ? 4 : 3);
3488 *  
3489 *  
3490 *   GridGenerator::hyper_cube(triangulation);
3491 *   global_Omega_diameter = GridTools::diameter(triangulation);
3492 *  
3493 *   triangulation.refine_global(initial_refinement);
3494 *  
3495 *   setup_dofs();
3496 *  
3497 *   unsigned int pre_refinement_step = 0;
3498 *  
3499 *   start_time_iteration:
3500 *  
3501 *   VectorTools::project(temperature_dof_handler,
3502 *   temperature_constraints,
3503 *   QGauss<dim>(temperature_degree + 2),
3504 *   EquationData::TemperatureInitialValues<dim>(),
3505 *   old_temperature_solution);
3506 *  
3507 *   timestep_number = 0;
3508 *   time_step = old_time_step = 0;
3509 *  
3510 *   double time = 0;
3511 *  
3512 *   do
3513 *   {
3514 *   std::cout << "Timestep " << timestep_number << ": t=" << time
3515 *   << std::endl;
3516 *  
3517 * @endcode
3518 *
3519 * The first steps in the time loop are all obvious &ndash; we
3520 * assemble the Stokes system, the preconditioner, the temperature
3521 * matrix (matrices and preconditioner do actually only change in case
3522 * we've remeshed before), and then do the solve. Before going on with
3523 * the next time step, we have to check whether we should first finish
3524 * the pre-refinement steps or if we should remesh (every fifth time
3525 * step), refining up to a level that is consistent with initial
3526 * refinement and pre-refinement steps. Last in the loop is to advance
3527 * the solutions, i.e., to copy the solutions to the next "older" time
3528 * level.
3529 *
3530 * @code
3531 *   assemble_stokes_system();
3532 *   build_stokes_preconditioner();
3533 *   assemble_temperature_matrix();
3534 *  
3535 *   solve();
3536 *  
3537 *   output_results();
3538 *  
3539 *   std::cout << std::endl;
3540 *  
3541 *   if ((timestep_number == 0) &&
3542 *   (pre_refinement_step < n_pre_refinement_steps))
3543 *   {
3544 *   refine_mesh(initial_refinement + n_pre_refinement_steps);
3545 *   ++pre_refinement_step;
3546 *   goto start_time_iteration;
3547 *   }
3548 *   else if ((timestep_number > 0) && (timestep_number % 5 == 0))
3549 *   refine_mesh(initial_refinement + n_pre_refinement_steps);
3550 *  
3551 *   time += time_step;
3552 *   ++timestep_number;
3553 *  
3554 *   old_stokes_solution = stokes_solution;
3555 *   old_old_temperature_solution = old_temperature_solution;
3556 *   old_temperature_solution = temperature_solution;
3557 *   }
3558 * @endcode
3559 *
3560 * Do all the above until we arrive at time 100.
3561 *
3562 * @code
3563 *   while (time <= 100);
3564 *   }
3565 *   } // namespace Step31
3566 *  
3567 *  
3568 *  
3569 * @endcode
3570 *
3571 *
3572 * <a name="step_31-Thecodemaincodefunction"></a>
3573 * <h3>The <code>main</code> function</h3>
3574 *
3575
3576 *
3577 * The main function looks almost the same as in all other programs.
3578 *
3579
3580 *
3581 * There is one difference we have to be careful about. This program uses
3582 * Trilinos and, typically, Trilinos is configured so that it can run in
3583 * %parallel using MPI. This doesn't mean that it <i>has</i> to run in
3584 * %parallel, and in fact this program (unlike @ref step_32 "step-32") makes no attempt at
3585 * all to do anything in %parallel using MPI. Nevertheless, Trilinos wants the
3586 * MPI system to be initialized. We do that be creating an object of type
3587 * Utilities::MPI::MPI_InitFinalize that initializes MPI (if available) using
3588 * the arguments given to main() (i.e., <code>argc</code> and
3589 * <code>argv</code>) and de-initializes it again when the object goes out of
3590 * scope.
3591 *
3592 * @code
3593 *   int main(int argc, char *argv[])
3594 *   {
3595 *   try
3596 *   {
3597 *   using namespace dealii;
3598 *   using namespace Step31;
3599 *  
3600 *   Utilities::MPI::MPI_InitFinalize mpi_initialization(
3601 *   argc, argv, numbers::invalid_unsigned_int);
3602 *  
3603 * @endcode
3604 *
3605 * This program can only be run in serial. Otherwise, throw an exception.
3606 *
3607 * @code
3608 *   AssertThrow(Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD) == 1,
3609 *   ExcMessage(
3610 *   "This program can only be run in serial, use ./step-31"));
3611 *  
3612 *   BoussinesqFlowProblem<2> flow_problem;
3613 *   flow_problem.run();
3614 *   }
3615 *   catch (std::exception &exc)
3616 *   {
3617 *   std::cerr << std::endl
3618 *   << std::endl
3619 *   << "----------------------------------------------------"
3620 *   << std::endl;
3621 *   std::cerr << "Exception on processing: " << std::endl
3622 *   << exc.what() << std::endl
3623 *   << "Aborting!" << std::endl
3624 *   << "----------------------------------------------------"
3625 *   << std::endl;
3626 *  
3627 *   return 1;
3628 *   }
3629 *   catch (...)
3630 *   {
3631 *   std::cerr << std::endl
3632 *   << std::endl
3633 *   << "----------------------------------------------------"
3634 *   << std::endl;
3635 *   std::cerr << "Unknown exception!" << std::endl
3636 *   << "Aborting!" << std::endl
3637 *   << "----------------------------------------------------"
3638 *   << std::endl;
3639 *   return 1;
3640 *   }
3641 *  
3642 *   return 0;
3643 *   }
3644 * @endcode
3645<a name="step_31-Results"></a><h1>Results</h1>
3646
3647
3648<a name="step_31-Resultsin2d"></a><h3> Results in 2d </h3>
3649
3650
3651When you run the program in 2d, the output will look something like
3652this:
3653<code>
3654<pre>
3655Number of active cells: 256 (on 5 levels)
3656Number of degrees of freedom: 3556 (2178+289+1089)
3657
3658Timestep 0: t=0
3659 Assembling...
3660 Rebuilding Stokes preconditioner...
3661 Solving...
3662 0 GMRES iterations for Stokes subsystem.
3663 Time step: 0.919118
3664 9 CG iterations for temperature.
3665 Temperature range: -0.16687 1.30011
3666
3667Number of active cells: 280 (on 6 levels)
3668Number of degrees of freedom: 4062 (2490+327+1245)
3669
3670Timestep 0: t=0
3671 Assembling...
3672 Rebuilding Stokes preconditioner...
3673 Solving...
3674 0 GMRES iterations for Stokes subsystem.
3675 Time step: 0.459559
3676 9 CG iterations for temperature.
3677 Temperature range: -0.0982971 0.598503
3678
3679Number of active cells: 520 (on 7 levels)
3680Number of degrees of freedom: 7432 (4562+589+2281)
3681
3682Timestep 0: t=0
3683 Assembling...
3684 Rebuilding Stokes preconditioner...
3685 Solving...
3686 0 GMRES iterations for Stokes subsystem.
3687 Time step: 0.229779
3688 9 CG iterations for temperature.
3689 Temperature range: -0.0551098 0.294493
3690
3691Number of active cells: 1072 (on 8 levels)
3692Number of degrees of freedom: 15294 (9398+1197+4699)
3693
3694Timestep 0: t=0
3695 Assembling...
3696 Rebuilding Stokes preconditioner...
3697 Solving...
3698 0 GMRES iterations for Stokes subsystem.
3699 Time step: 0.11489
3700 9 CG iterations for temperature.
3701 Temperature range: -0.0273524 0.156861
3702
3703Number of active cells: 2116 (on 9 levels)
3704Number of degrees of freedom: 30114 (18518+2337+9259)
3705
3706Timestep 0: t=0
3707 Assembling...
3708 Rebuilding Stokes preconditioner...
3709 Solving...
3710 0 GMRES iterations for Stokes subsystem.
3711 Time step: 0.0574449
3712 9 CG iterations for temperature.
3713 Temperature range: -0.014993 0.0738328
3714
3715Timestep 1: t=0.0574449
3716 Assembling...
3717 Solving...
3718 56 GMRES iterations for Stokes subsystem.
3719 Time step: 0.0574449
3720 9 CG iterations for temperature.
3721 Temperature range: -0.0273934 0.14488
3722
3723...
3724</pre>
3725</code>
3726
3727In the beginning we refine the mesh several times adaptively and
3728always return to time step zero to restart on the newly refined
3729mesh. Only then do we start the actual time iteration.
3730
3731The program runs for a while. The temperature field for time steps 0,
3732500, 1000, 1500, 2000, 3000, 4000, and 5000 looks like this (note that
3733the color scale used for the temperature is not always the same):
3734
3735<table align="center" class="doxtable">
3736 <tr>
3737 <td>
3738 <img src="https://www.dealii.org/images/steps/developer/step-31.2d.solution.00.png" alt="">
3739 </td>
3740 <td>
3741 <img src="https://www.dealii.org/images/steps/developer/step-31.2d.solution.01.png" alt="">
3742 </td>
3743 <td>
3744 <img src="https://www.dealii.org/images/steps/developer/step-31.2d.solution.02.png" alt="">
3745 </td>
3746 <td>
3747 <img src="https://www.dealii.org/images/steps/developer/step-31.2d.solution.03.png" alt="">
3748 </td>
3749 </tr>
3750 <tr>
3751 <td>
3752 <img src="https://www.dealii.org/images/steps/developer/step-31.2d.solution.04.png" alt="">
3753 </td>
3754 <td>
3755 <img src="https://www.dealii.org/images/steps/developer/step-31.2d.solution.05.png" alt="">
3756 </td>
3757 <td>
3758 <img src="https://www.dealii.org/images/steps/developer/step-31.2d.solution.06.png" alt="">
3759 </td>
3760 <td>
3761 <img src="https://www.dealii.org/images/steps/developer/step-31.2d.solution.07.png" alt="">
3762 </td>
3763 </tr>
3764</table>
3765
3766The visualizations shown here were generated using a version of the example
3767which did not enforce the constraints after transferring the mesh.
3768
3769As can be seen, we have three heat sources that heat fluid and
3770therefore produce a buoyancy effect that lets hots pockets of fluid
3771rise up and swirl around. By a chimney effect, the three streams are
3772pressed together by fluid that comes from the outside and wants to
3773join the updraft party. Note that because the fluid is initially at
3774rest, those parts of the fluid that were initially over the sources
3775receive a longer heating time than that fluid that is later dragged
3776over the source by the fully developed flow field. It is therefore
3777hotter, a fact that can be seen in the red tips of the three
3778plumes. Note also the relatively fine features of the flow field, a
3779result of the sophisticated transport stabilization of the temperature
3780equation we have chosen.
3781
3782In addition to the pictures above, the following ones show the
3783adaptive mesh and the flow field at the same time steps:
3784
3785<table align="center" class="doxtable">
3786 <tr>
3787 <td>
3788 <img src="https://www.dealii.org/images/steps/developer/step-31.2d.grid.00.png" alt="">
3789 </td>
3790 <td>
3791 <img src="https://www.dealii.org/images/steps/developer/step-31.2d.grid.01.png" alt="">
3792 </td>
3793 <td>
3794 <img src="https://www.dealii.org/images/steps/developer/step-31.2d.grid.02.png" alt="">
3795 </td>
3796 <td>
3797 <img src="https://www.dealii.org/images/steps/developer/step-31.2d.grid.03.png" alt="">
3798 </td>
3799 </tr>
3800 <tr>
3801 <td>
3802 <img src="https://www.dealii.org/images/steps/developer/step-31.2d.grid.04.png" alt="">
3803 </td>
3804 <td>
3805 <img src="https://www.dealii.org/images/steps/developer/step-31.2d.grid.05.png" alt="">
3806 </td>
3807 <td>
3808 <img src="https://www.dealii.org/images/steps/developer/step-31.2d.grid.06.png" alt="">
3809 </td>
3810 <td>
3811 <img src="https://www.dealii.org/images/steps/developer/step-31.2d.grid.07.png" alt="">
3812 </td>
3813 </tr>
3814</table>
3815
3816
3817<a name="step_31-Resultsin3d"></a><h3> Results in 3d </h3>
3818
3819
3820The same thing can of course be done in 3d by changing the template
3821parameter to the BoussinesqFlowProblem object in <code>main()</code>
3822from 2 to 3, so that the output now looks like follows:
3823
3824<code>
3825<pre>
3826Number of active cells: 64 (on 3 levels)
3827Number of degrees of freedom: 3041 (2187+125+729)
3828
3829Timestep 0: t=0
3830 Assembling...
3831 Rebuilding Stokes preconditioner...
3832 Solving...
3833 0 GMRES iterations for Stokes subsystem.
3834 Time step: 2.45098
3835 9 CG iterations for temperature.
3836 Temperature range: -0.675683 4.94725
3837
3838Number of active cells: 288 (on 4 levels)
3839Number of degrees of freedom: 12379 (8943+455+2981)
3840
3841Timestep 0: t=0
3842 Assembling...
3843 Rebuilding Stokes preconditioner...
3844 Solving...
3845 0 GMRES iterations for Stokes subsystem.
3846 Time step: 1.22549
3847 9 CG iterations for temperature.
3848 Temperature range: -0.527701 2.25764
3849
3850Number of active cells: 1296 (on 5 levels)
3851Number of degrees of freedom: 51497 (37305+1757+12435)
3852
3853Timestep 0: t=0
3854 Assembling...
3855 Rebuilding Stokes preconditioner...
3856 Solving...
3857 0 GMRES iterations for Stokes subsystem.
3858 Time step: 0.612745
3859 10 CG iterations for temperature.
3860 Temperature range: -0.496942 0.847395
3861
3862Number of active cells: 5048 (on 6 levels)
3863Number of degrees of freedom: 192425 (139569+6333+46523)
3864
3865Timestep 0: t=0
3866 Assembling...
3867 Rebuilding Stokes preconditioner...
3868 Solving...
3869 0 GMRES iterations for Stokes subsystem.
3870 Time step: 0.306373
3871 10 CG iterations for temperature.
3872 Temperature range: -0.267683 0.497739
3873
3874Timestep 1: t=0.306373
3875 Assembling...
3876 Solving...
3877 27 GMRES iterations for Stokes subsystem.
3878 Time step: 0.306373
3879 10 CG iterations for temperature.
3880 Temperature range: -0.461787 0.958679
3881
3882...
3883</pre>
3884</code>
3885
3886Visualizing the temperature isocontours at time steps 0,
388750, 100, 150, 200, 300, 400, 500, 600, 700, and 800 yields the
3888following plots:
3889
3890<table align="center" class="doxtable">
3891 <tr>
3892 <td>
3893 <img src="https://www.dealii.org/images/steps/developer/step-31.3d.solution.00.png" alt="">
3894 </td>
3895 <td>
3896 <img src="https://www.dealii.org/images/steps/developer/step-31.3d.solution.01.png" alt="">
3897 </td>
3898 <td>
3899 <img src="https://www.dealii.org/images/steps/developer/step-31.3d.solution.02.png" alt="">
3900 </td>
3901 <td>
3902 <img src="https://www.dealii.org/images/steps/developer/step-31.3d.solution.03.png" alt="">
3903 </td>
3904 </tr>
3905 <tr>
3906 <td>
3907 <img src="https://www.dealii.org/images/steps/developer/step-31.3d.solution.04.png" alt="">
3908 </td>
3909 <td>
3910 <img src="https://www.dealii.org/images/steps/developer/step-31.3d.solution.05.png" alt="">
3911 </td>
3912 <td>
3913 <img src="https://www.dealii.org/images/steps/developer/step-31.3d.solution.06.png" alt="">
3914 </td>
3915 <td>
3916 <img src="https://www.dealii.org/images/steps/developer/step-31.3d.solution.07.png" alt="">
3917 </td>
3918 </tr>
3919 <tr>
3920 <td>
3921 <img src="https://www.dealii.org/images/steps/developer/step-31.3d.solution.08.png" alt="">
3922 </td>
3923 <td>
3924 <img src="https://www.dealii.org/images/steps/developer/step-31.3d.solution.09.png" alt="">
3925 </td>
3926 <td>
3927 <img src="https://www.dealii.org/images/steps/developer/step-31.3d.solution.10.png" alt="">
3928 </td>
3929 <td>
3930 </td>
3931 </tr>
3932</table>
3933
3934That the first picture looks like three hedgehogs stems from the fact that our
3935scheme essentially projects the source times the first time step size onto the
3936mesh to obtain the temperature field in the first time step. Since the source
3937function is discontinuous, we need to expect over- and undershoots from this
3938project. This is in fact what happens (it's easier to check this in 2d) and
3939leads to the crumpled appearance of the isosurfaces. The visualizations shown
3940here were generated using a version of the example which did not enforce the
3941constraints after transferring the mesh.
3942
3943
3944
3945<a name="step_31-Numericalexperimentstodetermineoptimalparameters"></a><h3> Numerical experiments to determine optimal parameters </h3>
3946
3947
3948The program as is has three parameters that we don't have much of a
3949theoretical handle on how to choose in an optimal way. These are:
3950<ul>
3951 <li>The time step must satisfy a CFL condition
3952 @f$k\le \min_K \frac{c_kh_K}{\|\mathbf{u}\|_{L^\infty(K)}}@f$. Here, @f$c_k@f$ is
3953 dimensionless, but what is the right value?
3954 <li>In the computation of the artificial viscosity,
3955@f{eqnarray*}{
3956 \nu_\alpha(T)|_K
3957 =
3958 \beta
3959 \|\mathbf{u}\|_{L^\infty(K)}
3960 \min\left\{
3961 h_K,
3962 h_K^\alpha
3963 \frac{\|R_\alpha(T)\|_{L^\infty(K)}}{c(\mathbf{u},T)}
3964 \right\},
3965@f}
3966 with @f$c(\mathbf{u},T) =
3967 c_R\ \|\mathbf{u}\|_{L^\infty(\Omega)} \ \mathrm{var}(T)
3968 \ |\mathrm{diam}(\Omega)|^{\alpha-2}@f$.
3969 Here, the choice of the dimensionless %numbers @f$\beta,c_R@f$ is of
3970 interest.
3971</ul>
3972In all of these cases, we will have to expect that the correct choice of each
3973value depends on that of the others, and most likely also on the space
3974dimension and polynomial degree of the finite element used for the
3975temperature. Below we'll discuss a few numerical experiments to choose
3976constants @f$c_k@f$ and @f$\beta@f$.
3977
3978Below, we will not discuss the choice of @f$c_R@f$. In the program, we set
3979it to @f$c_R=2^{\frac{4-2\alpha}{d}}@f$. The reason for this value is a
3980bit complicated and has more to do with the history of the program
3981than reasoning: while the correct formula for the global scaling
3982parameter @f$c(\mathbf{u},T)@f$ is shown above, the program (including the
3983version shipped with deal.II 6.2) initially had a bug in that we
3984computed
3985@f$c(\mathbf{u},T) =
3986 \|\mathbf{u}\|_{L^\infty(\Omega)} \ \mathrm{var}(T)
3987 \ \frac{1}{|\mathrm{diam}(\Omega)|^{\alpha-2}}@f$ instead, where
3988we had set the scaling parameter to one. Since we only computed on the
3989unit square/cube where @f$\mathrm{diam}(\Omega)=2^{1/d}@f$, this was
3990entirely equivalent to using the correct formula with
3991@f$c_R=\left(2^{1/d}\right)^{4-2\alpha}=2^{\frac{4-2\alpha}{d}}@f$. Since
3992this value for @f$c_R@f$ appears to work just fine for the current
3993program, we corrected the formula in the program and set @f$c_R@f$ to a
3994value that reproduces exactly the results we had before. We will,
3995however, revisit this issue again in @ref step_32 "step-32".
3996
3997Now, however, back to the discussion of what values of @f$c_k@f$ and
3998@f$\beta@f$ to choose:
3999
4000
4001<a name="step_31-Choosingicsubksubiandbeta"></a><h4> Choosing <i>c<sub>k</sub></i> and beta </h4>
4002
4003
4004These two constants are definitely linked in some way. The reason is easy to
4005see: In the case of a pure advection problem,
4006@f$\frac{\partial T}{\partial t} + \mathbf{u}\cdot\nabla T = \gamma@f$, any
4007explicit scheme has to satisfy a CFL condition of the form
4008@f$k\le \min_K \frac{c_k^a h_K}{\|\mathbf{u}\|_{L^\infty(K)}}@f$. On the other hand,
4009for a pure diffusion problem,
4010@f$\frac{\partial T}{\partial t} + \nu \Delta T = \gamma@f$,
4011explicit schemes need to satisfy a condition
4012@f$k\le \min_K \frac{c_k^d h_K^2}{\nu}@f$. So given the form of @f$\nu@f$ above, an
4013advection diffusion problem like the one we have to solve here will result in
4014a condition of the form
4015@f$
4016k\le \min_K \min \left\{
4017 \frac{c_k^a h_K}{\|\mathbf{u}\|_{L^\infty(K)}},
4018 \frac{c_k^d h_K^2}{\beta \|\mathbf{u}\|_{L^\infty(K)} h_K}\right\}
4019 =
4020 \min_K \left( \min \left\{
4021 c_k^a,
4022 \frac{c_k^d}{\beta}\right\}
4023 \frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}} \right)
4024@f$.
4025It follows that we have to face the fact that we might want to choose @f$\beta@f$
4026larger to improve the stability of the numerical scheme (by increasing the
4027amount of artificial diffusion), but we have to pay a price in the form of
4028smaller, and consequently more time steps. In practice, one would therefore
4029like to choose @f$\beta@f$ as small as possible to keep the transport problem
4030sufficiently stabilized while at the same time trying to choose the time step
4031as large as possible to reduce the overall amount of work.
4032
4033The find the right balance, the only way is to do a few computational
4034experiments. Here's what we did: We modified the program slightly to allow
4035less mesh refinement (so we don't always have to wait that long) and to choose
4036@f$
4037 \nu(T)|_K
4038 =
4039 \beta
4040 \|\mathbf{u}\|_{L^\infty(K)} h_K
4041@f$ to eliminate the effect of the constant @f$c_R@f$ (we know that
4042solutions are stable by using this version of @f$\nu(T)@f$ as an artificial
4043viscosity, but that we can improve things -- i.e. make the solution
4044sharper -- by using the more complicated formula for this artificial
4045viscosity). We then run the program
4046for different values @f$c_k,\beta@f$ and observe maximal and minimal temperatures
4047in the domain. What we expect to see is this: If we choose the time step too
4048big (i.e. choose a @f$c_k@f$ bigger than theoretically allowed) then we will get
4049exponential growth of the temperature. If we choose @f$\beta@f$ too small, then
4050the transport stabilization becomes insufficient and the solution will show
4051significant oscillations but not exponential growth.
4052
4053
4054<a name="step_31-ResultsforQsub1subelements"></a><h5>Results for Q<sub>1</sub> elements</h5>
4055
4056
4057Here is what we get for
4058@f$\beta=0.01, \beta=0.1@f$, and @f$\beta=0.5@f$, different choices of @f$c_k@f$, and
4059bilinear elements (<code>temperature_degree=1</code>) in 2d:
4060
4061<table align="center" class="doxtable">
4062 <tr>
4063 <td>
4064 <img src="https://www.dealii.org/images/steps/developer/step-31.timestep.q1.beta=0.01.png" alt="">
4065 </td>
4066 <td>
4067 <img src="https://www.dealii.org/images/steps/developer/step-31.timestep.q1.beta=0.03.png" alt="">
4068 </td>
4069 </tr>
4070 <tr>
4071 <td>
4072 <img src="https://www.dealii.org/images/steps/developer/step-31.timestep.q1.beta=0.1.png" alt="">
4073 </td>
4074 <td>
4075 <img src="https://www.dealii.org/images/steps/developer/step-31.timestep.q1.beta=0.5.png" alt="">
4076 </td>
4077 </tr>
4078</table>
4079
4080The way to interpret these graphs goes like this: for @f$\beta=0.01@f$ and
4081@f$c_k=\frac 12,\frac 14@f$, we see exponential growth or at least large
4082variations, but if we choose
4083@f$k=\frac 18\frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}@f$
4084or smaller, then the scheme is
4085stable though a bit wobbly. For more artificial diffusion, we can choose
4086@f$k=\frac 14\frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}@f$
4087or smaller for @f$\beta=0.03@f$,
4088@f$k=\frac 13\frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}@f$
4089or smaller for @f$\beta=0.1@f$, and again need
4090@f$k=\frac 1{15}\frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}@f$
4091for @f$\beta=0.5@f$ (this time because much diffusion requires a small time
4092step).
4093
4094So how to choose? If we were simply interested in a large time step, then we
4095would go with @f$\beta=0.1@f$ and
4096@f$k=\frac 13\frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}@f$.
4097On the other hand, we're also interested in accuracy and here it may be of
4098interest to actually investigate what these curves show. To this end note that
4099we start with a zero temperature and that our sources are positive &mdash; so
4100we would intuitively expect that the temperature can never drop below
4101zero. But it does, a consequence of Gibb's phenomenon when using continuous
4102elements to approximate a discontinuous solution. We can therefore see that
4103choosing @f$\beta@f$ too small is bad: too little artificial diffusion leads to
4104over- and undershoots that aren't diffused away. On the other hand, for large
4105@f$\beta@f$, the minimum temperature drops below zero at the beginning but then
4106quickly diffuses back to zero.
4107
4108On the other hand, let's also look at the maximum temperature. Watching the
4109movie of the solution, we see that initially the fluid is at rest. The source
4110keeps heating the same volume of fluid whose temperature increases linearly at
4111the beginning until its buoyancy is able to move it upwards. The hottest part
4112of the fluid is therefore transported away from the solution and fluid taking
4113its place is heated for only a short time before being moved out of the source
4114region, therefore remaining cooler than the initial bubble. If @f$\kappa=0@f$
4115(in the program it is nonzero but very small) then the hottest part of the
4116fluid should be advected along with the flow with its temperature
4117constant. That's what we can see in the graphs with the smallest @f$\beta@f$: Once
4118the maximum temperature is reached, it hardly changes any more. On the other
4119hand, the larger the artificial diffusion, the more the hot spot is
4120diffused. Note that for this criterion, the time step size does not play a
4121significant role.
4122
4123So to sum up, likely the best choice would appear to be @f$\beta=0.03@f$
4124and @f$k=\frac 14\frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}@f$. The curve is
4125a bit wobbly, but overall pictures looks pretty reasonable with the
4126exception of some over and undershoots close to the start time due to
4127Gibb's phenomenon.
4128
4129
4130<a name="step_31-ResultsforQsub2subelements"></a><h5>Results for Q<sub>2</sub> elements</h5>
4131
4132
4133One can repeat the same sequence of experiments for higher order
4134elements as well. Here are the graphs for bi-quadratic shape functions
4135(<code>temperature_degree=2</code>) for the temperature, while we
4136retain the @f$Q_2/Q_1@f$ stable Taylor-Hood element for the Stokes system:
4137
4138<table align="center" class="doxtable">
4139 <tr>
4140 <td>
4141 <img src="https://www.dealii.org/images/steps/developer/step-31.timestep.q2.beta=0.01.png" alt="">
4142 </td>
4143 <td>
4144 <img src="https://www.dealii.org/images/steps/developer/step-31.timestep.q2.beta=0.03.png" alt="">
4145 </td>
4146 </tr>
4147 <tr>
4148 <td>
4149 <img src="https://www.dealii.org/images/steps/developer/step-31.timestep.q2.beta=0.1.png" alt="">
4150 </td>
4151 </tr>
4152</table>
4153
4154Again, small values of @f$\beta@f$ lead to less diffusion but we have to
4155choose the time step very small to keep things under control. Too
4156large values of @f$\beta@f$ make for more diffusion, but again require
4157small time steps. The best value would appear to be @f$\beta=0.03@f$, as
4158for the @f$Q_1@f$ element, and then we have to choose
4159@f$k=\frac 18\frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}@f$ &mdash; exactly
4160half the size for the @f$Q_1@f$ element, a fact that may not be surprising
4161if we state the CFL condition as the requirement that the time step be
4162small enough so that the distance transport advects in each time step
4163is no longer than one <i>grid point</i> away (which for @f$Q_1@f$ elements
4164is @f$h_K@f$, but for @f$Q_2@f$ elements is @f$h_K/2@f$). It turns out that @f$\beta@f$
4165needs to be slightly larger for obtaining stable results also late in
4166the simulation at times larger than 60, so we actually choose it as
4167@f$\beta = 0.034@f$ in the code.
4168
4169
4170<a name="step_31-Resultsfor3d"></a><h5>Results for 3d</h5>
4171
4172
4173One can repeat these experiments in 3d and find the optimal time step
4174for each value of @f$\beta@f$ and find the best value of @f$\beta@f$. What one
4175finds is that for the same @f$\beta@f$ already used in 2d, the time steps
4176needs to be a bit smaller, by around a factor of 1.2 or so. This is
4177easily explained: the time step restriction is
4178@f$k=\min_K \frac{ch_K}{\|\mathbf{u}\|_{L^\infty(K)}}@f$ where @f$h_K@f$ is
4179the <i>diameter</i> of the cell. However, what is really needed is the
4180distance between mesh points, which is @f$\frac{h_K}{\sqrt{d}}@f$. So a
4181more appropriate form would be
4182@f$k=\min_K \frac{ch_K}{\|\mathbf{u}\|_{L^\infty(K)}\sqrt{d}}@f$.
4183
4184The second find is that one needs to choose @f$\beta@f$ slightly bigger
4185(about @f$\beta=0.05@f$ or so). This then again reduces the time step we
4186can take.
4187
4188
4189
4190
4191<a name="step_31-Conclusions"></a><h5>Conclusions</h5>
4192
4193
4194Concluding, from the simple computations above, @f$\beta=0.034@f$ appears to be a
4195good choice for the stabilization parameter in 2d, and @f$\beta=0.05@f$ in 3d. In
4196a dimension independent way, we can model this as @f$\beta=0.017d@f$. If one does
4197longer computations (several thousand time steps) on finer meshes, one
4198realizes that the time step size is not quite small enough and that for
4199stability one will have to reduce the above values a bit more (by about a
4200factor of @f$\frac 78@f$).
4201
4202As a consequence, a formula that reconciles 2d, 3d, and variable polynomial
4203degree and takes all factors in account reads as follows:
4204@f{eqnarray*}{
4205 k =
4206 \frac 1{2 \cdot 1.7} \frac 1{\sqrt{d}}
4207 \frac 2d
4208 \frac 1{q_T}
4209 \frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}
4210 =
4211 \frac 1{1.7 d\sqrt{d}}
4212 \frac 1{q_T}
4213 \frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}.
4214@f}
4215In the first form (in the center of the equation), @f$\frac
42161{2 \cdot 1.7}@f$ is a universal constant, @f$\frac 1{\sqrt{d}}@f$
4217is the factor that accounts for the difference between cell diameter
4218and grid point separation,
4219@f$\frac 2d@f$ accounts for the increase in @f$\beta@f$ with space dimension,
4220@f$\frac 1{q_T}@f$ accounts for the distance between grid points for
4221higher order elements, and @f$\frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}@f$
4222for the local speed of transport relative to the cell size. This is
4223the formula that we use in the program.
4224
4225As for the question of whether to use @f$Q_1@f$ or @f$Q_2@f$ elements for the
4226temperature, the following considerations may be useful: First,
4227solving the temperature equation is hardly a factor in the overall
4228scheme since almost the entire compute time goes into solving the
4229Stokes system in each time step. Higher order elements for the
4230temperature equation are therefore not a significant drawback. On the
4231other hand, if one compares the size of the over- and undershoots the
4232solution produces due to the discontinuous source description, one
4233notices that for the choice of @f$\beta@f$ and @f$k@f$ as above, the @f$Q_1@f$
4234solution dips down to around @f$-0.47@f$, whereas the @f$Q_2@f$ solution only
4235goes to @f$-0.13@f$ (remember that the exact solution should never become
4236negative at all. This means that the @f$Q_2@f$ solution is significantly
4237more accurate; the program therefore uses these higher order elements,
4238despite the penalty we pay in terms of smaller time steps.
4239
4240
4241<a name="step_31-Possibilitiesforextensions"></a><h3> Possibilities for extensions </h3>
4242
4243
4244There are various ways to extend the current program. Of particular interest
4245is, of course, to make it faster and/or increase the resolution of the
4246program, in particular in 3d. This is the topic of the @ref step_32 "step-32"
4247tutorial program which will implement strategies to solve this problem in
4248%parallel on a cluster. It is also the basis of the much larger open
4249source code ASPECT (see https://aspect.geodynamics.org/ ) that can solve realistic
4250problems and that constitutes the further development of @ref step_32 "step-32".
4251
4252Another direction would be to make the fluid flow more realistic. The program
4253was initially written to simulate various cases simulating the convection of
4254material in the earth's mantle, i.e. the zone between the outer earth core and
4255the solid earth crust: there, material is heated from below and cooled from
4256above, leading to thermal convection. The physics of this fluid are much more
4257complicated than shown in this program, however: The viscosity of mantle
4258material is strongly dependent on the temperature, i.e. @f$\eta=\eta(T)@f$, with
4259the dependency frequently modeled as a viscosity that is reduced exponentially
4260with rising temperature. Secondly, much of the dynamics of the mantle is
4261determined by chemical reactions, primarily phase changes of the various
4262crystals that make up the mantle; the buoyancy term on the right hand side of
4263the Stokes equations then depends not only on the temperature, but also on the
4264chemical composition at a given location which is advected by the flow field
4265but also changes as a function of pressure and temperature. We will
4266investigate some of these effects in later tutorial programs as well.
4267 *
4268 *
4269<a name="step_31-PlainProg"></a>
4270<h1> The plain program</h1>
4271@include "step-31.cc"
4272*/
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, level_dof_access > > &cell)
Definition fe_q.h:550
const std::vector< Point< dim > > & get_unit_support_points() const
const unsigned int n_components
Definition function.h:163
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
virtual void vector_value(const Point< dim > &p, Vector< RangeNumberType > &values) const
static void estimate(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Quadrature< dim - 1 > &quadrature, const std::map< types::boundary_id, const Function< spacedim, Number > * > &neumann_bc, const ReadVector< Number > &solution, Vector< float > &error, const ComponentMask &component_mask={}, const Function< spacedim > *coefficients=nullptr, const unsigned int n_threads=numbers::invalid_unsigned_int, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id, const types::material_id material_id=numbers::invalid_material_id, const Strategy strategy=cell_diameter_over_24)
Definition point.h:111
float depth
Point< 3 > vertices[4]
Point< 2 > second
Definition grid_out.cc:4614
Point< 2 > first
Definition grid_out.cc:4613
unsigned int level
Definition grid_out.cc:4616
__global__ void set(Number *val, const Number s, const size_type N)
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcMessage(std::string arg1)
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(DOFINFO &, DOFINFO &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
Definition loop.h:442
@ update_hessians
Second derivatives of shape functions.
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
IndexSet complete_index_set(const IndexSet::size_type N)
Definition index_set.h:1193
const Event initial
Definition event.cc:64
void approximate(const SynchronousIterators< std::tuple< typename DoFHandler< dim, spacedim >::active_cell_iterator, Vector< float >::iterator > > &cell, const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof_handler, const InputVector &solution, const unsigned int component)
Expression sign(const Expression &x)
void downstream(DoFHandler< dim, spacedim > &dof_handler, const Tensor< 1, spacedim > &direction, const bool dof_wise_renumbering=false)
void interpolate(const DoFHandler< dim, spacedim > &dof1, const InVector &u1, const DoFHandler< dim, spacedim > &dof2, OutVector &u2)
void refine(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold, const unsigned int max_to_mark=numbers::invalid_unsigned_int)
void refine_and_coarsen_fixed_fraction(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double top_fraction, const double bottom_fraction, const unsigned int max_n_cells=std::numeric_limits< unsigned int >::max(), const VectorTools::NormType norm_type=VectorTools::L1_norm)
const std::vector< bool > & used
spacedim & mesh
Definition grid_tools.h:980
double diameter(const Triangulation< dim, spacedim > &tria)
if(marked_vertices.size() !=0) for(auto it
for(unsigned int j=best_vertex+1;j< vertices.size();++j) if(vertices_to_use[j])
@ matrix
Contents is actually a matrix.
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
Definition divergence.h:471
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)
constexpr ReturnType< rank, T >::value_type & extract(T &t, const ArrayType &indices)
spacedim const DoFHandler< dim, spacedim > const FullMatrix< double > & transfer
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
unsigned int n_active_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition tria.cc:15058
void copy(const T *begin, const T *end, U *dest)
int(&) functions(const void *v1, const void *v2)
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
Definition loop.h:70
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
const InputIterator OutputIterator out
Definition parallel.h:167
const Iterator const std_cxx20::type_identity_t< Iterator > & end
Definition parallel.h:610
const InputIterator OutputIterator const Function & function
Definition parallel.h:168
STL namespace.
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation