The step-32 tutorial program

Table of contents
  1. Introduction
  2. The commented program
  1. Results
  2. The plain program

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.

Introduction

This program does pretty much exactly what step-31 already does: it solves the Boussinesq equations that describe the motion of a fluid whose temperature is not in equilibrium. As such, all the equations we have described in step-31 still hold: we solve the same partial differential equation, using the same finite element scheme, the same time stepping algorithm, and the same stabilization method for the temperature advection-diffusion equation. As a consequence, you may first want to understand that program before you work on the current one.

The difference between step-31 and the current program is that here we want to do things in parallel, using both the availability of many machines in a cluster (with parallelization based on MPI) as well as many processor cores within a single machine (with parallelization based on threads). This program's main job is therefore to introduce the changes that are necessary to utilize the availability of these parallel compute resources.

In addition to these changes, we also use a slightly different preconditioner, and we will have to make a number of changes that have to do with the fact that we want to solve a realistic problem here, not a model problem. The latter, in particular, will require that we think about scaling issues as well as what all those parameters and coefficients in the equations under consideration actually mean. We will discuss first the issues that affect changes in the mathematical formulation and solver structure, then how to parallelize things, and finally the actual testcase we will consider.

Using the "right" pressure

In step-31, we used the following Stokes model for the velocity and pressure field:

\begin{eqnarray*} -\nabla \cdot (2 \eta \varepsilon ({\mathbf u})) + \nabla p &=& -\rho \; \beta \; T \mathbf{g}, \\ \nabla \cdot {\mathbf u} &=& 0. \end{eqnarray*}

The right hand side of the first equation appears a wee bit unmotivated. Here's how things should really be. On the right side, we need the external forces that act on the fluid, which we assume are given by gravity only. In the current case, we assume that the fluid does expand slightly for the purposes of this gravity force, but not enough that we need to modify the incompressibility condition (the second equation). What this means is that for the purpose of the right hand side, we can assume that $\rho=\rho(T)$. An assumption that may not be entirely justified is that we can assume that the changes of density as a function of temperature are small, leading to an expression of the form $\rho(T) = \rho_{\text{ref}} [1-\beta(T-T_{\text{ref}})]$, i.e. the density equals $\rho_{\text{ref}}$ at reference temperature and decreases linearly as the temperature increases (as the material expands). The force balance equation then looks properly written like this:

\begin{eqnarray*} -\nabla \cdot (2 \eta \varepsilon ({\mathbf u})) + \nabla p &=& \rho_{\text{ref}} [1-\beta(T-T_{\text{ref}})] \mathbf{g}. \end{eqnarray*}

Now note that the gravity force results from a gravity potential as $\mathbf g=-\nabla \varphi$, so that we can re-write this as follows:

\begin{eqnarray*} -\nabla \cdot (2 \eta \varepsilon ({\mathbf u})) + \nabla p &=& -\rho_{\text{ref}} \; \beta\; T\; \mathbf{g} -\rho_{\text{ref}} [1+\beta T_{\text{ref}}] \nabla\varphi. \end{eqnarray*}

The second term on the right is time independent, and so we could introduce a new "dynamic" pressure $p_{\text{dyn}}=p+\rho_{\text{ref}} [1+\beta T_{\text{ref}}] \varphi=p_{\text{total}}-p_{\text{static}}$ with which the Stokes equations would read:

\begin{eqnarray*} -\nabla \cdot (2 \eta \varepsilon ({\mathbf u})) + \nabla p_{\text{dyn}} &=& -\rho_{\text{ref}} \; \beta \; T \; \mathbf{g}, \\ \nabla \cdot {\mathbf u} &=& 0. \end{eqnarray*}

This is exactly the form we used in step-31, and it was appropriate to do so because all changes in the fluid flow are only driven by the dynamic pressure that results from temperature differences.

On the other hand, we will here use the form of the Stokes equations that considers the total pressure instead:

\begin{eqnarray*} -\nabla \cdot (2 \eta \varepsilon ({\mathbf u})) + \nabla p &=& \rho(T)\; \mathbf{g}, \\ \nabla \cdot {\mathbf u} &=& 0. \end{eqnarray*}

This way we can plot the pressure in our program in such a way that it actually shows the total pressure that includes the effects of temperature differences as well as the static pressure of the overlying rocks. Since the pressure does not appear any further in any of the other equations, whether to use one or the other is more a matter of taste than of correctness. The flow field is exactly the same, but we get a pressure that we can now compare with values that are given in geophysical books as those that hold at the bottom of the earth mantle, for example.

A second reason to do this is discussed in the results section and concerns possible extensions to the model we use here. It has to do with the fact that while the temperature equation we use here does not include a term that contains the pressure. It should, however: rock, like gas, heats up as you compress it. Consequently, material that rises up cools adiabatically, and cold material that sinks down heats adiabatically. We discuss this further below.

The scaling of discretized equations

Remember that we want to solve the following set of equations:

\begin{eqnarray*} -\nabla \cdot (2 \eta \varepsilon ({\mathbf u})) + \nabla p &=& \rho(T) \mathbf{g}, \\ \nabla \cdot {\mathbf u} &=& 0, \\ \frac{\partial T}{\partial t} + {\mathbf u} \cdot \nabla T - \nabla \cdot \kappa \nabla T &=& \gamma, \end{eqnarray*}

augmented by appropriate boundary and initial conditions. As discussed in step-31, we will solve this set of equations by solving for a Stokes problem first in each time step, and then moving the temperature equation forward by one time interval.

The problem under consideration in this current section is with the Stokes problem: if we discretize it as usual, we get a linear system

\begin{eqnarray*} M \; X = \left(\begin{array}{cc} A & B^T \\ B & 0 \end{array}\right) \left(\begin{array}{c} U \\ P \end{array}\right) = \left(\begin{array}{c} F_U \\ 0 \end{array}\right) = F \end{eqnarray*}

which in this program we will solve with a BiCGStab solver. BiCGStab iterates until the residual of these linear equations is below a certain tolerance, i.e. until

\[ \left\| \left(\begin{array}{c} F_U - A U^{(k)} - B P^{(k)} \\ B^T U^{(k)} \end{array}\right) \right\| < \text{Tol}. \]

This does not make any sense from the viewpoint of physical units: the quantities involved here have physical units so that the first part of the residual has units $\frac{\text{Pa}}{\text{m}} \text{m}^{\text{dim}}$ (most easily established by considering the term $(\nabla \cdot \mathbf v, p)_{\Omega}$ and considering that the pressure has units $\text{Pa}=\frac{\text{kg}}{\text{m\; s}^2}$ and the integration yields a factor of $\text{m}^{\text{dim}}$), whereas the second part of the residual has units $\frac{\text{m}^{\text{dim}}}{\text{s}}$. Taking the norm of this residual vector would yield a quantity with units $\text{m}^{\text{dim}-1} \sqrt{\left(\text{Pa}\right)^2 + \left(\frac{\text{m}}{\text{s}}\right)^2}$. This, quite obviously, does not make sense, and we should not be surprised that doing so is eventually going to come back hurting us.

So why is this an issue here, but not in step-31? The reason back there is that everything was nicely balanced: velocities were on the order of one, the pressure likewise, the viscosity was one, and the domain had a diameter of $\sqrt{2}$. As a result, while non-sensical, nothing bad happened. On the other hand, as we will explain below, things here will not be that simply scaled: $\eta$ will be around $10^{21}$, velocities on the order of $10^{-8}$, pressure around $10^8$, and the diameter of the domain is $10^7$. In other words, the order of magnitude for the first equation is going to be $\eta\text{div}\varepsilon(\mathbf u) \approx 10^{21} \frac{10^{-8}}{(10^7)^2} \approx 10^{-1}$, whereas the second equation will be around $\text{div}{\mathbf u}\approx \frac{10^{-8}}{10^7} \approx 10^{-15}$. Well, so what this will lead to is this: if BiCGStab wants to make the residual small, it will almost entirely focus on the first set of equations because they are so much bigger, and ignore the divergence equation that describes mass conservation. That's exactly what happens: unless we set the tolerance to extremely small values, the resulting flow field is definitely not divergence free. As an auxiliary problem, it turns out that it is difficult to find a tolerance that always works; in practice, one often ends up with a tolerance that requires 30 or 40 iterations for most time steps, and 10,000 for some others.

So what's a numerical analyst to do in a case like this? The answer is to start at the root and first make sure that everything is mathematically consistent first. In our case, this means that if we want to solve the system of Stokes equations jointly, we have to scale them so that they all have the same physical dimensions. In our case, this means multiplying the second equation by something that has units $\frac{\text{Pa\; s}}{\text{m}}$; one choice is to multiply with $\frac{\eta}{L}$ where $L$ is a typical lengthscale in our domain. Using the numbers above, this factor is around $10^{14}$, which just so happens to be the order of magnitude that would make the two equations numerically about the same. So, we now get this for the Stokes system:

\begin{eqnarray*} -\nabla \cdot (2 \eta \varepsilon ({\mathbf u})) + \nabla p &=& \rho(T) \; \mathbf{g}, \\ \frac{\eta}{L} \nabla \cdot {\mathbf u} &=& 0. \end{eqnarray*}

The trouble with this is that the result is not symmetric any more (we have $\frac{\eta}{L} \nabla \cdot$ at the bottom left, but not its transpose operator at the top right). This, however, can be cured by introducing a scaled pressure $\hat p = \frac{L}{\eta}p$, and we get the scaled equations

\begin{eqnarray*} -\nabla \cdot (2 \eta \varepsilon ({\mathbf u})) + \nabla \left(\frac{\eta}{L} \hat p\right) &=& \rho(T) \; \mathbf{g}, \\ \frac{\eta}{L} \nabla \cdot {\mathbf u} &=& 0. \end{eqnarray*}

This is now symmetric. Obviously, we can easily recover the original pressure $p$ from the scaled pressure $\hat p$ that we compute as a result of this procedure.

In the program below, we will introduce a factor EquationData::pressure_scaling that corresponds to $\frac{\eta}{L}$, and we will use this factor in the assembly of the system matrix and preconditioner. We will recover the unscaled pressure in the output_results function.

Changes to the Stokes preconditioner

In this tutorial program, we apply a variant of the preconditioner used in step-31. That preconditioner was built to operate on the system matrix M in block form such that the product matrix

\begin{eqnarray*} P^{-1} M = \left(\begin{array}{cc} A^{-1} & 0 \\ S^{-1} B A^{-1} & -S^{-1} \end{array}\right) \left(\begin{array}{cc} A & B^T \\ B & 0 \end{array}\right) \end{eqnarray*}

is of a form that Krylov-based iterative solvers like GMRES can solve in a few iterations. We then replaced the exact inverse of A by the action of an AMG preconditioner $\tilde{A}$ based on a vector Laplace matrix, approximated the Schur complement $S = B A^{-1} B^T$ by a mass matrix $M_p$ on the pressure space and wrote an InverseMatrix class for implementing the action of $M_p^{-1}\approx S^{-1}$ on vectors. In the InverseMatrix class, we used a CG solve with an incomplete Cholesky (IC) preconditioner for performing the inner solves.

An observation one can make is that we use just the action of a preconditioner for approximating the velocity inverse $A^{-1}$ (and the outer GMRES iteration takes care of the approximate character of the inverse), whereas we use a more or less exact inverse for $M_p^{-1}$, realized by a fully converged CG solve. What we change here is to skip that exact inverse matrix and replace it – as usual – by the action of a preconditioner only. This works, as we will demonstrate below. For efficiency reasons, we want to avoid increasing the number of iterations for the block solve. Keep in mind that most of the time in the solution of the matrix system is the application of the AMG preconditioner (about half the time of the total solve), and the application of matrix A (about one third of the total solve time). This means that we really do not want to do those operations more often when we remove the inner solve on the Schur complement approximation. It turns out that the Trilinos IC preconditioner would not fulfill this requirement, however, the Trilinos ILU preconditioner does. It does even better than so — it decreases the iteration count for large 3D problems. The reason for that decrease is that we avoid some errors that CG introduces: Even a converged solve has still some residual. That is a problem because that small error interferes with the outer iterative solver, probably because a CG solver does some nonlinear operations by weighting vectors by some inner products.

Except the simplification in the preconditioner, we replaced the GMRES solver by BiCGStab. This is merely to demonstrate that GMRES is not the only possible option for a saddle point system like the Stokes equations. BiCGStab harmonizes nicely with the ILU preconditioner on a pressure mass matrix as approximation for $S^{-1}$, so it is at least as good as GMRES in this example. Keep in mind the discussion in the results section of the step-22 tutorial program, where we observed that BiCGStab does not like inner solves with CG, which made us prefer GMRES in step-31.

As a final note, let us remark that in step-31 we computed the Schur complement $S=B A^{-1} B^T$ by approximating $-\text{div}(-\eta\Delta)^{-1}\nabla \approx \frac 1{\eta} \mathbf{1}$. Now, however, we have re-scaled the $B$ and $B^T$ operators. So $S$ should now approximate $-\frac{\eta}{L}\text{div}(-\eta\Delta)^{-1}\nabla \frac{\eta}{L} \approx \left(\frac{\eta}{L}\right)^2 \frac 1{\eta} \mathbf{1}$. This is exactly the operator we use to approximate $S$.

Changes to the artificial viscosity stabilization

As in step-31, we will use an artificial viscosity of the form

\begin{eqnarray*} \nu_\alpha(T)|_K = \beta \|\mathbf{u}\|_{L^\infty(K)} \min\left\{ h_K, h_K^\alpha \frac{\|R_\alpha(T)\|_{L^\infty(K)}}{c(\mathbf{u},T)} \right\} \end{eqnarray*}

in this problem, where $c(\mathbf{u},T) = c_R\ \|\mathbf{u}\|_{L^\infty(\Omega)} \ \mathrm{var}(T) \ |\mathrm{diam}(\Omega)|^{\alpha-2}$ (for the meaning of the various terms in these formulas, see step-31. In the results section of that program, we have discussed our choice for $c_R$ and how we arrived at the value used there mostly by accident, and in more detail how $\beta$ was chosen. For the current program, we want to go about this issue a bit more systematically for both parameters, using the same line of reasoning with which we chose two other parameters in our discretization, $c_k$ and $\beta$, in the results section of step-31. In particular, remember that we would like to make the artificial viscosity as small as possible while keeping it as large as necessary. To see what is happening, note that below we will impose boundary conditions for the temperature between 973 and 4273 Kelvin, and initial conditions are also chosen in this range; because there are no internal heat sources or sinks, the temperature should consequently always be in this range, barring any internal oscillations. If the minimal temperature drops below 973 Kelvin, then we need to add stabilization by either increasing $\beta$ or decreasing $c_R$.

As we did in step-31, we first determine an optimal value of $\beta$ by using the "traditional" formula

\begin{eqnarray*} \nu_\alpha(T)|_K = \beta \|\mathbf{u}\|_{L^\infty(K)} h_K, \end{eqnarray*}

which we know to be stable if only $\beta$ is large enough. Doing a couple hundred time steps (on a coarser mesh than the one shown in the program, and with a different viscosity that affects transport velocities and therefore time step sizes) in 2d will produce the following graph:

step-32.beta.2d.png

As can be seen, values $\beta \le 0.05$ are too small whereas $\beta=0.052$ appears to work, at least to the time horizon shown here. As a remark on the side, there are at least two questions one may wonder here: First, what happens at the time when the solution becomes unstable? Looking at the graphical output, we can see that with the unreasonably coarse mesh chosen for these experiments, around time $t=10^{15}$ seconds the plumes of hot material that have been rising towards the cold outer boundary and have then spread sideways are starting to get close to each other, squezzing out the cold material inbetween. This creates a layer of cells into which fluids flows from two opposite sides and flows out toward a third, apparently a scenario that then produce these instabilities without sufficient stabilization. Second: In step-31, we used $\beta=0.015\cdot\text{dim}$; why does this not work here? The answer to this is not entirely clear -- stabilization parameters are certainly known to depend on things like the shape of cells, for which we had square in step-31 but have trapezoids in the current program. Whatever the exact cause, we at least have a value of $\beta$, namely 0.052 for 2d, that works for the current program.

With this value fixed, we can go back to the original formula for the viscosity $\nu$ and play with the constant $c_R$, making it as large as possible in order to make $\nu$ as small as possible. This gives us a picture like this:

step-32.beta_cr.2d.png

Consequently, $c_R=0.1$ would appear to be the right value.

Parallelization on clusters

Parallelization of scientific codes across multiple machines in a cluster of computers is almost always done using the Message Passing Interface (MPI). This program is no exception to that, and it follows the step-17 and step-18 programs in this.

MPI is a rather awkward interface to program with, and so we usually try to not use it directly but through an interface layer that abstracts most of the MPI operations into a friendlier interface. In the two programs mentioned above, this was achieved by using the PETSc library that provides support for parallel linear algebra in a way that almost completely hides the MPI layer under it. PETSc is powerful, providing a large number of functions that deal with matrices, vectors, and iterative solvers and preconditioners, along with lots of other stuff, most of which runs quite well in parallel. It is, however, a few years old already, written in C, and generally not quite as easy to use as some other libraries. As a consequence, deal.II also has interfaces to Trilinos, a library similar to PETSc in its aims and with a lot of the same functionality. It is, however, a project that is several years younger, is written in C++ and by people who generally have put a significant emphasis on software design. We have already used Trilinos in step-31, and will do so again here, with the difference that we will use its parallel capabilities.

deal.II's Trilinos interfaces encapsulate pretty much everything Trilinos provides into wrapper classes (in namespace TrilinosWrappers) that make the Trilinos matrix, vector, solver and preconditioner classes look very much the same as deal.II's own implementations of this functionality. However, as opposed to deal.II's classes, they can be used in parallel if we give them the necessary information. As a consequence, there are two Trilinos classes that we have to deal with directly (rather than through wrappers), both of which are part of Trilinos' Epetra library of basic linear algebra and tool classes:

The only other things specific to programming using MPI that we will use in this program are the following facilities deal.II provides:

The rest of the program is almost completely agnostic about the fact that we don't store all objects completely locally. There will be a few points where we can not use certain programming techniques (though without making explicit reference to MPI or parallelization) or where we need access to all elements of a vector and therefore need to localize its elements (i.e. create a vector that has all its elements stored on the current machine), but we will comment on these locations as we get to them in the program code.

Parallelization within individual nodes of a cluster

The second strategy to parallelize a program is to make use of the fact that most computers today have more than one processor that all have access to the same memory. In other words, in this model, we don't explicitly have to say which pieces of data reside where -- all of the data we need is directly accessible and all we have to do is split processing this data between the available processors. We will then couple this with the MPI parallelization outlined above, i.e. we will have all the processors on a machine work together to, for example, assemble the local contributions to the global matrix for the cells that this machine actually "owns" but not for those cells that are owned by other machines. We will use this strategy for four kinds of operations we frequently do in this program: assembly of the Stokes and temperature matrices, assembly of the matrix that forms the Stokes preconditioner, and assembly of the right hand side of the temperature system.

All of these operations essentially look as follows: we need to loop over all cells for which cell->subdomain_id() equals the index our machine has within the communicator object used for all communication (i.e. essentially MPI_COMM_WORLD, as explained above), on each cell we need to assemble the local contributions to the global matrix or vector, and then we have to copy each cell's contribution into the global matrix or vector. Note that the first part of this (the loop) defines a range of iterators on which something has to happen. The second part, assembly of local contributions is something that takes the majority of CPU time in this sequence of steps, is a typical example of things that can be done in parallel: each cell's contribution is entirely independent of all other cells' contributions. The third part, copying into the global matrix, must not happen in parallel since we are modifying one object and so several threads can not at the same time read an existing matrix element, add their contribution, and write the sum back into memory without danger of producing a race condition.

deal.II has a class that is made for exactly this workflow: WorkStream. Its use is extensively documented in the module on Parallel computing with multiple processors accessing shared memory (in the section on the WorkStream class) and we won't repeat here the rationale and detailed instructions laid out there, though you will want to read through this module to understand the distinction between scratch space and per-cell data. Suffice it to mention that we need the following:

We will comment on a few more points in the actual code, but in general their structure should be clear from the discussion in Parallel computing with multiple processors accessing shared memory.

The underlying technology for WorkStream identifies "tasks" that need to be worked on (e.g. assembling local contributions on a cell) and schedules these tasks automatically to available processors. WorkStream creates these tasks automatically, by splitting the iterator range into suitable chunks, but as outlined in Parallel computing with multiple processors accessing shared memory, one can also create tasks explicitly. We use this in one place in the program, namely where we set up the Stokes system and preconditioner matrices as well as the temperature matrix. These are independent operations that, if enough processors are available, can be worked on in parallel (if not enough processors are available -- because the system has only one, or because the others are working on something else for us -- then these tasks will be worked on sequentially). Consequently, BoussinesqFlowProblem::setup_dofs creates tasks for the three calls to BoussinesqFlowProblem::setup_stokes_matrix, BoussinesqFlowProblem::setup_stokes_preconditioner, and BoussinesqFlowProblem::setup_temperature_matrices that are then scheduled to available resources. There is one problem with this, however: if we have more than one MPI process running in parallel, then all of these processes need to communicate in a certain order and that requires that the various setup_* functions can't run in parallel. To make things worse, even if there is only a single MPI process, we still have to make a few calls to the MPI runtime (for example to set up communicators for the various matrices and vectors; all of these communicators contain only a single processor, but we still have to call MPI functions for this). Unfortunately, most MPI implementations are not thread-safe, and we can't call MPI functions from multiple threads at once. For these setup operations, we will therefore only be able to make use of multiple processor cores within the same machine if we are not running under the control of MPI -- a condition we can check using the Utilities::System::job_supports_mpi() function.

Implementation details

One detail worth discussing before we show the actual implementation is how to deal with distributed vs. localized vectors. When we build the linear systems, the various matrices and right hand side vectors are always built distributed, i.e. each processor builds contributions for a subset of cells and adds them to the global matrix objects; Trilinos then makes sure that the right entries end up on the right processors so that each processor stores a part of these distributed objects.

When we then come around to solving linear systems, the solution vectors will also be distributed. On the other hand, in later steps, processors will need to be able to access random elements for various tasks. For example, processor zero needs to have access to all elements in order to produce output; all processors need to have access to temperature and Stokes solution data on the cells they own to build right hand side vectors for the next time step; and all processors need to access data on the cells they own as well as their neighbors (whether these neighbors are owned by the same process or another one) in order to compute the jump of the gradient across cell interfaces in the Kelly error estimator when computing refinement information.

In other words, for some operations we will have to exchange information between processors. We could try to be smart and really only exchange that data that we really need. This would probably be important if we were to run this program on thousands of processors, as then storing all elements of a solution vector may become impractical. On the other hand, deal.II currently has a number of other bottlenecks that make sure that thousands of processors is not possible anyway. Consequently, here we opt for the simplest choice: solve linear systems using distributed vectors, and immediately afterwards each processor requests a localized copy of the solution vector containing all elements that we will use henceforth for all operations. The distributed vector is no longer necessary at this point and is, in fact, deallocated again immediately.

The testcase

The setup for this program is mildly reminiscent of the problem we wanted to solve in the first place (see the introduction of step-31): convection in the earth mantle. As a consequence, we choose the following data, all of which appears in the program in units of meters and seconds (the SI system) even if we list them here in other units.

As a reminder, let us again state the equations we want to solve are these:

\begin{eqnarray*} -\nabla \cdot (2 \eta \varepsilon ({\mathbf u})) + \nabla \left( \frac{\eta}{L} \hat p\right) &=& \rho(T) \mathbf{g}, \\ \frac{\eta}{L} \nabla \cdot {\mathbf u} &=& 0, \\ \frac{\partial T}{\partial t} + {\mathbf u} \cdot \nabla T - \nabla \cdot \kappa \nabla T &=& \gamma, \end{eqnarray*}

augmented by boundary and initial conditions. We then have to choose data for the following quantities:

All of these pieces of equation data are defined in the program in the EquationData namespace. When run, the program produces long-term maximal velocities around 10-40 centimeters per year (see the results section below), approximately the physically correct order of magnitude. We will set the end time to 1 billion years.

The commented program

Include files

We 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/function.h>
 #include <base/utilities.h>
 #include <base/conditional_ostream.h>
 #include <base/work_stream.h>
 #include <base/timer.h>
 
 #include <lac/full_matrix.h>
 #include <lac/solver_bicgstab.h>
 #include <lac/solver_cg.h>
 #include <lac/constraint_matrix.h>
 #include <lac/block_sparsity_pattern.h>
 #include <lac/trilinos_block_vector.h>
 #include <lac/trilinos_sparse_matrix.h>
 #include <lac/trilinos_block_sparse_matrix.h>
 #include <lac/trilinos_precondition.h>
 
 #include <grid/tria.h>
 #include <grid/grid_generator.h>
 #include <grid/tria_accessor.h>
 #include <grid/tria_iterator.h>
 #include <grid/filtered_iterator.h>
 #include <grid/tria_boundary_lib.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 <fe/fe_q.h>
 #include <fe/fe_dgq.h>
 #include <fe/fe_system.h>
 #include <fe/fe_values.h>
 #include <fe/mapping_q1.h>
 
 #include <numerics/vectors.h>
 #include <numerics/matrices.h>
 #include <numerics/data_out.h>
 #include <numerics/error_estimator.h>
 #include <numerics/solution_transfer.h>
 
 #include <fstream>
 #include <iostream>
 #include <sstream>

This is the only include file that is new: We use an IndexSet to describe the parallel partitioning of vectors and matrices.

 #include <base/index_set.h>

Next, we import all deal.II names into global namespace:

 using namespace dealii;

Equation data

In the following namespace, we define the various pieces of equation data. All of these are exhaustively discussed in the description of the testcase in the introduction:

 namespace EquationData
 {
   const double eta                   = 1e21;
   const double kappa                 = 1e-6;
   const double reference_density     = 3300;
   const double reference_temperature = 293;
   const double expansion_coefficient = 2e-5;
 
   const double R0      = 6371000.-2890000.;
   const double R1      = 6371000.-  35000.;
 
   const double T0      = 4000+273;
   const double T1      =  700+273;
 
   const double year_in_seconds  = 60*60*24*365.2425;
   const double end_time         = 1e9 * year_in_seconds;
 
   const double pressure_scaling = eta / (R1-R0);
 
 
   double density (const double temperature)
   {
     return (reference_density *
             (1 - expansion_coefficient * (temperature -
                                           reference_temperature)));
   }
 
 
   template <int dim>
   Tensor<1,dim> gravity_vector (const Point<dim> &p)
   {
     return -9.81 * p / R1;
   }
 
 
   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>  &p,
                                         const unsigned int) const
   {
     const double r = p.norm();
     const double h = R1-R0;
 
     const double s = (r-R0)/h;
 
     return T1+(T0-T1)*((1-s)*(1-s));
   }
 
 
   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
   {
     return 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);
   }
 }

