doc/user/stokessolver.tex

\section{The Stokes Problem}
\label{STOKES PROBLEM} 
In this section we discuss how to solve the Stokes problem which is defined as follows:

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}
\begin{equation}\label{Stokes 1}
-\left(\eta(v\hackscore{i,j}+ v\hackscore{j,i})\right)\hackscore{,j}+p\hackscore{,i}=f\hackscore{i}-\sigma\hackscore{ij,j}
\end{equation}
where  $f\hackscore{i}$ defines an internal force \index{force, internal} and $\sigma\hackscore{ij}$ is an initial 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.

We assume an incompressible media:
\begin{equation}\label{Stokes 2}
-v\hackscore{i,i}=0
\end{equation}
Natural boundary conditions are taken in the form 
\begin{equation}\label{Stokes Boundary}
\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}
\end{equation}
which can be overwritten by constraints of the form 
\begin{equation}\label{Stokes Boundary0}
v\hackscore{i}(x)=v^D\hackscore{i}(x)
\end{equation}
at some locations $x$ at the boundary of the domain. $s\hackscore{i}$ defines a normal stress and 
$\alpha\ge 0$ the spring constant for restoring normal force.
The index $i$ may depend on the location $x$ on the boundary.
$v^D$ is a given function on the domain.

\subsection{Solution Method \label{STOKES SOLVE}}
If we assume that $\eta$ is independent from the velocity and pressure equations~\ref{Stokes 1} and~\ref{Stokes 2} 
can be written in the block form
\begin{equation}
\left[ \begin{array}{cc}
A     & B^{*} \\
B & 0 \\
\end{array} \right]
\left[ \begin{array}{c}
v \\
p \\
\end{array} \right]
=\left[ \begin{array}{c}
G \\
0 \\
\end{array} \right]
\label{STOKES}
\end{equation}
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).
For more details on the mathematics see references \cite{AAMIRBERKYAN2008,MBENZI2005}. 

If $v\hackscore{0}$ and $p\hackscore{0}$ are given initial guesses for
velocity and pressure we calculate a correction $dv$ for the velocity by solving the first
equation of equation~\ref{STOKES}
 \begin{equation}\label{STOKES ITER STEP 1}
 A dv\hackscore{1} = G - A v\hackscore{0} - B^{*} p\hackscore{0}
\end{equation}
We then insert the new approximation $v\hackscore{1}=v\hackscore{0}+dv\hackscore{1}$ to calculate a correction $dp\hackscore{2}$
for the pressure and an additional correction $dv\hackscore{2}$ for the velocity by solving
 \begin{equation}
 \begin{array}{rcl}
 B A^{-1} B^{*} dp\hackscore{2} & = & Bv\hackscore{1} \\
 A dv\hackscore{2} & = & B^{*} dp\hackscore{2} 
\end{array}
 \label{STOKES ITER STEP 2}
 \end{equation}
The new velocity and pressure are then given by $v\hackscore{2}=v\hackscore{1}-dv\hackscore{2}$ and
$p\hackscore{2}=p\hackscore{0}+dp\hackscore{2}$ which will fulfill the block system~\ref{STOKES}. 
This solution strategy is called the Uzawa scheme \index{Uzawa scheme}. 

There is a problem with this scheme: In practice we will use an iterative scheme
to solve any problem for operator $A$. So we will be unable to calculate the operator
$ B A^{-1} B^{*}$ required for $dp\hackscore{2}$ explicitly. In fact, we need to use another
iterative scheme to solve the first equation in~\ref{STOKES ITER STEP 2} where in each iteration step
an iterative solver for $A$ is applied. Another issue is the fact that the
viscosity $\eta$ may depend on velocity or pressure and so we need to iterate over the 
three equations~\ref{STOKES ITER STEP 1} and~\ref{STOKES ITER STEP 2}. 

