This program was contributed by Martin Kronbichler and Wolfgang Bangerth.
This material is based upon work partly supported by the National Science Foundation under Award No. EAR-0426271 and The California Institute of Technology. Any opinions, findings, and conclusions or recommendations expressed in this publication are those of the author and do not necessarily reflect the views of the National Science Foundation or of The California Institute of Technology.
This program deals with an interesting physical problem: how does a fluid (i.e. a liquid or gas) behave if it experiences differences in buoyancy caused by temperature differences? It is clear that those parts of the fluid that are hotter (and therefore lighter) are going to rise up and those that are cooler (and denser) are going to sink down with gravity.
In cases where the fluid moves slowly enough such that inertia effects can be neglected, the equations that describe such behavior are the Boussinesq equations that read as follows:
These equations fall into the class of vector-valued problems (a toplevel overview of this topic can be found in the Handling vector valued problems module). Here, u is the velocity field, p the pressure, and T the temperature of the fluid.
is the symmetric gradient of the velocity. As can be seen, velocity and pressure solve a Stokes equation describing the motion of an incompressible fluid, an equation we have previously considered in step-22; we will draw extensively on the experience we have gained in that program, in particular with regard to efficient linear Stokes solvers.
The forcing term of the fluid motion is the buoyancy of the fluid, expressed as the product of the Rayleigh number
, the temperature T and the gravity vector g. (A possibly more intuitive formulation would use
as right hand side where
is the average temperature, and the right hand side then describes the forces due to local deviations from the average density; this formulation is entirely equivalent if the gravity vector results from a gravity potential
, i.e.
, and yields the exact same solution except for the pressure which will now be
.)
While the first two equations describe how the fluid reacts to temperature differences by moving around, the third equation states how the fluid motion affects the temperature field: it is an advection diffusion equation, i.e. the temperature is attached to the fluid particles and advected along in the flow field, with an additional diffusion (heat conduction) term. In many applications, the diffusion coefficient is fairly small, and the temperature equation is in fact transport, not diffusion dominated and therefore in character more hyperbolic than elliptic; we will have to take this into account when developing a stable discretization.
In the equations above, the term
on the right hand side denotes the heat sources and may be a spatially and temporally varying function.
and
denote the viscosity and diffusivity coefficients, which we assume constant for this tutorial program. The more general case when
depends on the temperature is an important factor in physical applications: Most materials become more fluid as they get hotter (i.e.,
decreases with T); sometimes, as in the case of rock minerals at temperatures close to their melting point,
may change by orders of magnitude over the typical range of temperatures.
, called the Rayleigh number, is a dimensionless number that describes the ratio of heat transport due to convection induced by buoyancy changes from temperature differences, and of heat transport due to thermal diffusion. A small Rayleigh number implies that buoyancy is not strong relative to viscosity and fluid motion u is slow enough so that heat diffusion
is the dominant heat transport term. On the other hand, a fluid with a high Rayleigh number will show vigorous convection that dominates heat conduction.
For most fluids for which we are interested in computing thermal convection, the Rayleigh number is very large, often
or larger. From the structure of the equations, we see that this will lead to large pressure differences and large velocities. Consequently, the convection term in the convection-diffusion equation for T will also be very large and an accurate solution of this equation will require us to choose small time steps. Problems with large Rayleigh numbers are therefore hard to solve numerically for similar reasons that make solving the Navier-Stokes equations hard to solve when the Reynolds number
is large.
Note that a large Rayleigh number does not necessarily involve large velocities in absolute terms. For example, the Rayleigh number in the earth mantle has a Rayleigh number larger than
. Yet the velocities are small: the material is in fact solid rock but it is so hot and under pressure that it can flow very slowly, on the order of at most a few centimeters per year. Nevertheless, this can lead to mixing over time scales of many million years, a time scale much shorter than for the same amount of heat to be distributed by thermal conductivity and a time scale of relevance to affect the evolution of the earth's interior and surface structure.
Since the Boussinesq equations are derived under the assumption that inertia of the fluid's motion does not play a role, the flow field is at each time entirely determined by buoyancy difference at that time, not by the flow field at previous times. This is reflected by the fact that the first two equations above are the steady state Stokes equation that do not contain a time derivative. Consequently, we do not need initial conditions for either velocities or pressure. On the other hand, the temperature field does satisfy an equation with a time derivative, so we need initial conditions for T.
As for boundary conditions: if
then the temperature satisfies a second order differential equation that requires boundary data all around the boundary for all times. These can either be a prescribed boundary temperature
(Dirichlet boundary conditions), or a prescribed thermal flux
; in this program, we will use an insulated boundary condition, i.e. prescribe no thermal flux:
.
Similarly, the velocity field requires us to pose boundary conditions. These may be no-slip no-flux conditions u=0 on
if the fluid sticks to the boundary, or no normal flux conditions
if the fluid can flow along but not across the boundary, or any number of other conditions that are physically reasonable. In this program, we will use no normal flux conditions.
Like the equations solved in step-21, we here have a system of differential-algebraic equations (DAE): with respect to the time variable, only the temperature equation is a differential equation whereas the Stokes system for u and p has no time-derivatives and is therefore of the sort of an algebraic constraint that has to hold at each time instant. The main difference to step-21 is that the algebraic constraint there was a mixed Laplace system of the form
where now we have a Stokes system
where
is an operator similar to the Laplacian
applied to a vector field.
Given the similarity to what we have done in step-21, it may not come as a surprise that we choose a similar approach, although we will have to make adjustments for the change in operator in the top-left corner of the differential operator.
The structure of the problem as a DAE allows us to use the same strategy as we have already used in step-21, i.e. we use a time lag scheme: we first solve the temperature equation (using an extrapolated velocity field), and then insert the new temperature solution into the right hand side of the velocity equation. The way we implement this in our code looks at things from a slightly different perspective, though. We first solve the Stokes equations for velocity and pressure using the temperature field from the previous time step, which means that we get the velocity for the previous time step. In other words, we first solve the Stokes system for time step n-1 as
and then the temperature equation with an extrapolated velocity field to time n.
In contrast to step-21, we'll use a higher order time stepping scheme here, namely the Backward Differentiation Formula scheme of order 2 (BDF-2 in short) that replaces the time derivative
by the (one-sided) difference quotient
with k the time step size. This gives the discretized-in-time temperature equation
Note how the temperature equation is solved semi-explicitly: diffusion is treated implicitly whereas advection is treated explicitly using an extrapolation (or forward-projection) of temperature and velocity, including the just-computed velocity
. The forward-projection to the current time level n is derived from a Taylor expansion,
. We need this projection for maintaining the order of accuracy of the BDF-2 scheme. In other words, the temperature fields we use in the explicit right hand side are second order approximations of the current temperature field — not quite an explicit time stepping scheme, but by character not too far away either.
The introduction of the temperature extrapolation limits the time step by a Courant-Friedrichs-Lewy (CFL) condition just like it was in step-21. (We wouldn't have had that stability condition if we treated the advection term implicitly since the BDF-2 scheme is A-stable, at the price that we needed to build a new temperature matrix at each time step.) We will discuss the exact choice of time step in the results section, but for the moment of importance is that this CFL condition means that the time step size k may change from time step to time step, and that we have to modify the above formula slightly. If
are the time steps sizes of the current and previous time step, then we use the approximations
and
, and above equation is generalized as follows:
where
denotes the extrapolation of velocity u and temperature T to time level n, using the values at the two previous time steps. That's not an easy to read equation, but will provide us with the desired higher order accuracy. As a consistency check, it is easy to verify that it reduces to the same equation as above if
.
As a final remark we note that the choice of a higher order time stepping scheme of course forces us to keep more time steps in memory; in particular, we here will need to have
around, a vector that we could previously discard. This seems like a nuisance that we were able to avoid previously by using only a first order time stepping scheme, but as we will see below when discussing the topic of stabilization, we will need this vector anyway and so keeping it around for time discretization is essentially for free and gives us the opportunity to use a higher order scheme.
Like solving the mixed Laplace equations, solving the Stokes equations requires us to choose particular pairs of finite elements for velocities and pressure variables. Because this has already been discussed in step-22, we only cover this topic briefly: Here, we use the stable pair
. These are continuous elements, so we can form the weak form of the Stokes equation without problem by integrating by parts and substituting continuous functions by their discrete counterparts:
for all test functions
. The first term of the first equation is considered as the inner product between tensors, i.e.
. Because the second tensor in this product is symmetric, the anti-symmetric component of
plays no role and it leads to the entirely same form if we use the symmetric gradient of
instead. Consequently, the formulation we consider and that we implement is
This is exactly the same as what we already discussed in step-22 and there is not much more to say about this here.
The more interesting question is what to do with the temperature advection-diffusion equation. By default, not all discretizations of this equation are equally stable unless we either do something like upwinding, stabilization, or all of this. One way to achieve this is to use discontinuous elements (i.e. the FE_DGQ class that we used, for example, in the discretization of the transport equation in step-12, or in discretizing the pressure in step-20 and step-21) and to define a flux at the interface between cells that takes into account upwinding. If we had a pure advection problem this would probably be the simplest way to go. However, here we have some diffusion as well, and the discretization of the Laplace operator with discontinuous elements is cumbersome because of the significant number of additional terms that need to be integrated on each face between cells. Discontinuous elements also have the drawback that the use of numerical fluxes introduces an additional numerical diffusion that acts everywhere, whereas we would really like to minimize the effect of numerical diffusion to a minimum and only apply it where it is necessary to stabilize the scheme.
A better alternative is therefore to add some nonlinear viscosity to the model. Essentially, what this does is to transform the temperature equation from the form
to something like
where
is an addition viscosity (diffusion) term that only acts in the vicinity of shocks and other discontinuities.
is chosen in such a way that if T satisfies the original equations, the additional viscosity is zero.
To achieve this, the literature contains a number of approaches. We will here follow one developed by Guermond and Popov that builds on a suitably defined residual and a limiting procedure for the additional viscosity. To this end, let us define a residual
as follows:
where we will later choose the stabilization exponent
from within the range
. Note that
will be zero if
satisfies the temperature equation, since then the term in parentheses will be zero. Multiplying terms out, we get the following, entirely equivalent form:
With this residual, we can now define the artificial viscosity as a piecewise constant function defined on each cell
with diameter
separately as follows:
Here,
is a stabilization constant (a dimensional analysis reveals that it is unitless and therefore independent of scaling; we will discuss its choice in the results section) and
is a normalization constant that must have units
. We will choose it as
, where
is the range of present temperature values (remember that buoyancy is driven by temperature variations, not the absolute temperature) and
is a dimensionless constant. To understand why this method works consider this: If on a particular cell
the temperature field is smooth, then we expect the residual to be small there (in fact to be on the order of
) and the stabilization term that injects artificial diffusion will there be of size
— i.e. rather small, just as we hope it to be when no additional diffusion is necessary. On the other hand, if we are on or close to a discontinuity of the temperature field, then the residual will be large; the minimum operation in the definition of
will then ensure that the stabilization has size
— the optimal amount of artificial viscosity to ensure stability of the scheme.
It is certainly a good questions whether this scheme really works? Computations by Guermond and Popov have shown that this form of stabilization actually performs much better than most of the other stabilization schemes that are around (for example streamline diffusion, to name only the simplest one). Furthermore, for
they can even prove that it produces better convergence orders for the linear transport equation than for example streamline diffusion. For
, no theoretical results are currently available, but numerical tests indicate that the results are considerably better than for
.
A more practical question is how to introduce this artificial diffusion into the equations we would like to solve. Note that the numerical viscosity
is temperature-dependent, so the equation we want to solve is nonlinear in T — not what one desires from a simple method to stabilize an equation, and even less so if we realize that
is non-differentiable in T. However, there is no reason to despair: we still have to discretize in time and we can treat the term explicitly.
In the definition of the stabilization parameter, we approximate the time derivative by
. This approximation makes only use of available time data and this is the reason why we need to store data of two previous time steps (which enabled us to use the BDF-2 scheme without additional storage cost). We could now simply evaluate the rest of the terms at
, but then the discrete residual would be nothing else than a backward Euler approximation, which is only first order accurate. So, in case of smooth solutions, the residual would be still of the order h, despite the second order time accuracy in the outer BDF-2 scheme and the spatial FE discretization. This is certainly not what we want to have (in fact, we desired to have small residuals in regions where the solution behaves nicely), so a bit more care is needed. The key to this problem is to observe that the first derivative as we constructed it is actually centered at
. We get the desired second order accurate residual calculation if we evaluate all spatial terms at
by using the approximation
, which means that we calculate the nonlinear viscosity as a function of this intermediate temperature,
. Note that this evaluation of the residual is nothing else than a Crank-Nicholson scheme, so we can be sure that now everything is alright. One might wonder whether it is a problem that the numerical viscosity now is not evaluated at time n (as opposed to the rest of the equation). However, this offset is uncritical: For smooth solutions,
will vary continuously, so the error in time offset is k times smaller than the nonlinear viscosity itself, i.e., it is a small higher order contribution that is left out. That's fine because the term itself is already at the level of discretization error in smooth regions.
Using the BDF-2 scheme introduced above, this yields for the simpler case of uniform time steps of size k:
On the left side of this equation remains the term from the time derivative and the original (physical) diffusion which we treat implicitly (this is actually a nice term: the matrices that result from the left hand side are the mass matrix and a multiple of the Laplace matrix — both are positive definite and if the time step size k is small, the sum is simple to invert). On the right hand side, the terms in the first line result from the time derivative; in the second line is the artificial diffusion at time
; the third line contains the advection term, and the fourth the sources. Note that the artificial diffusion operates on the extrapolated temperature at the current time in the same way as we have discussed the advection works in the section on time stepping.
The form for non-uniform time steps that we will have to use in reality is a bit more complicated (which is why we showed the simpler form above first) and reads:
After settling all these issues, the weak form follows naturally from the strong form shown in the last equation, and we immediately arrive at the weak form of the discretized equations:
for all discrete test functions
. Here, the diffusion term has been integrated by parts, and we have used that we will impose no thermal flux,
.
This then results in a matrix equation of form
which given the structure of matrix on the left (the sum of two positive definite matrices) is easily solved using the Conjugate Gradient method.
As explained above, our approach to solving the joint system for velocities/pressure on the one hand and temperature on the other is to use an operator splitting where we first solve the Stokes system for the velocities and pressures using the old temperature field, and then solve for the new temperature field using the just computed velocity field.
Solving the linear equations coming from the Stokes system has been discussed in great detail in step-22. In particular, in the results section of that program, we have discussed a number of alternative linear solver strategies that turned out to be more efficient than the original approach. The best alternative identified there we to use a GMRES solver preconditioned by a block matrix involving the Schur complement. Specifically, the Stokes operator leads to a block structured matrix
and as discussed there a good preconditioner is
where S is the Schur complement of the Stokes operator
. Of course, this preconditioner is not useful because we can't form the various inverses of matrices, but we can use the following as a preconditioner:
where
are approximations to the inverse matrices. In particular, it turned out that S is spectrally equivalent to the mass matrix and consequently replacing
by a CG solver applied to the mass matrix on the pressure space was a good choice. In a small deviation from step-22, we here have a coefficient
in the momentum equation, and by the same derivation as there we should arrive at the conclusion that it is the weighted mass matrix with entries
that we should be using.
It was more complicated to come up with a good replacement
, which corresponds to the discretized symmetric Laplacian of the vector-valued velocity field, i.e.
. In step-22 we used a sparse LU decomposition (using the SparseDirectUMFPACK class) of A for
— the perfect preconditioner — in 2d, but for 3d memory and compute time is not usually sufficient to actually compute this decomposition; consequently, we only use an incomplete LU decomposition (ILU, using the SparseILU class) in 3d.
For this program, we would like to go a bit further. To this end, note that the symmetrized bilinear form on vector fields,
is not too far away from the nonsymmetrized version,
(note that the factor 2 has disappeared in this form). The latter, however, has the advantage that the dim vector components of the test functions are not coupled (well, almost, see below), i.e. the resulting matrix is block-diagonal: one block for each vector component, and each of these blocks is equal to the Laplace matrix for this vector component. So assuming we order degrees of freedom in such a way that first all x-components of the velocity are numbered, then the y-components, and then the z-components, then the matrix
that is associated with this slightly different bilinear form has the form
where
is a Laplace matrix of size equal to the number of shape functions associated with each component of the vector-valued velocity. With this matrix, one could be tempted to define our preconditioner for the velocity matrix A as follows:
where
is a preconditioner for the Laplace matrix — something where we know very well how to build good preconditioners!
In reality, the story is not quite as simple: To make the matrix
definite, we need to make the individual blocks
definite by applying boundary conditions. One can try to do so by applying Dirichlet boundary conditions all around the boundary, and then the so-defined preconditioner
turns out to be a good preconditioner for A if the latter matrix results from a Stokes problem where we also have Dirichlet boundary conditions on the velocity components all around the domain, i.e. if we enforce u=0.
Unfortunately, this "if" is an "if and only if": in the program below we will want to use no-flux boundary conditions of the form
(i.e. flow parallel to the boundary is allowed, but no flux through the boundary). In this case, it turns out that the block diagonal matrix defined above is not a good preconditioner because it neglects the coupling of components at the boundary. A better way to do things is therefore if we build the matrix
as the vector Laplace matrix
and then apply the same boundary condition as we applied to A. If this is Dirichlet boundary conditions all around the domain, the
will decouple to three diagonal blocks as above, and if the boundary conditions are of the form
then this will introduce a coupling of degrees of freedom at the boundary but only there. This, in fact, turns out to be a much better preconditioner than the one introduced above, and has almost all the benefits of what we hoped to get.
To sum this whole story up, we can observe:
and its preconditioner. While this means that we have to store an additional matrix we didn't need before, the preconditioner
is likely going to need much less memory than storing a preconditioner for the coupled matrix A. This is because the matrix
has only a third of the entries per row for all rows corresponding to interior degrees of freedom, and contains coupling between vector components only on those parts of the boundary where the boundary conditions introduce such a coupling. Storing the matrix is therefore comparatively cheap, and we can expect that computing and storing the preconditioner
will also be much cheaper compared to doing so for the fully coupled matrix.
This is the easy part: The matrix for the temperature equation has the form
, where
are mass and stiffness matrices on the temperature space, and
are constants related the time stepping scheme and the current and previous time step. This being the sum of a symmetric positive definite and a symmetric positive semidefinite matrix, the result is also symmetric positive definite. Furthermore,
is a number proportional to the time step, and so becomes small whenever the mesh is fine, damping the effect of the then ill-conditioned stiffness matrix.
As a consequence, inverting this matrix with the Conjugate Gradient algorithm, using a simple preconditioner, is trivial and very cheap compared to inverting the Stokes matrix.
One of the things worth explaining up front about the program below is the use of two different DoFHandler objects. If one looks at the structure of the equations above and the scheme for their solution, one realizes that there is little commonality that keeps the Stokes part and the temperature part together. In all previous tutorial programs in which we have discussed vector-valued problems we have always only used a single finite element with several vector components, and a single DoFHandler object. Sometimes, we have substructured the resulting matrix into blocks to facilitate particular solver schemes; this was, for example, the case in the step-22 program for the Stokes equations upon which the current program is based.
We could of course do the same here. The linear system that we would get would look like this:
The problem with this is: We never use the whole matrix at the same time. In fact, it never really exists at the same time: As explained above,
and
depend on the already computed solution
, in the first case through the time step (that depends on
because it has to satisfy a CFL condition). So we can only assemble it once we've already solved the top left
block Stokes system, and once we've moved on to the temperature equation we don't need the Stokes part any more. Furthermore, we don't actually build the matrix
: Because by the time we get to the temperature equation we already know
, and because we have to assemble the right hand side
at this time anyway, we simply move the term
to the right hand side and assemble it along with all the other terms there. What this means is that there does not remain a part of the matrix where temperature variables and Stokes variables couple, and so a global enumeration of all degrees of freedom is no longer important: It is enough if we have an enumeration of all Stokes degrees of freedom, and of all temperature degrees of freedom independently.
In essence, there is consequently not much use in putting everything into a block matrix (though there are of course the same good reasons to do so for the
Stokes part), or, for that matter, in putting everything into the same DoFHandler object.
But are there downsides to doing so? These exist, though they may not be obvious at first. The main problem is that if we need to create one global finite element that contains velocity, pressure, and temperature shape functions, and use this to initialize the DoFHandler. But we also use this finite element object to initialize all FEValues or FEFaceValues objects that we use. This may not appear to be that big a deal, but imagine what happens when, for example, we evaluate the residual
that we need to compute the artificial viscosity
. For this, we need the Laplacian of the temperature, which we compute using the tensor of second derivatives (Hessians) of the shape functions (we have to give the update_hessians flag to the FEValues object for this). Now, if we have a finite that contains the shape functions for velocities, pressures, and temperatures, that means that we have to compute the Hessians of all shape functions, including the many higher order shape functions for the velocities. That's a lot of computations that we don't need, and indeed if one were to do that (as we had in an early version of the program), assembling the right hand side took about a quarter of the overall compute time.
So what we will do is to use two different finite element objects, one for the Stokes components and one for the temperatures. With this come two different DoFHandlers, two sparsity patterns and two matrices for the Stokes and temperature parts, etc. And whenever we have to assemble something that contains both temperature and Stokes shape functions (in particular the right hand sides of Stokes and temperature equations), then we use two FEValues objects initialized with two cell iterators that we walk in parallel through the two DoFHandler objects associated with the same Triangulation object; for these two FEValues objects, we use of course the same quadrature objects so that we can iterate over the same set of quadrature points, but each FEValues object will get update flags only according to what it actually needs to compute. In particular, when we compute the residual as above, we only ask for the values of the Stokes shape functions, but also the Hessians of the temperature shape functions — much cheaper indeed, and as it turns out: assembling the right hand side of the temperature equation is now a component of the program that is hardly measurable.
With these changes, timing the program yields that only the following operations are relevant for the overall run time:
BoussinesqFlowProblem::setup_dofs: 7% of overall run time. In much the same way as we used PETSc to support our linear algebra needs in step-17 and step-18, we use interfaces to the Trilinos library (see Trilinos) in this program. Trilinos is a very large collection of everything that has to do with linear and nonlinear algebra, as well as all sorts of tools around that (and looks like it will grow in many other directions in the future as well).
The main reason for using Trilinos, similar to our exploring PETSc, is that it is a very powerful library that provides a lot more tools than deal.II's own linear algebra library. That includes, in particular, the ability to work in parallel on a cluster, using MPI, and a wider variety of preconditioners. In the latter class, one of the most interesting capabilities is the existence of the Trilinos ML package that implements an Algebraic Multigrid (AMG) method. We will use this preconditioner to precondition the second order operator part of the momentum equation. The ability to solve problems in parallel will be explored in step-32, using the same problem as discussed here.
The case we want to solve here is as follows: we solve the Boussinesq equations described above with
, i.e. a relatively slow moving fluid that has virtually no thermal diffusive conductivity and transports heat mainly through convection. On the boundary, we will require no-normal flux for the velocity (
) and for the temperature (
). This is one of the cases discussed in the introduction of step-22 and fixes one component of the velocity while allowing flow to be parallel to the boundary. There remain dim-1 components to be fixed, namely the tangential components of the normal stress; for these, we choose homogenous conditions which means that we do not have to anything special. Initial conditions are only necessary for the temperature field, and we choose it to be constant zero.
The evolution of the problem is then entirely driven by the right hand side
of the temperature equation, i.e. by heat sources and sinks. Here, we choose a setup invented in advance of a Christmas lecture: real candles are of course prohibited in U.S. class rooms, but virtual ones are allowed. We therefore choose three spherical heat sources unequally spaced close to the bottom of the domain, imitating three candles. The fluid located at these sources, initially at rest, is then heated up and as the temperature rises gains buoyancy, rising up; more fluid is dragged up and through the sources, leading to three hote plumes that rise up until they are captured by the recirculation of fluid that sinks down on the outside, replacing the air that rises due to heating.
The first step, as always, is to include the functionality of these well-known deal.II library files and some C++ header files.
#include <base/quadrature_lib.h> #include <base/logstream.h> #include <base/utilities.h> #include <lac/full_matrix.h> #include <lac/solver_gmres.h> #include <lac/solver_cg.h> #include <grid/tria.h> #include <grid/grid_generator.h> #include <grid/tria_accessor.h> #include <grid/tria_iterator.h> #include <grid/grid_tools.h> #include <grid/grid_refinement.h> #include <dofs/dof_handler.h> #include <dofs/dof_renumbering.h> #include <dofs/dof_accessor.h> #include <dofs/dof_tools.h> #include <dofs/dof_constraints.h> #include <fe/fe_q.h> #include <fe/fe_system.h> #include <fe/fe_values.h> #include <numerics/vectors.h> #include <numerics/data_out.h> #include <numerics/error_estimator.h> #include <numerics/solution_transfer.h>
Then we need to include some header files that provide vector, matrix, and preconditioner classes that implement interfaces to the respective Trilinos classes. In particular, we will need interfaces to the matrix and vector classes based on Trilinos as well as Trilinos preconditioners:
#include <lac/trilinos_sparse_matrix.h> #include <lac/trilinos_block_sparse_matrix.h> #include <lac/trilinos_vector.h> #include <lac/trilinos_block_vector.h> #include <lac/trilinos_precondition.h>
Finally, here are two C++ headers that haven't been included yet by one of the aforelisted header files:
#include <fstream> #include <sstream>
At the end of this top-matter, we import all deal.II names into the global namespace:
using namespace dealii;
Again, the next stage in the program is the definition of the equation data, that is, the various boundary conditions, the right hand sides and the initial condition (remember that we're about to solve a time-dependent system). The basic strategy for this definition is the same as in step-22. Regarding the details, though, there are some differences.
The first thing is that we don't set any non-homogenous boundary conditions on the velocity, since as is explained in the introduction we will use no-flux conditions
. So what is left are dim-1 conditions for the tangential part of the normal component of the stress tensor,
; we assume homogenous values for these components, i.e. a natural boundary condition that requires no specific action (it appears as a zero term in the right hand side of the weak form).
For the temperature T, we assume no thermal energy flux, i.e.
. This, again, is a boundary condition that does not require us to do anything in particular.
Secondly, we have to set initial conditions for the temperature (no initial conditions are required for the velocity and pressure, since the Stokes equations for the quasi-stationary case we consider here have no time derivatives of the velocity or pressure). Here, we choose a very simple test case, where the initial temperature is zero, and all dynamics are driven by the temperature right hand side.
Thirdly, we need to define the right hand side of the temperature equation. We choose it to be constant within three circles (or spheres in 3d) somewhere at the bottom of the domain, as explained in the introduction, and zero outside.
Finally, or maybe firstly, at the top of this namespace, we define the various material constants we need (
and the Rayleigh number
):
namespace EquationData { const double eta = 1; const double kappa = 1e-6; const double Rayleigh_number = 10; template <int dim> class TemperatureInitialValues : public Function<dim> { public: TemperatureInitialValues () : Function<dim>(1) {} virtual double value (const Point<dim> &p, const unsigned int component = 0) const; virtual void vector_value (const Point<dim> &p, Vector<double> &value) const; }; template <int dim> double TemperatureInitialValues<dim>::value (const Point<dim> &, const unsigned int) const { return 0; } template <int dim> void TemperatureInitialValues<dim>::vector_value (const Point<dim> &p, Vector<double> &values) const { for (unsigned int c=0; c<this->n_components; ++c) values(c) = TemperatureInitialValues<dim>::value (p, c); } template <int dim> class TemperatureRightHandSide : public Function<dim> { public: TemperatureRightHandSide () : Function<dim>(1) {} virtual double value (const Point<dim> &p, const unsigned int component = 0) const; virtual void vector_value (const Point<dim> &p, Vector<double> &value) const; }; template <int dim> double TemperatureRightHandSide<dim>::value (const Point<dim> &p, const unsigned int component) const { Assert (component == 0, ExcMessage ("Invalid operation for a scalar function.")); Assert ((dim==2) || (dim==3), ExcNotImplemented()); static const Point<dim> source_centers[3] = { (dim == 2 ? Point<dim>(.3,.1) : Point<dim>(.3,.5,.1)), (dim == 2 ? Point<dim>(.45,.1) : Point<dim>(.45,.5,.1)), (dim == 2 ? Point<dim>(.75,.1) : Point<dim>(.75,.5,.1)) }; static const double source_radius = (dim == 2 ? 1./32 : 1./8); return ((source_centers[0].distance (p) < source_radius) || (source_centers[1].distance (p) < source_radius) || (source_centers[2].distance (p) < source_radius) ? 1 : 0); } template <int dim> void TemperatureRightHandSide<dim>::vector_value (const Point<dim> &p, Vector<double> &values) const { for (unsigned int c=0; c<this->n_components; ++c) values(c) = TemperatureRightHandSide<dim>::value (p, c); } }
This section introduces some objects that are used for the solution of the linear equations of the Stokes system that we need to solve in each time step. Many of the ideas used here are the same as in step-20, where Schur complement based preconditioners and solvers have been introduced, with the actual interface taken from step-22 (in particular the discussion in the "Results" section of step-22, in which we introduce alternatives to the direct Schur complement approach). Note, however, that here we don't use the Schur complement to solve the Stokes equations, though an approximate Schur complement (the mass matrix on the pressure space) appears in the preconditioner.
namespace LinearSolvers
{
InverseMatrix class template
This class is an interface to calculate the action of an "inverted" matrix on a vector (using the vmult operation) in the same way as the corresponding class in step-22: when the product of an object of this class is requested, we solve a linear equation system with that matrix using the CG method, accelerated by a preconditioner of (templated) class Preconditioner.
In a minor deviation from the implementation of the same class in step-22 (and step-20), we make the vmult function take any kind of vector type (it will yield compiler errors, however, if the matrix does not allow a matrix-vector product with this kind of vector).
Secondly, we catch any exceptions that the solver may have thrown. The reason is as follows: When debugging a program like this one occasionally makes a mistake of passing an indefinite or non-symmetric matrix or preconditioner to the current class. The solver will, in that case, not converge and throw a run-time exception. If not caught here it will propagate up the call stack and may end up in main() where we output an error message that will say that the CG solver failed. The question then becomes: Which CG solver? The one that inverted the mass matrix? The one that inverted the top left block with the Laplace operator? Or a CG solver in one of the several other nested places where we use linear solvers in the current code? No indication about this is present in a run-time exception because it doesn't store the stack of calls through which we got to the place where the exception was generated.
So rather than letting the exception propagate freely up to main() we realize that there is little that an outer function can do if the inner solver fails and rather convert the run-time exception into an assertion that fails and triggers a call to abort(), allowing us to trace back in a debugger how we got to the current place.
template <class Matrix, class Preconditioner> class InverseMatrix : public Subscriptor { public: InverseMatrix (const Matrix &m, const Preconditioner &preconditioner); template <typename VectorType> void vmult (VectorType &dst, const VectorType &src) const; private: const SmartPointer<const Matrix> matrix; const Preconditioner &preconditioner; }; template <class Matrix, class Preconditioner> InverseMatrix<Matrix,Preconditioner>:: InverseMatrix (const Matrix &m, const Preconditioner &preconditioner) : matrix (&m), preconditioner (preconditioner) {} template <class Matrix, class Preconditioner> template <typename VectorType> void InverseMatrix<Matrix,Preconditioner>:: vmult (VectorType &dst, const VectorType &src) const { SolverControl solver_control (src.size(), 1e-7*src.l2_norm()); SolverCG<VectorType> cg (solver_control); dst = 0; try { cg.solve (*matrix, dst, src, preconditioner); } catch (std::exception &e) { Assert (false, ExcMessage(e.what())); } }
This is the implementation of the Schur complement preconditioner as described in detail in the introduction. As opposed to step-20 and step-22, we solve the block system all-at-once using GMRES, and use the Schur complement of the block structured matrix to build a good preconditioner instead.
Let's have a look at the ideal preconditioner matrix
described in the introduction. If we apply this matrix in the solution of a linear system, convergence of an iterative GMRES solver will be governed by the matrix
which indeed is very simple. A GMRES solver based on exact matrices would converge in one iteration, since all eigenvalues are equal (any Krylov method takes at most as many iterations as there are distinct eigenvalues). Such a preconditioner for the blocked Stokes system has been proposed by Silvester and Wathen ("Fast iterative solution of stabilised Stokes systems part II. Using general block preconditioners", SIAM J. Numer. Anal., 31 (1994), pp. 1352-1367).
Replacing P by
keeps that spirit alive: the product
will still be close to a matrix with eigenvalues 1 with a distribution that does not depend on the problem size. This lets us hope to be able to get a number of GMRES iterations that is problem-size independent.
The deal.II users who have already gone through the step-20 and step-22 tutorials can certainly imagine how we're going to implement this. We replace the exact inverse matrices in
by some approximate inverses built from the InverseMatrix class, and the inverse Schur complement will be approximated by the pressure mass matrix
(weighted by
as mentioned in the introduction). As pointed out in the results section of step-22, we can replace the exact inverse of A by just the application of a preconditioner, in this case on a vector Laplace matrix as was explained in the introduction. This does increase the number of (outer) GMRES iterations, but is still significantly cheaper than an exact inverse, which would require between 20 and 35 CG iterations for each outer solver step (using the AMG preconditioner).
Having the above explanations in mind, we define a preconditioner class with a vmult functionality, which is all we need for the interaction with the usual solver functions further below in the program code.
First the declarations. These are similar to the definition of the Schur complement in step-20, with the difference that we need some more preconditioners in the constructor and that the matrices we use here are built upon Trilinos:
template <class PreconditionerA, class PreconditionerMp> class BlockSchurPreconditioner : public Subscriptor { public: BlockSchurPreconditioner ( const TrilinosWrappers::BlockSparseMatrix &S, const InverseMatrix<TrilinosWrappers::SparseMatrix, PreconditionerMp> &Mpinv, const PreconditionerA &Apreconditioner); void vmult (TrilinosWrappers::BlockVector &dst, const TrilinosWrappers::BlockVector &src) const; private: const SmartPointer<const TrilinosWrappers::BlockSparseMatrix> stokes_matrix; const SmartPointer<const InverseMatrix<TrilinosWrappers::SparseMatrix, PreconditionerMp > > m_inverse; const PreconditionerA &a_preconditioner; mutable TrilinosWrappers::Vector tmp; }; template <class PreconditionerA, class PreconditionerMp> BlockSchurPreconditioner<PreconditionerA, PreconditionerMp>:: BlockSchurPreconditioner(const TrilinosWrappers::BlockSparseMatrix &S, const InverseMatrix<TrilinosWrappers::SparseMatrix, PreconditionerMp> &Mpinv, const PreconditionerA &Apreconditioner) : stokes_matrix (&S), m_inverse (&Mpinv), a_preconditioner (Apreconditioner), tmp (stokes_matrix->block(1,1).matrix->RowMap()) {}
Next is the vmult function. We implement the action of
as described above in three successive steps. In formulas, we want to compute
where
are both vectors with two block components.
The first step multiplies the velocity part of the vector by a preconditioner of the matrix A, i.e. we compute
. The resulting velocity vector is then multiplied by
and subtracted from the pressure, i.e. we want to compute
. This second step only acts on the pressure vector and is accomplished by the residual function of our matrix classes, except that the sign is wrong. Consequently, we change the sign in the temporary pressure vector and finally multiply by the inverse pressure mass matrix to get the final pressure vector, completing our work on the Stokes preconditioner:
template <class PreconditionerA, class PreconditionerMp> void BlockSchurPreconditioner<PreconditionerA, PreconditionerMp>::vmult ( TrilinosWrappers::BlockVector &dst, const TrilinosWrappers::BlockVector &src) const { a_preconditioner.vmult (dst.block(0), src.block(0)); stokes_matrix->block(1,0).residual(tmp, dst.block(0), src.block(1)); tmp *= -1; m_inverse->vmult (dst.block(1), tmp); } }
BoussinesqFlowProblem class templateThe definition of the class that defines the top-level logic of solving the time-dependent Boussinesq problem is mainly based on the step-22 tutorial program. The main differences are that now we also have to solve for the temperature equation, which forces us to have a second DoFHandler object for the temperature variable as well as matrices, right hand sides, and solution vectors for the current and previous time steps. As mentioned in the introduction, all linear algebra objects are going to use wrappers of the corresponding Trilinos functionality.
The member functions of this class are reminiscent of step-21, where we also used a staggered scheme that first solve the flow equations (here the Stokes equations, in step-21 Darcy flow) and then update the advected quantity (here the temperature, there the saturation). The functions that are new are mainly concerned with determining the time step, as well as the proper size of the artificial viscosity stabilization.
The last three variables indicate whether the various matrices or preconditioners need to be rebuilt the next time the corresponding build functions are called. This allows us to move the corresponding if into the respective function and thereby keeping our main run() function clean and easy to read.
template <int dim> class BoussinesqFlowProblem { public: BoussinesqFlowProblem (); void run (); private: void setup_dofs (); void assemble_stokes_preconditioner (); void build_stokes_preconditioner (); void assemble_stokes_system (); void assemble_temperature_system (); void assemble_temperature_matrix (); double get_maximal_velocity () const; std::pair<double,double> get_extrapolated_temperature_range () const; void solve (); void output_results () const; void refine_mesh (const unsigned int max_grid_level); static double compute_viscosity(const std::vector<double> &old_temperature, const std::vector<double> &old_old_temperature, const std::vector<Tensor<1,dim> > &old_temperature_grads, const std::vector<Tensor<1,dim> > &old_old_temperature_grads, const std::vector<Tensor<2,dim> > &old_temperature_hessians, const std::vector<Tensor<2,dim> > &old_old_temperature_hessians, const std::vector<Vector<double> > &old_stokes_values, const std::vector<Vector<double> > &old_old_stokes_values, const std::vector<double> &gamma_values, const double global_u_infty, const double global_T_variation, const double global_Omega_diameter, const double cell_diameter, const double old_time_step); Triangulation<dim> triangulation; const unsigned int stokes_degree; FESystem<dim> stokes_fe; DoFHandler<dim> stokes_dof_handler; ConstraintMatrix stokes_constraints; std::vector<unsigned int> stokes_block_sizes; TrilinosWrappers::BlockSparseMatrix stokes_matrix; TrilinosWrappers::BlockSparseMatrix stokes_preconditioner_matrix; TrilinosWrappers::BlockVector stokes_solution; TrilinosWrappers::BlockVector old_stokes_solution; TrilinosWrappers::BlockVector stokes_rhs; const unsigned int temperature_degree; FE_Q<dim> temperature_fe; DoFHandler<dim> temperature_dof_handler; ConstraintMatrix temperature_constraints; TrilinosWrappers::SparseMatrix temperature_mass_matrix; TrilinosWrappers::SparseMatrix temperature_stiffness_matrix; TrilinosWrappers::SparseMatrix temperature_matrix; TrilinosWrappers::Vector temperature_solution; TrilinosWrappers::Vector old_temperature_solution; TrilinosWrappers::Vector old_old_temperature_solution; TrilinosWrappers::Vector temperature_rhs; double time_step; double old_time_step; unsigned int timestep_number; boost::shared_ptr<TrilinosWrappers::PreconditionAMG> Amg_preconditioner; boost::shared_ptr<TrilinosWrappers::PreconditionIC> Mp_preconditioner; bool rebuild_stokes_matrix; bool rebuild_temperature_matrices; bool rebuild_stokes_preconditioner; };
The constructor of this class is an extension of the constructor in step-22. We need to add the various variables that concern the temperature. As discussed in the introduction, we are going to use
(Taylor-Hood) elements again for the Stokes part, and
elements for the temperature. However, by using variables that store the polynomial degree of the Stokes and temperature finite elements, it is easy to consistently modify the degree of the elements as well as all quadrature formulas used on them downstream. Moreover, we initialize the time stepping as well as the options for matrix assembly and preconditioning:
template <int dim> BoussinesqFlowProblem<dim>::BoussinesqFlowProblem () : triangulation (Triangulation<dim>::maximum_smoothing), stokes_degree (1), stokes_fe (FE_Q<dim>(stokes_degree+1), dim, FE_Q<dim>(stokes_degree), 1), stokes_dof_handler (triangulation), temperature_degree (2), temperature_fe (temperature_degree), temperature_dof_handler (triangulation), time_step (0), old_time_step (0), timestep_number (0), rebuild_stokes_matrix (true), rebuild_temperature_matrices (true), rebuild_stokes_preconditioner (true) {}
Starting the real functionality of this class is a helper function that determines the maximum (
) velocity in the domain (at the quadrature points, in fact). How it works should be relatively obvious to all who have gotten to this point of the tutorial.
The only point worth thinking about a bit is how to choose the quadrature points we use here. Since the goal of this function is to find the maximal velocity over a domain by looking at quadrature points on each cell. So we should ask how we should best choose these quadrature points on each cell. To this end, recall that if we had a single
field (rather than the vector-valued field of higher order) then the maximum would be attained at a vertex of the mesh. In other words, we should use the QTrapez class that has quadrature points only at the vertices of cells.
For higher order shape functions, the situation is more complicated: the maxima and minima may be attained at points between the support points of shape functions (for the usual
elements the support points are the equidistant Lagrange interpolation points); furthermore, since we are looking for the maximum magnitude of a vector-valued quantity, we can even less say with certainty where the set of potential maximal points are. Nevertheless, intuitively if not provably, the Lagrange interpolation points appear to be a better choice than the Gauss points.
There are now different methods to produce a quadrature formula with quadrature points equal to the interpolation points of the finite element. One option would be to use the FiniteElement::get_unit_support_points() function, reduce the output to a unique set of points to avoid duplicate function evaluations, and create a Quadrature object using these points. Another option, chosen here, is to use the QTrapez class and combine it with the QIterated class that repeats the QTrapez formula on a number of sub-cells in each coordinate direction. To cover all support points, we need to iterate it stokes_degree+1 times since this is the polynomial degree of the Stokes element in use:
template <int dim> double BoussinesqFlowProblem<dim>::get_maximal_velocity () const { const QIterated<dim> quadrature_formula (QTrapez<1>(), stokes_degree+1); const unsigned int n_q_points = quadrature_formula.size(); FEValues<dim> fe_values (stokes_fe, quadrature_formula, update_values); std::vector<Vector<double> > stokes_values(n_q_points, Vector<double>(dim+1)); double max_velocity = 0; typename DoFHandler<dim>::active_cell_iterator cell = stokes_dof_handler.begin_active(), endc = stokes_dof_handler.end(); for (; cell!=endc; ++cell) { fe_values.reinit (cell); fe_values.get_function_values (stokes_solution, stokes_values); for (unsigned int q=0; q<n_q_points; ++q) { Tensor<1,dim> velocity; for (unsigned int i=0; i<dim; ++i) velocity[i] = stokes_values[q](i); max_velocity = std::max (max_velocity, velocity.norm()); } } return max_velocity; }
Next a function that determines the minimum and maximum temperature at quadrature points inside
when extrapolated from the two previous time steps to the current one. We need this information in the computation of the artificial viscosity parameter
as discussed in the introduction.
The formula for the extrapolated temperature is
. The way to compute it is to loop over all quadrature points and update the maximum and minimum value if the current value is bigger/smaller than the previous one. We initialize the variables that store the max and min before the loop over all quadrature points by bounding
, where
is the set of the support points (i.e. nodal points, but note that the maximum of a finite element function can be attained at a point that's not a support point unless one is using
elements). So if we initialize the minimal value by this upper bound, and the maximum value by the negative of this upper bound, then we know for a fact that it is larger/smaller than the minimum/maximum and that the loop over all quadrature points is ultimately going to update the initial value with the correct one.
The only other complication worth mentioning here is that in the first time step,
is not yet available of course. In that case, we can only use
which we have from the initial temperature. As quadrature points, we use the same choice as in the previous function though with the difference that now the number of repetitions is determined by the polynomial degree of the temperature field.
template <int dim> std::pair<double,double> BoussinesqFlowProblem<dim>::get_extrapolated_temperature_range () const { const QIterated<dim> quadrature_formula (QTrapez<1>(), temperature_degree); const unsigned int n_q_points = quadrature_formula.size(); FEValues<dim> fe_values (temperature_fe, quadrature_formula, update_values); std::vector<double> old_temperature_values(n_q_points); std::vector<double> old_old_temperature_values(n_q_points); if (timestep_number != 0) { double min_temperature = (1. + time_step/old_time_step) * old_temperature_solution.linfty_norm() + time_step/old_time_step * old_old_temperature_solution.linfty_norm(), max_temperature = -min_temperature; typename DoFHandler<dim>::active_cell_iterator cell = temperature_dof_handler.begin_active(), endc = temperature_dof_handler.end(); for (; cell!=endc; ++cell) { fe_values.reinit (cell); fe_values.get_function_values (old_temperature_solution, old_temperature_values); fe_values.get_function_values (old_old_temperature_solution, old_old_temperature_values); for (unsigned int q=0; q<n_q_points; ++q) { const double temperature = (1. + time_step/old_time_step) * old_temperature_values[q]- time_step/old_time_step * old_old_temperature_values[q]; min_temperature = std::min (min_temperature, temperature); max_temperature = std::max (max_temperature, temperature); } } return std::make_pair(min_temperature, max_temperature); } else { double min_temperature = old_temperature_solution.linfty_norm(), max_temperature = -min_temperature; typename DoFHandler<dim>::active_cell_iterator cell = temperature_dof_handler.begin_active(), endc = temperature_dof_handler.end(); for (; cell!=endc; ++cell) { fe_values.reinit (cell); fe_values.get_function_values (old_temperature_solution, old_temperature_values); for (unsigned int q=0; q<n_q_points; ++q) { const double temperature = old_temperature_values[q]; min_temperature = std::min (min_temperature, temperature); max_temperature = std::max (max_temperature, temperature); } } return std::make_pair(min_temperature, max_temperature); } }
The last of the tool functions computes the artificial viscosity parameter
on a cell
as a function of the extrapolated temperature, its gradient and Hessian (second derivatives), the velocity, the right hand side
all on the quadrature points of the current cell, and various other parameters as described in detail in the introduction.
There are some universal constants worth mentioning here. First, we need to fix
; we choose
, a choice discussed in detail in the results section of this tutorial program. The second is the exponent
;
appears to work fine for the current program, even though some additional benefit might be expected from chosing
. Finally, there is one thing that requires special casing: In the first time step, the velocity equals zero, and the formula for
is not defined. In that case, we return
, a choice admittedly more motivated by heuristics than anything else (it is in the same order of magnitude, however, as the value returned for most cells on the second time step).
The rest of the function should be mostly obvious based on the material discussed in the introduction:
template <int dim> double BoussinesqFlowProblem<dim>:: compute_viscosity (const std::vector<double> &old_temperature, const std::vector<double> &old_old_temperature, const std::vector<Tensor<1,dim> > &old_temperature_grads, const std::vector<Tensor<1,dim> > &old_old_temperature_grads, const std::vector<Tensor<2,dim> > &old_temperature_hessians, const std::vector<Tensor<2,dim> > &old_old_temperature_hessians, const std::vector<Vector<double> > &old_stokes_values, const std::vector<Vector<double> > &old_old_stokes_values, const std::vector<double> &gamma_values, const double global_u_infty, const double global_T_variation, const double global_Omega_diameter, const double cell_diameter, const double old_time_step) { const double beta = 0.015 * dim; const double alpha = 1; if (global_u_infty == 0) return 5e-3 * cell_diameter; const unsigned int n_q_points = old_temperature.size(); double max_residual = 0; double max_velocity = 0; for (unsigned int q=0; q < n_q_points; ++q) { Tensor<1,dim> u; for (unsigned int d=0; d<dim; ++d) u[d] = (old_stokes_values[q](d) + old_old_stokes_values[q](d)) / 2; const double dT_dt = (old_temperature[q] - old_old_temperature[q]) / old_time_step; const double u_grad_T = u * (old_temperature_grads[q] + old_old_temperature_grads[q]) / 2; const double kappa_Delta_T = EquationData::kappa * (trace(old_temperature_hessians[q]) + trace(old_old_temperature_hessians[q])) / 2; const double residual = std::abs((dT_dt + u_grad_T - kappa_Delta_T - gamma_values[q]) * std::pow((old_temperature[q]+old_old_temperature[q]) / 2, alpha-1.)); max_residual = std::max (residual, max_residual); max_velocity = std::max (std::sqrt (u*u), max_velocity); } const double global_scaling = global_u_infty * global_T_variation / std::pow(global_Omega_diameter, alpha - 2.); return (beta * max_velocity * std::min (cell_diameter, std::pow(cell_diameter,alpha) * max_residual / global_scaling)); }
This is the function that sets up the DoFHandler objects we have here (one for the Stokes part and one for the temperature part) as well as set to the right sizes the various objects required for the linear algebra in this program. Its basic operations are similar to what we do in step-22.
The body of the function first enumerates all degrees of freedom for the Stokes and temperature systems. For the Stokes part, degrees of freedom are then sorted to ensure that velocities precede pressure DoFs so that we can partition the Stokes matrix into a
matrix. As a difference to step-22, we do not perform any additional DoF renumbering. In that program, it paid off since our solver was heavily dependent on ILU's, whereas we use AMG here which is not sensitive to the DoF numbering. The IC preconditioner for the inversion of the pressure mass matrix would of course take advantage of a Cuthill-McKee like renumbering, but its costs are low compared to the velocity portion, so the additional work does not pay off.
We then proceed with the generation of the hanging node constraints that arise from adaptive grid refinement for both DoFHandler objects. For the velocity, we impose no-flux boundary conditions
by adding constraints to the object that already stores the hanging node constraints matrix. The second parameter in the function describes the first of the velocity components in the total dof vector, which is zero here. The variable no_normal_flux_boundaries denotes the boundary indicators for which to set the no flux boundary conditions; here, this is boundary indicator zero.
After having done so, we count the number of degrees of freedom in the various blocks:
template <int dim> void BoussinesqFlowProblem<dim>::setup_dofs () { std::vector<unsigned int> stokes_sub_blocks (dim+1,0); stokes_sub_blocks[dim] = 1; { stokes_dof_handler.distribute_dofs (stokes_fe); DoFRenumbering::component_wise (stokes_dof_handler, stokes_sub_blocks); stokes_constraints.clear (); DoFTools::make_hanging_node_constraints (stokes_dof_handler, stokes_constraints); std::set<unsigned char> no_normal_flux_boundaries; no_normal_flux_boundaries.insert (0); VectorTools::compute_no_normal_flux_constraints (stokes_dof_handler, 0, no_normal_flux_boundaries, stokes_constraints); stokes_constraints.close (); } { temperature_dof_handler.distribute_dofs (temperature_fe); temperature_constraints.clear (); DoFTools::make_hanging_node_constraints (temperature_dof_handler, temperature_constraints); temperature_constraints.close (); } std::vector<unsigned int> stokes_dofs_per_block (2); DoFTools::count_dofs_per_block (stokes_dof_handler, stokes_dofs_per_block, stokes_sub_blocks); const unsigned int n_u = stokes_dofs_per_block[0], n_p = stokes_dofs_per_block[1], n_T = temperature_dof_handler.n_dofs(); std::cout << "Number of active cells: " << triangulation.n_active_cells() << " (on " << triangulation.n_levels() << " levels)" << std::endl << "Number of degrees of freedom: " << n_u + n_p + n_T << " (" << n_u << '+' << n_p << '+'<< n_T <<')' << std::endl << std::endl;
The next step is to create the sparsity pattern for the Stokes and temperature system matrices as well as the preconditioner matrix from which we build the Stokes preconditioner. As in step-22, we choose to create the pattern not as in the first few tutorial programs, but by using the blocked version of CompressedSimpleSparsityPattern. The reason for doing this is mainly memory, that is, the SparsityPattern class would consume too much memory when used in three spatial dimensions as we intend to do for this program.
So, we first release the memory stored in the matrices, then set up an object of type BlockCompressedSimpleSparsityPattern consisting of
blocks (for the Stokes system matrix and precondi