Linear solvers and preconditioners

In comparison to step-31, we did one change in the linear algebra of the problem: We exchange the InverseMatrix that previously held the approximation of the Schur complement by a preconditioner only (we will choose ILU in the application code below), as discussed in the introduction. This trick we already did for the velocity block - the idea of this is that the solver iterations on the block system will eventually also make the approximation for the Schur complement good. If the preconditioner we're using is good enough, there will be no increase in the outer iteration count compared to using converged solves for the inverse matrices of velocity and Schur complement. All we need to do for implementing that change is to give the respective variable in the BlockSchurPreconditioner class another name.

 namespace LinearSolvers
 {
   template <class PreconditionerA, class PreconditionerMp>
   class BlockSchurPreconditioner : public Subscriptor
   {
     public:
       BlockSchurPreconditioner (
         const TrilinosWrappers::BlockSparseMatrix  &S,
         const PreconditionerMp                     &Mppreconditioner,
         const PreconditionerA                      &Apreconditioner);
 
       void vmult (TrilinosWrappers::MPI::BlockVector       &dst,
                   const TrilinosWrappers::MPI::BlockVector &src) const;
 
     private:
       const SmartPointer<const TrilinosWrappers::BlockSparseMatrix> stokes_matrix;
       const PreconditionerMp &mp_preconditioner;
       const PreconditionerA  &a_preconditioner;
       mutable TrilinosWrappers::MPI::Vector tmp;
   };
 
 
 
   template <class PreconditionerA, class PreconditionerMp>
   BlockSchurPreconditioner<PreconditionerA, PreconditionerMp>::
   BlockSchurPreconditioner(const TrilinosWrappers::BlockSparseMatrix &S,
                            const PreconditionerMp                    &Mppreconditioner,
                            const PreconditionerA                     &Apreconditioner)
                   :
                   stokes_matrix     (&S),
                   mp_preconditioner (Mppreconditioner),
                   a_preconditioner  (Apreconditioner),
                   tmp               (stokes_matrix->block(1,1).row_partitioner())
   {}
 
 
 
   template <class PreconditionerA, class PreconditionerMp>
   void BlockSchurPreconditioner<PreconditionerA, PreconditionerMp>::vmult (
     TrilinosWrappers::MPI::BlockVector       &dst,
     const TrilinosWrappers::MPI::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;
     mp_preconditioner.vmult (dst.block(1), tmp);
   }
 }

Definition of assembly data structures

As described in the introduction, we will use the WorkStream mechanism discussed in the Parallel computing with multiple processors accessing shared memory module to parallelize operations among the processors of a single machine. The WorkStream class requires that data is passed around in two kinds of data structures, one for scratch data and one to pass data from the assembly function to the function that copies local contributions into global objects.

The following namespace (and the two sub-namespaces) contains a collection of data structures that serve this purpose, one pair for each of the four operations discussed in the introduction that we will want to parallelize. Each assembly routine gets two sets of data: a Scratch array that collects all the classes and arrays that are used for the calculation of the cell contribution, and a CopyData array that keeps local matrices and vectors which will be written into the global matrix. Whereas CopyData is a container for the final data that is written into the global matrices and vector (and, thus, absolutely necessary), the Scratch arrays are merely there for performance reasons — it would be much more expensive to set up a FEValues object on each cell, than creating it only once and updating some derivative data.

Using the program in step-31, we have four assembly routines. One for the preconditioner matrix of the Stokes system, one for the Stokes matrix and right hand side, one for the temperature matrices and one for the right hand side of the temperature equation. We organize the scratch arrays and a CopyData arrays for each of those four assembly components using a struct environment.

Regarding the Scratch array, each struct is equipped with a constructor that create an FEValues object for a finite element, a quadrature formula and some update flags. Moreover, we manually implement a copy constructor (since the FEValues class is not copyable by itself), and provide some additional vector fields that are used to improve performance of assembly.

 namespace Assembly
 {
   namespace Scratch
   {
     template <int dim>
     struct StokesPreconditioner
     {
         StokesPreconditioner (const FiniteElement<dim> &stokes_fe,
                               const Quadrature<dim>    &stokes_quadrature,
                               const UpdateFlags         update_flags);
         StokesPreconditioner (const StokesPreconditioner &data);
 
         FEValues<dim>               stokes_fe_values;
 
         std::vector<Tensor<2,dim> > grad_phi_u;
         std::vector<double>         phi_p;
     };
 
     template <int dim>
     StokesPreconditioner<dim>::
     StokesPreconditioner (const FiniteElement<dim> &stokes_fe,
                           const Quadrature<dim>    &stokes_quadrature,
                           const UpdateFlags         update_flags)
                     :
                     stokes_fe_values (stokes_fe, stokes_quadrature,
                                       update_flags),
                     grad_phi_u (stokes_fe.dofs_per_cell),
                     phi_p (stokes_fe.dofs_per_cell)
     {}
 
 
 
     template <int dim>
     StokesPreconditioner<dim>::
     StokesPreconditioner (const StokesPreconditioner &scratch)
                     :
                     stokes_fe_values (scratch.stokes_fe_values.get_fe(),
                                       scratch.stokes_fe_values.get_quadrature(),
                                       scratch.stokes_fe_values.get_update_flags()),
                     grad_phi_u (scratch.grad_phi_u),
                     phi_p (scratch.phi_p)
     {}

Observe that we derive the StokesSystem scratch array from the StokesPreconditioner array. We do this because all the objects that are necessary for the assembly of the preconditioner are also needed for the actual matrix system and right hand side, plus some extra data. This makes the program more compact. Note also that the assembly of the Stokes system and the temperature right hand side further down requires data from temperature and velocity, respectively, so we actually need two FEValues objects for those two cases.

     template <int dim>
     struct StokesSystem : public StokesPreconditioner<dim>
     {
         StokesSystem (const FiniteElement<dim> &stokes_fe,
                       const Quadrature<dim>    &stokes_quadrature,
                       const UpdateFlags         stokes_update_flags,
                       const FiniteElement<dim> &temperature_fe,
                       const UpdateFlags         temperature_update_flags);
 
         StokesSystem (const StokesSystem<dim> &data);
 
         FEValues<dim>  temperature_fe_values;
 
         std::vector<Tensor<1,dim> >          phi_u;
         std::vector<SymmetricTensor<2,dim> > grads_phi_u;
         std::vector<double>                  div_phi_u;
 
         std::vector<double>                  old_temperature_values;
     };
 
 
     template <int dim>
     StokesSystem<dim>::
     StokesSystem (const FiniteElement<dim> &stokes_fe,
                   const Quadrature<dim>    &stokes_quadrature,
                   const UpdateFlags         stokes_update_flags,
                   const FiniteElement<dim> &temperature_fe,
                   const UpdateFlags         temperature_update_flags)
                     :
                     StokesPreconditioner<dim> (stokes_fe, stokes_quadrature,
                                                stokes_update_flags),
                     temperature_fe_values (temperature_fe, stokes_quadrature,
                                            temperature_update_flags),
                     phi_u (stokes_fe.dofs_per_cell),
                     grads_phi_u (stokes_fe.dofs_per_cell),
                     div_phi_u (stokes_fe.dofs_per_cell),
                     old_temperature_values (stokes_quadrature.n_quadrature_points)
     {}
 
 
     template <int dim>
     StokesSystem<dim>::
     StokesSystem (const StokesSystem<dim> &scratch)
                     :
                     StokesPreconditioner<dim> (scratch),
                     temperature_fe_values (scratch.temperature_fe_values.get_fe(),
                                            scratch.temperature_fe_values.get_quadrature(),
                                            scratch.temperature_fe_values.get_update_flags()),
                     phi_u (scratch.phi_u),
                     grads_phi_u (scratch.grads_phi_u),
                     div_phi_u (scratch.div_phi_u),
                     old_temperature_values (scratch.old_temperature_values)
     {}
 
 
 
     template <int dim>
     struct TemperatureMatrix
     {
         TemperatureMatrix (const FiniteElement<dim> &temperature_fe,
                            const Quadrature<dim>    &temperature_quadrature);
         TemperatureMatrix (const TemperatureMatrix &data);
 
         FEValues<dim>               temperature_fe_values;
 
         std::vector<double>         phi_T;
         std::vector<Tensor<1,dim> > grad_phi_T;
     };
 
     template <int dim>
     TemperatureMatrix<dim>::
     TemperatureMatrix (const FiniteElement<dim> &temperature_fe,
                        const Quadrature<dim>    &temperature_quadrature)
                     :
                     temperature_fe_values (temperature_fe, temperature_quadrature,
                                            update_values    | update_gradients |
                                            update_JxW_values),
                     phi_T (temperature_fe.dofs_per_cell),
                     grad_phi_T (temperature_fe.dofs_per_cell)
     {}
 
 
     template <int dim>
     TemperatureMatrix<dim>::
     TemperatureMatrix (const TemperatureMatrix &scratch)
                     :
                     temperature_fe_values (scratch.temperature_fe_values.get_fe(),
                                            scratch.temperature_fe_values.get_quadrature(),
                                            scratch.temperature_fe_values.get_update_flags()),
                     phi_T (scratch.phi_T),
                     grad_phi_T (scratch.grad_phi_T)
     {}
 
 
     template <int dim>
     struct TemperatureRHS
     {
         TemperatureRHS (const FiniteElement<dim> &temperature_fe,
                         const FiniteElement<dim> &stokes_fe,
                         const Quadrature<dim>    &quadrature);
         TemperatureRHS (const TemperatureRHS &data);
 
         FEValues<dim>               temperature_fe_values;
         FEValues<dim>               stokes_fe_values;
 
         std::vector<double>         phi_T;
         std::vector<Tensor<1,dim> > grad_phi_T;
 
         std::vector<Tensor<1,dim> > old_velocity_values;
         std::vector<Tensor<1,dim> > old_old_velocity_values;
 
         std::vector<double>         old_temperature_values;
         std::vector<double>         old_old_temperature_values;
         std::vector<Tensor<1,dim> > old_temperature_grads;
         std::vector<Tensor<1,dim> > old_old_temperature_grads;
         std::vector<double>         old_temperature_laplacians;
         std::vector<double>         old_old_temperature_laplacians;
 
         std::vector<double>         gamma_values;
     };
 
     template <int dim>
     TemperatureRHS<dim>::
     TemperatureRHS (const FiniteElement<dim> &temperature_fe,
                     const FiniteElement<dim> &stokes_fe,
                     const Quadrature<dim>    &quadrature)
                     :
                     temperature_fe_values (temperature_fe, quadrature,
                                            update_values    |
                                            update_gradients |
                                            update_hessians  |
                                            update_quadrature_points |
                                            update_JxW_values),
                     stokes_fe_values (stokes_fe, quadrature,
                                       update_values),
                     phi_T (temperature_fe.dofs_per_cell),
                     grad_phi_T (temperature_fe.dofs_per_cell),
 
                     old_velocity_values (quadrature.n_quadrature_points),
                     old_old_velocity_values (quadrature.n_quadrature_points),
 
                     old_temperature_values (quadrature.n_quadrature_points),
                     old_old_temperature_values(quadrature.n_quadrature_points),
                     old_temperature_grads(quadrature.n_quadrature_points),
                     old_old_temperature_grads(quadrature.n_quadrature_points),
                     old_temperature_laplacians(quadrature.n_quadrature_points),
                     old_old_temperature_laplacians(quadrature.n_quadrature_points),
 
                     gamma_values (quadrature.n_quadrature_points)
     {}
 
 
     template <int dim>
     TemperatureRHS<dim>::
     TemperatureRHS (const TemperatureRHS &scratch)
                     :
                     temperature_fe_values (scratch.temperature_fe_values.get_fe(),
                                            scratch.temperature_fe_values.get_quadrature(),
                                            scratch.temperature_fe_values.get_update_flags()),
                     stokes_fe_values (scratch.stokes_fe_values.get_fe(),
                                       scratch.stokes_fe_values.get_quadrature(),
                                       scratch.stokes_fe_values.get_update_flags()),
                     phi_T (scratch.phi_T),
                     grad_phi_T (scratch.grad_phi_T),
 
                     old_velocity_values (scratch.old_velocity_values),
                     old_old_velocity_values (scratch.old_old_velocity_values),
 
                     old_temperature_values (scratch.old_temperature_values),
                     old_old_temperature_values (scratch.old_old_temperature_values),
                     old_temperature_grads (scratch.old_temperature_grads),
                     old_old_temperature_grads (scratch.old_old_temperature_grads),
                     old_temperature_laplacians (scratch.old_temperature_laplacians),
                     old_old_temperature_laplacians (scratch.old_old_temperature_laplacians),
 
                     gamma_values (scratch.gamma_values)
     {}
   }

The CopyData arrays are similar to the Scratch arrays. They provide a constructor, a copy operation, and some arrays for local matrix, local vectors and the relation between local and global degrees of freedom (a.k.a. local_dof_indices).

   namespace CopyData
   {
     template <int dim>
     struct StokesPreconditioner
     {
         StokesPreconditioner (const FiniteElement<dim> &stokes_fe);
         StokesPreconditioner (const StokesPreconditioner &data);
 
         FullMatrix<double>          local_matrix;
         std::vector<unsigned int>   local_dof_indices;
     };
 
     template <int dim>
     StokesPreconditioner<dim>::
     StokesPreconditioner (const FiniteElement<dim> &stokes_fe)
                     :
                     local_matrix (stokes_fe.dofs_per_cell,
                                   stokes_fe.dofs_per_cell),
                     local_dof_indices (stokes_fe.dofs_per_cell)
     {}
 
 
 
     template <int dim>
     StokesPreconditioner<dim>::
     StokesPreconditioner (const StokesPreconditioner &data)
                     :
                     local_matrix (data.local_matrix),
                     local_dof_indices (data.local_dof_indices)
     {}
 
 
 
     template <int dim>
     struct StokesSystem : public StokesPreconditioner<dim>
     {
         StokesSystem (const FiniteElement<dim> &stokes_fe);
         StokesSystem (const StokesSystem<dim> &data);
 
         Vector<double> local_rhs;
     };
 
 
     template <int dim>
     StokesSystem<dim>::
     StokesSystem (const FiniteElement<dim> &stokes_fe)
                     :
                     StokesPreconditioner<dim> (stokes_fe),
                     local_rhs (stokes_fe.dofs_per_cell)
     {}
 
 
     template <int dim>
     StokesSystem<dim>::
     StokesSystem (const StokesSystem<dim> &data)
                     :
                     StokesPreconditioner<dim> (data),
                     local_rhs (data.local_rhs)
     {}
 
 
 
     template <int dim>
     struct TemperatureMatrix
     {
         TemperatureMatrix (const FiniteElement<dim> &temperature_fe);
         TemperatureMatrix (const TemperatureMatrix &data);
 
         FullMatrix<double>          local_mass_matrix;
         FullMatrix<double>          local_stiffness_matrix;
         std::vector<unsigned int>   local_dof_indices;
     };
 
     template <int dim>
     TemperatureMatrix<dim>::
     TemperatureMatrix (const FiniteElement<dim> &temperature_fe)
                     :
                     local_mass_matrix (temperature_fe.dofs_per_cell,
                                        temperature_fe.dofs_per_cell),
                     local_stiffness_matrix (temperature_fe.dofs_per_cell,
                                             temperature_fe.dofs_per_cell),
                     local_dof_indices (temperature_fe.dofs_per_cell)
     {}
 
 
     template <int dim>
     TemperatureMatrix<dim>::
     TemperatureMatrix (const TemperatureMatrix &data)
                     :
                     local_mass_matrix (data.local_mass_matrix),
                     local_stiffness_matrix (data.local_stiffness_matrix),
                     local_dof_indices (data.local_dof_indices)
     {}
 
 
     template <int dim>
     struct TemperatureRHS
     {
         TemperatureRHS (const FiniteElement<dim> &temperature_fe);
         TemperatureRHS (const TemperatureRHS &data);
 
         Vector<double>              local_rhs;
         std::vector<unsigned int>   local_dof_indices;
         FullMatrix<double>          matrix_for_bc;
     };
 
     template <int dim>
     TemperatureRHS<dim>::
     TemperatureRHS (const FiniteElement<dim> &temperature_fe)
                     :
                     local_rhs (temperature_fe.dofs_per_cell),
                     local_dof_indices (temperature_fe.dofs_per_cell),
                     matrix_for_bc (temperature_fe.dofs_per_cell,
                                    temperature_fe.dofs_per_cell)
     {}
 
 
     template <int dim>
     TemperatureRHS<dim>::
     TemperatureRHS (const TemperatureRHS &data)
                     :
                     local_rhs (data.local_rhs),
                     local_dof_indices (data.local_dof_indices),
                     matrix_for_bc (data.matrix_for_bc)
     {}
   }
 }

The BoussinesqFlowProblem class template

This is the declaration of the main class. It is very similar to step-31. Following the task-based parallelization paradigm, we split all the assembly routines into two parts: a first part that can do all the calculations on a certain cell without taking care of other threads, and a second part (which is writing the local data into the global matrices and vectors) which can be entered by only one thread at a time. In order to implement that, we provide functions for each of those two steps for all the four assembly routines that we use in this program.

The pcout (for parallel std::cout) object is used to simplify writing output: each MPI process can use this to generate output as usual, but since each of these processes will produce the same output it will just be replicated many times over; with the ConditionalOStream class, only the output generated by one task will actually be printed to screen, whereas the output by all the other threads will simply be forgotten.

In a bit of naming confusion, you will notice below that some of the variables from namespace TrilinosWrappers are taken from namespace TrilinosWrappers::MPI (such as the right hand side vectors) whereas others are not (such as the various matrices). For the matrices, we happen to use the same class names for parallel and sequential data structures, i.e. all matrices will actually be considered parallel below. On the other hand, for vectors, only those from namespace TrilinosWrappers::MPI are actually distributed. In particular, we will frequently have to query velocities and temperatures at arbitrary quadrature points; consequently, rather than "localizing" a vector whenever we need a localized vector, we solve linear systems in parallel but then immediately localize the solution for further processing. The various *_solution vectors are therefore filled immediately after solving their respective linear system in parallel.

The only other new data member is computing_timer. Its class type, TimerOutput, can be used to conveniently account for compute time spent in certain "sections" of the code that are repeatedly entered. For example, we will enter (and leave) sections for Stokes matrix assembly and would like to accumulate the run time spent in this section over all time steps. At the end of the program, the destructor of the TimerOutput class will automatically produce a nice summary of the times spent in all the sections. For this output, one can choose whether wall clock or CPU times are to be printed, as well as whether we want to produce output every time we leave a section -- which would be quite a lot of additional output -- or just in the end of the program (this choice is made in the from this variable in the results section of this tutorial program.

 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_matrix ();
     void assemble_temperature_system (const double maximal_velocity);
     void project_temperature_field ();
     double get_maximal_velocity () const;
     std::pair<double,double> get_extrapolated_temperature_range () const;
     void solve ();
     void output_results ();
     void refine_mesh (const unsigned int max_grid_level);
 
     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<double>          &old_temperature_laplacians,
                       const std::vector<double>          &old_old_temperature_laplacians,
                       const std::vector<Tensor<1,dim> >  &old_velocity_values,
                       const std::vector<Tensor<1,dim> >  &old_old_velocity_values,
                       const std::vector<double>          &gamma_values,
                       const double                        global_u_infty,
                       const double                        global_T_variation,
                       const double                        cell_diameter) const;
 
 
     ConditionalOStream                  pcout;
 
     Triangulation<dim>                  triangulation;
     double                              global_Omega_diameter;
 
     const unsigned int                  stokes_degree;
     FESystem<dim>                       stokes_fe;
     DoFHandler<dim>                     stokes_dof_handler;
     ConstraintMatrix                    stokes_constraints;
 
     TrilinosWrappers::BlockSparseMatrix stokes_matrix;
     TrilinosWrappers::BlockSparseMatrix stokes_preconditioner_matrix;
 
     TrilinosWrappers::BlockVector       stokes_solution;
     TrilinosWrappers::BlockVector       old_stokes_solution;
     TrilinosWrappers::MPI::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::MPI::Vector       temperature_rhs;
 
 
     double time_step;
     double old_time_step;
     unsigned int timestep_number;
 
     std_cxx1x::shared_ptr<TrilinosWrappers::PreconditionAMG> Amg_preconditioner;
     std_cxx1x::shared_ptr<TrilinosWrappers::PreconditionILU> Mp_preconditioner;
     std_cxx1x::shared_ptr<TrilinosWrappers::PreconditionIC>  T_preconditioner;
 
     bool rebuild_stokes_matrix;
     bool rebuild_stokes_preconditioner;
     bool rebuild_temperature_matrices;
     bool rebuild_temperature_preconditioner;
 
     TimerOutput computing_timer;
 
     void setup_stokes_matrix (const std::vector<IndexSet> &stokes_partitioning);
     void setup_stokes_preconditioner (const std::vector<IndexSet> &stokes_partitioning);
     void setup_temperature_matrices (const IndexSet &temperature_partitioning);
 
     void
     local_assemble_stokes_preconditioner (const typename DoFHandler<dim>::active_cell_iterator &cell,
                                           Assembly::Scratch::StokesPreconditioner<dim> &scratch,
                                           Assembly::CopyData::StokesPreconditioner<dim> &data);
 
     void
     copy_local_to_global_stokes_preconditioner (const Assembly::CopyData::StokesPreconditioner<dim> &data);
 
 
     void
     local_assemble_stokes_system (const typename DoFHandler<dim>::active_cell_iterator &cell,
                                   Assembly::Scratch::StokesSystem<dim>  &scratch,
                                   Assembly::CopyData::StokesSystem<dim> &data);
 
     void
     copy_local_to_global_stokes_system (const Assembly::CopyData::StokesSystem<dim> &data);
 
 
     void
     local_assemble_temperature_matrix (const typename DoFHandler<dim>::active_cell_iterator &cell,
                                        Assembly::Scratch::TemperatureMatrix<dim>  &scratch,
                                        Assembly::CopyData::TemperatureMatrix<dim> &data);
 
     void
     copy_local_to_global_temperature_matrix (const Assembly::CopyData::TemperatureMatrix<dim> &data);
 
 
 
     void
     local_assemble_temperature_rhs (const std::pair<double,double> global_T_range,
                                     const double                   global_max_velocity,
                                     const typename DoFHandler<dim>::active_cell_iterator &cell,
                                     Assembly::Scratch::TemperatureRHS<dim> &scratch,
                                     Assembly::CopyData::TemperatureRHS<dim> &data);
 
     void
     copy_local_to_global_temperature_rhs (const Assembly::CopyData::TemperatureRHS<dim> &data);
 };

