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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2793 - (show annotations)
Tue Dec 1 06:10:10 2009 UTC (9 years, 11 months ago) by gross
File MIME type: application/x-tex
File size: 25596 byte(s)
some new util functions
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\hackscore{i,j}+ v\hackscore{i,j})\right)\hackscore{,j}+p\hackscore{,i}=f\hackscore{i}-\sigma\hackscore{ij,j}
8 \end{equation}
9 where $f\hackscore{i}$ defines an internal force \index{force, internal} and $\sigma\hackscore{ij}$ is an intial 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\hackscore{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\hackscore{i,j}+ v\hackscore{i,j})\right)n\hackscore{j}-n\hackscore{i}p=s\hackscore{i} - \alpha \cdot n\hackscore{i} n\hackscore{j} v\hackscore{j}+\sigma\hackscore{ij} n\hackscore{j}
18 \end{equation}
19 which can be overwritten by constraints of the form
20 \begin{equation}\label{Stokes Boundary0}
21 v\hackscore{i}(x)=v^D\hackscore{i}(x)
22 \end{equation}
23 at some locations $x$ at the boundary of the domain. $s\hackscore{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{SADDLEPOINT}
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\hackscore{0}$ and $p\hackscore{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{SADDLEPOINT}
52 \begin{equation}\label{SADDLEPOINT ITER STEP 1}
53 A dv\hackscore{1} = G - A v\hackscore{0} - B^{*} p\hackscore{0}
54 \end{equation}
55 We then insert the new approximation $v\hackscore{1}=v\hackscore{0}+dv\hackscore{1}$ to calculate a correction $dp\hackscore{2}$
56 for the pressure and an additional correction $dv\hackscore{2}$ for the velocity by solving
57 \begin{equation}
58 \begin{array}{rcl}
59 B A^{-1} B^{*} dp\hackscore{2} & = & Bv\hackscore{1} \\
60 A dv\hackscore{2} & = & B^{*} dp\hackscore{2}
61 \end{array}
62 \label{SADDLEPOINT ITER STEP 2}
63 \end{equation}
64 The new velocity and pressure are then given by $v=v\hackscore{1}-dv\hackscore{2}$ and
65 $p=p\hackscore{0}+dp\hackscore{2}$ which will fullfill the block system~\ref{SADDLEPOINT}.
66 This solution strategy is called the Uzawa scheme \index{Uzawa scheme}. There are
67 two problems with this scheme. Firstly, in practice we will use an iterative scheme
68 to solve the problem for operator $A$. So we will be unable to calculate the operator
69 $ B A^{-1} B^{*}$ required for $dp\hackscore{2}$. In fact, we need to use an iterative scheme
70 to solve the first equation in~\ref{SADDLEPOINT ITER STEP 2} where in each iteration step
71 an iterative solver for $A$ is applied. The second issue is that
72 viscosity $\eta$ may depend of velocity or pressure and so we need to iterate over the
73 three equations~\ref{SADDLEPOINT ITER STEP 1} and~\ref{SADDLEPOINT ITER STEP 2}. Using the
74 two norms
75 \begin{equation}
76 \|v\|\hackscore{1}^2 = \int\hackscore{\Omega} v\hackscore{j,k}v\hackscore{j,k} \; dx
77 \mbox{ and }
78 \|p\|\hackscore{0}^2= \int\hackscore{\Omega} p^2 \; dx.
79 \label{STOKES STOP}
80 \end{equation}
81 to define the stopping criterium
82 \begin{equation}
83 \max(\|Bv\hackscore{1}\|\hackscore{0},\|v-v\hackscore{0}\|\hackscore{1}) \le \tau \cdot \|v\|\hackscore{1}
84 \end{equation}
85 to terminate the overall iteration with a overall given tolerance $0<\tau<1$.
86 Notice that because of the first equation of
87 and~\ref{SADDLEPOINT ITER STEP 2} $\|Bv\hackscore{1}\|\hackscore{0}$ is equal to the
88 norm of $B A^{-1} B^{*} dp\hackscore{2}$ and consequently provides a norm for the pressure correction.
89
90 Let us first have a look at the solution process~\ref{SADDLEPOINT ITER STEP 1} which needs to resolve the
91 non-linearity of the viscosity coeficient. We will run a iterative process by resolving
92 the equation~\ref{SADDLEPOINT ITER STEP 1} with $v\hackscore{0} \leftarrow v\hackscore{1}$. The first observation is that there is no point start
93 the iteration process over the second set of equation if $\|Bv\hackscore{1}\|\hackscore{0}$
94 is smaller then $\|v\hackscore{1}-v\hackscore{0}\|\hackscore{1}$. So we will terminate the non-linear iteration
95 process of the first equation if
96 \begin{equation}
97 \theta \cdot \|v\hackscore{1}-v\hackscore{0}\| \le \|Bv\hackscore{1}\|\hackscore{0}
98 \mbox{ or }
99 \|v\hackscore{1}-v\hackscore{0}\|\hackscore{1} \le \tau \cdot \|v\hackscore{1}\|\hackscore{1}
100 \end{equation}
101 where $0<\theta<1$ is a given factor (typically $\theta=0.5$).
102
103 We need to think about appropriate stopping criteria when solving
104 \ref{SADDLEPOINT ITER STEP 1}. We would like to use very weak convergence criteria to reduce computational costs but they need to be tight enough to not interfere with the
105 convergence of the iteration process one level above.
106 We will use a sparse matrix solver so we have not full control on the norm $\|.\|\hackscore{s}$ used in the stopping criteria
107 \begin{equation}
108 \| G - A v\hackscore{1} - B^{*} p\hackscore{0} \|\hackscore{s} \le \tau\hackscore{1} \| G - A v\hackscore{0} - B^{*} p\hackscore{0} \|\hackscore{s}
109 \end{equation}
110 If $e\hackscore{1}$ is the error of the returned solution $dv\hackscore{1}$ this translates into the condition
111 \begin{equation}
112 \| e\hackscore{1} \|\hackscore{1} \le K \tau\hackscore{1} \| dv\hackscore{1} - e\hackscore{1} \|\hackscore{1}
113 \end{equation}
114 The constant $K$ represents some uncertainty combining a variety of unknown factors such as the
115 solver being used and the condition number of the stiffness matrix. This leads to the estimate
116 \begin{equation}
117 \| dv\hackscore{1} \|\hackscore{1} \le (1+K \tau\hackscore{1}) \chi \| dv\hackscore{1}^- \|\hackscore{1}
118 \label{DV COND}
119 \end{equation}
120 where $chi$ is the convergence rate of the \ref{SADDLEPOINT ITER STEP 1} iteration and $dv\hackscore{1}^-$ the correction from the last step. In order to make sure that the error $e\hackscore{1}$ can be neglected we need
121 to have $K \tau\hackscore{1} \le \chi$ which leads to the condition
122 \begin{equation}
123 \tau\hackscore{1} \le \frac{\min(\chi,\frac{1}{2})}{K}
124 \label{NEW TAU 1}
125 \end{equation}
126 We can estimate $K$ by comparing estimates for the convergence rate $\chi=\frac{\| dv\hackscore{1} \|\hackscore{1}}{\| dv\hackscore{1}^- \|\hackscore{1}}$ by checking our error prediction~\ref{DV COND}. If we have
127 access to the estimate of the convergence rate $\chi^-$ of the previous iteration step we
128 get
129 \begin{equation}
130 \frac{1}{K} = \frac{\chi^-}{\chi-\chi^{-}} \tau\hackscore{1}^-
131 \end{equation}
132 if $\chi \ge \chi^- (1 + \chi^-) $ which means that our prediction for a suitable tolerane based on the avaibale $K$ value was wrong.
133 So we start with the value $K=1$ in~\ref{NEW TAU 1} and recalculate
134 $K$ if we have underestimate $K$.
135
136
137
138
139
140
141
142 If $dv\hackscore{1}$
143 From Equation~\ref{SADDLEPOINT ITER STEP 1} we have
144 \begin{equation}
145 v\hackscore{1} = e\hackscore{1} + A^{-1} ( G - B^{*} p )
146 \end{equation}
147 This translates into the conditoon
148 \begin{equation}
149 \| e\hackscore{1} \|\hackscore{1} \le K \tau\hackscore{1} \| dv\hackscore{1} - e\hackscore{1} \|\hackscore{1}
150 \mbox{ therefore }
151 \| e\hackscore{1} \|\hackscore{1} \le \frac{K \tau\hackscore{1}}{1-K \tau\hackscore{1}} \| dv\hackscore{1} \|\hackscore{1}
152 \end{equation}
153 The constant $K$ represents some uncertainty combining a variety of unknown factors such as the
154 solver being used and the condition number of the stiffness matrix.
155
156
157
158
159
160 ===================================
161
162
163 We use iterative techniques to solve this problem: Given an approximation $v$ and $p$ for
164 velocity and pressure we perform the following steps in the Uzawa scheme \index{Uzawa scheme} style:
165 \begin{enumerate}
166 \item calculate viscosity $\eta(v,p)$ and assemble operator $A$ from $\eta$.
167 \item Solve for $dv$:
168 \begin{equation}
169 A dv = G - A v - B^{*} p \label{SADDLEPOINT ITER STEP 1}
170 \end{equation}
171 \item update $v\hackscore{1}= v+ dv$
172 \item if $\max(\|Bv\hackscore{1}\|\hackscore{0},\|dv\|\hackscore{1}) \le \tau \cdot \|v\hackscore{1}\|\hackscore{1}$: return v$\hackscore{1},p$
173 \item Solve for $dp$:
174 \begin{equation}
175 \begin{array}{rcl}
176 B A^{-1} B^{*} dp & = & Bv\hackscore{1} \\
177 A dv\hackscore{2} & = & B^{*} dp
178 \end{array}
179 \label{SADDLEPOINT ITER STEP 2}
180 \end{equation}
181 \item update $p\hackscore{2}\leftarrow p+ dp$ and $v\hackscore{2}= v\hackscore{1} - dv\hackscore{2}$
182 \item goto Step 0 with $v=v\hackscore{2}$ and $p=p\hackscore{2}$.
183 \end{enumerate}
184 where $\tau$ is the given tolerance. Notice that the operator $A$ if it is depending on $v$ or $p$ is updated once only.
185 In this algorithm
186 \begin{equation}
187 \|v\|\hackscore{1}^2 = \int\hackscore{\Omega} v\hackscore{j,k}v\hackscore{j,k} \; dx
188 \mbox{ and }
189 \|p\|\hackscore{0}^2= \int\hackscore{\Omega} p^2 \; dx.
190 \label{STOKES STOP}
191 \end{equation}
192 so the stopping criterion used is checking for convergence as well as for
193 the fact that the incompressiblity condition~\ref{Stokes 2} is fullfilled.
194
195 To solve the update step~\ref{SADDLEPOINT ITER STEP 2} we use an iterative solver with iteration
196 operator $B A^{-1} B^{*}$, eg. using generalized minimal residual method (GMRES) \index{linear solver!GMRES}\index{GMRES}, preconditioned conjugate gradient method (PCG) \index{linear solver!PCG}\index{PCG}. As suitable preconditioner \index{preconditioner} for this operator is $\frac{1}{\eta}$, ie
197 the evaluation of the preconditioner $P$ for a given pressure increment $q$ is the solution of
198 \begin{equation} \label{P PREC}
199 \frac{1}{\eta} Pq = q \; .
200 \end{equation}
201 Note that in each evaluation of the iteration operator $q=B A^{-1} B^{*} s$ one needs to solve
202 the problem
203 \begin{equation} \label{P OPERATOR}
204 A w = B^{*} s
205 \end{equation}
206 with sufficient accuracy to return $q=Bw$. Notice that the residual $r$ is given as
207 \begin{equation} \label{STOKES RES }
208 r= B (v\hackscore{1} - A^{-1} B^{*} dp) = B (v\hackscore{1} - A^{-1} B^{*} dp) = B (v\hackscore{1}-dv\hackscore{2}) = B v\hackscore{2}
209 \end{equation}
210 so in fact the residual $r$ is represented by the updated velocity $v\hackscore{2}$. This saves the recovery of
211 $dv\hackscore{2}$ in~\ref{SADDLEPOINT ITER STEP 2} after $dp$ has been calculated as iterative method such as PCG calculate the solution approximations along with the their residual. In PCG the iteration is terminated if
212 \begin{equation} \label{P OPERATOR}
213 \| P^{\frac{1}{2}}B v\hackscore{2} \|\hackscore{0} \le \tau\hackscore{2} \| P^{\frac{1}{2}}B v\hackscore{1} \|\hackscore{0}
214 \end{equation}
215 where $\tau\hackscore{2}$ is the given tolerane. The update step~\ref{P OPERATOR} involves the
216 solution of a sparse matrix problem. We use the tolerance $\tau\hackscore{2}^2$ in order to make sure that any
217 error from solving this problem does not interfere with the PCG iteration.
218
219
220
221 We need to think about appropriate stopping criteria when solving
222 \ref{SADDLEPOINT ITER STEP 1}, \ref{SADDLEPOINT ITER STEP 2} and~\ref{P OPERATOR}. We would like to use very weak convergence criteria to reduce computational costs but they need to be tight enough to not interfere with the
223 convergence of the iteration process one level above. From Equation~\ref{SADDLEPOINT ITER STEP 1} we have
224 \begin{equation}
225 v\hackscore{1} = e\hackscore{1} + A^{-1} ( G - B^{*} p )
226 \end{equation}
227 We will use a sparse matrix solver so we have not full control on the norm $\|.\|\hackscore{s}$ used in the stopping criteria
228 \begin{equation}
229 \| G - A v\hackscore{1} - B^{*} p \|\hackscore{s} \le \tau\hackscore{1} \| G - A v - B^{*} p \|\hackscore{s}
230 \end{equation}
231 This translates into the conditoon
232 \begin{equation}
233 \| e\hackscore{1} \|\hackscore{1} \le K \tau\hackscore{1} \| dv\hackscore{1} - e\hackscore{1} \|\hackscore{1}
234 \mbox{ therefore }
235 \| e\hackscore{1} \|\hackscore{1} \le \frac{K \tau\hackscore{1}}{1-K \tau\hackscore{1}} \| dv\hackscore{1} \|\hackscore{1}
236 \end{equation}
237 The constant $K$ represents some uncertainty combining a variety of unknown factors such as the
238 solver being used and the condition number of the stiffness matrix.
239
240 From the first equation of~\ref{SADDLEPOINT ITER STEP 2} we have
241 \begin{equation}
242 p\hackscore{2} = p + (B A^{-1} B^{*})^{-1} (e\hackscore{2} + Bv\hackscore{1} ) =
243 (B A^{-1} B^{*})^{-1} ( e\hackscore{2} + B e\hackscore{1} + B A^{-1} G)
244 \end{equation}
245 and simlar
246 \begin{equation}
247 v\hackscore{2} = v\hackscore{1} - dv\hackscore{2}
248 = v\hackscore{1} - A^{-1} B^{*}dp
249 = e\hackscore{1} + A^{-1} ( G - B^{*} p ) - A^{-1} B^{*} (p\hackscore{2}-p)
250 = e\hackscore{1} + A^{-1} ( G - B^{*} p\hackscore{2})
251 \end{equation}
252 This shows that we can write the iterative process as a fixed point iteration to solve (assume all errors are zero)
253 \begin{equation}
254 v = \Phi(v,p) \mbox{ and } p = \Psi(u,p)
255 \end{equation}
256 where
257 \begin{equation}
258 \begin{array}{rcl}
259 \Psi(v,p) & = & (B A^{-1} B^{*})^{-1} B A^{-1} G \\
260 \Phi(u,p) & = & A^{-1} ( G - B^{*} (B A^{-1} B^{*})^{-1} B A^{-1} G )
261 \end{array}
262 \end{equation}
263 Notice that if $A$ is independent from $v$ and $p$ the operators $\Phi(v,p)$ and $\Psi(u,p)$ are constant
264 and threfore the iteration will terminate - providing no termination errors in sub-iterations - after one step.
265 We also can give a formula for the error
266 \begin{equation}
267 \begin{array}{rcl}
268 \delta p & = & (B A^{-1} B^{*})^{-1} ( e\hackscore{2} + B e\hackscore{1} ) \\
269 \delta v & = & e\hackscore{1} - A^{-1} B^{*}\delta p \;.
270 \end{array}\label{STOKES ERRORS}
271 \end{equation}
272 Notice that $B\delta v = - e\hackscore{2}=-Bv\hackscore{2}$.
273
274 With this notation
275 \begin{equation}
276 \begin{array}{rcl}
277 v\hackscore{2} = \Phi(v,p) + \delta v \\
278 p\hackscore{2} = \Psi(v,p) + \delta p \\
279 \end{array}
280 \end{equation}
281 We use the values $dv=v\hackscore{2}-v$ and $Bv\hackscore{1}=B A^{-1} B^{*} dp$ to measure the error of the
282 current approximation $v\hackscore{2}$ and $v\hackscore{2}$ towards the exact solution.
283 Assuming that the iteration does converge with a convergence rate $\chi^{-}$ we have the estimate
284 \begin{equation}
285 \max(\|dv\|\hackscore{1}, \|Bv\hackscore{1}\|\hackscore{0})
286 \le \chi^{-} \max(\|dv^{-}\|\hackscore{1}, \|Bv\hackscore{1}^{-}\|\hackscore{0})
287 + \max(\|\delta v\|\hackscore{1}, \|(B A^{-1} B^{*}) \delta p\|\hackscore{0})
288 \end{equation}
289 were the upper index $-$ referes to the increments in the last step.
290
291 Now we will establish estimates for $\|\delta v\|$ and $\|Bv\hackscore{1}\|$ from
292 formulas~\ref{STOKES ERRORS} where neglect the $\delta p$ in the equation for $\delta v$ assuming that
293 $\delta p$ is controlled.
294 \begin{equation}
295 \|\delta v\|\hackscore{1} \approx \|e\hackscore{1}\|\hackscore{1} \le \frac{K \tau\hackscore{1}}{1-K \tau\hackscore{1}} \| dv\hackscore{1} \|\hackscore{1} \approx \frac{K \tau\hackscore{1}}{1-K \tau\hackscore{1}} \| dv\|\hackscore{1}
296 \end{equation}
297 and
298 \begin{equation}
299 \|(B A^{-1} B^{*}) \delta p\|\hackscore{0} \approx \| e\hackscore{2} \|\hackscore{0}
300 = \| B v\hackscore{2} \|\hackscore{0} \le M \tau\hackscore{2} \| B v\hackscore{1} \|\hackscore{0}
301 \end{equation}
302 which leads to
303 \begin{equation}
304 \max(\|dv\|\hackscore{1}, \|Bv\hackscore{1}\|\hackscore{0})
305 \le \frac{\chi^{-}}{1-max{(M \tau\hackscore{2}, \frac{K \tau\hackscore{1}}{1-K \tau\hackscore{1}})}} \max(\|dv^{-}\|\hackscore{1}, \|Bv\hackscore{1}^{-}\|\hackscore{0}) \label{STOKES ESTIMTED IMPROVEMENT}
306 \end{equation}
307 where the upper index $-$ refers to the previous iteration step.
308 If we allow require $max{(M \tau\hackscore{2}, \frac{K \tau\hackscore{1}}{1-K \tau\hackscore{1}})}\le \gamma$
309 with a given $0<\gamma<1$ we can set
310 \begin{equation} \label{STOKES SET TOL}
311 \tau\hackscore{2} = \frac{\gamma}{M} \mbox{ and } \tau\hackscore{1} = \frac{1}{K} \frac{\gamma}{1+\gamma}
312 \end{equation}
313 Problem is that we do not know $M$ and $K$ but can use the convergence rate
314 \begin{equation}
315 \chi := \frac{\max(\|dv\|\hackscore{1}, \|Bv\hackscore{1}\|\hackscore{0})}{\max(\|dv^{-}\|\hackscore{1}, \|Bv\hackscore{1}^{-}\|\hackscore{0})}
316 \end{equation}
317 to construct estimates of $M$ and $K$. We look at the two cases where our prediction and choice of
318 the tolerances was good or where we went wrong:
319 \begin{equation}
320 \chi \le \frac{\chi^{-}}{1-\gamma} \mbox{ or }
321 \chi = \frac{\chi^{-}}{1-\gamma\hackscore{0}} >\frac{\chi^{-}}{1-\gamma}.
322 \end{equation}
323 which translates to
324 \begin{equation}
325 \gamma\hackscore{0} \le \gamma \mbox{ or } \gamma\hackscore{0}=1-\frac{\chi^{-}}{ \chi}>\gamma>0
326 \end{equation}
327 In the second case use \ref{STOKES SET TOL} for $\gamma=\gamma\hackscore{0}$ to get new estimates $M^{+}$ and $K^{+}$
328 for $M$ and $K$:
329 \begin{equation} \label{TOKES CONST UPDATE}
330 M^{+} =
331 \frac{\gamma\hackscore{0}}{\gamma} M
332 \mbox{ and } K^{+} =
333 \frac{\gamma\hackscore{0}}{\gamma}
334 \frac{1+\gamma}{1+\gamma\hackscore{0}}
335 K
336 \mbox{ with } \gamma\hackscore{0}=\max(1-\frac{\chi^{-}}{ \chi},\gamma)
337 \end{equation}
338 With these updated constants we can now get better values for the tolerances via an updated value $\gamma^{+}$ for $\gamma$ and the corrected values $M^{+}$ and $K^{+}$ for $M$ and $K$. If we are in the case of convergence which we
339 meassure by
340 \begin{equation}
341 \chi < \chi\hackscore{max}
342 \end{equation}
343 where $\chi\hackscore{max}$ is given value with $0<\chi\hackscore{max}<1$. We then consider the following
344 criteria
345 \begin{itemize}
346 \item As we are in case of convergence we try to relax the tolerances by increasing $\gamma$ by a factor $\frac{1}{\omega}$ ($0<\omega<1$). So we would like to choose
347 $\gamma^{+} \ge \frac{\gamma}{\omega}$.
348 \item We need to make sure that the estimated convergence rate based on $\gamma^{+}$ will still maintain convergence
349 \begin{equation}
350 \frac{\chi}{1-\gamma^{+}} \le \chi\hackscore{max}
351 \end{equation}
352 which translates into
353 \begin{equation}
354 \gamma^{+} \le 1 - \frac{\chi}{\chi\hackscore{max}}
355 \end{equation}
356 Notice that the right hand side is positive.
357 \item We do not want to iterate more then necessary. So we would like reduce the error in the next just below the
358 desired tolerance:
359 \begin{equation}
360 \frac{\chi}{1-\gamma^{+}}\max(\|dv\|\hackscore{1}, \|Bv\hackscore{1}\|\hackscore{0}) \ge f \cdot \tau \|v\hackscore{2}\|\hackscore{1}
361 \end{equation}
362 where the left hand side provides an estimate for the error to prediced in the next step. $f$ is a
363 safty factor. This leads to
364 \begin{equation}
365 \gamma^{+} \le
366 1 -
367 \chi \frac{\max(\|dv\|\hackscore{1}, \|Bv\hackscore{1}\|\hackscore{0})} { f \cdot \tau \|v\hackscore{2}\|\hackscore{1}}
368 \end{equation}
369 \end{itemize}
370 Putting all these criteria together we get
371 \begin{equation}
372 \gamma^{+}=\min(\max(
373 \frac{\gamma}{\omega},
374 1 -
375 \chi \frac{\max(\|dv\|\hackscore{1}, \|Bv\hackscore{1}\|\hackscore{0})} { f \cdot \tau \|v\hackscore{2}\|\hackscore{1}}), 1 - \frac{\chi}{\chi\hackscore{max}})
376 \label{STOKES SET GAMMA PLUS}
377 \end{equation}
378 In case we cannot see convergence $\gamma$ is reduced by the factor $\omega$
379 \begin{equation}
380 \gamma^{+}=\max(\omega \cdot \gamma, \gamma\hackscore{min})\label{STOKES SET GAMMA PLUS2}
381 \end{equation}
382 where $\gamma\hackscore{min}$ is introduced to make sure that $\gamma^{+}$ stays away from zero.
383
384
385 Appling~\ref{STOKES SET TOL} for $\gamma^{+}$ and~\ref{TOKES CONST UPDATE} we get the update formula
386 \begin{equation} \label{STOKES CONST UPDATE}
387 \tau\hackscore{2}^{+} =
388 \frac{\gamma^{+}}{\gamma\hackscore{0}} \tau\hackscore{2}
389 \mbox{ and } \tau\hackscore{1}^{+} =
390 \frac{\gamma^{+}}{\gamma\hackscore{0}}
391 \frac{1+\gamma\hackscore{0}}{1+\gamma^{+}}
392 \tau\hackscore{1}
393 \end{equation}
394 to get the new tolerances $\tau\hackscore{1}^{+}$ and $\tau\hackscore{2}^{+}$ to be used in the next iteration step.
395 Notice that the coefficients $M$ and $K$ have been eliminated. The calculation of $\gamma\hackscore{0}$ requires
396 at least two iteration steps (or a previous time step).
397
398 Typical initial values are
399 \begin{equation} \label{STOKES VALUES}
400 \tau\hackscore{1}=\tau\hackscore{1}=\frac{1}{100} ; \; \omega=\frac{1}{2} \; \gamma=\frac{1}{10} ;\gamma\hackscore{min}=10^{-6}
401 \end{equation}
402 where after the first step we set $\gamma\hackscore{0}=\gamma$ as no $\chi^{-}$ is available.
403
404
405 \subsection{Functions}
406
407 \begin{classdesc}{StokesProblemCartesian}{domain}
408 opens the Stokes problem\index{Stokes problem} on the \Domain domain. The domain
409 needs to support LBB compliant elements for the Stokes problem, see~\cite{LBB} for detials~\index{LBB condition}.
410 For instance one can use second order polynomials for velocity and
411 first order polynomials for the pressure on the same element. Alternativly, one can use
412 macro elements\index{macro elements} using linear polynomial for both pressure and velocity bu with a subdivided
413 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
414 \begin{python}
415 velocity=Vector(0.0, Solution(mesh))
416 pressure=Scalar(0.0, ReducedSolution(mesh))
417 \end{python}
418 \end{classdesc}
419
420 \begin{methoddesc}[StokesProblemCartesian]{initialize}{\optional{f=Data(), \optional{fixed_u_mask=Data(), \optional{eta=1, \optional{surface_stress=Data(), \optional{stress=Data()}, \optional{
421 restoration_factor=0}}}}}}
422 assigns values to the model parameters. In any call all values must be set.
423 \var{f} defines the external force $f$, \var{eta} the viscosity $\eta$,
424 \var{surface_stress} the surface stress $s$ and \var{stress} the initial stress $\sigma$.
425 The locations and compontents where the velocity is fixed are set by
426 the values of \var{fixed_u_mask}. \var{restoration_factor} defines the restoring force factor $\alpha$.
427 The method will try to cast the given values to appropriate
428 \Data class objects.
429 \end{methoddesc}
430
431 \begin{methoddesc}[StokesProblemCartesian]{solve}{v,p
432 \optional{, max_iter=100 \optional{, verbose=False \optional{, usePCG=True }}}}
433 solves the problem and return approximations for velocity and pressure.
434 The arguments \var{v} and \var{p} define initial guess.
435 \var{v} must have function space \var{Solution(domain)} and
436 \var{p} must have function space \var{ReducedSolution(domain)}.
437 The values of \var{v} marked
438 by \var{fixed_u_mask} remain unchanged.
439 If \var{usePCG} is set to \True
440 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
441 the PCG scheme is more efficient.
442 \var{max_iter} defines the maximum number of iteration steps.
443
444 If \var{verbose} is set to \True informations on the progress of of the solver are printed.
445 \end{methoddesc}
446
447
448 \begin{methoddesc}[StokesProblemCartesian]{setTolerance}{\optional{tolerance=1.e-4}}
449 sets the tolerance in an appropriate norm relative to the right hand side. The tolerance must be non-negative and less than 1.
450 \end{methoddesc}
451 \begin{methoddesc}[StokesProblemCartesian]{getTolerance}{}
452 returns the current relative tolerance.
453 \end{methoddesc}
454 \begin{methoddesc}[StokesProblemCartesian]{setAbsoluteTolerance}{\optional{tolerance=0.}}
455 sets the absolute tolerance for the error in the relevant norm. The tolerance must be non-negative. Typically the
456 absolute talerance is set to 0.
457 \end{methoddesc}
458
459 \begin{methoddesc}[StokesProblemCartesian]{getAbsoluteTolerance}{}
460 sreturns the current absolute tolerance.
461 \end{methoddesc}
462
463 \begin{methoddesc}[StokesProblemCartesian]{getSolverOptionsVelocity}{}
464 returns the solver options used solve the equations~(\ref{V CALC}) for velocity.
465 \end{methoddesc}
466
467 \begin{methoddesc}[StokesProblemCartesian]{getSolverOptionsPressure}{}
468 returns the solver options used solve the equation~(\ref{P PREC}) for pressure.
469 \end{methoddesc}
470
471 \begin{methoddesc}[StokesProblemCartesian]{getSolverOptionsDiv}{}
472 set the solver options for solving the equation to project the divergence of the velocity onto the function space of pressure.
473 \end{methoddesc}
474
475
476 \subsection{Example: Lit Driven Cavity}
477 The following script \file{lit\hackscore driven\hackscore cavity.py}
478 \index{scripts!\file{helmholtz.py}} which is available in the \ExampleDirectory
479 illustrates the usage of the \class{StokesProblemCartesian} class to solve
480 the lit driven cavity problem:
481 \begin{python}
482 from esys.escript import *
483 from esys.finley import Rectangle
484 from esys.escript.models import StokesProblemCartesian
485 NE=25
486 dom = Rectangle(NE,NE,order=2)
487 x = dom.getX()
488 sc=StokesProblemCartesian(dom)
489 mask= (whereZero(x[0])*[1.,0]+whereZero(x[0]-1))*[1.,0] + \
490 (whereZero(x[1])*[0.,1.]+whereZero(x[1]-1))*[1.,1]
491 sc.initialize(eta=.1, fixed_u_mask= mask)
492 v=Vector(0.,Solution(dom))
493 v[0]+=whereZero(x[1]-1.)
494 p=Scalar(0.,ReducedSolution(dom))
495 v,p=sc.solve(v,p, verbose=True)
496 saveVTK("u.xml",velocity=v,pressure=p)
497 \end{python}

  ViewVC Help
Powered by ViewVC 1.1.26