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

revision 3324 by jfenwick, Fri Oct 22 01:56:02 2010 UTC revision 3325 by caltinay, Fri Oct 29 00:30:51 2010 UTC
# Line 1  Line 1
1
2    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3    %
4    % Copyright (c) 2003-2010 by University of Queensland
5    % Earth Systems Science Computational Center (ESSCC)
6    % http://www.uq.edu.au/esscc
7    %
8    % Primary Business: Queensland, Australia
11    %
12    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
13
14  \section{The Stokes Problem}  \section{The Stokes Problem}
15  \label{STOKES PROBLEM}  \label{STOKES PROBLEM}
16  In this section we discuss how to solve the Stokes problem which is defined as follows:  In this section we discuss how to solve the Stokes problem.
17    We want to calculate the velocity\index{velocity} field $v$ and pressure $p$
18  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}  of an incompressible fluid\index{incompressible fluid}.
19    They are given as the solution of the Stokes problem\index{Stokes problem}
20  \label{Stokes 1}  \label{Stokes 1}
21  -\left(\eta(v_{i,j}+ v_{j,i})\right)_{,j}+p_{,i}=f_{i}-\sigma_{ij,j}  -\left(\eta(v_{i,j}+ v_{j,i})\right)_{,j}+p_{,i}=f_{i}-\sigma_{ij,j}
22
23  where  $f_{i}$ defines an internal force \index{force, internal} and $\sigma_{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.  where  $f_{i}$ defines an internal force\index{force, internal} and
24    $\sigma_{ij}$ is an initial stress\index{stress, initial}.
25    The viscosity $\eta$ may weakly depend on pressure and velocity.
26    If relevant we will use the notation $\eta(v,p)$ to express this dependency.
27
28  We assume an incompressible media:  We assume an incompressible medium:
29  \label{Stokes 2}  \label{Stokes 2}
30  -v_{i,i}=0  -v_{i,i}=0
31
# Line 20  which can be overwritten by constraints Line 37  which can be overwritten by constraints
37  \label{Stokes Boundary0}  \label{Stokes Boundary0}
38  v_{i}(x)=v^D_{i}(x)  v_{i}(x)=v^D_{i}(x)
39
40  at some locations $x$ at the boundary of the domain. $s_{i}$ defines a normal stress and  at some locations $x$ at the boundary of the domain.
41  $\alpha\ge 0$ the spring constant for restoring normal force.  $s_{i}$ defines a normal stress and $\alpha\ge 0$ the spring constant for
42    restoring normal force.
43  The index $i$ may depend on the location $x$ on the boundary.  The index $i$ may depend on the location $x$ on the boundary.
44  $v^D$ is a given function on the domain.  $v^D$ is a given function on the domain.
45
46  \subsection{Solution Method \label{STOKES SOLVE}}  \subsection{Solution Method \label{STOKES SOLVE}}
47  If we assume that $\eta$ is independent from the velocity and pressure equations~\ref{Stokes 1} and~\ref{Stokes 2}  If we assume that $\eta$ is independent from the velocity and pressure,
48  can be written in the block form  equations~\ref{Stokes 1} and~\ref{Stokes 2} can be written in the block form
49
50  \left[ \begin{array}{cc}  \left[ \begin{array}{cc}
51  A     & B^{*} \\  A     & B^{*} \\
# Line 43  G \\ Line 61  G \\
61  \end{array} \right]  \end{array} \right]
62  \label{STOKES}  \label{STOKES}
63
64  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 a coercive, self-adjoint linear operator in a suitable Hilbert
65  For more details on the mathematics see references \cite{AAMIRBERKYAN2008,MBENZI2005}.  space, $B$ is the $(-1) \cdot$ divergence operator and $B^{*}$ is the adjoint
67  If $v_{0}$ and $p_{0}$ are given initial guesses for  For more details on the mathematics see references \cite{AAMIRBERKYAN2008,MBENZI2005}.
68  velocity and pressure we calculate a correction $dv$ for the velocity by solving the first
69  equation of equation~\ref{STOKES}  If $v_{0}$ and $p_{0}$ are given initial guesses for velocity and pressure we
70    calculate a correction $dv$ for the velocity by solving the first equation of
71    \eqn{STOKES}
72   \label{STOKES ITER STEP 1}   \label{STOKES ITER STEP 1}
73   A dv_{1} = G - A v_{0} - B^{*} p_{0}   A dv_{1} = G - A v_{0} - B^{*} p_{0}
74
75  We then insert the new approximation $v_{1}=v_{0}+dv_{1}$ to calculate a correction $dp_{2}$  We then insert the new approximation $v_{1}=v_{0}+dv_{1}$ to calculate a
76  for the pressure and an additional correction $dv_{2}$ for the velocity by solving  correction $dp_{2}$ for the pressure and an additional correction $dv_{2}$ for
77    the velocity by solving
78
79   \begin{array}{rcl}   \begin{array}{rcl}
80   B A^{-1} B^{*} dp_{2} & = & Bv_{1} \\   B A^{-1} B^{*} dp_{2} & = & Bv_{1} \\
# Line 63  for the pressure and an additional corre Line 84  for the pressure and an additional corre
84
85  The new velocity and pressure are then given by $v_{2}=v_{1}-dv_{2}$ and  The new velocity and pressure are then given by $v_{2}=v_{1}-dv_{2}$ and
86  $p_{2}=p_{0}+dp_{2}$ which will fulfill the block system~\ref{STOKES}.  $p_{2}=p_{0}+dp_{2}$ which will fulfill the block system~\ref{STOKES}.
87  This solution strategy is called the Uzawa scheme \index{Uzawa scheme}.  This solution strategy is called the Uzawa scheme\index{Uzawa scheme}.
88
89  There is a problem with this scheme: In practice we will use an iterative scheme  There is a problem with this scheme: in practice we will use an iterative
90  to solve any problem for operator $A$. So we will be unable to calculate the operator  scheme to solve any problem for operator $A$.
91  $B A^{-1} B^{*}$ required for $dp_{2}$ explicitly. In fact, we need to use another  So we will be unable to calculate the operator $B A^{-1} B^{*}$ required for
92  iterative scheme to solve the first equation in~\ref{STOKES ITER STEP 2} where in each iteration step  $dp_{2}$ explicitly. In fact, we need to use another iterative scheme to solve
93    the first equation in~\ref{STOKES ITER STEP 2} where in each iteration step
94  an iterative solver for $A$ is applied. Another issue is the fact that the  an iterative solver for $A$ is applied. Another issue is the fact that the
95  viscosity $\eta$ may depend on velocity or pressure and so we need to iterate over the  viscosity $\eta$ may depend on velocity or pressure and so we need to iterate
96  three equations~\ref{STOKES ITER STEP 1} and~\ref{STOKES ITER STEP 2}.  over the three equations~\ref{STOKES ITER STEP 1} and~\ref{STOKES ITER STEP 2}.
97
98  In the following we will use the two norms  In the following we will use the two norms
99
# Line 80  In the following we will use the two nor Line 102  In the following we will use the two nor
102  \|p\|_{0}^2= \int_{\Omega} p^2 \; dx.  \|p\|_{0}^2= \int_{\Omega} p^2 \; dx.
103  \label{STOKES STOP}  \label{STOKES STOP}
104
105  for velocity $v$ and pressure $p$. The iteration is terminated if the stopping criteria  for velocity $v$ and pressure $p$.
106    The iteration is terminated if the stopping criterion
107   \label{STOKES STOPPING CRITERIA}   \label{STOKES STOPPING CRITERIA}
108  \max(\|Bv_{1}\|_{0},\|v_{2}-v_{0}\|_{1}) \le \tau \cdot \|v_{2}\|_{1}  \max(\|Bv_{1}\|_{0},\|v_{2}-v_{0}\|_{1}) \le \tau \cdot \|v_{2}\|_{1}
109
110   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  for a given tolerance $0<\tau<1$ is met.
111  $\|Bv_{1}\|_{0}$ equals the  Notice that because of the first equation of~\ref{STOKES ITER STEP 2} we have
112  norm of $B A^{-1} B^{*} dp_{2}$ and consequently provides a norm for the pressure correction.  that $\|Bv_{1}\|_{0}$ equals the norm of $B A^{-1} B^{*} dp_{2}$ and
113    consequently provides a norm for the pressure correction.
114
115  We want to optimize the tolerance choice for solving~\ref{STOKES ITER STEP 1}  We want to optimize the tolerance choice for solving~\ref{STOKES ITER STEP 1}
116  and~\ref{STOKES ITER STEP 2}. To do this we write the iteration scheme as a fixed point problem. Here  and~\ref{STOKES ITER STEP 2}. To do this we write the iteration scheme as a
117  we consider the errors produced by the iterative solvers being used.  fixed point problem. Here we consider the errors produced by the iterative
118  From Equation~\ref{STOKES ITER STEP 1} we have  solvers being used.
119    From \eqn{STOKES ITER STEP 1} we have
120   \label{STOKES total V1}   \label{STOKES total V1}
121  v_{1} = e_{1} + v_{0} + A^{-1} ( G - Av_{0} - B^{*} p_{0} )  v_{1} = e_{1} + v_{0} + A^{-1} ( G - Av_{0} - B^{*} p_{0} )
122
123  where $e_{1}$ is the error when solving~\ref{STOKES ITER STEP 1}.    where $e_{1}$ is the error when solving~\ref{STOKES ITER STEP 1}.
124  We will use a sparse matrix solver so we have not full control on the norm $\|.\|_{s}$ used in the stopping criteria for this equation. In fact we will have a stopping criteria of the form  We will use a sparse matrix solver so we have not full control on the norm
125    $\|.\|_{s}$ used in the stopping criterion for this equation.
126    In fact we will have a stopping criterion of the form
127
128  \| A e_{1} \|_{s}  = \| G - A v_{1} - B^{*} p_{0} \|_{s} \le \tau_{1} \| G - A v_{0} - B^{*} p_{0} \|_{s}  \| A e_{1} \|_{s}  = \| G - A v_{1} - B^{*} p_{0} \|_{s} \le \tau_{1} \| G - A v_{0} - B^{*} p_{0} \|_{s}
129
130  where $\tau_{1}$ is the tolerance which we need to choose. This translates into the condition  where $\tau_{1}$ is the tolerance which we need to choose.
131    This translates into the condition
132
133  \| e_{1} \|_{1} \le K \tau_{1} \| dv_{1} - e_{1} \|_{1}  \| e_{1} \|_{1} \le K \tau_{1} \| dv_{1} - e_{1} \|_{1}
134
135  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
136  norm being used and the condition number of the stiffness matrix.    factors such as the norm being used and the condition number of the stiffness matrix.
137  From the first equation of~\ref{STOKES ITER STEP 2} we have  From the first equation of~\ref{STOKES ITER STEP 2} we have
138  \label{STOKES total P2}  \label{STOKES total P2}
139  p_{2} =  p_{0} + (B A^{-1} B^{*})^{-1} (e_{2} + Bv_{1} )  p_{2} =  p_{0} + (B A^{-1} B^{*})^{-1} (e_{2} + Bv_{1} )
140
141  where $e_{2}$ represents the error when solving~\ref{STOKES ITER STEP 2}.  where $e_{2}$ represents the error when solving~\ref{STOKES ITER STEP 2}.
142  We use an iterative preconditioned conjugate gradient method (PCG) \index{linear solver!PCG}\index{PCG} with iteration  We use an iterative preconditioned conjugate gradient method
143  operator $B A^{-1} B^{*}$ using the $\|.\|_{0}$ norm. As suitable preconditioner \index{preconditioner} for the iteration  (PCG)\index{linear solver!PCG}\index{PCG} with iteration operator
144  operator we use $\frac{1}{\eta}$ \cite{ELMAN}, ie  $B A^{-1} B^{*}$ using the $\|.\|_{0}$ norm.
145  the evaluation of the preconditioner $P$ for a given pressure increment $q$ is the solution of  As suitable preconditioner\index{preconditioner} for the iteration operator we
146    use $\frac{1}{\eta}$ \cite{ELMAN}, i.e. the evaluation of the preconditioner
147    $P$ for a given pressure increment $q$ is the solution of
148   \label{STOKES P PREC}   \label{STOKES P PREC}
149  \frac{1}{\eta} (Pq) = q \; .  \frac{1}{\eta} (Pq) = q \; .
150
151  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$
152  the problem  one needs to solve the problem
153   \label{STOKES P OPERATOR}   \label{STOKES P OPERATOR}
154  A w = B^{*} s  A w = B^{*} s
155
156  with sufficient accuracy to return $q=Bw$. We assume that the desired tolerance is  with sufficient accuracy to return $q=Bw$. We assume that the desired
157  sufficiently small, for instance one can take $\tau_{2}^2$  tolerance is sufficiently small, for instance one can take $\tau_{2}^2$ where
158  where $\tau_{2}$ is the tolerance for~\ref{STOKES ITER STEP 2}.  $\tau_{2}$ is the tolerance for~\ref{STOKES ITER STEP 2}.
159
160  In an implementation we use the fact that the residual $r$ is given as  In an implementation we use the fact that the residual $r$ is given as
161   \label{STOKES RES }   \label{STOKES RES }
162   r= B (v_{1} -  A^{-1} B^{*} dp) =  B (v_{1} - A^{-1} B^{*} dp) = B (v_{1}-dv_{2}) = B v_{2}   r= B (v_{1} -  A^{-1} B^{*} dp) =  B (v_{1} - A^{-1} B^{*} dp) = B (v_{1}-dv_{2}) = B v_{2}
163
164  In particular we have $e_{2} = B v_{2}$  In particular we have $e_{2} = B v_{2}$.
165  So the residual $r$ is represented by the updated velocity $v_{2}=v_{1}-dv_{2}$. In practice, if  So the residual $r$ is represented by the updated velocity $v_{2}=v_{1}-dv_{2}$.
166  one uses the velocity to represent the residual $r$ there is no need  In practice, if one uses the velocity to represent the residual $r$ there is
167  to recover $dv_{2}$ in~\ref{STOKES ITER STEP 2} after $dp_{2}$ has been calculated.  no need to recover $dv_{2}$ in~\ref{STOKES ITER STEP 2} after $dp_{2}$ has
168    been calculated.
169  In PCG the iteration is terminated if  In PCG the iteration is terminated if
170   \label{STOKES P OPERATOR ERROR}   \label{STOKES P OPERATOR ERROR}
171  \| P^{\frac{1}{2}}B v_{2} \|_{0} \le \tau_{2} \| P^{\frac{1}{2}}B v_{1} \|_{0}  \| P^{\frac{1}{2}}B v_{2} \|_{0} \le \tau_{2} \| P^{\frac{1}{2}}B v_{1} \|_{0}
# Line 143  where $\tau_{2}$ is the given tolerance. Line 174  where $\tau_{2}$ is the given tolerance.
174   \label{STOKES P OPERATOR ERROR 2}   \label{STOKES P OPERATOR ERROR 2}
175  \|e_{2}\|_{0} = \| B v_{2} \|_{0} \le M \tau_{2} \| B v_{1} \|_{0}  \|e_{2}\|_{0} = \| B v_{2} \|_{0} \le M \tau_{2} \| B v_{1} \|_{0}
176
177  where $M$ is taking care of the fact that $P^{\frac{1}{2}}$ is dropped.    where $M$ is taking care of the fact that $P^{\frac{1}{2}}$ is dropped.
178
179  As we assume that there is no significant error from solving with the operator $A$ we have  As we assume that there is no significant error from solving with the operator
180    $A$ we have
181   \label{STOKES total V2}   \label{STOKES total V2}
182  v_{2} =  v_{1} - dv_{2}  v_{2} =  v_{1} - dv_{2}
183  = v_{1}  - A^{-1} B^{*}dp  = v_{1}  - A^{-1} B^{*}dp
184
185  Combining the equations~\ref{STOKES total V1},~\ref{STOKES total P2} and~\ref{STOKES total V2} and  Combining the equations~\ref{STOKES total V1},~\ref{STOKES total P2} and~\ref{STOKES total V2}
186  setting the errors to zero we can write the solution process as a fix point problem  and setting the errors to zero we can write the solution process as a fix
187    point problem
188
189  v = \Phi(v,p) \mbox{ and } p = \Psi(u,p)  v = \Phi(v,p) \mbox{ and } p = \Psi(u,p)
190
191  with suitable functions $\Phi(v,p)$ and $\Psi(v,p)$ representing the iteration operator without  with suitable functions $\Phi(v,p)$ and $\Psi(v,p)$ representing the
192  errors. In fact for a linear problem,  $\Phi$ and $\Psi$ are constant. With this notation we can write  iteration operator without errors. In fact, for a linear problem, $\Phi$ and
193  the update step in the form $p_{2}= \delta p + \Psi(v_{0},p_{0})$ and  $\Psi$ are constant. With this notation we can write the update step in the
194  $v_{2}= \delta v + \Phi(v_{0},p_{0})$ where  form $p_{2}= \delta p + \Psi(v_{0},p_{0})$ and $v_{2}= \delta v + \Phi(v_{0},p_{0})$
195  the total error $\delta p$ and $\delta v$ are given as  where the total error $\delta p$ and $\delta v$ are given as
196
197   \begin{array}{rcl}   \begin{array}{rcl}
198  \delta p & = &  (B A^{-1} B^{*})^{-1} ( e_{2} + B e_{1} ) \\  \delta p & = &  (B A^{-1} B^{*})^{-1} ( e_{2} + B e_{1} ) \\
199  \delta v & = &  e_{1} -  A^{-1} B^{*}\delta p  \;.  \delta v & = &  e_{1} -  A^{-1} B^{*}\delta p  \;.
200  \end{array}\label{STOKES ERRORS}  \end{array}\label{STOKES ERRORS}
201
202  Notice that $B\delta v = - e_{2}=-Bv_{2}$. Our task is now to choose the tolerances  Notice that $B\delta v = - e_{2}=-Bv_{2}$.
203  $\tau_{1}$ and $\tau_{2}$ such that the global errors $\delta p$ and $\delta v$  Our task is now to choose the tolerances $\tau_{1}$ and $\tau_{2}$ such that
204  do not stop the convergence of the iteration process.  the global errors $\delta p$ and $\delta v$ do not stop the convergence of the
205    iteration process.
206
207  To measure convergence we use  To measure convergence we use
208
# Line 180  $v_{0}$ we will use the estimate Line 214  $v_{0}$ we will use the estimate
214
215  \epsilon = \max(\|v_{2}-v_{0}\|_{1}, \|B v_{1}\|_{0})  \epsilon = \max(\|v_{2}-v_{0}\|_{1}, \|B v_{1}\|_{0})
216
217  to estimate the progress of the iteration step after the step is completed.  to estimate the progress of the iteration step after the step is completed.
218  Note that the estimate of $\epsilon$    Note that the estimate of $\epsilon$ is used in the stopping
219  used in the stopping criteria~\ref{STOKES STOPPING CRITERIA}. If $\chi^{-}$ is the convergence rate assuming  criterion~\ref{STOKES STOPPING CRITERIA}.
220  exact calculations, i.e. $e_{1}=0$ and $e_{2}=0$, we are expecting  If $\chi^{-}$ is the convergence rate assuming exact calculations, i.e.
221  to maintain $\epsilon \le \chi^{-} \cdot \epsilon^{-}$. For the  $e_{1}=0$ and $e_{2}=0$, we are expecting to maintain $\epsilon \le \chi^{-} \cdot \epsilon^{-}$.
222  pressure increment we get:  For the pressure increment we get:
223   \label{STOKES EST 1}   \label{STOKES EST 1}
224  \begin{array}{rcl}  \begin{array}{rcl}
225  \|B A^{-1} B^{*} (p_{2}-p)\|_{0}  \|B A^{-1} B^{*} (p_{2}-p)\|_{0}
# Line 196  pressure increment we get: Line 230  pressure increment we get:
230  & \le & \chi^{-} \cdot \epsilon^{-} + M \tau_{2} \|B v_{1}\|_{0} \\    & \le & \chi^{-} \cdot \epsilon^{-} + M \tau_{2} \|B v_{1}\|_{0} \\
231  \end{array}  \end{array}
232
233  So we choose the value for $\tau_{2}$ from  So we choose the value for $\tau_{2}$ from
234   \label{STOKES TOL2}   \label{STOKES TOL2}
235   M \tau_{2} \|B v_{1}\|_{0}  \le (\chi^{-})^2 \epsilon^{-}   M \tau_{2} \|B v_{1}\|_{0}  \le (\chi^{-})^2 \epsilon^{-}
236
237  in order to make the perturbation for the termination of the pressure iteration a second order effect. We use a  in order to make the perturbation for the termination of the pressure
238  similar argument for the velocity:  iteration a second order effect. We use a similar argument for the velocity:
239  \label{STOKES EST 2}  \label{STOKES EST 2}
240  \begin{array}{rcl}  \begin{array}{rcl}
241  \|v_{2}-v\|_{1} & \le & \|v_{2}-\delta v-v\|_{1} + \| \delta v\|_{1} \\  \|v_{2}-v\|_{1} & \le & \|v_{2}-\delta v-v\|_{1} + \| \delta v\|_{1} \\
# Line 216  So we choose the value for $\tau_{1}$ fr Line 250  So we choose the value for $\tau_{1}$ fr
250   \label{STOKES TOL1}   \label{STOKES TOL1}
251  K \tau_{1} \le \chi^{-}  K \tau_{1} \le \chi^{-}
252
253  Assuming we have estimates for $M$ and $K$ (if none is available, we use the value $1$.)  Assuming we have estimates for $M$ and $K$\footnote{if no estimates are
254  we can use~\ref{STOKES TOL1} and~\ref{STOKES TOL2} to get appropriate values for the tolerances. After  available, we use the value $1$} we can use~\ref{STOKES TOL1} and
255  the step has been completed we can calculate a new convergence rate $\chi =\frac{\epsilon}{\epsilon^{-}}$.  \ref{STOKES TOL2} to get appropriate values for the tolerances.
256  For partial reasons we restrict $\chi$ to be less or equal a given maximum value $\chi_{max}\le 1$.  After the step has been completed we can calculate a new convergence rate
257  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  $\chi =\frac{\epsilon}{\epsilon^{-}}$.
258    For partial reasons we restrict $\chi$ to be less or equal a given maximum
259    value $\chi_{max}\le 1$.
260    If we see $\chi \le \chi^{-} (1+\chi^{-})$ our choices for the tolerances were
261    suitable. Otherwise, we need to adjust the values for $K$ and $M$.
262    From the estimates~\ref{STOKES EST 1} and~\ref{STOKES EST 2} we establish
263  \label{STOKES EST 3}  \label{STOKES EST 3}
264  \chi \le ( 1 + \max(M \frac{\tau_{2} \|B v_{1}\|_{0}}{\chi^{-} \epsilon^{-}},K \tau_{1}  ) ) \cdot \chi^{-}  \chi \le ( 1 + \max(M \frac{\tau_{2} \|B v_{1}\|_{0}}{\chi^{-} \epsilon^{-}},K \tau_{1}  ) ) \cdot \chi^{-}
265
266  If we assume that this inequality would be an equation if we would have chosen the right values  If we assume that this inequality would be an equation if we would have chosen
267  $M^{+}$ and $K^{+}$ then we get  the right values $M^{+}$ and $K^{+}$ then we get
268  \label{STOKES EST 3b}  \label{STOKES EST 3b}
269  \chi =  ( 1 + \max(M^{+} \frac{\chi^{-}}{M},K^{+} \frac{\chi^{-}}{K}) ) \cdot \chi^{-}  \chi =  ( 1 + \max(M^{+} \frac{\chi^{-}}{M},K^{+} \frac{\chi^{-}}{K}) ) \cdot \chi^{-}
270
271  From this equation we set for  From this equation we see if our choice for $K$ was not good enough.
272  than our choice for $K$ was not good enough. In this case we can calculate a new value  In this case we can calculate a new value
273
274  K^{+} =  \frac{\chi-\chi^{-}}{(\chi^{-})^2} K  K^{+} =  \frac{\chi-\chi^{-}}{(\chi^{-})^2} K
275
276  In practice we will use  In practice we will use
277
278  K^{+}  = \max(\frac{\chi-\chi^{-}}{(\chi^{-})^2} K,\frac{1}{2}K,1)  K^{+}  = \max(\frac{\chi-\chi^{-}}{(\chi^{-})^2} K,\frac{1}{2}K,1)
279
280  where the second term is used to reduce a potential overestimate of $K$.    where the second term is used to reduce a potential overestimate of $K$.
281  The same identity is used for to update $M$. The updated $M^{+}$ and $K^{+}$  The same identity is used for to update $M$. The updated $M^{+}$ and $K^{+}$
282  are then use in the next iteration step to control the tolerances.  are then use in the next iteration step to control the tolerances.
283
284  In some cases one can observe that there is a significant change  In some cases one can observe that there is a significant change in the
285  in the velocity but the new velocity $v_{1}$ has still a  velocity but the new velocity $v_{1}$ has still a small divergence, i.e.
286  small divergence, i.e. we have  we have $\|Bv_{1}\|_{0} \ll \|v_{1}-v_{0}\|_{1}$.
287  $\|Bv_{1}\|_{0} \ll \|v_{1}-v_{0}\|_{1}$.  In this case we will get a small pressure increment and consequently only very
288  In this case we will get a small pressure increment and consequently only very small changes to  small changes to the velocity as a result of the second update step which
289  the velocity as a result of the second update step which therefore can be skipped and  therefore can be skipped and we can directly repeat the first update step
290  we can directly repeat the first update step until the increment in velocity becomes  until the increment in velocity becomes significant relative to its divergence.
291  significant relative to its divergence. In practice we will ignore the second half of the iteration step  In practice we will ignore the second half of the iteration step as long as
as long as
292   \label{STOKES LARGE BV1}   \label{STOKES LARGE BV1}
293  \|Bv_{1}\|_{0} \le \theta \cdot \|v_{1}-v_{0}\|  \|Bv_{1}\|_{0} \le \theta \cdot \|v_{1}-v_{0}\|
294
295  where $0<\theta<1$ is a given factor. In this case we will also check the stopping criteria  where $0<\theta<1$ is a given factor. In this case we will also check the
296  with $v_{1}\rightarrow v_{2}$ but we will not correct $M$ in this case.  stopping criterion with $v_{1}\rightarrow v_{2}$ but we will not correct $M$
297    in this case.
298
299  Starting from an initial guess $v_{0}$ and $p_{0}$ for velocity and pressure  Starting from an initial guess $v_{0}$ and $p_{0}$ for velocity and pressure
300  the solution procedure is implemented as follows.  the solution procedure is implemented as follows:
301  \begin{enumerate}  \begin{enumerate}
302   \item calculate viscosity $\eta(v_{0},p)_{0}$ and assemble operator $A$ from $\eta$.   \item calculate viscosity $\eta(v_{0},p)_{0}$ and assemble operator $A$ from $\eta$
303   \item calculate the tolerance $\tau_{1}$ from~\ref{STOKES TOL1}.   \item calculate the tolerance $\tau_{1}$ from \eqn{STOKES TOL1}
304   \item Solve~\ref{STOKES ITER STEP 1} for $dv_{1}$ with tolerance $\tau_{1}$.   \item solve \eqn{STOKES ITER STEP 1} for $dv_{1}$ with tolerance $\tau_{1}$
305   \item update $v_{1}= v_{0}+ dv_{1}$   \item update $v_{1}= v_{0}+ dv_{1}$
306   \item if $Bv_{1}$ is large (see~\ref{STOKES LARGE BV1}):   \item if $Bv_{1}$ is large (see~\ref{STOKES LARGE BV1}):
307   \begin{enumerate}   \begin{enumerate}
308   \item calculate the tolerance $\tau_{2}$ from~\ref{STOKES TOL2}.     \item calculate the tolerance $\tau_{2}$ from~\ref{STOKES TOL2}
309   \item Solve~\ref{STOKES ITER STEP 2} for $dp_{2}$ and $v_{2}$ with tolerance $\tau_{2}$.     \item solve~\ref{STOKES ITER STEP 2} for $dp_{2}$ and $v_{2}$ with tolerance $\tau_{2}$
310   \item update $p_{2}\leftarrow p_{0}+ dp_{2}$     \item update $p_{2}\leftarrow p_{0}+ dp_{2}$
311   \end{enumerate}   \end{enumerate}
312   \item else:   \item else:
313    \begin{enumerate}    \begin{itemize}
314    \item update $p_{2}\leftarrow p$ and $v_{2}\leftarrow v_{1}$      \item update $p_{2}\leftarrow p$ and $v_{2}\leftarrow v_{1}$
315     \end{enumerate}    \end{itemize}
316     \item calculate convergence measure $\epsilon$ and convergence rate $\chi$   \item calculate convergence measure $\epsilon$ and convergence rate $\chi$
317  \item if stopping criteria~\ref{STOKES STOPPING CRITERIA} holds:   \item if stopping criterion~\ref{STOKES STOPPING CRITERIA} holds:
318   \begin{enumerate}   \begin{itemize}
319   \item return $v_{2}$ and $p_{2}$     \item return $v_{2}$ and $p_{2}$
320   \end{enumerate}   \end{itemize}
321   \item else:   \item else:
322   \begin{enumerate}   \begin{enumerate}
323       \item update $M$ and $K$
324       \item update $M$ and $K$     \item goto step 1 with $v_{0}\leftarrow v_{2}$ and $p_{0}\leftarrow p_{2}$.
325       \item goto step 1 with $v_{0}\leftarrow v_{2}$ and $p_{0}\leftarrow p_{2}$.   \end{enumerate}
\end{enumerate}
326  \end{enumerate}  \end{enumerate}
327
328  \subsection{Functions}  \subsection{Functions}
329
330  \begin{classdesc}{StokesProblemCartesian}{domain}  \begin{classdesc}{StokesProblemCartesian}{domain}
331  opens the Stokes problem\index{Stokes problem} on the \Domain domain. The domain  opens the Stokes problem\index{Stokes problem} on the \Domain domain.
332  needs to support LBB compliant elements for the Stokes problem, see~\cite{LBB} for details~\index{LBB condition}.  The domain needs to support LBB compliant elements for the Stokes problem, see~\cite{LBB} for details\index{LBB condition}.
333  For instance one can use second order polynomials for velocity and  For instance one can use second order polynomials for velocity and first order
334  first order polynomials for the pressure on the same element. Alternatively, one can use  polynomials for the pressure on the same element.
335  macro elements\index{macro elements} using linear polynomial for both pressure and velocity    with a subdivided  Alternatively, one can use macro elements\index{macro elements} using linear
336  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  polynomials for both pressure and velocity with a subdivided element for the
337    velocity. Typically, the macro element is more cost effective.
338    The fact that pressure and velocity are represented in different ways is
339    expressed by
340  \begin{python}  \begin{python}
341  velocity=Vector(0.0, Solution(mesh))    velocity=Vector(0.0, Solution(mesh))
342  pressure=Scalar(0.0, ReducedSolution(mesh))    pressure=Scalar(0.0, ReducedSolution(mesh))
343  \end{python}  \end{python}
344  \end{classdesc}  \end{classdesc}
345
346  \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(),
348        \optional{surface_stress=Data(), \optional{stress=Data()},%
349        \optional{restoration_factor=0}}}}}}
350  assigns values to the model parameters. In any call all values must be set.  assigns values to the model parameters. In any call all values must be set.
351  \var{f} defines the external force $f$, \var{eta} the viscosity $\eta$,  \var{f} defines the external force $f$, \var{eta} the viscosity $\eta$,
352  \var{surface_stress} the surface stress $s$ and \var{stress} the initial stress $\sigma$.  \var{surface_stress} the surface stress $s$ and \var{stress} the initial stress $\sigma$.
353  The locations and components where the velocity is fixed are set by  The locations and components where the velocity is fixed are set by the values
354  the values of \var{fixed_u_mask}. \var{restoration_factor} defines the restoring force factor $\alpha$.  of \var{fixed_u_mask}.
355  The method will try to cast the given values to appropriate  \var{restoration_factor} defines the restoring force factor $\alpha$.
356  \Data class objects.  The method will try to cast the given values to appropriate \Data class objects.
357  \end{methoddesc}  \end{methoddesc}
358
359  \begin{methoddesc}[StokesProblemCartesian]{solve}{v,p  \begin{methoddesc}[StokesProblemCartesian]{solve}{v,p
360  \optional{, max_iter=100 \optional{, verbose=False \optional{, usePCG=True }}}}  \optional{, max_iter=100 \optional{, verbose=False \optional{, usePCG=True }}}}
361  solves the problem and return approximations for velocity and pressure.  solves the problem and returns approximations for velocity and pressure.
362  The arguments \var{v} and \var{p} define initial guess.  The arguments \var{v} and \var{p} define initial guesses.
363  \var{v} must have function space \var{Solution(domain)} and  \var{v} must have function space \var{Solution(domain)} and \var{p} must have
364  \var{p} must have function space \var{ReducedSolution(domain)}.  function space \var{ReducedSolution(domain)}.
365  The values of \var{v} marked  The values of \var{v} marked by \var{fixed_u_mask} remain unchanged.
366  by \var{fixed_u_mask} remain unchanged.  If \var{usePCG} is set to \True then the preconditioned conjugate gradient
367  If \var{usePCG} is set to \True  method (PCG)\index{preconditioned conjugate gradient method!PCG} scheme is
368  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  used. Otherwise the problem is solved with the generalized minimal residual
369  the PCG scheme is more efficient.  method (GMRES)\index{generalized minimal residual method!GMRES}.
370  \var{max_iter} defines the maximum number of iteration steps.  In most cases the PCG scheme is more efficient.
371    \var{max_iter} defines the maximum number of iteration steps.
372  If \var{verbose} is set to \True informations on the progress of of the solver are printed.  If \var{verbose} is set to \True information on the progress of of the solver
373    is printed.
374  \end{methoddesc}  \end{methoddesc}
375

