# Diff of /trunk/doc/user/stokessolver.tex

revision 2793 by gross, Tue Dec 1 06:10:10 2009 UTC revision 2843 by gross, Thu Jan 14 05:51:28 2010 UTC
# Line 4  In this section we discuss how to solve Line 4  In this section we discuss how to solve
4
5  We want to calculate the velocity \index{velocity} field $v$ and pressure $p$ of an incompressible fluid \index{incompressible fluid}. They are given as the solution of the Stokes problem\index{Stokes problem}  We want to calculate the velocity \index{velocity} field $v$ and pressure $p$ of an incompressible fluid \index{incompressible fluid}. They are given as the solution of the Stokes problem\index{Stokes problem}
6  \label{Stokes 1}  \label{Stokes 1}
7  -\left(\eta(v\hackscore{i,j}+ v\hackscore{i,j})\right)\hackscore{,j}+p\hackscore{,i}=f\hackscore{i}-\sigma\hackscore{ij,j}  -\left(\eta(v\hackscore{i,j}+ v\hackscore{j,i})\right)\hackscore{,j}+p\hackscore{,i}=f\hackscore{i}-\sigma\hackscore{ij,j}
8
9  where  $f\hackscore{i}$ defines an internal force \index{force, internal} and $\sigma\hackscore{ij}$ is an intial stress \index{stress, initial}. The viscosity $\eta$ may weakly depend on pressure and velocity. If relevant we will use the notation $\eta(v,p)$ to express this dependency.  where  $f\hackscore{i}$ defines an internal force \index{force, internal} and $\sigma\hackscore{ij}$ is an intial stress \index{stress, initial}. The viscosity $\eta$ may weakly depend on pressure and velocity. If relevant we will use the notation $\eta(v,p)$ to express this dependency.
10
# Line 14  We assume an incompressible media: Line 14  We assume an incompressible media:
14
15  Natural boundary conditions are taken in the form  Natural boundary conditions are taken in the form
16  \label{Stokes Boundary}  \label{Stokes Boundary}
17  \left(\eta(v\hackscore{i,j}+ v\hackscore{i,j})\right)n\hackscore{j}-n\hackscore{i}p=s\hackscore{i} - \alpha \cdot n\hackscore{i} n\hackscore{j} v\hackscore{j}+\sigma\hackscore{ij} n\hackscore{j}  \left(\eta(v\hackscore{i,j}+ v\hackscore{j,i})\right)n\hackscore{j}-n\hackscore{i}p=s\hackscore{i} - \alpha \cdot n\hackscore{i} n\hackscore{j} v\hackscore{j}+\sigma\hackscore{ij} n\hackscore{j}
18
19  which can be overwritten by constraints of the form  which can be overwritten by constraints of the form
20  \label{Stokes Boundary0}  \label{Stokes Boundary0}
# Line 41  p \\ Line 41  p \\
41  G \\  G \\
42  0 \\  0 \\
43  \end{array} \right]  \end{array} \right]
45
46  where $A$ is coercive, self-adjoint linear operator in a suitable Hilbert space, $B$ is the $(-1) \cdot$ divergence operator and $B^{*}$ is it adjoint operator (=gradient operator).  where $A$ is coercive, self-adjoint linear operator in a suitable Hilbert space, $B$ is the $(-1) \cdot$ divergence operator and $B^{*}$ is it adjoint operator (=gradient operator).
47  For more details on the mathematics see references \cite{AAMIRBERKYAN2008,MBENZI2005}.  For more details on the mathematics see references \cite{AAMIRBERKYAN2008,MBENZI2005}.
48
49  If $v\hackscore{0}$ and $p\hackscore{0}$ are given initial guesses for  If $v\hackscore{0}$ and $p\hackscore{0}$ are given initial guesses for
50  velocity and pressure we calculate a correction $dv$ for the velocity by solving the first  velocity and pressure we calculate a correction $dv$ for the velocity by solving the first
51  equation of equation~\ref{SADDLEPOINT}  equation of equation~\ref{STOKES}
52   \label{SADDLEPOINT ITER STEP 1}   \label{STOKES ITER STEP 1}
53   A dv\hackscore{1} = G - A v\hackscore{0} - B^{*} p\hackscore{0}   A dv\hackscore{1} = G - A v\hackscore{0} - B^{*} p\hackscore{0}
54
55  We then insert the new approximation $v\hackscore{1}=v\hackscore{0}+dv\hackscore{1}$ to calculate a correction $dp\hackscore{2}$  We then insert the new approximation $v\hackscore{1}=v\hackscore{0}+dv\hackscore{1}$ to calculate a correction $dp\hackscore{2}$
# Line 59  for the pressure and an additional corre Line 59  for the pressure and an additional corre
59   B A^{-1} B^{*} dp\hackscore{2} & = & Bv\hackscore{1} \\   B A^{-1} B^{*} dp\hackscore{2} & = & Bv\hackscore{1} \\
60   A dv\hackscore{2} & = & B^{*} dp\hackscore{2}   A dv\hackscore{2} & = & B^{*} dp\hackscore{2}
61  \end{array}  \end{array}
62   \label{SADDLEPOINT ITER STEP 2}   \label{STOKES ITER STEP 2}
63
64  The new velocity and pressure are then given by $v=v\hackscore{1}-dv\hackscore{2}$ and  The new velocity and pressure are then given by $v\hackscore{2}=v\hackscore{1}-dv\hackscore{2}$ and
65  $p=p\hackscore{0}+dp\hackscore{2}$ which will fullfill the block system~\ref{SADDLEPOINT}.  $p\hackscore{2}=p\hackscore{0}+dp\hackscore{2}$ which will fullfill the block system~\ref{STOKES}.
66  This solution strategy is called the Uzawa scheme \index{Uzawa scheme}. There are  This solution strategy is called the Uzawa scheme \index{Uzawa scheme}.
67  two problems with this scheme. Firstly, in practice we will use an iterative scheme
68  to solve the problem for operator $A$. So we will be unable to calculate the operator  There is a problem with this scheme: In practice we will use an iterative scheme
69  $B A^{-1} B^{*}$ required for $dp\hackscore{2}$. In fact, we need to use an iterative scheme  to solve any problem for operator $A$. So we will be unable to calculate the operator
70  to solve the first equation in~\ref{SADDLEPOINT ITER STEP 2} where in each iteration step  $B A^{-1} B^{*}$ required for $dp\hackscore{2}$ explicitly. In fact, we need to use another
71  an iterative solver for $A$ is applied. The second issue is that  iterative scheme to solve the first equation in~\ref{STOKES ITER STEP 2} where in each iteration step
72  viscosity $\eta$ may depend of velocity or pressure and so we need to iterate over the  an iterative solver for $A$ is applied. Another issue is the fact that the
73  three equations~\ref{SADDLEPOINT ITER STEP 1} and~\ref{SADDLEPOINT ITER STEP 2}. Using the  viscosity $\eta$ may depend on velocity or pressure and so we need to iterate over the
74  two norms  three equations~\ref{STOKES ITER STEP 1} and~\ref{STOKES ITER STEP 2}.
75
76    In the following we will use the two norms
77
78  \|v\|\hackscore{1}^2 = \int\hackscore{\Omega} v\hackscore{j,k}v\hackscore{j,k} \; dx  \|v\|\hackscore{1}^2 = \int\hackscore{\Omega} v\hackscore{j,k}v\hackscore{j,k} \; dx
79  \mbox{ and }  \mbox{ and }
80  \|p\|\hackscore{0}^2= \int\hackscore{\Omega} p^2 \; dx.  \|p\|\hackscore{0}^2= \int\hackscore{\Omega} p^2 \; dx.
81  \label{STOKES STOP}  \label{STOKES STOP}
82
83  to define the stopping criterium  for velicity $v$ and pressure $p$. The iteration is terminated if the stopping criterium
84     \label{STOKES STOPPING CRITERIA}
85  \max(\|Bv\hackscore{1}\|\hackscore{0},\|v-v\hackscore{0}\|\hackscore{1}) \le \tau \cdot \|v\|\hackscore{1}  \max(\|Bv\hackscore{1}\|\hackscore{0},\|v\hackscore{2}-v\hackscore{0}\|\hackscore{1}) \le \tau \cdot \|v\hackscore{2}\|\hackscore{1}
86
87  to terminate the overall iteration with a overall given tolerance $0<\tau<1$.   for a given  given tolerance $0<\tau<1$ is meet. Notice that because of the first equation of~\ref{STOKES ITER STEP 2} we have that
88  Notice that because of the first equation of  $\|Bv\hackscore{1}\|\hackscore{0}$ equals the
89  and~\ref{SADDLEPOINT ITER STEP 2} $\|Bv\hackscore{1}\|\hackscore{0}$ is equal to the  norm of $B A^{-1} B^{*} dp\hackscore{2}$ and consequently provides a norm for the pressure correction.
90  norm of $B A^{-1} B^{*} dp\hackscore{2}$ and consequently provides a norm for the  pressure correction.
91    We want to optimize the tolerance choice for solving~\ref{STOKES ITER STEP 1}
92  Let us first have a look at the solution process~\ref{SADDLEPOINT ITER STEP 1} which needs to resolve the  and~\ref{STOKES ITER STEP 2}. To do this we write the iteration scheme as a fixed point problem. Here
93  non-linearity of the viscosity coeficient. We will run a iterative process by resolving  we consider the errors produced by the iterative solvers beeing used.
94  the equation~\ref{SADDLEPOINT ITER STEP 1} with $v\hackscore{0} \leftarrow v\hackscore{1}$. The first observation is that there is no point start  From Equation~\ref{STOKES ITER STEP 1} we have
95  the iteration process over the second set of equation if $\|Bv\hackscore{1}\|\hackscore{0}$   \label{STOKES total V1}
96  is smaller then $\|v\hackscore{1}-v\hackscore{0}\|\hackscore{1}$. So we will terminate the non-linear iteration  v\hackscore{1} = e\hackscore{1} + v\hackscore{0} + A^{-1} ( G - Av\hackscore{0} - B^{*} p\hackscore{0} )
process of the first equation if