In the following we will use the two norms
\begin{equation}
\|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}
\end{equation}
for velocity $v$ and pressure $p$. The iteration is terminated if the stopping criteria
 \begin{equation} \label{STOKES STOPPING CRITERIA}
\max(\|Bv\hackscore{1}\|\hackscore{0},\|v\hackscore{2}-v\hackscore{0}\|\hackscore{1}) \le \tau \cdot \|v\hackscore{2}\|\hackscore{1} 
 \end{equation}
 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
$\|Bv\hackscore{1}\|\hackscore{0}$ equals the
norm of $B A^{-1} B^{*} dp\hackscore{2}$ and consequently provides a norm for the pressure correction.

We want to optimize the tolerance choice for solving~\ref{STOKES ITER STEP 1}
and~\ref{STOKES ITER STEP 2}. To do this we write the iteration scheme as a fixed point problem. Here
we consider the errors produced by the iterative solvers being used. 
From Equation~\ref{STOKES ITER STEP 1} we have 
\begin{equation} \label{STOKES total V1}
v\hackscore{1} = e\hackscore{1} + v\hackscore{0} + A^{-1} ( G - Av\hackscore{0} - B^{*} p\hackscore{0} ) 
\end{equation}
where $ e\hackscore{1}$ is the error when solving~\ref{STOKES ITER STEP 1}.  
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 
\begin{equation} 
\| 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} 
\end{equation}
where $\tau\hackscore{1}$ is the tolerance which we need to choose. This translates into the condition
\begin{equation} 
\| e\hackscore{1} \|\hackscore{1} \le K \tau\hackscore{1} \| dv\hackscore{1} - e\hackscore{1} \|\hackscore{1} 
\end{equation}
The constant $K$ represents some uncertainty combining a variety of unknown factors such as the 
norm being used and the condition number of the stiffness matrix.  
From the first equation of~\ref{STOKES ITER STEP 2} we have
\begin{equation}\label{STOKES total P2}
p\hackscore{2} =  p\hackscore{0} + (B A^{-1} B^{*})^{-1} (e\hackscore{2} + Bv\hackscore{1} )
\end{equation}
where $e\hackscore{2}$ represents the error when solving~\ref{STOKES ITER STEP 2}.
We use an iterative preconditioned conjugate gradient method (PCG) \index{linear solver!PCG}\index{PCG} with iteration 
operator $B A^{-1} B^{*}$ using the $\|.\|\hackscore{0}$ norm. As suitable preconditioner \index{preconditioner} for the iteration
operator we use $\frac{1}{\eta}$ \cite{ELMAN}, ie 
the evaluation of the preconditioner $P$ for a given pressure increment $q$ is the solution of
\begin{equation} \label{STOKES P PREC}
\frac{1}{\eta} (Pq) = q \; . 
\end{equation}
Note that in each evaluation of the iteration operator $q=B A^{-1} B^{*} s$ one needs to solve
the problem
\begin{equation} \label{STOKES P OPERATOR}
A w = B^{*} s 
\end{equation}
with sufficient accuracy to return $q=Bw$. We assume that the desired tolerance is 
sufficiently small, for instance one can take $\tau\hackscore{2}^2$ 
where $\tau\hackscore{2}$ is the tolerance for~\ref{STOKES ITER STEP 2}.

In an implementation we use the fact that the residual $r$ is given as
\begin{equation} \label{STOKES RES }
 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}
\end{equation}
In particular we have $e\hackscore{2} = B v\hackscore{2}$
So the residual $r$ is represented by the updated velocity $v\hackscore{2}=v\hackscore{1}-dv\hackscore{2}$. In practice, if
one uses the velocity to represent the residual $r$ there is no need 
to recover $dv\hackscore{2}$ in~\ref{STOKES ITER STEP 2} after $dp\hackscore{2}$ has been calculated.
In PCG the iteration is terminated if
\begin{equation} \label{STOKES P OPERATOR ERROR}
\| P^{\frac{1}{2}}B v\hackscore{2} \|\hackscore{0} \le \tau\hackscore{2} \| P^{\frac{1}{2}}B v\hackscore{1} \|\hackscore{0}
\end{equation}
where $\tau\hackscore{2}$ is the given tolerance. This translates into
\begin{equation} \label{STOKES P OPERATOR ERROR 2}
\|e\hackscore{2}\|\hackscore{0} = \| B v\hackscore{2} \|\hackscore{0} \le M \tau\hackscore{2} \| B v\hackscore{1} \|\hackscore{0}
\end{equation}
where $M$ is taking care of the fact that $P^{\frac{1}{2}}$ is dropped.   

