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

Revision 6651 - (hide annotations)
Wed Feb 7 02:12:08 2018 UTC (2 years, 11 months ago) by jfenwick
File MIME type: application/x-tex
File size: 21102 byte(s)
Make everyone sad by touching all the files


 1 caltinay 3325 2 jfenwick 3989 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 3 jfenwick 6651 % Copyright (c) 2003-2018 by The University of Queensland 4 jfenwick 3989 5 caltinay 3325 % 6 % Primary Business: Queensland, Australia 7 jfenwick 6112 % Licensed under the Apache License, version 2.0 8 9 caltinay 3325 % 10 jfenwick 3989 % Development until 2012 by Earth Systems Science Computational Center (ESSCC) 11 jfenwick 4657 % Development 2012-2013 by School of Earth Sciences 12 % Development from 2014 by Centre for Geoscience Computing (GeoComp) 13 jfenwick 3989 % 14 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 15 caltinay 3325 16 gross 2719 \section{The Stokes Problem} 17 \label{STOKES PROBLEM} 18 caltinay 3325 In this section we discuss how to solve the Stokes problem. 19 We want to calculate the velocity\index{velocity} field $v$ and pressure $p$ 20 of an incompressible fluid\index{incompressible fluid}. 21 They are given as the solution of the Stokes problem\index{Stokes problem} 22 gross 2719 \begin{equation}\label{Stokes 1} 23 jfenwick 3295 -\left(\eta(v_{i,j}+ v_{j,i})\right)_{,j}+p_{,i}=f_{i}-\sigma_{ij,j} 24 gross 2719 \end{equation} 25 caltinay 3325 where $f_{i}$ defines an internal force\index{force, internal} and 26 $\sigma_{ij}$ is an initial stress\index{stress, initial}. 27 The viscosity $\eta$ may weakly depend on pressure and velocity. 28 If relevant we will use the notation $\eta(v,p)$ to express this dependency. 29 gross 2719 30 caltinay 3325 We assume an incompressible medium: 31 gross 2719 \begin{equation}\label{Stokes 2} 32 jfenwick 3295 -v_{i,i}=0 33 gross 2719 \end{equation} 34 Natural boundary conditions are taken in the form 35 \begin{equation}\label{Stokes Boundary} 36 jfenwick 3295 \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} 37 gross 2719 \end{equation} 38 which can be overwritten by constraints of the form 39 \begin{equation}\label{Stokes Boundary0} 40 jfenwick 3295 v_{i}(x)=v^D_{i}(x) 41 gross 2719 \end{equation} 42 caltinay 3325 at some locations $x$ at the boundary of the domain. 43 $s_{i}$ defines a normal stress and $\alpha\ge 0$ the spring constant for 44 restoring normal force. 45 gross 2719 The index $i$ may depend on the location $x$ on the boundary. 46 $v^D$ is a given function on the domain. 47 48 \subsection{Solution Method \label{STOKES SOLVE}} 49 caltinay 3325 If we assume that $\eta$ is independent from the velocity and pressure, 50 equations~\ref{Stokes 1} and~\ref{Stokes 2} can be written in the block form 51 gross 2719 \begin{equation} 52 \left[ \begin{array}{cc} 53 A & B^{*} \\ 54 B & 0 \\ 55 \end{array} \right] 56 \left[ \begin{array}{c} 57 v \\ 58 p \\ 59 \end{array} \right] 60 =\left[ \begin{array}{c} 61 G \\ 62 0 \\ 63 \end{array} \right] 64 gross 2843 \label{STOKES} 65 gross 2719 \end{equation} 66 caltinay 3325 where $A$ is a coercive, self-adjoint linear operator in a suitable Hilbert 67 space, $B$ is the $(-1) \cdot$ divergence operator and $B^{*}$ is the adjoint 68 operator (=gradient operator). 69 For more details on the mathematics see references \cite{AAMIRBERKYAN2008,MBENZI2005}. 70 gross 2793 71 caltinay 3325 If $v_{0}$ and $p_{0}$ are given initial guesses for velocity and pressure we 72 calculate a correction $dv$ for the velocity by solving the first equation of 73 \eqn{STOKES} 74 gross 2843 \begin{equation}\label{STOKES ITER STEP 1} 75 jfenwick 3295 A dv_{1} = G - A v_{0} - B^{*} p_{0} 76 gross 2793 \end{equation} 77 caltinay 3325 We then insert the new approximation $v_{1}=v_{0}+dv_{1}$ to calculate a 78 correction $dp_{2}$ for the pressure and an additional correction $dv_{2}$ for 79 the velocity by solving 80 gross 2793 \begin{equation} 81 \begin{array}{rcl} 82 jfenwick 3295 B A^{-1} B^{*} dp_{2} & = & Bv_{1} \\ 83 A dv_{2} & = & B^{*} dp_{2} 84 gross 2793 \end{array} 85 gross 2843 \label{STOKES ITER STEP 2} 86 gross 2793 \end{equation} 87 jfenwick 3295 The new velocity and pressure are then given by $v_{2}=v_{1}-dv_{2}$ and 88 $p_{2}=p_{0}+dp_{2}$ which will fulfill the block system~\ref{STOKES}. 89 caltinay 3325 This solution strategy is called the Uzawa scheme\index{Uzawa scheme}. 90 gross 2843 91 caltinay 3325 There is a problem with this scheme: in practice we will use an iterative 92 scheme to solve any problem for operator $A$. 93 So we will be unable to calculate the operator $B A^{-1} B^{*}$ required for 94 $dp_{2}$ explicitly. In fact, we need to use another iterative scheme to solve 95 the first equation in~\ref{STOKES ITER STEP 2} where in each iteration step 96 gross 2843 an iterative solver for $A$ is applied. Another issue is the fact that the 97 caltinay 3325 viscosity $\eta$ may depend on velocity or pressure and so we need to iterate 98 over the three equations~\ref{STOKES ITER STEP 1} and~\ref{STOKES ITER STEP 2}. 99 gross 2843 100 In the following we will use the two norms 101 gross 2793 \begin{equation} 102 jfenwick 3295 \|v\|_{1}^2 = \int_{\Omega} v_{j,k}v_{j,k} \; dx 103 gross 2793 \mbox{ and } 104 jfenwick 3295 \|p\|_{0}^2= \int_{\Omega} p^2 \; dx. 105 gross 2793 \label{STOKES STOP} 106 \end{equation} 107 caltinay 3325 for velocity $v$ and pressure $p$. 108 The iteration is terminated if the stopping criterion 109 gross 2843 \begin{equation} \label{STOKES STOPPING CRITERIA} 110 jfenwick 3295 \max(\|Bv_{1}\|_{0},\|v_{2}-v_{0}\|_{1}) \le \tau \cdot \|v_{2}\|_{1} 111 gross 2793 \end{equation} 112 caltinay 3325 for a given tolerance $0<\tau<1$ is met. 113 Notice that because of the first equation of~\ref{STOKES ITER STEP 2} we have 114 that $\|Bv_{1}\|_{0}$ equals the norm of $B A^{-1} B^{*} dp_{2}$ and 115 consequently provides a norm for the pressure correction. 116 gross 2793 117 gross 2843 We want to optimize the tolerance choice for solving~\ref{STOKES ITER STEP 1} 118 caltinay 3325 and~\ref{STOKES ITER STEP 2}. To do this we write the iteration scheme as a 119 fixed point problem. Here we consider the errors produced by the iterative 120 solvers being used. 121 From \eqn{STOKES ITER STEP 1} we have 122 gross 2843 \begin{equation} \label{STOKES total V1} 123 jfenwick 3295 v_{1} = e_{1} + v_{0} + A^{-1} ( G - Av_{0} - B^{*} p_{0} ) 124 gross 2793 \end{equation} 125 caltinay 3325 where $e_{1}$ is the error when solving~\ref{STOKES ITER STEP 1}. 126 We will use a sparse matrix solver so we have not full control on the norm 127 $\|.\|_{s}$ used in the stopping criterion for this equation. 128 In fact we will have a stopping criterion of the form 129 gross 2793 \begin{equation} 130 jfenwick 3295 \| A e_{1} \|_{s} = \| G - A v_{1} - B^{*} p_{0} \|_{s} \le \tau_{1} \| G - A v_{0} - B^{*} p_{0} \|_{s} 131 gross 2793 \end{equation} 132 caltinay 3325 where $\tau_{1}$ is the tolerance which we need to choose. 133 This translates into the condition 134 gross 2793 \begin{equation} 135 jfenwick 3295 \| e_{1} \|_{1} \le K \tau_{1} \| dv_{1} - e_{1} \|_{1} 136 gross 2793 \end{equation} 137 caltinay 3325 The constant $K$ represents some uncertainty combining a variety of unknown 138 factors such as the norm being used and the condition number of the stiffness matrix. 139 gross 2843 From the first equation of~\ref{STOKES ITER STEP 2} we have 140 \begin{equation}\label{STOKES total P2} 141 jfenwick 3295 p_{2} = p_{0} + (B A^{-1} B^{*})^{-1} (e_{2} + Bv_{1} ) 142 gross 2793 \end{equation} 143 jfenwick 3295 where $e_{2}$ represents the error when solving~\ref{STOKES ITER STEP 2}. 144 caltinay 3325 We use an iterative preconditioned conjugate gradient method 145 (PCG)\index{linear solver!PCG}\index{PCG} with iteration operator 146 $B A^{-1} B^{*}$ using the $\|.\|_{0}$ norm. 147 As suitable preconditioner\index{preconditioner} for the iteration operator we 148 use $\frac{1}{\eta}$ \cite{ELMAN}, i.e. the evaluation of the preconditioner 149 $P$ for a given pressure increment $q$ is the solution of 150 gross 2843 \begin{equation} \label{STOKES P PREC} 151 \frac{1}{\eta} (Pq) = q \; . 152 gross 2719 \end{equation} 153 caltinay 3325 Note that in each evaluation of the iteration operator $q=B A^{-1} B^{*} s$ 154 one needs to solve the problem 155 gross 2843 \begin{equation} \label{STOKES P OPERATOR} 156 gross 2719 A w = B^{*} s 157 \end{equation} 158 caltinay 3325 with sufficient accuracy to return $q=Bw$. We assume that the desired 159 tolerance is sufficiently small, for instance one can take $\tau_{2}^2$ where 160 $\tau_{2}$ is the tolerance for~\ref{STOKES ITER STEP 2}. 161 gross 2843 162 In an implementation we use the fact that the residual $r$ is given as 163 gross 2719 \begin{equation} \label{STOKES RES } 164 jfenwick 3295 r= B (v_{1} - A^{-1} B^{*} dp) = B (v_{1} - A^{-1} B^{*} dp) = B (v_{1}-dv_{2}) = B v_{2} 165 gross 2719 \end{equation} 166 caltinay 3325 In particular we have $e_{2} = B v_{2}$. 167 So the residual $r$ is represented by the updated velocity $v_{2}=v_{1}-dv_{2}$. 168 In practice, if one uses the velocity to represent the residual $r$ there is 169 no need to recover $dv_{2}$ in~\ref{STOKES ITER STEP 2} after $dp_{2}$ has 170 been calculated. 171 gross 2843 In PCG the iteration is terminated if 172 \begin{equation} \label{STOKES P OPERATOR ERROR} 173 jfenwick 3295 \| P^{\frac{1}{2}}B v_{2} \|_{0} \le \tau_{2} \| P^{\frac{1}{2}}B v_{1} \|_{0} 174 gross 2719 \end{equation} 175 jfenwick 3295 where $\tau_{2}$ is the given tolerance. This translates into 176 gross 2843 \begin{equation} \label{STOKES P OPERATOR ERROR 2} 177 jfenwick 3295 \|e_{2}\|_{0} = \| B v_{2} \|_{0} \le M \tau_{2} \| B v_{1} \|_{0} 178 gross 2719 \end{equation} 179 caltinay 3325 where $M$ is taking care of the fact that $P^{\frac{1}{2}}$ is dropped. 180 gross 2719 181 caltinay 3325 As we assume that there is no significant error from solving with the operator 182 $A$ we have 183 gross 2843 \begin{equation} \label{STOKES total V2} 184 jfenwick 3295 v_{2} = v_{1} - dv_{2} 185 = v_{1} - A^{-1} B^{*}dp 186 gross 2719 \end{equation} 187 caltinay 3325 Combining the equations~\ref{STOKES total V1},~\ref{STOKES total P2} and~\ref{STOKES total V2} 188 and setting the errors to zero we can write the solution process as a fix 189 point problem 190 gross 2719 \begin{equation} 191 v = \Phi(v,p) \mbox{ and } p = \Psi(u,p) 192 \end{equation} 193 caltinay 3325 with suitable functions $\Phi(v,p)$ and $\Psi(v,p)$ representing the 194 iteration operator without errors. In fact, for a linear problem, $\Phi$ and 195 $\Psi$ are constant. With this notation we can write the update step in the 196 form $p_{2}= \delta p + \Psi(v_{0},p_{0})$ and $v_{2}= \delta v + \Phi(v_{0},p_{0})$ 197 where the total error $\delta p$ and $\delta v$ are given as 198 gross 2719 \begin{equation} 199 \begin{array}{rcl} 200 jfenwick 3295 \delta p & = & (B A^{-1} B^{*})^{-1} ( e_{2} + B e_{1} ) \\ 201 \delta v & = & e_{1} - A^{-1} B^{*}\delta p \;. 202 gross 2719 \end{array}\label{STOKES ERRORS} 203 \end{equation} 204 caltinay 3325 Notice that $B\delta v = - e_{2}=-Bv_{2}$. 205 Our task is now to choose the tolerances $\tau_{1}$ and $\tau_{2}$ such that 206 the global errors $\delta p$ and $\delta v$ do not stop the convergence of the 207 iteration process. 208 gross 2793 209 gross 2843 To measure convergence we use 210 gross 2719 \begin{equation} 211 jfenwick 3295 \epsilon = \max(\|v_{2}-v\|_{1}, \|B A^{-1} B^{*} (p_{2}-p)\|_{0}) 212 gross 2719 \end{equation} 213 jfenwick 3295 In practice using the fact that $B A^{-1} B^{*} (p_{2}-p_{0}) = B v_{1}$ 214 and assuming that $v_{2}$ gives a better approximation to the true $v$ than 215 $v_{0}$ we will use the estimate 216 gross 2719 \begin{equation} 217 jfenwick 3295 \epsilon = \max(\|v_{2}-v_{0}\|_{1}, \|B v_{1}\|_{0}) 218 gross 2719 \end{equation} 219 caltinay 3325 to estimate the progress of the iteration step after the step is completed. 220 Note that the estimate of $\epsilon$ is used in the stopping 221 criterion~\ref{STOKES STOPPING CRITERIA}. 222 If $\chi^{-}$ is the convergence rate assuming exact calculations, i.e. 223 $e_{1}=0$ and $e_{2}=0$, we are expecting to maintain $\epsilon \le \chi^{-} \cdot \epsilon^{-}$. 224 For the pressure increment we get: 225 gross 2843 \begin{equation} \label{STOKES EST 1} 226 \begin{array}{rcl} 227 jfenwick 3295 \|B A^{-1} B^{*} (p_{2}-p)\|_{0} 228 & \le & \|B A^{-1} B^{*} (p_{2}-\delta p-p)\|_{0} + 229 \|B A^{-1} B^{*} \delta p\|_{0} \\ 230 & = & \chi^{-} \cdot \epsilon^{-} + \|e_{2} + B e_{1}\|_{0} \\ 231 & \approx & \chi^{-} \cdot \epsilon^{-} + \|e_{2}\|_{0} \\ 232 & \le & \chi^{-} \cdot \epsilon^{-} + M \tau_{2} \|B v_{1}\|_{0} \\ 233 gross 2843 \end{array} 234 gross 2719 \end{equation} 235 caltinay 3325 So we choose the value for $\tau_{2}$ from 236 gross 2843 \begin{equation} \label{STOKES TOL2} 237 jfenwick 3295 M \tau_{2} \|B v_{1}\|_{0} \le (\chi^{-})^2 \epsilon^{-} 238 gross 2719 \end{equation} 239 caltinay 3325 in order to make the perturbation for the termination of the pressure 240 iteration a second order effect. We use a similar argument for the velocity: 241 gross 2843 \begin{equation}\label{STOKES EST 2} 242 \begin{array}{rcl} 243 jfenwick 3295 \|v_{2}-v\|_{1} & \le & \|v_{2}-\delta v-v\|_{1} + \| \delta v\|_{1} \\ 244 & \le & \chi^{-} \cdot \epsilon^{-} + \| e_{1} - A^{-1} B^{*}\delta p \|_{1} \\ 245 & \approx & \chi^{-} \cdot \epsilon^{-} + \| e_{1} \|_{1} \\ 246 & \le & \chi^{-} \cdot \epsilon^{-} + K \tau_{1} \| dv_{1} - e_{1} \|_{1} 247 gross 2843 \\ 248 jfenwick 3295 & \le & ( 1 + K \tau_{1}) \chi^{-} \cdot \epsilon^{-} 249 gross 2843 \end{array} 250 gross 2719 \end{equation} 251 jfenwick 3295 So we choose the value for $\tau_{1}$ from 252 gross 2843 \begin{equation} \label{STOKES TOL1} 253 jfenwick 3295 K \tau_{1} \le \chi^{-} 254 gross 2719 \end{equation} 255 caltinay 3325 Assuming we have estimates for $M$ and $K$\footnote{if no estimates are 256 available, we use the value $1$} we can use~\ref{STOKES TOL1} and 257 \ref{STOKES TOL2} to get appropriate values for the tolerances. 258 After the step has been completed we can calculate a new convergence rate 259 $\chi =\frac{\epsilon}{\epsilon^{-}}$. 260 For partial reasons we restrict $\chi$ to be less or equal a given maximum 261 value $\chi_{max}\le 1$. 262 If we see $\chi \le \chi^{-} (1+\chi^{-})$ our choices for the tolerances were 263 suitable. Otherwise, we need to adjust the values for $K$ and $M$. 264 From the estimates~\ref{STOKES EST 1} and~\ref{STOKES EST 2} we establish 265 gross 2843 \begin{equation}\label{STOKES EST 3} 266 jfenwick 3295 \chi \le ( 1 + \max(M \frac{\tau_{2} \|B v_{1}\|_{0}}{\chi^{-} \epsilon^{-}},K \tau_{1} ) ) \cdot \chi^{-} 267 gross 2719 \end{equation} 268 caltinay 3325 If we assume that this inequality would be an equation if we would have chosen 269 the right values $M^{+}$ and $K^{+}$ then we get 270 gross 2843 \begin{equation}\label{STOKES EST 3b} 271 \chi = ( 1 + \max(M^{+} \frac{\chi^{-}}{M},K^{+} \frac{\chi^{-}}{K}) ) \cdot \chi^{-} 272 gross 2719 \end{equation} 273 caltinay 3325 From this equation we see if our choice for $K$ was not good enough. 274 In this case we can calculate a new value 275 gross 2843 \begin{equation} 276 K^{+} = \frac{\chi-\chi^{-}}{(\chi^{-})^2} K 277 gross 2719 \end{equation} 278 caltinay 3325 In practice we will use 279 gross 2843 \begin{equation} 280 K^{+} = \max(\frac{\chi-\chi^{-}}{(\chi^{-})^2} K,\frac{1}{2}K,1) 281 gross 2719 \end{equation} 282 caltinay 3325 where the second term is used to reduce a potential overestimate of $K$. 283 The same identity is used for to update $M$. The updated $M^{+}$ and $K^{+}$ 284 are then use in the next iteration step to control the tolerances. 285 gross 2719 286 caltinay 3325 In some cases one can observe that there is a significant change in the 287 velocity but the new velocity $v_{1}$ has still a small divergence, i.e. 288 we have $\|Bv_{1}\|_{0} \ll \|v_{1}-v_{0}\|_{1}$. 289 In this case we will get a small pressure increment and consequently only very 290 small changes to the velocity as a result of the second update step which 291 therefore can be skipped and we can directly repeat the first update step 292 until the increment in velocity becomes significant relative to its divergence. 293 In practice we will ignore the second half of the iteration step as long as 294 gross 2843 \begin{equation}\label{STOKES LARGE BV1} 295 jfenwick 3295 \|Bv_{1}\|_{0} \le \theta \cdot \|v_{1}-v_{0}\| 296 gross 2719 \end{equation} 297 caltinay 3325 where $0<\theta<1$ is a given factor. In this case we will also check the 298 stopping criterion with $v_{1}\rightarrow v_{2}$ but we will not correct $M$ 299 in this case. 300 gross 2719 301 jfenwick 3295 Starting from an initial guess $v_{0}$ and $p_{0}$ for velocity and pressure 302 caltinay 3325 the solution procedure is implemented as follows: 303 gross 2843 \begin{enumerate} 304 caltinay 3325 \item calculate viscosity $\eta(v_{0},p)_{0}$ and assemble operator $A$ from $\eta$ 305 \item calculate the tolerance $\tau_{1}$ from \eqn{STOKES TOL1} 306 \item solve \eqn{STOKES ITER STEP 1} for $dv_{1}$ with tolerance $\tau_{1}$ 307 jfenwick 3295 \item update $v_{1}= v_{0}+ dv_{1}$ 308 \item if $Bv_{1}$ is large (see~\ref{STOKES LARGE BV1}): 309 gross 2843 \begin{enumerate} 310 caltinay 3325 \item calculate the tolerance $\tau_{2}$ from~\ref{STOKES TOL2} 311 \item solve~\ref{STOKES ITER STEP 2} for $dp_{2}$ and $v_{2}$ with tolerance $\tau_{2}$ 312 \item update $p_{2}\leftarrow p_{0}+ dp_{2}$ 313 gross 2843 \end{enumerate} 314 \item else: 315 caltinay 3325 \begin{itemize} 316 \item update $p_{2}\leftarrow p$ and $v_{2}\leftarrow v_{1}$ 317 \end{itemize} 318 \item calculate convergence measure $\epsilon$ and convergence rate $\chi$ 319 \item if stopping criterion~\ref{STOKES STOPPING CRITERIA} holds: 320 \begin{itemize} 321 \item return $v_{2}$ and $p_{2}$ 322 \end{itemize} 323 \item else: 324 gross 2843 \begin{enumerate} 325 caltinay 3325 \item update $M$ and $K$ 326 \item goto step 1 with $v_{0}\leftarrow v_{2}$ and $p_{0}\leftarrow p_{2}$. 327 gross 2843 \end{enumerate} 328 \end{enumerate} 329 gross 2719 330 \subsection{Functions} 331 332 \begin{classdesc}{StokesProblemCartesian}{domain} 333 caltinay 3325 opens the Stokes problem\index{Stokes problem} on the \Domain domain. 334 The domain needs to support LBB compliant elements for the Stokes problem, see~\cite{LBB} for details\index{LBB condition}. 335 For instance one can use second order polynomials for velocity and first order 336 polynomials for the pressure on the same element. 337 Alternatively, one can use macro elements\index{macro elements} using linear 338 polynomials for both pressure and velocity with a subdivided element for the 339 velocity. Typically, the macro element is more cost effective. 340 The fact that pressure and velocity are represented in different ways is 341 expressed by 342 gross 2793 \begin{python} 343 caltinay 3325 velocity=Vector(0.0, Solution(mesh)) 344 pressure=Scalar(0.0, ReducedSolution(mesh)) 345 gross 2793 \end{python} 346 gross 2719 \end{classdesc} 347 348 caltinay 3325 \begin{methoddesc}[StokesProblemCartesian]{initialize}{\optional{f=Data(), 349 \optional{fixed_u_mask=Data(), \optional{eta=1,% 350 \optional{surface_stress=Data(), \optional{stress=Data()},% 351 \optional{restoration_factor=0}}}}}} 352 gross 2719 assigns values to the model parameters. In any call all values must be set. 353 \var{f} defines the external force $f$, \var{eta} the viscosity $\eta$, 354 \var{surface_stress} the surface stress $s$ and \var{stress} the initial stress $\sigma$. 355 caltinay 3325 The locations and components where the velocity is fixed are set by the values 356 of \var{fixed_u_mask}. 357 \var{restoration_factor} defines the restoring force factor $\alpha$. 358 The method will try to cast the given values to appropriate \Data class objects. 359 gross 2719 \end{methoddesc} 360 361 \begin{methoddesc}[StokesProblemCartesian]{solve}{v,p 362 \optional{, max_iter=100 \optional{, verbose=False \optional{, usePCG=True }}}} 363 caltinay 3325 solves the problem and returns approximations for velocity and pressure. 364 The arguments \var{v} and \var{p} define initial guesses. 365 \var{v} must have function space \var{Solution(domain)} and \var{p} must have 366 function space \var{ReducedSolution(domain)}. 367 The values of \var{v} marked by \var{fixed_u_mask} remain unchanged. 368 If \var{usePCG} is set to \True then the preconditioned conjugate gradient 369 method (PCG)\index{preconditioned conjugate gradient method!PCG} scheme is 370 used. Otherwise the problem is solved with the generalized minimal residual 371 method (GMRES)\index{generalized minimal residual method!GMRES}. 372 In most cases the PCG scheme is more efficient. 373 \var{max_iter} defines the maximum number of iteration steps. 374 If \var{verbose} is set to \True information on the progress of of the solver 375 is printed. 376 gross 2719 \end{methoddesc} 377 378 \begin{methoddesc}[StokesProblemCartesian]{setTolerance}{\optional{tolerance=1.e-4}} 379 caltinay 3325 sets the tolerance in an appropriate norm relative to the right hand side. 380 The tolerance must be non-negative and less than 1. 381 gross 2719 \end{methoddesc} 382 caltinay 3325 383 gross 2719 \begin{methoddesc}[StokesProblemCartesian]{getTolerance}{} 384 returns the current relative tolerance. 385 \end{methoddesc} 386 caltinay 3325 387 gross 2719 \begin{methoddesc}[StokesProblemCartesian]{setAbsoluteTolerance}{\optional{tolerance=0.}} 388 caltinay 3325 sets the absolute tolerance for the error in the relevant norm. 389 The tolerance must be non-negative. 390 Typically the absolute tolerance is set to 0. 391 gross 2719 \end{methoddesc} 392 393 \begin{methoddesc}[StokesProblemCartesian]{getAbsoluteTolerance}{} 394 gross 2851 returns the current absolute tolerance. 395 gross 2719 \end{methoddesc} 396 397 \begin{methoddesc}[StokesProblemCartesian]{getSolverOptionsVelocity}{} 398 caltinay 3325 returns the solver options used to solve \eqn{STOKES ITER STEP 1} and \eqn{STOKES P OPERATOR}) for velocity. 399 gross 2719 \end{methoddesc} 400 401 \begin{methoddesc}[StokesProblemCartesian]{getSolverOptionsPressure}{} 402 caltinay 3325 returns the solver options used to solve the preconditioner \eqn{STOKES P PREC} for pressure. 403 gross 2719 \end{methoddesc} 404 405 \begin{methoddesc}[StokesProblemCartesian]{getSolverOptionsDiv}{} 406 caltinay 3325 sets the solver options for solving the equation to project the divergence of 407 the velocity onto the function space of pressure. 408 gross 2719 \end{methoddesc} 409 410 gross 3582 \begin{methoddesc}[StokesProblemCartesian]{setStokesEquation}{\optional{f=None, 411 \optional{fixed_u_mask=None, \optional{eta=None,% 412 \optional{surface_stress=None, \optional{stress=None},% 413 \optional{restoration_factor=None}}}}}} 414 assigns new values to the model parameters, see method \function{initialize} for description of the 415 parameter list. In contrast to \function{initialize} only values given in the argument list are set. 416 Typically this method is called to update parameters such as viscosity $\eta$ within a time integration scheme 417 after the problem has been set up by an initial call of method \function{initialize}. 418 \end{methoddesc} 419 420 \begin{methoddesc}[StokesProblemCartesian]{updateStokesEquation}{v, p} 421 this method is called by the solver to allow for updating the model parameter during the iteration process for 422 incompressibility. In order to implement a non-linear dependence of model parameters such as viscosity $\eta$ 423 from the current velocity \var{v} or pressure field \var{p} this function can be overwritten in the following way: 424 \begin{python} 425 class MyStokesProblem(StokesProblemCartesian): 426 def updateStokesEquation(self, v, p): 427 my_eta= 428 self.setStokesEquation(eta=my_eta) 429 \end{python} 430 Note that \function{setStokesEquation} to update the model. \warning{It is not guaranteed that the iteration converges if model parameters are altered.} 431 \end{methoddesc} 432 433 caltinay 3325 \subsection{Example: Lid-driven Cavity} 434 The following script \file{lid_driven_cavity.py}\index{scripts!\file{lid_driven_cavity.py}} 435 which is available in the \ExampleDirectory illustrates the usage of the 436 \class{StokesProblemCartesian} class to solve the lid-driven cavity problem: 437 gross 2719 \begin{python} 438 caltinay 3325 from esys.escript import * 439 from esys.finley import Rectangle 440 caltinay 3348 from esys.weipa import saveVTK 441 caltinay 3325 from esys.escript.models import StokesProblemCartesian 442 NE=25 443 dom = Rectangle(NE,NE,order=2) 444 x = dom.getX() 445 sc=StokesProblemCartesian(dom) 446 mask= (whereZero(x[0])*[1.,0]+whereZero(x[0]-1))*[1.,0] + \ 447 (whereZero(x[1])*[0.,1.]+whereZero(x[1]-1))*[1.,1] 448 sc.initialize(eta=.1, fixed_u_mask=mask) 449 v=Vector(0., Solution(dom)) 450 v[0]+=whereZero(x[1]-1.) 451 p=Scalar(0.,ReducedSolution(dom)) 452 v,p=sc.solve(v, p, verbose=True) 453 saveVTK("u.vtu", velocity=v, pressure=p) 454 gross 2719 \end{python} 455 caltinay 3325