\theta \cdot \|v\hackscore{1}-v\hackscore{0}\| \le \|Bv\hackscore{1}\|\hackscore{0}
\mbox{ or }
\|v\hackscore{1}-v\hackscore{0}\|\hackscore{1} \le \tau \cdot \|v\hackscore{1}\|\hackscore{1}

where $0<\theta<1$ is a given factor (typically $\theta=0.5$).

We need to think about appropriate stopping criteria when solving
\ref{SADDLEPOINT ITER STEP 1}. We would like to use very weak convergence criteria to reduce computational costs but they need to be tight enough to not interfere with the
convergence of the iteration process one level above.
We will use a sparse matrix solver so we have not full control on the norm $\|.\|\hackscore{s}$ used in the stopping criteria

\| G - A v\hackscore{1} - B^{*} p\hackscore{0}  \|\hackscore{s} \le \tau\hackscore{1} \| G - A v\hackscore{0}  - B^{*} p\hackscore{0} \|\hackscore{s}
97
98  If $e\hackscore{1}$ is the error of the returned solution $dv\hackscore{1}$ this translates into the condition  where $e\hackscore{1}$ is the error when solving~\ref{STOKES ITER STEP 1}.
99    We will use a sparse matrix solver so we have not full control on the norm $\|.\|\hackscore{s}$ used in the stopping criteria for this equation. In fact we will have a stopping criteria of the form
100
101  \| e\hackscore{1} \|\hackscore{1} \le K \tau\hackscore{1} \| dv\hackscore{1} - e\hackscore{1} \|\hackscore{1}  \| A e\hackscore{1} \|\hackscore{s}  = \| G - A v\hackscore{1} - B^{*} p\hackscore{0} \|\hackscore{s} \le \tau\hackscore{1} \| G - A v\hackscore{0} - B^{*} p\hackscore{0} \|\hackscore{s}
102
103  The constant $K$ represents some uncertainty combining a variety of unknown factors such as the  where $\tau\hackscore{1}$ is the tolerance which we need to choose. This translates into the condition
solver being used and the condition number of the stiffness matrix.  This leads to the estimate

