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

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


 1 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 3 % Copyright (c) 2003-2018 by The University of Queensland 4 5 % 6 % Primary Business: Queensland, Australia 7 % Licensed under the Apache License, version 2.0 8 9 % 10 % Development until 2012 by Earth Systems Science Computational Center (ESSCC) 11 % Development 2012-2013 by School of Earth Sciences 12 % Development from 2014 by Centre for Geoscience Computing (GeoComp) 13 % 14 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 15 16 \section{The Stokes Problem} 17 \label{STOKES PROBLEM} 18 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 \begin{equation}\label{Stokes 1} 23 -\left(\eta(v_{i,j}+ v_{j,i})\right)_{,j}+p_{,i}=f_{i}-\sigma_{ij,j} 24 \end{equation} 25 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 30 We assume an incompressible medium: 31 \begin{equation}\label{Stokes 2} 32 -v_{i,i}=0 33 \end{equation} 34 Natural boundary conditions are taken in the form 35 \begin{equation}\label{Stokes Boundary} 36 \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 \end{equation} 38 which can be overwritten by constraints of the form 39 \begin{equation}\label{Stokes Boundary0} 40 v_{i}(x)=v^D_{i}(x) 41 \end{equation} 42 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 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 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 \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 \label{STOKES} 65 \end{equation} 66 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 71 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 \begin{equation}\label{STOKES ITER STEP 1} 75 A dv_{1} = G - A v_{0} - B^{*} p_{0} 76 \end{equation} 77 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 \begin{equation} 81 \begin{array}{rcl} 82 B A^{-1} B^{*} dp_{2} & = & Bv_{1} \\ 83 A dv_{2} & = & B^{*} dp_{2} 84 \end{array} 85 \label{STOKES ITER STEP 2} 86 \end{equation} 87 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 This solution strategy is called the Uzawa scheme\index{Uzawa scheme}. 90 91 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 an iterative solver for $A$ is applied. Another issue is the fact that the 97 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 100 In the following we will use the two norms 101 \begin{equation} 102 \|v\|_{1}^2 = \int_{\Omega} v_{j,k}v_{j,k} \; dx 103 \mbox{ and } 104 \|p\|_{0}^2= \int_{\Omega} p^2 \; dx. 105 \label{STOKES STOP} 106 \end{equation} 107 for velocity $v$ and pressure $p$. 108 The iteration is terminated if the stopping criterion 109 \begin{equation} \label{STOKES STOPPING CRITERIA} 110 \max(\|Bv_{1}\|_{0},\|v_{2}-v_{0}\|_{1}) \le \tau \cdot \|v_{2}\|_{1} 111 \end{equation} 112 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 117 We want to optimize the tolerance choice for solving~\ref{STOKES ITER STEP 1} 118 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 \begin{equation} \label{STOKES total V1} 123 v_{1} = e_{1} + v_{0} + A^{-1} ( G - Av_{0} - B^{*} p_{0} ) 124 \end{equation} 125 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 \begin{equation} 130 \| A e_{1} \|_{s} = \| G - A v_{1} - B^{*} p_{0} \|_{s} \le \tau_{1} \| G - A v_{0} - B^{*} p_{0} \|_{s} 131 \end{equation} 132 where $\tau_{1}$ is the tolerance which we need to choose. 133 This translates into the condition 134 \begin{equation} 135 \| e_{1} \|_{1} \le K \tau_{1} \| dv_{1} - e_{1} \|_{1} 136 \end{equation} 137 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 From the first equation of~\ref{STOKES ITER STEP 2} we have 140 \begin{equation}\label{STOKES total P2} 141 p_{2} = p_{0} + (B A^{-1} B^{*})^{-1} (e_{2} + Bv_{1} ) 142 \end{equation} 143 where $e_{2}$ represents the error when solving~\ref{STOKES ITER STEP 2}. 144 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 \begin{equation} \label{STOKES P PREC} 151 \frac{1}{\eta} (Pq) = q \; . 152 \end{equation} 153 Note that in each evaluation of the iteration operator $q=B A^{-1} B^{*} s$ 154 one needs to solve the problem 155 \begin{equation} \label{STOKES P OPERATOR} 156 A w = B^{*} s 157 \end{equation} 158 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 162 In an implementation we use the fact that the residual $r$ is given as 163 \begin{equation} \label{STOKES RES } 164 r= B (v_{1} - A^{-1} B^{*} dp) = B (v_{1} - A^{-1} B^{*} dp) = B (v_{1}-dv_{2}) = B v_{2} 165 \end{equation} 166 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 In PCG the iteration is terminated if 172 \begin{equation} \label{STOKES P OPERATOR ERROR} 173 \| P^{\frac{1}{2}}B v_{2} \|_{0} \le \tau_{2} \| P^{\frac{1}{2}}B v_{1} \|_{0} 174 \end{equation} 175 where $\tau_{2}$ is the given tolerance. This translates into 176 \begin{equation} \label{STOKES P OPERATOR ERROR 2} 177 \|e_{2}\|_{0} = \| B v_{2} \|_{0} \le M \tau_{2} \| B v_{1} \|_{0} 178 \end{equation} 179 where $M$ is taking care of the fact that $P^{\frac{1}{2}}$ is dropped. 180 181 As we assume that there is no significant error from solving with the operator 182 $A$ we have 183 \begin{equation} \label{STOKES total V2} 184 v_{2} = v_{1} - dv_{2} 185 = v_{1} - A^{-1} B^{*}dp 186 \end{equation} 187 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 \begin{equation} 191 v = \Phi(v,p) \mbox{ and } p = \Psi(u,p) 192 \end{equation} 193 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 \begin{equation} 199 \begin{array}{rcl} 200 \delta p & = & (B A^{-1} B^{*})^{-1} ( e_{2} + B e_{1} ) \\ 201 \delta v & = & e_{1} - A^{-1} B^{*}\delta p \;. 202 \end{array}\label{STOKES ERRORS} 203 \end{equation} 204 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 209 To measure convergence we use 210 \begin{equation} 211 \epsilon = \max(\|v_{2}-v\|_{1}, \|B A^{-1} B^{*} (p_{2}-p)\|_{0}) 212 \end{equation} 213 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 \begin{equation} 217 \epsilon = \max(\|v_{2}-v_{0}\|_{1}, \|B v_{1}\|_{0}) 218 \end{equation} 219 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 \begin{equation} \label{STOKES EST 1} 226 \begin{array}{rcl} 227 \|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 \end{array} 234 \end{equation} 235 So we choose the value for $\tau_{2}$ from 236 \begin{equation} \label{STOKES TOL2} 237 M \tau_{2} \|B v_{1}\|_{0} \le (\chi^{-})^2 \epsilon^{-} 238 \end{equation} 239 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 \begin{equation}\label{STOKES EST 2} 242 \begin{array}{rcl} 243 \|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 \\ 248 & \le & ( 1 + K \tau_{1}) \chi^{-} \cdot \epsilon^{-} 249 \end{array} 250 \end{equation} 251 So we choose the value for $\tau_{1}$ from 252 \begin{equation} \label{STOKES TOL1} 253 K \tau_{1} \le \chi^{-} 254 \end{equation} 255 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 \begin{equation}\label{STOKES EST 3} 266 \chi \le ( 1 + \max(M \frac{\tau_{2} \|B v_{1}\|_{0}}{\chi^{-} \epsilon^{-}},K \tau_{1} ) ) \cdot \chi^{-} 267 \end{equation} 268 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 \begin{equation}\label{STOKES EST 3b} 271 \chi = ( 1 + \max(M^{+} \frac{\chi^{-}}{M},K^{+} \frac{\chi^{-}}{K}) ) \cdot \chi^{-} 272 \end{equation} 273 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 \begin{equation} 276 K^{+} = \frac{\chi-\chi^{-}}{(\chi^{-})^2} K 277 \end{equation} 278 In practice we will use 279 \begin{equation} 280 K^{+} = \max(\frac{\chi-\chi^{-}}{(\chi^{-})^2} K,\frac{1}{2}K,1) 281 \end{equation} 282 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 286 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 \begin{equation}\label{STOKES LARGE BV1} 295 \|Bv_{1}\|_{0} \le \theta \cdot \|v_{1}-v_{0}\| 296 \end{equation} 297 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 301 Starting from an initial guess $v_{0}$ and $p_{0}$ for velocity and pressure 302 the solution procedure is implemented as follows: 303 \begin{enumerate} 304 \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 \item update $v_{1}= v_{0}+ dv_{1}$ 308 \item if $Bv_{1}$ is large (see~\ref{STOKES LARGE BV1}): 309 \begin{enumerate} 310 \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 \end{enumerate} 314 \item else: 315 \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 \begin{enumerate} 325 \item update $M$ and $K$ 326 \item goto step 1 with $v_{0}\leftarrow v_{2}$ and $p_{0}\leftarrow p_{2}$. 327 \end{enumerate} 328 \end{enumerate} 329 330 \subsection{Functions} 331 332 \begin{classdesc}{StokesProblemCartesian}{domain} 333 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 \begin{python} 343 velocity=Vector(0.0, Solution(mesh)) 344 pressure=Scalar(0.0, ReducedSolution(mesh)) 345 \end{python} 346 \end{classdesc} 347 348 \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 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 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 \end{methoddesc} 360 361 \begin{methoddesc}[StokesProblemCartesian]{solve}{v,p 362 \optional{, max_iter=100 \optional{, verbose=False \optional{, usePCG=True }}}} 363 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 \end{methoddesc} 377 378 \begin{methoddesc}[StokesProblemCartesian]{setTolerance}{\optional{tolerance=1.e-4}} 379 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 \end{methoddesc} 382 383 \begin{methoddesc}[StokesProblemCartesian]{getTolerance}{} 384 returns the current relative tolerance. 385 \end{methoddesc} 386 387 \begin{methoddesc}[StokesProblemCartesian]{setAbsoluteTolerance}{\optional{tolerance=0.}} 388 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 \end{methoddesc} 392 393 \begin{methoddesc}[StokesProblemCartesian]{getAbsoluteTolerance}{} 394 returns the current absolute tolerance. 395 \end{methoddesc} 396 397 \begin{methoddesc}[StokesProblemCartesian]{getSolverOptionsVelocity}{} 398 returns the solver options used to solve \eqn{STOKES ITER STEP 1} and \eqn{STOKES P OPERATOR}) for velocity. 399 \end{methoddesc} 400 401 \begin{methoddesc}[StokesProblemCartesian]{getSolverOptionsPressure}{} 402 returns the solver options used to solve the preconditioner \eqn{STOKES P PREC} for pressure. 403 \end{methoddesc} 404 405 \begin{methoddesc}[StokesProblemCartesian]{getSolverOptionsDiv}{} 406 sets the solver options for solving the equation to project the divergence of 407 the velocity onto the function space of pressure. 408 \end{methoddesc} 409 410 \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 \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 \begin{python} 438 from esys.escript import * 439 from esys.finley import Rectangle 440 from esys.weipa import saveVTK 441 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)*[1.,0]+whereZero(x-1))*[1.,0] + \ 447 (whereZero(x)*[0.,1.]+whereZero(x-1))*[1.,1] 448 sc.initialize(eta=.1, fixed_u_mask=mask) 449 v=Vector(0., Solution(dom)) 450 v+=whereZero(x-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 \end{python} 455