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

Revision 3295 - (show annotations)
Fri Oct 22 01:56:02 2010 UTC (10 years, 3 months ago) by jfenwick
File MIME type: application/x-tex
File size: 19088 byte(s)

 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_{i,j}+ v_{j,i})\right)_{,j}+p_{,i}=f_{i}-\sigma_{ij,j} 8 \end{equation} 9 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. 10 11 We assume an incompressible media: 12 \begin{equation}\label{Stokes 2} 13 -v_{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_{i,j}+ v_{j,i})\right)n_{j}-n_{i}p=s_{i} - \alpha \cdot n_{i} n_{j} v_{j}+\sigma_{ij} n_{j} 18 \end{equation} 19 which can be overwritten by constraints of the form 20 \begin{equation}\label{Stokes Boundary0} 21 v_{i}(x)=v^D_{i}(x) 22 \end{equation} 23 at some locations $x$ at the boundary of the domain. $s_{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{STOKES} 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_{0}$ and $p_{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{STOKES} 52 \begin{equation}\label{STOKES ITER STEP 1} 53 A dv_{1} = G - A v_{0} - B^{*} p_{0} 54 \end{equation} 55 We then insert the new approximation $v_{1}=v_{0}+dv_{1}$ to calculate a correction $dp_{2}$ 56 for the pressure and an additional correction $dv_{2}$ for the velocity by solving 57 \begin{equation} 58 \begin{array}{rcl} 59 B A^{-1} B^{*} dp_{2} & = & Bv_{1} \\ 60 A dv_{2} & = & B^{*} dp_{2} 61 \end{array} 62 \label{STOKES ITER STEP 2} 63 \end{equation} 64 The new velocity and pressure are then given by $v_{2}=v_{1}-dv_{2}$ and 65 $p_{2}=p_{0}+dp_{2}$ which will fulfill the block system~\ref{STOKES}. 66 This solution strategy is called the Uzawa scheme \index{Uzawa scheme}. 67 68 There is a problem with this scheme: In practice we will use an iterative scheme 69 to solve any problem for operator $A$. So we will be unable to calculate the operator 70 $B A^{-1} B^{*}$ required for $dp_{2}$ explicitly. In fact, we need to use another 71 iterative scheme to solve the first equation in~\ref{STOKES ITER STEP 2} where in each iteration step 72 an iterative solver for $A$ is applied. Another issue is the fact that the 73 viscosity $\eta$ may depend on velocity or pressure and so we need to iterate over the 74 three equations~\ref{STOKES ITER STEP 1} and~\ref{STOKES ITER STEP 2}. 75 76 In the following we will use the two norms 77 \begin{equation} 78 \|v\|_{1}^2 = \int_{\Omega} v_{j,k}v_{j,k} \; dx 79 \mbox{ and } 80 \|p\|_{0}^2= \int_{\Omega} p^2 \; dx. 81 \label{STOKES STOP} 82 \end{equation} 83 for velocity $v$ and pressure $p$. The iteration is terminated if the stopping criteria 84 \begin{equation} \label{STOKES STOPPING CRITERIA} 85 \max(\|Bv_{1}\|_{0},\|v_{2}-v_{0}\|_{1}) \le \tau \cdot \|v_{2}\|_{1} 86 \end{equation} 87 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 88 $\|Bv_{1}\|_{0}$ equals the 89 norm of $B A^{-1} B^{*} dp_{2}$ and consequently provides a norm for the pressure correction. 90 91 We want to optimize the tolerance choice for solving~\ref{STOKES ITER STEP 1} 92 and~\ref{STOKES ITER STEP 2}. To do this we write the iteration scheme as a fixed point problem. Here 93 we consider the errors produced by the iterative solvers being used. 94 From Equation~\ref{STOKES ITER STEP 1} we have 95 \begin{equation} \label{STOKES total V1} 96 v_{1} = e_{1} + v_{0} + A^{-1} ( G - Av_{0} - B^{*} p_{0} ) 97 \end{equation} 98 where $e_{1}$ is the error when solving~\ref{STOKES ITER STEP 1}. 99 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 100 \begin{equation} 101 \| A e_{1} \|_{s} = \| G - A v_{1} - B^{*} p_{0} \|_{s} \le \tau_{1} \| G - A v_{0} - B^{*} p_{0} \|_{s} 102 \end{equation} 103 where $\tau_{1}$ is the tolerance which we need to choose. This translates into the condition 104 \begin{equation} 105 \| e_{1} \|_{1} \le K \tau_{1} \| dv_{1} - e_{1} \|_{1} 106 \end{equation} 107 The constant $K$ represents some uncertainty combining a variety of unknown factors such as the 108 norm being used and the condition number of the stiffness matrix. 109 From the first equation of~\ref{STOKES ITER STEP 2} we have 110 \begin{equation}\label{STOKES total P2} 111 p_{2} = p_{0} + (B A^{-1} B^{*})^{-1} (e_{2} + Bv_{1} ) 112 \end{equation} 113 where $e_{2}$ represents the error when solving~\ref{STOKES ITER STEP 2}. 114 We use an iterative preconditioned conjugate gradient method (PCG) \index{linear solver!PCG}\index{PCG} with iteration 115 operator $B A^{-1} B^{*}$ using the $\|.\|_{0}$ norm. As suitable preconditioner \index{preconditioner} for the iteration 116 operator we use $\frac{1}{\eta}$ \cite{ELMAN}, ie 117 the evaluation of the preconditioner $P$ for a given pressure increment $q$ is the solution of 118 \begin{equation} \label{STOKES P PREC} 119 \frac{1}{\eta} (Pq) = q \; . 120 \end{equation} 121 Note that in each evaluation of the iteration operator $q=B A^{-1} B^{*} s$ one needs to solve 122 the problem 123 \begin{equation} \label{STOKES P OPERATOR} 124 A w = B^{*} s 125 \end{equation} 126 with sufficient accuracy to return $q=Bw$. We assume that the desired tolerance is 127 sufficiently small, for instance one can take $\tau_{2}^2$ 128 where $\tau_{2}$ is the tolerance for~\ref{STOKES ITER STEP 2}. 129 130 In an implementation we use the fact that the residual $r$ is given as 131 \begin{equation} \label{STOKES RES } 132 r= B (v_{1} - A^{-1} B^{*} dp) = B (v_{1} - A^{-1} B^{*} dp) = B (v_{1}-dv_{2}) = B v_{2} 133 \end{equation} 134 In particular we have $e_{2} = B v_{2}$ 135 So the residual $r$ is represented by the updated velocity $v_{2}=v_{1}-dv_{2}$. In practice, if 136 one uses the velocity to represent the residual $r$ there is no need 137 to recover $dv_{2}$ in~\ref{STOKES ITER STEP 2} after $dp_{2}$ has been calculated. 138 In PCG the iteration is terminated if 139 \begin{equation} \label{STOKES P OPERATOR ERROR} 140 \| P^{\frac{1}{2}}B v_{2} \|_{0} \le \tau_{2} \| P^{\frac{1}{2}}B v_{1} \|_{0} 141 \end{equation} 142 where $\tau_{2}$ is the given tolerance. This translates into 143 \begin{equation} \label{STOKES P OPERATOR ERROR 2} 144 \|e_{2}\|_{0} = \| B v_{2} \|_{0} \le M \tau_{2} \| B v_{1} \|_{0} 145 \end{equation} 146 where $M$ is taking care of the fact that $P^{\frac{1}{2}}$ is dropped. 147 148 As we assume that there is no significant error from solving with the operator $A$ we have 149 \begin{equation} \label{STOKES total V2} 150 v_{2} = v_{1} - dv_{2} 151 = v_{1} - A^{-1} B^{*}dp 152 \end{equation} 153 Combining the equations~\ref{STOKES total V1},~\ref{STOKES total P2} and~\ref{STOKES total V2} and 154 setting the errors to zero we can write the solution process as a fix point problem 155 \begin{equation} 156 v = \Phi(v,p) \mbox{ and } p = \Psi(u,p) 157 \end{equation} 158 with suitable functions $\Phi(v,p)$ and $\Psi(v,p)$ representing the iteration operator without 159 errors. In fact for a linear problem, $\Phi$ and $\Psi$ are constant. With this notation we can write 160 the update step in the form $p_{2}= \delta p + \Psi(v_{0},p_{0})$ and 161 $v_{2}= \delta v + \Phi(v_{0},p_{0})$ where 162 the total error $\delta p$ and $\delta v$ are given as 163 \begin{equation} 164 \begin{array}{rcl} 165 \delta p & = & (B A^{-1} B^{*})^{-1} ( e_{2} + B e_{1} ) \\ 166 \delta v & = & e_{1} - A^{-1} B^{*}\delta p \;. 167 \end{array}\label{STOKES ERRORS} 168 \end{equation} 169 Notice that $B\delta v = - e_{2}=-Bv_{2}$. Our task is now to choose the tolerances 170 $\tau_{1}$ and $\tau_{2}$ such that the global errors $\delta p$ and $\delta v$ 171 do not stop the convergence of the iteration process. 172 173 To measure convergence we use 174 \begin{equation} 175 \epsilon = \max(\|v_{2}-v\|_{1}, \|B A^{-1} B^{*} (p_{2}-p)\|_{0}) 176 \end{equation} 177 In practice using the fact that $B A^{-1} B^{*} (p_{2}-p_{0}) = B v_{1}$ 178 and assuming that $v_{2}$ gives a better approximation to the true $v$ than 179 $v_{0}$ we will use the estimate 180 \begin{equation} 181 \epsilon = \max(\|v_{2}-v_{0}\|_{1}, \|B v_{1}\|_{0}) 182 \end{equation} 183 to estimate the progress of the iteration step after the step is completed. 184 Note that the estimate of $\epsilon$ 185 used in the stopping criteria~\ref{STOKES STOPPING CRITERIA}. If $\chi^{-}$ is the convergence rate assuming 186 exact calculations, i.e. $e_{1}=0$ and $e_{2}=0$, we are expecting 187 to maintain $\epsilon \le \chi^{-} \cdot \epsilon^{-}$. For the 188 pressure increment we get: 189 \begin{equation} \label{STOKES EST 1} 190 \begin{array}{rcl} 191 \|B A^{-1} B^{*} (p_{2}-p)\|_{0} 192 & \le & \|B A^{-1} B^{*} (p_{2}-\delta p-p)\|_{0} + 193 \|B A^{-1} B^{*} \delta p\|_{0} \\ 194 & = & \chi^{-} \cdot \epsilon^{-} + \|e_{2} + B e_{1}\|_{0} \\ 195 & \approx & \chi^{-} \cdot \epsilon^{-} + \|e_{2}\|_{0} \\ 196 & \le & \chi^{-} \cdot \epsilon^{-} + M \tau_{2} \|B v_{1}\|_{0} \\ 197 \end{array} 198 \end{equation} 199 So we choose the value for $\tau_{2}$ from 200 \begin{equation} \label{STOKES TOL2} 201 M \tau_{2} \|B v_{1}\|_{0} \le (\chi^{-})^2 \epsilon^{-} 202 \end{equation} 203 in order to make the perturbation for the termination of the pressure iteration a second order effect. We use a 204 similar argument for the velocity: 205 \begin{equation}\label{STOKES EST 2} 206 \begin{array}{rcl} 207 \|v_{2}-v\|_{1} & \le & \|v_{2}-\delta v-v\|_{1} + \| \delta v\|_{1} \\ 208 & \le & \chi^{-} \cdot \epsilon^{-} + \| e_{1} - A^{-1} B^{*}\delta p \|_{1} \\ 209 & \approx & \chi^{-} \cdot \epsilon^{-} + \| e_{1} \|_{1} \\ 210 & \le & \chi^{-} \cdot \epsilon^{-} + K \tau_{1} \| dv_{1} - e_{1} \|_{1} 211 \\ 212 & \le & ( 1 + K \tau_{1}) \chi^{-} \cdot \epsilon^{-} 213 \end{array} 214 \end{equation} 215 So we choose the value for $\tau_{1}$ from 216 \begin{equation} \label{STOKES TOL1} 217 K \tau_{1} \le \chi^{-} 218 \end{equation} 219 Assuming we have estimates for $M$ and $K$ (if none is available, we use the value $1$.) 220 we can use~\ref{STOKES TOL1} and~\ref{STOKES TOL2} to get appropriate values for the tolerances. After 221 the step has been completed we can calculate a new convergence rate $\chi =\frac{\epsilon}{\epsilon^{-}}$. 222 For partial reasons we restrict $\chi$ to be less or equal a given maximum value $\chi_{max}\le 1$. 223 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 224 \begin{equation}\label{STOKES EST 3} 225 \chi \le ( 1 + \max(M \frac{\tau_{2} \|B v_{1}\|_{0}}{\chi^{-} \epsilon^{-}},K \tau_{1} ) ) \cdot \chi^{-} 226 \end{equation} 227 If we assume that this inequality would be an equation if we would have chosen the right values 228 $M^{+}$ and $K^{+}$ then we get 229 \begin{equation}\label{STOKES EST 3b} 230 \chi = ( 1 + \max(M^{+} \frac{\chi^{-}}{M},K^{+} \frac{\chi^{-}}{K}) ) \cdot \chi^{-} 231 \end{equation} 232 From this equation we set for 233 than our choice for $K$ was not good enough. In this case we can calculate a new value 234 \begin{equation} 235 K^{+} = \frac{\chi-\chi^{-}}{(\chi^{-})^2} K 236 \end{equation} 237 In practice we will use 238 \begin{equation} 239 K^{+} = \max(\frac{\chi-\chi^{-}}{(\chi^{-})^2} K,\frac{1}{2}K,1) 240 \end{equation} 241 where the second term is used to reduce a potential overestimate of $K$. 242 The same identity is used for to update $M$. The updated $M^{+}$ and $K^{+}$ 243 are then use in the next iteration step to control the tolerances. 244 245 In some cases one can observe that there is a significant change 246 in the velocity but the new velocity $v_{1}$ has still a 247 small divergence, i.e. we have 248 $\|Bv_{1}\|_{0} \ll \|v_{1}-v_{0}\|_{1}$. 249 In this case we will get a small pressure increment and consequently only very small changes to 250 the velocity as a result of the second update step which therefore can be skipped and 251 we can directly repeat the first update step until the increment in velocity becomes 252 significant relative to its divergence. In practice we will ignore the second half of the iteration step 253 as long as 254 \begin{equation}\label{STOKES LARGE BV1} 255 \|Bv_{1}\|_{0} \le \theta \cdot \|v_{1}-v_{0}\| 256 \end{equation} 257 where $0<\theta<1$ is a given factor. In this case we will also check the stopping criteria 258 with $v_{1}\rightarrow v_{2}$ but we will not correct $M$ in this case. 259 260 Starting from an initial guess $v_{0}$ and $p_{0}$ for velocity and pressure 261 the solution procedure is implemented as follows. 262 \begin{enumerate} 263 \item calculate viscosity $\eta(v_{0},p)_{0}$ and assemble operator $A$ from $\eta$. 264 \item calculate the tolerance $\tau_{1}$ from~\ref{STOKES TOL1}. 265 \item Solve~\ref{STOKES ITER STEP 1} for $dv_{1}$ with tolerance $\tau_{1}$. 266 \item update $v_{1}= v_{0}+ dv_{1}$ 267 \item if $Bv_{1}$ is large (see~\ref{STOKES LARGE BV1}): 268 \begin{enumerate} 269 \item calculate the tolerance $\tau_{2}$ from~\ref{STOKES TOL2}. 270 \item Solve~\ref{STOKES ITER STEP 2} for $dp_{2}$ and $v_{2}$ with tolerance $\tau_{2}$. 271 \item update $p_{2}\leftarrow p_{0}+ dp_{2}$ 272 \end{enumerate} 273 \item else: 274 \begin{enumerate} 275 \item update $p_{2}\leftarrow p$ and $v_{2}\leftarrow v_{1}$ 276 \end{enumerate} 277 \item calculate convergence measure $\epsilon$ and convergence rate $\chi$ 278 \item if stopping criteria~\ref{STOKES STOPPING CRITERIA} holds: 279 \begin{enumerate} 280 \item return $v_{2}$ and $p_{2}$ 281 \end{enumerate} 282 \item else: 283 \begin{enumerate} 284 285 \item update $M$ and $K$ 286 \item goto step 1 with $v_{0}\leftarrow v_{2}$ and $p_{0}\leftarrow p_{2}$. 287 \end{enumerate} 288 \end{enumerate} 289 290 \subsection{Functions} 291 292 \begin{classdesc}{StokesProblemCartesian}{domain} 293 opens the Stokes problem\index{Stokes problem} on the \Domain domain. The domain 294 needs to support LBB compliant elements for the Stokes problem, see~\cite{LBB} for details~\index{LBB condition}. 295 For instance one can use second order polynomials for velocity and 296 first order polynomials for the pressure on the same element. Alternatively, one can use 297 macro elements\index{macro elements} using linear polynomial for both pressure and velocity with a subdivided 298 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 299 \begin{python} 300 velocity=Vector(0.0, Solution(mesh)) 301 pressure=Scalar(0.0, ReducedSolution(mesh)) 302 \end{python} 303 \end{classdesc} 304 305 \begin{methoddesc}[StokesProblemCartesian]{initialize}{\optional{f=Data(), \optional{fixed_u_mask=Data(), \optional{eta=1, \optional{surface_stress=Data(), \optional{stress=Data()}, \optional{ 306 restoration_factor=0}}}}}} 307 assigns values to the model parameters. In any call all values must be set. 308 \var{f} defines the external force $f$, \var{eta} the viscosity $\eta$, 309 \var{surface_stress} the surface stress $s$ and \var{stress} the initial stress $\sigma$. 310 The locations and components where the velocity is fixed are set by 311 the values of \var{fixed_u_mask}. \var{restoration_factor} defines the restoring force factor $\alpha$. 312 The method will try to cast the given values to appropriate 313 \Data class objects. 314 \end{methoddesc} 315 316 \begin{methoddesc}[StokesProblemCartesian]{solve}{v,p 317 \optional{, max_iter=100 \optional{, verbose=False \optional{, usePCG=True }}}} 318 solves the problem and return approximations for velocity and pressure. 319 The arguments \var{v} and \var{p} define initial guess. 320 \var{v} must have function space \var{Solution(domain)} and 321 \var{p} must have function space \var{ReducedSolution(domain)}. 322 The values of \var{v} marked 323 by \var{fixed_u_mask} remain unchanged. 324 If \var{usePCG} is set to \True 325 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 326 the PCG scheme is more efficient. 327 \var{max_iter} defines the maximum number of iteration steps. 328 329 If \var{verbose} is set to \True informations on the progress of of the solver are printed. 330 \end{methoddesc} 331 332 333 \begin{methoddesc}[StokesProblemCartesian]{setTolerance}{\optional{tolerance=1.e-4}} 334 sets the tolerance in an appropriate norm relative to the right hand side. The tolerance must be non-negative and less than 1. 335 \end{methoddesc} 336 \begin{methoddesc}[StokesProblemCartesian]{getTolerance}{} 337 returns the current relative tolerance. 338 \end{methoddesc} 339 \begin{methoddesc}[StokesProblemCartesian]{setAbsoluteTolerance}{\optional{tolerance=0.}} 340 sets the absolute tolerance for the error in the relevant norm. The tolerance must be non-negative. Typically the 341 absolute tolerance is set to 0. 342 \end{methoddesc} 343 344 \begin{methoddesc}[StokesProblemCartesian]{getAbsoluteTolerance}{} 345 returns the current absolute tolerance. 346 \end{methoddesc} 347 348 \begin{methoddesc}[StokesProblemCartesian]{getSolverOptionsVelocity}{} 349 returns the solver options used solve the equations~(\ref{STOKES ITER STEP 1}) and~(\ref{STOKES P OPERATOR}) for velocity. 350 \end{methoddesc} 351 352 \begin{methoddesc}[StokesProblemCartesian]{getSolverOptionsPressure}{} 353 returns the solver options used solve the preconditioner equation~(\ref{STOKES P PREC}) for pressure. 354 \end{methoddesc} 355 356 \begin{methoddesc}[StokesProblemCartesian]{getSolverOptionsDiv}{} 357 set the solver options for solving the equation to project the divergence of the velocity onto the function space of pressure. 358 \end{methoddesc} 359 360 361 \subsection{Example: Lit Driven Cavity} 362 The following script \file{lit_driven_cavity.py} 363 \index{scripts!\file{helmholtz.py}} which is available in the \ExampleDirectory 364 illustrates the usage of the \class{StokesProblemCartesian} class to solve 365 the lit driven cavity problem: 366 \begin{python} 367 from esys.escript import * 368 from esys.finley import Rectangle 369 from esys.escript.models import StokesProblemCartesian 370 NE=25 371 dom = Rectangle(NE,NE,order=2) 372 x = dom.getX() 373 sc=StokesProblemCartesian(dom) 374 mask= (whereZero(x[0])*[1.,0]+whereZero(x[0]-1))*[1.,0] + \ 375 (whereZero(x[1])*[0.,1.]+whereZero(x[1]-1))*[1.,1] 376 sc.initialize(eta=.1, fixed_u_mask= mask) 377 v=Vector(0.,Solution(dom)) 378 v[0]+=whereZero(x[1]-1.) 379 p=Scalar(0.,ReducedSolution(dom)) 380 v,p=sc.solve(v,p, verbose=True) 381 saveVTK("u.xml",velocity=v,pressure=p) 382 \end{python}