\| dv\hackscore{1} \|\hackscore{1} \le (1+K \tau\hackscore{1}) \chi \| dv\hackscore{1}^- \|\hackscore{1}
\label{DV COND}

where $chi$ is the convergence rate of the \ref{SADDLEPOINT ITER STEP 1} iteration and $dv\hackscore{1}^-$ the correction from the last step. In order to make sure that the error $e\hackscore{1}$ can be neglected we need
to have $K \tau\hackscore{1} \le \chi$ which leads to the condition

\tau\hackscore{1} \le \frac{\min(\chi,\frac{1}{2})}{K}
\label{NEW TAU 1}

We can estimate $K$ by comparing estimates for the convergence rate $\chi=\frac{\| dv\hackscore{1} \|\hackscore{1}}{\| dv\hackscore{1}^- \|\hackscore{1}}$ by checking our error prediction~\ref{DV COND}. If we have
access to the estimate of the convergence rate $\chi^-$ of the previous iteration step we
get

\frac{1}{K} = \frac{\chi^-}{\chi-\chi^{-}} \tau\hackscore{1}^-

if $\chi \ge \chi^- (1 + \chi^-)$ which means that our prediction for a suitable tolerane based on the avaibale $K$ value was wrong.
So we start with the value $K=1$ in~\ref{NEW TAU 1} and recalculate
$K$ if we have underestimate $K$.

If $dv\hackscore{1}$
From Equation~\ref{SADDLEPOINT ITER STEP 1} we have

v\hackscore{1} = e\hackscore{1} + A^{-1} ( G - B^{*} p )

This translates into the conditoon
104
105  \| e\hackscore{1} \|\hackscore{1} \le K \tau\hackscore{1} \| dv\hackscore{1} - e\hackscore{1} \|\hackscore{1}  \| e\hackscore{1} \|\hackscore{1} \le K \tau\hackscore{1} \| dv\hackscore{1} - e\hackscore{1} \|\hackscore{1}
\mbox{ therefore }
\| e\hackscore{1} \|\hackscore{1} \le \frac{K \tau\hackscore{1}}{1-K \tau\hackscore{1}}  \| dv\hackscore{1} \|\hackscore{1}
106
107  The constant $K$ represents some uncertainty combining a variety of unknown factors such as the  The constant $K$ represents some uncertainty combining a variety of unknown factors such as the
108  solver being used and the condition number of the stiffness matrix.    norm being used and the condition number of the stiffness matrix.
109    From the first equation of~\ref{STOKES ITER STEP 2} we have
110    \label{STOKES total P2}
111    p\hackscore{2} =  p\hackscore{0} + (B A^{-1} B^{*})^{-1} (e\hackscore{2} + Bv\hackscore{1} )
112
113    where $e\hackscore{2}$ represents the error when solving~\ref{STOKES ITER STEP 2}.
114  ===================================  We use an iterative preconditioned conjugate gradient method (PCG) \index{linear solver!PCG}\index{PCG} with iteration
115    operator $B A^{-1} B^{*}$ using the $\|.\|\hackscore{0}$ norm. As suitable preconditioner \index{preconditioner} for the iteration
116    operator we use $\frac{1}{\eta}$ \cite{ELMAN}, ie
We use iterative techniques to solve this problem: Given an approximation $v$ and $p$ for
velocity and pressure we perform the following steps in the Uzawa scheme \index{Uzawa scheme} style:
\begin{enumerate}
\item calculate viscosity $\eta(v,p)$ and assemble operator $A$ from $\eta$.
\item Solve for $dv$:

A dv = G - A v - B^{*} p  \label{SADDLEPOINT ITER STEP 1}

\item update $v\hackscore{1}= v+ dv$
\item if $\max(\|Bv\hackscore{1}\|\hackscore{0},\|dv\|\hackscore{1}) \le \tau \cdot \|v\hackscore{1}\|\hackscore{1}$: return v$\hackscore{1},p$
\item Solve for $dp$:

\begin{array}{rcl}
B A^{-1} B^{*} dp & = & Bv\hackscore{1} \\
A dv\hackscore{2} & = & B^{*} dp
\end{array}

