1269 *
const unsigned int = 0)
const override
1274 *
virtual void vector_value(
const Point<dim> &p,
1277 *
for (
unsigned int c = 0; c < this->n_components; ++c)
1278 *
value(c) = TemperatureInitialValues<dim>::value(p, c);
1284 *
template <
int dim>
1285 *
class TemperatureRightHandSide :
public Function<dim>
1288 *
TemperatureRightHandSide()
1293 *
const unsigned int component = 0) const override
1297 *
ExcMessage(
"Invalid operation for a scalar function."));
1301 *
static const Point<dim> source_centers[3] = {
1303 *
(dim == 2 ?
Point<dim>(.45, .1) :
Point<dim>(.45, .5, .1)),
1304 *
(dim == 2 ?
Point<dim>(.75, .1) :
Point<dim>(.75, .5, .1))};
1305 *
static const double source_radius = (dim == 2 ? 1. / 32 : 1. / 8);
1307 *
return ((source_centers[0].distance(p) < source_radius) ||
1308 *
(source_centers[1].distance(p) < source_radius) ||
1309 *
(source_centers[2].distance(p) < source_radius) ?
1317 *
for (
unsigned int c = 0; c < this->
n_components; ++c)
1318 *
value(c) = TemperatureRightHandSide<dim>::value(p, c);
1328 * <a name=
"step_31-Linearsolversandpreconditioners"></a>
1329 * <h3>Linear solvers and preconditioners</h3>
1333 * This section introduces some objects that are used
for the solution of
1334 * the linear equations of the Stokes system that we need to solve in each
1335 * time step. Many of the ideas used here are the same as in @ref step_20
"step-20", where
1336 * Schur complement based preconditioners and solvers have been introduced,
1337 * with the actual
interface taken from @ref step_22
"step-22" (in particular the
1338 * discussion in the
"Results" section of @ref step_22
"step-22", in which we introduce
1339 * alternatives to the direct Schur complement approach). Note, however,
1340 * that here we don
't use the Schur complement to solve the Stokes
1341 * equations, though an approximate Schur complement (the mass matrix on the
1342 * pressure space) appears in the preconditioner.
1345 * namespace LinearSolvers
1350 * <a name="step_31-ThecodeInverseMatrixcodeclasstemplate"></a>
1351 * <h4>The <code>InverseMatrix</code> class template</h4>
1355 * This class is an interface to calculate the action of an "inverted"
1356 * matrix on a vector (using the <code>vmult</code> operation) in the same
1357 * way as the corresponding class in @ref step_22 "step-22": when the product of an
1358 * object of this class is requested, we solve a linear equation system
1359 * with that matrix using the CG method, accelerated by a preconditioner
1360 * of (templated) class <code>PreconditionerType</code>.
1364 * In a minor deviation from the implementation of the same class in
1365 * @ref step_22 "step-22", we make the <code>vmult</code> function take any
1366 * kind of vector type (it will yield compiler errors, however, if the
1367 * matrix does not allow a matrix-vector product with this kind of
1372 * Secondly, we catch any exceptions that the solver may have thrown. The
1373 * reason is as follows: When debugging a program like this one
1374 * occasionally makes a mistake of passing an indefinite or nonsymmetric
1375 * matrix or preconditioner to the current class. The solver will, in that
1376 * case, not converge and throw a run-time exception. If not caught here
1377 * it will propagate up the call stack and may end up in
1378 * <code>main()</code> where we output an error message that will say that
1379 * the CG solver failed. The question then becomes: Which CG solver? The
1380 * one that inverted the mass matrix? The one that inverted the top left
1381 * block with the Laplace operator? Or a CG solver in one of the several
1382 * other nested places where we use linear solvers in the current code? No
1383 * indication about this is present in a run-time exception because it
1384 * doesn't store the stack of calls through which we got to the place
1385 * where the exception was generated.
1389 * So rather than letting the exception propagate freely up to
1390 * <code>
main()</code>, we acknowledge that in the current context
1391 * there is little that an outer function can
do if the inner
1392 * solver fails. As a consequence, we deal with the situation by
1393 * catching the exception and letting the program fail by
1394 * triggering an assertion with a `
false` condition (which of
1395 * course always fails) and that uses the error message associated
1396 * with the caught exception (returned by calling `
e.what()`) as
1397 * error text. In other words, instead of letting the error
1398 * message be produced by `
main()`, we rather report it here where
1399 * we can abort the program preserving information about where the
1403 *
template <class MatrixType, class PreconditionerType>
1407 *
InverseMatrix(
const MatrixType &m,
1408 *
const PreconditionerType &preconditioner);
1411 *
template <
typename VectorType>
1412 *
void vmult(VectorType &dst,
const VectorType &src)
const;
1416 *
const PreconditionerType &preconditioner;
1420 *
template <
class MatrixType,
class PreconditionerType>
1421 *
InverseMatrix<MatrixType, PreconditionerType>::InverseMatrix(
1422 *
const MatrixType &m,
1423 *
const PreconditionerType &preconditioner)
1425 *
, preconditioner(preconditioner)
1430 *
template <
class MatrixType,
class PreconditionerType>
1431 *
template <
typename VectorType>
1432 *
void InverseMatrix<MatrixType, PreconditionerType>::vmult(
1434 *
const VectorType &src)
const
1436 *
SolverControl solver_control(src.size(), 1e-7 * src.l2_norm());
1443 *
cg.solve(*matrix, dst, src, preconditioner);
1445 *
catch (std::exception &e)
1447 *
Assert(
false, ExcMessage(
e.what()));
1454 * <a name=
"step_31-Schurcomplementpreconditioner"></a>
1455 * <h4>Schur complement preconditioner</h4>
1459 * This is the implementation of the Schur complement preconditioner as
1460 * described in detail in the introduction. As opposed to @ref step_20
"step-20" and
1461 * @ref step_22
"step-22", we solve the block system all-at-once
using GMRES, and use the
1462 * Schur complement of the block structured
matrix to build a good
1463 * preconditioner instead.
1467 * Let
's have a look at the ideal preconditioner matrix
1468 * @f$P=\left(\begin{array}{cc} A & 0 \\ B & -S \end{array}\right)@f$
1469 * described in the introduction. If we apply this matrix in the solution
1470 * of a linear system, convergence of an iterative GMRES solver will be
1471 * governed by the matrix @f{eqnarray*} P^{-1}\left(\begin{array}{cc} A &
1472 * B^T \\ B & 0 \end{array}\right) = \left(\begin{array}{cc} I & A^{-1}
1473 * B^T \\ 0 & I \end{array}\right), @f} which indeed is very simple. A
1474 * GMRES solver based on exact matrices would converge in one iteration,
1475 * since all eigenvalues are equal (any Krylov method takes at most as
1476 * many iterations as there are distinct eigenvalues). Such a
1477 * preconditioner for the blocked Stokes system has been proposed by
1478 * Silvester and Wathen ("Fast iterative solution of stabilised Stokes
1479 * systems part II. Using general block preconditioners", SIAM
1480 * J. Numer. Anal., 31 (1994), pp. 1352-1367).
1484 * Replacing @f$P@f$ by @f$\tilde{P}@f$ keeps that spirit alive: the product
1485 * @f$P^{-1} A@f$ will still be close to a matrix with eigenvalues 1 with a
1486 * distribution that does not depend on the problem size. This lets us
1487 * hope to be able to get a number of GMRES iterations that is
1488 * problem-size independent.
1492 * The deal.II users who have already gone through the @ref step_20 "step-20" and @ref step_22 "step-22"
1493 * tutorials can certainly imagine how we're going to implement this. We
1494 * replace the exact inverse matrices in @f$P^{-1}@f$ by some
approximate
1495 * inverses built from the InverseMatrix
class, and the inverse Schur
1496 * complement will be approximated by the pressure mass
matrix @f$M_p@f$
1497 * (weighted by @f$\eta^{-1}@f$ as mentioned in the introduction). As pointed
1498 * out in the results section of @ref step_22
"step-22", we can replace the exact inverse
1499 * of @f$A@f$ by just the application of a preconditioner, in
this case
1500 * on a vector Laplace matrix as was explained in the introduction. This
1501 * does increase the number of (outer) GMRES iterations, but is still
1502 * significantly cheaper than an exact inverse, which would require
1503 * between 20 and 35 CG iterations for <em>each</em> outer solver step
1504 * (
using the AMG preconditioner).
1508 * Having the above explanations in mind, we define a preconditioner
class
1509 * with a <code>vmult</code> functionality, which is all we need
for the
1510 * interaction with the usual solver functions further below in the
1515 * First the declarations. These are similar to the definition of the
1516 * Schur complement in @ref step_20
"step-20", with the difference that we need some more
1517 * preconditioners in the constructor and that the matrices we use here
1518 * are built upon Trilinos:
1521 *
template <class PreconditionerTypeA, class PreconditionerTypeMp>
1525 *
BlockSchurPreconditioner(
1528 *
PreconditionerTypeMp> &Mpinv,
1529 *
const PreconditionerTypeA &Apreconditioner);
1538 *
PreconditionerTypeMp>>
1540 *
const PreconditionerTypeA &a_preconditioner;
1555 * an
IndexSet where every valid index is part of the
set. Note that this
1556 * program can only be run sequentially and will throw an exception if used
1560 *
template <class PreconditionerTypeA, class PreconditionerTypeMp>
1561 *
BlockSchurPreconditioner<PreconditionerTypeA, PreconditionerTypeMp>::
1562 *
BlockSchurPreconditioner(
1565 *
PreconditionerTypeMp> &Mpinv,
1566 *
const PreconditionerTypeA &Apreconditioner)
1567 *
: stokes_matrix(&S)
1568 *
, m_inverse(&Mpinv)
1569 *
, a_preconditioner(Apreconditioner)
1576 * Next is the <code>vmult</code> function. We implement the action of
1577 * @f$P^{-1}@f$ as described above in three successive steps. In formulas, we
1578 * want to compute @f$Y=P^{-1}X@f$ where @f$X,Y@f$ are both vectors with two block
1583 * The
first step multiplies the velocity part of the vector by a
1584 * preconditioner of the
matrix @f$A@f$, i.e., we compute @f$Y_0={\tilde
1585 *
A}^{-1}X_0@f$. The resulting velocity vector is then multiplied by @f$B@f$
1586 * and subtracted from the pressure, i.e., we want to compute @f$X_1-BY_0@f$.
1587 * This
second step only acts on the pressure vector and is accomplished
1588 * by the residual function of our
matrix classes, except that the
sign is
1589 * wrong. Consequently, we change the
sign in the temporary pressure
1590 * vector and
finally multiply by the inverse pressure mass
matrix to get
1591 * the
final pressure vector, completing our work on the Stokes
1595 *
template <
class PreconditionerTypeA,
class PreconditionerTypeMp>
1597 *
BlockSchurPreconditioner<PreconditionerTypeA, PreconditionerTypeMp>::vmult(
1601 *
a_preconditioner.vmult(dst.block(0), src.block(0));
1602 *
stokes_matrix->block(1, 0).residual(tmp, dst.block(0), src.block(1));
1604 *
m_inverse->vmult(dst.block(1), tmp);
1613 * <a name=
"step_31-ThecodeBoussinesqFlowProblemcodeclasstemplate"></a>
1614 * <h3>The <code>BoussinesqFlowProblem</code>
class template</h3>
1618 * The definition of the
class that defines the top-
level logic of solving
1619 * the time-dependent Boussinesq problem is mainly based on the @ref step_22
"step-22"
1620 * tutorial program. The
main differences are that now we also have to solve
1622 *
object for the
temperature variable as well as matrices, right hand
1623 * sides, and solution vectors
for the current and previous time steps. As
1624 * mentioned in the introduction, all linear algebra objects are going to
1625 * use wrappers of the corresponding Trilinos functionality.
1629 * The member
functions of
this class are reminiscent of @ref step_21
"step-21", where we
1630 * also used a staggered scheme that
first solve the flow equations (here
1631 * the Stokes equations, in @ref step_21
"step-21" Darcy flow) and then update the advected
1632 * quantity (here the temperature, there the saturation). The
functions that
1633 * are
new are mainly concerned with determining the time step, as well as
1634 * the proper
size of the artificial viscosity stabilization.
1638 * The last three variables indicate whether the various matrices or
1639 * preconditioners need to be rebuilt the next time the corresponding build
1640 *
functions are called. This allows us to move the corresponding
1641 * <code>
if</code> into the respective function and thereby keeping our
main
1642 * <code>
run()</code> function clean and easy to read.
1645 *
template <
int dim>
1646 *
class BoussinesqFlowProblem
1649 *
BoussinesqFlowProblem();
1653 *
void setup_dofs();
1654 *
void assemble_stokes_preconditioner();
1655 *
void build_stokes_preconditioner();
1656 *
void assemble_stokes_system();
1657 *
void assemble_temperature_system(
const double maximal_velocity);
1658 *
void assemble_temperature_matrix();
1659 *
double get_maximal_velocity()
const;
1660 *
std::pair<double, double> get_extrapolated_temperature_range()
const;
1662 *
void output_results()
const;
1663 *
void refine_mesh(
const unsigned int max_grid_level);
1665 *
double compute_viscosity(
1666 *
const std::vector<double> &old_temperature,
1667 *
const std::vector<double> &old_old_temperature,
1670 *
const std::vector<double> &old_temperature_laplacians,
1671 *
const std::vector<double> &old_old_temperature_laplacians,
1674 *
const std::vector<double> &gamma_values,
1675 *
const double global_u_infty,
1676 *
const double global_T_variation,
1677 *
const double cell_diameter)
const;
1681 *
double global_Omega_diameter;
1683 *
const unsigned int stokes_degree;
1688 *
std::vector<IndexSet> stokes_partitioning;
1697 *
const unsigned int temperature_degree;
1713 *
double old_time_step;
1714 *
unsigned int timestep_number;
1716 *
std::shared_ptr<TrilinosWrappers::PreconditionAMG> Amg_preconditioner;
1720 * Next, we select the type used
for the preconditioner. For older versions
1721 * of Trilinos, which have the historical Epetra sub-package, we can use an
1722 * incomplete Cholesky (IC) preconditioner. For newer Trilinos versions that
1723 * no longer have Epetra, we use a Jacobi preconditioner.
1726 *
#ifdef DEAL_II_TRILINOS_WITH_EPETRA
1731 *
std::shared_ptr<MpPreconditionType> Mp_preconditioner;
1733 *
bool rebuild_stokes_matrix;
1734 *
bool rebuild_temperature_matrices;
1735 *
bool rebuild_stokes_preconditioner;
1742 * <a name=
"step_31-BoussinesqFlowProblemclassimplementation"></a>
1743 * <h3>BoussinesqFlowProblem
class implementation</h3>
1748 * <a name=
"step_31-BoussinesqFlowProblemBoussinesqFlowProblem"></a>
1749 * <h4>BoussinesqFlowProblem::BoussinesqFlowProblem</h4>
1753 * The constructor of
this class is an extension of the constructor in
1754 * @ref step_22
"step-22". We need to add the various variables that concern the
1755 *
temperature. As discussed in the introduction, we are going to use
1756 * @f$Q_2^d\times Q_1@f$ (Taylor-Hood) elements again
for the Stokes part, and
1757 * @f$Q_2@f$ elements
for the
temperature. However, by
using variables that
1758 * store the polynomial degree of the Stokes and
temperature finite
1759 * elements, it is easy to consistently modify the degree of the elements as
1760 * well as all quadrature formulas used on them
downstream. Moreover, we
1761 * initialize the time stepping as well as the options
for matrix assembly
1762 * and preconditioning:
1765 *
template <
int dim>
1766 *
BoussinesqFlowProblem<dim>::BoussinesqFlowProblem()
1768 *
, global_Omega_diameter(std::numeric_limits<double>::quiet_NaN())
1769 *
, stokes_degree(1)
1771 *
, stokes_dof_handler(triangulation)
1774 *
temperature_degree(2)
1775 *
, temperature_fe(temperature_degree)
1776 *
, temperature_dof_handler(triangulation)
1780 *
, old_time_step(0)
1781 *
, timestep_number(0)
1782 *
, rebuild_stokes_matrix(
true)
1783 *
, rebuild_temperature_matrices(
true)
1784 *
, rebuild_stokes_preconditioner(
true)
1792 * <a name=
"step_31-BoussinesqFlowProblemget_maximal_velocity"></a>
1793 * <h4>BoussinesqFlowProblem::get_maximal_velocity</h4>
1797 * Starting the real functionality of
this class is a helper function that
1798 * determines the maximum (@f$L_\infty@f$) velocity in the domain (at the
1799 * quadrature points, in fact). How it works should be relatively obvious to
1800 * all who have gotten to
this point of the tutorial. Note that since we are
1801 * only interested in the velocity, rather than
using
1802 * <code>stokes_fe_values.get_function_values</code> to get the
values of
1803 * the entire Stokes solution (velocities and pressures) we use
1804 * <code>stokes_fe_values[velocities].get_function_values</code> to
extract
1805 * only the velocities part. This has the additional benefit that we get it
1807 * allowing us to process it right away
using the <code>
norm()</code>
1808 * function to get the magnitude of the velocity.
1812 * The only
point worth thinking about a bit is how to choose the quadrature
1813 * points we use here. Since the goal of
this function is to find the
1814 * maximal velocity over a domain by looking at quadrature points on each
1815 * cell. So we should ask how we should best choose these quadrature points
1816 * on each cell. To
this end, recall that
if we had a single @f$Q_1@f$ field
1817 * (rather than the vector-valued field of higher order) then the maximum
1818 * would be attained at a vertex of the mesh. In other words, we should use
1819 * the
QTrapezoid class that has quadrature points only at the vertices of
1824 * For higher order shape functions, the situation is more complicated: the
1825 * maxima and minima may be attained at points between the support points of
1826 * shape
functions (
for the usual @f$Q_p@f$ elements the support points are the
1827 * equidistant Lagrange interpolation points); furthermore, since we are
1828 * looking
for the maximum magnitude of a vector-valued quantity, we can
1829 * even less say with certainty where the
set of potential maximal points
1830 * are. Nevertheless, intuitively
if not provably, the Lagrange
1831 * interpolation points appear to be a better choice than the Gauss points.
1835 * There are now different methods to produce a quadrature formula with
1836 * quadrature points
equal to the interpolation points of the finite
1837 * element. One option would be to use the
1839 * unique
set of points to avoid duplicate function evaluations, and create
1840 * a
Quadrature object using these points. Another option, chosen here, is
1842 * repeats the
QTrapezoid formula on a number of sub-cells in each coordinate
1843 * direction. To cover all support points, we need to iterate it
1844 * <code>stokes_degree+1</code> times since this is the polynomial degree of
1845 * the Stokes element in use:
1848 *
template <
int dim>
1849 *
double BoussinesqFlowProblem<dim>::get_maximal_velocity() const
1852 *
const unsigned int n_q_points = quadrature_formula.size();
1855 *
std::vector<Tensor<1, dim>> velocity_values(n_q_points);
1856 *
double max_velocity = 0;
1860 *
for (
const auto &cell : stokes_dof_handler.active_cell_iterators())
1862 *
fe_values.
reinit(cell);
1863 *
fe_values[velocities].get_function_values(stokes_solution,
1866 *
for (
unsigned int q = 0; q < n_q_points; ++q)
1867 *
max_velocity =
std::max(max_velocity, velocity_values[q].
norm());
1870 *
return max_velocity;
1878 * <a name=
"step_31-BoussinesqFlowProblemget_extrapolated_temperature_range"></a>
1879 * <h4>BoussinesqFlowProblem::get_extrapolated_temperature_range</h4>
1883 * Next a function that determines the minimum and maximum
temperature at
1884 * quadrature points
inside @f$\Omega@f$ when extrapolated from the two previous
1885 * time steps to the current
one. We need
this information in the
1886 * computation of the artificial viscosity parameter @f$\nu@f$ as discussed in
1892 * @f$\left(1+\frac{k_n}{k_{n-1}} \right)T^{n-1} + \frac{k_n}{k_{n-1}}
1893 *
T^{n-2}@f$. The way to compute it is to
loop over all quadrature points and
1894 * update the maximum and minimum
value if the current
value is
1895 * bigger/smaller than the previous
one. We initialize the variables that
1896 * store the
max and
min before the
loop over all quadrature points by the
1897 * smallest and the largest number representable as a
double. Then we know
1898 *
for a fact that it is larger/smaller than the minimum/maximum and that
1899 * the
loop over all quadrature points is ultimately going to update the
1904 * The only other complication worth mentioning here is that in the
first
1905 * time step, @f$T^{k-2}@f$ is not yet available of course. In that
case, we can
1907 * quadrature points, we use the same choice as in the previous function
1908 * though with the difference that now the number of repetitions is
1909 * determined by the polynomial degree of the
temperature field.
1912 *
template <
int dim>
1913 *
std::pair<double, double>
1914 *
BoussinesqFlowProblem<dim>::get_extrapolated_temperature_range() const
1917 *
temperature_degree);
1918 *
const unsigned int n_q_points = quadrature_formula.size();
1921 *
std::vector<double> old_temperature_values(n_q_points);
1922 *
std::vector<double> old_old_temperature_values(n_q_points);
1924 *
if (timestep_number != 0)
1926 *
double min_temperature = std::numeric_limits<double>::max(),
1927 *
max_temperature = std::numeric_limits<double>::lowest();
1929 *
for (
const auto &cell : temperature_dof_handler.active_cell_iterators())
1931 *
fe_values.
reinit(cell);
1932 *
fe_values.get_function_values(old_temperature_solution,
1933 *
old_temperature_values);
1934 *
fe_values.get_function_values(old_old_temperature_solution,
1935 *
old_old_temperature_values);
1937 *
for (
unsigned int q = 0; q < n_q_points; ++q)
1940 *
(1. + time_step / old_time_step) * old_temperature_values[q] -
1941 *
time_step / old_time_step * old_old_temperature_values[q];
1943 *
min_temperature =
std::min(min_temperature, temperature);
1944 *
max_temperature =
std::max(max_temperature, temperature);
1948 *
return std::make_pair(min_temperature, max_temperature);
1952 *
double min_temperature = std::numeric_limits<double>::max(),
1953 *
max_temperature = std::numeric_limits<double>::lowest();
1955 *
for (
const auto &cell : temperature_dof_handler.active_cell_iterators())
1957 *
fe_values.
reinit(cell);
1958 *
fe_values.get_function_values(old_temperature_solution,
1959 *
old_temperature_values);
1961 *
for (
unsigned int q = 0; q < n_q_points; ++q)
1963 *
const double temperature = old_temperature_values[q];
1965 *
min_temperature =
std::min(min_temperature, temperature);
1966 *
max_temperature =
std::max(max_temperature, temperature);
1970 *
return std::make_pair(min_temperature, max_temperature);
1979 * <a name=
"step_31-BoussinesqFlowProblemcompute_viscosity"></a>
1980 * <h4>BoussinesqFlowProblem::compute_viscosity</h4>
1984 * The last of the tool
functions computes the artificial viscosity
1985 * parameter @f$\nu|_K@f$ on a cell @f$K@f$ as a function of the extrapolated
1987 * the right hand side @f$\gamma@f$ all on the quadrature points of the current
1988 * cell, and various other parameters as described in detail in the
1993 * There are some universal constants worth mentioning here. First, we need
1994 * to fix @f$\beta@f$; we choose @f$\beta=0.017\cdot dim@f$, a choice discussed in
1995 * detail in the results section of
this tutorial program. The
second is the
1996 * exponent @f$\alpha@f$; @f$\alpha=1@f$ appears to work fine
for the current
1997 * program, even though some additional benefit might be expected from
1998 * choosing @f$\alpha = 2@f$. Finally, there is
one thing that
requires special
1999 * casing: In the
first time step, the velocity equals
zero, and the formula
2000 *
for @f$\nu|_K@f$ is not defined. In that
case, we
return @f$\nu|_K=5\cdot 10^3
2001 * \cdot h_K@f$, a choice admittedly more motivated by heuristics than
2002 * anything
else (it is in the same order of magnitude, however, as the
2003 * value returned
for most cells on the
second time step).
2007 * The rest of the function should be mostly obvious based on the material
2008 * discussed in the introduction:
2011 *
template <
int dim>
2012 *
double BoussinesqFlowProblem<dim>::compute_viscosity(
2013 *
const std::vector<double> &old_temperature,
2014 *
const std::vector<double> &old_old_temperature,
2017 *
const std::vector<double> &old_temperature_laplacians,
2018 *
const std::vector<double> &old_old_temperature_laplacians,
2021 *
const std::vector<double> &gamma_values,
2022 *
const double global_u_infty,
2023 *
const double global_T_variation,
2024 *
const double cell_diameter)
const
2026 *
constexpr double beta = 0.017 * dim;
2027 *
constexpr double alpha = 1.0;
2029 *
if (global_u_infty == 0)
2030 *
return 5
e-3 * cell_diameter;
2032 *
const unsigned int n_q_points = old_temperature.size();
2034 *
double max_residual = 0;
2035 *
double max_velocity = 0;
2037 *
for (
unsigned int q = 0; q < n_q_points; ++q)
2040 *
(old_velocity_values[q] + old_old_velocity_values[q]) / 2;
2042 *
const double dT_dt =
2043 *
(old_temperature[q] - old_old_temperature[q]) / old_time_step;
2044 *
const double u_grad_T =
2045 *
u * (old_temperature_grads[q] + old_old_temperature_grads[q]) / 2;
2047 *
const double kappa_Delta_T =
2048 *
EquationData::kappa *
2049 *
(old_temperature_laplacians[q] + old_old_temperature_laplacians[q]) /
2052 *
const double residual =
2053 *
std::abs((dT_dt + u_grad_T - kappa_Delta_T - gamma_values[q]) *
2054 *
std::pow((old_temperature[q] + old_old_temperature[q]) / 2,
2057 *
max_residual =
std::max(residual, max_residual);
2061 *
const double c_R =
std::pow(2., (4. - 2 * alpha) / dim);
2062 *
const double global_scaling = c_R * global_u_infty * global_T_variation *
2063 *
std::pow(global_Omega_diameter, alpha - 2.);
2066 *
beta * max_velocity *
2068 *
std::pow(cell_diameter, alpha) * max_residual / global_scaling));
2076 * <a name=
"step_31-BoussinesqFlowProblemsetup_dofs"></a>
2077 * <h4>BoussinesqFlowProblem::setup_dofs</h4>
2081 * This is the function that sets up the
DoFHandler objects we have here
2083 * to the right sizes the various objects required
for the linear algebra in
2084 *
this program. Its basic operations are similar to what we
do in @ref step_22
"step-22".
2088 * The body of the function
first enumerates all degrees of freedom
for the
2089 * Stokes and temperature systems. For the Stokes part, degrees of freedom
2090 * are then sorted to ensure that velocities precede pressure DoFs so that
2091 * we can partition the Stokes matrix into a @f$2\times 2@f$
matrix. As a
2092 * difference to @ref step_22
"step-22", we
do not perform any additional DoF
2093 * renumbering. In that program, it paid off since our solver was heavily
2094 * dependent on ILU
's, whereas we use AMG here which is not sensitive to the
2095 * DoF numbering. The IC preconditioner for the inversion of the pressure
2096 * mass matrix would of course take advantage of a Cuthill-McKee like
2097 * renumbering, but its costs are low compared to the velocity portion, so
2098 * the additional work does not pay off.
2102 * We then proceed with the generation of the hanging node constraints that
2103 * arise from adaptive grid refinement for both DoFHandler objects. For the
2104 * velocity, we impose no-flux boundary conditions @f$\mathbf{u}\cdot
2105 * \mathbf{n}=0@f$ by adding constraints to the object that already stores the
2106 * hanging node constraints matrix. The second parameter in the function
2107 * describes the first of the velocity components in the total dof vector,
2108 * which is zero here. The variable <code>no_normal_flux_boundaries</code>
2109 * denotes the boundary indicators for which to set the no flux boundary
2110 * conditions; here, this is boundary indicator zero.
2114 * After having done so, we count the number of degrees of freedom in the
2118 * template <int dim>
2119 * void BoussinesqFlowProblem<dim>::setup_dofs()
2121 * std::vector<unsigned int> stokes_sub_blocks(dim + 1, 0);
2122 * stokes_sub_blocks[dim] = 1;
2125 * stokes_dof_handler.distribute_dofs(stokes_fe);
2126 * DoFRenumbering::component_wise(stokes_dof_handler, stokes_sub_blocks);
2128 * stokes_constraints.clear();
2129 * DoFTools::make_hanging_node_constraints(stokes_dof_handler,
2130 * stokes_constraints);
2131 * const std::set<types::boundary_id> no_normal_flux_boundaries = {0};
2132 * VectorTools::compute_no_normal_flux_constraints(stokes_dof_handler,
2134 * no_normal_flux_boundaries,
2135 * stokes_constraints);
2136 * stokes_constraints.close();
2139 * temperature_dof_handler.distribute_dofs(temperature_fe);
2141 * temperature_constraints.clear();
2142 * DoFTools::make_hanging_node_constraints(temperature_dof_handler,
2143 * temperature_constraints);
2144 * temperature_constraints.close();
2147 * const std::vector<types::global_dof_index> stokes_dofs_per_block =
2148 * DoFTools::count_dofs_per_fe_block(stokes_dof_handler, stokes_sub_blocks);
2150 * const types::global_dof_index n_u = stokes_dofs_per_block[0],
2151 * n_p = stokes_dofs_per_block[1],
2152 * n_T = temperature_dof_handler.n_dofs();
2154 * std::cout << "Number of active cells: " << triangulation.n_active_cells()
2155 * << " (on " << triangulation.n_levels() << " levels)" << std::endl
2156 * << "Number of degrees of freedom: " << n_u + n_p + n_T << " ("
2157 * << n_u << '+
' << n_p << '+
' << n_T << ')
' << std::endl
2162 * The next step is to create the sparsity pattern for the Stokes and
2163 * temperature system matrices as well as the preconditioner matrix from
2164 * which we build the Stokes preconditioner. As in @ref step_22 "step-22", we choose to
2165 * create the pattern by
2166 * using the blocked version of DynamicSparsityPattern.
2170 * So, we first release the memory stored in the matrices, then set up an
2171 * object of type BlockDynamicSparsityPattern consisting of
2172 * @f$2\times 2@f$ blocks (for the Stokes system matrix and preconditioner) or
2173 * DynamicSparsityPattern (for the temperature part). We then
2174 * fill these objects with the nonzero pattern, taking into account that
2175 * for the Stokes system matrix, there are no entries in the
2176 * pressure-pressure block (but all velocity vector components couple with
2177 * each other and with the pressure). Similarly, in the Stokes
2178 * preconditioner matrix, only the diagonal blocks are nonzero, since we
2179 * use the vector Laplacian as discussed in the introduction. This
2180 * operator only couples each vector component of the Laplacian with
2181 * itself, but not with the other vector components. (Application of the
2182 * constraints resulting from the no-flux boundary conditions will couple
2183 * vector components at the boundary again, however.)
2187 * When generating the sparsity pattern, we directly apply the constraints
2188 * from hanging nodes and no-flux boundary conditions. This approach was
2189 * already used in @ref step_27 "step-27", but is different from the one in early
2190 * tutorial programs where we first built the original sparsity pattern
2191 * and only then added the entries resulting from constraints. The reason
2192 * for doing so is that later during assembly we are going to distribute
2193 * the constraints immediately when transferring local to global
2194 * dofs. Consequently, there will be no data written at positions of
2195 * constrained degrees of freedom, so we can let the
2196 * DoFTools::make_sparsity_pattern function omit these entries by setting
2197 * the last Boolean flag to <code>false</code>. Once the sparsity pattern
2198 * is ready, we can use it to initialize the Trilinos matrices. Since the
2199 * Trilinos matrices store the sparsity pattern internally, there is no
2200 * need to keep the sparsity pattern around after the initialization of
2204 * stokes_partitioning.resize(2);
2205 * stokes_partitioning[0] = complete_index_set(n_u);
2206 * stokes_partitioning[1] = complete_index_set(n_p);
2208 * stokes_matrix.clear();
2210 * BlockDynamicSparsityPattern dsp(stokes_dofs_per_block,
2211 * stokes_dofs_per_block);
2213 * Table<2, DoFTools::Coupling> coupling(dim + 1, dim + 1);
2215 * for (unsigned int c = 0; c < dim + 1; ++c)
2216 * for (unsigned int d = 0; d < dim + 1; ++d)
2217 * if (!((c == dim) && (d == dim)))
2218 * coupling[c][d] = DoFTools::always;
2220 * coupling[c][d] = DoFTools::none;
2222 * DoFTools::make_sparsity_pattern(
2223 * stokes_dof_handler, coupling, dsp, stokes_constraints, false);
2225 * stokes_matrix.reinit(dsp);
2229 * Amg_preconditioner.reset();
2230 * Mp_preconditioner.reset();
2231 * stokes_preconditioner_matrix.clear();
2233 * BlockDynamicSparsityPattern dsp(stokes_dofs_per_block,
2234 * stokes_dofs_per_block);
2236 * Table<2, DoFTools::Coupling> coupling(dim + 1, dim + 1);
2237 * for (unsigned int c = 0; c < dim + 1; ++c)
2238 * for (unsigned int d = 0; d < dim + 1; ++d)
2240 * coupling[c][d] = DoFTools::always;
2242 * coupling[c][d] = DoFTools::none;
2244 * DoFTools::make_sparsity_pattern(
2245 * stokes_dof_handler, coupling, dsp, stokes_constraints, false);
2247 * stokes_preconditioner_matrix.reinit(dsp);
2252 * The creation of the temperature matrix (or, rather, matrices, since we
2253 * provide a temperature mass matrix and a temperature @ref GlossStiffnessMatrix "stiffness matrix",
2254 * that will be added together for time discretization) follows the
2255 * generation of the Stokes matrix – except that it is much easier
2256 * here since we do not need to take care of any blocks or coupling
2257 * between components. Note how we initialize the three temperature
2258 * matrices: We only use the sparsity pattern for reinitialization of the
2259 * first matrix, whereas we use the previously generated matrix for the
2260 * two remaining reinits. The reason for doing so is that reinitialization
2261 * from an already generated matrix allows Trilinos to reuse the sparsity
2262 * pattern instead of generating a new one for each copy. This saves both
2263 * some time and memory.
2267 * temperature_mass_matrix.clear();
2268 * temperature_stiffness_matrix.clear();
2269 * temperature_matrix.clear();
2271 * DynamicSparsityPattern dsp(n_T, n_T);
2272 * DoFTools::make_sparsity_pattern(temperature_dof_handler,
2274 * temperature_constraints,
2277 * temperature_matrix.reinit(dsp);
2278 * temperature_mass_matrix.reinit(temperature_matrix);
2279 * temperature_stiffness_matrix.reinit(temperature_matrix);
2284 * Lastly, we set the vectors for the Stokes solutions @f$\mathbf u^{n-1}@f$
2285 * and @f$\mathbf u^{n-2}@f$, as well as for the temperatures @f$T^{n}@f$,
2286 * @f$T^{n-1}@f$ and @f$T^{n-2}@f$ (required for time stepping) and all the system
2287 * right hand sides to their correct sizes and block structure:
2290 * IndexSet temperature_partitioning = complete_index_set(n_T);
2291 * stokes_solution.reinit(stokes_partitioning, MPI_COMM_WORLD);
2292 * old_stokes_solution.reinit(stokes_partitioning, MPI_COMM_WORLD);
2293 * stokes_rhs.reinit(stokes_partitioning, MPI_COMM_WORLD);
2295 * temperature_solution.reinit(temperature_partitioning, MPI_COMM_WORLD);
2296 * old_temperature_solution.reinit(temperature_partitioning, MPI_COMM_WORLD);
2297 * old_old_temperature_solution.reinit(temperature_partitioning,
2300 * temperature_rhs.reinit(temperature_partitioning, MPI_COMM_WORLD);
2308 * <a name="step_31-BoussinesqFlowProblemassemble_stokes_preconditioner"></a>
2309 * <h4>BoussinesqFlowProblem::assemble_stokes_preconditioner</h4>
2313 * This function assembles the matrix we use for preconditioning the Stokes
2314 * system. What we need are a vector Laplace matrix on the velocity
2315 * components and a mass matrix weighted by @f$\eta^{-1}@f$ on the pressure
2316 * component. We start by generating a quadrature object of appropriate
2317 * order, the FEValues object that can give values and gradients at the
2318 * quadrature points (together with quadrature weights). Next we create data
2319 * structures for the cell matrix and the relation between local and global
2320 * DoFs. The vectors <code>grad_phi_u</code> and <code>phi_p</code> are
2321 * going to hold the values of the basis functions in order to faster build
2322 * up the local matrices, as was already done in @ref step_22 "step-22". Before we start
2323 * the loop over all active cells, we have to specify which components are
2324 * pressure and which are velocity.
2327 * template <int dim>
2328 * void BoussinesqFlowProblem<dim>::assemble_stokes_preconditioner()
2330 * stokes_preconditioner_matrix = 0;
2332 * const QGauss<dim> quadrature_formula(stokes_degree + 2);
2333 * FEValues<dim> stokes_fe_values(stokes_fe,
2334 * quadrature_formula,
2335 * update_JxW_values | update_values |
2336 * update_gradients);
2338 * const unsigned int dofs_per_cell = stokes_fe.n_dofs_per_cell();
2339 * const unsigned int n_q_points = quadrature_formula.size();
2341 * FullMatrix<double> local_matrix(dofs_per_cell, dofs_per_cell);
2342 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
2344 * std::vector<Tensor<2, dim>> grad_phi_u(dofs_per_cell);
2345 * std::vector<double> phi_p(dofs_per_cell);
2347 * const FEValuesExtractors::Vector velocities(0);
2348 * const FEValuesExtractors::Scalar pressure(dim);
2350 * for (const auto &cell : stokes_dof_handler.active_cell_iterators())
2352 * stokes_fe_values.reinit(cell);
2357 * The creation of the local matrix is rather simple. There are only a
2358 * Laplace term (on the velocity) and a mass matrix weighted by
2359 * @f$\eta^{-1}@f$ to be generated, so the creation of the local matrix is
2360 * done in two lines. Once the local matrix is ready (loop over rows
2361 * and columns in the local matrix on each quadrature point), we get
2362 * the local DoF indices and write the local information into the
2363 * global matrix. We do this as in @ref step_27 "step-27", i.e., we directly apply the
2364 * constraints from hanging nodes locally. By doing so, we don't have
2365 * to
do that afterwards, and we don
't also write into entries of the
2366 * matrix that will actually be set to zero again later when
2367 * eliminating constraints.
2370 * for (unsigned int q = 0; q < n_q_points; ++q)
2372 * for (unsigned int k = 0; k < dofs_per_cell; ++k)
2374 * grad_phi_u[k] = stokes_fe_values[velocities].gradient(k, q);
2375 * phi_p[k] = stokes_fe_values[pressure].value(k, q);
2378 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
2379 * for (unsigned int j = 0; j < dofs_per_cell; ++j)
2380 * local_matrix(i, j) +=
2381 * (EquationData::eta *
2382 * scalar_product(grad_phi_u[i], grad_phi_u[j]) +
2383 * (1. / EquationData::eta) * phi_p[i] * phi_p[j]) *
2384 * stokes_fe_values.JxW(q);
2387 * cell->get_dof_indices(local_dof_indices);
2388 * stokes_constraints.distribute_local_to_global(
2389 * local_matrix, local_dof_indices, stokes_preconditioner_matrix);
2392 * stokes_preconditioner_matrix.compress(VectorOperation::add);
2400 * <a name="step_31-BoussinesqFlowProblembuild_stokes_preconditioner"></a>
2401 * <h4>BoussinesqFlowProblem::build_stokes_preconditioner</h4>
2405 * This function generates the inner preconditioners that are going to be
2406 * used for the Schur complement block preconditioner. Since the
2407 * preconditioners need only to be regenerated when the matrices change,
2408 * this function does not have to do anything in case the matrices have not
2409 * changed (i.e., the flag <code>rebuild_stokes_preconditioner</code> has
2410 * the value <code>false</code>). Otherwise its first task is to call
2411 * <code>assemble_stokes_preconditioner</code> to generate the
2412 * preconditioner matrices.
2416 * Next, we set up the preconditioner for the velocity-velocity matrix
2417 * @f$A@f$. As explained in the introduction, we are going to use an AMG
2418 * preconditioner based on a vector Laplace matrix @f$\hat{A}@f$ (which is
2419 * spectrally close to the Stokes matrix @f$A@f$). Usually, the
2420 * TrilinosWrappers::PreconditionAMG class can be seen as a good black-box
2421 * preconditioner which does not need any special knowledge. In this case,
2422 * however, we have to be careful: since we build an AMG for a vector
2423 * problem, we have to tell the preconditioner setup which dofs belong to
2424 * which vector component. We do this using the function
2425 * DoFTools::extract_constant_modes, a function that generates a set of
2426 * <code>dim</code> vectors, where each one has ones in the respective
2427 * component of the vector problem and zeros elsewhere. Hence, these are the
2428 * constant modes on each component, which explains the name of the
2432 * template <int dim>
2433 * void BoussinesqFlowProblem<dim>::build_stokes_preconditioner()
2435 * if (rebuild_stokes_preconditioner == false)
2438 * std::cout << " Rebuilding Stokes preconditioner..." << std::flush;
2440 * assemble_stokes_preconditioner();
2442 * Amg_preconditioner = std::make_shared<TrilinosWrappers::PreconditionAMG>();
2444 * const FEValuesExtractors::Vector velocity_components(0);
2445 * const std::vector<std::vector<bool>> constant_modes =
2446 * DoFTools::extract_constant_modes(
2447 * stokes_dof_handler, stokes_fe.component_mask(velocity_components));
2448 * TrilinosWrappers::PreconditionAMG::AdditionalData amg_data;
2449 * amg_data.constant_modes = constant_modes;
2453 * Next, we set some more options of the AMG preconditioner. In
2454 * particular, we need to tell the AMG setup that we use quadratic basis
2455 * functions for the velocity matrix (this implies more nonzero elements
2456 * in the matrix, so that a more robust algorithm needs to be chosen
2457 * internally). Moreover, we want to be able to control how the coarsening
2458 * structure is build up. The way the Trilinos smoothed aggregation AMG
2459 * does this is to look which matrix entries are of similar size as the
2460 * diagonal entry in order to algebraically build a coarse-grid
2461 * structure. By setting the parameter <code>aggregation_threshold</code>
2462 * to 0.02, we specify that all entries that are more than two percent of
2463 * size of some diagonal pivots in that row should form one coarse grid
2464 * point. This parameter is rather ad hoc, and some fine-tuning of it can
2465 * influence the performance of the preconditioner. As a rule of thumb,
2466 * larger values of <code>aggregation_threshold</code> will decrease the
2467 * number of iterations, but increase the costs per iteration. A look at
2468 * the Trilinos documentation will provide more information on these
2469 * parameters. With this data set, we then initialize the preconditioner
2470 * with the matrix we want it to apply to.
2474 * Finally, we also initialize the preconditioner for the inversion of the
2475 * pressure mass matrix. This matrix is symmetric and well-behaved, so we
2476 * can chose a simple preconditioner. We stick with an incomplete Cholesky
2477 * (IC) factorization preconditioner, which is designed for symmetric
2478 * matrices. We could have also chosen an SSOR preconditioner with
2479 * relaxation factor around 1.2, but IC is cheaper for our example. We
2480 * wrap the preconditioners into a <code>std::shared_ptr</code>
2481 * pointer, which makes it easier to recreate the preconditioner next time
2482 * around since we do not have to care about destroying the previously
2486 * amg_data.elliptic = true;
2487 * #ifdef DEAL_II_TRILINOS_WITH_EPETRA
2488 * amg_data.higher_order_elements = true;
2490 * amg_data.smoother_sweeps = 2;
2491 * amg_data.aggregation_threshold = 0.02;
2492 * Amg_preconditioner->initialize(stokes_preconditioner_matrix.block(0, 0),
2495 * Mp_preconditioner = std::make_shared<MpPreconditionType>();
2496 * Mp_preconditioner->initialize(stokes_preconditioner_matrix.block(1, 1));
2498 * std::cout << std::endl;
2500 * rebuild_stokes_preconditioner = false;
2508 * <a name="step_31-BoussinesqFlowProblemassemble_stokes_system"></a>
2509 * <h4>BoussinesqFlowProblem::assemble_stokes_system</h4>
2513 * The time lag scheme we use for advancing the coupled Stokes-temperature
2514 * system forces us to split up the assembly (and the solution of linear
2515 * systems) into two step. The first one is to create the Stokes system
2516 * matrix and right hand side, and the second is to create matrix and right
2517 * hand sides for the temperature dofs, which depends on the result of the
2518 * linear system for the velocity.
2522 * This function is called at the beginning of each time step. In the first
2523 * time step or if the mesh has changed, indicated by the
2524 * <code>rebuild_stokes_matrix</code>, we need to assemble the Stokes
2525 * matrix; on the other hand, if the mesh hasn't changed and the
matrix is
2526 * already available,
this is not necessary and all we need to
do is
2527 *
assemble the right hand side vector which changes in each time step.
2531 * Regarding the technical details of implementation, not much has changed
2532 * from @ref step_22
"step-22". We reset
matrix and vector, create a quadrature formula on
2533 * the cells, and then create the respective
FEValues object. For the update
2534 * flags, we require basis function derivatives only in
case of a full
2535 * assembly, since they are not needed
for the right hand side; as
always,
2536 * choosing the minimal
set of flags depending on what is currently needed
2542 * There is
one thing that needs to be commented – since we have a
2544 * generate a
second FEValues object for the proper evaluation of the
2545 *
temperature solution. This isn
't too complicated to realize here: just
2546 * use the temperature structures and set an update flag for the basis
2547 * function values which we need for evaluation of the temperature
2548 * solution. The only important part to remember here is that the same
2549 * quadrature formula is used for both FEValues objects to ensure that we
2550 * get matching information when we loop over the quadrature points of the
2555 * The declarations proceed with some shortcuts for array sizes, the
2556 * creation of the local matrix and right hand side as well as the vector
2557 * for the indices of the local dofs compared to the global system.
2560 * template <int dim>
2561 * void BoussinesqFlowProblem<dim>::assemble_stokes_system()
2563 * std::cout << " Assembling..." << std::flush;
2565 * if (rebuild_stokes_matrix == true)
2566 * stokes_matrix = 0;
2570 * const QGauss<dim> quadrature_formula(stokes_degree + 2);
2571 * FEValues<dim> stokes_fe_values(
2573 * quadrature_formula,
2574 * update_values | update_quadrature_points | update_JxW_values |
2575 * (rebuild_stokes_matrix == true ? update_gradients : UpdateFlags(0)));
2577 * FEValues<dim> temperature_fe_values(temperature_fe,
2578 * quadrature_formula,
2581 * const unsigned int dofs_per_cell = stokes_fe.n_dofs_per_cell();
2582 * const unsigned int n_q_points = quadrature_formula.size();
2584 * FullMatrix<double> local_matrix(dofs_per_cell, dofs_per_cell);
2585 * Vector<double> local_rhs(dofs_per_cell);
2587 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
2591 * Next we need a vector that will contain the values of the temperature
2592 * solution at the previous time level at the quadrature points to
2593 * assemble the source term in the right hand side of the momentum
2594 * equation. Let's call
this vector <code>old_solution_values</code>.
2598 * The
set of vectors we create next hold the evaluations of the basis
2600 * will be used
for creating the matrices. Putting these into their own
2601 * arrays rather than asking the
FEValues object for this information each
2602 * time it is needed is an optimization to accelerate the assembly
2603 * process, see @ref step_22
"step-22" for details.
2607 * The last two declarations are used to
extract the individual blocks
2608 * (velocity, pressure,
temperature) from the total FE system.
2611 *
std::vector<double> old_temperature_values(n_q_points);
2613 *
std::vector<Tensor<1, dim>> phi_u(dofs_per_cell);
2614 *
std::vector<SymmetricTensor<2, dim>> grads_phi_u(dofs_per_cell);
2615 *
std::vector<double> div_phi_u(dofs_per_cell);
2616 *
std::vector<double> phi_p(dofs_per_cell);
2623 * Now start the
loop over all cells in the problem. We are working on two
2624 * different DoFHandlers
for this assembly routine, so we must have two
2625 * different cell iterators
for the two objects in use. This might seem a
2626 * bit peculiar, since both the Stokes system and the
temperature system
2627 * use the same grid, but that
's the only way to keep degrees of freedom
2628 * in sync. The first statements within the loop are again all very
2629 * familiar, doing the update of the finite element data as specified by
2630 * the update flags, zeroing out the local arrays and getting the values
2631 * of the old solution at the quadrature points. Then we are ready to loop
2632 * over the quadrature points on the cell.
2635 * auto cell = stokes_dof_handler.begin_active();
2636 * const auto endc = stokes_dof_handler.end();
2637 * auto temperature_cell = temperature_dof_handler.begin_active();
2639 * for (; cell != endc; ++cell, ++temperature_cell)
2641 * stokes_fe_values.reinit(cell);
2642 * temperature_fe_values.reinit(temperature_cell);
2647 * temperature_fe_values.get_function_values(old_temperature_solution,
2648 * old_temperature_values);
2650 * for (unsigned int q = 0; q < n_q_points; ++q)
2652 * const double old_temperature = old_temperature_values[q];
2656 * Next we extract the values and gradients of basis functions
2657 * relevant to the terms in the inner products. As shown in
2658 * @ref step_22 "step-22" this helps accelerate assembly.
2662 * Once this is done, we start the loop over the rows and columns
2663 * of the local matrix and feed the matrix with the relevant
2664 * products. The right hand side is filled with the forcing term
2665 * driven by temperature in direction of gravity (which is
2666 * vertical in our example). Note that the right hand side term
2667 * is always generated, whereas the matrix contributions are only
2668 * updated when it is requested by the
2669 * <code>rebuild_matrices</code> flag.
2672 * for (unsigned int k = 0; k < dofs_per_cell; ++k)
2674 * phi_u[k] = stokes_fe_values[velocities].value(k, q);
2675 * if (rebuild_stokes_matrix)
2678 * stokes_fe_values[velocities].symmetric_gradient(k, q);
2680 * stokes_fe_values[velocities].divergence(k, q);
2681 * phi_p[k] = stokes_fe_values[pressure].value(k, q);
2685 * if (rebuild_stokes_matrix)
2686 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
2687 * for (unsigned int j = 0; j < dofs_per_cell; ++j)
2688 * local_matrix(i, j) +=
2689 * (EquationData::eta * 2 * (grads_phi_u[i] * grads_phi_u[j]) -
2690 * div_phi_u[i] * phi_p[j] - phi_p[i] * div_phi_u[j]) *
2691 * stokes_fe_values.JxW(q);
2693 * const Point<dim> gravity =
2694 * -((dim == 2) ? (Point<dim>(0, 1)) : (Point<dim>(0, 0, 1)));
2695 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
2696 * local_rhs(i) += (-EquationData::density * EquationData::beta *
2697 * gravity * phi_u[i] * old_temperature) *
2698 * stokes_fe_values.JxW(q);
2703 * The last step in the loop over all cells is to enter the local
2704 * contributions into the global matrix and vector structures to the
2705 * positions specified in <code>local_dof_indices</code>. Again, we
2706 * let the AffineConstraints class do the insertion of the cell
2707 * matrix elements to the global matrix, which already condenses the
2708 * hanging node constraints.
2711 * cell->get_dof_indices(local_dof_indices);
2713 * if (rebuild_stokes_matrix == true)
2714 * stokes_constraints.distribute_local_to_global(local_matrix,
2716 * local_dof_indices,
2720 * stokes_constraints.distribute_local_to_global(local_rhs,
2721 * local_dof_indices,
2725 * stokes_matrix.compress(VectorOperation::add);
2726 * stokes_rhs.compress(VectorOperation::add);
2728 * rebuild_stokes_matrix = false;
2730 * std::cout << std::endl;
2738 * <a name="step_31-BoussinesqFlowProblemassemble_temperature_matrix"></a>
2739 * <h4>BoussinesqFlowProblem::assemble_temperature_matrix</h4>
2743 * This function assembles the matrix in the temperature equation. The
2744 * temperature matrix consists of two parts, a mass matrix and the time step
2745 * size times a stiffness matrix given by a Laplace term times the amount of
2746 * diffusion. Since the matrix depends on the time step size (which varies
2747 * from one step to another), the temperature matrix needs to be updated
2748 * every time step. We could simply regenerate the matrices in every time
2749 * step, but this is not really efficient since mass and Laplace matrix do
2750 * only change when we change the mesh. Hence, we do this more efficiently
2751 * by generating two separate matrices in this function, one for the mass
2752 * matrix and one for the stiffness (diffusion) matrix. We will then sum up
2753 * the matrix plus the stiffness matrix times the time step size once we
2754 * know the actual time step.
2758 * So the details for this first step are very simple. In case we need to
2759 * rebuild the matrix (i.e., the mesh has changed), we zero the data
2760 * structures, get a quadrature formula and a FEValues object, and create
2761 * local matrices, local dof indices and evaluation structures for the basis
2765 * template <int dim>
2766 * void BoussinesqFlowProblem<dim>::assemble_temperature_matrix()
2768 * if (rebuild_temperature_matrices == false)
2771 * temperature_mass_matrix = 0;
2772 * temperature_stiffness_matrix = 0;
2774 * const QGauss<dim> quadrature_formula(temperature_degree + 2);
2775 * FEValues<dim> temperature_fe_values(temperature_fe,
2776 * quadrature_formula,
2777 * update_values | update_gradients |
2778 * update_JxW_values);
2780 * const unsigned int dofs_per_cell = temperature_fe.n_dofs_per_cell();
2781 * const unsigned int n_q_points = quadrature_formula.size();
2783 * FullMatrix<double> local_mass_matrix(dofs_per_cell, dofs_per_cell);
2784 * FullMatrix<double> local_stiffness_matrix(dofs_per_cell, dofs_per_cell);
2786 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
2788 * std::vector<double> phi_T(dofs_per_cell);
2789 * std::vector<Tensor<1, dim>> grad_phi_T(dofs_per_cell);
2793 * Now, let's start the
loop over all cells in the triangulation. We need
2794 * to
zero out the local matrices, update the finite element evaluations,
2795 * and then
loop over the rows and columns of the matrices on each
2796 * quadrature
point, where we then create the mass
matrix and the
2797 * stiffness
matrix (Laplace terms times the diffusion
2798 * <code>EquationData::kappa</code>. Finally, we let the constraints
2799 *
object insert these values into the global matrix, and directly
2800 * condense the constraints into the
matrix.
2803 *
for (
const auto &cell : temperature_dof_handler.active_cell_iterators())
2805 *
local_mass_matrix = 0;
2806 *
local_stiffness_matrix = 0;
2808 *
temperature_fe_values.reinit(cell);
2810 *
for (
unsigned int q = 0; q < n_q_points; ++q)
2812 *
for (
unsigned int k = 0; k < dofs_per_cell; ++k)
2814 *
grad_phi_T[k] = temperature_fe_values.shape_grad(k, q);
2815 *
phi_T[k] = temperature_fe_values.shape_value(k, q);
2818 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
2819 *
for (
unsigned int j = 0; j < dofs_per_cell; ++j)
2821 *
local_mass_matrix(i, j) +=
2822 *
(phi_T[i] * phi_T[j] * temperature_fe_values.JxW(q));
2823 *
local_stiffness_matrix(i, j) +=
2824 *
(EquationData::kappa * grad_phi_T[i] * grad_phi_T[j] *
2825 *
temperature_fe_values.JxW(q));
2829 *
cell->get_dof_indices(local_dof_indices);
2831 *
temperature_constraints.distribute_local_to_global(
2832 *
local_mass_matrix, local_dof_indices, temperature_mass_matrix);
2833 *
temperature_constraints.distribute_local_to_global(
2834 *
local_stiffness_matrix,
2835 *
local_dof_indices,
2836 *
temperature_stiffness_matrix);
2841 *
rebuild_temperature_matrices =
false;
2849 * <a name=
"step_31-BoussinesqFlowProblemassemble_temperature_system"></a>
2850 * <h4>BoussinesqFlowProblem::assemble_temperature_system</h4>
2854 * This function does the
second part of the assembly work on the
2856 *
matrix (where the time step
size comes into play), as well as the
2857 * creation of the velocity-dependent right hand side. The declarations
for
2858 * the right hand side assembly in
this function are pretty much the same as
2859 * the ones used in the other assembly routines, except that we restrict
2860 * ourselves to vectors
this time. We are going to calculate residuals on
2862 * derivatives, specified by the update flag <code>
update_hessians</code>.
2866 * The
temperature equation is coupled to the Stokes system by means of the
2867 * fluid velocity. These two parts of the solution are associated with
2868 * different DoFHandlers, so we again need to create a
second FEValues
2869 *
object for the evaluation of the velocity at the quadrature points.
2872 *
template <
int dim>
2873 *
void BoussinesqFlowProblem<dim>::assemble_temperature_system(
2874 *
const double maximal_velocity)
2876 *
const bool use_bdf2_scheme = (timestep_number != 0);
2878 *
if (use_bdf2_scheme ==
true)
2880 *
temperature_matrix.copy_from(temperature_mass_matrix);
2881 *
temperature_matrix *=
2882 *
(2 * time_step + old_time_step) / (time_step + old_time_step);
2883 *
temperature_matrix.add(time_step, temperature_stiffness_matrix);
2887 *
temperature_matrix.copy_from(temperature_mass_matrix);
2888 *
temperature_matrix.add(time_step, temperature_stiffness_matrix);
2891 *
temperature_rhs = 0;
2893 *
const QGauss<dim> quadrature_formula(temperature_degree + 2);
2895 *
quadrature_formula,
2901 *
quadrature_formula,
2904 *
const unsigned int dofs_per_cell = temperature_fe.n_dofs_per_cell();
2905 *
const unsigned int n_q_points = quadrature_formula.size();
2909 *
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
2913 * Next comes the declaration of vectors to hold the old and older
2914 * solution
values (as a notation
for time levels @f$n-1@f$ and
2915 * @f$n-2@f$, respectively) and
gradients at quadrature points of the
2916 * current cell. We also declare an
object to hold the
temperature right
2917 * hand side
values (<code>gamma_values</code>), and we again use
2920 * domain which will be used
for the definition of the stabilization
2921 * parameter (we got the maximal velocity as an input to
this function).
2924 *
std::vector<Tensor<1, dim>> old_velocity_values(n_q_points);
2925 *
std::vector<Tensor<1, dim>> old_old_velocity_values(n_q_points);
2926 *
std::vector<double> old_temperature_values(n_q_points);
2927 *
std::vector<double> old_old_temperature_values(n_q_points);
2928 *
std::vector<Tensor<1, dim>> old_temperature_grads(n_q_points);
2929 *
std::vector<Tensor<1, dim>> old_old_temperature_grads(n_q_points);
2930 *
std::vector<double> old_temperature_laplacians(n_q_points);
2931 *
std::vector<double> old_old_temperature_laplacians(n_q_points);
2933 *
EquationData::TemperatureRightHandSide<dim> temperature_right_hand_side;
2934 *
std::vector<double> gamma_values(n_q_points);
2936 *
std::vector<double> phi_T(dofs_per_cell);
2937 *
std::vector<Tensor<1, dim>> grad_phi_T(dofs_per_cell);
2939 *
const std::pair<double, double> global_T_range =
2940 *
get_extrapolated_temperature_range();
2946 * Now, let
's start the loop over all cells in the triangulation. Again,
2947 * we need two cell iterators that walk in parallel through the cells of
2948 * the two involved DoFHandler objects for the Stokes and temperature
2949 * part. Within the loop, we first set the local rhs to zero, and then get
2950 * the values and derivatives of the old solution functions at the
2951 * quadrature points, since they are going to be needed for the definition
2952 * of the stabilization parameters and as coefficients in the equation,
2953 * respectively. Note that since the temperature has its own DoFHandler
2954 * and FEValues object we get the entire solution at the quadrature point
2955 * (which is the scalar temperature field only anyway) whereas for the
2956 * Stokes part we restrict ourselves to extracting the velocity part (and
2957 * ignoring the pressure part) by using
2958 * <code>stokes_fe_values[velocities].get_function_values</code>.
2961 * auto cell = temperature_dof_handler.begin_active();
2962 * const auto endc = temperature_dof_handler.end();
2963 * auto stokes_cell = stokes_dof_handler.begin_active();
2965 * for (; cell != endc; ++cell, ++stokes_cell)
2969 * temperature_fe_values.reinit(cell);
2970 * stokes_fe_values.reinit(stokes_cell);
2972 * temperature_fe_values.get_function_values(old_temperature_solution,
2973 * old_temperature_values);
2974 * temperature_fe_values.get_function_values(old_old_temperature_solution,
2975 * old_old_temperature_values);
2977 * temperature_fe_values.get_function_gradients(old_temperature_solution,
2978 * old_temperature_grads);
2979 * temperature_fe_values.get_function_gradients(
2980 * old_old_temperature_solution, old_old_temperature_grads);
2982 * temperature_fe_values.get_function_laplacians(
2983 * old_temperature_solution, old_temperature_laplacians);
2984 * temperature_fe_values.get_function_laplacians(
2985 * old_old_temperature_solution, old_old_temperature_laplacians);
2987 * temperature_right_hand_side.value_list(
2988 * temperature_fe_values.get_quadrature_points(), gamma_values);
2990 * stokes_fe_values[velocities].get_function_values(stokes_solution,
2991 * old_velocity_values);
2992 * stokes_fe_values[velocities].get_function_values(
2993 * old_stokes_solution, old_old_velocity_values);
2997 * Next, we calculate the artificial viscosity for stabilization
2998 * according to the discussion in the introduction using the dedicated
2999 * function. With that at hand, we can get into the loop over
3000 * quadrature points and local rhs vector components. The terms here
3001 * are quite lengthy, but their definition follows the time-discrete
3002 * system developed in the introduction of this program. The BDF-2
3003 * scheme needs one more term from the old time step (and involves
3004 * more complicated factors) than the backward Euler scheme that is
3005 * used for the first time step. When all this is done, we distribute
3006 * the local vector into the global one (including hanging node
3011 * compute_viscosity(old_temperature_values,
3012 * old_old_temperature_values,
3013 * old_temperature_grads,
3014 * old_old_temperature_grads,
3015 * old_temperature_laplacians,
3016 * old_old_temperature_laplacians,
3017 * old_velocity_values,
3018 * old_old_velocity_values,
3021 * global_T_range.second - global_T_range.first,
3022 * cell->diameter());
3024 * for (unsigned int q = 0; q < n_q_points; ++q)
3026 * for (unsigned int k = 0; k < dofs_per_cell; ++k)
3028 * grad_phi_T[k] = temperature_fe_values.shape_grad(k, q);
3029 * phi_T[k] = temperature_fe_values.shape_value(k, q);
3032 * const double T_term_for_rhs =
3033 * (use_bdf2_scheme ?
3034 * (old_temperature_values[q] * (1 + time_step / old_time_step) -
3035 * old_old_temperature_values[q] * (time_step * time_step) /
3036 * (old_time_step * (time_step + old_time_step))) :
3037 * old_temperature_values[q]);
3039 * const Tensor<1, dim> ext_grad_T =
3040 * (use_bdf2_scheme ?
3041 * (old_temperature_grads[q] * (1 + time_step / old_time_step) -
3042 * old_old_temperature_grads[q] * time_step / old_time_step) :
3043 * old_temperature_grads[q]);
3045 * const Tensor<1, dim> extrapolated_u =
3046 * (use_bdf2_scheme ?
3047 * (old_velocity_values[q] * (1 + time_step / old_time_step) -
3048 * old_old_velocity_values[q] * time_step / old_time_step) :
3049 * old_velocity_values[q]);
3051 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
3053 * (T_term_for_rhs * phi_T[i] -
3054 * time_step * extrapolated_u * ext_grad_T * phi_T[i] -
3055 * time_step * nu * ext_grad_T * grad_phi_T[i] +
3056 * time_step * gamma_values[q] * phi_T[i]) *
3057 * temperature_fe_values.JxW(q);
3060 * cell->get_dof_indices(local_dof_indices);
3061 * temperature_constraints.distribute_local_to_global(local_rhs,
3062 * local_dof_indices,
3072 * <a name="step_31-BoussinesqFlowProblemsolve"></a>
3073 * <h4>BoussinesqFlowProblem::solve</h4>
3077 * This function solves the linear systems of equations. Following the
3078 * introduction, we start with the Stokes system, where we need to generate
3079 * our block Schur preconditioner. Since all the relevant actions are
3080 * implemented in the class <code>BlockSchurPreconditioner</code>, all we
3081 * have to do is to initialize the class appropriately. What we need to pass
3082 * down is an <code>InverseMatrix</code> object for the pressure mass
3083 * matrix, which we set up using the respective class together with the IC
3084 * preconditioner we already generated, and the AMG preconditioner for the
3085 * velocity-velocity matrix. Note that both <code>Mp_preconditioner</code>
3086 * and <code>Amg_preconditioner</code> are only pointers, so we use
3087 * <code>*</code> to pass down the actual preconditioner objects.
3091 * Once the preconditioner is ready, we create a GMRES solver for the block
3092 * system. Since we are working with Trilinos data structures, we have to
3093 * set the respective template argument in the solver. GMRES needs to
3094 * internally store temporary vectors for each iteration (see the discussion
3095 * in the results section of @ref step_22 "step-22") – the more vectors it can use,
3096 * the better it will generally perform. To keep memory demands in check, we
3097 * set the number of vectors to 100. This means that up to 100 solver
3098 * iterations, every temporary vector can be stored. If the solver needs to
3099 * iterate more often to get the specified tolerance, it will work on a
3100 * reduced set of vectors by restarting at every 100 iterations.
3104 * With this all set up, we solve the system and distribute the constraints
3105 * in the Stokes system, i.e., hanging nodes and no-flux boundary condition,
3106 * in order to have the appropriate solution values even at constrained
3107 * dofs. Finally, we write the number of iterations to the screen.
3110 * template <int dim>
3111 * void BoussinesqFlowProblem<dim>::solve()
3113 * std::cout << " Solving..." << std::endl;
3116 * const LinearSolvers::InverseMatrix<TrilinosWrappers::SparseMatrix,
3117 * MpPreconditionType>
3118 * mp_inverse(stokes_preconditioner_matrix.block(1, 1),
3119 * *Mp_preconditioner);
3121 * const LinearSolvers::BlockSchurPreconditioner<
3122 * TrilinosWrappers::PreconditionAMG,
3123 * MpPreconditionType>
3124 * preconditioner(stokes_matrix, mp_inverse, *Amg_preconditioner);
3126 * SolverControl solver_control(stokes_matrix.m(),
3127 * 1e-6 * stokes_rhs.l2_norm());
3129 * SolverGMRES<TrilinosWrappers::MPI::BlockVector> gmres(
3131 * SolverGMRES<TrilinosWrappers::MPI::BlockVector>::AdditionalData(100));
3133 * for (unsigned int i = 0; i < stokes_solution.size(); ++i)
3134 * if (stokes_constraints.is_constrained(i))
3135 * stokes_solution(i) = 0;
3137 * gmres.solve(stokes_matrix, stokes_solution, stokes_rhs, preconditioner);
3139 * stokes_constraints.distribute(stokes_solution);
3141 * std::cout << " " << solver_control.last_step()
3142 * << " GMRES iterations for Stokes subsystem." << std::endl;
3147 * Once we know the Stokes solution, we can determine the new time step
3148 * from the maximal velocity. We have to do this to satisfy the CFL
3149 * condition since convection terms are treated explicitly in the
3150 * temperature equation, as discussed in the introduction. The exact form
3151 * of the formula used here for the time step is discussed in the results
3152 * section of this program.
3156 * There is a snatch here. The formula contains a division by the maximum
3157 * value of the velocity. However, at the start of the computation, we
3158 * have a constant temperature field (we start with a constant
3159 * temperature, and it will be nonconstant only after the first time step
3160 * during which the source acts). Constant temperature means that no
3161 * buoyancy acts, and so the velocity is zero. Dividing by it will not
3162 * likely lead to anything good.
3166 * To avoid the resulting infinite time step, we ask whether the maximal
3167 * velocity is very small (in particular smaller than the values we
3168 * encounter during any of the following time steps) and if so rather than
3169 * dividing by zero we just divide by a small value, resulting in a large
3170 * but finite time step.
3173 * old_time_step = time_step;
3174 * const double maximal_velocity = get_maximal_velocity();
3176 * if (maximal_velocity >= 0.01)
3177 * time_step = 1. / (1.7 * dim * std::sqrt(1. * dim)) / temperature_degree *
3178 * GridTools::minimal_cell_diameter(triangulation) /
3181 * time_step = 1. / (1.7 * dim * std::sqrt(1. * dim)) / temperature_degree *
3182 * GridTools::minimal_cell_diameter(triangulation) / .01;
3185 * << "Time step: " << time_step << std::endl;
3187 * temperature_solution = old_temperature_solution;
3191 * Next we set up the temperature system and the right hand side using the
3192 * function <code>assemble_temperature_system()</code>. Knowing the
3193 * matrix and right hand side of the temperature equation, we set up a
3194 * preconditioner and a solver. The temperature matrix is a mass matrix
3195 * (with eigenvalues around one) plus a Laplace matrix (with eigenvalues
3196 * between zero and @f$ch^{-2}@f$) times a small number proportional to the
3197 * time step @f$k_n@f$. Hence, the resulting symmetric and positive definite
3198 * matrix has eigenvalues in the range @f$[1,1+k_nh^{-2}]@f$ (up to
3199 * constants). This matrix is only moderately ill conditioned even for
3200 * small mesh sizes and we get a reasonably good preconditioner by simple
3201 * means, for example with an incomplete Cholesky decomposition
3202 * preconditioner (IC) as we also use for preconditioning the pressure
3203 * mass matrix solver. As a solver, we choose the conjugate gradient
3204 * method CG. As before, we tell the solver to use Trilinos vectors via
3205 * the template argument <code>TrilinosWrappers::MPI::Vector</code>.
3206 * Finally, we solve, distribute the hanging node constraints and write out
3207 * the number of iterations.
3210 * assemble_temperature_system(maximal_velocity);
3212 * SolverControl solver_control(temperature_matrix.m(),
3213 * 1e-8 * temperature_rhs.l2_norm());
3214 * SolverCG<TrilinosWrappers::MPI::Vector> cg(solver_control);
3216 * MpPreconditionType preconditioner;
3217 * preconditioner.initialize(temperature_matrix);
3219 * cg.solve(temperature_matrix,
3220 * temperature_solution,
3224 * temperature_constraints.distribute(temperature_solution);
3226 * std::cout << " " << solver_control.last_step()
3227 * << " CG iterations for temperature." << std::endl;
3231 * At the end of this function, we step through the vector and read out
3232 * the maximum and minimum temperature value, which we also want to
3233 * output. This will come in handy when determining the correct constant
3234 * in the choice of time step as discuss in the results section of this
3238 * double min_temperature = temperature_solution(0),
3239 * max_temperature = temperature_solution(0);
3240 * for (unsigned int i = 0; i < temperature_solution.size(); ++i)
3243 * std::min<double>(min_temperature, temperature_solution(i));
3245 * std::max<double>(max_temperature, temperature_solution(i));
3248 * std::cout << " Temperature range: " << min_temperature << ' '
3249 * << max_temperature << std::endl;
3258 * <a name="step_31-BoussinesqFlowProblemoutput_results"></a>
3259 * <h4>BoussinesqFlowProblem::output_results</h4>
3263 * This function writes the solution to a VTK output file for visualization,
3264 * which is done every tenth time step. This is usually quite a simple task,
3265 * since the deal.II library provides functions that do almost all the job
3266 * for us. There is one new function compared to previous examples: We want
3267 * to visualize both the Stokes solution and the temperature as one data
3268 * set, but we have done all the calculations based on two different
3269 * DoFHandler objects. Luckily, the DataOut class is prepared to deal with
3270 * it. All we have to do is to not attach one single DoFHandler at the
3271 * beginning and then use that for all added vector, but specify the
3272 * DoFHandler to each vector separately. The rest is done as in @ref step_22 "step-22". We
3273 * create solution names (that are going to appear in the visualization
3274 * program for the individual components). The first <code>dim</code>
3275 * components are the vector velocity, and then we have pressure for the
3276 * Stokes part, whereas temperature is scalar. This information is read out
3277 * using the DataComponentInterpretation helper class. Next, we actually
3278 * attach the data vectors with their DoFHandler objects, build patches
3279 * according to the degree of freedom, which are (sub-) elements that
3280 * describe the data for visualization programs. Finally, we open a file
3281 * (that includes the time step number) and write the vtk data into it.
3284 * template <int dim>
3285 * void BoussinesqFlowProblem<dim>::output_results() const
3287 * if (timestep_number % 10 != 0)
3290 * std::vector<std::string> stokes_names(dim, "velocity");
3291 * stokes_names.emplace_back("p");
3292 * std::vector<DataComponentInterpretation::DataComponentInterpretation>
3293 * stokes_component_interpretation(
3294 * dim + 1, DataComponentInterpretation::component_is_scalar);
3295 * for (unsigned int i = 0; i < dim; ++i)
3296 * stokes_component_interpretation[i] =
3297 * DataComponentInterpretation::component_is_part_of_vector;
3299 * DataOut<dim> data_out;
3300 * data_out.add_data_vector(stokes_dof_handler,
3303 * stokes_component_interpretation);
3304 * data_out.add_data_vector(temperature_dof_handler,
3305 * temperature_solution,
3307 * data_out.build_patches(std::min(stokes_degree, temperature_degree));
3309 * std::ofstream output("solution-" +
3310 * Utilities::int_to_string(timestep_number, 4) + ".vtk");
3311 * data_out.write_vtk(output);
3319 * <a name="step_31-BoussinesqFlowProblemrefine_mesh"></a>
3320 * <h4>BoussinesqFlowProblem::refine_mesh</h4>
3324 * This function takes care of the adaptive mesh refinement. The three tasks
3325 * this function performs is to first find out which cells to
3326 * refine/coarsen, then to actually do the refinement and eventually
3327 * transfer the solution vectors between the two different grids. The first
3328 * task is simply achieved by using the well-established Kelly error
3329 * estimator on the temperature (it is the temperature we're mainly
3330 * interested in
for this program, and we need to be accurate in regions of
3332 * diffusion). The
second task is to actually
do the remeshing. That
3333 * involves only basic functions as well, such as the
3334 * <code>refine_and_coarsen_fixed_fraction</code> that refines those cells
3335 * with the largest estimated error that together make up 80 per cent of the
3336 * error, and coarsens those cells with the smallest error that make up
for
3337 * a combined 10 per cent of the error.
3341 * If implemented like
this, we would get a program that will not make much
3342 * progress: Remember that we expect
temperature fields that are nearly
3343 * discontinuous (the diffusivity @f$\kappa@f$ is very small after all) and
3344 * consequently we can expect that a freely adapted mesh will
refine further
3345 * and further into the areas of large
gradients. This decrease in mesh
size
3346 * will then be accompanied by a decrease in time step, requiring an
3347 * exceedingly large number of time steps to solve to a given final time. It
3348 * will also lead to meshes that are much better at resolving
3349 * discontinuities after several mesh refinement cycles than in the
3354 * In particular to prevent the decrease in time step
size and the
3355 * correspondingly large number of time steps, we limit the maximal
3356 * refinement depth of the mesh. To this
end, after the refinement indicator
3357 * has been applied to the cells, we simply
loop over all cells on the
3358 * finest
level and unselect them from refinement
if they would result in
3359 * too high a mesh
level.
3362 *
template <
int dim>
3364 *
BoussinesqFlowProblem<dim>::refine_mesh(
const unsigned
int max_grid_level)
3371 *
temperature_solution,
3372 *
estimated_error_per_cell);
3375 *
estimated_error_per_cell,
3378 *
if (triangulation.n_levels() > max_grid_level)
3380 *
triangulation.active_cell_iterators_on_level(max_grid_level))
3381 *
cell->clear_refine_flag();
3385 * As part of mesh refinement we need to transfer the solution vectors
3386 * from the old mesh to the
new one. To
this end we use the
3388 * should be transferred to the
new grid (we will lose the old grid once
3389 * we have done the refinement so the transfer has to happen concurrently
3390 * with refinement). What we definitely need are the current and the old
3391 *
temperature (BDF-2 time stepping
requires two old solutions). Since the
3394 * structure. Moreover, we choose to transfer the Stokes solution, too,
3395 * since we need the velocity at two previous time steps, of which only
3396 *
one is calculated on the fly.
3402 * handlers. With
this at place, we can prepare the triangulation and the
3403 *
data vectors
for refinement (in
this order).
3406 *
const std::vector<TrilinosWrappers::MPI::Vector> x_temperature = {
3407 *
temperature_solution, old_temperature_solution};
3411 *
temperature_dof_handler);
3413 *
stokes_dof_handler);
3415 *
triangulation.prepare_coarsening_and_refinement();
3416 *
temperature_trans.prepare_for_coarsening_and_refinement(x_temperature);
3417 *
stokes_trans.prepare_for_coarsening_and_refinement(x_stokes);
3421 * Now everything is ready, so
do the refinement and recreate the dof
3422 * structure on the
new grid, and initialize the
matrix structures and the
3423 *
new vectors in the <code>setup_dofs</code> function. Next, we actually
3424 * perform the interpolation of the solutions between the grids. We create
3425 * another
copy of temporary vectors
for temperature (now corresponding to
3426 * the
new grid), and let the
interpolate function
do the job. Then, the
3427 * resulting array of vectors is written into the respective vector member
3432 * Remember that the
set of constraints will be updated
for the
new
3433 * triangulation in the setup_dofs() call.
3436 *
triangulation.execute_coarsening_and_refinement();
3442 *
temperature_trans.interpolate(tmp);
3444 *
temperature_solution = tmp[0];
3445 *
old_temperature_solution = tmp[1];
3449 * After the solution has been transferred we then enforce the constraints
3450 * on the transferred solution.
3453 *
temperature_constraints.distribute(temperature_solution);
3454 *
temperature_constraints.distribute(old_temperature_solution);
3458 * For the Stokes vector, everything is just the same – except that
3459 * we
do not need another temporary vector since we just
interpolate a
3460 * single vector. In the
end, we have to tell the program that the matrices
3461 * and preconditioners need to be regenerated, since the mesh has changed.
3464 *
stokes_trans.interpolate(stokes_solution);
3466 *
stokes_constraints.distribute(stokes_solution);
3468 *
rebuild_stokes_matrix =
true;
3469 *
rebuild_temperature_matrices =
true;
3470 *
rebuild_stokes_preconditioner =
true;
3478 * <a name=
"step_31-BoussinesqFlowProblemrun"></a>
3479 * <h4>BoussinesqFlowProblem::run</h4>
3483 * This function performs all the essential steps in the Boussinesq
3484 * program. It starts by setting up a grid (depending on the spatial
3485 * dimension, we choose some different
level of initial refinement and
3486 * additional adaptive refinement steps, and then create a cube in
3487 * <code>dim</code> dimensions and
set up the dofs
for the
first time. Since
3488 * we want to start the time stepping already with an adaptively refined
3489 * grid, we perform some pre-refinement steps, consisting of all assembly,
3490 * solution and refinement, but without actually advancing in time. Rather,
3491 * we use the vilified <code>goto</code> statement to jump out of the time
3492 * loop right after mesh refinement to start all over again on the
new mesh
3493 * beginning at the <code>start_time_iteration</code> label. (The use of the
3494 * <code>goto</code> is discussed in @ref step_26
"step-26".)
3498 * Before we start, we project the initial values to the grid and obtain the
3499 *
first data for the <code>old_temperature_solution</code> vector. Then, we
3500 * initialize time step number and time step and start the time
loop.
3503 *
template <int dim>
3504 *
void BoussinesqFlowProblem<dim>::run()
3506 * const unsigned int initial_refinement = (dim == 2 ? 4 : 2);
3507 * const unsigned int n_pre_refinement_steps = (dim == 2 ? 4 : 3);
3510 * GridGenerator::hyper_cube(triangulation);
3511 * global_Omega_diameter = GridTools::diameter(triangulation);
3513 * triangulation.refine_global(initial_refinement);
3517 * unsigned int pre_refinement_step = 0;
3519 * start_time_iteration:
3521 * VectorTools::project(temperature_dof_handler,
3522 * temperature_constraints,
3523 * QGauss<dim>(temperature_degree + 2),
3524 * EquationData::TemperatureInitialValues<dim>(),
3525 * old_temperature_solution);
3527 * timestep_number = 0;
3528 * time_step = old_time_step = 0;
3534 * std::cout <<
"Timestep " << timestep_number <<
": t=" << time
3539 * The first steps in the time loop are all obvious – we
3540 * assemble the Stokes system, the preconditioner, the temperature
3541 * matrix (matrices and preconditioner do actually only change in case
3542 * we
've remeshed before), and then do the solve. Before going on with
3543 * the next time step, we have to check whether we should first finish
3544 * the pre-refinement steps or if we should remesh (every fifth time
3545 * step), refining up to a level that is consistent with initial
3546 * refinement and pre-refinement steps. Last in the loop is to advance
3547 * the solutions, i.e., to copy the solutions to the next "older" time
3551 * assemble_stokes_system();
3552 * build_stokes_preconditioner();
3553 * assemble_temperature_matrix();
3559 * std::cout << std::endl;
3561 * if ((timestep_number == 0) &&
3562 * (pre_refinement_step < n_pre_refinement_steps))
3564 * refine_mesh(initial_refinement + n_pre_refinement_steps);
3565 * ++pre_refinement_step;
3566 * goto start_time_iteration;
3568 * else if ((timestep_number > 0) && (timestep_number % 5 == 0))
3569 * refine_mesh(initial_refinement + n_pre_refinement_steps);
3571 * time += time_step;
3572 * ++timestep_number;
3574 * old_stokes_solution = stokes_solution;
3575 * old_old_temperature_solution = old_temperature_solution;
3576 * old_temperature_solution = temperature_solution;
3580 * Do all the above until we arrive at time 100.
3583 * while (time <= 100);
3585 * } // namespace Step31
3592 * <a name="step_31-Thecodemaincodefunction"></a>
3593 * <h3>The <code>main</code> function</h3>
3597 * The main function looks almost the same as in all other programs.
3601 * There is one difference we have to be careful about. This program uses
3602 * Trilinos and, typically, Trilinos is configured so that it can run in
3603 * %parallel using MPI. This doesn't mean that it <i>has</i> to run in
3604 * %parallel, and in fact this program (unlike @ref step_32
"step-32") makes no attempt at
3605 * all to do anything in %parallel using MPI. Nevertheless, Trilinos wants the
3606 * MPI system to be initialized. We do that be creating an object of type
3607 * Utilities::MPI::MPI_InitFinalize that initializes MPI (if available) using
3608 * the arguments given to main() (i.e., <code>argc</code> and
3609 * <code>argv</code>) and de-initializes it again when the object goes out of
3613 * int main(int argc, char *argv[])
3617 * using namespace dealii;
3618 * using namespace Step31;
3620 * Utilities::MPI::MPI_InitFinalize mpi_initialization(
3621 * argc, argv, numbers::invalid_unsigned_int);
3625 * This program can only be run in serial. Otherwise, throw an exception.
3628 * AssertThrow(Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD) == 1,
3630 *
"This program can only be run in serial, use ./step-31"));
3632 * BoussinesqFlowProblem<2> flow_problem;
3633 * flow_problem.run();
3635 * catch (std::exception &exc)
3637 * std::cerr << std::endl
3639 * <<
"----------------------------------------------------"
3641 * std::cerr <<
"Exception on processing: " << std::endl
3642 * << exc.what() << std::endl
3643 * <<
"Aborting!" << std::endl
3644 * <<
"----------------------------------------------------"
3651 * std::cerr << std::endl
3653 * <<
"----------------------------------------------------"
3655 * std::cerr <<
"Unknown exception!" << std::endl
3656 * <<
"Aborting!" << std::endl
3657 * <<
"----------------------------------------------------"
3665@anchor step_31-ResultsSection
3666<a name=
"step_31-Results"></a><h1>Results</h1>
3669<a name=
"step_31-Resultsin2d"></a><h3> Results in 2d </h3>
3672When you run the program in 2d, the output will look something like
3676Number of active cells: 256 (on 5 levels)
3677Number of degrees of freedom: 3556 (2178+289+1089)
3681 Rebuilding Stokes preconditioner...
3683 0 GMRES iterations for Stokes subsystem.
3685 9 CG iterations for temperature.
3686 Temperature range: -0.16687 1.30011
3688Number of active cells: 280 (on 6 levels)
3689Number of degrees of freedom: 4062 (2490+327+1245)
3693 Rebuilding Stokes preconditioner...
3695 0 GMRES iterations for Stokes subsystem.
3697 9 CG iterations for temperature.
3698 Temperature range: -0.0982971 0.598503
3700Number of active cells: 520 (on 7 levels)
3701Number of degrees of freedom: 7432 (4562+589+2281)
3705 Rebuilding Stokes preconditioner...
3707 0 GMRES iterations for Stokes subsystem.
3709 9 CG iterations for temperature.
3710 Temperature range: -0.0551098 0.294493
3712Number of active cells: 1072 (on 8 levels)
3713Number of degrees of freedom: 15294 (9398+1197+4699)
3717 Rebuilding Stokes preconditioner...
3719 0 GMRES iterations for Stokes subsystem.
3721 9 CG iterations for temperature.
3722 Temperature range: -0.0273524 0.156861
3724Number of active cells: 2116 (on 9 levels)
3725Number of degrees of freedom: 30114 (18518+2337+9259)
3729 Rebuilding Stokes preconditioner...
3731 0 GMRES iterations for Stokes subsystem.
3732 Time step: 0.0574449
3733 9 CG iterations for temperature.
3734 Temperature range: -0.014993 0.0738328
3736Timestep 1: t=0.0574449
3739 56 GMRES iterations for Stokes subsystem.
3740 Time step: 0.0574449
3741 9 CG iterations for temperature.
3742 Temperature range: -0.0273934 0.14488
3748In the beginning we refine the mesh several times adaptively and
3749always return to time step zero to restart on the newly refined
3750mesh. Only then do we start the actual time iteration.
3752The program runs for a while. The temperature field for time steps 0,
3753500, 1000, 1500, 2000, 3000, 4000, and 5000 looks like this (note that
3754the color scale used for the temperature is not always the same):
3756<table align=
"center" class=
"doxtable">
3759 <img src=
"https://dealii.org/images/steps/developer/step-31.2d.solution.00.png" alt=
"">
3762 <img src=
"https://dealii.org/images/steps/developer/step-31.2d.solution.01.png" alt=
"">
3765 <img src=
"https://dealii.org/images/steps/developer/step-31.2d.solution.02.png" alt=
"">
3768 <img src=
"https://dealii.org/images/steps/developer/step-31.2d.solution.03.png" alt=
"">
3773 <img src=
"https://dealii.org/images/steps/developer/step-31.2d.solution.04.png" alt=
"">
3776 <img src=
"https://dealii.org/images/steps/developer/step-31.2d.solution.05.png" alt=
"">
3779 <img src=
"https://dealii.org/images/steps/developer/step-31.2d.solution.06.png" alt=
"">
3782 <img src=
"https://dealii.org/images/steps/developer/step-31.2d.solution.07.png" alt=
"">
3787The visualizations shown here were generated using a version of the example
3788which did not enforce the constraints after transferring the mesh.
3790As can be seen, we have three heat sources that heat fluid and
3791therefore produce a buoyancy effect that lets hots pockets of fluid
3792rise up and swirl around. By a chimney effect, the three streams are
3793pressed together by fluid that comes from the outside and wants to
3794join the updraft party. Note that because the fluid is initially at
3795rest, those parts of the fluid that were initially over the sources
3796receive a longer heating time than that fluid that is later dragged
3797over the source by the fully developed flow field. It is therefore
3798hotter, a fact that can be seen in the red tips of the three
3799plumes. Note also the relatively fine features of the flow field, a
3800result of the sophisticated transport stabilization of the temperature
3801equation we have chosen.
3803In addition to the pictures above, the following ones show the
3804adaptive mesh and the flow field at the same time steps:
3806<table align=
"center" class=
"doxtable">
3809 <img src=
"https://dealii.org/images/steps/developer/step-31.2d.grid.00.png" alt=
"">
3812 <img src=
"https://dealii.org/images/steps/developer/step-31.2d.grid.01.png" alt=
"">
3815 <img src=
"https://dealii.org/images/steps/developer/step-31.2d.grid.02.png" alt=
"">
3818 <img src=
"https://dealii.org/images/steps/developer/step-31.2d.grid.03.png" alt=
"">
3823 <img src=
"https://dealii.org/images/steps/developer/step-31.2d.grid.04.png" alt=
"">
3826 <img src=
"https://dealii.org/images/steps/developer/step-31.2d.grid.05.png" alt=
"">
3829 <img src=
"https://dealii.org/images/steps/developer/step-31.2d.grid.06.png" alt=
"">
3832 <img src=
"https://dealii.org/images/steps/developer/step-31.2d.grid.07.png" alt=
"">
3838<a name=
"step_31-Resultsin3d"></a><h3> Results in 3d </h3>
3841The same thing can of course be done in 3d by changing the template
3842parameter to the BoussinesqFlowProblem object in <code>main()</code>
3843from 2 to 3, so that the output now looks like follows:
3847Number of active cells: 64 (on 3 levels)
3848Number of degrees of freedom: 3041 (2187+125+729)
3852 Rebuilding Stokes preconditioner...
3854 0 GMRES iterations for Stokes subsystem.
3856 9 CG iterations for temperature.
3857 Temperature range: -0.675683 4.94725
3859Number of active cells: 288 (on 4 levels)
3860Number of degrees of freedom: 12379 (8943+455+2981)
3864 Rebuilding Stokes preconditioner...
3866 0 GMRES iterations for Stokes subsystem.
3868 9 CG iterations for temperature.
3869 Temperature range: -0.527701 2.25764
3871Number of active cells: 1296 (on 5 levels)
3872Number of degrees of freedom: 51497 (37305+1757+12435)
3876 Rebuilding Stokes preconditioner...
3878 0 GMRES iterations for Stokes subsystem.
3880 10 CG iterations for temperature.
3881 Temperature range: -0.496942 0.847395
3883Number of active cells: 5048 (on 6 levels)
3884Number of degrees of freedom: 192425 (139569+6333+46523)
3888 Rebuilding Stokes preconditioner...
3890 0 GMRES iterations for Stokes subsystem.
3892 10 CG iterations for temperature.
3893 Temperature range: -0.267683 0.497739
3895Timestep 1: t=0.306373
3898 27 GMRES iterations for Stokes subsystem.
3900 10 CG iterations for temperature.
3901 Temperature range: -0.461787 0.958679
3907Visualizing the temperature isocontours at time steps 0,
390850, 100, 150, 200, 300, 400, 500, 600, 700, and 800 yields the
3911<table align=
"center" class=
"doxtable">
3914 <img src=
"https://dealii.org/images/steps/developer/step-31.3d.solution.00.png" alt=
"">
3917 <img src=
"https://dealii.org/images/steps/developer/step-31.3d.solution.01.png" alt=
"">
3920 <img src=
"https://dealii.org/images/steps/developer/step-31.3d.solution.02.png" alt=
"">
3923 <img src=
"https://dealii.org/images/steps/developer/step-31.3d.solution.03.png" alt=
"">
3928 <img src=
"https://dealii.org/images/steps/developer/step-31.3d.solution.04.png" alt=
"">
3931 <img src=
"https://dealii.org/images/steps/developer/step-31.3d.solution.05.png" alt=
"">
3934 <img src=
"https://dealii.org/images/steps/developer/step-31.3d.solution.06.png" alt=
"">
3937 <img src=
"https://dealii.org/images/steps/developer/step-31.3d.solution.07.png" alt=
"">
3942 <img src=
"https://dealii.org/images/steps/developer/step-31.3d.solution.08.png" alt=
"">
3945 <img src=
"https://dealii.org/images/steps/developer/step-31.3d.solution.09.png" alt=
"">
3948 <img src=
"https://dealii.org/images/steps/developer/step-31.3d.solution.10.png" alt=
"">
3955That the first picture looks like three hedgehogs stems from the fact that our
3956scheme essentially projects the source times the first time step size onto the
3957mesh to obtain the temperature field in the first time step. Since the source
3958function is discontinuous, we need to expect over- and undershoots from this
3959project. This is in fact what happens (it
's easier to check this in 2d) and
3960leads to the crumpled appearance of the isosurfaces. The visualizations shown
3961here were generated using a version of the example which did not enforce the
3962constraints after transferring the mesh.
3966<a name="step_31-Numericalexperimentstodetermineoptimalparameters"></a><h3> Numerical experiments to determine optimal parameters </h3>
3969The program as is has three parameters that we don't have much of a
3970theoretical handle on how to choose in an optimal way. These are:
3972 <li>The time step must satisfy a CFL condition
3973 @f$k\le \min_K \frac{c_kh_K}{\|\mathbf{u}\|_{L^\infty(K)}}@f$. Here, @f$c_k@f$ is
3974 dimensionless, but what is the right value?
3975 <li>In the computation of the artificial viscosity,
3980 \|\mathbf{u}\|_{L^\infty(K)}
3984 \frac{\|R_\alpha(T)\|_{L^\infty(K)}}{c(\mathbf{u},T)}
3987 with @f$c(\mathbf{u},T) =
3988 c_R\ \|\mathbf{u}\|_{L^\infty(\Omega)} \ \mathrm{var}(T)
3989 \ |\mathrm{diam}(\Omega)|^{\alpha-2}@f$.
3990 Here, the choice of the dimensionless %numbers @f$\beta,c_R@f$ is of
3993In all of these cases, we will have to expect that the correct choice of each
3994value depends on that of the others, and most likely also on the space
3995dimension and polynomial degree of the finite element used for the
3996temperature. Below we
'll discuss a few numerical experiments to choose
3997constants @f$c_k@f$ and @f$\beta@f$.
3999Below, we will not discuss the choice of @f$c_R@f$. In the program, we set
4000it to @f$c_R=2^{\frac{4-2\alpha}{d}}@f$. The reason for this value is a
4001bit complicated and has more to do with the history of the program
4002than reasoning: while the correct formula for the global scaling
4003parameter @f$c(\mathbf{u},T)@f$ is shown above, the program (including the
4004version shipped with deal.II 6.2) initially had a bug in that we
4007 \|\mathbf{u}\|_{L^\infty(\Omega)} \ \mathrm{var}(T)
4008 \ \frac{1}{|\mathrm{diam}(\Omega)|^{\alpha-2}}@f$ instead, where
4009we had set the scaling parameter to one. Since we only computed on the
4010unit square/cube where @f$\mathrm{diam}(\Omega)=2^{1/d}@f$, this was
4011entirely equivalent to using the correct formula with
4012@f$c_R=\left(2^{1/d}\right)^{4-2\alpha}=2^{\frac{4-2\alpha}{d}}@f$. Since
4013this value for @f$c_R@f$ appears to work just fine for the current
4014program, we corrected the formula in the program and set @f$c_R@f$ to a
4015value that reproduces exactly the results we had before. We will,
4016however, revisit this issue again in @ref step_32 "step-32".
4018Now, however, back to the discussion of what values of @f$c_k@f$ and
4019@f$\beta@f$ to choose:
4022<a name="step_31-Choosingicsubksubiandbeta"></a><h4> Choosing <i>c<sub>k</sub></i> and beta </h4>
4025These two constants are definitely linked in some way. The reason is easy to
4026see: In the case of a pure advection problem,
4027@f$\frac{\partial T}{\partial t} + \mathbf{u}\cdot\nabla T = \gamma@f$, any
4028explicit scheme has to satisfy a CFL condition of the form
4029@f$k\le \min_K \frac{c_k^a h_K}{\|\mathbf{u}\|_{L^\infty(K)}}@f$. On the other hand,
4030for a pure diffusion problem,
4031@f$\frac{\partial T}{\partial t} + \nu \Delta T = \gamma@f$,
4032explicit schemes need to satisfy a condition
4033@f$k\le \min_K \frac{c_k^d h_K^2}{\nu}@f$. So given the form of @f$\nu@f$ above, an
4034advection diffusion problem like the one we have to solve here will result in
4035a condition of the form
4037k\le \min_K \min \left\{
4038 \frac{c_k^a h_K}{\|\mathbf{u}\|_{L^\infty(K)}},
4039 \frac{c_k^d h_K^2}{\beta \|\mathbf{u}\|_{L^\infty(K)} h_K}\right\}
4041 \min_K \left( \min \left\{
4043 \frac{c_k^d}{\beta}\right\}
4044 \frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}} \right)
4046It follows that we have to face the fact that we might want to choose @f$\beta@f$
4047larger to improve the stability of the numerical scheme (by increasing the
4048amount of artificial diffusion), but we have to pay a price in the form of
4049smaller, and consequently more time steps. In practice, one would therefore
4050like to choose @f$\beta@f$ as small as possible to keep the transport problem
4051sufficiently stabilized while at the same time trying to choose the time step
4052as large as possible to reduce the overall amount of work.
4054The find the right balance, the only way is to do a few computational
4055experiments. Here's what we did: We modified the program slightly to allow
4056less mesh refinement (so we don
't always have to wait that long) and to choose
4061 \|\mathbf{u}\|_{L^\infty(K)} h_K
4062@f$ to eliminate the effect of the constant @f$c_R@f$ (we know that
4063solutions are stable by using this version of @f$\nu(T)@f$ as an artificial
4064viscosity, but that we can improve things -- i.e. make the solution
4065sharper -- by using the more complicated formula for this artificial
4066viscosity). We then run the program
4067for different values @f$c_k,\beta@f$ and observe maximal and minimal temperatures
4068in the domain. What we expect to see is this: If we choose the time step too
4069big (i.e. choose a @f$c_k@f$ bigger than theoretically allowed) then we will get
4070exponential growth of the temperature. If we choose @f$\beta@f$ too small, then
4071the transport stabilization becomes insufficient and the solution will show
4072significant oscillations but not exponential growth.
4075<a name="step_31-ResultsforQsub1subelements"></a><h5>Results for Q<sub>1</sub> elements</h5>
4078Here is what we get for
4079@f$\beta=0.01, \beta=0.1@f$, and @f$\beta=0.5@f$, different choices of @f$c_k@f$, and
4080bilinear elements (<code>temperature_degree=1</code>) in 2d:
4082<table align="center" class="doxtable">
4085 <img src="https://dealii.org/images/steps/developer/step-31.timestep.q1.beta=0.01.png" alt="">
4088 <img src="https://dealii.org/images/steps/developer/step-31.timestep.q1.beta=0.03.png" alt="">
4093 <img src="https://dealii.org/images/steps/developer/step-31.timestep.q1.beta=0.1.png" alt="">
4096 <img src="https://dealii.org/images/steps/developer/step-31.timestep.q1.beta=0.5.png" alt="">
4101The way to interpret these graphs goes like this: for @f$\beta=0.01@f$ and
4102@f$c_k=\frac 12,\frac 14@f$, we see exponential growth or at least large
4103variations, but if we choose
4104@f$k=\frac 18\frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}@f$
4105or smaller, then the scheme is
4106stable though a bit wobbly. For more artificial diffusion, we can choose
4107@f$k=\frac 14\frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}@f$
4108or smaller for @f$\beta=0.03@f$,
4109@f$k=\frac 13\frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}@f$
4110or smaller for @f$\beta=0.1@f$, and again need
4111@f$k=\frac 1{15}\frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}@f$
4112for @f$\beta=0.5@f$ (this time because much diffusion requires a small time
4115So how to choose? If we were simply interested in a large time step, then we
4116would go with @f$\beta=0.1@f$ and
4117@f$k=\frac 13\frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}@f$.
4118On the other hand, we're also interested in accuracy and here it may be of
4119interest to actually investigate what these curves show. To this end note that
4120we start with a zero temperature and that our sources are positive — so
4121we would intuitively expect that the temperature can never drop below
4122zero. But it does, a consequence of Gibb
's phenomenon when using continuous
4123elements to approximate a discontinuous solution. We can therefore see that
4124choosing @f$\beta@f$ too small is bad: too little artificial diffusion leads to
4125over- and undershoots that aren't diffused away. On the other hand, for large
4126@f$\beta@f$, the minimum temperature drops below zero at the beginning but then
4127quickly diffuses back to zero.
4129On the other hand, let
's also look at the maximum temperature. Watching the
4130movie of the solution, we see that initially the fluid is at rest. The source
4131keeps heating the same volume of fluid whose temperature increases linearly at
4132the beginning until its buoyancy is able to move it upwards. The hottest part
4133of the fluid is therefore transported away from the solution and fluid taking
4134its place is heated for only a short time before being moved out of the source
4135region, therefore remaining cooler than the initial bubble. If @f$\kappa=0@f$
4136(in the program it is nonzero but very small) then the hottest part of the
4137fluid should be advected along with the flow with its temperature
4138constant. That's what we can see in the graphs with the smallest @f$\beta@f$: Once
4139the maximum temperature is reached, it hardly changes any more. On the other
4140hand, the larger the artificial diffusion, the more the hot spot is
4141diffused. Note that for this criterion, the time step size does not play a
4144So to sum up, likely the best choice would appear to be @f$\beta=0.03@f$
4145and @f$k=\frac 14\frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}@f$. The curve is
4146a bit wobbly, but overall pictures looks pretty reasonable with the
4147exception of some over and undershoots close to the start time due to
4151<a name="step_31-ResultsforQsub2subelements"></a><h5>Results for Q<sub>2</sub> elements</h5>
4154One can repeat the same sequence of experiments for higher order
4155elements as well. Here are the graphs for bi-quadratic shape functions
4156(<code>temperature_degree=2</code>) for the temperature, while we
4157retain the @f$Q_2/Q_1@f$ stable Taylor-Hood element for the Stokes system:
4159<table align="center" class="doxtable">
4162 <img src="https://dealii.org/images/steps/developer/step-31.timestep.q2.beta=0.01.png" alt="">
4165 <img src="https://dealii.org/images/steps/developer/step-31.timestep.q2.beta=0.03.png" alt="">
4170 <img src="https://dealii.org/images/steps/developer/step-31.timestep.q2.beta=0.1.png" alt="">
4175Again, small values of @f$\beta@f$ lead to less diffusion but we have to
4176choose the time step very small to keep things under control. Too
4177large values of @f$\beta@f$ make for more diffusion, but again require
4178small time steps. The best value would appear to be @f$\beta=0.03@f$, as
4179for the @f$Q_1@f$ element, and then we have to choose
4180@f$k=\frac 18\frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}@f$ — exactly
4181half the size for the @f$Q_1@f$ element, a fact that may not be surprising
4182if we state the CFL condition as the requirement that the time step be
4183small enough so that the distance transport advects in each time step
4184is no longer than one <i>grid point</i> away (which for @f$Q_1@f$ elements
4185is @f$h_K@f$, but for @f$Q_2@f$ elements is @f$h_K/2@f$). It turns out that @f$\beta@f$
4186needs to be slightly larger for obtaining stable results also late in
4187the simulation at times larger than 60, so we actually choose it as
4188@f$\beta = 0.034@f$ in the code.
4191<a name="step_31-Resultsfor3d"></a><h5>Results for 3d</h5>
4194One can repeat these experiments in 3d and find the optimal time step
4195for each value of @f$\beta@f$ and find the best value of @f$\beta@f$. What one
4196finds is that for the same @f$\beta@f$ already used in 2d, the time steps
4197needs to be a bit smaller, by around a factor of 1.2 or so. This is
4198easily explained: the time step restriction is
4199@f$k=\min_K \frac{ch_K}{\|\mathbf{u}\|_{L^\infty(K)}}@f$ where @f$h_K@f$ is
4200the <i>diameter</i> of the cell. However, what is really needed is the
4201distance between mesh points, which is @f$\frac{h_K}{\sqrt{d}}@f$. So a
4202more appropriate form would be
4203@f$k=\min_K \frac{ch_K}{\|\mathbf{u}\|_{L^\infty(K)}\sqrt{d}}@f$.
4205The second find is that one needs to choose @f$\beta@f$ slightly bigger
4206(about @f$\beta=0.05@f$ or so). This then again reduces the time step we
4212<a name="step_31-Conclusions"></a><h5>Conclusions</h5>
4215Concluding, from the simple computations above, @f$\beta=0.034@f$ appears to be a
4216good choice for the stabilization parameter in 2d, and @f$\beta=0.05@f$ in 3d. In
4217a dimension independent way, we can model this as @f$\beta=0.017d@f$. If one does
4218longer computations (several thousand time steps) on finer meshes, one
4219realizes that the time step size is not quite small enough and that for
4220stability one will have to reduce the above values a bit more (by about a
4221factor of @f$\frac 78@f$).
4223As a consequence, a formula that reconciles 2d, 3d, and variable polynomial
4224degree and takes all factors in account reads as follows:
4227 \frac 1{2 \cdot 1.7} \frac 1{\sqrt{d}}
4230 \frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}
4232 \frac 1{1.7 d\sqrt{d}}
4234 \frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}.
4236In the first form (in the center of the equation), @f$\frac
42371{2 \cdot 1.7}@f$ is a universal constant, @f$\frac 1{\sqrt{d}}@f$
4238is the factor that accounts for the difference between cell diameter
4239and grid point separation,
4240@f$\frac 2d@f$ accounts for the increase in @f$\beta@f$ with space dimension,
4241@f$\frac 1{q_T}@f$ accounts for the distance between grid points for
4242higher order elements, and @f$\frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}@f$
4243for the local speed of transport relative to the cell size. This is
4244the formula that we use in the program.
4246As for the question of whether to use @f$Q_1@f$ or @f$Q_2@f$ elements for the
4247temperature, the following considerations may be useful: First,
4248solving the temperature equation is hardly a factor in the overall
4249scheme since almost the entire compute time goes into solving the
4250Stokes system in each time step. Higher order elements for the
4251temperature equation are therefore not a significant drawback. On the
4252other hand, if one compares the size of the over- and undershoots the
4253solution produces due to the discontinuous source description, one
4254notices that for the choice of @f$\beta@f$ and @f$k@f$ as above, the @f$Q_1@f$
4255solution dips down to around @f$-0.47@f$, whereas the @f$Q_2@f$ solution only
4256goes to @f$-0.13@f$ (remember that the exact solution should never become
4257negative at all. This means that the @f$Q_2@f$ solution is significantly
4258more accurate; the program therefore uses these higher order elements,
4259despite the penalty we pay in terms of smaller time steps.
4262<a name="step_31-Possibilitiesforextensions"></a><h3> Possibilities for extensions </h3>
4265There are various ways to extend the current program. Of particular interest
4266is, of course, to make it faster and/or increase the resolution of the
4267program, in particular in 3d. This is the topic of the @ref step_32 "step-32"
4268tutorial program which will implement strategies to solve this problem in
4269%parallel on a cluster. It is also the basis of the much larger open
4270source code ASPECT (see https://aspect.geodynamics.org/ ) that can solve realistic
4271problems and that constitutes the further development of @ref step_32 "step-32".
4273Another direction would be to make the fluid flow more realistic. The program
4274was initially written to simulate various cases simulating the convection of
4275material in the earth's mantle, i.e. the zone between the outer earth core and
4276the solid earth crust: there, material is heated from below and cooled from
4277above, leading to thermal convection. The physics of this fluid are much more
4278complicated than shown in this program, however: The viscosity of mantle
4279material is strongly dependent on the temperature, i.e. @f$\eta=\eta(T)@f$, with
4280the dependency frequently modeled as a viscosity that is reduced exponentially
4281with rising temperature. Secondly, much of the dynamics of the mantle is
4282determined by chemical reactions, primarily phase changes of the various
4283crystals that make up the mantle; the buoyancy term on the right hand side of
4284the Stokes equations then depends not only on the temperature, but also on the
4285chemical composition at a given location which is advected by the flow field
4286but also changes as a function of pressure and temperature. We will
4287investigate some of these effects in later tutorial programs as well.
4290<a name=
"step_31-PlainProg"></a>
4291<h1> The plain program</h1>
4292@include
"step-31.cc"
* * for(const auto &cell :triangulation.active_cell_iterators())
* * int main(int argc, char **argv)
* x_component_mask set(0, true)
* * * struct InterferenceTaperTransform *
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, level_dof_access > > &cell)
const std::vector< Point< dim > > & get_unit_support_points() const
const unsigned int n_components
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
virtual void vector_value(const Point< dim > &p, Vector< RangeNumberType > &values) const
static void estimate(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Quadrature< dim - 1 > &quadrature, const std::map< types::boundary_id, const Function< spacedim, Number > * > &neumann_bc, const ReadVector< Number > &solution, Vector< float > &error, const ComponentMask &component_mask={}, const Function< spacedim > *coefficients=nullptr, const unsigned int n_threads=numbers::invalid_unsigned_int, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id, const types::material_id material_id=numbers::invalid_material_id, const Strategy strategy=cell_diameter_over_24)
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcMessage(std::string arg1)
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
@ update_hessians
Second derivatives of shape functions.
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
IndexSet complete_index_set(const IndexSet::size_type N)
std::vector< index_type > data
void approximate(const SynchronousIterators< std::tuple< typename DoFHandler< dim, spacedim >::active_cell_iterator, Vector< float >::iterator > > &cell, const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof_handler, const InputVector &solution, const unsigned int component)
Expression sign(const Expression &x)
void downstream(DoFHandler< dim, spacedim > &dof_handler, const Tensor< 1, spacedim > &direction, const bool dof_wise_renumbering=false)
void refine(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold, const unsigned int max_to_mark=numbers::invalid_unsigned_int)
void refine_and_coarsen_fixed_fraction(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double top_fraction, const double bottom_fraction, const unsigned int max_n_cells=std::numeric_limits< unsigned int >::max(), const VectorTools::NormType norm_type=VectorTools::L1_norm)
@ matrix
Contents is actually a matrix.
constexpr types::blas_int zero
constexpr types::blas_int one
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
* * * ScaleZFunction< dim, Number, components >::ScaleZFunction * component(component)
* * if(update_pressure &update_flags) * compute_pressure(constitutive_request
* const Number temperature
* * * * std::vector< Number > ThermoPlasticMaterial< dim, ViscoplasticYieldLaw, Number >::get_state_parameters const
constexpr ReturnType< rank, T >::value_type & extract(T &t, const ArrayType &indices)
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
unsigned int n_active_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
void copy(const T *begin, const T *end, U *dest)
int(&) functions(const void *v1, const void *v2)
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
1.9.8