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

Revision 3325 - (show annotations)
Fri Oct 29 00:30:51 2010 UTC (11 years, 8 months ago) by caltinay
File MIME type: application/x-tex
File size: 19440 byte(s)
Section 6.1.


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