\item update $p\hackscore{2}\leftarrow p+ dp$ and $v\hackscore{2}= v\hackscore{1} - dv\hackscore{2}$
\item goto Step 0 with $v=v\hackscore{2}$ and $p=p\hackscore{2}$.
\end{enumerate}
where $\tau$ is the given tolerance. Notice that the operator $A$ if it is depending on $v$ or $p$ is updated once only.
In this algorithm

\|v\|\hackscore{1}^2 = \int\hackscore{\Omega} v\hackscore{j,k}v\hackscore{j,k} \; dx
\mbox{ and }
\|p\|\hackscore{0}^2= \int\hackscore{\Omega} p^2 \; dx.
\label{STOKES STOP}

so the stopping criterion used is checking for convergence as well as for
the fact that the incompressiblity condition~\ref{Stokes 2} is fullfilled.

To solve the update step~\ref{SADDLEPOINT ITER STEP 2} we use an iterative solver with iteration
operator $B A^{-1} B^{*}$, eg. using generalized minimal residual method (GMRES) \index{linear solver!GMRES}\index{GMRES}, preconditioned conjugate gradient method (PCG) \index{linear solver!PCG}\index{PCG}. As suitable preconditioner \index{preconditioner} for this operator is $\frac{1}{\eta}$, ie
117  the evaluation of the preconditioner $P$ for a given pressure increment $q$ is the solution of  the evaluation of the preconditioner $P$ for a given pressure increment $q$ is the solution of
118   \label{P PREC}   \label{STOKES P PREC}
119  \frac{1}{\eta} Pq = q \; .  \frac{1}{\eta} (Pq) = q \; .
120
121  Note that in each evaluation of the iteration operator $q=B A^{-1} B^{*} s$ one needs to solve  Note that in each evaluation of the iteration operator $q=B A^{-1} B^{*} s$ one needs to solve
122  the problem  the problem
123   \label{P OPERATOR}   \label{STOKES P OPERATOR}
124  A w = B^{*} s  A w = B^{*} s
125
126  with sufficient accuracy to return $q=Bw$. Notice that the residual $r$ is given as  with sufficient accuracy to return $q=Bw$. We assume that the choicen tolerance is
127    sufficiently small, for instane one can take $\tau\hackscore{2}^2$
128    where $\tau\hackscore{2}$ is the tolerance for~\ref{STOKES ITER STEP 2}.
129
130    In an implementation we use the fact that the residual $r$ is given as
131   \label{STOKES RES }   \label{STOKES RES }
132   r= B (v\hackscore{1} -  A^{-1} B^{*} dp) =  B (v\hackscore{1} - A^{-1} B^{*} dp) = B (v\hackscore{1}-dv\hackscore{2}) = B v\hackscore{2}   r= B (v\hackscore{1} -  A^{-1} B^{*} dp) =  B (v\hackscore{1} - A^{-1} B^{*} dp) = B (v\hackscore{1}-dv\hackscore{2}) = B v\hackscore{2}
133
134  so in fact the residual $r$ is represented by the updated velocity $v\hackscore{2}$. This saves the recovery of  In particular we have $e\hackscore{2} = B v\hackscore{2}$
135  $dv\hackscore{2}$ in~\ref{SADDLEPOINT ITER STEP 2} after $dp$ has been calculated as iterative method such as PCG calculate the solution approximations along with the their residual. In PCG the iteration is terminated if  So the residual $r$ is represented by the updated velocity $v\hackscore{2}=v\hackscore{1}-dv\hackscore{2}$. In practice, if
136   \label{P OPERATOR}  one uses the velocity to represent the residual $r$ there is no need
137    to recover $dv\hackscore{2}$ in~\ref{STOKES ITER STEP 2} after $dp\hackscore{2}$ has been calculated.
138    In PCG the iteration is terminated if
139     \label{STOKES P OPERATOR ERROR}
140  \| P^{\frac{1}{2}}B v\hackscore{2} \|\hackscore{0} \le \tau\hackscore{2} \| P^{\frac{1}{2}}B v\hackscore{1} \|\hackscore{0}  \| P^{\frac{1}{2}}B v\hackscore{2} \|\hackscore{0} \le \tau\hackscore{2} \| P^{\frac{1}{2}}B v\hackscore{1} \|\hackscore{0}
141
142  where $\tau\hackscore{2}$ is the given tolerane. The update step~\ref{P OPERATOR} involves the  where $\tau\hackscore{2}$ is the given tolerance. This translates into
143  solution of a sparse matrix problem. We use the tolerance $\tau\hackscore{2}^2$ in order to make sure that any   \label{STOKES P OPERATOR ERROR 2}
144  error from solving this problem does not interfere with the PCG iteration.  \|e\hackscore{2}\|\hackscore{0} = \| B v\hackscore{2} \|\hackscore{0} \le M \tau\hackscore{2} \| B v\hackscore{1} \|\hackscore{0}

We need to think about appropriate stopping criteria when solving
\ref{SADDLEPOINT ITER STEP 1}, \ref{SADDLEPOINT ITER STEP 2} and~\ref{P OPERATOR}. We would like to use very weak convergence criteria to reduce computational costs but they need to be tight enough to not interfere with the
convergence of the iteration process one level above. From Equation~\ref{SADDLEPOINT ITER STEP 1} we have

v\hackscore{1} = e\hackscore{1} + A^{-1} ( G - B^{*} p )

We will use a sparse matrix solver so we have not full control on the norm $\|.\|\hackscore{s}$ used in the stopping criteria