BoussinesqFlowProblem class implementation

BoussinesqFlowProblem::BoussinesqFlowProblem

The constructor of the problem is very similar to the constructor in step-31. What is different is the parallel communication: Trilinos uses a message passing interface (MPI) for data distribution. When entering the BoussinesqFlowProblem class, we have to decide how the parallization is to be done. We choose a rather simple strategy and let all processors that are running the program work together, specified by the communicator comm_world(). Next, we create some modified output stream as we already did in step-18. In MPI, all the processors run the same program individually (they simply operate on different chunks of data and exchange some part of that data from time to time). Next, we need to initialize the pcout object in order to print the user information only on one processor. The implementation of this idea is to check the process number when pcout gets a true argument, and it uses the std::cout stream for output. If we are one processor five, for instance, then we will give a false argument to pcout, which means that the output of that processor will not be printed anywhere.

Finally, we enter the preferred options for the TimerOutput object to its constructor. We restrict the output to the pcout stream (processor 0), and then we specify that we want to get a summary table in the end of the program which shows us wallclock times (as opposed to CPU times).

 template <int dim>
 BoussinesqFlowProblem<dim>::BoussinesqFlowProblem ()
                 :
                 pcout (std::cout,
                        (Utilities::System::
                         get_this_mpi_process(MPI_COMM_WORLD)
                         == 0)),
 
                 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_stokes_preconditioner (true),
                 rebuild_temperature_matrices (true),
                 rebuild_temperature_preconditioner (true),
 
                 computing_timer (pcout, TimerOutput::summary,
                                  TimerOutput::wall_times)
 {}

The BoussinesqFlowProblem helper functions

Except two small details, this function is the very same as in step-31. The first detail is actually common to all functions that implement loop over all cells in the triangulation: When operating in parallel, each processor only works on a chunk of cells. This chunk of cells is identified via a so-called subdomain_id, as we also did in step-18. All we need to change is hence to perform the cell-related operations only on the process with the correct ID. The second difference is the way we calculate the maximum value. Before, we could simply have a double variable that we checked against on each quadrature point for each cell. Now, we have to be a bit more careful since each processor only operates on a subset of cells. What we do is to first let each processor calculate the maximum among its cells, and then do a global communication operation called MaxAll that searches for the maximum value among all the maximum values of the individual processors. MPI provides such a call, but it's even simpler to use the respective function of the MPI communicator object since that will do the right thing even if we work without MPI and on a single machine only. The call to MaxAll needs three arguments, namely the local maximum (input), a field for the global maximum (output), and an integer value one that says that we only work on one double.

 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<Tensor<1,dim> > velocity_values(n_q_points);
 
   const FEValuesExtractors::Vector velocities (0);
 
   double max_local_velocity = 0;
 
   typename DoFHandler<dim>::active_cell_iterator
     cell = stokes_dof_handler.begin_active(),
     endc = stokes_dof_handler.end();
   for (; cell!=endc; ++cell)
     if (cell->subdomain_id() ==
         Utilities::System::get_this_mpi_process(MPI_COMM_WORLD))
       {
         fe_values.reinit (cell);
         fe_values[velocities].get_function_values (stokes_solution,
                                                    velocity_values);
 
         for (unsigned int q=0; q<n_q_points; ++q)
           max_local_velocity = std::max (max_local_velocity,
                                          velocity_values[q].norm());
       }
 
   double max_velocity = 0.;
 #ifdef DEAL_II_COMPILER_SUPPORTS_MPI
   MPI_Allreduce (&max_local_velocity, &max_velocity, 1, MPI_DOUBLE,
                  MPI_MAX, MPI_COMM_WORLD);
 #else
   max_velocity = max_local_velocity;
 #endif
 
   return max_velocity;
 }

Again, this is only a slightly modified version of the respective function in step-31. What is new is that each processor works on its partition of cells, and gets a minimum and maximum temperature on that partition. Two global communication steps synchronize the data among the processors.

 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_local_temperature = (1. + time_step/old_time_step) *
                                      old_temperature_solution.linfty_norm()
                                      +
                                      time_step/old_time_step *
                                      old_old_temperature_solution.linfty_norm(),
              max_local_temperature = -min_local_temperature;
 
       typename DoFHandler<dim>::active_cell_iterator
         cell = temperature_dof_handler.begin_active(),
         endc = temperature_dof_handler.end();
       for (; cell!=endc; ++cell)
         if (cell->subdomain_id() ==
             Utilities::System::get_this_mpi_process(MPI_COMM_WORLD))
           {
             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_local_temperature = std::min (min_local_temperature,
                                                   temperature);
                 max_local_temperature = std::max (max_local_temperature,
                                                   temperature);
               }
           }
 
       double min_temperature, max_temperature;
 #ifdef DEAL_II_COMPILER_SUPPORTS_MPI
       MPI_Allreduce (&max_local_temperature, &max_temperature, 1, MPI_DOUBLE,
                      MPI_MAX, MPI_COMM_WORLD);
       MPI_Allreduce (&min_local_temperature, &min_temperature, 1, MPI_DOUBLE,
                      MPI_MIN, MPI_COMM_WORLD);
 #else
       min_temperature = min_local_temperature;
       max_temperature = max_local_temperature;
 #endif
 
       return std::make_pair(min_temperature, max_temperature);
     }
   else
     {
       double min_local_temperature = old_temperature_solution.linfty_norm(),
              max_local_temperature = -min_local_temperature;
 
       typename DoFHandler<dim>::active_cell_iterator
         cell = temperature_dof_handler.begin_active(),
         endc = temperature_dof_handler.end();
       for (; cell!=endc; ++cell)
         if (cell->subdomain_id() ==
             Utilities::System::get_this_mpi_process(MPI_COMM_WORLD))
           {
             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_local_temperature = std::min (min_local_temperature,
                                                   temperature);
                 max_local_temperature = std::max (max_local_temperature,
                                                   temperature);
               }
           }
 
       double min_temperature, max_temperature;
 #ifdef DEAL_II_COMPILER_SUPPORTS_MPI
       MPI_Allreduce (&max_local_temperature, &max_temperature, 1, MPI_DOUBLE,
                      MPI_MAX, MPI_COMM_WORLD);
       MPI_Allreduce (&min_local_temperature, &min_temperature, 1, MPI_DOUBLE,
                      MPI_MIN, MPI_COMM_WORLD);
 #else
       min_temperature = min_local_temperature;
       max_temperature = max_local_temperature;
 #endif
 
       return std::make_pair(min_temperature, max_temperature);
     }
 }

The function that calculates the viscosity is purely local, so this is the same code as in step-31.

 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<double>          &old_temperature_laplacians,
                    const std::vector<double>          &old_old_temperature_laplacians,
                    const std::vector<Tensor<1,dim> >  &old_velocity_values,
                    const std::vector<Tensor<1,dim> >  &old_old_velocity_values,
                    const std::vector<double>          &gamma_values,
                    const double                        global_u_infty,
                    const double                        global_T_variation,
                    const double                        cell_diameter) const
 {
   const double beta = 0.026 * 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)
     {
       const Tensor<1,dim> u = (old_velocity_values[q] +
                                old_old_velocity_values[q]) / 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
                                    * (old_temperature_laplacians[q] +
                                       old_old_temperature_laplacians[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);
     }
 
   if (timestep_number == 0)
     return beta * max_velocity * cell_diameter;
   else
     {
       Assert (old_time_step > 0, ExcInternalError());
 
       const double c_R = 0.11;
       const double global_scaling = c_R * 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 function is new compared to step-31. What is does is to re-implement the library function VectorTools::project() for an MPI-based parallelization, a function we used for generating an initial vector for temperature based on some initial function. The library function only works with shared memory but doesn't know how to utilize multiple machines coupled through MPI to compute the projected solution. If run with more than one MPI process, this would mean that each processor projects the whole field, which is clearly not very efficient. The details of a project() function are not very difficult. All we do is to use a mass matrix and put the evaluation of the initial value function on the right hand side. The mass matrix for temperature we can simply generate using the respective assembly function, so all we need to do here is to create the right hand side and do a CG solve. The assembly function does a loop over all cells and evaluates the function in the EquationData namespace, and does this only on cells pertaining to the respective processor. The implementation of this assembly differs from the assembly we do for the principal assembly functions further down (which include thread-based parallelization with the WorkStream concept). Here we chose to keep things simple (keeping in mind that this function is also only called at the beginning of the program, not every time step), and generating that right hand side is cheap anyway so we won't even notice that this part is not parallized by threads.

Regarding the implementation of inhomogeneous Dirichlet boundary conditions: Since we use the temperature ConstraintMatrix, we can apply the boundary conditions directly when building the respective matrix and right hand side. In this case, the boundary conditions are inhomogeneous, which makes this procedure somewhat tricky. Remember that we get the matrix from some other function. However, the correct imposition of boundary conditions needs the matrix data we work on plus the right hand side simultaneously, since the right hand side is created by Gaussian elimination on the matrix rows. In order to not introduce the matrix assembly at this place, but still having the matrix data available, we choose to create a dummy matrix matrix_for_bc that we only fill with data when we need it for imposing boundary conditions. These positions are exactly those where we have an inhomogeneous entry in the ConstraintMatrix. There are only a few such positions (on the boundary dofs), so it is still much cheaper to use this function than to create the full matrix here. To implement this, we ask the constraint matrix whether the dof under consideration is inhomogeneously constraint. In that case, we generate the respective matrix column that we need for creating the correct right hand side. Note that this (manually generated) matrix entry needs to be exactly the entry that we would fill the matrix with — otherwise, this will not work.

 template <int dim>
 void BoussinesqFlowProblem<dim>::project_temperature_field ()
 {
   assemble_temperature_matrix ();
 
   QGauss<dim> quadrature(temperature_degree+2);
   UpdateFlags update_flags = UpdateFlags(update_values   |
                                          update_quadrature_points |
                                          update_JxW_values);
   FEValues<dim> fe_values (temperature_fe, quadrature, update_flags);
 
   const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
                      n_q_points    = fe_values.n_quadrature_points;
 
   std::vector<unsigned int> local_dof_indices (dofs_per_cell);
   Vector<double> cell_vector (dofs_per_cell);
   FullMatrix<double> matrix_for_bc (dofs_per_cell, dofs_per_cell);
 
   typename DoFHandler<dim>::active_cell_iterator
     cell = temperature_dof_handler.begin_active(),
     endc = temperature_dof_handler.end();
 
   std::vector<double> rhs_values(n_q_points);
 
   TrilinosWrappers::MPI::Vector
     rhs (temperature_mass_matrix.row_partitioner()),
     solution (temperature_mass_matrix.row_partitioner());
 
   for (; cell!=endc; ++cell)
     if (cell->subdomain_id() ==
         Utilities::System::get_this_mpi_process(MPI_COMM_WORLD))
       {
         cell->get_dof_indices (local_dof_indices);
         fe_values.reinit(cell);
 
         EquationData::TemperatureInitialValues<dim>().value_list
           (fe_values.get_quadrature_points(), rhs_values);
 
         cell_vector = 0;
         matrix_for_bc = 0;
         for (unsigned int point=0; point<n_q_points; ++point)
           for (unsigned int i=0; i<dofs_per_cell; ++i)
             {
               cell_vector(i) += rhs_values[point] *
                                 fe_values.shape_value(i,point) *
                                 fe_values.JxW(point);
               if (temperature_constraints.is_inhomogeneously_constrained(local_dof_indices[i]))
                 {
                   for (unsigned int j=0; j<dofs_per_cell; ++j)
                     matrix_for_bc(j,i) += fe_values.shape_value(i,point) *
                                           fe_values.shape_value(j,point) *
                                           fe_values.JxW(point);
                 }
             }
 
         temperature_constraints.distribute_local_to_global (cell_vector,
                                                             local_dof_indices,
                                                             rhs,
                                                             matrix_for_bc);
       }
 
   rhs.compress ();
 
   SolverControl solver_control(5*rhs.size(), 1e-12*rhs.l2_norm());
   SolverCG<TrilinosWrappers::MPI::Vector> cg(solver_control);
 
   TrilinosWrappers::PreconditionSSOR preconditioner_mass;
   preconditioner_mass.initialize(temperature_mass_matrix, 1.3);
 
   cg.solve (temperature_mass_matrix, solution, rhs, preconditioner_mass);
 
   old_temperature_solution = solution;
   temperature_constraints.distribute (old_temperature_solution);
 }

The BoussinesqFlowProblem setup functions

The following three functions set up the Stokes matrix, the matrix used for the Stokes preconditioner, and the temperature matrix. The code is mostly the same as in step-31, but it has been broken out into three functions of their own for simplicity, but also so that they can easily be run in parallel on multiple threads (unless we are running with MPI, in which case this is not possible, as explained in the introduction).

The main functional difference between the code here and that in step-31 is that the matrices we want to set up are distributed across multiple processors. Since we still want to build up the sparsity pattern first for efficiency reasons, we could continue to build the entire sparsity pattern as a BlockCompressedSimpleSparsityPattern, as we did in step-31. However, that would be inefficient: every processor would build the same sparsity pattern, but only initialize a small part of the matrix using it.

Rather, we use an object of type TrilinosWrappers::BlockSparsityPattern, which is (obviously) a wrapper around a sparsity pattern object provided by Trilinos. The advantage is that the Trilinos sparsity pattern class can communicate across multiple processors: if this processor fills in all the nonzero entries that result from the cells it owns, and every other processor does so as well, then at the end after some MPI communication initiated by the compress() call, we will have the globally assembled sparsity pattern available with which the global matrix can be initialized.

The only other change we need to make is to tell the DoFTools::make_sparsity_pattern function that it is only supposed to work on a subset of cells, namely the ones whose subdomain_id equals the number of the current processor, and to ignore all other cells.

This strategy is replicated across all three of the following functions.

Note that Trilinos matrices store the information contained in the sparsity patterns, so we can safely release the sp variable once the matrix has been given the sparsity structure.

 template <int dim>
 void BoussinesqFlowProblem<dim>::
   setup_stokes_matrix (const std::vector<IndexSet> &stokes_partitioning)
 {
   stokes_matrix.clear ();
 
   TrilinosWrappers::BlockSparsityPattern sp (stokes_partitioning,
                                              MPI_COMM_WORLD);
 
   Table<2,DoFTools::Coupling> coupling (dim+1, dim+1);
 
   for (unsigned int c=0; c<dim+1; ++c)
     for (unsigned int d=0; d<dim+1; ++d)
       if (! ((c==dim) && (d==dim)))
         coupling[c][d] = DoFTools::always;
       else
         coupling[c][d] = DoFTools::none;
 
   DoFTools::make_sparsity_pattern (stokes_dof_handler, coupling, sp,
                                    stokes_constraints, false,
                                    Utilities::System::
                                    get_this_mpi_process(MPI_COMM_WORLD));
   sp.compress();
 
   stokes_matrix.reinit (sp);
 }
 
 
 
 template <int dim>
 void BoussinesqFlowProblem<dim>::
   setup_stokes_preconditioner (const std::vector<IndexSet> &stokes_partitioning)
 {
   Amg_preconditioner.reset ();
   Mp_preconditioner.reset ();
 
   stokes_preconditioner_matrix.clear ();
 
   TrilinosWrappers::BlockSparsityPattern sp (stokes_partitioning,
                                              MPI_COMM_WORLD);
 
   Table<2,DoFTools::Coupling> coupling (dim+1, dim+1);
   for (unsigned int c=0; c<dim+1; ++c)
     for (unsigned int d=0; d<dim+1; ++d)
       if (c == d)
         coupling[c][d] = DoFTools::always;
       else
         coupling[c][d] = DoFTools::none;
 
   DoFTools::make_sparsity_pattern (stokes_dof_handler, coupling, sp,
                                    stokes_constraints, false,
                                    Utilities::System::
                                    get_this_mpi_process(MPI_COMM_WORLD));
   sp.compress();
 
   stokes_preconditioner_matrix.reinit (sp);
 }
 
 
 template <int dim>
 void BoussinesqFlowProblem<dim>::
   setup_temperature_matrices (const IndexSet &temperature_partitioner)
 {
   T_preconditioner.reset ();
   temperature_mass_matrix.clear ();
   temperature_stiffness_matrix.clear ();
   temperature_matrix.clear ();
 
   TrilinosWrappers::SparsityPattern sp (temperature_partitioner,
                                         MPI_COMM_WORLD);
   DoFTools::make_sparsity_pattern (temperature_dof_handler, sp,
                                    temperature_constraints, false,
                                    Utilities::System::
                                    get_this_mpi_process(MPI_COMM_WORLD));
   sp.compress();
 
   temperature_matrix.reinit (sp);
   temperature_mass_matrix.reinit (sp);
   temperature_stiffness_matrix.reinit (sp);
 }

The remainder of the setup function (after splitting out the three functions above) mostly has to deal with the things we need to do for parallelization across processors. In particular, at the top it calls GridTools::partition_triangulation to subdivide all cells into subdomains of roughly equal size and roughly minimal interface length (using METIS). We then distribute degrees of freedom for Stokes and temperature DoFHandler objects, and re-sort them in such a way that all degrees of freedom associated with subdomain zero come before all those associated with subdomain one, etc. For the Stokes part, this entails, however, that velocities and pressures become intermixed, but this is trivially solved by sorting again by blocks; it is worth noting that this latter operation leaves the relative ordering of all velocities and pressures alone, i.e. within the velocity block we will still have all those associated with subdomain zero before all velocities associated with subdomain one, etc. This is important since we store each of the blocks of this matrix distributed across all processors and want this to be done in such a way that each processor stores that part of the matrix that is roughly equal to the degrees of freedom located on those cells that it will actually work on. Note how we set boundary conditions on the temperature by using the ConstraintMatrix object.

After this, we have to set up the various partitioners (of type IndexSet, see the introduction) that describe which parts of each matrix or vector will be stored where, then call the functions that actually set up the matrices (concurrently if not using MPI but sequentially otherwise, as explained in the introduction), and at the end also resize the various vectors we keep around in this program. We given those vectors the correct size using the aforementioned Epetra_Map. Most of the vectors are actually localized, i.e., they store all dofs in the problem on each processor. In that case, the only information that is used is the global size. This is different for the two right hand side vectors, which are distributed ones, see also the class declaration.

Note how this function enters and leaves a timed section so that we can get a time report at the end of the program. Note also the use of the pcout variable: to every process it looks like we can write to screen, but only the output of the first processor actually ends up somewhere. We could of course have achieved the same effect by writing to std::cout but would then have had to guard every access to that stream by something like if (Utilities:: System:: get_this_mpi_process (MPI_COMM_WORLD) == 0), hardly a pretty solution.

 template <int dim>
 void BoussinesqFlowProblem<dim>::setup_dofs ()
 {
   computing_timer.enter_section("Setup dof systems");
   std::vector<unsigned int> stokes_sub_blocks (dim+1,0);
   stokes_sub_blocks[dim] = 1;
 
   GridTools::partition_triangulation (Utilities::System::
                                       get_n_mpi_processes(MPI_COMM_WORLD),
                                       triangulation);
 
   {
     stokes_dof_handler.distribute_dofs (stokes_fe);
     DoFRenumbering::subdomain_wise (stokes_dof_handler);
     DoFRenumbering::component_wise (stokes_dof_handler, stokes_sub_blocks);
 
     stokes_constraints.clear ();
     DoFTools::make_hanging_node_constraints (stokes_dof_handler,
                                              stokes_constraints);
 
     std::vector<bool> velocity_mask (dim+1, true);
     velocity_mask[dim] = false;
     VectorTools::interpolate_boundary_values (stokes_dof_handler,
                                               0,
                                               ZeroFunction<dim>(dim+1),
                                               stokes_constraints,
                                               velocity_mask);
 
     std::set<unsigned char> no_normal_flux_boundaries;
     no_normal_flux_boundaries.insert (1);
     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);
     DoFRenumbering::subdomain_wise (temperature_dof_handler);
 
     temperature_constraints.clear ();
     VectorTools::interpolate_boundary_values (temperature_dof_handler,
                                               0,
                                               EquationData::TemperatureInitialValues<dim>(),
                                               temperature_constraints);
     VectorTools::interpolate_boundary_values (temperature_dof_handler,
                                               1,
                                               EquationData::TemperatureInitialValues<dim>(),
                                               temperature_constraints);
     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();
 
   pcout << "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;
 
   std::vector<IndexSet> stokes_partitioning;
   IndexSet temperature_partitioning (n_T);
   {
     const unsigned int my_id =
       Utilities::System::get_this_mpi_process(MPI_COMM_WORLD);
     IndexSet stokes_index_set =
       DoFTools::dof_indices_with_subdomain_association(stokes_dof_handler,
                                                        my_id);
     stokes_partitioning.push_back(stokes_index_set.get_view(0,n_u));
     stokes_partitioning.push_back(stokes_index_set.get_view(n_u,n_u+n_p));
 
     temperature_partitioning =
       DoFTools::dof_indices_with_subdomain_association(temperature_dof_handler,
                                                        my_id);
   }
 
   if (Utilities::System::job_supports_mpi() == false)
     {
       Threads::TaskGroup<> tasks;
       tasks += Threads::new_task (&BoussinesqFlowProblem<dim>::setup_stokes_matrix,
                                   *this,
                                   stokes_partitioning);
       tasks += Threads::new_task (&BoussinesqFlowProblem<dim>::setup_stokes_preconditioner,
                                   *this,
                                   stokes_partitioning);
       tasks += Threads::new_task (&BoussinesqFlowProblem<dim>::setup_temperature_matrices,
                                   *this,
                                   temperature_partitioning);
       tasks.join_all ();
     }
   else
     {
       setup_stokes_matrix (stokes_partitioning);
       setup_stokes_preconditioner (stokes_partitioning);
       setup_temperature_matrices (temperature_partitioning);
     }
 
   stokes_rhs.reinit (stokes_partitioning, MPI_COMM_WORLD);
   stokes_solution.reinit (stokes_rhs);
   old_stokes_solution.reinit (stokes_solution);
 
   temperature_rhs.reinit (temperature_partitioning, MPI_COMM_WORLD);
   temperature_solution.reinit (temperature_rhs);
   old_temperature_solution.reinit (temperature_solution);
   old_old_temperature_solution.reinit (temperature_solution);
 
   computing_timer.exit_section();
 }

The BoussinesqFlowProblem assembly functions

Following the discussion in the introduction and in the Parallel computing with multiple processors accessing shared memory module, we split the assembly functions into different parts:

Stokes preconditioner assembly

Let us start with the functions that builds the Stokes preconditioner. The first two of these are pretty trivial, given the discussion above. Note in particular that the main point in using the scratch data object is that we want to avoid allocating any objects on the free space each time we visit a new cell. As a consequence, the assembly function below only has automatic local variables, and everything else is accessed through the scratch data object, which is allocated only once before we start the loop over all cells:

 template <int dim>
 void
 BoussinesqFlowProblem<dim>::
 local_assemble_stokes_preconditioner (const typename DoFHandler<dim>::active_cell_iterator &cell,
                                       Assembly::Scratch::StokesPreconditioner<dim> &scratch,
                                       Assembly::CopyData::StokesPreconditioner<dim> &data)
 {
   const unsigned int   dofs_per_cell   = stokes_fe.dofs_per_cell;
   const unsigned int   n_q_points      = scratch.stokes_fe_values.n_quadrature_points;
 
   const FEValuesExtractors::Vector velocities (0);
   const FEValuesExtractors::Scalar pressure (dim);
 
   scratch.stokes_fe_values.reinit (cell);
   cell->get_dof_indices (data.local_dof_indices);
 
   data.local_matrix = 0;
 
   for (unsigned int q=0; q<n_q_points; ++q)
     {
       for (unsigned int k=0; k<dofs_per_cell; ++k)
         {
           scratch.grad_phi_u[k] = scratch.stokes_fe_values[velocities].gradient(k,q);
           scratch.phi_p[k]      = scratch.stokes_fe_values[pressure].value (k, q);
         }
 
       for (unsigned int i=0; i<dofs_per_cell; ++i)
         for (unsigned int j=0; j<dofs_per_cell; ++j)
           data.local_matrix(i,j) += (EquationData::eta *
                                      scalar_product (scratch.grad_phi_u[i],
                                                      scratch.grad_phi_u[j])
                                      +
                                      (1./EquationData::eta) *
                                      EquationData::pressure_scaling *
                                      EquationData::pressure_scaling *
                                      scratch.phi_p[i] * scratch.phi_p[j])
                                     * scratch.stokes_fe_values.JxW(q);
     }
 }
 
 
 
 template <int dim>
 void
 BoussinesqFlowProblem<dim>::
 copy_local_to_global_stokes_preconditioner (const Assembly::CopyData::StokesPreconditioner<dim> &data)
 {
   stokes_constraints.distribute_local_to_global (data.local_matrix,
                                                  data.local_dof_indices,
                                                  stokes_preconditioner_matrix);
 }

When we create the WorkStream, we modify the start and end iterator into a so-called SubdomainFilter that tells the individual processes which cells to work on. This is exactly the case discussed in the introduction. Note how we use the construct std_cxx1x::bind to create a function object that is compatible with the WorkStream class. It uses placeholders _1, _2, _3 for the local assembly function that specify cell, scratch data, and copy data, as well as the placeholder _1 for the copy function that expects the data to be written into the global matrix. On the other hand, the implicit zeroth argument of member functions (namely the this pointer of the object on which that member function is to operate on) is bound to the this pointer of the current function. The WorkStream class, as a consequence, does not need to know anything about the object these functions work on.

When the WorkStream is executed, it will create several local assembly routines of the first kind for several cells and let some available processors work on them. The function that needs to be synchronized, i.e., the write operation into the global matrix, however, is executed by only one thread at a time in the prescribed order. Of course, this only holds for the parallelization on a single MPI process. Different MPI processes will have their own WorkStream objects and do that work completely independently. In a distributed calculation, some data will accumulate at degrees of freedom that are not owned by the respective processor. It would be inefficient to send data around every time we encounter such a dof. What happens instead is that the Trilinos sparse matrix will keep that data and send it to the owner at the end of assembly, by calling the compress() command.

 template <int dim>
 void
 BoussinesqFlowProblem<dim>::assemble_stokes_preconditioner ()
 {
   stokes_preconditioner_matrix = 0;
 
   const QGauss<dim> quadrature_formula(stokes_degree+2);
 
   typedef
     FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
     SubdomainFilter;
 
   WorkStream::
     run (SubdomainFilter (IteratorFilters::SubdomainEqualTo
                           (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)),
                           stokes_dof_handler.begin_active()),
          SubdomainFilter (IteratorFilters::SubdomainEqualTo
                           (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)),
                           stokes_dof_handler.end()),
          std_cxx1x::bind (&BoussinesqFlowProblem<dim>::
                           local_assemble_stokes_preconditioner,
                           this,
                           _1,
                           _2,
                           _3),
          std_cxx1x::bind (&BoussinesqFlowProblem<dim>::
                           copy_local_to_global_stokes_preconditioner,
                           this,
                           _1),
          Assembly::Scratch::
          StokesPreconditioner<dim> (stokes_fe, quadrature_formula,
                                     update_JxW_values |
                                     update_values |
                                     update_gradients),
          Assembly::CopyData::
          StokesPreconditioner<dim> (stokes_fe));
 
   stokes_preconditioner_matrix.compress();
 }

