The step-22 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 deals with the Stokes system of equations which reads as follows in non-dimensionalized form:

\begin{eqnarray*} -\textrm{div}\; \varepsilon(\textbf{u}) + \nabla p &=& \textbf{f}, \\ -\textrm{div}\; \textbf{u} &=& 0, \end{eqnarray*}

where $\textbf u$ denotes the velocity of a fluid, $p$ is its pressure, $\textbf f$ are external forces, and $\varepsilon(\textbf{u})= \nabla^s{\textbf{u}}= \frac 12 \left[ (\nabla \textbf{u}) + (\nabla \textbf{u})^T\right]$ is the rank-2 tensor of symmetrized gradients; a component-wise definition of it is $\varepsilon(\textbf{u})_{ij}=\frac 12\left(\frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i}\right)$.

The Stokes equations describe the steady-state motion of a slow-moving, viscous fluid such as honey, rocks in the earth mantle, or other cases where inertia does not play a significant role. If a fluid is moving fast enough that inertia forces are significant compared to viscous friction, the Stokes equations are no longer valid; taking into account inertia effects then leads to the nonlinear Navier-Stokes equations. However, in this tutorial program, we will focus on the simpler Stokes system.

To be well-posed, we will have to add boundary conditions to the equations. What boundary conditions are readily possible here will become clear once we discuss the weak form of the equations.

The equations covered here fall into the class of vector-valued problems. A toplevel overview of this topic can be found in the Handling vector valued problems module.

Weak form

The weak form of the equations is obtained by writing it in vector form as

\begin{eqnarray*} \left( {-\textrm{div}\; \varepsilon(\textbf{u}) + \nabla p} \atop {-\textrm{div}\; \textbf{u}} \right) = \left( {\textbf{f}} \atop 0 \right), \end{eqnarray*}

forming the dot product from the left with a vector-valued test function $\phi = \left({\textbf v \atop q}\right)$ and integrating over the domain $\Omega$, yielding the following set of equations:

\begin{eqnarray*} (\mathrm v, -\textrm{div}\; \varepsilon(\textbf{u}) + \nabla p)_{\Omega} - (q,\textrm{div}\; \textbf{u})_{\Omega} = (\textbf{v}, \textbf{f})_\Omega, \end{eqnarray*}

which has to hold for all test functions $\phi = \left({\textbf v \atop q}\right)$.

In practice, one wants to impose as little regularity on the pressure variable as possible; consequently, we integrate by parts the second term:

\begin{eqnarray*} (\mathrm v, -\textrm{div}\; \varepsilon(\textbf{u}))_{\Omega} - (\textrm{div}\; \textbf{v}, p)_{\Omega} + (\textbf{n}\cdot\textbf{v}, p)_{\partial\Omega} - (q,\textrm{div}\; \textbf{u})_{\Omega} = (\textbf{v}, \textbf{f})_\Omega. \end{eqnarray*}

Likewise, we integrate by parts the first term to obtain

\begin{eqnarray*} (\nabla \mathrm v,\varepsilon(\textbf{u}))_{\Omega} - (\textbf{n} \otimes \mathrm v,\varepsilon(\textbf{u}))_{\partial\Omega} - (\textrm{div}\; \textbf{v}, p)_{\Omega} + (\textbf{n}\cdot\textbf{v}, p)_{\partial\Omega} - (q,\textrm{div}\; \textbf{u})_{\Omega} = (\textbf{v}, \textbf{f})_\Omega, \end{eqnarray*}

where the scalar product between two tensor-valued quantities is here defined as

\begin{eqnarray*} (\nabla \mathrm v,\varepsilon(\textbf{u}))_{\Omega} = \int_\Omega \sum_{i,j=1}^d \frac{\partial v_j}{\partial x_i} \varepsilon(\textbf{u})_{ij} \ dx. \end{eqnarray*}

Because the scalar product between a general tensor like $\nabla\mathrm v$ and a symmetric tensor like $\varepsilon(\textbf{u})$ equals the scalar product between the symmetrized forms of the two, we can also write the bilinear form above as follows:

\begin{eqnarray*} (\varepsilon(\mathrm v),\varepsilon(\textbf{u}))_{\Omega} - (\textbf{n} \otimes \mathrm v,\varepsilon(\textbf{u}))_{\partial\Omega} - (\textrm{div}\; \textbf{v}, p)_{\Omega} + (\textbf{n}\cdot\textbf{v}, p)_{\partial\Omega} - (q,\textrm{div}\; \textbf{u})_{\Omega} = (\textbf{v}, \textbf{f})_\Omega, \end{eqnarray*}

We will deal with the boundary terms in the next section, but it is already clear from the domain terms

\begin{eqnarray*} (\varepsilon(\mathrm v),\varepsilon(\textbf{u}))_{\Omega} - (\textrm{div}\; \textbf{v}, p)_{\Omega} - (q,\textrm{div}\; \textbf{u})_{\Omega} \end{eqnarray*}

of the bilinear form that the Stokes equations yield a symmetric bilinear form, and consequently a symmetric (if indefinite) system matrix.

Boundary conditions

The weak form just derived immediately presents us with different possibilities for imposing boundary conditions:

  1. Dirichlet velocity boundary conditions: On a part $\Gamma_D\subset\partial\Omega$ we may impose Dirichlet conditions on the velocity $\textbf u$:

    \begin{eqnarray*} \textbf u = \textbf g_D \qquad\qquad \textrm{on}\ \Gamma_D. \end{eqnarray*}

    Because test functions $\textbf v$ come from the tangent space of the solution variable, we have that $\textbf v=0$ on $\Gamma_D$ and consequently that

    \begin{eqnarray*} -(\textbf{n} \otimes \mathrm v,\varepsilon(\textbf{u}))_{\Gamma_D} + (\textbf{n}\cdot\textbf{v}, p)_{\Gamma_D} = 0. \end{eqnarray*}

    In other words, as usual, strongly imposed boundary values do not appear in the weak form.

    It is noteworthy that if we impose Dirichlet boundary values on the entire boundary, then the pressure is only determined up to a constant. An algorithmic realization of that would use similar tools as have been seen in step-11.

  2. Neumann-type or natural boundary conditions: On the rest of the boundary $\Gamma_N=\partial\Omega\backslash\Gamma_D$, let us re-write the boundary terms as follows:

    \begin{eqnarray*} -(\textbf{n} \otimes \mathrm v,\varepsilon(\textbf{u}))_{\Gamma_N} + (\textbf{n}\cdot\textbf{v}, p)_{\Gamma_N} &=& \sum_{i,j=1}^d -(n_i v_j,\varepsilon(\textbf{u})_{ij})_{\Gamma_N} + \sum_{i=1}^d (n_i v_i, p)_{\Gamma_N} \\ &=& \sum_{i,j=1}^d -(n_i v_j,\varepsilon(\textbf{u})_{ij})_{\Gamma_N} + \sum_{i,j=1}^d (n_i v_j, p \delta_{ij})_{\Gamma_N} \\ &=& \sum_{i,j=1}^d (n_i v_j,p \delta_{ij} - \varepsilon(\textbf{u})_{ij})_{\Gamma_N} \\ &=& (\textbf{n} \otimes \mathrm v, p \textbf{1} - \varepsilon(\textbf{u}))_{\Gamma_N}. \\ &=& (\mathrm v, \textbf{n}\cdot [p \textbf{1} - \varepsilon(\textbf{u})])_{\Gamma_N}. \end{eqnarray*}

    In other words, on the Neumann part of the boundary we can prescribe values for the total stress:

    \begin{eqnarray*} \textbf{n}\cdot [p \textbf{1} - \varepsilon(\textbf{u})] = \textbf g_N \qquad\qquad \textrm{on}\ \Gamma_N. \end{eqnarray*}

    If the boundary is subdivided into Dirichlet and Neumann parts $\Gamma_D,\Gamma_N$, this then leads to the following weak form:

    \begin{eqnarray*} (\varepsilon(\mathrm v),\varepsilon(\textbf{u}))_{\Omega} - (\textrm{div}\; \textbf{v}, p)_{\Omega} - (q,\textrm{div}\; \textbf{u})_{\Omega} = (\textbf{v}, \textbf{f})_\Omega - (\textbf{v}, \textbf g_N)_{\Gamma_N}. \end{eqnarray*}

  3. Robin-type boundary conditions: Robin boundary conditions are a mixture of Dirichlet and Neumann boundary conditions. They would read

    \begin{eqnarray*} \textbf{n}\cdot [p \textbf{1} - \varepsilon(\textbf{u})] = \textbf S \textbf u \qquad\qquad \textrm{on}\ \Gamma_R, \end{eqnarray*}

    with a rank-2 tensor (matrix) $\textbf S$. The associated weak form is

    \begin{eqnarray*} (\varepsilon(\mathrm v),\varepsilon(\textbf{u}))_{\Omega} - (\textrm{div}\; \textbf{v}, p)_{\Omega} - (q,\textrm{div}\; \textbf{u})_{\Omega} + (\textbf S \textbf u, \textbf{v})_{\Gamma_R} = (\textbf{v}, \textbf{f})_\Omega. \end{eqnarray*}

  4. Partial boundary conditions: It is possible to combine Dirichlet and Neumann boundary conditions by only enforcing each of them for certain components of the velocity. For example, one way to impose artificial boundary conditions is to require that the flow is perpendicular to the boundary, i.e. the tangential component $\textbf u_{\textbf t}=(\textbf 1-\textbf n\otimes\textbf n)\textbf u$ be zero, thereby constraining dim-1 components of the velocity. The remaining component can be constrained by requiring that the normal stress be zero, yielding the following set of boundary conditions:

    \begin{eqnarray*} \textbf u_{\textbf t} &=& 0, \\ \textbf n \cdot \left(\textbf{n}\cdot [p \textbf{1} - \varepsilon(\textbf{u})] \right) &=& 0. \end{eqnarray*}

