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

revision 2718 by jfenwick, Mon Sep 7 03:39:45 2009 UTC revision 2719 by gross, Wed Oct 14 06:38:03 2009 UTC
# Line 17  Line 17
17
18  The following sections give a breif overview of the model classes and their corresponding methods.  The following sections give a breif overview of the model classes and their corresponding methods.
19
20  \section{Stokes Problem}  \input{stokessolver}
The velocity \index{velocity} field $v$ and pressure $p$ of an incompressible fluid \index{incompressible fluid} is given as the solution of the Stokes problem\index{Stokes problem}
\label{Stokes 1}
-\left(\eta(v\hackscore{i,j}+ v\hackscore{i,j})\right)\hackscore{,j}+p\hackscore{,i}=f\hackscore{i}-\sigma\hackscore{ij,j}

where $\eta$ is the viscosity, $F\hackscore{i}$ defines an internal force \index{force, internal} and $\sigma\hackscore{ij}$ is an intial stress \index{stress, initial}. We assume an incompressible media:
\label{Stokes 2}
-v\hackscore{i,i}=0

Natural boundary conditions are taken in the form
\label{Stokes Boundary}
\left(\eta(v\hackscore{i,j}+ v\hackscore{i,j})\right)n\hackscore{j}-n\hackscore{i}p=s\hackscore{i}+\sigma\hackscore{ij} n\hackscore{j}

which can be overwritten by constraints of the form
\label{Stokes Boundary0}
v\hackscore{i}(x)=v^D\hackscore{i}(x)

at some locations $x$ at the boundary of the domain. The index $i$ may depend on the location $x$ on the boundary.
$v^D$ is a given function on the domain.
21
\subsection{Solution Method \label{STOKES SOLVE}}
In block form equation equations~\ref{Stokes 1} and~\ref{Stokes 2} takes the form of a saddle point problem

\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]

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}.
We use iterative techniques to solve this problem. To make sure that the incomressibilty condition holds
with sufficient accuracy we check for

\|v\hackscore{k,k}\| \hackscore \le  \epsilon
\|\sqrt{v\hackscore{j,k}v\hackscore{j,k}}\|
\label{STOKES STOP}

where $\epsilon$ is the desired relative accuracy and

\|p\|^2= \int\hackscore{\Omega} p^2 \; dx
\label{PRESSURE NORM}

defines the $L^2$-norm. We use the Uzawa scheme \index{Uzawa scheme} to solve the problem.

In fact the first equation in~\ref{SADDLEPOINT} gives for a known pressure

v=A^{-1}(G-B^{*}p)
\label{V CALC}

which is inserted into the second equation leading to

S p =  B A^{-1} G

with the Schur complement \index{Schur complement} $S=BA^{-1}B^{*}$. This problem can be solved iteratively
with the preconditioner $\hat{S}$ defined as $q=\hat{S}^{-1}p$ by solving
\label{P PREC}
\frac{1}{\eta}q = p

see~\cite{ELMAN} for more details. Note that the residual for the current approximation $p$ is given as

r=B A^{-1} (G - B^* p) = Bv

where $v$ is given by~\ref{V CALC}.

If one uses the generalized minimal residual method (GMRES) \index{generalized minimal residual method!GMRES}
the method is directly applied to the preconditioned system

\hat{S}^{-1} S p =  \hat{S}^{-1} B A^{-1} G

We use the norm

\|p\|\hackscore{GMRES} = \|\hat{S} p \|

Notice that for the residual $\hat{r}=\hat{S}^{-1} r$ one has

\

If $p^{0}$ provides an initial guess for the pressure we use~\ref{V CALC} to get a first initial guess for the
velocity $v^{0}$ which we use to set an absolute tolerance $ATOL =\epsilon \|\sqrt{v^{0}\hackscore{j,k}v^{0}\hackscore{j,k}}\|$.
The GMRES is terminated when

\|\hat{r}\|\hackscore{GMRES} \le ATOL

Notice that $\|\hat{r}\|\hackscore{GMRES}= \|r \| = \|Bv\| = \|v\hackscore{k,k}\|$ so we we can expect that
the target stopping criterion~\ref{STOKES STOP} is fullfilled. However, if $v$ is very different from the
initial choice of $v^{0}$ the value of $ATOL$ is corrected and GMRES is restarted with a new tolerance. For time dependend problems this apprach works well as value for $p$ form a previous time step provides a good initial guess.

Alternatively, as $S$ is symmetric and positive definite one can apply the preconditioned conjugate gradient method (PCG) \index{preconditioned conjugate gradient method!PCG}. PCG use the norm

\|r\|\hackscore{PCG}^2 = \int\hackscore{\Omega} r \hat{S}^{-1}r \; dx = \int\hackscore{\Omega}  \eta r^2 \; dx

To take the extra factor $\eta$ into consideration when checking the stopping criterion we use the following
definition for $ATOL$:

ATOL = \epsilon \frac{\|\sqrt{v^{0}\hackscore{j,k}v^{0}\hackscore{j,k}}\|  }{\|v^{0}\hackscore{k,k}\|}
\|v^{0}\hackscore{k,k}\|\hackscore{PCG}

\subsection{Functions}

opens the Stokes problem\index{Stokes problem} on the \Domain domain. The approximation
order needs to be two.
the tolerances for all subproblems are set automatically.
\end{classdesc}

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 compontents where the velocity is fixed are set by
the values of \var{fixed_u_mask}. 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. The values of \var{v} marked
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 talerance is set to 0.
\end{methoddesc}

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

\begin{methoddesc}[StokesProblemCartesian]{getSolverOptionsVelocity}{}
returns the solver options used  solve the equations~(\ref{V CALC}) for velocity.
\end{methoddesc}

\begin{methoddesc}[StokesProblemCartesian]{getSolverOptionsPressure}{}
returns the solver options used  solve the equation~(\ref{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)
(whereZero(x[1])*[0.,1.]+whereZero(x[1]-1))*[1.,1]
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}
22
23  \section{Darcy Flux}  \section{Darcy Flux}
24  We want to calculate the velocity $u$ and pressure $p$ on a domain $\Omega$ solving  We want to calculate the velocity $u$ and pressure $p$ on a domain $\Omega$ solving

Legend:
 Removed from v.2718 changed lines Added in v.2719