As we assume that there is no significant error from solving with the operator $A$ we have 
\begin{equation} \label{STOKES total V2}
v\hackscore{2} =  v\hackscore{1} - dv\hackscore{2} 
= v\hackscore{1}  - A^{-1} B^{*}dp 
\end{equation}
Combining the equations~\ref{STOKES total V1},~\ref{STOKES total P2} and~\ref{STOKES total V2} and
setting the errors to zero we can write the solution process as a fix point problem 
\begin{equation} 
v = \Phi(v,p) \mbox{ and } p = \Psi(u,p) 
\end{equation}
with suitable functions $\Phi(v,p)$ and $ \Psi(v,p)$ representing the iteration operator without 
errors. In fact for a linear problem,  $\Phi$ and $\Psi$ are constant. With this notation we can write
the update step in the form $p\hackscore{2}= \delta p + \Psi(v\hackscore{0},p\hackscore{0})$ and 
$v\hackscore{2}= \delta v + \Phi(v\hackscore{0},p\hackscore{0})$ where 
the total error $\delta p$ and $\delta v$ are given as
\begin{equation} 
 \begin{array}{rcl}
\delta p & = &  (B A^{-1} B^{*})^{-1} ( e\hackscore{2} + B e\hackscore{1} ) \\
\delta v & = &  e\hackscore{1} -  A^{-1} B^{*}\delta p  \;.
\end{array}\label{STOKES ERRORS}
\end{equation}
Notice that $B\delta v = - e\hackscore{2}=-Bv\hackscore{2}$. Our task is now to choose the tolerances
$\tau\hackscore{1}$ and $\tau\hackscore{2}$ such that the global errors $\delta p$ and $\delta v$
do not stop the convergence of the iteration process. 

To measure convergence we use
\begin{equation} 
\epsilon = \max(\|v\hackscore{2}-v\|\hackscore{1}, \|B A^{-1} B^{*} (p\hackscore{2}-p)\|\hackscore{0})
\end{equation}
In practice using the fact that $B A^{-1} B^{*} (p\hackscore{2}-p\hackscore{0}) = B v\hackscore{1}$
and assuming that $v\hackscore{2}$ gives a better approximation to the true $v$ than
$v\hackscore{0}$ we will use the estimate
\begin{equation} 
\epsilon = \max(\|v\hackscore{2}-v\hackscore{0}\|\hackscore{1}, \|B v\hackscore{1}\|\hackscore{0})
\end{equation}
to estimate the progress of the iteration step after the step is completed. 
Note that the estimate of $\epsilon$   
used in the stopping criteria~\ref{STOKES STOPPING CRITERIA}. If $\chi^{-}$ is the convergence rate assuming
exact calculations, i.e. $e\hackscore{1}=0$ and $e\hackscore{2}=0$, we are expecting 
to maintain $\epsilon \le \chi^{-} \cdot \epsilon^{-}$. For the 
pressure increment we get:
\begin{equation} \label{STOKES EST 1}
\begin{array}{rcl}
\|B A^{-1} B^{*} (p\hackscore{2}-p)\|\hackscore{0}
 & \le & \|B A^{-1} B^{*} (p\hackscore{2}-\delta p-p)\|\hackscore{0} +
