<a name="Intro"></a>
<h1>Introduction</h1>
+(This program was contributed by Martin Kronbichler and Wolfgang Bangerth.)
+
This program deals with the Stokes system of equations which reads as
follows in their non-dimensionalized form:
@f{eqnarray*}
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
-rankd-2 tensor of symmetrized gradients; a component-wise definition
+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)$.
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 interia effects then leads to the
+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.
\atop q}\right)$.
In practice, one wants to impose as little regularity on the pressure
-variable; consequently, we integrate by parts the second term:
+variable as possible; consequently, we integrate by parts the second term:
@f{eqnarray*}
(\mathrm v, -\textrm{div}\; \varepsilon(\textbf{u}))_{\Omega}
- (\textrm{div}\; \textbf{v}, p)_{\Omega}
=
(\textbf{v}, \textbf{f})_\Omega,
@f}
+We will deal with the boundary terms in the next section, but it is already
+clear from the domain terms
+@f{eqnarray*}
+ (\varepsilon(\mathrm v),\varepsilon(\textbf{u}))_{\Omega}
+ - (\textrm{div}\; \textbf{v}, p)_{\Omega}
+ -
+ (q,\textrm{div}\; \textbf{u})_{\Omega}
+@f}
+of the bilinear form that the Stokes equations yield a symmetric bilinear
+form, and consequently a symmetric (if indefinite) system matrix.
<h2>%Boundary conditions</h2>
-The weak form just presented immediately presents us with different
+The weak form just derived immediately presents us with different
possibilities for imposing boundary conditions:
<ol>
<li>Dirichlet velocity boundary conditions: On a part
In other words, as usual, strongly imposed boundary values do not
appear in the weak form.
-<li>Neumann-type boundary conditions: On the rest of the boundary
+ It is noteworthy that if we impose Dirichlet boundary values on the entire
+ boundary, then the pressure is only determined up to a constant.
+
+<li>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:
@f{eqnarray*}
-
(q,\textrm{div}\; \textbf{u})_{\Omega}
+
- (\textbf S \textbf u, \textbf{v})_{\Gamma_N}
+ (\textbf S \textbf u, \textbf{v})_{\Gamma_R}
=
(\textbf{v}, \textbf{f})_\Omega.
@f}
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
+ <code>dim</code>-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:
@f{eqnarray*}
\textbf u_{\textbf t} &=& 0,
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,Q$ also satisfy the LBB
+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 @ref step_20 "step-20".
For the Stokes equations, there are a number of possible choices to ensure
-
(\textbf{v}_h, \textbf g_N)_{\Gamma_N}
@f}
-for all test functions $\textbf v_h, q_h$.
+for all test functions $\textbf v_h, q_h$. Assembling the linear system
+associated with this problem follows the same lines used in @ref step_20
+"step-20", @ref step_21 "step-21", and explained in detail in the @ref
+vector_valued module.
\left(\begin{array}{cc}
A & B^T \\ B & 0
\end{array}\right)
- \left(\begin{array}{cc}
+ \left(\begin{array}{c}
U \\ P
\end{array}\right)
=
- \left(\begin{array}{cc}
+ \left(\begin{array}{c}
F \\ G
\end{array}\right),
@f}
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} 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.
+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.
<li>
Second, every time we multiplied with $B^TM^{-1}B$ we had to solve with the
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} \eta \nabla^s$ (where $\nabla^s$ is the symmetric gradient)
+$-\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
<ol>
<li>
-First, it makes the outer preconditioner simple: the
-Schur complement corresponds to the operator $-\textrm{div} (-\textrm{div}
-\eta \nabla^s)^{-1} \nabla$ on the pressure space; forgetting about the
-viscosity $\eta$ and ignoring 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. It turns out that it isn't
-easy to solve this Schur complement in a straight forward 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 <a
-href="#Results">results section</a> of this program to see that indeed
-the number of CG iterations does not change as we refine the mesh).
+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, however:
+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 <a href="#Results">results section</a> 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
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)$ where
-$\nabla^s$ is the symmetrized gradient of a vector field). Thus,
+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$
href="../../readme.html#optional-software">ReadMe file</a> 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 factor can be computed in a matter of a few
-seconds.
+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
unavoidable.
</ol>
-
-<h2>Implementation</h2>
-
-Note anything here that would be particular about this program.
+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:
+@code
+template <int dim>
+struct InnerPreconditioner;
+
+template <>
+struct InnerPreconditioner<2>
+{
+ typedef SparseDirectUMFPACK type;
+};
+
+template <>
+struct InnerPreconditioner<3>
+{
+ typedef SparseILU<double> type;
+};
+@endcode
+
+From hereon, we can refer to the type <code>typename
+InnerPreconditioner@<dim@>::type</code> 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.
<h2>The testcase</h2>
-Describe the testcase we are dealing with here.
+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 from below
+with magma. 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]$:
+@f{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,
+@f}
+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 <a href="#Results">results section</a>, 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.
+
projects / dealii-svn.git / commitdiff