Despite this wealth of possibilities, we will only use Dirichlet and (homogenous) Neumann boundary conditions in this tutorial program.

Discretization

As developed above, the weak form of the equations with Dirichlet and Neumann boundary conditions on $\Gamma_D$ and $\Gamma_N$ reads like this: find $\textbf u\in \textbf V_g = \{\varphi \in H^1(\Omega)^d: \varphi_{\Gamma_D}=\textbf g_D\}, p\in Q=L^2(\Omega)$ so that

\begin{eqnarray*} (\varepsilon(\mathrm v),\varepsilon(\textbf{u}))_{\Omega} - (\textrm{div}\; \textbf{v}, p)_{\Omega} - (q,\textrm{div}\; \textbf{u})_{\Omega} = (\textbf{v}, \textbf{f})_\Omega - (\textbf{v}, \textbf g_N)_{\Gamma_N} \end{eqnarray*}

for all test functions $\textbf v\in \textbf V_0 = \{\varphi \in H^1(\Omega)^d: \varphi_{\Gamma_D}=0\},q\in Q$.

These equations represent a symmetric saddle point problem. It is well known that then a solution only exists if the function spaces in which we search for a solution have to satisfy certain conditions, typically referred to as the Babuska-Brezzi or Ladyzhenskaya-Babuska-Brezzi (LBB) conditions. The continuous function spaces above satisfy them. However, when we discretize the equations by replacing the continuous variables and test functions by finite element functions in finite dimensional spaces $\textbf V_{g,h}\subset \textbf V_g, Q_h\subset Q$, we have to make sure that $\textbf V_h,Q_h$ also satisfy the LBB conditions. This is similar to what we had to do in step-20.

For the Stokes equations, there are a number of possible choices to ensure that the finite element spaces are compatible with the LBB condition. A simple and accurate choice that we will use here is $\textbf u_h\in Q_{p+1}^d, p_h\in Q_p$, i.e. use elements one order higher for the velocities than for the pressures.

This then leads to the following discrete problem: find $\textbf u_h,p_h$ so that

\begin{eqnarray*} (\varepsilon(\mathrm v_h),\varepsilon(\textbf u_h))_{\Omega} - (\textrm{div}\; \textbf{v}_h, p_h)_{\Omega} - (q_h,\textrm{div}\; \textbf{u}_h)_{\Omega} = (\textbf{v}_h, \textbf{f})_\Omega - (\textbf{v}_h, \textbf g_N)_{\Gamma_N} \end{eqnarray*}

for all test functions $\textbf v_h, q_h$. Assembling the linear system associated with this problem follows the same lines used in step-20, step-21, and explained in detail in the Handling vector valued problems module.

Linear solver and preconditioning issues

The weak form of the discrete equations naturally leads to the following linear system for the nodal values of the velocity and pressure fields:

\begin{eqnarray*} \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 \\ G \end{array}\right), \end{eqnarray*}

Like in step-20 and step-21, we will solve this system of equations by forming the Schur complement, i.e. we will first find the solution $P$ of

\begin{eqnarray*} BA^{-1}B^T P &=& BA^{-1} F - G, \\ \end{eqnarray*}

and then

\begin{eqnarray*} AU &=& F - B^TP. \end{eqnarray*}

The way we do this is pretty much exactly like we did in these previous tutorial programs, i.e. we use the same classes SchurComplement and InverseMatrix again. There are two significant differences, however:

  1. First, in the mixed Laplace equation we had to deal with the question of how to precondition the Schur complement $B^TM^{-1}B$, which was spectrally equivalent to the Laplace operator on the pressure space (because $B$ represents the gradient operator, $B^T$ its adjoint $-\textrm{div}$, and $M$ the identity (up to the material parameter $K^{-1}$), so $B^TM^{-1}B$ is something like $-\textrm{div} \mathbf 1 \nabla = -\Delta$). Consequently, the matrix is badly conditioned for small mesh sizes and we had to come up with an elaborate preconditioning scheme for the Schur complement.

  2. Second, every time we multiplied with $B^TM^{-1}B$ we had to solve with the mass matrix $M$. This wasn't particularly difficult, however, since the mass matrix is always well conditioned and so simple to invert using CG and a little bit of preconditioning.
In other words, preconditioning the inner solver for $M$ was simple whereas preconditioning the outer solver for $B^TM^{-1}B$ was complicated.