\|B A^{-1} B^{*} \delta p\|\hackscore{0} \\
 & = & \chi^{-} \cdot \epsilon^{-} + \|e\hackscore{2} + B e\hackscore{1}\|\hackscore{0}  \\
& \approx & \chi^{-} \cdot \epsilon^{-} + \|e\hackscore{2}\|\hackscore{0} \\
& \le & \chi^{-} \cdot \epsilon^{-} + M \tau\hackscore{2} \|B v\hackscore{1}\|\hackscore{0} \\  
\end{array}
\end{equation}
So we choose the value for $\tau\hackscore{2}$ from 
\begin{equation} \label{STOKES TOL2}
 M \tau\hackscore{2} \|B v\hackscore{1}\|\hackscore{0}  \le (\chi^{-})^2 \epsilon^{-}
\end{equation}
in order to make the perturbation for the termination of the pressure iteration a second order effect. We use a
similar argument for the velocity:
\begin{equation}\label{STOKES EST 2}
\begin{array}{rcl}
\|v\hackscore{2}-v\|\hackscore{1} & \le & \|v\hackscore{2}-\delta v-v\|\hackscore{1} + \| \delta v\|\hackscore{1} \\
& \le &  \chi^{-} \cdot \epsilon^{-}  + \| e\hackscore{1} -  A^{-1} B^{*}\delta p \|\hackscore{1} \\
& \approx &  \chi^{-} \cdot \epsilon^{-}  + \| e\hackscore{1} \|\hackscore{1} \\
& \le &  \chi^{-} \cdot \epsilon^{-}  +  K \tau\hackscore{1} \| dv\hackscore{1} - e\hackscore{1} \|\hackscore{1}
\\
& \le &  ( 1  + K \tau\hackscore{1}) \chi^{-} \cdot \epsilon^{-}
\end{array}
\end{equation}
So we choose the value for $\tau\hackscore{1}$ from
\begin{equation} \label{STOKES TOL1}
K \tau\hackscore{1} \le \chi^{-}
\end{equation}
Assuming we have estimates for $M$ and $K$ (if none is available, we use the value $1$.)
we can use~\ref{STOKES TOL1} and~\ref{STOKES TOL2} to get appropriate values for the tolerances. After
the step has been completed we can calculate a new convergence rate $\chi =\frac{\epsilon}{\epsilon^{-}}$. 
For partial reasons we restrict $\chi$ to be less or equal a given maximum value $\chi\hackscore{max}\le 1$.
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
\begin{equation}\label{STOKES EST 3}
\chi \le ( 1 + \max(M \frac{\tau\hackscore{2} \|B v\hackscore{1}\|\hackscore{0}}{\chi^{-} \epsilon^{-}},K \tau\hackscore{1}  ) ) \cdot \chi^{-} 
\end{equation}
If we assume that this inequality would be an equation if we would have chosen the right values
$M^{+}$ and $K^{+}$ then we get 
\begin{equation}\label{STOKES EST 3b}
\chi =  ( 1 + \max(M^{+} \frac{\chi^{-}}{M},K^{+} \frac{\chi^{-}}{K}) ) \cdot \chi^{-} 
\end{equation}
From this equation we set for 
than our choice for $K$ was not good enough. In this case we can calculate a new value
 \begin{equation}
K^{+} =  \frac{\chi-\chi^{-}}{(\chi^{-})^2} K
\end{equation}
In practice we will use 
 \begin{equation}
K^{+}  = \max(\frac{\chi-\chi^{-}}{(\chi^{-})^2} K,\frac{1}{2}K,1)
\end{equation}
where the second term is used to reduce a potential overestimate of $K$.  
The same identity is used for to update $M$. The updated $M^{+}$ and $K^{+}$ 
are then use in the next iteration step to control the tolerances. 

