/[escript]/trunk/doc/user/stokessolver.tex
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log


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}

  ViewVC Help
Powered by ViewVC 1.1.26