\| G - A v\hackscore{1} - B^{*} p \|\hackscore{s} \le \tau\hackscore{1} \| G - A v - B^{*} p \|\hackscore{s}
145
146  This translates into the conditoon  where $M$ is taking care of the fact that $P^{\frac{1}{2}}$ is dropped.

\| e\hackscore{1} \|\hackscore{1} \le K \tau\hackscore{1} \| dv\hackscore{1} - e\hackscore{1} \|\hackscore{1}
\mbox{ therefore }
\| e\hackscore{1} \|\hackscore{1} \le \frac{K \tau\hackscore{1}}{1-K \tau\hackscore{1}}  \| dv\hackscore{1} \|\hackscore{1}

The constant $K$ represents some uncertainty combining a variety of unknown factors such as the
solver being used and the condition number of the stiffness matrix.
147
148  From the first equation of~\ref{SADDLEPOINT ITER STEP 2} we have  As we assume that there is no significant error from solving with the operator $A$ we have
149     \label{STOKES total V2}
p\hackscore{2} =  p + (B A^{-1} B^{*})^{-1} (e\hackscore{2} + Bv\hackscore{1} ) =
(B A^{-1} B^{*})^{-1} ( e\hackscore{2} + B e\hackscore{1} + B A^{-1} G)

and simlar

150  v\hackscore{2} =  v\hackscore{1} - dv\hackscore{2}  v\hackscore{2} =  v\hackscore{1} - dv\hackscore{2}
151  = v\hackscore{1}  - A^{-1} B^{*}dp  = v\hackscore{1}  - A^{-1} B^{*}dp
= e\hackscore{1}  + A^{-1} ( G - B^{*} p ) - A^{-1} B^{*} (p\hackscore{2}-p)
= e\hackscore{1}  + A^{-1} ( G - B^{*} p\hackscore{2})
152
153  This shows that we can write the iterative process as a fixed point iteration to solve (assume all errors are zero)  Combining the equations~\ref{STOKES total V1},~\ref{STOKES total P2} and~\ref{STOKES total V2} and
154    setting the errors to zero we can write the solution process as a fix point problem
155
156  v = \Phi(v,p) \mbox{ and } p = \Psi(u,p)  v = \Phi(v,p) \mbox{ and } p = \Psi(u,p)
157
158  where  with suitable functions $\Phi(v,p)$ and $\Psi(v,p)$ representing the itartion operator without
159    errors. In fact for a linear problem,  $\Phi$ and $\Psi$ are constant. With this notation we can write
160   \begin{array}{rcl}  the update step in the form $p\hackscore{2}= \delta p + \Psi(v\hackscore{0},p\hackscore{0})$ and
161  \Psi(v,p) & = &  (B A^{-1} B^{*})^{-1} B A^{-1} G \\  $v\hackscore{2}= \delta v + \Phi(v\hackscore{0},p\hackscore{0})$ where
162  \Phi(u,p) & = & A^{-1} ( G - B^{*} (B A^{-1} B^{*})^{-1} B A^{-1} G )  the total error $\delta p$ and $\delta v$ are given as
\end{array}