The final function in this block initiates assemble of the Stokes preconditioner matrix and then builds the Stokes preconditioner. It is mostly the same as in the serial case. The only difference to step-31 is that we use an ILU preconditioner for the pressure mass matrix instead of IC, as discussed in the introduction.

 template <int dim>
 void
 BoussinesqFlowProblem<dim>::build_stokes_preconditioner ()
 {
   if (rebuild_stokes_preconditioner == false)
     return;
 
   computing_timer.enter_section ("   Build Stokes preconditioner");
   pcout << "   Rebuilding Stokes preconditioner..." << std::flush;
 
   assemble_stokes_preconditioner ();
 
   std::vector<std::vector<bool> > constant_modes;
   std::vector<bool>  velocity_components (dim+1,true);
   velocity_components[dim] = false;
   DoFTools::extract_constant_modes (stokes_dof_handler, velocity_components,
                                     constant_modes);
 
   Mp_preconditioner  = std_cxx1x::shared_ptr<TrilinosWrappers::PreconditionILU>
                        (new TrilinosWrappers::PreconditionILU());
   Amg_preconditioner = std_cxx1x::shared_ptr<TrilinosWrappers::PreconditionAMG>
                        (new TrilinosWrappers::PreconditionAMG());
 
   TrilinosWrappers::PreconditionAMG::AdditionalData Amg_data;
   Amg_data.constant_modes = constant_modes;
   Amg_data.elliptic = true;
   Amg_data.higher_order_elements = true;
   Amg_data.smoother_sweeps = 2;
   Amg_data.aggregation_threshold = 0.02;
 
   Mp_preconditioner->initialize (stokes_preconditioner_matrix.block(1,1));
   Amg_preconditioner->initialize (stokes_preconditioner_matrix.block(0,0),
                                   Amg_data);
 
   rebuild_stokes_preconditioner = false;
 
   pcout << std::endl;
   computing_timer.exit_section();
 }

Stokes system assembly

The next three functions implement the assembly of the Stokes system, again split up into a part performing local calculations, one for writing the local data into the global matrix and vector, and one for actually running the loop over all cells with the help of the WorkStream class. Note that the assembly of the Stokes matrix needs only to be done in case we have changed the mesh. Otherwise, just the (temperature-dependent) right hand side needs to be calculated here. Since we are working with distributed matrices and vectors, we have to call the respective compress() functions in the end of the assembly in order to send non-local data to the owner process.

 template <int dim>
 void
 BoussinesqFlowProblem<dim>::
 local_assemble_stokes_system (const typename DoFHandler<dim>::active_cell_iterator &cell,
                               Assembly::Scratch::StokesSystem<dim> &scratch,
                               Assembly::CopyData::StokesSystem<dim> &data)
 {
   const unsigned int dofs_per_cell = scratch.stokes_fe_values.get_fe().dofs_per_cell;
   const unsigned int n_q_points    = scratch.stokes_fe_values.n_quadrature_points;
 
   const FEValuesExtractors::Vector velocities (0);
   const FEValuesExtractors::Scalar pressure (dim);
 
   scratch.stokes_fe_values.reinit (cell);
 
   typename DoFHandler<dim>::active_cell_iterator
     temperature_cell (&triangulation,
                       cell->level(),
                       cell->index(),
                       &temperature_dof_handler);
   scratch.temperature_fe_values.reinit (temperature_cell);
 
   if (rebuild_stokes_matrix)
     data.local_matrix = 0;
   data.local_rhs = 0;
 
   scratch.temperature_fe_values.get_function_values (old_temperature_solution,
                                                      scratch.old_temperature_values);
 
   for (unsigned int q=0; q<n_q_points; ++q)
     {
       const double old_temperature = scratch.old_temperature_values[q];
 
       for (unsigned int k=0; k<dofs_per_cell; ++k)
         {
           scratch.phi_u[k] = scratch.stokes_fe_values[velocities].value (k,q);
           if (rebuild_stokes_matrix)
             {
               scratch.grads_phi_u[k] = scratch.stokes_fe_values[velocities].symmetric_gradient(k,q);
               scratch.div_phi_u[k]   = scratch.stokes_fe_values[velocities].divergence (k, q);
               scratch.phi_p[k]       = scratch.stokes_fe_values[pressure].value (k, q);
             }
         }
 
       if (rebuild_stokes_matrix)
         for (unsigned int i=0; i<dofs_per_cell; ++i)
           for (unsigned int j=0; j<dofs_per_cell; ++j)
             data.local_matrix(i,j) += (EquationData::eta * 2 *
                                        (scratch.grads_phi_u[i] * scratch.grads_phi_u[j])
                                        - (EquationData::pressure_scaling *
                                           scratch.div_phi_u[i] * scratch.phi_p[j])
                                        - (EquationData::pressure_scaling *
                                           scratch.phi_p[i] * scratch.div_phi_u[j]))
                                       * scratch.stokes_fe_values.JxW(q);
 
       const Tensor<1,dim>
         gravity = EquationData::gravity_vector (scratch.stokes_fe_values
                                                 .quadrature_point(q));
 
       for (unsigned int i=0; i<dofs_per_cell; ++i)
         data.local_rhs(i) += (EquationData::density(old_temperature) *
                               gravity  *
                               scratch.phi_u[i]) *
                              scratch.stokes_fe_values.JxW(q);
     }
 
   cell->get_dof_indices (data.local_dof_indices);
 }
 
 
 
 template <int dim>
 void
 BoussinesqFlowProblem<dim>::
 copy_local_to_global_stokes_system (const Assembly::CopyData::StokesSystem<dim> &data)
 {
   if (rebuild_stokes_matrix == true)
     stokes_constraints.distribute_local_to_global (data.local_matrix,
                                                    data.local_rhs,
                                                    data.local_dof_indices,
                                                    stokes_matrix,
                                                    stokes_rhs);
   else
     stokes_constraints.distribute_local_to_global (data.local_rhs,
                                                    data.local_dof_indices,
                                                    stokes_rhs);
 }
 
 
 
 template <int dim>
 void BoussinesqFlowProblem<dim>::assemble_stokes_system ()
 {
   pcout << "   Assembling..." << std::flush;
 
   computing_timer.enter_section ("   Assemble Stokes system");
 
   if (rebuild_stokes_matrix == true)
     stokes_matrix=0;
 
   stokes_rhs=0;
 
   const QGauss<dim> quadrature_formula(stokes_degree+2);
 
   typedef
     FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
     SubdomainFilter;
 
   WorkStream::
     run (SubdomainFilter (IteratorFilters::SubdomainEqualTo
                           (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)),
                           stokes_dof_handler.begin_active()),
          SubdomainFilter (IteratorFilters::SubdomainEqualTo
                           (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)),
                           stokes_dof_handler.end()),
          std_cxx1x::bind (&BoussinesqFlowProblem<dim>::
                           local_assemble_stokes_system,
                           this,
                           _1,
                           _2,
                           _3),
          std_cxx1x::bind (&BoussinesqFlowProblem<dim>::
                           copy_local_to_global_stokes_system,
                           this,
                           _1),
          Assembly::Scratch::
          StokesSystem<dim> (stokes_fe, quadrature_formula,
                             (update_values    |
                              update_quadrature_points  |
                              update_JxW_values |
                              (rebuild_stokes_matrix == true
                               ?
                               update_gradients
                               :
                               UpdateFlags(0))),
                             temperature_fe,
                             update_values),
          Assembly::CopyData::
          StokesSystem<dim> (stokes_fe));
 
   stokes_matrix.compress();
   stokes_rhs.compress();
 
   rebuild_stokes_matrix = false;
 
   pcout << std::endl;
   computing_timer.exit_section();
 }

Temperature matrix assembly

The task to be performed by the next three functions is to calculate a mass matrix and a Laplace matrix on the temperature system. These will be combined in order to yield the semi-implicit time stepping matrix that consists of the mass matrix plus a time step weight times the Laplace matrix. This function is again essentially the body of the loop over all cells from step-31.

The two following functions perform similar services as the ones above.

 template <int dim>
 void BoussinesqFlowProblem<dim>::
 local_assemble_temperature_matrix (const typename DoFHandler<dim>::active_cell_iterator &cell,
                                    Assembly::Scratch::TemperatureMatrix<dim> &scratch,
                                    Assembly::CopyData::TemperatureMatrix<dim> &data)
 {
   const unsigned int dofs_per_cell = scratch.temperature_fe_values.get_fe().dofs_per_cell;
   const unsigned int n_q_points    = scratch.temperature_fe_values.n_quadrature_points;
 
   scratch.temperature_fe_values.reinit (cell);
   cell->get_dof_indices (data.local_dof_indices);
 
   data.local_mass_matrix = 0;
   data.local_stiffness_matrix = 0;
 
   for (unsigned int q=0; q<n_q_points; ++q)
     {
       for (unsigned int k=0; k<dofs_per_cell; ++k)
         {
           scratch.grad_phi_T[k] = scratch.temperature_fe_values.shape_grad (k,q);
           scratch.phi_T[k]      = scratch.temperature_fe_values.shape_value (k, q);
         }
 
       for (unsigned int i=0; i<dofs_per_cell; ++i)
         for (unsigned int j=0; j<dofs_per_cell; ++j)
           {
             data.local_mass_matrix(i,j)
               += (scratch.phi_T[i] * scratch.phi_T[j]
                   *
                   scratch.temperature_fe_values.JxW(q));
             data.local_stiffness_matrix(i,j)
               += (EquationData::kappa * scratch.grad_phi_T[i] * scratch.grad_phi_T[j]
                   *
                   scratch.temperature_fe_values.JxW(q));
           }
     }
 }
 
 
 
 template <int dim>
 void
 BoussinesqFlowProblem<dim>::
 copy_local_to_global_temperature_matrix (const Assembly::CopyData::TemperatureMatrix<dim> &data)
 {
   temperature_constraints.distribute_local_to_global (data.local_mass_matrix,
                                                       data.local_dof_indices,
                                                       temperature_mass_matrix);
   temperature_constraints.distribute_local_to_global (data.local_stiffness_matrix,
                                                       data.local_dof_indices,
                                                       temperature_stiffness_matrix);
 }
 
 
 template <int dim>
 void BoussinesqFlowProblem<dim>::assemble_temperature_matrix ()
 {
   if (rebuild_temperature_matrices == false)
     return;
 
   computing_timer.enter_section ("   Assemble temperature matrices");
   temperature_mass_matrix = 0;
   temperature_stiffness_matrix = 0;
 
   const QGauss<dim> quadrature_formula(temperature_degree+2);
 
   typedef
     FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
     SubdomainFilter;
 
   WorkStream::
     run (SubdomainFilter (IteratorFilters::SubdomainEqualTo
                           (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)),
                           temperature_dof_handler.begin_active()),
          SubdomainFilter (IteratorFilters::SubdomainEqualTo
                           (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)),
                           temperature_dof_handler.end()),
          std_cxx1x::bind (&BoussinesqFlowProblem<dim>::
                           local_assemble_temperature_matrix,
                           this,
                           _1,
                           _2,
                           _3),
          std_cxx1x::bind (&BoussinesqFlowProblem<dim>::
                           copy_local_to_global_temperature_matrix,
                           this,
                           _1),
          Assembly::Scratch::
          TemperatureMatrix<dim> (temperature_fe, quadrature_formula),
          Assembly::CopyData::
          TemperatureMatrix<dim> (temperature_fe));
 
   temperature_mass_matrix.compress();
   temperature_stiffness_matrix.compress();
 
   rebuild_temperature_matrices = false;
   rebuild_temperature_preconditioner = true;
 
   computing_timer.exit_section();
 }

Temperature right hand side assembly

This is the last assembly function. It calculates the right hand side of the temperature system, which includes the convection and the stabilization terms. It includes a lot of evaluations of old solutions at the quadrature points (which are necessary for calculating the artificial viscosity of stabilization), but is otherwise similar to the other assembly functions. Notice, once again, how we resolve the dilemma of having inhomogeneous boundary conditions, but just making a right hand side at this point (compare the comments for the project function): We create some matrix columns with exactly the values that would be entered for the temperature stiffness matrix, in case we have inhomogeneously constrained dofs. That will account for the correct balance of the right hand side vector with the matrix system of temperature.

 template <int dim>
 void BoussinesqFlowProblem<dim>::
 local_assemble_temperature_rhs (const std::pair<double,double> global_T_range,
                                 const double                   global_max_velocity,
                                 const typename DoFHandler<dim>::active_cell_iterator &cell,
                                 Assembly::Scratch::TemperatureRHS<dim> &scratch,
                                 Assembly::CopyData::TemperatureRHS<dim> &data)
 {
   const bool use_bdf2_scheme = (timestep_number != 0);
 
   const unsigned int dofs_per_cell = scratch.temperature_fe_values.get_fe().dofs_per_cell;
   const unsigned int n_q_points    = scratch.temperature_fe_values.n_quadrature_points;
 
   EquationData::TemperatureRightHandSide<dim>  temperature_right_hand_side;
 
   const FEValuesExtractors::Vector velocities (0);
 
   data.local_rhs = 0;
   data.matrix_for_bc = 0;
   cell->get_dof_indices (data.local_dof_indices);
 
   scratch.temperature_fe_values.reinit (cell);
 
   typename DoFHandler<dim>::active_cell_iterator
     stokes_cell (&triangulation,
                  cell->level(),
                  cell->index(),
                  &stokes_dof_handler);
   scratch.stokes_fe_values.reinit (stokes_cell);
 
   scratch.temperature_fe_values.get_function_values (old_temperature_solution,
                                                      scratch.old_temperature_values);
   scratch.temperature_fe_values.get_function_values (old_old_temperature_solution,
                                                      scratch.old_old_temperature_values);
 
   scratch.temperature_fe_values.get_function_gradients (old_temperature_solution,
                                                         scratch.old_temperature_grads);
   scratch.temperature_fe_values.get_function_gradients (old_old_temperature_solution,
                                                         scratch.old_old_temperature_grads);
 
   scratch.temperature_fe_values.get_function_laplacians (old_temperature_solution,
                                                          scratch.old_temperature_laplacians);
   scratch.temperature_fe_values.get_function_laplacians (old_old_temperature_solution,
                                                          scratch.old_old_temperature_laplacians);
 
   temperature_right_hand_side.value_list (scratch.temperature_fe_values.get_quadrature_points(),
                                           scratch.gamma_values);
 
   scratch.stokes_fe_values[velocities].get_function_values (stokes_solution,
                                                             scratch.old_velocity_values);
   scratch.stokes_fe_values[velocities].get_function_values (old_stokes_solution,
                                                             scratch.old_old_velocity_values);
 
   const double nu
     = compute_viscosity (scratch.old_temperature_values,
                          scratch.old_old_temperature_values,
                          scratch.old_temperature_grads,
                          scratch.old_old_temperature_grads,
                          scratch.old_temperature_laplacians,
                          scratch.old_old_temperature_laplacians,
                          scratch.old_velocity_values,
                          scratch.old_old_velocity_values,
                          scratch.gamma_values,
                          global_max_velocity,
                          global_T_range.second - global_T_range.first,
                          cell->diameter());
 
   for (unsigned int q=0; q<n_q_points; ++q)
     {
       for (unsigned int k=0; k<dofs_per_cell; ++k)
         {
           scratch.phi_T[k]      = scratch.temperature_fe_values.shape_value (k, q);
           scratch.grad_phi_T[k] = scratch.temperature_fe_values.shape_grad (k,q);
         }
 
 
       const double old_Ts
         = (use_bdf2_scheme ?
            (scratch.old_temperature_values[q] *
             (time_step + old_time_step) / old_time_step
             -
             scratch.old_old_temperature_values[q] *
             (time_step * time_step) /
             (old_time_step * (time_step + old_time_step)))
            :
            scratch.old_temperature_values[q]);
 
       const Tensor<1,dim> ext_grad_T
         = (use_bdf2_scheme ?
            (scratch.old_temperature_grads[q] *
             (1+time_step/old_time_step)
             -
             scratch.old_old_temperature_grads[q] *
             time_step / old_time_step)
            :
            scratch.old_temperature_grads[q]);
 
       const Tensor<1,dim> extrapolated_u
         = (use_bdf2_scheme ?
            (scratch.old_velocity_values[q] * (1+time_step/old_time_step) -
             scratch.old_old_velocity_values[q] * time_step/old_time_step)
            :
            scratch.old_velocity_values[q]);
 
       for (unsigned int i=0; i<dofs_per_cell; ++i)
         {
           data.local_rhs(i) += (old_Ts * scratch.phi_T[i]
                                 -
                                 time_step *
                                 extrapolated_u * ext_grad_T * scratch.phi_T[i]
                                 -
                                 time_step *
                                 nu * ext_grad_T * scratch.grad_phi_T[i]
                                 +
                                 time_step *
                                 scratch.gamma_values[q] * scratch.phi_T[i])
                                *
                                scratch.temperature_fe_values.JxW(q);
 
           if (temperature_constraints.is_inhomogeneously_constrained(data.local_dof_indices[i]))
             {
               for (unsigned int j=0; j<dofs_per_cell; ++j)
                 data.matrix_for_bc(j,i) += (scratch.phi_T[i] * scratch.phi_T[j] *
                                             (use_bdf2_scheme ?
                                              ((2*time_step + old_time_step) /
                                               (time_step + old_time_step)) : 1.)
                                             +
                                             scratch.grad_phi_T[i] *
                                             scratch.grad_phi_T[j] *
                                             EquationData::kappa *
                                             time_step)
                                            *
                                            scratch.temperature_fe_values.JxW(q);
             }
         }
     }
 }
 
 
 template <int dim>
 void
 BoussinesqFlowProblem<dim>::
 copy_local_to_global_temperature_rhs (const Assembly::CopyData::TemperatureRHS<dim> &data)
 {
   temperature_constraints.distribute_local_to_global (data.local_rhs,
                                                       data.local_dof_indices,
                                                       temperature_rhs,
                                                       data.matrix_for_bc);
 }

In the function that runs the WorkStream for actually calculating the right hand side, we also generate the final matrix. As mentioned above, it is a sum of the mass matrix and the Laplace matrix, times some time step weight. This weight is specified by the BDF-2 time integration scheme, see the introduction in step-31. What is new in this tutorial program (in addition to the use of MPI parallelization and the WorkStream class), is that we now precompute the temperature preconditioner as well. The reason is that the setup of the IC preconditioner takes a noticable time compared to the solver because we usually only need between 10 and 20 iterations for solving the temperature system. Hence, it is more efficient to precompute the preconditioner, even though the matrix entries may slightly change because the time step might change. This is not too big a problem because we remesh every fifth time step (and regenerate the preconditioner then).

 template <int dim>
 void BoussinesqFlowProblem<dim>::assemble_temperature_system (const double maximal_velocity)
 {
   const bool use_bdf2_scheme = (timestep_number != 0);
 
   if (use_bdf2_scheme == true)
     {
       temperature_matrix.copy_from (temperature_mass_matrix);
       temperature_matrix *= (2*time_step + old_time_step) /
                             (time_step + old_time_step);
       temperature_matrix.add (time_step, temperature_stiffness_matrix);
     }
   else
     {
       temperature_matrix.copy_from (temperature_mass_matrix);
       temperature_matrix.add (time_step, temperature_stiffness_matrix);
     }
   temperature_matrix.compress();
 
   if (rebuild_temperature_preconditioner == true)
     {
       T_preconditioner =  std_cxx1x::shared_ptr<TrilinosWrappers::PreconditionIC>
                           (new TrilinosWrappers::PreconditionIC());
       T_preconditioner->initialize (temperature_matrix);
       rebuild_temperature_preconditioner = false;
     }
 
   temperature_rhs = 0;
 
   const QGauss<dim> quadrature_formula(temperature_degree+2);
   const std::pair<double,double>
     global_T_range = get_extrapolated_temperature_range();
 
   typedef
     FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
     SubdomainFilter;
 
   WorkStream::
     run (SubdomainFilter (IteratorFilters::SubdomainEqualTo
                           (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)),
                           temperature_dof_handler.begin_active()),
          SubdomainFilter (IteratorFilters::SubdomainEqualTo
                           (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)),
                           temperature_dof_handler.end()),
          std_cxx1x::bind (&BoussinesqFlowProblem<dim>::
                           local_assemble_temperature_rhs,
                           this,
                           global_T_range,
                           maximal_velocity,
                           _1,
                           _2,
                           _3),
          std_cxx1x::bind (&BoussinesqFlowProblem<dim>::
                           copy_local_to_global_temperature_rhs,
                           this,
                           _1),
          Assembly::Scratch::
          TemperatureRHS<dim> (temperature_fe, stokes_fe, quadrature_formula),
          Assembly::CopyData::
          TemperatureRHS<dim> (temperature_fe));
 
   temperature_rhs.compress();
 }

BoussinesqFlowProblem::solve

This function solves the linear systems in each time step of the Boussinesq problem. First, we work on the Stokes system and then on the temperature system. In essence, it does the same things as the respective function in step-31. However, there are a few things that we need to pay some attention to. The first thing is, as mentioned in the introduction, the way we store our solution: we keep the full vector with all degrees of freedom on each MPI node. When we enter a solver which is supposed to perform matrix-vector products with a distributed matrix, this is not the appropriate form, though. There, we will want to have the solution vector to be distributed in the same way as the matrix. So what we do first (after initializing the Schur-complement based preconditioner) is to generate a distributed vector called distributed_stokes_solution and put only the locally owned dofs into that, which is neatly done by the operator= of the Trilinos vector. Next, we need to set the pressure values at hanging nodes to zero. This we also did in step-31 in order not to disturb the Schur complement by some vector entries that actually are irrelevant during the solve stage. As a difference to step-31, here we do it only for the locally owned pressure dofs. After solving for the Stokes solution, each processor copies distributed solution back into the solution vector for which every element is locally owned.