Here, the situation is pretty much exactly the opposite. The difference stems from the fact that the matrix at the heart of the Schur complement does not stem from the identity operator but from a variant of the Laplace operator, $-\textrm{div} \nabla^s$ (where $\nabla^s$ is the symmetric gradient) acting on a vector field. In the investigation of this issue we largely follow the paper D. Silvester and A. Wathen: "Fast iterative solution of stabilised Stokes systems part II. Using general block preconditioners." (SIAM J. Numer. Anal., 31 (1994), pp. 1352-1367), which is available online here. Principally, the difference in the matrix at the heart of the Schur complement has two consequences:

  1. First, it makes the outer preconditioner simple: the Schur complement corresponds to the operator $-\textrm{div} (-\textrm{div} \nabla^s)^{-1} \nabla$ on the pressure space; forgetting about the fact that we deal with symmetric gradients instead of the regular one, the Schur complement is something like $-\textrm{div} (-\textrm{div} \nabla)^{-1} \nabla = -\textrm{div} (-\Delta)^{-1} \nabla$, which, even if not mathematically entirely concise, is spectrally equivalent to the identity operator (a heuristic argument would be to commute the operators into $-\textrm{div}(-\Delta)^{-1} \nabla = -\textrm{div}\nabla(-\Delta)^{-1} = -\Delta(-\Delta)^{-1} = \mathbf 1$). It turns out that it isn't easy to solve this Schur complement in a straightforward way with the CG method: using no preconditioner, the condition number of the Schur complement matrix depends on the size ratios of the largest to the smallest cells, and one still needs on the order of 50-100 CG iterations. However, there is a simple cure: precondition with the mass matrix on the pressure space and we get down to a number between 5-15 CG iterations, pretty much independently of the structure of the mesh (take a look at the results section of this program to see that indeed the number of CG iterations does not change as we refine the mesh).

    So all we need in addition to what we already have is the mass matrix on the pressure variables. We could do that by building this matrix on the side in a separate data structure. However, it is worth remembering that although we build the system matrix

    \begin{eqnarray*} \left(\begin{array}{cc} A & B^T \\ B & 0 \end{array}\right) \end{eqnarray*}

    as one object (of type BlockSparseMatrix), we never actually do matrix-vector products with this matrix, or any other operations that consider the entire matrix. Rather, we only build it in this form for convenience (because it reflects the structure of the FESystem finite element and associated DoFHandler object) but later only operate on the $(0,0),(0,1)$, and $(1,0)$ blocks of this matrix. In other words, our algorithm so far entirely ignores the $(1,1)$ (pressure-pressure) block as it is empty anyway.

    Now, as mentioned, we need a pressure mass matrix to precondition the Schur complement and that conveniently the pressure-pressure block of the matrix we build anyway is currently empty and ignored. So what we will do is to assemble the needed mass matrix in this space; this does change the global system matrix but since our algorithm never operates on the global matrix and instead only considers individual blocks, this fact does not affect what we actually compute. Later, when solving, we then precondition the Schur complement with $M_p^{-1}$ by doing a few CG iterations on the well-conditioned pressure mass matrix $M_p$ stored in the $(1,1)$ block.

  2. While the outer preconditioner has become simpler compared to the mixed Laplace case discussed in step-20, the issue of the inner solver has become more complicated. In the mixed Laplace discretization, the Schur complement has the form $B^TM^{-1}B$. Thus, every time we multiplied with the Schur complement, we had to solve a linear system $M_uz=y$; this isn't too complicated there, however, since the mass matrix $M_u$ on the pressure space is well-conditioned.

    On the other hand, for the Stokes equation we consider here, the Schur complement is $B^TA^{-1}B$ where the matrix $A$ is related to the Laplace operator (it is, in fact, the matrix corresponding to the bilinear form $(\nabla^s \varphi_i, \nabla^s\varphi_j)$). Thus, solving with $A$ is a lot more complicated: the matrix is badly conditioned and we know that we need many iterations unless we have a very good preconditioner. What is worse, we have to solve with $A$ every time we multiply with the Schur complement, which is 5-15 times using the preconditioner described above.

    Because we have to solve with $A$ several times, it pays off to spend a bit more time once to create a good preconditioner for this matrix. So here's what we're going to do: if in 2d, we use the ultimate preconditioner, namely a direct sparse LU decomposition of the matrix. This is implemented using the SparseDirectUMFPACK class that uses the UMFPACK direct solver to compute the decomposition. To use it, you will have to specify the --enable-umfpack switch when configuring the deal.II library, see the ReadMe file for instructions. With this, the inner solver converges in one iteration.

    In 2d, we can do this sort of thing because even reasonably large problems rarely have more than a few 100,000 unknowns with relatively few nonzero entries per row. Furthermore, the bandwidth of matrices in 2d is ${\cal O}(\sqrt{N})$ and therefore moderate. For such matrices, sparse factors can be computed in a matter of a few seconds. (As a point of reference, computing the sparse factors of a matrix of size $N$ and bandwidth $B$ takes ${\cal O}(NB^2)$ operations. In 2d, this is ${\cal O}(N^2)$; though this is a higher complexity than, for example, assembling the linear system which takes ${\cal O}(N)$, the constant for computing the decomposition is so small that it doesn't become the dominating factor in the entire program until we get to very large numbers of unknowns in the high 100,000s or more.)

    The situation changes in 3d, because there we quickly have many more unknowns and the bandwidth of matrices (which determines the number of nonzero entries in sparse LU factors) is ${\cal O}(N^{2/3})$, and there are many more entries per row as well. This makes using a sparse direct solver such as UMFPACK inefficient: only for problem sizes of a few 10,000 to maybe 100,000 unknowns can a sparse decomposition be computed using reasonable time and memory resources.

    What we do in that case is to use an incomplete LU decomposition (ILU) as a preconditioner, rather than actually computing complete LU factors. As it so happens, deal.II has a class that does this: SparseILU. Computing the ILU takes a time that only depends on the number of nonzero entries in the sparse matrix (or that we are willing to fill in the LU factors, if these should be more than the ones in the matrix), but is independent of the bandwidth of the matrix. It is therefore an operation that can efficiently also be computed in 3d. On the other hand, an incomplete LU decomposition, by definition, does not represent an exact inverse of the matrix $A$. Consequently, preconditioning with the ILU will still require more than one iteration, unlike preconditioning with the sparse direct solver. The inner solver will therefore take more time when multiplying with the Schur complement, a tradeoff unavoidable.

In the program below, we will make use of the fact that the SparseILU and SparseDirectUMFPACK classes have a very similar interface and can be used interchangeably. All that we need is a switch class that, depending on the dimension, provides a type that is either of the two classes mentioned above. This is how we do that:

template <int dim>
struct InnerPreconditioner;

template <>
struct InnerPreconditioner<2> 
{
    typedef SparseDirectUMFPACK type;
};

template <>
struct InnerPreconditioner<3> 
{
    typedef SparseILU<double> type;
};

From hereon, we can refer to the type typename InnerPreconditioner<dim>type and automatically get the correct preconditioner class. Because of the similarity of the interfaces of the two classes, we will be able to use them interchangeably using the same syntax in all places.

The testcase

The domain, right hand side and boundary conditions we implement below relate to a problem in geophysics: there, one wants to compute the flow field of magma in the earth's interior under a mid-ocean rift. Rifts are places where two continental plates are very slowly drifting apart (a few centimeters per year at most), leaving a crack in the earth crust that is filled with magma from below. Without trying to be entirely realistic, we model this situation by solving the following set of equations and boundary conditions on the domain $\Omega=[-2,2]\times[0,1]\times[-1,0]$:

\begin{eqnarray*} -\textrm{div}\; \varepsilon(\textbf{u}) + \nabla p &=& 0, \\ -\textrm{div}\; \textbf{u} &=& 0, \\ \mathbf u &=& \left(\begin{array}{c} -1 \\ 0 \\0 \end{array}\right) \qquad\qquad \textrm{at}\ z=0, x<0, \\ \mathbf u &=& \left(\begin{array}{c} +1 \\ 0 \\0 \end{array}\right) \qquad\qquad \textrm{at}\ z=0, x>0, \\ \mathbf u &=& \left(\begin{array}{c} 0 \\ 0 \\0 \end{array}\right) \qquad\qquad \textrm{at}\ z=0, x=0, \end{eqnarray*}

and using natural boundary conditions $\textbf{n}\cdot [p \textbf{1} - \varepsilon(\textbf{u})] = 0$ everywhere else. In other words, at the left part of the top surface we prescribe that the fluid moves with the continental plate to the left at speed $-1$, that it moves to the right on the right part of the top surface, and impose natural flow conditions everywhere else. If we are in 2d, the description is essentially the same, with the exception that we omit the second component of all vectors stated above.

As will become apparent in the results section, the flow field will pull material from below and move it to the left and right ends of the domain, as expected. The discontinuity of velocity boundary conditions will produce a singularity in the pressure at the center of the top surface that sucks material all the way to the top surface to fill the gap left by the outward motion of material at this location.

Implementation

The program developed below has seen a lot of TLC. We have run it over and over under profiling tools (mainly valgrind's cachegrind and callgrind tools, as well as the KDE KCachegrind program for visualization) to see where the bottlenecks are. This has paid off: through this effort, the program has become almost twice as fast when considering the runtime of the refinement cycles zero through three, reducing the overall number of CPU instructions executed from 869,574,060,348 to 474,507,755,764. For higher refinement levels, the gain is probably even larger since some algorithms that are not ${\cal O}(N)$ have been eliminated.

Essentially, there are currently two algorithms in the program that do not scale linearly with the number of degrees of freedom: renumbering of degrees of freedom, and the linear solver. As for the first, while reordering degrees of freedom may not scale linearly, it is an indispensible part of the overall algorithm as it greatly improves the quality of the sparse ILU, easily making up for the time spent on computing the renumbering; graphs and timings to demonstrate this are shown in the documentation of the DoFRenumbering namespace, also underlining the choice of the King reordering algorithm chosen below.

As for the linear solver: as mentioned above, our implementation here uses a Schur complement formulation. This is not necessarily the very best choice but demonstrates various important techniques available in deal.II. The question of which solver is best is again discussed in the section on improved solvers in the results part of this program, along with code showing alternative solvers and a comparison of their results.

Apart from this, many other algorithms have been tested and improved during the creation of this program. For example, in building the sparsity pattern, we originally used a BlockCompressedSparsityPattern object; however, its data structures are poorly adapted for the large numbers of nonzero entries per row created by our discretization in 3d, leading to a quadratic behavior. Replacing it with a BlockCompressedSetSparsityPattern removed this bottleneck at the price of some more memory consumption, by using a better adapted data structure. Likewise, the implementation of the decomposition step in the SparseILU class was very inefficient and has been replaced by one that is about 10 times faster. Small improvements were applied here and there.

A profile of how many CPU instructions are spent at the various different places in the program during refinement cycles zero through three in 3d is shown here:

step-22.profile-3.png

As can be seen, at this refinement level approximately one sixth of the instruction count is spent on matrix assembly and sparse ILU computation (third quarter from left to right), about one half on the actual solver (the SparseILU::vmult calls on the left), and the rest on other things. Since floating point operations such as in the SparseILU::vmult calls typically take much longer than many of the logical operations and table lookups in matrix assembly, the fraction of the run time taken up by matrix assembly is actually significantly less than the fraction of instructions, as will become apparent in the comparison we make in the results section.

For higher refinement levels, the boxes representing the solver as well as the beige box at the top right representing the stemming from reordering algorithm are going to grow at the expense of the other parts of the program, since they don't scale linearly. The fact that at this moderate refinement level (3168 cells and 93176 degrees of freedom) the linear solver already makes up more than half the instructions is a good sign that most of the algorithms used in this program are well-tuned and that major improvements in speeding up the program are most likely not to come from hand-optimizing individual aspects but by changing solver algorithms. We will address this point in the discussion of results below as well.

As a final point, and as a point of reference, the following picture also shows how the profile looked at an early stage of optimizing this program:

step-22.profile-3.original.png

As mentioned above, the runtime of this version was about twice as long as for the first profile, with the SparseILU decomposition taking up about 30% of the instruction count, and operations on the ill-suited CompressedSparsityPattern about 10%. Both these bottlenecks have since been completely removed with the exception of the use of the CompressedSparsityPattern through the DoFRenumbering::Cuthill_McKee algorithm. (The latter at least applies as of the time of writing of this program; the problem may be solved in the future.)

The commented program

Include files

As usual, we start by including some well-known files:

 #include <base/quadrature_lib.h>
 #include <base/logstream.h>
 #include <base/function.h>
 #include <base/utilities.h>
 
 #include <lac/block_vector.h>
 #include <lac/full_matrix.h>
 #include <lac/block_sparse_matrix.h>
 #include <lac/solver_cg.h>
 #include <lac/precondition.h>
 
 #include <grid/tria.h>
 #include <grid/grid_generator.h>
 #include <grid/tria_accessor.h>
 #include <grid/tria_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 <dofs/dof_constraints.h>
 
 #include <fe/fe_q.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>

Then we need to include the header file for the sparse direct solver UMFPACK:

 #include <lac/sparse_direct.h>

This includes the libary for the incomplete LU factorization that will be used as a preconditioner in 3D:

 #include <lac/sparse_ilu.h>

This is C++:

 #include <fstream>
 #include <sstream>

As in all programs, the namespace dealii in included:

 using namespace dealii;

Defining the inner preconditioner type

As explained in the introduction, we are going to use different preconditioners for two and three space dimensions, respectively. We differentiate between them by the use of the spatial dimension as a template parameter. See step-4 for details on templates. We are not going to create any preconditioner object here, all we do is to create class that holds a local typedef determining the preconditioner class so we can write our program in a dimension-independent way.

 template <int dim>
 struct InnerPreconditioner;

In 2D, we are going to use a sparse direct solver as preconditioner:

 template <>
 struct InnerPreconditioner<2> 
 {
     typedef SparseDirectUMFPACK type;
 };

And the ILU preconditioning in 3D, called by SparseILU:

 template <>
 struct InnerPreconditioner<3> 
 {
     typedef SparseILU<double> type;
 };

The StokesProblem class template

This is an adaptation of step-20, so the main class and the data types are the same as used there. In this example we also use adaptive grid refinement, which is handled in complete analogy to step-6:

 template <int dim>
 class StokesProblem 
 {
   public:
     StokesProblem (const unsigned int degree);
     void run ();
     
   private:
     void setup_dofs ();
     void assemble_system ();
     void solve ();
     void output_results (const unsigned int refinement_cycle) const;
     void refine_mesh ();
     
     const unsigned int   degree;
     
     Triangulation<dim>   triangulation;
     FESystem<dim>        fe;
     DoFHandler<dim>      dof_handler;
 
     ConstraintMatrix     hanging_node_constraints;
     
     BlockSparsityPattern      sparsity_pattern;
     BlockSparseMatrix<double> system_matrix;
 
     BlockVector<double> solution;
     BlockVector<double> system_rhs;

This one is new: We shall use a so-called shared pointer structure to access the preconditioner. Shared pointers are essentially just a convenient form of pointers. Several shared pointers can point to the same object (just like regular pointers), but when the last shared pointer object to point to a preconditioner object is deleted (for example if a shared pointer object goes out of scope, if the class of which it is a member is destroyed, or if the pointer is assigned a different preconditioner object) then the preconditioner object pointed to is also destroyed. This ensures that we don't have to manually track in how many places a preconditioner object is still referenced, it can never create a memory leak, and can never produce a dangling pointer to an already destroyed object:

     boost::shared_ptr<typename InnerPreconditioner<dim>::type> A_preconditioner;
 };

Boundary values and right hand side

As in step-20 and most other example programs, the next task is to define the data for the PDE: For the Stokes problem, we are going to use natural boundary values on parts of the boundary (i.e. homogenous Neumann-type) for which we won't have to do anything special (the homogeneity implies that the corresponding terms in the weak form are simply zero), and boundary conditions on the velocity (Dirichlet-type) on the rest of the boundary, as described in the introduction.

In order to enforce the Dirichlet boundary values on the velocity, we will use the VectorTools::interpolate_boundary_values function as usual which requires us to write a function object with as many components as the finite element has. In other words, we have to define the function on the $(u,p)$-space, but we are going to filter out the pressure component when interpolating the boundary values.

The following function object is a representation of the boundary values described in the introduction:

 template <int dim>
 class BoundaryValues : public Function<dim> 
 {
   public:
     BoundaryValues () : Function<dim>(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
 BoundaryValues<dim>::value (const Point<dim>  &p,
                             const unsigned int component) const 
 {
   Assert (component < this->n_components,
           ExcIndexRange (component, 0, this->n_components));
   
   if (component == 0)
     return (p[0] < 0 ? -1 : (p[0] > 0 ? 1 : 0));
   return 0;
 }
 
 
 template <int dim>
 void
 BoundaryValues<dim>::vector_value (const Point<dim> &p,
                                    Vector<double>   &values) const 
 {
   for (unsigned int c=0; c<this->n_components; ++c)
     values(c) = BoundaryValues<dim>::value (p, c);
 }

We implement similar functions for the right hand side which for the current example is simply zero:

 template <int dim>
 class RightHandSide : public Function<dim> 
 {
   public:
     RightHandSide () : Function<dim>(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
 RightHandSide<dim>::value (const Point<dim>  &/ *p* /,
                            const unsigned int / *component* /) const 
 {
   return 0;
 }
 
 
 template <int dim>
 void
 RightHandSide<dim>::vector_value (const Point<dim> &p,
                                   Vector<double>   &values) const 
 {
   for (unsigned int c=0; c<this->n_components; ++c)
     values(c) = RightHandSide<dim>::value (p, c);
 }

Linear solvers and preconditioners

The linear solvers and preconditioners are discussed extensively in the introduction. Here, we create the respective objects that will be used.

The InverseMatrix class template

The InverseMatrix class represents the data structure for an inverse matrix. It is derived from the one in step-20. The only difference is that we now do include a preconditioner to the matrix since we will apply this class to different kinds of matrices that will require different preconditioners (in step-20 we did not use a preconditioner in this class at all). The types of matrix and preconditioner are passed to this class via template parameters, and matrix and preconditioner objects of these types will then be passed to the constructor when an InverseMatrix object is created. The member function vmult is, as in step-20, a multiplication with a vector, obtained by solving a linear system:

 template <class Matrix, class Preconditioner>
 class InverseMatrix : public Subscriptor
 {
   public:
     InverseMatrix (const Matrix         &m,
                    const Preconditioner &preconditioner);
 
     void vmult (Vector<double>       &dst,
                 const Vector<double> &src) const;
 
   private:
     const SmartPointer<const Matrix> matrix;
     const SmartPointer<const Preconditioner> preconditioner;
 };
 
 
 template <class Matrix, class Preconditioner>
 InverseMatrix<Matrix,Preconditioner>::InverseMatrix (const Matrix &m,
                                                      const Preconditioner &preconditioner)
                 :
                 matrix (&m),
                 preconditioner (&preconditioner)
 {}

This is the implementation of the vmult function.

Note that we use a rather large tolerance for the solver control. The reason for this is that the function is used very frequently, and hence, any additional effort to make the residual in the CG solve smaller makes the solution more expensive. Note that we do not only use this class as a preconditioner for the Schur complement, but also when forming the inverse of the Laplace matrix - which needs to be accurate in order to obtain a solution to the right problem.

 template <class Matrix, class Preconditioner>
 void InverseMatrix<Matrix,Preconditioner>::vmult (Vector<double>       &dst,
                                                   const Vector<double> &src) const
 {
   SolverControl solver_control (src.size(), 1e-6*src.l2_norm());
   SolverCG<>    cg (solver_control);
 
   dst = 0;
 
   cg.solve (*matrix, dst, src, *preconditioner);
 }

The SchurComplement class template

This class implements the Schur complement discussed in the introduction. It is in analogy to step-20. Though, we now call it with a template parameter Preconditioner in order to access that when specifying the respective type of the inverse matrix class. As a consequence of the definition above, the declaration InverseMatrix now contains the second template parameter for a preconditioner class as above, which affects the SmartPointer object m_inverse as well.

 template <class Preconditioner>
 class SchurComplement : public Subscriptor
 {
   public:
     SchurComplement (const BlockSparseMatrix<double> &system_matrix,
                      const InverseMatrix<SparseMatrix<double>, Preconditioner> &A_inverse);
 
     void vmult (Vector<double>       &dst,
                 const Vector<double> &src) const;
 
   private:
     const SmartPointer<const BlockSparseMatrix<double> > system_matrix;
     const SmartPointer<const InverseMatrix<SparseMatrix<double>, Preconditioner> > A_inverse;
     
     mutable Vector<double> tmp1, tmp2;
 };
 
 
 
 template <class Preconditioner>
 SchurComplement<Preconditioner>::
 SchurComplement (const BlockSparseMatrix<double> &system_matrix,
                  const InverseMatrix<SparseMatrix<double>,Preconditioner> &A_inverse)
                 :
                 system_matrix (&system_matrix),
                 A_inverse (&A_inverse),
                 tmp1 (system_matrix.block(0,0).m()),
                 tmp2 (system_matrix.block(0,0).m())
 {}
 
 
 template <class Preconditioner>
 void SchurComplement<Preconditioner>::vmult (Vector<double>       &dst,
                                              const Vector<double> &src) const
 {
   system_matrix->block(0,1).vmult (tmp1, src);
   A_inverse->vmult (tmp2, tmp1);
   system_matrix->block(1,0).vmult (dst, tmp2);
 }

StokesProblem class implementation

StokesProblem::StokesProblem

The constructor of this class looks very similar to the one of step-20. The constructor initializes the variables for the polynomial degree, triangulation, finite element system and the dof handler. The underlying polynomial functions are of order degree+1 for the vector-valued velocity components and of order degree for the pressure. This gives the LBB-stable element pair $Q_{degree+1}^d\times Q_{degree}$, often referred to as the Taylor-Hood element.

Note that we initialize the triangulation with a MeshSmoothing argument, which ensures that the refinement of cells is done in a way that the approximation of the PDE solution remains well-behaved (problems arise if grids are too unstructered), see the documentation of Triangulation::MeshSmoothing for details.

 template <int dim>
 StokesProblem<dim>::StokesProblem (const unsigned int degree)
                 :
                 degree (degree),
                 triangulation (Triangulation<dim>::maximum_smoothing),
                 fe (FE_Q<dim>(degree+1), dim,
                     FE_Q<dim>(degree), 1),
                 dof_handler (triangulation)
 {}

StokesProblem::setup_dofs

Given a mesh, this function associates the degrees of freedom with it and creates the corresponding matrices and vectors. At the beginning it also releases the pointer to the preconditioner object (if the shared pointer pointed at anything at all at this point) since it will definitely not be needed any more after this point and will have to be re-computed after assembling the matrix, and unties the sparse matrix from its sparsity pattern object.

We then proceed with distributing degrees of freedom and renumbering them: In order to make the ILU preconditioner (in 3D) work efficiently, it is important to enumerate the degrees of freedom in such a way that it reduces the bandwidth of the matrix, or maybe more importantly: in such a way that the ILU is as close as possible to a real LU decomposition. On the other hand, we need to preserve the block structure of velocity and pressure already seen in in step-20 and step-21. This is done in two steps: First, all dofs are renumbered to improve the ILU and then we renumber once again by components. Since DoFRenumbering::component_wise does not touch the renumbering within the individual blocks, the basic renumbering from the first step remains. As for how the renumber degrees of freedom to improve the ILU: deal.II has a number of algorithms that attempt to find orderings to improve ILUs, or reduce the bandwidth of matrices, or optimize some other aspect. The DoFRenumbering namespace shows a comparison of the results we obtain with several of these algorithms based on the testcase discussed here in this tutorial program. Here, we will use the traditional Cuthill-McKee algorithm already used in some of the previous tutorial programs. In the section on improved ILU we're going to discuss this issue in more detail.

There is one more change compared to previous tutorial programs: There is no reason in sorting the dim velocity components individually. In fact, rather than first enumerating all $x$-velocities, then all $y$-velocities, etc, we would like to keep all velocities at the same location together and only separate between velocities (all components) and pressures. By default, this is not what the DoFRenumbering::component_wise function does: it treats each vector component separately; what we have to do is group several components into "blocks" and pass this block structure to that function. Consequently, we allocate a vector block_component with as many elements as there are components and describe all velocity components to correspond to block 0, while the pressure component will form block 1:

 template <int dim>
 void StokesProblem<dim>::setup_dofs ()
 {
   A_preconditioner.reset ();
   system_matrix.clear ();
   
   dof_handler.distribute_dofs (fe);  
   DoFRenumbering::Cuthill_McKee (dof_handler);
 
   std::vector<unsigned int> block_component (dim+1,0);
   block_component[dim] = 1;
   DoFRenumbering::component_wise (dof_handler, block_component);

Since we use adaptively refined grids the constraint matrix for hanging node constraints is generated from the DoF handler.

   hanging_node_constraints.clear ();
   DoFTools::make_hanging_node_constraints (dof_handler,
                                            hanging_node_constraints);
   hanging_node_constraints.close ();

In analogy to step-20, we count the dofs in the individual components. We could do this in the same way as there, but we want to operate on the block structure we used already for the renumbering: The function DoFTools::count_dofs_per_block does the same as DoFTools::count_dofs_per_component, but now grouped as velocity and pressure block via block_component.

   std::vector<unsigned int> dofs_per_block (2);
   DoFTools::count_dofs_per_block (dof_handler, dofs_per_block, block_component);  
   const unsigned int n_u = dofs_per_block[0],
                      n_p = dofs_per_block[1];
 
   std::cout << "   Number of active cells: "
             << triangulation.n_active_cells()
             << std::endl
             << "   Number of degrees of freedom: "
             << dof_handler.n_dofs()
             << " (" << n_u << '+' << n_p << ')'
             << std::endl;

The next task is to allocate a sparsity pattern for the system matrix we will create. We could do this in the same way as in step-20, i.e. directly build an object of type SparsityPattern through DoFTools::make_sparsity_pattern. However, there is a major reason not to do so: In 3D, the function DoFTools::max_couplings_between_dofs yields a conservative but rather large number for the coupling between the individual dofs, so that the memory initially provided for the creation of the sparsity pattern of the matrix is far too much -- so much actually that the initial sparsity pattern won't even fit into the physical memory of most systems already for moderately-sized 3D problems, see also the discussion in step-18. Instead, we first build a temporary object that uses a different data structure that doesn't require allocating more memory than necessary but isn't suitable for use as a basis of SparseMatrix or BlockSparseMatrix objects; in a second step we then copy this object into an object of BlockSparsityPattern. This is entirely analgous to what we already did in step-11 and step-18.

There is one snag again here, though: just as in step-27, it turns out that using the CompressedSparsityPattern (or the block version BlockCompressedSparsityPattern we would use here) has a bottleneck that makes the algorithm to build the sparsity pattern be quadratic in the number of degrees of freedom. This doesn't become noticable until we get well into the range of several 100,000 degrees of freedom, but eventually dominates the setup of the linear system when we get to more than a million degrees of freedom. This is due to the data structures used in the CompressedSparsityPattern class, nothing that can easily be changed. Fortunately, there is an easy solution, as already pointed out in step-27: the CompressedSetSparsityPattern class (and its block variant BlockCompressedSetSparsityPattern) has exactly the same interface, uses a different internal data structure, is slightly slower for smaller numbers of degrees of freedom (but there we don't care that much anyway) but is linear in the number of degrees of freedom and therefore much more efficient for large problems.

Consequently, this is the class that we will use for our intermediate sparsity representation. All this is done inside a new scope, which means that the memory of csp will be released once the information has been copied to sparsity_pattern.

   {
     BlockCompressedSetSparsityPattern csp (2,2);
 
     csp.block(0,0).reinit (n_u, n_u);
     csp.block(1,0).reinit (n_p, n_u);
     csp.block(0,1).reinit (n_u, n_p);
     csp.block(1,1).reinit (n_p, n_p);
   
     csp.collect_sizes();    
   
     DoFTools::make_sparsity_pattern (dof_handler, csp);
     hanging_node_constraints.condense (csp);
     sparsity_pattern.copy_from (csp);
   }

Finally, the system matrix, solution and right hand side are created from the block structure as in step-20:

   system_matrix.reinit (sparsity_pattern);
                                    
   solution.reinit (2);
   solution.block(0).reinit (n_u);
   solution.block(1).reinit (n_p);
   solution.collect_sizes ();
   
   system_rhs.reinit (2);
   system_rhs.block(0).reinit (n_u);
   system_rhs.block(1).reinit (n_p);
   system_rhs.collect_sizes ();
 }

StokesProblem::assemble_system

The assembly process follows the discussion in step-20 and in the introduction. We use the well-known abbreviations for the data structures that hold the local matrix, right hand side, and global numbers of the degrees of freedom for the present cell.

 template <int dim>
 void StokesProblem<dim>::assemble_system () 
 {
   system_matrix=0;
   system_rhs=0;
   
   QGauss<dim>   quadrature_formula(degree+2); 
   QGauss<dim-1> face_quadrature_formula(degree+2);
 
   FEValues<dim> fe_values (fe, quadrature_formula,
                            update_values    |
                            update_quadrature_points  |
                            update_JxW_values |
                            update_gradients);
   FEFaceValues<dim> fe_face_values (fe, face_quadrature_formula, 
                                     update_values    | 
                                     update_normal_vectors |
                                     update_quadrature_points  |
                                     update_JxW_values);
 
   const unsigned int   dofs_per_cell   = fe.dofs_per_cell;
   
   const unsigned int   n_q_points      = quadrature_formula.size();
 
   FullMatrix<double>   local_matrix (dofs_per_cell, dofs_per_cell);
   Vector<double>       local_rhs (dofs_per_cell);
 
   std::vector<unsigned int> local_dof_indices (dofs_per_cell);
   
   const RightHandSide<dim>          right_hand_side;
   std::vector<Vector<double> >      rhs_values (n_q_points,
                                                 Vector<double>(dim+1));

Next, we need two objects that work as extractors for the FEValues object. Their use is explained in detail in the report on Handling vector valued problems :

   const FEValuesExtractors::Vector velocities (0);
   const FEValuesExtractors::Scalar pressure (dim);

As an extension over step-20 and step-21, we include a few optimizations that make assembly faster for this particular problem. The improvements are based on the observation that we do a few calculations too many times when we do as in step-20: The symmetric gradient actually has dofs_per_cell different values per quadrature point, but we calculate it dofs_per_cell*dofs_per_cell times - for both the loop over i and the loop over j. So what we're going to do here is to avoid such double calculations by getting a vector of rank-2 tensors (and similarly for the divergence and the basis function value on pressure) at the quadrature point prior to starting the loop over the dofs on the cell. First, we create the respective objects that will hold the respective values. Then, we start the loop over all cells and the loop over the quadrature points, where we first extract these values. There is one more optimization we implement here: the local matrix (as well as the global one) is going to be symmetric, since the all the operations involved are symmetric with respect to $i$ and $j$. This is implemented by simply running the inner loop not to dofs_per_cell, but only up to i, the index of the outer loop.

   std::vector<SymmetricTensor<2,dim> > phi_grads_u (dofs_per_cell);
   std::vector<double>                  div_phi_u   (dofs_per_cell);
   std::vector<double>                  phi_p       (dofs_per_cell);
                                    
   typename DoFHandler<dim>::active_cell_iterator
     cell = dof_handler.begin_active(),
     endc = dof_handler.end();
   for (; cell!=endc; ++cell)
     { 
       fe_values.reinit (cell);
       local_matrix = 0;
       local_rhs = 0;
       
       right_hand_side.vector_value_list(fe_values.get_quadrature_points(),
                                         rhs_values);
       
       for (unsigned int q=0; q<n_q_points; ++q)
         {
           for (unsigned int k=0; k<dofs_per_cell; ++k)
             {
               phi_grads_u[k] = fe_values[velocities].symmetric_gradient (k, q);
               div_phi_u[k]   = fe_values[velocities].divergence (k, q);
               phi_p[k]       = fe_values[pressure].value (k, q);
             }
           for (unsigned int i=0; i<dofs_per_cell; ++i)
             {
               for (unsigned int j=0; j<=i; ++j)
                 {
                   local_matrix(i,j) += (phi_grads_u[i] * phi_grads_u[j]
                                         - div_phi_u[i] * phi_p[j]
                                         - phi_p[i] * div_phi_u[j]
                                         + phi_p[i] * phi_p[j])
                                        * fe_values.JxW(q);     
 
                 }
               const unsigned int component_i =
                 fe.system_to_component_index(i).first;
               local_rhs(i) += fe_values.shape_value(i,q) * 
                               rhs_values[q](component_i) *
                               fe_values.JxW(q);
             }
         }

Note that in the above computation of the local matrix contribution we added the term phi_p[i] * phi_p[j] , yielding a pressure mass matrix in the $(1,1)$ block of the matrix as discussed in the introduction. That this term only ends up in the $(1,1)$ block stems from the fact that both of the factors in phi_p[i] * phi_p[j] are only non-zero when all the other terms vanish (and the other way around).

Note also that operator* is overloaded for symmetric tensors, yielding the scalar product between the two tensors in the first line of the local matrix contribution.

The final step is, as usual, the transfer of the local contributions to the global system matrix. This works in the case of block vectors and matrices just the way we are used it to be, and also the terms constituting the pressure mass matrix are written into the correct position without any further interaction. We have to be careful about one thing, though. We have only build up half of the local matrix because of symmetry, but we're going to save the full system matrix in order to use the standard functions for solution. This is done by flipping the indices in case we are pointing into the empty part of the local matrix.

       cell->get_dof_indices (local_dof_indices);
 
       for (unsigned int i=0; i<dofs_per_cell; ++i)
         for (unsigned int j=0; j<dofs_per_cell; ++j)
           {
             if (j <= i)
               system_matrix.add (local_dof_indices[i],
                                             local_dof_indices[j],
                                  local_matrix(i,j));
             else
               system_matrix.add (local_dof_indices[i],
                                             local_dof_indices[j],
                                  local_matrix(j,i));
           }
       
       for (unsigned int i=0; i<dofs_per_cell; ++i)
         system_rhs(local_dof_indices[i]) += local_rhs(i);
     }

After the addition of the local contributions, we have to condense the hanging node constraints and interpolate Dirichlet boundary conditions. Further down below where we set up the mesh, we will associate the top boundary where we impose Dirichlet boundary conditions with boundary indicator 1. We will have to pass this boundary indicator as second argument to the function below interpolating boundary values. There is one more thing, though. The function describing the Dirichlet conditions was defined for all components, both velocity and pressure. However, the Dirichlet conditions are to be set for the velocity only. To this end, we use a component_mask that filters out the pressure component, so that the condensation is performed on velocity degrees of freedom only:

   hanging_node_constraints.condense (system_matrix);
   hanging_node_constraints.condense (system_rhs);  
 
   {
     std::map<unsigned int,double> boundary_values;
     std::vector<bool> component_mask (dim+1, true);
     component_mask[dim] = false;
     VectorTools::interpolate_boundary_values (dof_handler,
                                               1,
                                               BoundaryValues<dim>(),
                                               boundary_values,
                                               component_mask);
 
     MatrixTools::apply_boundary_values (boundary_values,
                                         system_matrix,
                                         solution,
                                         system_rhs);
   }

Before we're going to solve this linear system, we generate a preconditioner for the velocity-velocity matrix, i.e., block(0,0) in the system matrix. As mentioned above, this depends on the spatial dimension. Since the two classes described by the InnerPreconditioner<dim> :: type typedef have the same interface, we do not have to do anything different whether we want to use a sparse direct solver or an ILU:

   std::cout << "   Computing preconditioner..." << std::endl << std::flush;
       
   A_preconditioner
     = boost::shared_ptr<typename InnerPreconditioner<dim>::type>(new typename InnerPreconditioner<dim>::type());
   A_preconditioner->initialize (system_matrix.block(0,0),
                                 typename InnerPreconditioner<dim>::type::AdditionalData());
 
 }

StokesProblem::solve

After the discussion in the introduction and the definition of the respective classes above, the implementation of the solve function is rather straigt-forward and done in a similar way as in step-20. To start with, we need an object of the InverseMatrix class that represents the inverse of the matrix A. As described in the introduction, the inverse is generated with the help of an inner preconditioner of type InnerPreconditioner<dim>type.

 template <int dim>
 void StokesProblem<dim>::solve () 
 {
   const InverseMatrix<SparseMatrix<double>,
                       typename InnerPreconditioner<dim>::type>
     A_inverse (system_matrix.block(0,0), *A_preconditioner);
   Vector<double> tmp (solution.block(0).size());

This is as in step-20. We generate the right hand side $B A^{-1} F - G$ for the Schur complement and an object that represents the respective linear operation $B A^{-1} B^T$, now with a template parameter indicating the preconditioner - in accordance with the definition of the class.

   {
     Vector<double> schur_rhs (solution.block(1).size());
     A_inverse.vmult (tmp, system_rhs.block(0));
     system_matrix.block(1,0).vmult (schur_rhs, tmp);
     schur_rhs -= system_rhs.block(1);
   
     SchurComplement<typename InnerPreconditioner<dim>::type>
       schur_complement (system_matrix, A_inverse);

The usual control structures for the solver call are created...

     SolverControl solver_control (system_matrix.block(0,0).m(),
                                   1e-6*schur_rhs.l2_norm());
     SolverCG<>    cg (solver_control);

Now to the preconditioner to the Schur complement. As explained in the introduction, the preconditioning is done by a mass matrix in the pressure variable. It is stored in the $(1,1)$ block of the system matrix (that is not used anywhere else but in preconditioning).

Actually, the solver needs to have the preconditioner in the form $P^{-1}$, so we need to create an inverse operation. Once again, we use an object of the class InverseMatrix, which implements the vmult operation that is needed by the solver. In this case, we have to invert the pressure mass matrix. As it already turned out in earlier tutorial programs, the inversion of a mass matrix is a rather cheap and straight-forward operation (compared to, e.g., a Laplace matrix). The CG method with ILU preconditioning converges in 5-10 steps, independently on the mesh size. This is precisely what we do here: We choose another ILU preconditioner and take it along to the InverseMatrix object via the corresponding template parameter. A CG solver is then called within the vmult operation of the inverse matrix.

An alternative that is cheaper to build, but needs more iterations afterwards, would be to choose a SSOR preconditioner with factor 1.2. It needs about twice the number of iterations, but the costs for its generation are almost neglible.

     SparseILU<double> preconditioner;
     preconditioner.initialize (system_matrix.block(1,1), 
       SparseILU<double>::AdditionalData());
   
     InverseMatrix<SparseMatrix<double>,SparseILU<double> >
       m_inverse (system_matrix.block(1,1), preconditioner);

With the Schur complement and an efficient preconditioner at hand, we can solve the respective equation for the pressure (i.e. block 0 in the solution vector) in the usual way:

     cg.solve (schur_complement, solution.block(1), schur_rhs,
               m_inverse);

After this first solution step, the hanging node constraints have to be distributed to the solution in order to achieve a consistent pressure field.

     hanging_node_constraints.distribute (solution);
   
     std::cout << "  "
               << solver_control.last_step()
               << " outer CG Schur complement iterations for pressure"
               << std::flush
               << std::endl;    
   }

As in step-20, we finally need to solve for the velocity equation where we plug in the solution to the pressure equation. This involves only objects we already know - so we simply multiply $p$ by $B^T$, subtract the right hand side and multiply by the inverse of $A$. At the end, we need to distribute the constraints from hanging nodes in order to obtain a constistent flow field:

   {
     system_matrix.block(0,1).vmult (tmp, solution.block(1));
     tmp *= -1;
     tmp += system_rhs.block(0);
   
     A_inverse.vmult (solution.block(0), tmp);
 
     hanging_node_constraints.distribute (solution);
   }
 }

StokesProblem::output_results

The next function generates graphical output. In this example, we are going to use the VTK file format. We attach names to the individual variables in the problem: velocity to the dim components of velocity and pressure to the pressure.

Not all visualization programs have the ability to group individual vector components into a vector to provide vector plots; in particular, this holds for some VTK-based visualization programs. In this case, the logical grouping of components into vectors should already be described in the file containing the data. In other words, what we need to do is provide our output writers with a way to know which of the components of the finite element logically form a vector (with $d$ components in $d$ space dimensions) rather than letting them assume that we simply have a bunch of scalar fields. This is achieved using the members of the DataComponentInterpretation namespace: as with the filename, we create a vector in which the first dim components refer to the velocities and are given the tag DataComponentInterpretation::component_is_part_of_vector; we finally push one tag DataComponentInterpretation::component_is_scalar to describe the grouping of the pressure variable.

The rest of the function is then the same as in step-20.

 template <int dim>
 void
 StokesProblem<dim>::output_results (const unsigned int refinement_cycle)  const
 {
   std::vector<std::string> solution_names (dim, "velocity");
   solution_names.push_back ("pressure");
   
   std::vector<DataComponentInterpretation::DataComponentInterpretation>
     data_component_interpretation
     (dim, DataComponentInterpretation::component_is_part_of_vector);
   data_component_interpretation
     .push_back (DataComponentInterpretation::component_is_scalar);
       
   DataOut<dim> data_out;
   data_out.attach_dof_handler (dof_handler);  
   data_out.add_data_vector (solution, solution_names,
                             DataOut<dim>::type_dof_data,
                             data_component_interpretation);
   data_out.build_patches ();
   
   std::ostringstream filename;
   filename << "solution-"
            << Utilities::int_to_string (refinement_cycle, 2)
            << ".vtk";
 
   std::ofstream output (filename.str().c_str());
   data_out.write_vtk (output);
 }

StokesProblem::refine_mesh

This is the last interesting function of the StokesProblem class. As indicated by its name, it takes the solution to the problem and refines the mesh where this is needed. The procedure is the same as in the respective step in step-6, with the exception that we base the refinement only on the change in pressure, i.e., we call the Kelly error estimator with a mask object. Additionally, we do not coarsen the grid again:

 template <int dim>
 void
 StokesProblem<dim>::refine_mesh () 
 {
   Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
 
   std::vector<bool> component_mask (dim+1, false);
   component_mask[dim] = true;
   KellyErrorEstimator<dim>::estimate (dof_handler,
                                       QGauss<dim-1>(degree+1),
                                       typename FunctionMap<dim>::type(),
                                       solution,
                                       estimated_error_per_cell,
                                       component_mask);
 
   GridRefinement::refine_and_coarsen_fixed_number (triangulation,
                                                    estimated_error_per_cell,
                                                    0.3, 0.0);
   triangulation.execute_coarsening_and_refinement ();
 }

StokesProblem::run

The last step in the Stokes class is, as usual, the function that generates the initial grid and calls the other functions in the respective order.

We start off with a rectangle of size $4 \times 1$ (in 2d) or $4 \times 1 \times 1$ (in 3d), placed in $R^2/R^3$ as $(-2,2)\times(-1,0)$ or $(-2,2)\times(0,1)\times(-1,1)$, respectively. It is natural to start with equal mesh size in each direction, so we subdivide the initial rectangle four times in the first coordinate direction. To limit the scope of the variables involved in the creation of the mesh to the range where we actually need them, we put the entire block between a pair of braces:

 template <int dim>
 void StokesProblem<dim>::run () 
 {
   {
     std::vector<unsigned int> subdivisions (dim, 1);
     subdivisions[0] = 4;
 
     const Point<dim> bottom_left = (dim == 2 ?
                                     Point<dim>(-2,-1) :
                                     Point<dim>(-2,0,-1));
     const Point<dim> top_right   = (dim == 2 ?
                                     Point<dim>(2,0) :
                                     Point<dim>(2,1,0));
     
     GridGenerator::subdivided_hyper_rectangle (triangulation,
                                                subdivisions,
                                                bottom_left,
                                                top_right);
   }

A boundary indicator of 1 is set to all boundaries that are subject to Dirichlet boundary conditions, i.e. to faces that are located at 0 in the last coordinate direction. See the example description above for details.

   for (typename Triangulation<dim>::active_cell_iterator
          cell = triangulation.begin_active();
        cell != triangulation.end(); ++cell)
     for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
       if (cell->face(f)->center()[dim-1] == 0)
         cell->face(f)->set_all_boundary_indicators(1);

We then apply an initial refinement before solving for the first time. In 3D, there are going to be more degrees of freedom, so we refine less there:

   triangulation.refine_global (4-dim);

As first seen in step-6, we cycle over the different refinement levels and refine (except for the first cycle), setup the degrees of freedom and matrices, assemble, solve and create output:

   for (unsigned int refinement_cycle = 0; refinement_cycle<7;
        ++refinement_cycle)
     {
       std::cout << "Refinement cycle " << refinement_cycle << std::endl;
       
       if (refinement_cycle > 0)
         refine_mesh ();
       
       setup_dofs ();
 
       std::cout << "   Assembling..." << std::endl << std::flush;
       assemble_system ();      
 
       std::cout << "   Solving..." << std::flush;
       solve ();
       
       output_results (refinement_cycle);
 
       std::cout << std::endl;
     }
 }

The main function

The main function is the same as in step-20. We pass the element degree as a parameter and choose the space dimension at the well-known template slot.

 int main () 
 {
   try
     {
       deallog.depth_console (0);
 
       StokesProblem<3> flow_problem(1);
       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

Output of the program and graphical visualization

2D calculations

Running the program with the space dimension set to 2 in main() yields the following output:

examples/step-22> make run
============================ Remaking Makefile.dep
==============optimized===== step-22.cc
============================ Linking step-22
============================ Running step-22
Refinement cycle 0
   Number of active cells: 64
   Number of degrees of freedom: 679 (594+85)
   Assembling...
   Computing preconditioner...
   Solving...  11 outer CG Schur complement iterations for pressure

Refinement cycle 1
   Number of active cells: 160
   Number of degrees of freedom: 1683 (1482+201)
   Assembling...
   Computing preconditioner...
   Solving...  11 outer CG Schur complement iterations for pressure

Refinement cycle 2
   Number of active cells: 376
   Number of degrees of freedom: 3813 (3370+443)
   Assembling...
   Computing preconditioner...
   Solving...  11 outer CG Schur complement iterations for pressure

Refinement cycle 3
   Number of active cells: 880
   Number of degrees of freedom: 8723 (7722+1001)
   Assembling...
   Computing preconditioner...
   Solving...  11 outer CG Schur complement iterations for pressure

Refinement cycle 4
   Number of active cells: 2008
   Number of degrees of freedom: 19383 (17186+2197)
   Assembling...
   Computing preconditioner...
   Solving...  11 outer CG Schur complement iterations for pressure

Refinement cycle 5
   Number of active cells: 4288
   Number of degrees of freedom: 40855 (36250+4605)
   Assembling...
   Computing preconditioner...
   Solving...  11 outer CG Schur complement iterations for pressure

Refinement cycle 6
   Number of active cells: 8896
   Number of degrees of freedom: 83885 (74474+9411)
   Assembling...
   Computing preconditioner...
   Solving...  11 outer CG Schur complement iterations for pressure

The entire computation above takes about 30 seconds on a reasonably quick (for 2007 standards) machine.

What we see immediately from this is that the number of (outer) iterations does not increase as we refine the mesh. This confirms the statement in the introduction that preconditioning the Schur complement with the mass matrix indeed yields a matrix spectrally equivalent to the identity matrix (i.e. with eigenvalues bounded above and below independently of the mesh size or the relative sizes of cells). In other words, the mass matrix and the Schur complement are spectrally equivalent.

In the images below, we show the grids for the first six refinement steps in the program. Observe how the grid is refined in regions where the solution rapidly changes: On the upper boundary, we have Dirichlet boundary conditions that are -1 in the left half of the line and 1 in the right one, so there is an aprupt change at $x=0$. Likewise, there are changes from Dirichlet to Neumann data in the two upper corners, so there is need for refinement there as well:

step-22.2d.mesh-0.png

step-22.2d.mesh-1.png

step-22.2d.mesh-2.png

step-22.2d.mesh-3.png

step-22.2d.mesh-4.png

step-22.2d.mesh-5.png

Finally, following is a plot of the flow field. It shows fluid transported along with the moving upper boundary and being replaced by material coming from below:

step-22.2d.solution.png

This plot uses the capability of VTK-based visualization programs (in this case of VisIt) to show vector data; this is the result of us declaring the velocity components of the finite element in use to be a set of vector components, rather than independent scalar components in the StokesProblem<dim>output_results function of this tutorial program.

3D calculations

In 3d, the screen output of the program looks like this:

Refinement cycle 0
   Number of active cells: 32
   Number of degrees of freedom: 1356 (1275+81)
   Assembling...
   Computing preconditioner...
   Solving...  13 outer CG Schur complement iterations for pressure.

Refinement cycle 1
   Number of active cells: 144
   Number of degrees of freedom: 5088 (4827+261)
   Assembling...
   Computing preconditioner...
   Solving...  14 outer CG Schur complement iterations for pressure.

Refinement cycle 2
   Number of active cells: 704
   Number of degrees of freedom: 22406 (21351+1055)
   Assembling...
   Computing preconditioner...
   Solving...  14 outer CG Schur complement iterations for pressure.

Refinement cycle 3
   Number of active cells: 3168
   Number of degrees of freedom: 93176 (89043+4133)
   Assembling...
   Computing preconditioner...
   Solving...  15 outer CG Schur complement iterations for pressure.

Refinement cycle 4
   Number of active cells: 11456
   Number of degrees of freedom: 327808 (313659+14149)
   Assembling...
   Computing preconditioner...
   Solving...  15 outer CG Schur complement iterations for pressure.

Refinement cycle 5
   Number of active cells: 45056
   Number of degrees of freedom: 1254464 (1201371+53093)
   Assembling...
   Computing preconditioner...
   Solving...  14 outer CG Schur complement iterations for pressure.

Again, we see that the number of outer iterations does not increase as we refine the mesh. Nevertheless, the compute time increases significantly: for each of the iterations above separately, it takes a few seconds, a few seconds, 1min, 5min, 21min, and 1h35. (One level more can be run on machines with more memory, though one may have to replace the BlockCompressedSetSparsityPattern class by the BlockCompressedSparsityPattern class as the former needs significantly more memory and becomes the bottleneck; the latter class, however, has a superlinear runtime complexity.) This overall superlinear (in the number of unknowns) increase in runtime is due to the fact that our inner solver is not ${\cal O}(N)$: a simple experiment shows that as we keep refining the mesh, the average number of ILU-preconditioned CG iterations to invert the velocity-velocity block $A$ increases.

We will address the question of how possibly to improve our solver below.

As for the graphical output, the grids generated during the solution look as follow:

step-22.3d.mesh-0.png

step-22.3d.mesh-1.png

step-22.3d.mesh-2.png

step-22.3d.mesh-3.png

step-22.3d.mesh-4.png

step-22.3d.mesh-5.png

step-22.3d.mesh-6.png

Again, they show essentially the location of singularities introduced by boundary conditions. The vector field computed makes for an interesting graph:

step-22.3d.solution.png

The isocountours shown here as well are those of the pressure variable, showing the singularity at the point of discontinuous velocity boundary conditions.

Sparsity pattern

As explained during the generation of the sparsity pattern, it is important to have the numbering of degrees of freedom in mind when using preconditioners like incomplete LU decompositions. This is most conveniently visualized using the distribution of nonzero elements in the stiffness matrix.

If we don't do anything special to renumber degrees of freedom (i.e., without using DoFRenumbering::Cuthill_McKee, but with using DoFRenumbering::component_wise to ensure that degrees of freedom are appropriately sorted into their corresponding blocks of the matrix and vector, then we get the following image after the first adaptive refinement in two dimensions:

step-22.2d.sparsity-nor.png

In order to generate such a graph, you have to insert a piece of code like the following to the end of the setup step.

  {
    std::ofstream out ("sparsity_pattern.gpl");
    sparsity_pattern.print_gnuplot(out);
  }

It is clearly visible that the nonzero entries are spread over almost the whole matrix. This makes preconditioning by ILU inefficient: ILU generates a Gaussian elimination (LU decomposition) without fill-in elements, which means that more tentative fill-ins left out will result in a worse approximation of the complete decomposition.

In this program, we have thus chosen a more advanced renumbering of components. The renumbering with DoFRenumbering::Cuthill_McKee and grouping the components into velocity and pressure yields the following output:

step-22.2d.sparsity-ren.png

It is apparent that the situation has improved a lot. Most of the elements are now concentrated around the diagonal in the (0,0) block in the matrix. Similar effects are also visible for the other blocks. In this case, the ILU decomposition will be much closer to the full LU decomposition, which improves the quality of the preconditioner. (It may be interesting to note that the sparse direct solver UMFPACK does some internal renumbering of the equations before actually generating a sparse LU decomposition; that procedure leads to a very similar pattern to the one we got from the Cuthill-McKee algorithm.)

Finally, we want to have a closer look at a sparsity pattern in 3D. We show only the (0,0) block of the matrix, again after one adaptive refinement. Apart from the fact that the matrix size has increased, it is also visible that there are many more entries in the matrix. Moreover, even for the optimized renumbering, there will be a considerable amount of tentative fill-in elements. This illustrates why UMFPACK is not a good choice in 3D - a full decomposition needs many new entries that eventually won't fit into the physical memory (RAM):

step-22.3d.sparsity_uu-ren.png

Possible Extensions

Improved linear solver in 3D

We have seen in the section of computational results that the number of outer iterations does not depend on the mesh size, which is optimal in a sense of scalability. This does, however, not apply to the solver as a whole, as mentioned above: We did not look at the number of inner iterations when generating the inverse of the matrix $A$ and the mass matrix $M_p$. Of course, this is unproblematic in the 2D case where we precondition $A$ with a direct solver and the vmult operation of the inverse matrix structure will converge in one single CG step, but this changes in 3D where we only use an ILU preconditioner. There, the number of required preconditioned CG steps to invert $A$ increases as the mesh is refined, and each vmult operation involves on average approximately 14, 23, 36, 59, 75 and 101 inner CG iterations in the refinement steps shown above. (On the other hand, the number of iterations for applying the inverse pressure mass matrix is always around five, both in two and three dimensions.) To summarize, most work is spent on solving linear systems with the same matrix $A$ over and over again. What makes this look even worse is the fact that we actually invert a matrix that is about 95 precent the size of the total system matrix and stands for 85 precent of the non-zero entries in the sparsity pattern. Hence, the natural question is whether it is reasonable to solve a linear system with matrix $A$ for about 15 times when calculating the solution to the block system.

The answer is, of course, that we can do that in a few other (most of the time better) ways. Nevertheless, it has to be remarked that an indefinite system as the one at hand puts indeed much higher demands on the linear algebra than standard elliptic problems as we have seen in the early tutorial programs. The improvements are still rather unsatisfactory, if one compares with an elliptic problem of similar size.

Better ILU decomposition by smart