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

revision 2748 by gross, Tue Nov 17 07:32:59 2009 UTC revision 2793 by gross, Tue Dec 1 06:10:10 2009 UTC
# Line 26  The index $i$ may depend on the location Line 26  The index $i$ may depend on the location
26  $v^D$ is a given function on the domain.  $v^D$ is a given function on the domain.
27
28  \subsection{Solution Method \label{STOKES SOLVE}}  \subsection{Solution Method \label{STOKES SOLVE}}
29  In block form equation equations~\ref{Stokes 1} and~\ref{Stokes 2} takes the form of a saddle point problem  If we assume that $\eta$ is independent from the velocity and pressure equations~\ref{Stokes 1} and~\ref{Stokes 2}
30  \index{saddle point problem}  can be written in the block form
31
32  \left[ \begin{array}{cc}  \left[ \begin{array}{cc}
33  A     & B^{*} \\  A     & B^{*} \\
# Line 43  G \\ Line 43  G \\
43  \end{array} \right]  \end{array} \right]
45
46  where $A$ is coercive (assuming $A$ is not depending on $v$ or $p$), 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
50    velocity and pressure we calculate a correction $dv$ for the velocity by solving the first
53     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}$
56    for the pressure and an additional correction $dv\hackscore{2}$ for the velocity by solving
57
58     \begin{array}{rcl}
59     B A^{-1} B^{*} dp\hackscore{2} & = & Bv\hackscore{1} \\
60     A dv\hackscore{2} & = & B^{*} dp\hackscore{2}
61    \end{array}
63
64    The new velocity and pressure are then given by $v=v\hackscore{1}-dv\hackscore{2}$ and
65    $p=p\hackscore{0}+dp\hackscore{2}$ which will fullfill the block system~\ref{SADDLEPOINT}.
66    This solution strategy is called the Uzawa scheme \index{Uzawa scheme}. There are
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
69    $B A^{-1} B^{*}$ required for $dp\hackscore{2}$. In fact, we need to use an iterative scheme
70    to solve the first equation in~\ref{SADDLEPOINT ITER STEP 2} where in each iteration step
71    an iterative solver for $A$ is applied. The second issue is that
72    viscosity $\eta$ may depend of velocity or pressure and so we need to iterate over the
73    three equations~\ref{SADDLEPOINT ITER STEP 1} and~\ref{SADDLEPOINT ITER STEP 2}. Using the
74    two norms
75
76    \|v\|\hackscore{1}^2 = \int\hackscore{\Omega} v\hackscore{j,k}v\hackscore{j,k} \; dx
77    \mbox{ and }
78    \|p\|\hackscore{0}^2= \int\hackscore{\Omega} p^2 \; dx.
79    \label{STOKES STOP}
80
81    to define the stopping criterium
82
83    \max(\|Bv\hackscore{1}\|\hackscore{0},\|v-v\hackscore{0}\|\hackscore{1}) \le \tau \cdot \|v\|\hackscore{1}
84
85    to terminate the overall iteration with a overall given tolerance $0<\tau<1$.
86    Notice that because of the first equation of
87    and~\ref{SADDLEPOINT ITER STEP 2} $\|Bv\hackscore{1}\|\hackscore{0}$ is equal to the
88    norm of $B A^{-1} B^{*} dp\hackscore{2}$ and consequently provides a norm for the  pressure correction.
89
90    Let us first have a look at the solution process~\ref{SADDLEPOINT ITER STEP 1} which needs to resolve the
91    non-linearity of the viscosity coeficient. We will run a iterative process by resolving
92    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
93    the iteration process over the second set of equation if $\|Bv\hackscore{1}\|\hackscore{0}$
94    is smaller then $\|v\hackscore{1}-v\hackscore{0}\|\hackscore{1}$. So we will terminate the non-linear iteration
95    process of the first equation if
96
97    \theta \cdot \|v\hackscore{1}-v\hackscore{0}\| \le \|Bv\hackscore{1}\|\hackscore{0}
98    \mbox{ or }
99    \|v\hackscore{1}-v\hackscore{0}\|\hackscore{1} \le \tau \cdot \|v\hackscore{1}\|\hackscore{1}
100
101    where $0<\theta<1$ is a given factor (typically $\theta=0.5$).
102
103    We need to think about appropriate stopping criteria when solving
104    \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
105    convergence of the iteration process one level above.
106    We will use a sparse matrix solver so we have not full control on the norm $\|.\|\hackscore{s}$ used in the stopping criteria
107
108    \| 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}
109
110    If $e\hackscore{1}$ is the error of the returned solution $dv\hackscore{1}$ this translates into the condition
111
112    \| e\hackscore{1} \|\hackscore{1} \le K \tau\hackscore{1} \| dv\hackscore{1} - e\hackscore{1} \|\hackscore{1}
113
114    The constant $K$ represents some uncertainty combining a variety of unknown factors such as the
115    solver being used and the condition number of the stiffness matrix.  This leads to the estimate
116
117    \| dv\hackscore{1} \|\hackscore{1} \le (1+K \tau\hackscore{1}) \chi \| dv\hackscore{1}^- \|\hackscore{1}
118    \label{DV COND}
119
120    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
121    to have $K \tau\hackscore{1} \le \chi$ which leads to the condition
122
123     \tau\hackscore{1} \le \frac{\min(\chi,\frac{1}{2})}{K}
124     \label{NEW TAU 1}
125
126    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
127    access to the estimate of the convergence rate $\chi^-$ of the previous iteration step we
128    get
129
130     \frac{1}{K} = \frac{\chi^-}{\chi-\chi^{-}} \tau\hackscore{1}^-
131
132    if $\chi \ge \chi^- (1 + \chi^-)$ which means that our prediction for a suitable tolerane based on the avaibale $K$ value was wrong.
133    So we start with the value $K=1$ in~\ref{NEW TAU 1} and recalculate
134    $K$ if we have underestimate $K$.
135
136
137
138
139
140
141
142    If $dv\hackscore{1}$
143    From Equation~\ref{SADDLEPOINT ITER STEP 1} we have
144
145    v\hackscore{1} = e\hackscore{1} + A^{-1} ( G - B^{*} p )
146
147    This translates into the conditoon
148
149    \| e\hackscore{1} \|\hackscore{1} \le K \tau\hackscore{1} \| dv\hackscore{1} - e\hackscore{1} \|\hackscore{1}
150    \mbox{ therefore }
151    \| e\hackscore{1} \|\hackscore{1} \le \frac{K \tau\hackscore{1}}{1-K \tau\hackscore{1}}  \| dv\hackscore{1} \|\hackscore{1}
152
153    The constant $K$ represents some uncertainty combining a variety of unknown factors such as the
154    solver being used and the condition number of the stiffness matrix.
155
156
157
158
159
160    ===================================
161
162
163  We use iterative techniques to solve this problem: Given an approximation $v$ and $p$ for  We use iterative techniques to solve this problem: Given an approximation $v$ and $p$ for
164  velocity and pressure we perform the following steps in the Uzawa scheme \index{Uzawa scheme} style:  velocity and pressure we perform the following steps in the Uzawa scheme \index{Uzawa scheme} style:
165  \begin{enumerate}  \begin{enumerate}
# Line 66  velocity and pressure we perform the fol Line 181  velocity and pressure we perform the fol
181   \item update $p\hackscore{2}\leftarrow p+ dp$ and $v\hackscore{2}= v\hackscore{1} - dv\hackscore{2}$   \item update $p\hackscore{2}\leftarrow p+ dp$ and $v\hackscore{2}= v\hackscore{1} - dv\hackscore{2}$
182   \item goto Step 0 with $v=v\hackscore{2}$ and $p=p\hackscore{2}$.   \item goto Step 0 with $v=v\hackscore{2}$ and $p=p\hackscore{2}$.
183  \end{enumerate}  \end{enumerate}
184  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  where $\tau$ is the given tolerance. Notice that the operator $A$ if it is depending on $v$ or $p$ is updated once only.
185    In this algorithm
186
187  \|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
188  \mbox{ and }  \mbox{ and }
# Line 150  We also can give a formula for the error Line 266  We also can give a formula for the error
266
267   \begin{array}{rcl}   \begin{array}{rcl}
268  \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} ) \\
269  \delta v & = &   e\hackscore{1} -  A^{-1} B^{*}\delta p  \delta v & = &   e\hackscore{1} -  A^{-1} B^{*}\delta p  \;.
270  \end{array}\label{STOKES ERRORS}  \end{array}\label{STOKES ERRORS}
271
272  Notice that $B\delta v = - e\hackscore{2}=-Bv\hackscore{2}$  Notice that $B\delta v = - e\hackscore{2}=-Bv\hackscore{2}$.
273
274  With this notation  With this notation
275
276   \begin{array}{rcl}   \begin{array}{rcl}
# Line 161  v\hackscore{2} = \Phi(v,p) + \delta v \\ Line 278  v\hackscore{2} = \Phi(v,p) + \delta v \\
278  p\hackscore{2} = \Psi(v,p) + \delta p \\  p\hackscore{2} = \Psi(v,p) + \delta p \\
279  \end{array}  \end{array}
280
281  We use the $dv=v\hackscore{2}-v$ and $Bv\hackscore{1}=B A^{-1} B^{*} dp$ to measure the error of the  We use the values $dv=v\hackscore{2}-v$ and $Bv\hackscore{1}=B A^{-1} B^{*} dp$ to measure the error of the
282  current approximation $v\hackscore{2}$ and $v\hackscore{2}$ towards the exact solution.  current approximation $v\hackscore{2}$ and $v\hackscore{2}$ towards the exact solution.
283  Assuming that the iteration does converge with a convergence rate $\chi^{-}$ we have the estimate  Assuming that the iteration does converge with a convergence rate $\chi^{-}$ we have the estimate
284
285  \max(\|dv\|\hackscore{1}, \|Bv\hackscore{1}\|\hackscore{0})  \max(\|dv\|\hackscore{1}, \|Bv\hackscore{1}\|\hackscore{0})
286  \le \chi^{-} \max(\|dv^{-}\|\hackscore{1}, \|Bv\hackscore{1}^{-}\|\hackscore{0})  \le \chi^{-} \max(\|dv^{-}\|\hackscore{1}, \|Bv\hackscore{1}^{-}\|\hackscore{0})
287  + \max(\|\delta v\|\hackscore{1} + \|(B A^{-1} B^{*}) \delta p\|\hackscore{0})  + \max(\|\delta v\|\hackscore{1}, \|(B A^{-1} B^{*}) \delta p\|\hackscore{0})
288
289  were the upper index $-$ referes to the increments in the last step.  were the upper index $-$ referes to the increments in the last step.
290  Now we try to establish estimates for $\|\delta v\|$ and $\|Bv\hackscore{1}\|$ from
291    Now we will establish estimates for $\|\delta v\|$ and $\|Bv\hackscore{1}\|$ from
292  formulas~\ref{STOKES ERRORS} where neglect the $\delta p$ in the equation for $\delta v$ assuming that  formulas~\ref{STOKES ERRORS} where neglect the $\delta p$ in the equation for $\delta v$ assuming that
293  $\delta p$ is controlled.  $\delta p$ is controlled.
294
# Line 287  where after the first step we set $\gamm Line 405 where after the first step we set$\gamm
405  \subsection{Functions}  \subsection{Functions}
406
407  \begin{classdesc}{StokesProblemCartesian}{domain}  \begin{classdesc}{StokesProblemCartesian}{domain}
408  opens the Stokes problem\index{Stokes problem} on the \Domain domain. The approximation  opens the Stokes problem\index{Stokes problem} on the \Domain domain. The domain
409  order needs to be two. Order two will be used to approximate the velocity and order one  needs to support LBB compliant elements for the Stokes problem, see~\cite{LBB} for detials~\index{LBB condition}.
410  to approximate the pressure.  For instance one can use second order polynomials for velocity and
411    first order polynomials for the pressure on the same element. Alternativly, one can use
412    macro elements\index{macro elements} using linear polynomial for both pressure and velocity bu with a subdivided
413    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
414    \begin{python}
415    velocity=Vector(0.0, Solution(mesh))
416    pressure=Scalar(0.0, ReducedSolution(mesh))
417    \end{python}
418  \end{classdesc}  \end{classdesc}
419
420  \begin{methoddesc}[StokesProblemCartesian]{initialize}{\optional{f=Data(), \optional{fixed_u_mask=Data(), \optional{eta=1, \optional{surface_stress=Data(), \optional{stress=Data()}, \optional{  \begin{methoddesc}[StokesProblemCartesian]{initialize}{\optional{f=Data(), \optional{fixed_u_mask=Data(), \optional{eta=1, \optional{surface_stress=Data(), \optional{stress=Data()}, \optional{
# Line 306  The method will try to cast the given va Line 431  The method will try to cast the given va
431  \begin{methoddesc}[StokesProblemCartesian]{solve}{v,p  \begin{methoddesc}[StokesProblemCartesian]{solve}{v,p
432  \optional{, max_iter=100 \optional{, verbose=False \optional{, usePCG=True }}}}  \optional{, max_iter=100 \optional{, verbose=False \optional{, usePCG=True }}}}
433  solves the problem and return approximations for velocity and pressure.  solves the problem and return approximations for velocity and pressure.
434  The arguments \var{v} and \var{p} define initial guess. The values of \var{v} marked  The arguments \var{v} and \var{p} define initial guess.
435    \var{v} must have function space \var{Solution(domain)} and
436    \var{p} must have function space \var{ReducedSolution(domain)}.
437    The values of \var{v} marked