376  \begin{methoddesc}[StokesProblemCartesian]{setTolerance}{\optional{tolerance=1.e-4}}  \begin{methoddesc}[StokesProblemCartesian]{setTolerance}{\optional{tolerance=1.e-4}}
377  sets the tolerance in an appropriate norm relative to the right hand side. The tolerance must be non-negative and less than 1.  sets the tolerance in an appropriate norm relative to the right hand side.
378    The tolerance must be non-negative and less than 1.
379  \end{methoddesc}  \end{methoddesc}
380
381  \begin{methoddesc}[StokesProblemCartesian]{getTolerance}{}  \begin{methoddesc}[StokesProblemCartesian]{getTolerance}{}
382  returns the current relative tolerance.  returns the current relative tolerance.
383  \end{methoddesc}  \end{methoddesc}
384
385  \begin{methoddesc}[StokesProblemCartesian]{setAbsoluteTolerance}{\optional{tolerance=0.}}  \begin{methoddesc}[StokesProblemCartesian]{setAbsoluteTolerance}{\optional{tolerance=0.}}
386  sets the absolute tolerance for the error in the relevant norm. The tolerance must be non-negative. Typically the  sets the absolute tolerance for the error in the relevant norm.
387  absolute tolerance is set to 0.  The tolerance must be non-negative.
388    Typically the absolute tolerance is set to 0.
389  \end{methoddesc}  \end{methoddesc}
390
391  \begin{methoddesc}[StokesProblemCartesian]{getAbsoluteTolerance}{}  \begin{methoddesc}[StokesProblemCartesian]{getAbsoluteTolerance}{}
# Line 346  returns the current absolute tolerance. Line 393  returns the current absolute tolerance.
393  \end{methoddesc}  \end{methoddesc}
394
395  \begin{methoddesc}[StokesProblemCartesian]{getSolverOptionsVelocity}{}  \begin{methoddesc}[StokesProblemCartesian]{getSolverOptionsVelocity}{}
396  returns the solver options used  solve the equations~(\ref{STOKES ITER STEP 1}) and~(\ref{STOKES P OPERATOR}) for velocity.  returns the solver options used to solve \eqn{STOKES ITER STEP 1} and \eqn{STOKES P OPERATOR}) for velocity.
397  \end{methoddesc}  \end{methoddesc}
398
399  \begin{methoddesc}[StokesProblemCartesian]{getSolverOptionsPressure}{}  \begin{methoddesc}[StokesProblemCartesian]{getSolverOptionsPressure}{}
400  returns the solver options used solve the preconditioner equation~(\ref{STOKES P PREC}) for pressure.  returns the solver options used to solve the preconditioner \eqn{STOKES P PREC} for pressure.
401  \end{methoddesc}  \end{methoddesc}
402
403  \begin{methoddesc}[StokesProblemCartesian]{getSolverOptionsDiv}{}  \begin{methoddesc}[StokesProblemCartesian]{getSolverOptionsDiv}{}
404  set the solver options for solving the equation to project the divergence of the velocity onto the function space of pressure.  sets the solver options for solving the equation to project the divergence of
405    the velocity onto the function space of pressure.
406  \end{methoddesc}  \end{methoddesc}
407
408    \subsection{Example: Lid-driven Cavity}
409  \subsection{Example: Lit Driven Cavity}  The following script \file{lid_driven_cavity.py}\index{scripts!\file{lid_driven_cavity.py}}
410   The following script \file{lit_driven_cavity.py}  which is available in the \ExampleDirectory illustrates the usage of the
411  \index{scripts!\file{helmholtz.py}} which is available in the \ExampleDirectory  \class{StokesProblemCartesian} class to solve the lid-driven cavity problem:
illustrates the usage of the \class{StokesProblemCartesian} class to solve
the lit driven cavity problem:
412  \begin{python}  \begin{python}
413  from esys.escript import *    from esys.escript import *
414  from esys.finley import Rectangle    from esys.finley import Rectangle
415  from esys.escript.models import StokesProblemCartesian    from esys.escript.models import StokesProblemCartesian
416  NE=25    NE=25
417  dom = Rectangle(NE,NE,order=2)    dom = Rectangle(NE,NE,order=2)
418  x = dom.getX()    x = dom.getX()
419  sc=StokesProblemCartesian(dom)    sc=StokesProblemCartesian(dom)
421        (whereZero(x[1])*[0.,1.]+whereZero(x[1]-1))*[1.,1]          (whereZero(x[1])*[0.,1.]+whereZero(x[1]-1))*[1.,1]