Notice that if $A$ is independent from $v$ and $p$ the operators $\Phi(v,p)$ and $\Psi(u,p)$ are constant
and threfore the iteration will terminate - providing no termination errors in sub-iterations - after one step.
We also can give a formula for the error
163
164   \begin{array}{rcl}   \begin{array}{rcl}
165  \delta p & = &  (B A^{-1} B^{*})^{-1} ( e\hackscore{2} + B e\hackscore{1} ) \\  \delta p & = &  (B A^{-1} B^{*})^{-1} ( e\hackscore{2} + B e\hackscore{1} ) \\
166  \delta v & = &   e\hackscore{1} -  A^{-1} B^{*}\delta p  \;.  \delta v & = &  e\hackscore{1} -  A^{-1} B^{*}\delta p  \;.
167  \end{array}\label{STOKES ERRORS}  \end{array}\label{STOKES ERRORS}
168
169  Notice that $B\delta v = - e\hackscore{2}=-Bv\hackscore{2}$.  Notice that $B\delta v = - e\hackscore{2}=-Bv\hackscore{2}$. Our task is now to choose the tolerances
170    $\tau\hackscore{1}$ and $\tau\hackscore{2}$ such that the global errors $\delta p$ and $\delta v$
171  With this notation  do not stop the convergence of the iteration process.
172
173   \begin{array}{rcl}  To measure convergence we use
174  v\hackscore{2} = \Phi(v,p) + \delta v \\
175  p\hackscore{2} = \Psi(v,p) + \delta p \\  \epsilon = \max(\|v\hackscore{2}-v\|\hackscore{1}, \|B A^{-1} B^{*} (p\hackscore{2}-p)\|\hackscore{0})
176
177    In practice using the fact that $B A^{-1} B^{*} (p\hackscore{2}-p\hackscore{0}) = B v\hackscore{1}$
178    and assuming that $v\hackscore{2}$ gives a better approximation to the true $v$ than
179    $v\hackscore{0}$ we will use the estimate
180
181    \epsilon = \max(\|v\hackscore{2}-v\hackscore{0}\|\hackscore{1}, \|B v\hackscore{1}\|\hackscore{0})
182
183    to estimate the progress of the iteration step after the step is completed.
184    Note that the estimate of $\epsilon$
185    used in the stopping criteria~\ref{STOKES STOPPING CRITERIA}. If $\chi^{-}$ is the convergence rate assuming
186    exact calculations, i.e. $e\hackscore{1}=0$ and $e\hackscore{2}=0$, we are expecting
187    to maintain $\epsilon \le \chi^{-} \cdot \epsilon^{-}$. For the
188    pressure increment we get:
189     \label{STOKES EST 1}
190    \begin{array}{rcl}
191    \|B A^{-1} B^{*} (p\hackscore{2}-p)\|\hackscore{0}
192     & \le & \|B A^{-1} B^{*} (p\hackscore{2}-\delta p-p)\|\hackscore{0} +
193    \|B A^{-1} B^{*} \delta p\|\hackscore{0} \\
194     & = & \chi^{-} \cdot \epsilon^{-} + \|e\hackscore{2} + B e\hackscore{1}\|\hackscore{0}  \\
195    & \approx & \chi^{-} \cdot \epsilon^{-} + \|e\hackscore{2}\|\hackscore{0} \\
196    & \le & \chi^{-} \cdot \epsilon^{-} + M \tau\hackscore{2} \|B v\hackscore{1}\|\hackscore{0} \\
197  \end{array}  \end{array}
198
199  We use the values $dv=v\hackscore{2}-v$ and $Bv\hackscore{1}=B A^{-1} B^{*} dp$ to measure the error of the  So we choose the value for $\tau\hackscore{2}$ from
200  current approximation $v\hackscore{2}$ and $v\hackscore{2}$ towards the exact solution.   \label{STOKES TOL2}
201  Assuming that the iteration does converge with a convergence rate $\chi^{-}$ we have the estimate   M \tau\hackscore{2} \|B v\hackscore{1}\|\hackscore{0}  \le (\chi^{-})^2 \epsilon^{-}
202
203  \max(\|dv\|\hackscore{1}, \|Bv\hackscore{1}\|\hackscore{0})  in order to make the perturbation for the termination of the pressure iteration a second order effect. We use a
204  \le \chi^{-} \max(\|dv^{-}\|\hackscore{1}, \|Bv\hackscore{1}^{-}\|\hackscore{0})  similar argument for the velcity:
205  + \max(\|\delta v\|\hackscore{1}, \|(B A^{-1} B^{*}) \delta p\|\hackscore{0})  \label{STOKES EST 2}
206    \begin{array}{rcl}
207  were the upper index $-$ referes to the increments in the last step.  \|v\hackscore{2}-v\|\hackscore{1} & \le & \|v\hackscore{2}-\delta v-v\|\hackscore{1} + \| \delta v\|\hackscore{1} \\
208    & \le &  \chi^{-} \cdot \epsilon^{-}  + \| e\hackscore{1} -  A^{-1} B^{*}\delta p \|\hackscore{1} \\
209  Now we will establish estimates for $\|\delta v\|$ and $\|Bv\hackscore{1}\|$ from  & \approx &  \chi^{-} \cdot \epsilon^{-}  + \| e\hackscore{1} \|\hackscore{1} \\
210  formulas~\ref{STOKES ERRORS} where neglect the $\delta p$ in the equation for $\delta v$ assuming that  & \le &  \chi^{-} \cdot \epsilon^{-}  +  K \tau\hackscore{1} \| dv\hackscore{1} - e\hackscore{1} \|\hackscore{1}
211  $\delta p$ is controlled.  \\
212    & \le &  ( 1  + K \tau\hackscore{1}) \chi^{-} \cdot \epsilon^{-}
213  \|\delta v\|\hackscore{1} \approx \|e\hackscore{1}\|\hackscore{1} \le \frac{K \tau\hackscore{1}}{1-K \tau\hackscore{1}}  \| dv\hackscore{1} \|\hackscore{1} \approx \frac{K \tau\hackscore{1}}{1-K \tau\hackscore{1}} \| dv\|\hackscore{1}  \end{array}

and

\|(B A^{-1} B^{*}) \delta p\|\hackscore{0}  \approx \| e\hackscore{2} \|\hackscore{0}
= \| B v\hackscore{2} \|\hackscore{0} \le M \tau\hackscore{2} \| B v\hackscore{1} \|\hackscore{0}

\max(\|dv\|\hackscore{1}, \|Bv\hackscore{1}\|\hackscore{0})
\le \frac{\chi^{-}}{1-max{(M \tau\hackscore{2}, \frac{K \tau\hackscore{1}}{1-K \tau\hackscore{1}})}} \max(\|dv^{-}\|\hackscore{1}, \|Bv\hackscore{1}^{-}\|\hackscore{0}) \label{STOKES ESTIMTED IMPROVEMENT}

where the upper index $-$ refers to the previous iteration step.
If we allow require $max{(M \tau\hackscore{2}, \frac{K \tau\hackscore{1}}{1-K \tau\hackscore{1}})}\le \gamma$
with a given $0<\gamma<1$ we can set
\label{STOKES SET TOL}
\tau\hackscore{2} = \frac{\gamma}{M} \mbox{ and } \tau\hackscore{1} = \frac{1}{K} \frac{\gamma}{1+\gamma}

Problem is that we do not know $M$ and $K$ but can use the convergence rate

\chi := \frac{\max(\|dv\|\hackscore{1}, \|Bv\hackscore{1}\|\hackscore{0})}{\max(\|dv^{-}\|\hackscore{1}, \|Bv\hackscore{1}^{-}\|\hackscore{0})}

to construct estimates of $M$ and $K$. We look at the two cases where our prediction and choice of
the tolerances was good or where we went wrong:

\chi \le \frac{\chi^{-}}{1-\gamma} \mbox{ or }
\chi = \frac{\chi^{-}}{1-\gamma\hackscore{0}} >\frac{\chi^{-}}{1-\gamma}.

which translates to