Apart from these two changes, everything is the same as in step-31, so we don't need to further comment on it.

 template <int dim>
 void BoussinesqFlowProblem<dim>::solve ()
 {
   computing_timer.enter_section ("   Solve Stokes system");
   pcout << "   Solving..." << std::endl;
 
   {
     const LinearSolvers::BlockSchurPreconditioner<TrilinosWrappers::PreconditionAMG,
                                                   TrilinosWrappers::PreconditionILU>
       preconditioner (stokes_matrix, *Mp_preconditioner, *Amg_preconditioner);
 
     TrilinosWrappers::MPI::BlockVector
       distributed_stokes_solution (stokes_rhs);
     distributed_stokes_solution = stokes_solution;
 
     const unsigned int
       start = (distributed_stokes_solution.block(0).size() +
                distributed_stokes_solution.block(1).local_range().first),
       end   = (distributed_stokes_solution.block(0).size() +
                distributed_stokes_solution.block(1).local_range().second);
     for (unsigned int i=start; i<end; ++i)
       if (stokes_constraints.is_constrained (i))
         distributed_stokes_solution(i) = 0;
 
     SolverControl solver_control (stokes_matrix.m(), 1e-22*stokes_rhs.l2_norm());
     SolverBicgstab<TrilinosWrappers::MPI::BlockVector>
       bicgstab (solver_control, false);
 
     bicgstab.solve(stokes_matrix, distributed_stokes_solution, stokes_rhs,
                    preconditioner);
 
     stokes_solution = distributed_stokes_solution;
 
     pcout << "   "
           << solver_control.last_step()
           << " BiCGStab iterations for Stokes subsystem."
           << std::endl;
 
     stokes_constraints.distribute (stokes_solution);
   }
   computing_timer.exit_section();
 
 
   computing_timer.enter_section ("   Assemble temperature rhs");
 
   old_time_step = time_step;
   const double maximal_velocity = get_maximal_velocity();
   if (maximal_velocity > 1e-10)
     time_step = 1./(1.6*dim*std::sqrt(1.*dim)) /
                 temperature_degree *
                 GridTools::minimal_cell_diameter(triangulation) /
                 maximal_velocity;
   else
     time_step = 1./(1.6*dim*std::sqrt(1.*dim)) /
                 temperature_degree *
                 GridTools::minimal_cell_diameter(triangulation) /
                 1e-10;
 
   pcout << "   Maximal velocity: "
         << maximal_velocity * EquationData::year_in_seconds * 100
         << " cm/year"
         << std::endl;
   pcout << "   " << "Time step: "
         << time_step/EquationData::year_in_seconds
         << " years"
         << std::endl;
 
   temperature_solution = old_temperature_solution;
 
   assemble_temperature_system (maximal_velocity);
 
   computing_timer.exit_section ();
 
   computing_timer.enter_section ("   Solve temperature system");
 
   {
     SolverControl solver_control (temperature_matrix.m(),
                                   1e-12*temperature_rhs.l2_norm());
     SolverCG<TrilinosWrappers::MPI::Vector>   cg (solver_control);
 
     TrilinosWrappers::MPI::Vector
       distributed_temperature_solution (temperature_rhs);
     distributed_temperature_solution = temperature_solution;
 
     cg.solve (temperature_matrix, distributed_temperature_solution,
               temperature_rhs, *T_preconditioner);
 
     temperature_solution = distributed_temperature_solution;
     temperature_constraints.distribute (temperature_solution);
 
     pcout << "   "
           << solver_control.last_step()
           << " CG iterations for temperature" << std::endl;
     computing_timer.exit_section();
 
     double min_temperature = temperature_solution(0),
            max_temperature = temperature_solution(0);
     for (unsigned int i=1; i<temperature_solution.size(); ++i)
       {
         min_temperature = std::min<double> (min_temperature,
                                             temperature_solution(i));
         max_temperature = std::max<double> (max_temperature,
                                             temperature_solution(i));
       }
 
     pcout << "   Temperature range: "
           << min_temperature << ' ' << max_temperature
           << std::endl;
   }
 }

BoussinesqFlowProblem::output_results

This function does mostly what the corresponding one did in to step-31, in particular merging data from the two DoFHandler objects (for the Stokes and the temperature parts of the problem) into one is the same. There are three minor changes: we make sure that only a single processor actually does some work here; take care of scaling variables in a useful way; and in addition to the Stokes and temperature parts in the joint_fe finite element, we also add a piecewise constant field that denotes the subdomain id a cell corresponds to. This allows us to visualize the partitioning of the domain. As a consequence, we also have to change the assertion about the number of degrees of freedom in the joint DoFHandler object (which is now equal to the number of Stokes degrees of freedom plus the temperature degrees of freedom plus the number of active cells as that is the number of partition variables we want to add), and adjust the number of elements in the arrays we use to name the components of the joint solution vector and to identify which of these components are scalars or parts of dim-dimensional vectors.

As for scaling: as mentioned in the introduction, to keep the Stokes equations properly scaled and symmetric, we introduced a new pressure $\hat p = \frac{L}{\eta}p$. What we really wanted, however, was the original pressure $p$, so while copying data from the Stokes DoFHandler into the joint one, we undo this scaling. While we're at it messing with the results of the simulation, we do two more things: First, the pressure is only defined up to a constant. To make it more easily comparable, we compute the minimal value of the pressure computed and shift all values up by that amount -- in essence making all pressure variables positive or zero. Secondly, let's also take care of the awkward units we use for the velocity: it is computed in SI units of meters per second, which of course is a very small number in the earth mantle. We therefore rescale things into centimeters per year, the unit commonly used in geophysics.

 template <int dim>
 void BoussinesqFlowProblem<dim>::output_results ()
 {
   if (timestep_number % 25 != 0)
     return;
 
   computing_timer.enter_section ("Postprocessing");
 
   if (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD) == 0)
     {
 
       const FESystem<dim> joint_fe (stokes_fe, 1,
                                     temperature_fe, 1,
                                     FE_DGQ<dim>(0), 1);
       DoFHandler<dim> joint_dof_handler (triangulation);
       joint_dof_handler.distribute_dofs (joint_fe);
       Assert (joint_dof_handler.n_dofs() ==
               stokes_dof_handler.n_dofs() +
               temperature_dof_handler.n_dofs() +
               triangulation.n_active_cells(),
               ExcInternalError());
 
       Vector<double> joint_solution (joint_dof_handler.n_dofs());
 
       {
         double minimal_pressure = stokes_solution.block(1)(0);
         for (unsigned int i=0; i<stokes_solution.block(1).size(); ++i)
           minimal_pressure = std::min<double> (stokes_solution.block(1)(i),
                                                minimal_pressure);
 
         std::vector<unsigned int> local_joint_dof_indices (joint_fe.dofs_per_cell);
         std::vector<unsigned int> local_stokes_dof_indices (stokes_fe.dofs_per_cell);
         std::vector<unsigned int> local_temperature_dof_indices (temperature_fe.dofs_per_cell);
 
         typename DoFHandler<dim>::active_cell_iterator
           joint_cell       = joint_dof_handler.begin_active(),
           joint_endc       = joint_dof_handler.end(),
           stokes_cell      = stokes_dof_handler.begin_active(),
           temperature_cell = temperature_dof_handler.begin_active();
         for (; joint_cell!=joint_endc; ++joint_cell, ++stokes_cell, ++temperature_cell)
           {
             joint_cell->get_dof_indices (local_joint_dof_indices);
             stokes_cell->get_dof_indices (local_stokes_dof_indices);
             temperature_cell->get_dof_indices (local_temperature_dof_indices);
 
             for (unsigned int i=0; i<joint_fe.dofs_per_cell; ++i)
               if (joint_fe.system_to_base_index(i).first.first == 0)
                 {
                   Assert (joint_fe.system_to_base_index(i).second
                           <
                           local_stokes_dof_indices.size(),
                           ExcInternalError());
 
                   const unsigned int index_in_stokes_fe
                     = joint_fe.system_to_base_index(i).second;
 
                   if (stokes_fe.system_to_component_index(index_in_stokes_fe).first
                       < dim)
                     {
                       joint_solution(local_joint_dof_indices[i])
                         = (stokes_solution(local_stokes_dof_indices
                                            [joint_fe.system_to_base_index(i).second])
                            *
                            EquationData::year_in_seconds
                            *
                            100);
                     }
                   else
                     {
                       Assert (stokes_fe.system_to_component_index(index_in_stokes_fe).first
                               ==
                               dim,
                               ExcInternalError());
 
                       joint_solution(local_joint_dof_indices[i])
                         = ((stokes_solution(local_stokes_dof_indices
                                             [joint_fe.system_to_base_index(i).second])
                             -
                             minimal_pressure)
                            *
                            EquationData::pressure_scaling);
                     }
                 }
               else if (joint_fe.system_to_base_index(i).first.first == 1)
                 {
                   Assert (joint_fe.system_to_base_index(i).second
                           <
                           local_temperature_dof_indices.size(),
                           ExcInternalError());
                   joint_solution(local_joint_dof_indices[i])
                     = temperature_solution(local_temperature_dof_indices
                                            [joint_fe.system_to_base_index(i).second]);
                 }
               else
                 {
                   Assert (joint_fe.system_to_base_index(i).first.first == 2,
                           ExcInternalError());
                   Assert (joint_fe.system_to_base_index(i).second
                           == 0,
                           ExcInternalError());
                   joint_solution(local_joint_dof_indices[i])
                     = joint_cell->subdomain_id();
                 }
           }
       }
 
       std::vector<std::string> joint_solution_names (dim, "velocity");
       joint_solution_names.push_back ("p");
       joint_solution_names.push_back ("T");
       joint_solution_names.push_back ("partition");
 
       DataOut<dim> data_out;
 
       data_out.attach_dof_handler (joint_dof_handler);
 
       std::vector<DataComponentInterpretation::DataComponentInterpretation>
         data_component_interpretation
         (dim+3, DataComponentInterpretation::component_is_scalar);
       for (unsigned int i=0; i<dim; ++i)
         data_component_interpretation[i]
           = DataComponentInterpretation::component_is_part_of_vector;
 
       data_out.add_data_vector (joint_solution, joint_solution_names,
                                 DataOut<dim>::type_dof_data,
                                 data_component_interpretation);
       data_out.build_patches (std::min(stokes_degree, temperature_degree));
 
       std::ostringstream filename;
       filename << "solution-" << Utilities::int_to_string(timestep_number, 5)
                << ".vtk";
 
       std::ofstream output (filename.str().c_str());
       data_out.write_vtk (output);
     }
 
   computing_timer.exit_section ();
 }

BoussinesqFlowProblem::refine_mesh

This function isn't really new either. Since the setup_dofs function that we call in the middle has its own timer section, we split timing this function into two sections. It will also allow us to easily identify which of the two is more expensive.

One thing of note, however, is that we don't want to compute all error indicators on all cells, of course. Rather, it would be nice if each processor could only compute the error indicators for those cells it actually owns. However, in order for mesh refinement to proceed in the same way on all processors, all processors would have to exchange their refinement indicators. We do so in two steps: first, we call the KellyErrorEstimator::estimate function with an argument (usually defaulted, but explicitly given here) thatindicates the subdomain id of all those cells that we want to work on; note that this means that we also have to specify values for all those default arguments that lie before the one we want to give.

Secondly, we need to exchange the data. To do this, we could add up the refinement indicators from all processors, since they all only worked on a disjoint subset of the elements of the vector that holds these indicators. We could set up a distributed Trilinos vector for this, but that appears unnecessarily complicated because we would have to specify a partition of this vector, and none appears immediately obvious. Rather, we want to use the Trilinos communicator class to this for us, taking the local indicators as a collection of floating point values rather than a linear algebra vector. Unfortunately, the Trilinos communicator class doesn't appear to have function that wraps around the MPI add function; it has one that computes the maximum of a bunch of values, though, which in our case is equally good -- maybe even better, in case two processors should compute values for the same cell (which they shouldn't of course, unless we have made a mistake in specifying the arguments to the estimate function below). There is little snag again, however, that makes this a bit awkward: the Trilinos communicator class can take the maximum over all processors for each element of a vector, but only if the vector contains doubles. The vector returned by the KellyErrorEstimator::estimate function, on the other hand, has floats as its data type. An ugly, if workable way, is therefore to compute the indicators as floats, convert the vector to doubles, and form the maximum of that.

At the end of this chain of events, every processors has the complete set of refinement indicators, and the rest of the function proceeds as before.

 template <int dim>
 void BoussinesqFlowProblem<dim>::refine_mesh (const unsigned int max_grid_level)
 {
   computing_timer.enter_section ("Refine mesh structure, part 1");
   Vector<float> local_estimated_error_per_cell (triangulation.n_active_cells());
 
   KellyErrorEstimator<dim>::estimate (temperature_dof_handler,
                                       QGauss<dim-1>(temperature_degree+1),
                                       typename FunctionMap<dim>::type(),
                                       temperature_solution,
                                        local_estimated_error_per_cell,
                                       std::vector<bool>(),
                                       0,
                                       0,
                                       Utilities::System::get_this_mpi_process(MPI_COMM_WORLD));
 
   Vector<double> x_local_estimated_error_per_cell (triangulation.n_active_cells());
   x_local_estimated_error_per_cell = local_estimated_error_per_cell;
   Vector<double> estimated_error_per_cell (triangulation.n_active_cells());
 #ifdef DEAL_II_COMPILER_SUPPORTS_MPI
   MPI_Allreduce (&x_local_estimated_error_per_cell(0),
                  &estimated_error_per_cell(0),
                  triangulation.n_active_cells(), MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
 #else
   estimated_error_per_cell = x_local_estimated_error_per_cell;
 #endif
 
   GridRefinement::refine_and_coarsen_fixed_fraction (triangulation,
                                                      estimated_error_per_cell,
                                                      0.6, 0.2);
   if (triangulation.n_levels() > max_grid_level)
     for (typename Triangulation<dim>::active_cell_iterator
            cell = triangulation.begin_active(max_grid_level);
          cell != triangulation.end(); ++cell)
       cell->clear_refine_flag ();
 
   std::vector<TrilinosWrappers::Vector> x_temperature (2);
   x_temperature[0] = temperature_solution;
   x_temperature[1] = old_temperature_solution;
   TrilinosWrappers::BlockVector x_stokes = stokes_solution;
 
   SolutionTransfer<dim,TrilinosWrappers::Vector>
     temperature_trans(temperature_dof_handler);
   SolutionTransfer<dim,TrilinosWrappers::BlockVector>
     stokes_trans(stokes_dof_handler);
 
   triangulation.prepare_coarsening_and_refinement();
   temperature_trans.prepare_for_coarsening_and_refinement(x_temperature);
   stokes_trans.prepare_for_coarsening_and_refinement(x_stokes);
 
   triangulation.execute_coarsening_and_refinement ();
   computing_timer.exit_section();
 
   setup_dofs ();
 
   computing_timer.enter_section ("Refine mesh structure, part 2");
 
   std::vector<TrilinosWrappers::Vector> tmp (2);
   tmp[0].reinit (temperature_solution);
   tmp[1].reinit (temperature_solution);
   temperature_trans.interpolate(x_temperature, tmp);
 
   temperature_solution = tmp[0];
   old_temperature_solution = tmp[1];
   temperature_constraints.distribute(temperature_solution);
   temperature_constraints.distribute(old_temperature_solution);
 
   TrilinosWrappers::BlockVector x_stokes_new = stokes_solution;
   stokes_trans.interpolate (x_stokes, x_stokes_new);
   stokes_solution = x_stokes_new;
 
   rebuild_stokes_matrix              = true;
   rebuild_stokes_preconditioner      = true;
   rebuild_temperature_matrices       = true;
   rebuild_temperature_preconditioner = true;
 
   computing_timer.exit_section();
 }

BoussinesqFlowProblem::run

This is the final function in this class. It actually runs the program. It is, once more, very similar to step-31. The only thing that really changed is that we use the project_temperature_field() function instead of the library function VectorTools::project, the rest is as before.

 template <int dim>
 void BoussinesqFlowProblem<dim>::run ()
 {
   const unsigned int initial_refinement = (dim == 2 ? 5 : 2);
   const unsigned int n_pre_refinement_steps = (dim == 2 ? 2 : 2);
 
   GridGenerator::hyper_shell (triangulation,
                               Point<dim>(),
                               EquationData::R0,
                               EquationData::R1,
                               12,
                               true);
 
   static HyperShellBoundary<dim> boundary;
   triangulation.set_boundary (0, boundary);
   triangulation.set_boundary (1, boundary);
 
   global_Omega_diameter = GridTools::diameter (triangulation);
 
   triangulation.refine_global (initial_refinement);
 
   setup_dofs();
 
   unsigned int pre_refinement_step = 0;
 
   start_time_iteration:
 
   project_temperature_field ();
 
   timestep_number           = 0;
   time_step = old_time_step = 0;
 
   double time = 0;
 
   do
     {
       pcout << "Timestep " << timestep_number
             << ":  t=" << time/EquationData::year_in_seconds
             << " years"
             << std::endl;
 
       assemble_stokes_system ();
       build_stokes_preconditioner ();
       assemble_temperature_matrix ();
 
       solve ();
 
       pcout << std::endl;
 
       if ((timestep_number == 0) &&
           (pre_refinement_step < n_pre_refinement_steps))
         {
           refine_mesh (initial_refinement + n_pre_refinement_steps);
           ++pre_refinement_step;
           goto start_time_iteration;
         }
       else
         if ((timestep_number > 0) && (timestep_number % 10 == 0))
           refine_mesh (initial_refinement + n_pre_refinement_steps);
 
       output_results ();
 
       time += time_step;
       ++timestep_number;
 
       old_stokes_solution          = stokes_solution;
       old_old_temperature_solution = old_temperature_solution;
       old_temperature_solution     = temperature_solution;
     }
   while (time <= EquationData::end_time);
 }

The main function

This is copied verbatim from step-31:

 int main (int argc, char *argv[])
 {
   try
     {
       deallog.depth_console (0);
 
       Utilities::System::MPI_InitFinalize mpi_initialization(argc, argv);
 
       BoussinesqFlowProblem<2> flow_problem;
       flow_problem.run ();
     }
   catch (std::exception &exc)
     {
       std::cerr << std::endl << std::endl
                 << "----------------------------------------------------"
                 << std::endl;
       std::cerr << "Exception on processing: " << std::endl
                 << exc.what() << std::endl
                 << "Aborting!" << std::endl
                 << "----------------------------------------------------"
                 << std::endl;
 
       return 1;
     }
   catch (...)
     {
       std::cerr << std::endl << std::endl
                 << "----------------------------------------------------"
                 << std::endl;
       std::cerr << "Unknown exception!" << std::endl
                 << "Aborting!" << std::endl
                 << "----------------------------------------------------"
                 << std::endl;
       return 1;
     }
 
   return 0;
 }

Results

When run, the program simulates convection in 3d in much the same way as step-31 did, though with an entirely different testcase.

Comparison of results with step-31

Before we go to this testcase, however, let us show a few results from a slightly earlier version of this program that was solving exactly the testcase we used in step-31, just that we now solve it in parallel and with much higher resolution. We show these results mainly for comparison.

Here are two images that show this higher resolution if we choose a 3d computation in main() and if we set initial_refinement=3 and n_pre_refinement_steps=4. At the time steps shown, the meshes had around 72,000 and 236,000 cells, for a total of 2,680,000 and 8,250,000 degrees of freedom, respectively, more than an order of magnitude more than we had available in step-31:

step-32.3d.cube.0.png
step-32.3d.cube.1.png

The computation was done on a subset of 50 processors of the Brazos cluster at Texas A&M University.

Results for a 2d circular shell testcase

If we run the program as shown above, the output will look roughly like this, producing the final part of the output after some 2 days when run on 10 processors:

Number of active cells: 12288 (on 6 levels)
Number of degrees of freedom: 162432 (99840+12672+49920)

Timestep 0:  t=0 years
   Assembling...
   Rebuilding Stokes preconditioner...
   Solving...
   14 BiCGStab iterations for Stokes subsystem.
   Maximal velocity: 1.80187 cm/year
   Time step: 607937 years
   15 CG iterations for temperature
   Temperature range: 973 4273.27

Number of active cells: 17988 (on 7 levels)
Number of degrees of freedom: 240504 (147888+18672+73944)

Timestep 0:  t=0 years
   Assembling...
   Rebuilding Stokes preconditioner...
   Solving...
   3 BiCGStab iterations for Stokes subsystem.
   Maximal velocity: 1.52763 cm/year
   Time step: 357565 years
   17 CG iterations for temperature
   Temperature range: 973 4274.08

Number of active cells: 34716 (on 8 levels)
Number of degrees of freedom: 474550 (291896+36706+145948)

Timestep 0:  t=0 years
   Assembling...
   Rebuilding Stokes preconditioner...
   Solving...
   3 BiCGStab iterations for Stokes subsystem.
   Maximal velocity: 2.97931 cm/year
   Time step: 91574.3 years
   17 CG iterations for temperature
   Temperature range: 973 4277.31

Timestep 1:  t=91574.3 years
   Assembling...
   Solving...
   2 BiCGStab iterations for Stokes subsystem.
   Maximal velocity: 2.86431 cm/year
   Time step: 95251.1 years
   18 CG iterations for temperature
   Temperature range: 973 4277.31

[...]

Timestep 53906:  t=9.99985e+08 years
   Assembling...
   Solving...
   4 BiCGStab iterations for Stokes subsystem.
   Maximal velocity: 12.2898 cm/year
   Time step: 22199.6 years
   18 CG iterations for temperature
   Temperature range: 973 4273.07

+---------------------------------------------+------------+------------+
| Total wallclock time elapsed since start    |  1.87e+05s |            |
|                                             |            |            |
| Section                         | no. calls |  wall time | % of total |
+---------------------------------+-----------+------------+------------+
| Assemble Stokes system          |     53909 |  1.06e+04s |       5.7% |
| Assemble temperature matrices   |      5393 |       913s |      0.49% |
| Assemble temperature rhs        |     53909 |  1.24e+04s |       6.6% |
| Build Stokes preconditioner     |      5393 |  4.52e+03s |       2.4% |
| Solve Stokes system             |     53909 |  6.94e+04s |        37% |
| Solve temperature system        |     53909 |  2.34e+04s |        13% |
| Postprocessing                  |      2157 |  5.51e+03s |       2.9% |
| Refine mesh structure, part 1   |      5392 |  4.71e+04s |        25% |
| Refine mesh structure, part 2   |      5392 |  1.02e+03s |      0.55% |
| Setup dof systems               |      5393 |  1.11e+04s |       5.9% |
+---------------------------------+-----------+------------+------------+

As can be seen here, we spend most of the compute time in assembling linear systems, refining the mesh, and in particular in solving the Stokes and temperature linear systems.

The 50% spent on solving the linear systems are affected in large part because the Brazos cluster has a relatively slow ethernet interconnect. A cluster with a faster interconnect, for example using infiniband, should do better in this regard. On the other hand, there is little hope to do better with assembling the linear systems, though one could do significantly better with estimating the error by making sure that each processor only estimates the error on those cells it owns.

The program writes output every 25th time step, but we won't show all 2100 or so images this produces. Rather, let us only show the output from every 2500th time step here, even though this does, of course, not do full justice to the dynamics that are going on:
step-32.2d.temperature.0000.png

step-32.2d.temperature.0100.png

step-32.2d.temperature.0200.png

step-32.2d.temperature.0300.png

step-32.2d.temperature.0400.png

step-32.2d.temperature.0500.png

step-32.2d.temperature.0600.png

step-32.2d.temperature.0700.png

step-32.2d.temperature.0800.png

step-32.2d.temperature.0900.png

step-32.2d.temperature.1000.png

step-32.2d.temperature.1100.png

step-32.2d.temperature.1200.png

step-32.2d.temperature.1300.png

step-32.2d.temperature.1400.png

step-32.2d.temperature.1500.png

step-32.2d.temperature.1600.png

step-32.2d.temperature.1700.png

step-32.2d.temperature.1800.png

step-32.2d.temperature.1900.png

step-32.2d.temperature.2000.png

step-32.2d.temperature.2100.png

step-32.2d.grid.2100.png

step-32.2d.partition.2100.png

The last two images show the grid as well as the partitioning of the mesh for the last timestep shown into the 10 subdomains used for this computation. The full dynamics are really only visible by looking at an animation. <a href="http://www.math.tamu.edu/~bangerth/images/pictures/convection-outward/step-32.2d.convection.gif">At this site is such an animation. Beware that this animation is about 20MB large, though it is well worth watching due to its almost artistic quality.

If you watch the movie, you'll see that the convection pattern goes through several stages: First, it gets rid of the instable temperature layering with the hot material overlaid by the dense cold material. After this great driver is removed and we have a sort of stable situation, a few blobs start to separate from the hot boundary layer at the inner ring and rise up, with a few cold fingers also dropping down from the outer ring. During this phase, the solution remains mostly symmetric, reflecting the 12-fold symmetry of the original mesh. In a final phase, the fluid enters vigorous chaotic stirring in which all symmetries are lost. This is a pattern that then continues to dominate flow.

These different phases can also be identified if we look at the maximal velocity as a function of time in the simulation:

step-32.2d.t_vs_vmax.png

Here, the velocity (shown in centimeters per year) becomes very large, to the order of several meters per year) at the beginning when the temperature layering is instable. It then calms down to relatively small values before picking up again in the chaotic stirring regime. There, it remains in the range of 10-40 centimeters per year, quite within the physically expected region.

