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

revision 2851 by gross, Fri Jan 15 07:37:47 2010 UTC revision 3295 by jfenwick, Fri Oct 22 01:56:02 2010 UTC
# Line 4  In this section we discuss how to solve Line 4  In this section we discuss how to solve
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}  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}  \begin{equation}\label{Stokes 1}
7  -\left(\eta(v\hackscore{i,j}+ v\hackscore{j,i})\right)\hackscore{,j}+p\hackscore{,i}=f\hackscore{i}-\sigma\hackscore{ij,j}  -\left(\eta(v_{i,j}+ v_{j,i})\right)_{,j}+p_{,i}=f_{i}-\sigma_{ij,j}
8  \end{equation}  \end{equation}
9  where  $f\hackscore{i}$ defines an internal force \index{force, internal} and $\sigma\hackscore{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.  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:  We assume an incompressible media:
12  \begin{equation}\label{Stokes 2}  \begin{equation}\label{Stokes 2}
13  -v\hackscore{i,i}=0  -v_{i,i}=0
14  \end{equation}  \end{equation}
15  Natural boundary conditions are taken in the form  Natural boundary conditions are taken in the form
16  \begin{equation}\label{Stokes Boundary}  \begin{equation}\label{Stokes Boundary}
17  \left(\eta(v\hackscore{i,j}+ v\hackscore{j,i})\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}  \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}  \end{equation}
19  which can be overwritten by constraints of the form  which can be overwritten by constraints of the form
20  \begin{equation}\label{Stokes Boundary0}  \begin{equation}\label{Stokes Boundary0}
21  v\hackscore{i}(x)=v^D\hackscore{i}(x)  v_{i}(x)=v^D_{i}(x)
22  \end{equation}  \end{equation}
23  at some locations $x$ at the boundary of the domain. $s\hackscore{i}$ defines a normal stress and  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.  $\alpha\ge 0$ the spring constant for restoring normal force.
25  The index $i$ may depend on the location $x$ on the boundary.  The index $i$ may depend on the location $x$ on the boundary.
26  $v^D$ is a given function on the domain.  $v^D$ is a given function on the domain.
# Line 46  G \\ Line 46  G \\
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).  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}.  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  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  velocity and pressure we calculate a correction $dv$ for the velocity by solving the first
51  equation of equation~\ref{STOKES}  equation of equation~\ref{STOKES}
52   \begin{equation}\label{STOKES ITER STEP 1}   \begin{equation}\label{STOKES ITER STEP 1}
53   A dv\hackscore{1} = G - A v\hackscore{0} - B^{*} p\hackscore{0}   A dv_{1} = G - A v_{0} - B^{*} p_{0}
54  \end{equation}  \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}$  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\hackscore{2}$ for the velocity by solving  for the pressure and an additional correction $dv_{2}$ for the velocity by solving
57   \begin{equation}   \begin{equation}
58   \begin{array}{rcl}   \begin{array}{rcl}
59   B A^{-1} B^{*} dp\hackscore{2} & = & Bv\hackscore{1} \\   B A^{-1} B^{*} dp_{2} & = & Bv_{1} \\
60   A dv\hackscore{2} & = & B^{*} dp\hackscore{2}   A dv_{2} & = & B^{*} dp_{2}
61  \end{array}  \end{array}
62   \label{STOKES ITER STEP 2}   \label{STOKES ITER STEP 2}
63   \end{equation}   \end{equation}
64  The new velocity and pressure are then given by $v\hackscore{2}=v\hackscore{1}-dv\hackscore{2}$ and  The new velocity and pressure are then given by $v_{2}=v_{1}-dv_{2}$ and
65  $p\hackscore{2}=p\hackscore{0}+dp\hackscore{2}$ which will fulfill the block system~\ref{STOKES}.  $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}.  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  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  to solve any problem for operator $A$. So we will be unable to calculate the operator
70  $B A^{-1} B^{*}$ required for $dp\hackscore{2}$ explicitly. In fact, we need to use another  $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  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  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  viscosity $\eta$ may depend on velocity or pressure and so we need to iterate over the
# Line 75  three equations~\ref{STOKES ITER STEP 1} Line 75  three equations~\ref{STOKES ITER STEP 1}
75
76  In the following we will use the two norms  In the following we will use the two norms
77  \begin{equation}  \begin{equation}
78  \|v\|\hackscore{1}^2 = \int\hackscore{\Omega} v\hackscore{j,k}v\hackscore{j,k} \; dx  \|v\|_{1}^2 = \int_{\Omega} v_{j,k}v_{j,k} \; dx
79  \mbox{ and }  \mbox{ and }
80  \|p\|\hackscore{0}^2= \int\hackscore{\Omega} p^2 \; dx.  \|p\|_{0}^2= \int_{\Omega} p^2 \; dx.
81  \label{STOKES STOP}  \label{STOKES STOP}
82  \end{equation}  \end{equation}
83  for velocity $v$ and pressure $p$. The iteration is terminated if the stopping criteria  for velocity $v$ and pressure $p$. The iteration is terminated if the stopping criteria
84   \begin{equation} \label{STOKES STOPPING CRITERIA}   \begin{equation} \label{STOKES STOPPING CRITERIA}
85  \max(\|Bv\hackscore{1}\|\hackscore{0},\|v\hackscore{2}-v\hackscore{0}\|\hackscore{1}) \le \tau \cdot \|v\hackscore{2}\|\hackscore{1}  \max(\|Bv_{1}\|_{0},\|v_{2}-v_{0}\|_{1}) \le \tau \cdot \|v_{2}\|_{1}
86   \end{equation}   \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   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\hackscore{1}\|\hackscore{0}$ equals the  $\|Bv_{1}\|_{0}$ equals the
89  norm of $B A^{-1} B^{*} dp\hackscore{2}$ and consequently provides a norm for the pressure correction.  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}  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  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.  we consider the errors produced by the iterative solvers being used.
94  From Equation~\ref{STOKES ITER STEP 1} we have  From Equation~\ref{STOKES ITER STEP 1} we have
95  \begin{equation} \label{STOKES total V1}  \begin{equation} \label{STOKES total V1}
96  v\hackscore{1} = e\hackscore{1} + v\hackscore{0} + A^{-1} ( G - Av\hackscore{0} - B^{*} p\hackscore{0} )  v_{1} = e_{1} + v_{0} + A^{-1} ( G - Av_{0} - B^{*} p_{0} )
97  \end{equation}  \end{equation}
98  where $e\hackscore{1}$ is the error when solving~\ref{STOKES ITER STEP 1}.    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 $\|.\|\hackscore{s}$ used in the stopping criteria for this equation. In fact we will have a stopping criteria of the form  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}  \begin{equation}
101  \| A e\hackscore{1} \|\hackscore{s}  = \| 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}  \| 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}  \end{equation}
103  where $\tau\hackscore{1}$ is the tolerance which we need to choose. This translates into the condition  where $\tau_{1}$ is the tolerance which we need to choose. This translates into the condition
104  \begin{equation}  \begin{equation}
105  \| e\hackscore{1} \|\hackscore{1} \le K \tau\hackscore{1} \| dv\hackscore{1} - e\hackscore{1} \|\hackscore{1}  \| e_{1} \|_{1} \le K \tau_{1} \| dv_{1} - e_{1} \|_{1}
106  \end{equation}  \end{equation}
107  The constant $K$ represents some uncertainty combining a variety of unknown factors such as the  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.    norm being used and the condition number of the stiffness matrix.
109  From the first equation of~\ref{STOKES ITER STEP 2} we have  From the first equation of~\ref{STOKES ITER STEP 2} we have
110  \begin{equation}\label{STOKES total P2}  \begin{equation}\label{STOKES total P2}
111  p\hackscore{2} =  p\hackscore{0} + (B A^{-1} B^{*})^{-1} (e\hackscore{2} + Bv\hackscore{1} )  p_{2} =  p_{0} + (B A^{-1} B^{*})^{-1} (e_{2} + Bv_{1} )
112  \end{equation}  \end{equation}
113  where $e\hackscore{2}$ represents the error when solving~\ref{STOKES ITER STEP 2}.  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  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 $\|.\|\hackscore{0}$ norm. As suitable preconditioner \index{preconditioner} for the iteration  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  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  the evaluation of the preconditioner $P$ for a given pressure increment $q$ is the solution of
118  \begin{equation} \label{STOKES P PREC}  \begin{equation} \label{STOKES P PREC}
# Line 124  the problem Line 124  the problem
124  A w = B^{*} s  A w = B^{*} s
125  \end{equation}  \end{equation}
126  with sufficient accuracy to return $q=Bw$. We assume that the desired tolerance is  with sufficient accuracy to return $q=Bw$. We assume that the desired tolerance is
127  sufficiently small, for instance one can take $\tau\hackscore{2}^2$  sufficiently small, for instance one can take $\tau_{2}^2$
128  where $\tau\hackscore{2}$ is the tolerance for~\ref{STOKES ITER STEP 2}.  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  In an implementation we use the fact that the residual $r$ is given as
131  \begin{equation} \label{STOKES RES }  \begin{equation} \label{STOKES RES }
132   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}   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}  \end{equation}
134  In particular we have $e\hackscore{2} = B v\hackscore{2}$  In particular we have $e_{2} = B v_{2}$
135  So the residual $r$ is represented by the updated velocity $v\hackscore{2}=v\hackscore{1}-dv\hackscore{2}$. In practice, if  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  one uses the velocity to represent the residual $r$ there is no need
137  to recover $dv\hackscore{2}$ in~\ref{STOKES ITER STEP 2} after $dp\hackscore{2}$ has been calculated.  to recover $dv_{2}$ in~\ref{STOKES ITER STEP 2} after $dp_{2}$ has been calculated.
138  In PCG the iteration is terminated if  In PCG the iteration is terminated if
139  \begin{equation} \label{STOKES P OPERATOR ERROR}  \begin{equation} \label{STOKES P OPERATOR ERROR}
140  \| P^{\frac{1}{2}}B v\hackscore{2} \|\hackscore{0} \le \tau\hackscore{2} \| P^{\frac{1}{2}}B v\hackscore{1} \|\hackscore{0}  \| P^{\frac{1}{2}}B v_{2} \|_{0} \le \tau_{2} \| P^{\frac{1}{2}}B v_{1} \|_{0}
141  \end{equation}  \end{equation}
142  where $\tau\hackscore{2}$ is the given tolerance. This translates into  where $\tau_{2}$ is the given tolerance. This translates into
143  \begin{equation} \label{STOKES P OPERATOR ERROR 2}  \begin{equation} \label{STOKES P OPERATOR ERROR 2}
144  \|e\hackscore{2}\|\hackscore{0} = \| B v\hackscore{2} \|\hackscore{0} \le M \tau\hackscore{2} \| B v\hackscore{1} \|\hackscore{0}  \|e_{2}\|_{0} = \| B v_{2} \|_{0} \le M \tau_{2} \| B v_{1} \|_{0}
145  \end{equation}  \end{equation}
146  where $M$ is taking care of the fact that $P^{\frac{1}{2}}$ is dropped.    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  As we assume that there is no significant error from solving with the operator $A$ we have
149  \begin{equation} \label{STOKES total V2}  \begin{equation} \label{STOKES total V2}
150  v\hackscore{2} =  v\hackscore{1} - dv\hackscore{2}  v_{2} =  v_{1} - dv_{2}
151  = v\hackscore{1}  - A^{-1} B^{*}dp  = v_{1}  - A^{-1} B^{*}dp
152  \end{equation}  \end{equation}
153  Combining the equations~\ref{STOKES total V1},~\ref{STOKES total P2} and~\ref{STOKES total V2} and  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  setting the errors to zero we can write the solution process as a fix point problem
# Line 157  v = \Phi(v,p) \mbox{ and } p = \Psi(u,p) Line 157  v = \Phi(v,p) \mbox{ and } p = \Psi(u,p)
157  \end{equation}  \end{equation}
158  with suitable functions $\Phi(v,p)$ and $\Psi(v,p)$ representing the iteration operator without  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  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\hackscore{2}= \delta p + \Psi(v\hackscore{0},p\hackscore{0})$ and  the update step in the form $p_{2}= \delta p + \Psi(v_{0},p_{0})$ and
161  $v\hackscore{2}= \delta v + \Phi(v\hackscore{0},p\hackscore{0})$ where  $v_{2}= \delta v + \Phi(v_{0},p_{0})$ where
162  the total error $\delta p$ and $\delta v$ are given as  the total error $\delta p$ and $\delta v$ are given as
163  \begin{equation}  \begin{equation}
164   \begin{array}{rcl}   \begin{array}{rcl}
165  \delta p & = &  (B A^{-1} B^{*})^{-1} ( e\hackscore{2} + B e\hackscore{1} ) \\  \delta p & = &  (B A^{-1} B^{*})^{-1} ( e_{2} + B e_{1} ) \\
166  \delta v & = &  e\hackscore{1} -  A^{-1} B^{*}\delta p  \;.  \delta v & = &  e_{1} -  A^{-1} B^{*}\delta p  \;.
167  \end{array}\label{STOKES ERRORS}  \end{array}\label{STOKES ERRORS}
168  \end{equation}  \end{equation}
169  Notice that $B\delta v = - e\hackscore{2}=-Bv\hackscore{2}$. Our task is now to choose the tolerances  Notice that $B\delta v = - e_{2}=-Bv_{2}$. Our task is now to choose the tolerances
170  $\tau\hackscore{1}$ and $\tau\hackscore{2}$ such that the global errors $\delta p$ and $\delta v$  $\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.  do not stop the convergence of the iteration process.
172
173  To measure convergence we use  To measure convergence we use
174  \begin{equation}  \begin{equation}
175  \epsilon = \max(\|v\hackscore{2}-v\|\hackscore{1}, \|B A^{-1} B^{*} (p\hackscore{2}-p)\|\hackscore{0})  \epsilon = \max(\|v_{2}-v\|_{1}, \|B A^{-1} B^{*} (p_{2}-p)\|_{0})
176  \end{equation}  \end{equation}
177  In practice using the fact that $B A^{-1} B^{*} (p\hackscore{2}-p\hackscore{0}) = B v\hackscore{1}$  In practice using the fact that $B A^{-1} B^{*} (p_{2}-p_{0}) = B v_{1}$
178  and assuming that $v\hackscore{2}$ gives a better approximation to the true $v$ than  and assuming that $v_{2}$ gives a better approximation to the true $v$ than
179  $v\hackscore{0}$ we will use the estimate  $v_{0}$ we will use the estimate
180  \begin{equation}  \begin{equation}
181  \epsilon = \max(\|v\hackscore{2}-v\hackscore{0}\|\hackscore{1}, \|B v\hackscore{1}\|\hackscore{0})  \epsilon = \max(\|v_{2}-v_{0}\|_{1}, \|B v_{1}\|_{0})
182  \end{equation}  \end{equation}
183  to estimate the progress of the iteration step after the step is completed.  to estimate the progress of the iteration step after the step is completed.
184  Note that the estimate of $\epsilon$    Note that the estimate of $\epsilon$
185  used in the stopping criteria~\ref{STOKES STOPPING CRITERIA}. If $\chi^{-}$ is the convergence rate assuming  used in the stopping criteria~\ref{STOKES STOPPING CRITERIA}. If $\chi^{-}$ is the convergence rate assuming
186  exact calculations, i.e. $e\hackscore{1}=0$ and $e\hackscore{2}=0$, we are expecting  exact calculations, i.e. $e_{1}=0$ and $e_{2}=0$, we are expecting
187  to maintain $\epsilon \le \chi^{-} \cdot \epsilon^{-}$. For the  to maintain $\epsilon \le \chi^{-} \cdot \epsilon^{-}$. For the
188  pressure increment we get:  pressure increment we get:
189  \begin{equation} \label{STOKES EST 1}  \begin{equation} \label{STOKES EST 1}
190  \begin{array}{rcl}  \begin{array}{rcl}
191  \|B A^{-1} B^{*} (p\hackscore{2}-p)\|\hackscore{0}  \|B A^{-1} B^{*} (p_{2}-p)\|_{0}
192   & \le & \|B A^{-1} B^{*} (p\hackscore{2}-\delta p-p)\|\hackscore{0} +   & \le & \|B A^{-1} B^{*} (p_{2}-\delta p-p)\|_{0} +
193  \|B A^{-1} B^{*} \delta p\|\hackscore{0} \\  \|B A^{-1} B^{*} \delta p\|_{0} \\
194   & = & \chi^{-} \cdot \epsilon^{-} + \|e\hackscore{2} + B e\hackscore{1}\|\hackscore{0}  \\   & = & \chi^{-} \cdot \epsilon^{-} + \|e_{2} + B e_{1}\|_{0}  \\
195  & \approx & \chi^{-} \cdot \epsilon^{-} + \|e\hackscore{2}\|\hackscore{0} \\  & \approx & \chi^{-} \cdot \epsilon^{-} + \|e_{2}\|_{0} \\
196  & \le & \chi^{-} \cdot \epsilon^{-} + M \tau\hackscore{2} \|B v\hackscore{1}\|\hackscore{0} \\    & \le & \chi^{-} \cdot \epsilon^{-} + M \tau_{2} \|B v_{1}\|_{0} \\
197  \end{array}  \end{array}
198  \end{equation}  \end{equation}
199  So we choose the value for $\tau\hackscore{2}$ from  So we choose the value for $\tau_{2}$ from
200  \begin{equation} \label{STOKES TOL2}  \begin{equation} \label{STOKES TOL2}
201   M \tau\hackscore{2} \|B v\hackscore{1}\|\hackscore{0}  \le (\chi^{-})^2 \epsilon^{-}   M \tau_{2} \|B v_{1}\|_{0}  \le (\chi^{-})^2 \epsilon^{-}
202  \end{equation}  \end{equation}
203  in order to make the perturbation for the termination of the pressure iteration a second order effect. We use a  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:  similar argument for the velocity:
205  \begin{equation}\label{STOKES EST 2}  \begin{equation}\label{STOKES EST 2}
206  \begin{array}{rcl}  \begin{array}{rcl}
207  \|v\hackscore{2}-v\|\hackscore{1} & \le & \|v\hackscore{2}-\delta v-v\|\hackscore{1} + \| \delta v\|\hackscore{1} \\  \|v_{2}-v\|_{1} & \le & \|v_{2}-\delta v-v\|_{1} + \| \delta v\|_{1} \\
208  & \le &  \chi^{-} \cdot \epsilon^{-}  + \| e\hackscore{1} -  A^{-1} B^{*}\delta p \|\hackscore{1} \\  & \le &  \chi^{-} \cdot \epsilon^{-}  + \| e_{1} -  A^{-1} B^{*}\delta p \|_{1} \\
209  & \approx &  \chi^{-} \cdot \epsilon^{-}  + \| e\hackscore{1} \|\hackscore{1} \\  & \approx &  \chi^{-} \cdot \epsilon^{-}  + \| e_{1} \|_{1} \\
210  & \le &  \chi^{-} \cdot \epsilon^{-}  +  K \tau\hackscore{1} \| dv\hackscore{1} - e\hackscore{1} \|\hackscore{1}  & \le &  \chi^{-} \cdot \epsilon^{-}  +  K \tau_{1} \| dv_{1} - e_{1} \|_{1}
211  \\  \\
212  & \le &  ( 1  + K \tau\hackscore{1}) \chi^{-} \cdot \epsilon^{-}  & \le &  ( 1  + K \tau_{1}) \chi^{-} \cdot \epsilon^{-}
213  \end{array}  \end{array}
214  \end{equation}  \end{equation}
215  So we choose the value for $\tau\hackscore{1}$ from  So we choose the value for $\tau_{1}$ from
216  \begin{equation} \label{STOKES TOL1}  \begin{equation} \label{STOKES TOL1}
217  K \tau\hackscore{1} \le \chi^{-}  K \tau_{1} \le \chi^{-}
218  \end{equation}  \end{equation}
219  Assuming we have estimates for $M$ and $K$ (if none is available, we use the value $1$.)  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  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^{-}}$.  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\hackscore{max}\le 1$.  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  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}  \begin{equation}\label{STOKES EST 3}
225  \chi \le ( 1 + \max(M \frac{\tau\hackscore{2} \|B v\hackscore{1}\|\hackscore{0}}{\chi^{-} \epsilon^{-}},K \tau\hackscore{1}  ) ) \cdot \chi^{-}  \chi \le ( 1 + \max(M \frac{\tau_{2} \|B v_{1}\|_{0}}{\chi^{-} \epsilon^{-}},K \tau_{1}  ) ) \cdot \chi^{-}
226  \end{equation}  \end{equation}
227  If we assume that this inequality would be an equation if we would have chosen the right values  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  $M^{+}$ and $K^{+}$ then we get
# Line 243  The same identity is used for to update Line 243  The same identity is used for to update
243  are then use in the next iteration step to control the tolerances.  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  In some cases one can observe that there is a significant change
246  in the velocity but the new velocity $v\hackscore{1}$ has still a  in the velocity but the new velocity $v_{1}$ has still a
247  small divergence, i.e. we have  small divergence, i.e. we have
248  $\|Bv\hackscore{1}\|\hackscore{0} \ll \|v\hackscore{1}-v\hackscore{0}\|\hackscore{1}$.  $\|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  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  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  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  significant relative to its divergence. In practice we will ignore the second half of the iteration step
253  as long as  as long as
254   \begin{equation}\label{STOKES LARGE BV1}   \begin{equation}\label{STOKES LARGE BV1}
255  \|Bv\hackscore{1}\|\hackscore{0} \le \theta \cdot \|v\hackscore{1}-v\hackscore{0}\|  \|Bv_{1}\|_{0} \le \theta \cdot \|v_{1}-v_{0}\|
256  \end{equation}  \end{equation}
257  where $0<\theta<1$ is a given factor. In this case we will also check the stopping criteria  where $0<\theta<1$ is a given factor. In this case we will also check the stopping criteria
258  with $v\hackscore{1}\rightarrow v\hackscore{2}$ but we will not correct $M$ in this case.  with $v_{1}\rightarrow v_{2}$ but we will not correct $M$ in this case.
259
260  Starting from an initial guess $v\hackscore{0}$ and $p\hackscore{0}$ for velocity and pressure  Starting from an initial guess $v_{0}$ and $p_{0}$ for velocity and pressure
261  the solution procedure is implemented as follows.  the solution procedure is implemented as follows.
262  \begin{enumerate}  \begin{enumerate}
263   \item calculate viscosity $\eta(v\hackscore{0},p)\hackscore{0}$ and assemble operator $A$ from $\eta$.   \item calculate viscosity $\eta(v_{0},p)_{0}$ and assemble operator $A$ from $\eta$.
264   \item calculate the tolerance $\tau\hackscore{1}$ from~\ref{STOKES TOL1}.   \item calculate the tolerance $\tau_{1}$ from~\ref{STOKES TOL1}.
265   \item Solve~\ref{STOKES ITER STEP 1} for $dv\hackscore{1}$ with tolerance $\tau\hackscore{1}$.   \item Solve~\ref{STOKES ITER STEP 1} for $dv_{1}$ with tolerance $\tau_{1}$.
266   \item update $v\hackscore{1}= v\hackscore{0}+ dv\hackscore{1}$   \item update $v_{1}= v_{0}+ dv_{1}$
267   \item if $Bv\hackscore{1}$ is large (see~\ref{STOKES LARGE BV1}):   \item if $Bv_{1}$ is large (see~\ref{STOKES LARGE BV1}):
268   \begin{enumerate}   \begin{enumerate}
269   \item calculate the tolerance $\tau\hackscore{2}$ from~\ref{STOKES TOL2}.   \item calculate the tolerance $\tau_{2}$ from~\ref{STOKES TOL2}.
270   \item Solve~\ref{STOKES ITER STEP 2} for $dp\hackscore{2}$ and $v\hackscore{2}$ with tolerance $\tau\hackscore{2}$.   \item Solve~\ref{STOKES ITER STEP 2} for $dp_{2}$ and $v_{2}$ with tolerance $\tau_{2}$.
271   \item update $p\hackscore{2}\leftarrow p\hackscore{0}+ dp\hackscore{2}$   \item update $p_{2}\leftarrow p_{0}+ dp_{2}$
272   \end{enumerate}   \end{enumerate}
273   \item else:   \item else:
274    \begin{enumerate}    \begin{enumerate}
275    \item update $p\hackscore{2}\leftarrow p$ and $v\hackscore{2}\leftarrow v\hackscore{1}$    \item update $p_{2}\leftarrow p$ and $v_{2}\leftarrow v_{1}$
276     \end{enumerate}     \end{enumerate}
277     \item calculate convergence measure $\epsilon$ and convergence rate $\chi$     \item calculate convergence measure $\epsilon$ and convergence rate $\chi$
278  \item if stopping criteria~\ref{STOKES STOPPING CRITERIA} holds:  \item if stopping criteria~\ref{STOKES STOPPING CRITERIA} holds:
279   \begin{enumerate}   \begin{enumerate}
280   \item return $v\hackscore{2}$ and $p\hackscore{2}$   \item return $v_{2}$ and $p_{2}$
281   \end{enumerate}   \end{enumerate}
282   \item else:   \item else:
283   \begin{enumerate}   \begin{enumerate}
284
285       \item update $M$ and $K$       \item update $M$ and $K$
286       \item goto step 1 with $v\hackscore{0}\leftarrow v\hackscore{2}$ and $p\hackscore{0}\leftarrow p\hackscore{2}$.       \item goto step 1 with $v_{0}\leftarrow v_{2}$ and $p_{0}\leftarrow p_{2}$.
287  \end{enumerate}  \end{enumerate}
288  \end{enumerate}  \end{enumerate}
289
# Line 359  set the solver options for solving the e Line 359  set the solver options for solving the e
359
360
361  \subsection{Example: Lit Driven Cavity}  \subsection{Example: Lit Driven Cavity}
362   The following script \file{lit\hackscore driven\hackscore cavity.py}   The following script \file{lit_driven_cavity.py}
363  \index{scripts!\file{helmholtz.py}} which is available in the \ExampleDirectory  \index{scripts!\file{helmholtz.py}} which is available in the \ExampleDirectory
364  illustrates the usage of the \class{StokesProblemCartesian} class to solve  illustrates the usage of the \class{StokesProblemCartesian} class to solve
365  the lit driven cavity problem:  the lit driven cavity problem:

Legend:
 Removed from v.2851 changed lines Added in v.3295