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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Copyright dates update

1
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 % Copyright (c) 2003-2018 by The University of Queensland
4 % http://www.uq.edu.au
5 %
6 % Primary Business: Queensland, Australia
7 % Licensed under the Apache License, version 2.0
8 % http://www.apache.org/licenses/LICENSE-2.0
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=<eta derived from v and p>
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[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 \end{python}
455

  ViewVC Help
Powered by ViewVC 1.1.26