Possibilities for extensions

Apart from the various possibilities for extensions already outlined in the step-31, here are a few more ideas:

The plain program

(If you are looking at a locally installed deal.II version, then the program can be found at /u /bangerth /tmp /homepage /deal.II /examples /step-32 /step-32.cc . Otherwise, this is only the path on some remote server.)

 / * Subversion Id:  @ref step_32 "step-32".cc 20544 2010-02-10 20:51:42Z bangerth  * /
 / * Author: Martin Kronbichler, Uppsala University,
            Wolfgang Bangerth, Texas A&M University 2008, 2009 * /
 / *                                                                * /
 / *    Copyright (C) 2008, 2009, 2010 by the deal.II authors * /
 / *                                                                * /
 / *    This file is subject to QPL and may not be  distributed     * /
 / *    without copyright and license information. Please refer     * /
 / *    to the file deal.II/doc/license.html for the  text  and     * /
 / *    further information on this license.                        * /
 #include <base/quadrature_lib.h>
 #include <base/logstream.h>
 #include <base/function.h>
 #include <base/utilities.h>
 #include <base/conditional_ostream.h>
 #include <base/work_stream.h>
 #include <base/timer.h>
 
 #include <lac/full_matrix.h>
 #include <lac/solver_bicgstab.h>
 #include <lac/solver_cg.h>
 #include <lac/constraint_matrix.h>
 #include <lac/block_sparsity_pattern.h>
 #include <lac/trilinos_block_vector.h>
 #include <lac/trilinos_sparse_matrix.h>
 #include <lac/trilinos_block_sparse_matrix.h>
 #include <lac/trilinos_precondition.h>
 
 #include <grid/tria.h>
 #include <grid/grid_generator.h>
 #include <grid/tria_accessor.h>
 #include <grid/tria_iterator.h>
 #include <grid/filtered_iterator.h>
 #include <grid/tria_boundary_lib.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 <fe/fe_q.h>
 #include <fe/fe_dgq.h>
 #include <fe/fe_system.h>
 #include <fe/fe_values.h>
 #include <fe/mapping_q1.h>
 
 #include <numerics/vectors.h>
 #include <numerics/matrices.h>
 #include <numerics/data_out.h>
 #include <numerics/error_estimator.h>
 #include <numerics/solution_transfer.h>
 
 #include <fstream>
 #include <iostream>
 #include <sstream>
 
 #include <base/index_set.h>
 
 
 using namespace dealii;
 namespace EquationData
 {
   const double eta                   = 1e21;
   const double kappa                 = 1e-6;
   const double reference_density     = 3300;
   const double reference_temperature = 293;
   const double expansion_coefficient = 2e-5;
 
   const double R0      = 6371000.-2890000.;
   const double R1      = 6371000.-  35000.;
 
   const double T0      = 4000+273;
   const double T1      =  700+273;
 
   const double year_in_seconds  = 60*60*24*365.2425;
   const double end_time         = 1e9 * year_in_seconds;
 
   const double pressure_scaling = eta / (R1-R0);
 
 
   double density (const double temperature)
   {
     return (reference_density *
            (1 - expansion_coefficient * (temperature -
                                          reference_temperature)));
   }
 
 
   template <int dim>
   Tensor<1,dim> gravity_vector (const Point<dim> &p)
   {
     return -9.81 * p / R1;
   }
 
 
   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>  &p,
                                        const unsigned int) const
   {
     const double r = p.norm();
     const double h = R1-R0;
 
     const double s = (r-R0)/h;
 
     return T1+(T0-T1)*((1-s)*(1-s));
   }
 
 
   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
   {
     return 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);
   }
 }
 namespace LinearSolvers
 {
   template <class PreconditionerA, class PreconditionerMp>
   class BlockSchurPreconditioner : public Subscriptor
   {
     public:
       BlockSchurPreconditioner (
        const TrilinosWrappers::BlockSparseMatrix  &S,
        const PreconditionerMp                     &Mppreconditioner,
        const PreconditionerA                      &Apreconditioner);
 
       void vmult (TrilinosWrappers::MPI::BlockVector       &dst,
                  const TrilinosWrappers::MPI::BlockVector &src) const;
 
     private:
       const SmartPointer<const TrilinosWrappers::BlockSparseMatrix> stokes_matrix;
       const PreconditionerMp &mp_preconditioner;
       const PreconditionerA  &a_preconditioner;
       mutable TrilinosWrappers::MPI::Vector tmp;
   };
 
 
 
   template <class PreconditionerA, class PreconditionerMp>
   BlockSchurPreconditioner<PreconditionerA, PreconditionerMp>::
   BlockSchurPreconditioner(const TrilinosWrappers::BlockSparseMatrix &S,
                           const PreconditionerMp                    &Mppreconditioner,
                           const PreconditionerA                     &Apreconditioner)
                  :
                  stokes_matrix     (&S),
                  mp_preconditioner (Mppreconditioner),
                  a_preconditioner  (Apreconditioner),
                  tmp               (stokes_matrix->block(1,1).row_partitioner())
   {}
 
 
 
   template <class PreconditionerA, class PreconditionerMp>
   void BlockSchurPreconditioner<PreconditionerA, PreconditionerMp>::vmult (
     TrilinosWrappers::MPI::BlockVector       &dst,
     const TrilinosWrappers::MPI::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;
     mp_preconditioner.vmult (dst.block(1), tmp);
   }
 }
 namespace Assembly
 {
   namespace Scratch
   {
     template <int dim>
     struct StokesPreconditioner
     {
        StokesPreconditioner (const FiniteElement<dim> &stokes_fe,
                              const Quadrature<dim>    &stokes_quadrature,
                              const UpdateFlags         update_flags);
        StokesPreconditioner (const StokesPreconditioner &data);
 
        FEValues<dim>               stokes_fe_values;
 
        std::vector<Tensor<2,dim> > grad_phi_u;
        std::vector<double>         phi_p;
     };
 
     template <int dim>
     StokesPreconditioner<dim>::
     StokesPreconditioner (const FiniteElement<dim> &stokes_fe,
                          const Quadrature<dim>    &stokes_quadrature,
                          const UpdateFlags         update_flags)
                    :
                    stokes_fe_values (stokes_fe, stokes_quadrature,
                                      update_flags),
                    grad_phi_u (stokes_fe.dofs_per_cell),
                    phi_p (stokes_fe.dofs_per_cell)
     {}
 
 
 
     template <int dim>
     StokesPreconditioner<dim>::
     StokesPreconditioner (const StokesPreconditioner &scratch)
                    :
                    stokes_fe_values (scratch.stokes_fe_values.get_fe(),
                                      scratch.stokes_fe_values.get_quadrature(),
                                      scratch.stokes_fe_values.get_update_flags()),
                    grad_phi_u (scratch.grad_phi_u),
                    phi_p (scratch.phi_p)
     {}
 
 
 
     template <int dim>
     struct StokesSystem : public StokesPreconditioner<dim>
     {
        StokesSystem (const FiniteElement<dim> &stokes_fe,
                      const Quadrature<dim>    &stokes_quadrature,
                      const UpdateFlags         stokes_update_flags,
                      const FiniteElement<dim> &temperature_fe,
                      const UpdateFlags         temperature_update_flags);
 
        StokesSystem (const StokesSystem<dim> &data);
 
        FEValues<dim>  temperature_fe_values;
 
        std::vector<Tensor<1,dim> >          phi_u;
        std::vector<SymmetricTensor<2,dim> > grads_phi_u;
        std::vector<double>                  div_phi_u;
 
        std::vector<double>                  old_temperature_values;
     };
 
 
     template <int dim>
     StokesSystem<dim>::
     StokesSystem (const FiniteElement<dim> &stokes_fe,
                  const Quadrature<dim>    &stokes_quadrature,
                  const UpdateFlags         stokes_update_flags,
                  const FiniteElement<dim> &temperature_fe,
                  const UpdateFlags         temperature_update_flags)
                    :
                    StokesPreconditioner<dim> (stokes_fe, stokes_quadrature,
                                               stokes_update_flags),
                    temperature_fe_values (temperature_fe, stokes_quadrature,
                                           temperature_update_flags),
                    phi_u (stokes_fe.dofs_per_cell),
                    grads_phi_u (stokes_fe.dofs_per_cell),
                    div_phi_u (stokes_fe.dofs_per_cell),
                    old_temperature_values (stokes_quadrature.n_quadrature_points)
     {}
 
 
     template <int dim>
     StokesSystem<dim>::
     StokesSystem (const StokesSystem<dim> &scratch)
                    :
                    StokesPreconditioner<dim> (scratch),
                    temperature_fe_values (scratch.temperature_fe_values.get_fe(),
                                           scratch.temperature_fe_values.get_quadrature(),
                                           scratch.temperature_fe_values.get_update_flags()),
                    phi_u (scratch.phi_u),
                    grads_phi_u (scratch.grads_phi_u),
                    div_phi_u (scratch.div_phi_u),
                    old_temperature_values (scratch.old_temperature_values)
     {}
 
 
 
     template <int dim>
     struct TemperatureMatrix
     {
        TemperatureMatrix (const FiniteElement<dim> &temperature_fe,
                           const Quadrature<dim>    &temperature_quadrature);
        TemperatureMatrix (const TemperatureMatrix &data);
 
        FEValues<dim>               temperature_fe_values;
 
        std::vector<double>         phi_T;
        std::vector<Tensor<1,dim> > grad_phi_T;
     };
 
     template <int dim>
     TemperatureMatrix<dim>::
     TemperatureMatrix (const FiniteElement<dim> &temperature_fe,
                       const Quadrature<dim>    &temperature_quadrature)
                    :
                    temperature_fe_values (temperature_fe, temperature_quadrature,
                                           update_values    | update_gradients |
                                           update_JxW_values),
                    phi_T (temperature_fe.dofs_per_cell),
                    grad_phi_T (temperature_fe.dofs_per_cell)
     {}
 
 
     template <int dim>
     TemperatureMatrix<dim>::
     TemperatureMatrix (const TemperatureMatrix &scratch)
                    :
                    temperature_fe_values (scratch.temperature_fe_values.get_fe(),
                                           scratch.temperature_fe_values.get_quadrature(),
                                           scratch.temperature_fe_values.get_update_flags()),
                    phi_T (scratch.phi_T),
                    grad_phi_T (scratch.grad_phi_T)
     {}
 
 
     template <int dim>
     struct TemperatureRHS
     {
        TemperatureRHS (const FiniteElement<dim> &temperature_fe,
                        const FiniteElement<dim> &stokes_fe,
                        const Quadrature<dim>    &quadrature);
        TemperatureRHS (const TemperatureRHS &data);
 
        FEValues<dim>               temperature_fe_values;
        FEValues<dim>               stokes_fe_values;
 
        std::vector<double>         phi_T;
        std::vector<Tensor<1,dim> > grad_phi_T;
 
        std::vector<Tensor<1,dim> > old_velocity_values;
        std::vector<Tensor<1,dim> > old_old_velocity_values;
 
        std::vector<double>         old_temperature_values;
        std::vector<double>         old_old_temperature_values;
        std::vector<Tensor<1,dim> > old_temperature_grads;
        std::vector<Tensor<1,dim> > old_old_temperature_grads;
        std::vector<double>         old_temperature_laplacians;
        std::vector<double>         old_old_temperature_laplacians;
 
        std::vector<double>         gamma_values;
     };
 
     template <int dim>
     TemperatureRHS<dim>::
     TemperatureRHS (const FiniteElement<dim> &temperature_fe,
                    const FiniteElement<dim> &stokes_fe,
                    const Quadrature<dim>    &quadrature)
                    :
                    temperature_fe_values (temperature_fe, quadrature,
                                           update_values    |
                                           update_gradients |
                                           update_hessians  |
                                           update_quadrature_points |
                                           update_JxW_values),
                    stokes_fe_values (stokes_fe, quadrature,
                                      update_values),
                    phi_T (temperature_fe.dofs_per_cell),
                    grad_phi_T (temperature_fe.dofs_per_cell),
 
                    old_velocity_values (quadrature.n_quadrature_points),
                    old_old_velocity_values (quadrature.n_quadrature_points),
 
                    old_temperature_values (quadrature.n_quadrature_points),
                    old_old_temperature_values(quadrature.n_quadrature_points),
                    old_temperature_grads(quadrature.n_quadrature_points),
                    old_old_temperature_grads(quadrature.n_quadrature_points),
                    old_temperature_laplacians(quadrature.n_quadrature_points),
                    old_old_temperature_laplacians(quadrature.n_quadrature_points),
 
                    gamma_values (quadrature.n_quadrature_points)
     {}
 
 
     template <int dim>
     TemperatureRHS<dim>::
     TemperatureRHS (const TemperatureRHS &scratch)
                    :
                    temperature_fe_values (scratch.temperature_fe_values.get_fe(),
                                           scratch.temperature_fe_values.get_quadrature(),
                                           scratch.temperature_fe_values.get_update_flags()),
                    stokes_fe_values (scratch.stokes_fe_values.get_fe(),
                                      scratch.stokes_fe_values.get_quadrature(),
                                      scratch.stokes_fe_values.get_update_flags()),
                    phi_T (scratch.phi_T),
                    grad_phi_T (scratch.grad_phi_T),
 
                    old_velocity_values (scratch.old_velocity_values),
                    old_old_velocity_values (scratch.old_old_velocity_values),
 
                    old_temperature_values (scratch.old_temperature_values),
                    old_old_temperature_values (scratch.old_old_temperature_values),
                    old_temperature_grads (scratch.old_temperature_grads),
                    old_old_temperature_grads (scratch.old_old_temperature_grads),
                    old_temperature_laplacians (scratch.old_temperature_laplacians),
                    old_old_temperature_laplacians (scratch.old_old_temperature_laplacians),
 
                    gamma_values (scratch.gamma_values)
     {}
   }
 
 
   namespace CopyData
   {
     template <int dim>
     struct StokesPreconditioner
     {
        StokesPreconditioner (const FiniteElement<dim> &stokes_fe);
        StokesPreconditioner (const StokesPreconditioner &data);
 
        FullMatrix<double>          local_matrix;
        std::vector<unsigned int>   local_dof_indices;
     };
 
     template <int dim>
     StokesPreconditioner<dim>::
     StokesPreconditioner (const FiniteElement<dim> &stokes_fe)
                    :
                    local_matrix (stokes_fe.dofs_per_cell,
                                  stokes_fe.dofs_per_cell),
                    local_dof_indices (stokes_fe.dofs_per_cell)
     {}
 
 
 
     template <int dim>
     StokesPreconditioner<dim>::
     StokesPreconditioner (const StokesPreconditioner &data)
                    :
                    local_matrix (data.local_matrix),
                    local_dof_indices (data.local_dof_indices)
     {}
 
 
 
     template <int dim>
     struct StokesSystem : public StokesPreconditioner<dim>
     {
        StokesSystem (const FiniteElement<dim> &stokes_fe);
        StokesSystem (const StokesSystem<dim> &data);
 
        Vector<double> local_rhs;
     };
 
 
     template <int dim>
     StokesSystem<dim>::
     StokesSystem (const FiniteElement<dim> &stokes_fe)
                    :
                    StokesPreconditioner<dim> (stokes_fe),
                    local_rhs (stokes_fe.dofs_per_cell)
     {}
 
 
     template <int dim>
     StokesSystem<dim>::
     StokesSystem (const StokesSystem<dim> &data)
                    :
                    StokesPreconditioner<dim> (data),
                    local_rhs (data.local_rhs)
     {}
 
 
 
     template <int dim>
     struct TemperatureMatrix
     {
        TemperatureMatrix (const FiniteElement<dim> &temperature_fe);
        TemperatureMatrix (const TemperatureMatrix &data);
 
        FullMatrix<double>          local_mass_matrix;
        FullMatrix<double>          local_stiffness_matrix;
        std::vector<unsigned int>   local_dof_indices;
     };
 
     template <int dim>
     TemperatureMatrix<dim>::
     TemperatureMatrix (const FiniteElement<dim> &temperature_fe)
                    :
                    local_mass_matrix (temperature_fe.dofs_per_cell,
                                       temperature_fe.dofs_per_cell),
                    local_stiffness_matrix (temperature_fe.dofs_per_cell,
                                            temperature_fe.dofs_per_cell),
                    local_dof_indices (temperature_fe.dofs_per_cell)
     {}
 
 
     template <int dim>
     TemperatureMatrix<dim>::
     TemperatureMatrix (const TemperatureMatrix &data)
                    :
                    local_mass_matrix (data.local_mass_matrix),
                    local_stiffness_matrix (data.local_stiffness_matrix),
                    local_dof_indices (data.local_dof_indices)
     {}
 
 
     template <int dim>
     struct TemperatureRHS
     {
        TemperatureRHS (const FiniteElement<dim> &temperature_fe);
        TemperatureRHS (const TemperatureRHS &data);
 
        Vector<double>              local_rhs;
        std::vector<unsigned int>   local_dof_indices;
         FullMatrix<double>          matrix_for_bc;
     };
 
     template <int dim>
     TemperatureRHS<dim>::
     TemperatureRHS (const FiniteElement<dim> &temperature_fe)
                    :
                    local_rhs (temperature_fe.dofs_per_cell),
                    local_dof_indices (temperature_fe.dofs_per_cell),
                    matrix_for_bc (temperature_fe.dofs_per_cell,
                                   temperature_fe.dofs_per_cell)
     {}
 
 
     template <int dim>
     TemperatureRHS<dim>::
     TemperatureRHS (const TemperatureRHS &data)
                    :
                    local_rhs (data.local_rhs),
                    local_dof_indices (data.local_dof_indices),
                    matrix_for_bc (data.matrix_for_bc)
     {}
   }
 }
 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_matrix ();
     void assemble_temperature_system (const double maximal_velocity);
     void project_temperature_field ();
     double get_maximal_velocity () const;
     std::pair<double,double> get_extrapolated_temperature_range () const;
     void solve ();
     void output_results ();
     void refine_mesh (const unsigned int max_grid_level);
 
     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<double>          &old_temperature_laplacians,
                      const std::vector<double>          &old_old_temperature_laplacians,
                      const std::vector<Tensor<1,dim> >  &old_velocity_values,
                      const std::vector<Tensor<1,dim> >  &old_old_velocity_values,
                      const std::vector<double>          &gamma_values,
                      const double                        global_u_infty,
                      const double                        global_T_variation,
                      const double                        cell_diameter) const;
 
 
     ConditionalOStream                  pcout;
 
     Triangulation<dim>                  triangulation;
     double                              global_Omega_diameter;
 
     const unsigned int                  stokes_degree;
     FESystem<dim>                       stokes_fe;
     DoFHandler<dim>                     stokes_dof_handler;
     ConstraintMatrix                    stokes_constraints;
 
     TrilinosWrappers::BlockSparseMatrix stokes_matrix;
     TrilinosWrappers::BlockSparseMatrix stokes_preconditioner_matrix;
 
     TrilinosWrappers::BlockVector       stokes_solution;
     TrilinosWrappers::BlockVector       old_stokes_solution;
     TrilinosWrappers::MPI::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::MPI::Vector       temperature_rhs;
 
 
     double time_step;
     double old_time_step;
     unsigned int timestep_number;
 
     std_cxx1x::shared_ptr<TrilinosWrappers::PreconditionAMG> Amg_preconditioner;
     std_cxx1x::shared_ptr<TrilinosWrappers::PreconditionILU> Mp_preconditioner;
     std_cxx1x::shared_ptr<TrilinosWrappers::PreconditionIC>  T_preconditioner;
 
     bool rebuild_stokes_matrix;
     bool rebuild_stokes_preconditioner;
     bool rebuild_temperature_matrices;
     bool rebuild_temperature_preconditioner;
 
     TimerOutput computing_timer;
 
     void setup_stokes_matrix (const std::vector<IndexSet> &stokes_partitioning);
     void setup_stokes_preconditioner (const std::vector<IndexSet> &stokes_partitioning);
     void setup_temperature_matrices (const IndexSet &temperature_partitioning);
 
     void
     local_assemble_stokes_preconditioner (const typename DoFHandler<dim>::active_cell_iterator &cell,
                                          Assembly::Scratch::StokesPreconditioner<dim> &scratch,
                                          Assembly::CopyData::StokesPreconditioner<dim> &data);
 
     void
     copy_local_to_global_stokes_preconditioner (const Assembly::CopyData::StokesPreconditioner<dim> &data);
 
 
     void
     local_assemble_stokes_system (const typename DoFHandler<dim>::active_cell_iterator &cell,
                                  Assembly::Scratch::StokesSystem<dim>  &scratch,
                                  Assembly::CopyData::StokesSystem<dim> &data);
 
     void
     copy_local_to_global_stokes_system (const Assembly::CopyData::StokesSystem<dim> &data);
 
 
     void
     local_assemble_temperature_matrix (const typename DoFHandler<dim>::active_cell_iterator &cell,
                                       Assembly::Scratch::TemperatureMatrix<dim>  &scratch,
                                       Assembly::CopyData::TemperatureMatrix<dim> &data);
 
     void
     copy_local_to_global_temperature_matrix (const Assembly::CopyData::TemperatureMatrix<dim> &data);
 
 
 
     void
     local_assemble_temperature_rhs (const std::pair<double,double> global_T_range,
                                    const double                   global_max_velocity,
                                    const typename DoFHandler<dim>::active_cell_iterator &cell,
                                    Assembly::Scratch::TemperatureRHS<dim> &scratch,
                                    Assembly::CopyData::TemperatureRHS<dim> &data);
 
     void
     copy_local_to_global_temperature_rhs (const Assembly::CopyData::TemperatureRHS<dim> &data);
 };
 template <int dim>
 BoussinesqFlowProblem<dim>::BoussinesqFlowProblem ()
                 :
                pcout (std::cout,
                       (Utilities::System::
                        get_this_mpi_process(MPI_COMM_WORLD)
                        == 0)),
 
                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_stokes_preconditioner (true),
                rebuild_temperature_matrices (true),
                rebuild_temperature_preconditioner (true),
 
                computing_timer (pcout, TimerOutput::summary,
                                 TimerOutput::wall_times)
 {}
 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<Tensor<1,dim> > velocity_values(n_q_points);
 
   const FEValuesExtractors::Vector velocities (0);
 
   double max_local_velocity = 0;
 
   typename DoFHandler<dim>::active_cell_iterator
     cell = stokes_dof_handler.begin_active(),
     endc = stokes_dof_handler.end();
   for (; cell!=endc; ++cell)
     if (cell->subdomain_id() ==
        Utilities::System::get_this_mpi_process(MPI_COMM_WORLD))
       {
        fe_values.reinit (cell);
        fe_values[velocities].get_function_values (stokes_solution,
                                                   velocity_values);
 
        for (unsigned int q=0; q<n_q_points; ++q)
          max_local_velocity = std::max (max_local_velocity,
                                         velocity_values[q].norm());
       }
 
   double max_velocity = 0.;
 #ifdef DEAL_II_COMPILER_SUPPORTS_MPI
   MPI_Allreduce (&max_local_velocity, &max_velocity, 1, MPI_DOUBLE,
                 MPI_MAX, MPI_COMM_WORLD);
 #else
   max_velocity = max_local_velocity;
 #endif
 
   return max_velocity;
 }
 
 
 
 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_local_temperature = (1. + time_step/old_time_step) *
                                     old_temperature_solution.linfty_norm()
                                     +
                                     time_step/old_time_step *
                                     old_old_temperature_solution.linfty_norm(),
             max_local_temperature = -min_local_temperature;
 
       typename DoFHandler<dim>::active_cell_iterator
        cell = temperature_dof_handler.begin_active(),
        endc = temperature_dof_handler.end();
       for (; cell!=endc; ++cell)
        if (cell->subdomain_id() ==
            Utilities::System::get_this_mpi_process(MPI_COMM_WORLD))
          {
            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_local_temperature = std::min (min_local_temperature,
                                                  temperature);
                max_local_temperature = std::max (max_local_temperature,
                                                  temperature);
              }
          }
 
       double min_temperature, max_temperature;
 #ifdef DEAL_II_COMPILER_SUPPORTS_MPI
       MPI_Allreduce (&max_local_temperature, &max_temperature, 1, MPI_DOUBLE,
                     MPI_MAX, MPI_COMM_WORLD);
       MPI_Allreduce (&min_local_temperature, &min_temperature, 1, MPI_DOUBLE,
                     MPI_MIN, MPI_COMM_WORLD);
 #else
       min_temperature = min_local_temperature;
       max_temperature = max_local_temperature;
 #endif
 
       return std::make_pair(min_temperature, max_temperature);
     }
   else
     {
       double min_local_temperature = old_temperature_solution.linfty_norm(),
             max_local_temperature = -min_local_temperature;
 
       typename DoFHandler<dim>::active_cell_iterator
        cell = temperature_dof_handler.begin_active(),
        endc = temperature_dof_handler.end();
       for (; cell!=endc; ++cell)
        if (cell->subdomain_id() ==
            Utilities::System::get_this_mpi_process(MPI_COMM_WORLD))
          {
            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_local_temperature = std::min (min_local_temperature,
                                                  temperature);
                max_local_temperature = std::max (max_local_temperature,
                                                  temperature);
              }
          }
 
       double min_temperature, max_temperature;
 #ifdef DEAL_II_COMPILER_SUPPORTS_MPI
       MPI_Allreduce (&max_local_temperature, &max_temperature, 1, MPI_DOUBLE,
                     MPI_MAX, MPI_COMM_WORLD);
       MPI_Allreduce (&min_local_temperature, &min_temperature, 1, MPI_DOUBLE,
                     MPI_MIN, MPI_COMM_WORLD);
 #else
       min_temperature = min_local_temperature;
       max_temperature = max_local_temperature;
 #endif
 
       return std::make_pair(min_temperature, max_temperature);
     }
 }
 
 
 
 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<double>          &old_temperature_laplacians,
                   const std::vector<double>          &old_old_temperature_laplacians,
                   const std::vector<Tensor<1,dim> >  &old_velocity_values,
                   const std::vector<Tensor<1,dim> >  &old_old_velocity_values,
                   const std::vector<double>          &gamma_values,
                   const double                        global_u_infty,
                   const double                        global_T_variation,
                   const double                        cell_diameter) const
 {
   const double beta = 0.026 * 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)
     {
       const Tensor<1,dim> u = (old_velocity_values[q] +
                               old_old_velocity_values[q]) / 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
                                   * (old_temperature_laplacians[q] +
                                      old_old_temperature_laplacians[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);
     }
 
   if (timestep_number == 0)
     return beta * max_velocity * cell_diameter;
   else
     {
       Assert (old_time_step > 0, ExcInternalError());
 
       const double c_R = 0.11;
       const double global_scaling = c_R * 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));
     }
 }
 
 
 
 template <int dim>
 void BoussinesqFlowProblem<dim>::project_temperature_field ()
 {
   assemble_temperature_matrix ();
 
   QGauss<dim> quadrature(temperature_degree+2);
   UpdateFlags update_flags = UpdateFlags(update_values   |
                                         update_quadrature_points |
                                         update_JxW_values);
   FEValues<dim> fe_values (temperature_fe, quadrature, update_flags);
 
   const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
                     n_q_points    = fe_values.n_quadrature_points;
 
   std::vector<unsigned int> local_dof_indices (dofs_per_cell);
   Vector<double> cell_vector (dofs_per_cell);
   FullMatrix<double> matrix_for_bc (dofs_per_cell, dofs_per_cell);
 
   typename DoFHandler<dim>::active_cell_iterator
     cell = temperature_dof_handler.begin_active(),
     endc = temperature_dof_handler.end();
 
   std::vector<double> rhs_values(n_q_points);
 
   TrilinosWrappers::MPI::Vector
     rhs (temperature_mass_matrix.row_partitioner()),
     solution (temperature_mass_matrix.row_partitioner());
 
   for (; cell!=endc; ++cell)
     if (cell->subdomain_id() ==
        Utilities::System::get_this_mpi_process(MPI_COMM_WORLD))
       {
        cell->get_dof_indices (local_dof_indices);
        fe_values.reinit(cell);
 
        EquationData::TemperatureInitialValues<dim>().value_list
          (fe_values.get_quadrature_points(), rhs_values);
 
        cell_vector = 0;
        matrix_for_bc = 0;
        for (unsigned int point=0; point<n_q_points; ++point)
          for (unsigned int i=0; i<dofs_per_cell; ++i)
            {
              cell_vector(i) += rhs_values[point] *
                                fe_values.shape_value(i,point) *
                                fe_values.JxW(point);
              if (temperature_constraints.is_inhomogeneously_constrained(local_dof_indices[i]))
                {
                  for (unsigned int j=0; j<dofs_per_cell; ++j)
                    matrix_for_bc(j,i) += fe_values.shape_value(i,point) *
                                          fe_values.shape_value(j,point) *
                                          fe_values.JxW(point);
                }
            }
 
        temperature_constraints.distribute_local_to_global (cell_vector,
                                                            local_dof_indices,
                                                            rhs,
                                                            matrix_for_bc);
       }
 
   rhs.compress ();
 
   SolverControl solver_control(5*rhs.size(), 1e-12*rhs.l2_norm());
   SolverCG<TrilinosWrappers::MPI::Vector> cg(solver_control);
 
   TrilinosWrappers::PreconditionSSOR preconditioner_mass;
   preconditioner_mass.initialize(temperature_mass_matrix, 1.3);
 
   cg.solve (temperature_mass_matrix, solution, rhs, preconditioner_mass);
 
   old_temperature_solution = solution;
   temperature_constraints.distribute (old_temperature_solution);
 }
 template <int dim>
 void BoussinesqFlowProblem<dim>::
   setup_stokes_matrix (const std::vector<IndexSet> &stokes_partitioning)
 {
   stokes_matrix.clear ();
 
   TrilinosWrappers::BlockSparsityPattern sp (stokes_partitioning,
                                             MPI_COMM_WORLD);
 
   Table<2,DoFTools::Coupling> coupling (dim+1, dim+1);
 
   for (unsigned int c=0; c<dim+1; ++c)
     for (unsigned int d=0; d<dim+1; ++d)
       if (! ((c==dim) && (d==dim)))
        coupling[c][d] = DoFTools::always;
       else
        coupling[c][d] = DoFTools::none;
 
   DoFTools::make_sparsity_pattern (stokes_dof_handler, coupling, sp,
                                   stokes_constraints, false,
                                   Utilities::System::
                                   get_this_mpi_process(MPI_COMM_WORLD));
   sp.compress();
 
   stokes_matrix.reinit (sp);
 }
 
 
 
 template <int dim>
 void BoussinesqFlowProblem<dim>::
   setup_stokes_preconditioner (const std::vector<IndexSet> &stokes_partitioning)
 {
   Amg_preconditioner.reset ();
   Mp_preconditioner.reset ();
 
   stokes_preconditioner_matrix.clear ();
 
   TrilinosWrappers::BlockSparsityPattern sp (stokes_partitioning,
                                             MPI_COMM_WORLD);
 
   Table<2,DoFTools::Coupling> coupling (dim+1, dim+1);
   for (unsigned int c=0; c<dim+1; ++c)
     for (unsigned int d=0; d<dim+1; ++d)
       if (c == d)
        coupling[c][d] = DoFTools::always;
       else
        coupling[c][d] = DoFTools::none;
 
   DoFTools::make_sparsity_pattern (stokes_dof_handler, coupling, sp,
                                   stokes_constraints, false,
                                   Utilities::System::
                                   get_this_mpi_process(MPI_COMM_WORLD));
   sp.compress();
 
   stokes_preconditioner_matrix.reinit (sp);
 }
 
 
 template <int dim>
 void BoussinesqFlowProblem<dim>::
   setup_temperature_matrices (const IndexSet &temperature_partitioner)
 {
   T_preconditioner.reset ();
   temperature_mass_matrix.clear ();
   temperature_stiffness_matrix.clear ();
   temperature_matrix.clear ();
 
   TrilinosWrappers::SparsityPattern sp (temperature_partitioner,
                                        MPI_COMM_WORLD);
   DoFTools::make_sparsity_pattern (temperature_dof_handler, sp,
                                   temperature_constraints, false,
                                   Utilities::System::
                                   get_this_mpi_process(MPI_COMM_WORLD));
   sp.compress();
 
   temperature_matrix.reinit (sp);
   temperature_mass_matrix.reinit (sp);
   temperature_stiffness_matrix.reinit (sp);
 }
 
 
 
 template <int dim>
 void BoussinesqFlowProblem<dim>::setup_dofs ()
 {
   computing_timer.enter_section("Setup dof systems");
   std::vector<unsigned int> stokes_sub_blocks (dim+1,0);
   stokes_sub_blocks[dim] = 1;
 
   GridTools::partition_triangulation (Utilities::System::
                                      get_n_mpi_processes(MPI_COMM_WORLD),
                                      triangulation);
 
   {
     stokes_dof_handler.distribute_dofs (stokes_fe);
     DoFRenumbering::subdomain_wise (stokes_dof_handler);
     DoFRenumbering::component_wise (stokes_dof_handler, stokes_sub_blocks);
 
     stokes_constraints.clear ();
     DoFTools::make_hanging_node_constraints (stokes_dof_handler,
                                             stokes_constraints);
 
     std::vector<bool> velocity_mask (dim+1, true);
     velocity_mask[dim] = false;
     VectorTools::interpolate_boundary_values (stokes_dof_handler,
                                              0,
                                              ZeroFunction<dim>(dim+1),
                                              stokes_constraints,
                                              velocity_mask);
 
     std::set<unsigned char> no_normal_flux_boundaries;
     no_normal_flux_boundaries.insert (1);
     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);
     DoFRenumbering::subdomain_wise (temperature_dof_handler);
 
     temperature_constraints.clear ();
     VectorTools::interpolate_boundary_values (temperature_dof_handler,
                                              0,
                                              EquationData::TemperatureInitialValues<dim>(),
                                              temperature_constraints);
     VectorTools::interpolate_boundary_values (temperature_dof_handler,
                                              1,
                                              EquationData::TemperatureInitialValues<dim>(),
                                              temperature_constraints);
     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();
 
   pcout << "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;
 
   std::vector<IndexSet> stokes_partitioning;
   IndexSet temperature_partitioning (n_T);
   {
     const unsigned int my_id =
       Utilities::System::get_this_mpi_process(MPI_COMM_WORLD);
     IndexSet stokes_index_set =
       DoFTools::dof_indices_with_subdomain_association(stokes_dof_handler,
                                                       my_id);
     stokes_partitioning.push_back(stokes_index_set.get_view(0,n_u));
     stokes_partitioning.push_back(stokes_index_set.get_view(n_u,n_u+n_p));
 
     temperature_partitioning =
       DoFTools::dof_indices_with_subdomain_association(temperature_dof_handler,
                                                       my_id);
   }
 
   if (Utilities::System::job_supports_mpi() == false)
     {
       Threads::TaskGroup<> tasks;
       tasks += Threads::new_task (&BoussinesqFlowProblem<dim>::setup_stokes_matrix,
                                  *this,
                                  stokes_partitioning);
       tasks += Threads::new_task (&BoussinesqFlowProblem<dim>::setup_stokes_preconditioner,
                                  *this,
                                  stokes_partitioning);
       tasks += Threads::new_task (&BoussinesqFlowProblem<dim>::setup_temperature_matrices,
                                  *this,
                                  temperature_partitioning);
       tasks.join_all ();
     }
   else
     {
       setup_stokes_matrix (stokes_partitioning);
       setup_stokes_preconditioner (stokes_partitioning);
       setup_temperature_matrices (temperature_partitioning);
     }
 
   stokes_rhs.reinit (stokes_partitioning, MPI_COMM_WORLD);
   stokes_solution.reinit (stokes_rhs);
   old_stokes_solution.reinit (stokes_solution);
 
   temperature_rhs.reinit (temperature_partitioning, MPI_COMM_WORLD);
   temperature_solution.reinit (temperature_rhs);
   old_temperature_solution.reinit (temperature_solution);
   old_old_temperature_solution.reinit (temperature_solution);
 
   computing_timer.exit_section();
 }
 template <int dim>
 void
 BoussinesqFlowProblem<dim>::
 local_assemble_stokes_preconditioner (const typename DoFHandler<dim>::active_cell_iterator &cell,
                                      Assembly::Scratch::StokesPreconditioner<dim> &scratch,
                                      Assembly::CopyData::StokesPreconditioner<dim> &data)
 {
   const unsigned int   dofs_per_cell   = stokes_fe.dofs_per_cell;
   const unsigned int   n_q_points      = scratch.stokes_fe_values.n_quadrature_points;
 
   const FEValuesExtractors::Vector velocities (0);
   const FEValuesExtractors::Scalar pressure (dim);
 
   scratch.stokes_fe_values.reinit (cell);
   cell->get_dof_indices (data.local_dof_indices);
 
   data.local_matrix = 0;
 
   for (unsigned int q=0; q<n_q_points; ++q)
     {
       for (unsigned int k=0; k<dofs_per_cell; ++k)
        {
          scratch.grad_phi_u[k] = scratch.stokes_fe_values[velocities].gradient(k,q);
          scratch.phi_p[k]      = scratch.stokes_fe_values[pressure].value (k, q);
        }
 
       for (unsigned int i=0; i<dofs_per_cell; ++i)
        for (unsigned int j=0; j<dofs_per_cell; ++j)
          data.local_matrix(i,j) += (EquationData::eta *
                                     scalar_product (scratch.grad_phi_u[i],
                                                     scratch.grad_phi_u[j])
                                     +
                                     (1./EquationData::eta) *
                                     EquationData::pressure_scaling *
                                     EquationData::pressure_scaling *
                                     scratch.phi_p[i] * scratch.phi_p[j])
                                    * scratch.stokes_fe_values.JxW(q);
     }
 }
 
 
 
 template <int dim>
 void
 BoussinesqFlowProblem<dim>::
 copy_local_to_global_stokes_preconditioner (const Assembly::CopyData::StokesPreconditioner<dim> &data)
 {
   stokes_constraints.distribute_local_to_global (data.local_matrix,
                                                 data.local_dof_indices,
                                                 stokes_preconditioner_matrix);
 }
 
 
 
 template <int dim>
 void
 BoussinesqFlowProblem<dim>::assemble_stokes_preconditioner ()
 {
   stokes_preconditioner_matrix = 0;
 
   const QGauss<dim> quadrature_formula(stokes_degree+2);
 
   typedef
     FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
     SubdomainFilter;
 
   WorkStream::
     run (SubdomainFilter (IteratorFilters::SubdomainEqualTo
                          (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)),
                          stokes_dof_handler.begin_active()),
         SubdomainFilter (IteratorFilters::SubdomainEqualTo
                          (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)),
                          stokes_dof_handler.end()),
         std_cxx1x::bind (&BoussinesqFlowProblem<dim>::
                          local_assemble_stokes_preconditioner,
                          this,
                          _1,
                          _2,
                          _3),
         std_cxx1x::bind (&BoussinesqFlowProblem<dim>::
                          copy_local_to_global_stokes_preconditioner,
                          this,
                          _1),
         Assembly::Scratch::
         StokesPreconditioner<dim> (stokes_fe, quadrature_formula,
                                    update_JxW_values |
                                    update_values |
                                    update_gradients),
         Assembly::CopyData::
         StokesPreconditioner<dim> (stokes_fe));
 
   stokes_preconditioner_matrix.compress();
 }
 
 
 
 template <int dim>
 void
 BoussinesqFlowProblem<dim>::build_stokes_preconditioner ()
 {
   if (rebuild_stokes_preconditioner == false)
     return;
 
   computing_timer.enter_section ("   Build Stokes preconditioner");
   pcout << "   Rebuilding Stokes preconditioner..." << std::flush;
 
   assemble_stokes_preconditioner ();
 
   std::vector<std::vector<bool> > constant_modes;
   std::vector<bool>  velocity_components (dim+1,true);
   velocity_components[dim] = false;
   DoFTools::extract_constant_modes (stokes_dof_handler, velocity_components,
                                    constant_modes);
 
   Mp_preconditioner  = std_cxx1x::shared_ptr<TrilinosWrappers::PreconditionILU>
                       (new TrilinosWrappers::PreconditionILU());
   Amg_preconditioner = std_cxx1x::shared_ptr<TrilinosWrappers::PreconditionAMG>
                       (new TrilinosWrappers::PreconditionAMG());
 
   TrilinosWrappers::PreconditionAMG::AdditionalData Amg_data;
   Amg_data.constant_modes = constant_modes;
   Amg_data.elliptic = true;
   Amg_data.higher_order_elements = true;
   Amg_data.smoother_sweeps = 2;
   Amg_data.aggregation_threshold = 0.02;
 
   Mp_preconditioner->initialize (stokes_preconditioner_matrix.block(1,1));
   Amg_preconditioner->initialize (stokes_preconditioner_matrix.block(0,0),
                                  Amg_data);
 
   rebuild_stokes_preconditioner = false;
 
   pcout << std::endl;
   computing_timer.exit_section();
 }
 template <int dim>
 void
 BoussinesqFlowProblem<dim>::
 local_assemble_stokes_system (const typename DoFHandler<dim>::active_cell_iterator &cell,
                              Assembly::Scratch::StokesSystem<dim> &scratch,
                              Assembly::CopyData::StokesSystem<dim> &data)
 {
   const unsigned int dofs_per_cell = scratch.stokes_fe_values.get_fe().dofs_per_cell;
   const unsigned int n_q_points    = scratch.stokes_fe_values.n_quadrature_points;
 
   const FEValuesExtractors::Vector velocities (0);
   const FEValuesExtractors::Scalar pressure (dim);
 
   scratch.stokes_fe_values.reinit (cell);
 
   typename DoFHandler<dim>::active_cell_iterator
     temperature_cell (&triangulation,
                      cell->level(),
                      cell->index(),
                      &temperature_dof_handler);
   scratch.temperature_fe_values.reinit (temperature_cell);
 
   if (rebuild_stokes_matrix)
     data.local_matrix = 0;
   data.local_rhs = 0;
 
   scratch.temperature_fe_values.get_function_values (old_temperature_solution,
                                                     scratch.old_temperature_values);
 
   for (unsigned int q=0; q<n_q_points; ++q)
     {
       const double old_temperature = scratch.old_temperature_values[q];
 
       for (unsigned int k=0; k<dofs_per_cell; ++k)
        {
          scratch.phi_u[k] = scratch.stokes_fe_values[velocities].value (k,q);
          if (rebuild_stokes_matrix)
            {
              scratch.grads_phi_u[k] = scratch.stokes_fe_values[velocities].symmetric_gradient(k,q);
              scratch.div_phi_u[k]   = scratch.stokes_fe_values[velocities].divergence (k, q);
              scratch.phi_p[k]       = scratch.stokes_fe_values[pressure].value (k, q);
            }
        }
 
       if (rebuild_stokes_matrix)
        for (unsigned int i=0; i<dofs_per_cell; ++i)
          for (unsigned int j=0; j<dofs_per_cell; ++j)
            data.local_matrix(i,j) += (EquationData::eta * 2 *
                                       (scratch.grads_phi_u[i] * scratch.grads_phi_u[j])
                                       - (EquationData::pressure_scaling *
                                          scratch.div_phi_u[i] * scratch.phi_p[j])
                                       - (EquationData::pressure_scaling *
                                          scratch.phi_p[i] * scratch.div_phi_u[j]))
                                      * scratch.stokes_fe_values.JxW(q);
 
       const Tensor<1,dim>
        gravity = EquationData::gravity_vector (scratch.stokes_fe_values
                                                .quadrature_point(q));
 
       for (unsigned int i=0; i<dofs_per_cell; ++i)
        data.local_rhs(i) += (EquationData::density(old_temperature) *
                              gravity  *
                              scratch.phi_u[i]) *
                             scratch.stokes_fe_values.JxW(q);
     }
 
   cell->get_dof_indices (data.local_dof_indices);
 }
 
 
 
 template <int dim>
 void
 BoussinesqFlowProblem<dim>::
 copy_local_to_global_stokes_system (const Assembly::CopyData::StokesSystem<dim> &data)
 {
   if (rebuild_stokes_matrix == true)
     stokes_constraints.distribute_local_to_global (data.local_matrix,
                                                   data.local_rhs,
                                                   data.local_dof_indices,
                                                   stokes_matrix,
                                                   stokes_rhs);
   else
     stokes_constraints.distribute_local_to_global (data.local_rhs,
                                                   data.local_dof_indices,
                                                   stokes_rhs);
 }
 
 
 
 template <int dim>
 void BoussinesqFlowProblem<dim>::assemble_stokes_system ()
 {
   pcout << "   Assembling..." << std::flush;
 
   computing_timer.enter_section ("   Assemble Stokes system");
 
   if (rebuild_stokes_matrix == true)
     stokes_matrix=0;
 
   stokes_rhs=0;
 
   const QGauss<dim> quadrature_formula(stokes_degree+2);
 
   typedef
     FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
     SubdomainFilter;
 
   WorkStream::
     run (SubdomainFilter (IteratorFilters::SubdomainEqualTo
                          (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)),
                          stokes_dof_handler.begin_active()),
         SubdomainFilter (IteratorFilters::SubdomainEqualTo
                          (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)),
                          stokes_dof_handler.end()),
         std_cxx1x::bind (&BoussinesqFlowProblem<dim>::
                          local_assemble_stokes_system,
                          this,
                          _1,
                          _2,
                          _3),
         std_cxx1x::bind (&BoussinesqFlowProblem<dim>::
                          copy_local_to_global_stokes_system,
                          this,
                          _1),
         Assembly::Scratch::
         StokesSystem<dim> (stokes_fe, quadrature_formula,
                            (update_values    |
                             update_quadrature_points  |
                             update_JxW_values |
                             (rebuild_stokes_matrix == true
                              ?
                              update_gradients
                              :
                              UpdateFlags(0))),
                            temperature_fe,
                            update_values),
         Assembly::CopyData::
         StokesSystem<dim> (stokes_fe));
 
   stokes_matrix.compress();
   stokes_rhs.compress();
 
   rebuild_stokes_matrix = false;
 
   pcout << std::endl;
   computing_timer.exit_section();
 }
 template <int dim>
 void BoussinesqFlowProblem<dim>::
 local_assemble_temperature_matrix (const typename DoFHandler<dim>::active_cell_iterator &cell,
                                   Assembly::Scratch::TemperatureMatrix<dim> &scratch,
                                   Assembly::CopyData::TemperatureMatrix<dim> &data)
 {
   const unsigned int dofs_per_cell = scratch.temperature_fe_values.get_fe().dofs_per_cell;
   const unsigned int n_q_points    = scratch.temperature_fe_values.n_quadrature_points;
 
   scratch.temperature_fe_values.reinit (cell);
   cell->get_dof_indices (data.local_dof_indices);
 
   data.local_mass_matrix = 0;
   data.local_stiffness_matrix = 0;
 
   for (unsigned int q=0; q<n_q_points; ++q)
     {
       for (unsigned int k=0; k<dofs_per_cell; ++k)
        {
          scratch.grad_phi_T[k] = scratch.temperature_fe_values.shape_grad (k,q);
          scratch.phi_T[k]      = scratch.temperature_fe_values.shape_value (k, q);
        }
 
       for (unsigned int i=0; i<dofs_per_cell; ++i)
        for (unsigned int j=0; j<dofs_per_cell; ++j)
          {
            data.local_mass_matrix(i,j)
              += (scratch.phi_T[i] * scratch.phi_T[j]
                  *
                  scratch.temperature_fe_values.JxW(q));
            data.local_stiffness_matrix(i,j)
              += (EquationData::kappa * scratch.grad_phi_T[i] * scratch.grad_phi_T[j]
                  *
                  scratch.temperature_fe_values.JxW(q));
          }
     }
 }
 
 
 
 template <int dim>
 void
 BoussinesqFlowProblem<dim>::
 copy_local_to_global_temperature_matrix (const Assembly::CopyData::TemperatureMatrix<dim> &data)
 {
   temperature_constraints.distribute_local_to_global (data.local_mass_matrix,
                                                      data.local_dof_indices,
                                                      temperature_mass_matrix);
   temperature_constraints.distribute_local_to_global (data.local_stiffness_matrix,
                                                      data.local_dof_indices,
                                                      temperature_stiffness_matrix);
 }
 
 
 template <int dim>
 void BoussinesqFlowProblem<dim>::assemble_temperature_matrix ()
 {
   if (rebuild_temperature_matrices == false)
     return;
 
   computing_timer.enter_section ("   Assemble temperature matrices");
   temperature_mass_matrix = 0;
   temperature_stiffness_matrix = 0;
 
   const QGauss<dim> quadrature_formula(temperature_degree+2);
 
   typedef
     FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
     SubdomainFilter;
 
   WorkStream::
     run (SubdomainFilter (IteratorFilters::SubdomainEqualTo
                          (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)),
                          temperature_dof_handler.begin_active()),
         SubdomainFilter (IteratorFilters::SubdomainEqualTo
                          (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)),
                          temperature_dof_handler.end()),
         std_cxx1x::bind (&BoussinesqFlowProblem<dim>::
                          local_assemble_temperature_matrix,
                          this,
                          _1,
                          _2,
                          _3),
         std_cxx1x::bind (&BoussinesqFlowProblem<dim>::
                          copy_local_to_global_temperature_matrix,
                          this,
                          _1),
         Assembly::Scratch::
         TemperatureMatrix<dim> (temperature_fe, quadrature_formula),
         Assembly::CopyData::
         TemperatureMatrix<dim> (temperature_fe));
 
   temperature_mass_matrix.compress();
   temperature_stiffness_matrix.compress();
 
   rebuild_temperature_matrices = false;
   rebuild_temperature_preconditioner = true;
 
   computing_timer.exit_section();
 }
 template <int dim>
 void BoussinesqFlowProblem<dim>::
 local_assemble_temperature_rhs (const std::pair<double,double> global_T_range,
                                const double                   global_max_velocity,
                                const typename DoFHandler<dim>::active_cell_iterator &cell,
                                Assembly::Scratch::TemperatureRHS<dim> &scratch,
                                Assembly::CopyData::TemperatureRHS<dim> &data)
 {
   const bool use_bdf2_scheme = (timestep_number != 0);
 
   const unsigned int dofs_per_cell = scratch.temperature_fe_values.get_fe().dofs_per_cell;
   const unsigned int n_q_points    = scratch.temperature_fe_values.n_quadrature_points;
 
   EquationData::TemperatureRightHandSide<dim>  temperature_right_hand_side;
 
   const FEValuesExtractors::Vector velocities (0);
 
   data.local_rhs = 0;
   data.matrix_for_bc = 0;
   cell->get_dof_indices (data.local_dof_indices);
 
   scratch.temperature_fe_values.reinit (cell);
 
   typename DoFHandler<dim>::active_cell_iterator
     stokes_cell (&triangulation,
                 cell->level(),
                 cell->index(),
                 &stokes_dof_handler);
   scratch.stokes_fe_values.reinit (stokes_cell);
 
   scratch.temperature_fe_values.get_function_values (old_temperature_solution,
                                                     scratch.old_temperature_values);
   scratch.temperature_fe_values.get_function_values (old_old_temperature_solution,
                                                     scratch.old_old_temperature_values);
 
   scratch.temperature_fe_values.get_function_gradients (old_temperature_solution,
                                                        scratch.old_temperature_grads);
   scratch.temperature_fe_values.get_function_gradients (old_old_temperature_solution,
                                                        scratch.old_old_temperature_grads);
 
   scratch.temperature_fe_values.get_function_laplacians (old_temperature_solution,
                                                         scratch.old_temperature_laplacians);
   scratch.temperature_fe_values.get_function_laplacians (old_old_temperature_solution,
                                                         scratch.old_old_temperature_laplacians);
 
   temperature_right_hand_side.value_list (scratch.temperature_fe_values.get_quadrature_points(),
                                          scratch.gamma_values);
 
   scratch.stokes_fe_values[velocities].get_function_values (stokes_solution,
                                                            scratch.old_velocity_values);
   scratch.stokes_fe_values[velocities].get_function_values (old_stokes_solution,
                                                            scratch.old_old_velocity_values);
 
   const double nu
     = compute_viscosity (scratch.old_temperature_values,
                         scratch.old_old_temperature_values,
                         scratch.old_temperature_grads,
                         scratch.old_old_temperature_grads,
                         scratch.old_temperature_laplacians,
                         scratch.old_old_temperature_laplacians,
                         scratch.old_velocity_values,
                         scratch.old_old_velocity_values,
                         scratch.gamma_values,
                         global_max_velocity,
                         global_T_range.second - global_T_range.first,
                         cell->diameter());
 
   for (unsigned int q=0; q<n_q_points; ++q)
     {
       for (unsigned int k=0; k<dofs_per_cell; ++k)
        {
          scratch.phi_T[k]      = scratch.temperature_fe_values.shape_value (k, q);
          scratch.grad_phi_T[k] = scratch.temperature_fe_values.shape_grad (k,q);
        }
 
 
       const double old_Ts
        = (use_bdf2_scheme ?
           (scratch.old_temperature_values[q] *
            (time_step + old_time_step) / old_time_step
            -
            scratch.old_old_temperature_values[q] *
            (time_step * time_step) /
            (old_time_step * (time_step + old_time_step)))
           :
           scratch.old_temperature_values[q]);
 
       const Tensor<1,dim> ext_grad_T
        = (use_bdf2_scheme ?
           (scratch.old_temperature_grads[q] *
            (1+time_step/old_time_step)
            -
            scratch.old_old_temperature_grads[q] *
            time_step / old_time_step)
           :
           scratch.old_temperature_grads[q]);
 
       const Tensor<1,dim> extrapolated_u
        = (use_bdf2_scheme ?
           (scratch.old_velocity_values[q] * (1+time_step/old_time_step) -
            scratch.old_old_velocity_values[q] * time_step/old_time_step)
           :
           scratch.old_velocity_values[q]);
 
       for (unsigned int i=0; i<dofs_per_cell; ++i)
        {
          data.local_rhs(i) += (old_Ts * scratch.phi_T[i]
                                -
                                time_step *
                                extrapolated_u * ext_grad_T * scratch.phi_T[i]
                                -
                                time_step *
                                nu * ext_grad_T * scratch.grad_phi_T[i]
                                +
                                time_step *
                                scratch.gamma_values[q] * scratch.phi_T[i])
                               *
                               scratch.temperature_fe_values.JxW(q);
 
          if (temperature_constraints.is_inhomogeneously_constrained(data.local_dof_indices[i]))
            {
              for (unsigned int j=0; j<dofs_per_cell; ++j)
                data.matrix_for_bc(j,i) += (scratch.phi_T[i] * scratch.phi_T[j] *
                                            (use_bdf2_scheme ?
                                             ((2*time_step + old_time_step) /
                                              (time_step + old_time_step)) : 1.)
                                            +
                                            scratch.grad_phi_T[i] *
                                            scratch.grad_phi_T[j] *
                                            EquationData::kappa *
                                            time_step)
                                           *
                                           scratch.temperature_fe_values.JxW(q);
            }
        }
     }
 }
 
 
 template <int dim>
 void
 BoussinesqFlowProblem<dim>::
 copy_local_to_global_temperature_rhs (const Assembly::CopyData::TemperatureRHS<dim> &data)
 {
   temperature_constraints.distribute_local_to_global (data.local_rhs,
                                                      data.local_dof_indices,
                                                      temperature_rhs,
                                                      data.matrix_for_bc);
 }
 
 
 
 template <int dim>
 void BoussinesqFlowProblem<dim>::assemble_temperature_system (const double maximal_velocity)
 {
   const bool use_bdf2_scheme = (timestep_number != 0);
 
   if (use_bdf2_scheme == true)
     {
       temperature_matrix.copy_from (temperature_mass_matrix);
       temperature_matrix *= (2*time_step + old_time_step) /
                            (time_step + old_time_step);
       temperature_matrix.add (time_step, temperature_stiffness_matrix);
     }
   else
     {
       temperature_matrix.copy_from (temperature_mass_matrix);
       temperature_matrix.add (time_step, temperature_stiffness_matrix);
     }
   temperature_matrix.compress();
 
   if (rebuild_temperature_preconditioner == true)
     {
       T_preconditioner =  std_cxx1x::shared_ptr<TrilinosWrappers::PreconditionIC>
                          (new TrilinosWrappers::PreconditionIC());
       T_preconditioner->initialize (temperature_matrix);
       rebuild_temperature_preconditioner = false;
     }
 
   temperature_rhs = 0;
 
   const QGauss<dim> quadrature_formula(temperature_degree+2);
   const std::pair<double,double>
     global_T_range = get_extrapolated_temperature_range();
 
   typedef
     FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
     SubdomainFilter;
 
   WorkStream::
     run (SubdomainFilter (IteratorFilters::SubdomainEqualTo
                          (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)),
                          temperature_dof_handler.begin_active()),
         SubdomainFilter (IteratorFilters::SubdomainEqualTo
                          (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)),
                          temperature_dof_handler.end()),
         std_cxx1x::bind (&BoussinesqFlowProblem<dim>::
                          local_assemble_temperature_rhs,
                          this,
                          global_T_range,
                          maximal_velocity,
                          _1,
                          _2,
                          _3),
         std_cxx1x::bind (&BoussinesqFlowProblem<dim>::
                          copy_local_to_global_temperature_rhs,
                          this,
                          _1),
         Assembly::Scratch::
         TemperatureRHS<dim> (temperature_fe, stokes_fe, quadrature_formula),
         Assembly::CopyData::
         TemperatureRHS<dim> (temperature_fe));
 
   temperature_rhs.compress();
 }
 template <int dim>
 void BoussinesqFlowProblem<dim>::solve ()
 {
   computing_timer.enter_section ("   Solve Stokes system");
   pcout << "   Solving..." << std::endl;
 
   {
     const LinearSolvers::BlockSchurPreconditioner<TrilinosWrappers::PreconditionAMG,
                                                   TrilinosWrappers::PreconditionILU>
       preconditioner (stokes_matrix, *Mp_preconditioner, *Amg_preconditioner);
 
     TrilinosWrappers::MPI::BlockVector
       distributed_stokes_solution (stokes_rhs);
     distributed_stokes_solution = stokes_solution;
 
     const unsigned int
       start = (distributed_stokes_solution.block(0).size() +
               distributed_stokes_solution.block(1).local_range().first),
       end   = (distributed_stokes_solution.block(0).size() +
               distributed_stokes_solution.block(1).local_range().second);
     for (unsigned int i=start; i<end; ++i)
       if (stokes_constraints.is_constrained (i))
        distributed_stokes_solution(i) = 0;
 
     SolverControl solver_control (stokes_matrix.m(), 1e-22*stokes_rhs.l2_norm());
     SolverBicgstab<TrilinosWrappers::MPI::BlockVector>
       bicgstab (solver_control, false);
 
     bicgstab.solve(stokes_matrix, distributed_stokes_solution, stokes_rhs,
                   preconditioner);
 
     stokes_solution = distributed_stokes_solution;
 
     pcout << "   "
          << solver_control.last_step()
          << " BiCGStab iterations for Stokes subsystem."
          << std::endl;
 
     stokes_constraints.distribute (stokes_solution);
   }
   computing_timer.exit_section();
 
 
   computing_timer.enter_section ("   Assemble temperature rhs");
 
   old_time_step = time_step;
   const double maximal_velocity = get_maximal_velocity();
   if (maximal_velocity > 1e-10)
     time_step = 1./(1.6*dim*std::sqrt(1.*dim)) /
                temperature_degree *
                GridTools::minimal_cell_diameter(triangulation) /
                maximal_velocity;
   else
     time_step = 1./(1.6*dim*std::sqrt(1.*dim)) /
                temperature_degree *
                GridTools::minimal_cell_diameter(triangulation) /
                1e-10;
 
   pcout << "   Maximal velocity: "
        << maximal_velocity * EquationData::year_in_seconds * 100
        << " cm/year"
        << std::endl;
   pcout << "   " << "Time step: "
        << time_step/EquationData::year_in_seconds
        << " years"
        << std::endl;
 
   temperature_solution = old_temperature_solution;
 
   assemble_temperature_system (maximal_velocity);
 
   computing_timer.exit_section ();
 
   computing_timer.enter_section ("   Solve temperature system");
 
   {
     SolverControl solver_control (temperature_matrix.m(),
                                  1e-12*temperature_rhs.l2_norm());
     SolverCG<TrilinosWrappers::MPI::Vector>   cg (solver_control);
 
     TrilinosWrappers::MPI::Vector
       distributed_temperature_solution (temperature_rhs);
     distributed_temperature_solution = temperature_solution;
 
     cg.solve (temperature_matrix, distributed_temperature_solution,
              temperature_rhs, *T_preconditioner);
 
     temperature_solution = distributed_temperature_solution;
     temperature_constraints.distribute (temperature_solution);
 
     pcout << "   "
          << solver_control.last_step()
          << " CG iterations for temperature" << std::endl;
     computing_timer.exit_section();
 
     double min_temperature = temperature_solution(0),
           max_temperature = temperature_solution(0);
     for (unsigned int i=1; i<temperature_solution.size(); ++i)
       {
        min_temperature = std::min<double> (min_temperature,
                                            temperature_solution(i));
        max_temperature = std::max<double> (max_temperature,
                                            temperature_solution(i));
       }
 
     pcout << "   Temperature range: "
          << min_temperature << ' ' << max_temperature
          << std::endl;
   }
 }
 template <int dim>
 void BoussinesqFlowProblem<dim>::output_results ()
 {
   if (timestep_number % 25 != 0)
     return;
 
   computing_timer.enter_section ("Postprocessing");
 
   if (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD) == 0)
     {
 
       const FESystem<dim> joint_fe (stokes_fe, 1,
                                    temperature_fe, 1,
                                    FE_DGQ<dim>(0), 1);
       DoFHandler<dim> joint_dof_handler (triangulation);
       joint_dof_handler.distribute_dofs (joint_fe);
       Assert (joint_dof_handler.n_dofs() ==
              stokes_dof_handler.n_dofs() +
              temperature_dof_handler.n_dofs() +
              triangulation.n_active_cells(),
              ExcInternalError());
 
       Vector<double> joint_solution (joint_dof_handler.n_dofs());
 
       {
        double minimal_pressure = stokes_solution.block(1)(0);
        for (unsigned int i=0; i<stokes_solution.block(1).size(); ++i)
          minimal_pressure = std::min<double> (stokes_solution.block(1)(i),
                                               minimal_pressure);
 
        std::vector<unsigned int> local_joint_dof_indices (joint_fe.dofs_per_cell);
        std::vector<unsigned int> local_stokes_dof_indices (stokes_fe.dofs_per_cell);
        std::vector<unsigned int> local_temperature_dof_indices (temperature_fe.dofs_per_cell);
 
        typename DoFHandler<dim>::active_cell_iterator
          joint_cell       = joint_dof_handler.begin_active(),
          joint_endc       = joint_dof_handler.end(),
          stokes_cell      = stokes_dof_handler.begin_active(),
          temperature_cell = temperature_dof_handler.begin_active();
        for (; joint_cell!=joint_endc; ++joint_cell, ++stokes_cell, ++temperature_cell)
          {
            joint_cell->get_dof_indices (local_joint_dof_indices);
            stokes_cell->get_dof_indices (local_stokes_dof_indices);
            temperature_cell->get_dof_indices (local_temperature_dof_indices);
 
            for (unsigned int i=0; i<joint_fe.dofs_per_cell; ++i)
              if (joint_fe.system_to_base_index(i).first.first == 0)
                {
                  Assert (joint_fe.system_to_base_index(i).second
                          <
                          local_stokes_dof_indices.size(),
                          ExcInternalError());
 
                  const unsigned int index_in_stokes_fe
                    = joint_fe.system_to_base_index(i).second;
 
                  if (stokes_fe.system_to_component_index(index_in_stokes_fe).first
                      < dim)
                    {
                      joint_solution(local_joint_dof_indices[i])
                        = (stokes_solution(local_stokes_dof_indices
                                           [joint_fe.system_to_base_index(i).second])
                           *
                           EquationData::year_in_seconds
                           *
                           100);
                    }
                  else
                    {
                      Assert (stokes_fe.system_to_component_index(index_in_stokes_fe).first
                              ==
                              dim,
                              ExcInternalError());
 
                      joint_solution(local_joint_dof_indices[i])
                        = ((stokes_solution(local_stokes_dof_indices
                                            [joint_fe.system_to_base_index(i).second])
                            -
                            minimal_pressure)
                           *
                           EquationData::pressure_scaling);
                    }
                }
              else if (joint_fe.system_to_base_index(i).first.first == 1)
                {
                  Assert (joint_fe.system_to_base_index(i).second
                          <
                          local_temperature_dof_indices.size(),
                          ExcInternalError());
                  joint_solution(local_joint_dof_indices[i])
                    = temperature_solution(local_temperature_dof_indices
                                           [joint_fe.system_to_base_index(i).second]);
                }
              else
                {
                  Assert (joint_fe.system_to_base_index(i).first.first == 2,
                          ExcInternalError());
                  Assert (joint_fe.system_to_base_index(i).second
                          == 0,
                          ExcInternalError());
                  joint_solution(local_joint_dof_indices[i])
                    = joint_cell->subdomain_id();
                }
          }
       }
 
       std::vector<std::string> joint_solution_names (dim, "velocity");
       joint_solution_names.push_back ("p");
       joint_solution_names.push_back ("T");
       joint_solution_names.push_back ("partition");
 
       DataOut<dim> data_out;
 
       data_out.attach_dof_handler (joint_dof_handler);
 
       std::vector<DataComponentInterpretation::DataComponentInterpretation>
        data_component_interpretation
        (dim+3, DataComponentInterpretation::component_is_scalar);
       for (unsigned int i=0; i<dim; ++i)
        data_component_interpretation[i]
          = DataComponentInterpretation::component_is_part_of_vector;
 
       data_out.add_data_vector (joint_solution, joint_solution_names,
                                DataOut<dim>::type_dof_data,
                                data_component_interpretation);
       data_out.build_patches (std::min(stokes_degree, temperature_degree));
 
       std::ostringstream filename;
       filename << "solution-" << Utilities::int_to_string(timestep_number, 5)
               << ".vtk";
 
       std::ofstream output (filename.str().c_str());
       data_out.write_vtk (output);
     }
 
   computing_timer.exit_section ();
 }
 template <int dim>
 void BoussinesqFlowProblem<dim>::refine_mesh (const unsigned int max_grid_level)
 {
   computing_timer.enter_section ("Refine mesh structure, part 1");
   Vector<float> local_estimated_error_per_cell (triangulation.n_active_cells());
 
   KellyErrorEstimator<dim>::estimate (temperature_dof_handler,
                                      QGauss<dim-1>(temperature_degree+1),
                                      typename FunctionMap<dim>::type(),
                                      temperature_solution,
                                      local_estimated_error_per_cell,
                                      std::vector<bool>(),
                                      0,
                                      0,
                                      Utilities::System::get_this_mpi_process(MPI_COMM_WORLD));
 
   Vector<double> x_local_estimated_error_per_cell (triangulation.n_active_cells());
   x_local_estimated_error_per_cell = local_estimated_error_per_cell;
   Vector<double> estimated_error_per_cell (triangulation.n_active_cells());
 #ifdef DEAL_II_COMPILER_SUPPORTS_MPI
   MPI_Allreduce (&x_local_estimated_error_per_cell(0),
                 &estimated_error_per_cell(0),
                 triangulation.n_active_cells(), MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
 #else
   estimated_error_per_cell = x_local_estimated_error_per_cell;
 #endif
 
   GridRefinement::refine_and_coarsen_fixed_fraction (triangulation,
                                                     estimated_error_per_cell,
                                                     0.6, 0.2);
   if (triangulation.n_levels() > max_grid_level)
     for (typename Triangulation<dim>::active_cell_iterator
           cell = triangulation.begin_active(max_grid_level);
         cell != triangulation.end(); ++cell)
       cell->clear_refine_flag ();
 
   std::vector<TrilinosWrappers::Vector> x_temperature (2);
   x_temperature[0] = temperature_solution;
   x_temperature[1] = old_temperature_solution;
   TrilinosWrappers::BlockVector x_stokes = stokes_solution;
 
   SolutionTransfer<dim,TrilinosWrappers::Vector>
     temperature_trans(temperature_dof_handler);
   SolutionTransfer<dim,TrilinosWrappers::BlockVector>
     stokes_trans(stokes_dof_handler);
 
   triangulation.prepare_coarsening_and_refinement();
   temperature_trans.prepare_for_coarsening_and_refinement(x_temperature);
   stokes_trans.prepare_for_coarsening_and_refinement(x_stokes);
 
   triangulation.execute_coarsening_and_refinement ();
   computing_timer.exit_section();
 
   setup_dofs ();
 
   computing_timer.enter_section ("Refine mesh structure, part 2");
 
   std::vector<TrilinosWrappers::Vector> tmp (2);
   tmp[0].reinit (temperature_solution);
   tmp[1].reinit (temperature_solution);
   temperature_trans.interpolate(x_temperature, tmp);
 
   temperature_solution = tmp[0];
   old_temperature_solution = tmp[1];
   temperature_constraints.distribute(temperature_solution);
   temperature_constraints.distribute(old_temperature_solution);
 
   TrilinosWrappers::BlockVector x_stokes_new = stokes_solution;
   stokes_trans.interpolate (x_stokes, x_stokes_new);
   stokes_solution = x_stokes_new;
 
   rebuild_stokes_matrix              = true;
   rebuild_stokes_preconditioner      = true;
   rebuild_temperature_matrices       = true;
   rebuild_temperature_preconditioner = true;
 
   computing_timer.exit_section();
 }
 template <int dim>
 void BoussinesqFlowProblem<dim>::run ()
 {
   const unsigned int initial_refinement = (dim == 2 ? 5 : 2);
   const unsigned int n_pre_refinement_steps = (dim == 2 ? 2 : 2);
 
   GridGenerator::hyper_shell (triangulation,
                              Point<dim>(),
                              EquationData::R0,
                              EquationData::R1,
                              12,
                              true);
 
   static HyperShellBoundary<dim> boundary;
   triangulation.set_boundary (0, boundary);
   triangulation.set_boundary (1, boundary);
 
   global_Omega_diameter = GridTools::diameter (triangulation);
 
   triangulation.refine_global (initial_refinement);
 
   setup_dofs();
 
   unsigned int pre_refinement_step = 0;
 
   start_time_iteration:
 
   project_temperature_field ();
 
   timestep_number           = 0;
   time_step = old_time_step = 0;
 
   double time = 0;
 
   do
     {
       pcout << "Timestep " << timestep_number
            << ":  t=" << time/EquationData::year_in_seconds
            << " years"
            << std::endl;
 
       assemble_stokes_system ();
       build_stokes_preconditioner ();
       assemble_temperature_matrix ();
 
       solve ();
 
       pcout << std::endl;
 
       if ((timestep_number == 0) &&
          (pre_refinement_step < n_pre_refinement_steps))
        {
          refine_mesh (initial_refinement + n_pre_refinement_steps);
          ++pre_refinement_step;
          goto start_time_iteration;
        }
       else
        if ((timestep_number > 0) && (timestep_number % 10 == 0))
          refine_mesh (initial_refinement + n_pre_refinement_steps);
 
       output_results ();
 
       time += time_step;
       ++timestep_number;
 
       old_stokes_solution          = stokes_solution;
       old_old_temperature_solution = old_temperature_solution;
       old_temperature_solution     = temperature_solution;
     }
   while (time <= EquationData::end_time);
 }
 int main (int argc, char *argv[])
 {
   try
     {
       deallog.depth_console (0);
 
       Utilities::System::MPI_InitFinalize mpi_initialization(argc, argv);
 
       BoussinesqFlowProblem<2> flow_problem;
       flow_problem.run ();
     }
   catch (std::exception &exc)
     {
       std::cerr << std::endl << std::endl
                 << "----------------------------------------------------"
                 << std::endl;
       std::cerr << "Exception on processing: " << std::endl
                 << exc.what() << std::endl
                 << "Aborting!" << std::endl
                 << "----------------------------------------------------"
                 << std::endl;
 
       return 1;
     }
   catch (...)
     {
       std::cerr << std::endl << std::endl
                 << "----------------------------------------------------"
                 << std::endl;
       std::cerr << "Unknown exception!" << std::endl
                 << "Aborting!" << std::endl
                 << "----------------------------------------------------"
                 << std::endl;
       return 1;
     }
 
   return 0;
 }

deal.II documentation generated on Wed Jul 28 23:05:55 2010 by doxygen 1.5.6