In some cases one can observe that there is a significant change 
in the velocity but the new velocity $v\hackscore{1}$ has still a 
small divergence, i.e. we have
$\|Bv\hackscore{1}\|\hackscore{0} \ll \|v\hackscore{1}-v\hackscore{0}\|\hackscore{1}$. 
In this case we will get a small pressure increment and consequently only very small changes to
the velocity as a result of the second update step which therefore can be skipped and
we can directly repeat the first update step until the increment in velocity becomes
significant relative to its divergence. In practice we will ignore the second half of the iteration step
as long as 
 \begin{equation}\label{STOKES LARGE BV1}
\|Bv\hackscore{1}\|\hackscore{0} \le \theta \cdot \|v\hackscore{1}-v\hackscore{0}\| 
\end{equation}
where $0<\theta<1$ is a given factor. In this case we will also check the stopping criteria 
with $v\hackscore{1}\rightarrow v\hackscore{2}$ but we will not correct $M$ in this case.

Starting from an initial guess $v\hackscore{0}$ and $p\hackscore{0}$ for velocity and pressure
the solution procedure is implemented as follows.
\begin{enumerate}
 \item calculate viscosity $\eta(v\hackscore{0},p)\hackscore{0}$ and assemble operator $A$ from $\eta$.
 \item calculate the tolerance $\tau\hackscore{1}$ from~\ref{STOKES TOL1}.
 \item Solve~\ref{STOKES ITER STEP 1} for $dv\hackscore{1}$ with tolerance $\tau\hackscore{1}$.
 \item update $v\hackscore{1}= v\hackscore{0}+ dv\hackscore{1}$
 \item if $Bv\hackscore{1}$ is large (see~\ref{STOKES LARGE BV1}):
 \begin{enumerate}
 \item calculate the tolerance $\tau\hackscore{2}$ from~\ref{STOKES TOL2}.
 \item Solve~\ref{STOKES ITER STEP 2} for $dp\hackscore{2}$ and $v\hackscore{2}$ with tolerance $\tau\hackscore{2}$.
 \item update $p\hackscore{2}\leftarrow p\hackscore{0}+ dp\hackscore{2}$
 \end{enumerate}
 \item else:
  \begin{enumerate}
  \item update $p\hackscore{2}\leftarrow p$ and $v\hackscore{2}\leftarrow v\hackscore{1}$
   \end{enumerate}
   \item calculate convergence measure $\epsilon$ and convergence rate $\chi$
\item if stopping criteria~\ref{STOKES STOPPING CRITERIA} holds:
 \begin{enumerate}
 \item return $v\hackscore{2}$ and $p\hackscore{2}$
 \end{enumerate}
 \item else:
 \begin{enumerate}

     \item update $M$ and $K$
     \item goto step 1 with $v\hackscore{0}\leftarrow v\hackscore{2}$ and $p\hackscore{0}\leftarrow p\hackscore{2}$.
\end{enumerate}
\end{enumerate}

\subsection{Functions}

\begin{classdesc}{StokesProblemCartesian}{domain}
opens the Stokes problem\index{Stokes problem} on the \Domain domain. The domain
needs to support LBB compliant elements for the Stokes problem, see~\cite{LBB} for details~\index{LBB condition}.
For instance one can use second order polynomials for velocity and 
first order polynomials for the pressure on the same element. Alternatively, one can use 
macro elements\index{macro elements} using linear polynomial for both pressure and velocity    with a subdivided
element for the velocity. Typically, the macro element is more cost effective. The fact that pressure and velocity are represented in different way is expressed by
\begin{python}
velocity=Vector(0.0, Solution(mesh))
pressure=Scalar(0.0, ReducedSolution(mesh))
\end{python}
\end{classdesc}

\begin{methoddesc}[StokesProblemCartesian]{initialize}{\optional{f=Data(), \optional{fixed_u_mask=Data(), \optional{eta=1, \optional{surface_stress=Data(), \optional{stress=Data()}, \optional{
restoration_factor=0}}}}}}
assigns values to the model parameters. In any call all values must be set.
\var{f} defines the external force $f$, \var{eta} the viscosity $\eta$,
\var{surface_stress} the surface stress $s$ and \var{stress} the initial stress $\sigma$.
The locations and components where the velocity is fixed are set by 
the values of \var{fixed_u_mask}. \var{restoration_factor} defines the restoring force factor $\alpha$. 
The method will try to cast the given values to appropriate 
\Data class objects.
\end{methoddesc}

