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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26