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

Revision 2793 - (show annotations)
Tue Dec 1 06:10:10 2009 UTC (9 years, 11 months ago) by gross
File MIME type: application/x-tex
File size: 25596 byte(s)
some new util functions

 1 \section{The Stokes Problem} 2 \label{STOKES PROBLEM} 3 In this section we discuss how to solve the Stokes problem which is defined as follows: 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} 6 \begin{equation}\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} 8 \end{equation} 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. 10 11 We assume an incompressible media: 12 \begin{equation}\label{Stokes 2} 13 -v\hackscore{i,i}=0 14 \end{equation} 15 Natural boundary conditions are taken in the form 16 \begin{equation}\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} 18 \end{equation} 19 which can be overwritten by constraints of the form 20 \begin{equation}\label{Stokes Boundary0} 21 v\hackscore{i}(x)=v^D\hackscore{i}(x) 22 \end{equation} 23 at some locations $x$ at the boundary of the domain. $s\hackscore{i}$ defines a normal stress and 24 $\alpha\ge 0$ the spring constant for restoring normal force. 25 The index $i$ may depend on the location $x$ on the boundary. 26 $v^D$ is a given function on the domain. 27 28 \subsection{Solution Method \label{STOKES SOLVE}} 29 If we assume that $\eta$ is independent from the velocity and pressure equations~\ref{Stokes 1} and~\ref{Stokes 2} 30 can be written in the block form 31 \begin{equation} 32 \left[ \begin{array}{cc} 33 A & B^{*} \\ 34 B & 0 \\ 35 \end{array} \right] 36 \left[ \begin{array}{c} 37 v \\ 38 p \\ 39 \end{array} \right] 40 =\left[ \begin{array}{c} 41 G \\ 42 0 \\ 43 \end{array} \right] 44 \label{SADDLEPOINT} 45 \end{equation} 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). 47 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 51 equation of equation~\ref{SADDLEPOINT} 52 \begin{equation}\label{SADDLEPOINT ITER STEP 1} 53 A dv\hackscore{1} = G - A v\hackscore{0} - B^{*} p\hackscore{0} 54 \end{equation} 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 \begin{equation} 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} 62 \label{SADDLEPOINT ITER STEP 2} 63 \end{equation} 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 \begin{equation} 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 \end{equation} 81 to define the stopping criterium 82 \begin{equation} 83 \max(\|Bv\hackscore{1}\|\hackscore{0},\|v-v\hackscore{0}\|\hackscore{1}) \le \tau \cdot \|v\|\hackscore{1} 84 \end{equation} 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 \begin{equation} 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 \end{equation} 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 \begin{equation} 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 \end{equation} 110 If $e\hackscore{1}$ is the error of the returned solution $dv\hackscore{1}$ this translates into the condition 111 \begin{equation} 112 \| e\hackscore{1} \|\hackscore{1} \le K \tau\hackscore{1} \| dv\hackscore{1} - e\hackscore{1} \|\hackscore{1} 113 \end{equation} 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 \begin{equation} 117 \| dv\hackscore{1} \|\hackscore{1} \le (1+K \tau\hackscore{1}) \chi \| dv\hackscore{1}^- \|\hackscore{1} 118 \label{DV COND} 119 \end{equation} 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 \begin{equation} 123 \tau\hackscore{1} \le \frac{\min(\chi,\frac{1}{2})}{K} 124 \label{NEW TAU 1} 125 \end{equation} 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 \begin{equation} 130 \frac{1}{K} = \frac{\chi^-}{\chi-\chi^{-}} \tau\hackscore{1}^- 131 \end{equation} 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 \begin{equation} 145 v\hackscore{1} = e\hackscore{1} + A^{-1} ( G - B^{*} p ) 146 \end{equation} 147 This translates into the conditoon 148 \begin{equation} 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 \end{equation} 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 164 velocity and pressure we perform the following steps in the Uzawa scheme \index{Uzawa scheme} style: 165 \begin{enumerate} 166 \item calculate viscosity $\eta(v,p)$ and assemble operator $A$ from $\eta$. 167 \item Solve for $dv$: 168 \begin{equation} 169 A dv = G - A v - B^{*} p \label{SADDLEPOINT ITER STEP 1} 170 \end{equation} 171 \item update $v\hackscore{1}= v+ dv$ 172 \item if $\max(\|Bv\hackscore{1}\|\hackscore{0},\|dv\|\hackscore{1}) \le \tau \cdot \|v\hackscore{1}\|\hackscore{1}$: return v$\hackscore{1},p$ 173 \item Solve for $dp$: 174 \begin{equation} 175 \begin{array}{rcl} 176 B A^{-1} B^{*} dp & = & Bv\hackscore{1} \\ 177 A dv\hackscore{2} & = & B^{*} dp 178 \end{array} 179 \label{SADDLEPOINT ITER STEP 2} 180 \end{equation} 181 \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}$. 183 \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. 185 In this algorithm 186 \begin{equation} 187 \|v\|\hackscore{1}^2 = \int\hackscore{\Omega} v\hackscore{j,k}v\hackscore{j,k} \; dx 188 \mbox{ and } 189 \|p\|\hackscore{0}^2= \int\hackscore{\Omega} p^2 \; dx. 190 \label{STOKES STOP} 191 \end{equation} 192 so the stopping criterion used is checking for convergence as well as for 193 the fact that the incompressiblity condition~\ref{Stokes 2} is fullfilled. 194 195 To solve the update step~\ref{SADDLEPOINT ITER STEP 2} we use an iterative solver with iteration 196 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 197 the evaluation of the preconditioner $P$ for a given pressure increment $q$ is the solution of 198 \begin{equation} \label{P PREC} 199 \frac{1}{\eta} Pq = q \; . 200 \end{equation} 201 Note that in each evaluation of the iteration operator $q=B A^{-1} B^{*} s$ one needs to solve 202 the problem 203 \begin{equation} \label{P OPERATOR} 204 A w = B^{*} s 205 \end{equation} 206 with sufficient accuracy to return $q=Bw$. Notice that the residual $r$ is given as 207 \begin{equation} \label{STOKES RES } 208 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} 209 \end{equation} 210 so in fact the residual $r$ is represented by the updated velocity $v\hackscore{2}$. This saves the recovery of 211 $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 212 \begin{equation} \label{P OPERATOR} 213 \| P^{\frac{1}{2}}B v\hackscore{2} \|\hackscore{0} \le \tau\hackscore{2} \| P^{\frac{1}{2}}B v\hackscore{1} \|\hackscore{0} 214 \end{equation} 215 where $\tau\hackscore{2}$ is the given tolerane. The update step~\ref{P OPERATOR} involves the 216 solution of a sparse matrix problem. We use the tolerance $\tau\hackscore{2}^2$ in order to make sure that any 217 error from solving this problem does not interfere with the PCG iteration. 218 219 220 221 We need to think about appropriate stopping criteria when solving 222 \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 223 convergence of the iteration process one level above. From Equation~\ref{SADDLEPOINT ITER STEP 1} we have 224 \begin{equation} 225 v\hackscore{1} = e\hackscore{1} + A^{-1} ( G - B^{*} p ) 226 \end{equation} 227 We will use a sparse matrix solver so we have not full control on the norm $\|.\|\hackscore{s}$ used in the stopping criteria 228 \begin{equation} 229 \| G - A v\hackscore{1} - B^{*} p \|\hackscore{s} \le \tau\hackscore{1} \| G - A v - B^{*} p \|\hackscore{s} 230 \end{equation} 231 This translates into the conditoon 232 \begin{equation} 233 \| e\hackscore{1} \|\hackscore{1} \le K \tau\hackscore{1} \| dv\hackscore{1} - e\hackscore{1} \|\hackscore{1} 234 \mbox{ therefore } 235 \| e\hackscore{1} \|\hackscore{1} \le \frac{K \tau\hackscore{1}}{1-K \tau\hackscore{1}} \| dv\hackscore{1} \|\hackscore{1} 236 \end{equation} 237 The constant $K$ represents some uncertainty combining a variety of unknown factors such as the 238 solver being used and the condition number of the stiffness matrix. 239 240 From the first equation of~\ref{SADDLEPOINT ITER STEP 2} we have 241 \begin{equation} 242 p\hackscore{2} = p + (B A^{-1} B^{*})^{-1} (e\hackscore{2} + Bv\hackscore{1} ) = 243 (B A^{-1} B^{*})^{-1} ( e\hackscore{2} + B e\hackscore{1} + B A^{-1} G) 244 \end{equation} 245 and simlar 246 \begin{equation} 247 v\hackscore{2} = v\hackscore{1} - dv\hackscore{2} 248 = v\hackscore{1} - A^{-1} B^{*}dp 249 = e\hackscore{1} + A^{-1} ( G - B^{*} p ) - A^{-1} B^{*} (p\hackscore{2}-p) 250 = e\hackscore{1} + A^{-1} ( G - B^{*} p\hackscore{2}) 251 \end{equation} 252 This shows that we can write the iterative process as a fixed point iteration to solve (assume all errors are zero) 253 \begin{equation} 254 v = \Phi(v,p) \mbox{ and } p = \Psi(u,p) 255 \end{equation} 256 where 257 \begin{equation} 258 \begin{array}{rcl} 259 \Psi(v,p) & = & (B A^{-1} B^{*})^{-1} B A^{-1} G \\ 260 \Phi(u,p) & = & A^{-1} ( G - B^{*} (B A^{-1} B^{*})^{-1} B A^{-1} G ) 261 \end{array} 262 \end{equation} 263 Notice that if $A$ is independent from $v$ and $p$ the operators $\Phi(v,p)$ and $\Psi(u,p)$ are constant 264 and threfore the iteration will terminate - providing no termination errors in sub-iterations - after one step. 265 We also can give a formula for the error 266 \begin{equation} 267 \begin{array}{rcl} 268 \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 \;. 270 \end{array}\label{STOKES ERRORS} 271 \end{equation} 272 Notice that $B\delta v = - e\hackscore{2}=-Bv\hackscore{2}$. 273 274 With this notation 275 \begin{equation} 276 \begin{array}{rcl} 277 v\hackscore{2} = \Phi(v,p) + \delta v \\ 278 p\hackscore{2} = \Psi(v,p) + \delta p \\ 279 \end{array} 280 \end{equation} 281 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. 283 Assuming that the iteration does converge with a convergence rate $\chi^{-}$ we have the estimate 284 \begin{equation} 285 \max(\|dv\|\hackscore{1}, \|Bv\hackscore{1}\|\hackscore{0}) 286 \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}) 288 \end{equation} 289 were the upper index $-$ referes to the increments in the last step. 290 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 293 $\delta p$ is controlled. 294 \begin{equation} 295 \|\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} 296 \end{equation} 297 and 298 \begin{equation} 299 \|(B A^{-1} B^{*}) \delta p\|\hackscore{0} \approx \| e\hackscore{2} \|\hackscore{0} 300 = \| B v\hackscore{2} \|\hackscore{0} \le M \tau\hackscore{2} \| B v\hackscore{1} \|\hackscore{0} 301 \end{equation} 302 which leads to 303 \begin{equation} 304 \max(\|dv\|\hackscore{1}, \|Bv\hackscore{1}\|\hackscore{0}) 305 \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} 306 \end{equation} 307 where the upper index $-$ refers to the previous iteration step. 308 If we allow require $max{(M \tau\hackscore{2}, \frac{K \tau\hackscore{1}}{1-K \tau\hackscore{1}})}\le \gamma$ 309 with a given $0<\gamma<1$ we can set 310 \begin{equation} \label{STOKES SET TOL} 311 \tau\hackscore{2} = \frac{\gamma}{M} \mbox{ and } \tau\hackscore{1} = \frac{1}{K} \frac{\gamma}{1+\gamma} 312 \end{equation} 313 Problem is that we do not know $M$ and $K$ but can use the convergence rate 314 \begin{equation} 315 \chi := \frac{\max(\|dv\|\hackscore{1}, \|Bv\hackscore{1}\|\hackscore{0})}{\max(\|dv^{-}\|\hackscore{1}, \|Bv\hackscore{1}^{-}\|\hackscore{0})} 316 \end{equation} 317 to construct estimates of $M$ and $K$. We look at the two cases where our prediction and choice of 318 the tolerances was good or where we went wrong: 319 \begin{equation} 320 \chi \le \frac{\chi^{-}}{1-\gamma} \mbox{ or } 321 \chi = \frac{\chi^{-}}{1-\gamma\hackscore{0}} >\frac{\chi^{-}}{1-\gamma}. 322 \end{equation} 323 which translates to 324 \begin{equation} 325 \gamma\hackscore{0} \le \gamma \mbox{ or } \gamma\hackscore{0}=1-\frac{\chi^{-}}{ \chi}>\gamma>0 326 \end{equation} 327 In the second case use \ref{STOKES SET TOL} for $\gamma=\gamma\hackscore{0}$ to get new estimates $M^{+}$ and $K^{+}$ 328 for $M$ and $K$: 329 \begin{equation} \label{TOKES CONST UPDATE} 330 M^{+} = 331 \frac{\gamma\hackscore{0}}{\gamma} M 332 \mbox{ and } K^{+} = 333 \frac{\gamma\hackscore{0}}{\gamma} 334 \frac{1+\gamma}{1+\gamma\hackscore{0}} 335 K 336 \mbox{ with } \gamma\hackscore{0}=\max(1-\frac{\chi^{-}}{ \chi},\gamma) 337 \end{equation} 338 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 339 meassure by 340 \begin{equation} 341 \chi < \chi\hackscore{max} 342 \end{equation} 343 where $\chi\hackscore{max}$ is given value with $0<\chi\hackscore{max}<1$. We then consider the following 344 criteria 345 \begin{itemize} 346 \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 347 $\gamma^{+} \ge \frac{\gamma}{\omega}$. 348 \item We need to make sure that the estimated convergence rate based on $\gamma^{+}$ will still maintain convergence 349 \begin{equation} 350 \frac{\chi}{1-\gamma^{+}} \le \chi\hackscore{max} 351 \end{equation} 352 which translates into 353 \begin{equation} 354 \gamma^{+} \le 1 - \frac{\chi}{\chi\hackscore{max}} 355 \end{equation} 356 Notice that the right hand side is positive. 357 \item We do not want to iterate more then necessary. So we would like reduce the error in the next just below the 358 desired tolerance: 359 \begin{equation} 360 \frac{\chi}{1-\gamma^{+}}\max(\|dv\|\hackscore{1}, \|Bv\hackscore{1}\|\hackscore{0}) \ge f \cdot \tau \|v\hackscore{2}\|\hackscore{1} 361 \end{equation} 362 where the left hand side provides an estimate for the error to prediced in the next step. $f$ is a 363 safty factor. This leads to 364 \begin{equation} 365 \gamma^{+} \le 366 1 - 367 \chi \frac{\max(\|dv\|\hackscore{1}, \|Bv\hackscore{1}\|\hackscore{0})} { f \cdot \tau \|v\hackscore{2}\|\hackscore{1}} 368 \end{equation} 369 \end{itemize} 370 Putting all these criteria together we get 371 \begin{equation} 372 \gamma^{+}=\min(\max( 373 \frac{\gamma}{\omega}, 374 1 - 375 \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}}) 376 \label{STOKES SET GAMMA PLUS} 377 \end{equation} 378 In case we cannot see convergence $\gamma$ is reduced by the factor $\omega$ 379 \begin{equation} 380 \gamma^{+}=\max(\omega \cdot \gamma, \gamma\hackscore{min})\label{STOKES SET GAMMA PLUS2} 381 \end{equation} 382 where $\gamma\hackscore{min}$ is introduced to make sure that $\gamma^{+}$ stays away from zero. 383 384 385 Appling~\ref{STOKES SET TOL} for $\gamma^{+}$ and~\ref{TOKES CONST UPDATE} we get the update formula 386 \begin{equation} \label{STOKES CONST UPDATE} 387 \tau\hackscore{2}^{+} = 388 \frac{\gamma^{+}}{\gamma\hackscore{0}} \tau\hackscore{2} 389 \mbox{ and } \tau\hackscore{1}^{+} = 390 \frac{\gamma^{+}}{\gamma\hackscore{0}} 391 \frac{1+\gamma\hackscore{0}}{1+\gamma^{+}} 392 \tau\hackscore{1} 393 \end{equation} 394 to get the new tolerances $\tau\hackscore{1}^{+}$ and $\tau\hackscore{2}^{+}$ to be used in the next iteration step. 395 Notice that the coefficients $M$ and $K$ have been eliminated. The calculation of $\gamma\hackscore{0}$ requires 396 at least two iteration steps (or a previous time step). 397 398 Typical initial values are 399 \begin{equation} \label{STOKES VALUES} 400 \tau\hackscore{1}=\tau\hackscore{1}=\frac{1}{100} ; \; \omega=\frac{1}{2} \; \gamma=\frac{1}{10} ;\gamma\hackscore{min}=10^{-6} 401 \end{equation} 402 where after the first step we set $\gamma\hackscore{0}=\gamma$ as no $\chi^{-}$ is available. 403 404 405 \subsection{Functions} 406 407 \begin{classdesc}{StokesProblemCartesian}{domain} 408 opens the Stokes problem\index{Stokes problem} on the \Domain domain. The domain 409 needs to support LBB compliant elements for the Stokes problem, see~\cite{LBB} for detials~\index{LBB condition}. 410 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} 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{ 421 restoration_factor=0}}}}}} 422 assigns values to the model parameters. In any call all values must be set. 423 \var{f} defines the external force $f$, \var{eta} the viscosity $\eta$, 424 \var{surface_stress} the surface stress $s$ and \var{stress} the initial stress $\sigma$. 425 The locations and compontents where the velocity is fixed are set by 426 the values of \var{fixed_u_mask}. \var{restoration_factor} defines the restoring force factor $\alpha$. 427 The method will try to cast the given values to appropriate 428 \Data class objects. 429 \end{methoddesc} 430 431 \begin{methoddesc}[StokesProblemCartesian]{solve}{v,p 432 \optional{, max_iter=100 \optional{, verbose=False \optional{, usePCG=True }}}} 433 solves the problem and return approximations for velocity and pressure. 434 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 438 by \var{fixed_u_mask} remain unchanged. 439 If \var{usePCG} is set to \True 440 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 441 the PCG scheme is more efficient. 442 \var{max_iter} defines the maximum number of iteration steps. 443 444 If \var{verbose} is set to \True informations on the progress of of the solver are printed. 445 \end{methoddesc} 446 447 448 \begin{methoddesc}[StokesProblemCartesian]{setTolerance}{\optional{tolerance=1.e-4}} 449 sets the tolerance in an appropriate norm relative to the right hand side. The tolerance must be non-negative and less than 1. 450 \end{methoddesc} 451 \begin{methoddesc}[StokesProblemCartesian]{getTolerance}{} 452 returns the current relative tolerance. 453 \end{methoddesc} 454 \begin{methoddesc}[StokesProblemCartesian]{setAbsoluteTolerance}{\optional{tolerance=0.}} 455 sets the absolute tolerance for the error in the relevant norm. The tolerance must be non-negative. Typically the 456 absolute talerance is set to 0. 457 \end{methoddesc} 458 459 \begin{methoddesc}[StokesProblemCartesian]{getAbsoluteTolerance}{} 460 sreturns the current absolute tolerance. 461 \end{methoddesc} 462 463 \begin{methoddesc}[StokesProblemCartesian]{getSolverOptionsVelocity}{} 464 returns the solver options used solve the equations~(\ref{V CALC}) for velocity. 465 \end{methoddesc} 466 467 \begin{methoddesc}[StokesProblemCartesian]{getSolverOptionsPressure}{} 468 returns the solver options used solve the equation~(\ref{P PREC}) for pressure. 469 \end{methoddesc} 470 471 \begin{methoddesc}[StokesProblemCartesian]{getSolverOptionsDiv}{} 472 set the solver options for solving the equation to project the divergence of the velocity onto the function space of pressure. 473 \end{methoddesc} 474 475 476 \subsection{Example: Lit Driven Cavity} 477 The following script \file{lit\hackscore driven\hackscore cavity.py} 478 \index{scripts!\file{helmholtz.py}} which is available in the \ExampleDirectory 479 illustrates the usage of the \class{StokesProblemCartesian} class to solve 480 the lit driven cavity problem: 481 \begin{python} 482 from esys.escript import * 483 from esys.finley import Rectangle 484 from esys.escript.models import StokesProblemCartesian 485 NE=25 486 dom = Rectangle(NE,NE,order=2) 487 x = dom.getX() 488 sc=StokesProblemCartesian(dom) 489 mask= (whereZero(x)*[1.,0]+whereZero(x-1))*[1.,0] + \ 490 (whereZero(x)*[0.,1.]+whereZero(x-1))*[1.,1] 491 sc.initialize(eta=.1, fixed_u_mask= mask) 492 v=Vector(0.,Solution(dom)) 493 v+=whereZero(x-1.) 494 p=Scalar(0.,ReducedSolution(dom)) 495 v,p=sc.solve(v,p, verbose=True) 496 saveVTK("u.xml",velocity=v,pressure=p) 497 \end{python}