\begin{methoddesc}[StokesProblemCartesian]{solve}{v,p
\optional{, max_iter=100 \optional{, verbose=False \optional{, usePCG=True }}}}
solves the problem and return approximations for velocity and pressure. 
The arguments \var{v} and \var{p} define initial guess.
\var{v} must have function space \var{Solution(domain)} and
\var{p} must have function space \var{ReducedSolution(domain)}.
The values of \var{v} marked
by \var{fixed_u_mask} remain unchanged. 
If \var{usePCG} is set to \True 
reconditioned conjugate gradient method (PCG) \index{preconditioned conjugate gradient method!PCG}  scheme is used. Otherwise the problem is solved generalized minimal residual method (GMRES) \index{generalized minimal residual method!GMRES}. In most cases 
the PCG scheme is more efficient.
\var{max_iter} defines the maximum number of iteration steps. 

If \var{verbose} is set to \True informations on the progress of of the solver are printed.
\end{methoddesc}


\begin{methoddesc}[StokesProblemCartesian]{setTolerance}{\optional{tolerance=1.e-4}}
sets the tolerance in an appropriate norm relative to the right hand side. The tolerance must be non-negative and less than 1.
\end{methoddesc}
\begin{methoddesc}[StokesProblemCartesian]{getTolerance}{}
returns the current relative tolerance.
\end{methoddesc}
\begin{methoddesc}[StokesProblemCartesian]{setAbsoluteTolerance}{\optional{tolerance=0.}}
sets the absolute tolerance for the error in the relevant norm. The tolerance must be non-negative. Typically the
absolute tolerance is set to 0.
\end{methoddesc}

\begin{methoddesc}[StokesProblemCartesian]{getAbsoluteTolerance}{}
returns the current absolute tolerance.
\end{methoddesc}

\begin{methoddesc}[StokesProblemCartesian]{getSolverOptionsVelocity}{}
returns the solver options used  solve the equations~(\ref{STOKES ITER STEP 1}) and~(\ref{STOKES P OPERATOR}) for velocity.
\end{methoddesc}

\begin{methoddesc}[StokesProblemCartesian]{getSolverOptionsPressure}{}
returns the solver options used solve the preconditioner equation~(\ref{STOKES P PREC}) for pressure.
\end{methoddesc}

\begin{methoddesc}[StokesProblemCartesian]{getSolverOptionsDiv}{}
set the solver options for solving the equation to project the divergence of the velocity onto the function space of pressure.
\end{methoddesc}


\subsection{Example: Lit Driven Cavity}
 The following script \file{lit\hackscore driven\hackscore cavity.py} 
\index{scripts!\file{helmholtz.py}} which is available in the \ExampleDirectory
illustrates the usage of the \class{StokesProblemCartesian} class to solve
the lit driven cavity problem:
\begin{python}
from esys.escript import *
from esys.finley import Rectangle
from esys.escript.models import StokesProblemCartesian
NE=25
dom = Rectangle(NE,NE,order=2)
x = dom.getX()
sc=StokesProblemCartesian(dom)
mask= (whereZero(x[0])*[1.,0]+whereZero(x[0]-1))*[1.,0] + \
      (whereZero(x[1])*[0.,1.]+whereZero(x[1]-1))*[1.,1]
sc.initialize(eta=.1, fixed_u_mask= mask)
v=Vector(0.,Solution(dom))
v[0]+=whereZero(x[1]-1.)
p=Scalar(0.,ReducedSolution(dom))
v,p=sc.solve(v,p, verbose=True)
saveVTK("u.xml",velocity=v,pressure=p)
\end{python}