\gamma\hackscore{0} \le \gamma \mbox{ or } \gamma\hackscore{0}=1-\frac{\chi^{-}}{ \chi}>\gamma>0

In the second case use \ref{STOKES SET TOL} for $\gamma=\gamma\hackscore{0}$ to get new estimates $M^{+}$ and $K^{+}$
for $M$ and $K$:
\label{TOKES CONST UPDATE}
M^{+} =
\frac{\gamma\hackscore{0}}{\gamma} M
\mbox{ and } K^{+} =
\frac{\gamma\hackscore{0}}{\gamma}
\frac{1+\gamma}{1+\gamma\hackscore{0}}
K
\mbox{ with } \gamma\hackscore{0}=\max(1-\frac{\chi^{-}}{ \chi},\gamma)

With these updated constants we can now get better values for the tolerances via an updated value $\gamma^{+}$ for $\gamma$ and the corrected values $M^{+}$ and $K^{+}$ for $M$ and $K$. If we are in the case of convergence which we
meassure by

\chi < \chi\hackscore{max}

where $\chi\hackscore{max}$ is given value with $0<\chi\hackscore{max}<1$. We then consider the following
criteria
\begin{itemize}
\item As we are in case of convergence we try to relax the tolerances by increasing $\gamma$ by a factor $\frac{1}{\omega}$ ($0<\omega<1$). So we would like to choose
$\gamma^{+} \ge \frac{\gamma}{\omega}$.
\item We need to make sure that the estimated convergence rate based on $\gamma^{+}$ will still maintain convergence

\frac{\chi}{1-\gamma^{+}} \le \chi\hackscore{max}

which translates into

\gamma^{+} \le 1 - \frac{\chi}{\chi\hackscore{max}}

Notice that the right hand side is positive.
\item We do not want to iterate more then necessary. So we would like reduce the error in the next just below the
desired tolerance:

\frac{\chi}{1-\gamma^{+}}\max(\|dv\|\hackscore{1}, \|Bv\hackscore{1}\|\hackscore{0}) \ge f \cdot \tau \|v\hackscore{2}\|\hackscore{1}
214
215  where the left hand side provides an estimate for the error to prediced in the next step. $f$ is a  So we choose the value for $\tau\hackscore{1}$ from
216  safty factor. This leads to   \label{STOKES TOL1}
217    K \tau\hackscore{1} \le \chi^{-}
218  \gamma^{+} \le
219  1 -  Assuming we have etsimates for $M$ and $K$ (if none is available, we use the value $1$.)
220  \chi \frac{\max(\|dv\|\hackscore{1}, \|Bv\hackscore{1}\|\hackscore{0})} { f \cdot \tau \|v\hackscore{2}\|\hackscore{1}}  we can use~\ref{STOKES TOL1} and~\ref{STOKES TOL2} to get appropriate values for the tolerances. After
221    the step has been completed we can calculate a new convergence rate $\chi =\frac{\epsilon}{\epsilon^{-}}$.
222    For partical reasons we restrict $\chi$ to be less or equal a given maximum value $\chi\hackscore{max}\le 1$.
223    If we see $\chi \le \chi^{-} (1+\chi^{-})$ our choices for the tolerances was suitable. Otherwise, we need to adjust the values for $K$ and $M$. From the estimates~\ref{STOKES EST 1} and~\ref{STOKES EST 2} we establish
224    \label{STOKES EST 3}
225    \chi \le ( 1 + \max(M \frac{\tau\hackscore{2} \|B v\hackscore{1}\|\hackscore{0}}{\chi^{-} \epsilon^{-}},K \tau\hackscore{1}  ) ) \cdot \chi^{-}
226
227    If we assume that this inequality would be an equation if we would have chosen the right values
228    $M^{+}$ and $K^{+}$ then we get
229    \label{STOKES EST 3b}
230    \chi =  ( 1 + \max(M^{+} \frac{\chi^{-}}{M},K^{+} \frac{\chi^{-}}{K}) ) \cdot \chi^{-}
231
232  \end{itemize}  From this equation we set for
233  Putting all these criteria together we get  than our choice for $K$ was not good anough. In this case we can calaculate a new value
234
235  \gamma^{+}=\min(\max(  K^{+} =  \frac{\chi-\chi^{-}}{(\chi^{-})^2} K
\frac{\gamma}{\omega},
1 -
\chi \frac{\max(\|dv\|\hackscore{1}, \|Bv\hackscore{1}\|\hackscore{0})} { f \cdot \tau \|v\hackscore{2}\|\hackscore{1}}), 1 - \frac{\chi}{\chi\hackscore{max}})
\label{STOKES SET GAMMA PLUS}
236
237  In case we cannot see convergence $\gamma$ is reduced by the factor $\omega$  In practice we will use
238
239  \gamma^{+}=\max(\omega \cdot \gamma, \gamma\hackscore{min})\label{STOKES SET GAMMA PLUS2}  K^{+}  = \max(\frac{\chi-\chi^{-}}{(\chi^{-})^2} K,\frac{1}{2}K,1)
240
241  where $\gamma\hackscore{min}$ is introduced to make sure that $\gamma^{+}$ stays away from zero.  where the second term is used to reduce a potential overestimate of $K$.
242    The same identity is used for to update $M$. The updated $M^{+}$ and $K^{+}$
243    are then use in the next iteration step to control the tolerances.
244  Appling~\ref{STOKES SET TOL} for $\gamma^{+}$ and~\ref{TOKES CONST UPDATE} we get the update formula
245   \label{STOKES CONST UPDATE}  In some cases one can observe that there is a significant change
246  \tau\hackscore{2}^{+} =  in the velocity but the new velocity $v\hackscore{1}$ has still a
247  \frac{\gamma^{+}}{\gamma\hackscore{0}} \tau\hackscore{2}  small divergence, i.e. we have
248  \mbox{ and } \tau\hackscore{1}^{+} =  $\|Bv\hackscore{1}\|\hackscore{0} \ll \|v\hackscore{1}-v\hackscore{0}\|\hackscore{1}$.
249  \frac{\gamma^{+}}{\gamma\hackscore{0}}  In this case we will get a small pressure increment and consequently only very small changes to
250  \frac{1+\gamma\hackscore{0}}{1+\gamma^{+}}  the velocity as a result of the second update step which therefore can be skipped and
251  \tau\hackscore{1}  we can directly repeat the first update step until the increment in velocity becomes
252    significant relative to its divergence. In practice we will ignore the second half of the iteration step
253    as long as
254     \label{STOKES LARGE BV1}
255    \|Bv\hackscore{1}\|\hackscore{0} \le \theta \cdot \|v\hackscore{1}-v\hackscore{0}\|
256
257  to get the new tolerances $\tau\hackscore{1}^{+}$ and $\tau\hackscore{2}^{+}$ to be used in the next iteration step.  where $0<\theta<1$ is a given factor. In this case we will also check the stopping criteria
258  Notice that the coefficients $M$ and $K$ have been eliminated. The calculation of $\gamma\hackscore{0}$ requires  with $v\hackscore{1}\rightarrow v\hackscore{2}$ but we will not correct $M$ in this case.
at least two iteration steps (or a previous time step).
259
260  Typical initial values are  Starting from an initial guess $v\hackscore{0}$ and $p\hackscore{0}$ for velocity and pressure
261   \label{STOKES VALUES}  the solution procedure is implemented as follows.
262  \tau\hackscore{1}=\tau\hackscore{1}=\frac{1}{100} ; \; \omega=\frac{1}{2} \; \gamma=\frac{1}{10} ;\gamma\hackscore{min}=10^{-6}  \begin{enumerate}
263     \item calculate viscosity $\eta(v\hackscore{0},p)\hackscore{0}$ and assemble operator $A$ from $\eta$.
264  where after the first step we set $\gamma\hackscore{0}=\gamma$ as no $\chi^{-}$ is available.   \item calculate the tolerance $\tau\hackscore{1}$ from~\ref{STOKES TOL1}.
265     \item Solve~\ref{STOKES ITER STEP 1} for $dv\hackscore{1}$ with tolerance $\tau\hackscore{1}$.
266     \item update $v\hackscore{1}= v\hackscore{0}+ dv\hackscore{1}$
267     \item if $Bv\hackscore{1}$ is large (see~\ref{STOKES LARGE BV1}):
268     \begin{enumerate}
269     \item calculate the tolerance $\tau\hackscore{2}$ from~\ref{STOKES TOL2}.
270     \item Solve~\ref{STOKES ITER STEP 2} for $dp\hackscore{2}$ and $v\hackscore{2}$ with tolerance $\tau\hackscore{2}$.
271     \item update $p\hackscore{2}\leftarrow p\hackscore{0}+ dp\hackscore{2}$
272     \end{enumerate}
273     \item else:
274      \begin{enumerate}
275      \item update $p\hackscore{2}\leftarrow p$ and $v\hackscore{2}\leftarrow v\hackscore{1}$
276       \end{enumerate}
277       \item calculate convergence messure $\epsilon$ and convergence rate $\chi$
278    \item if stopping criterium~\ref{STOKES STOPPING CRITERIA} holds:
279     \begin{enumerate}
280     \item return $v\hackscore{2}$ and $p\hackscore{2}$
281     \end{enumerate}
282     \item else:
283     \begin{enumerate}
284
285         \item update $M$ and $K$
286         \item goto step 1 with $v\hackscore{0}\leftarrow v\hackscore{2}$ and $p\hackscore{0}\leftarrow p\hackscore{2}$.
287    \end{enumerate}
288    \end{enumerate}
289
290  \subsection{Functions}  \subsection{Functions}
291
# Line 461  sreturns the current absolute tolerance. Line 346  sreturns the current absolute tolerance.
346  \end{methoddesc}  \end{methoddesc}
347
348  \begin{methoddesc}[StokesProblemCartesian]{getSolverOptionsVelocity}{}  \begin{methoddesc}[StokesProblemCartesian]{getSolverOptionsVelocity}{}
349  returns the solver options used  solve the equations~(\ref{V CALC}) for velocity.  returns the solver options used  solve the equations~(\ref{STOKES ITER STEP 1}) and~(\ref{STOKES P OPERATOR}) for velocity.
350  \end{methoddesc}  \end{methoddesc}
351
352  \begin{methoddesc}[StokesProblemCartesian]{getSolverOptionsPressure}{}  \begin{methoddesc}[StokesProblemCartesian]{getSolverOptionsPressure}{}
353  returns the solver options used  solve the equation~(\ref{P PREC}) for pressure.  returns the solver options used solve the preconditioner equation~(\ref{STOKES P PREC}) for pressure.
354  \end{methoddesc}  \end{methoddesc}
355
356  \begin{methoddesc}[StokesProblemCartesian]{getSolverOptionsDiv}{}  \begin{methoddesc}[StokesProblemCartesian]{getSolverOptionsDiv}{}

Legend:
 Removed from v.2793 changed lines Added in v.2843