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