1 |
|
|
2 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
3 |
|
% |
4 |
|
% Copyright (c) 2003-2010 by University of Queensland |
5 |
|
% Earth Systems Science Computational Center (ESSCC) |
6 |
|
% http://www.uq.edu.au/esscc |
7 |
|
% |
8 |
|
% Primary Business: Queensland, Australia |
9 |
|
% Licensed under the Open Software License version 3.0 |
10 |
|
% http://www.opensource.org/licenses/osl-3.0.php |
11 |
|
% |
12 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
13 |
|
|
14 |
\section{The Stokes Problem} |
\section{The Stokes Problem} |
15 |
\label{STOKES PROBLEM} |
\label{STOKES PROBLEM} |
16 |
In this section we discuss how to solve the Stokes problem which is defined as follows: |
In this section we discuss how to solve the Stokes problem. |
17 |
|
We want to calculate the velocity\index{velocity} field $v$ and pressure $p$ |
18 |
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} |
of an incompressible fluid\index{incompressible fluid}. |
19 |
|
They are given as the solution of the Stokes problem\index{Stokes problem} |
20 |
\begin{equation}\label{Stokes 1} |
\begin{equation}\label{Stokes 1} |
21 |
-\left(\eta(v_{i,j}+ v_{j,i})\right)_{,j}+p_{,i}=f_{i}-\sigma_{ij,j} |
-\left(\eta(v_{i,j}+ v_{j,i})\right)_{,j}+p_{,i}=f_{i}-\sigma_{ij,j} |
22 |
\end{equation} |
\end{equation} |
23 |
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. |
where $f_{i}$ defines an internal force\index{force, internal} and |
24 |
|
$\sigma_{ij}$ is an initial stress\index{stress, initial}. |
25 |
|
The viscosity $\eta$ may weakly depend on pressure and velocity. |
26 |
|
If relevant we will use the notation $\eta(v,p)$ to express this dependency. |
27 |
|
|
28 |
We assume an incompressible media: |
We assume an incompressible medium: |
29 |
\begin{equation}\label{Stokes 2} |
\begin{equation}\label{Stokes 2} |
30 |
-v_{i,i}=0 |
-v_{i,i}=0 |
31 |
\end{equation} |
\end{equation} |
37 |
\begin{equation}\label{Stokes Boundary0} |
\begin{equation}\label{Stokes Boundary0} |
38 |
v_{i}(x)=v^D_{i}(x) |
v_{i}(x)=v^D_{i}(x) |
39 |
\end{equation} |
\end{equation} |
40 |
at some locations $x$ at the boundary of the domain. $s_{i}$ defines a normal stress and |
at some locations $x$ at the boundary of the domain. |
41 |
$\alpha\ge 0$ the spring constant for restoring normal force. |
$s_{i}$ defines a normal stress and $\alpha\ge 0$ the spring constant for |
42 |
|
restoring normal force. |
43 |
The index $i$ may depend on the location $x$ on the boundary. |
The index $i$ may depend on the location $x$ on the boundary. |
44 |
$v^D$ is a given function on the domain. |
$v^D$ is a given function on the domain. |
45 |
|
|
46 |
\subsection{Solution Method \label{STOKES SOLVE}} |
\subsection{Solution Method \label{STOKES SOLVE}} |
47 |
If we assume that $\eta$ is independent from the velocity and pressure equations~\ref{Stokes 1} and~\ref{Stokes 2} |
If we assume that $\eta$ is independent from the velocity and pressure, |
48 |
can be written in the block form |
equations~\ref{Stokes 1} and~\ref{Stokes 2} can be written in the block form |
49 |
\begin{equation} |
\begin{equation} |
50 |
\left[ \begin{array}{cc} |
\left[ \begin{array}{cc} |
51 |
A & B^{*} \\ |
A & B^{*} \\ |
61 |
\end{array} \right] |
\end{array} \right] |
62 |
\label{STOKES} |
\label{STOKES} |
63 |
\end{equation} |
\end{equation} |
64 |
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 a coercive, self-adjoint linear operator in a suitable Hilbert |
65 |
For more details on the mathematics see references \cite{AAMIRBERKYAN2008,MBENZI2005}. |
space, $B$ is the $(-1) \cdot$ divergence operator and $B^{*}$ is the adjoint |
66 |
|
operator (=gradient operator). |
67 |
If $v_{0}$ and $p_{0}$ are given initial guesses for |
For more details on the mathematics see references \cite{AAMIRBERKYAN2008,MBENZI2005}. |
68 |
velocity and pressure we calculate a correction $dv$ for the velocity by solving the first |
|
69 |
equation of equation~\ref{STOKES} |
If $v_{0}$ and $p_{0}$ are given initial guesses for velocity and pressure we |
70 |
|
calculate a correction $dv$ for the velocity by solving the first equation of |
71 |
|
\eqn{STOKES} |
72 |
\begin{equation}\label{STOKES ITER STEP 1} |
\begin{equation}\label{STOKES ITER STEP 1} |
73 |
A dv_{1} = G - A v_{0} - B^{*} p_{0} |
A dv_{1} = G - A v_{0} - B^{*} p_{0} |
74 |
\end{equation} |
\end{equation} |
75 |
We then insert the new approximation $v_{1}=v_{0}+dv_{1}$ to calculate a correction $dp_{2}$ |
We then insert the new approximation $v_{1}=v_{0}+dv_{1}$ to calculate a |
76 |
for the pressure and an additional correction $dv_{2}$ for the velocity by solving |
correction $dp_{2}$ for the pressure and an additional correction $dv_{2}$ for |
77 |
|
the velocity by solving |
78 |
\begin{equation} |
\begin{equation} |
79 |
\begin{array}{rcl} |
\begin{array}{rcl} |
80 |
B A^{-1} B^{*} dp_{2} & = & Bv_{1} \\ |
B A^{-1} B^{*} dp_{2} & = & Bv_{1} \\ |
84 |
\end{equation} |
\end{equation} |
85 |
The new velocity and pressure are then given by $v_{2}=v_{1}-dv_{2}$ and |
The new velocity and pressure are then given by $v_{2}=v_{1}-dv_{2}$ and |
86 |
$p_{2}=p_{0}+dp_{2}$ which will fulfill the block system~\ref{STOKES}. |
$p_{2}=p_{0}+dp_{2}$ which will fulfill the block system~\ref{STOKES}. |
87 |
This solution strategy is called the Uzawa scheme \index{Uzawa scheme}. |
This solution strategy is called the Uzawa scheme\index{Uzawa scheme}. |
88 |
|
|
89 |
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 |
90 |
to solve any problem for operator $A$. So we will be unable to calculate the operator |
scheme to solve any problem for operator $A$. |
91 |
$ B A^{-1} B^{*}$ required for $dp_{2}$ explicitly. In fact, we need to use another |
So we will be unable to calculate the operator $ B A^{-1} B^{*}$ required for |
92 |
iterative scheme to solve the first equation in~\ref{STOKES ITER STEP 2} where in each iteration step |
$dp_{2}$ explicitly. In fact, we need to use another iterative scheme to solve |
93 |
|
the first equation in~\ref{STOKES ITER STEP 2} where in each iteration step |
94 |
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 |
95 |
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 |
96 |
three equations~\ref{STOKES ITER STEP 1} and~\ref{STOKES ITER STEP 2}. |
over the three equations~\ref{STOKES ITER STEP 1} and~\ref{STOKES ITER STEP 2}. |
97 |
|
|
98 |
In the following we will use the two norms |
In the following we will use the two norms |
99 |
\begin{equation} |
\begin{equation} |
102 |
\|p\|_{0}^2= \int_{\Omega} p^2 \; dx. |
\|p\|_{0}^2= \int_{\Omega} p^2 \; dx. |
103 |
\label{STOKES STOP} |
\label{STOKES STOP} |
104 |
\end{equation} |
\end{equation} |
105 |
for velocity $v$ and pressure $p$. The iteration is terminated if the stopping criteria |
for velocity $v$ and pressure $p$. |
106 |
|
The iteration is terminated if the stopping criterion |
107 |
\begin{equation} \label{STOKES STOPPING CRITERIA} |
\begin{equation} \label{STOKES STOPPING CRITERIA} |
108 |
\max(\|Bv_{1}\|_{0},\|v_{2}-v_{0}\|_{1}) \le \tau \cdot \|v_{2}\|_{1} |
\max(\|Bv_{1}\|_{0},\|v_{2}-v_{0}\|_{1}) \le \tau \cdot \|v_{2}\|_{1} |
109 |
\end{equation} |
\end{equation} |
110 |
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 tolerance $0<\tau<1$ is met. |
111 |
$\|Bv_{1}\|_{0}$ equals the |
Notice that because of the first equation of~\ref{STOKES ITER STEP 2} we have |
112 |
norm of $B A^{-1} B^{*} dp_{2}$ and consequently provides a norm for the pressure correction. |
that $\|Bv_{1}\|_{0}$ equals the norm of $B A^{-1} B^{*} dp_{2}$ and |
113 |
|
consequently provides a norm for the pressure correction. |
114 |
|
|
115 |
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} |
116 |
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 |
117 |
we consider the errors produced by the iterative solvers being used. |
fixed point problem. Here we consider the errors produced by the iterative |
118 |
From Equation~\ref{STOKES ITER STEP 1} we have |
solvers being used. |
119 |
|
From \eqn{STOKES ITER STEP 1} we have |
120 |
\begin{equation} \label{STOKES total V1} |
\begin{equation} \label{STOKES total V1} |
121 |
v_{1} = e_{1} + v_{0} + A^{-1} ( G - Av_{0} - B^{*} p_{0} ) |
v_{1} = e_{1} + v_{0} + A^{-1} ( G - Av_{0} - B^{*} p_{0} ) |
122 |
\end{equation} |
\end{equation} |
123 |
where $ e_{1}$ is the error when solving~\ref{STOKES ITER STEP 1}. |
where $e_{1}$ is the error when solving~\ref{STOKES ITER STEP 1}. |
124 |
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 |
We will use a sparse matrix solver so we have not full control on the norm |
125 |
|
$\|.\|_{s}$ used in the stopping criterion for this equation. |
126 |
|
In fact we will have a stopping criterion of the form |
127 |
\begin{equation} |
\begin{equation} |
128 |
\| A e_{1} \|_{s} = \| G - A v_{1} - B^{*} p_{0} \|_{s} \le \tau_{1} \| G - A v_{0} - B^{*} p_{0} \|_{s} |
\| A e_{1} \|_{s} = \| G - A v_{1} - B^{*} p_{0} \|_{s} \le \tau_{1} \| G - A v_{0} - B^{*} p_{0} \|_{s} |
129 |
\end{equation} |
\end{equation} |
130 |
where $\tau_{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. |
131 |
|
This translates into the condition |
132 |
\begin{equation} |
\begin{equation} |
133 |
\| e_{1} \|_{1} \le K \tau_{1} \| dv_{1} - e_{1} \|_{1} |
\| e_{1} \|_{1} \le K \tau_{1} \| dv_{1} - e_{1} \|_{1} |
134 |
\end{equation} |
\end{equation} |
135 |
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 |
136 |
norm being used and the condition number of the stiffness matrix. |
factors such as the norm being used and the condition number of the stiffness matrix. |
137 |
From the first equation of~\ref{STOKES ITER STEP 2} we have |
From the first equation of~\ref{STOKES ITER STEP 2} we have |
138 |
\begin{equation}\label{STOKES total P2} |
\begin{equation}\label{STOKES total P2} |
139 |
p_{2} = p_{0} + (B A^{-1} B^{*})^{-1} (e_{2} + Bv_{1} ) |
p_{2} = p_{0} + (B A^{-1} B^{*})^{-1} (e_{2} + Bv_{1} ) |
140 |
\end{equation} |
\end{equation} |
141 |
where $e_{2}$ represents the error when solving~\ref{STOKES ITER STEP 2}. |
where $e_{2}$ represents the error when solving~\ref{STOKES ITER STEP 2}. |
142 |
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 |
143 |
operator $B A^{-1} B^{*}$ using the $\|.\|_{0}$ norm. As suitable preconditioner \index{preconditioner} for the iteration |
(PCG)\index{linear solver!PCG}\index{PCG} with iteration operator |
144 |
operator we use $\frac{1}{\eta}$ \cite{ELMAN}, ie |
$B A^{-1} B^{*}$ using the $\|.\|_{0}$ norm. |
145 |
the evaluation of the preconditioner $P$ for a given pressure increment $q$ is the solution of |
As suitable preconditioner\index{preconditioner} for the iteration operator we |
146 |
|
use $\frac{1}{\eta}$ \cite{ELMAN}, i.e. the evaluation of the preconditioner |
147 |
|
$P$ for a given pressure increment $q$ is the solution of |
148 |
\begin{equation} \label{STOKES P PREC} |
\begin{equation} \label{STOKES P PREC} |
149 |
\frac{1}{\eta} (Pq) = q \; . |
\frac{1}{\eta} (Pq) = q \; . |
150 |
\end{equation} |
\end{equation} |
151 |
Note that in each evaluation of the iteration operator $q=B A^{-1} B^{*} s$ one needs to solve |
Note that in each evaluation of the iteration operator $q=B A^{-1} B^{*} s$ |
152 |
the problem |
one needs to solve the problem |
153 |
\begin{equation} \label{STOKES P OPERATOR} |
\begin{equation} \label{STOKES P OPERATOR} |
154 |
A w = B^{*} s |
A w = B^{*} s |
155 |
\end{equation} |
\end{equation} |
156 |
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 |
157 |
sufficiently small, for instance one can take $\tau_{2}^2$ |
tolerance is sufficiently small, for instance one can take $\tau_{2}^2$ where |
158 |
where $\tau_{2}$ is the tolerance for~\ref{STOKES ITER STEP 2}. |
$\tau_{2}$ is the tolerance for~\ref{STOKES ITER STEP 2}. |
159 |
|
|
160 |
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 |
161 |
\begin{equation} \label{STOKES RES } |
\begin{equation} \label{STOKES RES } |
162 |
r= B (v_{1} - A^{-1} B^{*} dp) = B (v_{1} - A^{-1} B^{*} dp) = B (v_{1}-dv_{2}) = B v_{2} |
r= B (v_{1} - A^{-1} B^{*} dp) = B (v_{1} - A^{-1} B^{*} dp) = B (v_{1}-dv_{2}) = B v_{2} |
163 |
\end{equation} |
\end{equation} |
164 |
In particular we have $e_{2} = B v_{2}$ |
In particular we have $e_{2} = B v_{2}$. |
165 |
So the residual $r$ is represented by the updated velocity $v_{2}=v_{1}-dv_{2}$. In practice, if |
So the residual $r$ is represented by the updated velocity $v_{2}=v_{1}-dv_{2}$. |
166 |
one uses the velocity to represent the residual $r$ there is no need |
In practice, if one uses the velocity to represent the residual $r$ there is |
167 |
to recover $dv_{2}$ in~\ref{STOKES ITER STEP 2} after $dp_{2}$ has been calculated. |
no need to recover $dv_{2}$ in~\ref{STOKES ITER STEP 2} after $dp_{2}$ has |
168 |
|
been calculated. |
169 |
In PCG the iteration is terminated if |
In PCG the iteration is terminated if |
170 |
\begin{equation} \label{STOKES P OPERATOR ERROR} |
\begin{equation} \label{STOKES P OPERATOR ERROR} |
171 |
\| P^{\frac{1}{2}}B v_{2} \|_{0} \le \tau_{2} \| P^{\frac{1}{2}}B v_{1} \|_{0} |
\| P^{\frac{1}{2}}B v_{2} \|_{0} \le \tau_{2} \| P^{\frac{1}{2}}B v_{1} \|_{0} |
174 |
\begin{equation} \label{STOKES P OPERATOR ERROR 2} |
\begin{equation} \label{STOKES P OPERATOR ERROR 2} |
175 |
\|e_{2}\|_{0} = \| B v_{2} \|_{0} \le M \tau_{2} \| B v_{1} \|_{0} |
\|e_{2}\|_{0} = \| B v_{2} \|_{0} \le M \tau_{2} \| B v_{1} \|_{0} |
176 |
\end{equation} |
\end{equation} |
177 |
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. |
178 |
|
|
179 |
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 |
180 |
|
$A$ we have |
181 |
\begin{equation} \label{STOKES total V2} |
\begin{equation} \label{STOKES total V2} |
182 |
v_{2} = v_{1} - dv_{2} |
v_{2} = v_{1} - dv_{2} |
183 |
= v_{1} - A^{-1} B^{*}dp |
= v_{1} - A^{-1} B^{*}dp |
184 |
\end{equation} |
\end{equation} |
185 |
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} |
186 |
setting the errors to zero we can write the solution process as a fix point problem |
and setting the errors to zero we can write the solution process as a fix |
187 |
|
point problem |
188 |
\begin{equation} |
\begin{equation} |
189 |
v = \Phi(v,p) \mbox{ and } p = \Psi(u,p) |
v = \Phi(v,p) \mbox{ and } p = \Psi(u,p) |
190 |
\end{equation} |
\end{equation} |
191 |
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 |
192 |
errors. In fact for a linear problem, $\Phi$ and $\Psi$ are constant. With this notation we can write |
iteration operator without errors. In fact, for a linear problem, $\Phi$ and |
193 |
the update step in the form $p_{2}= \delta p + \Psi(v_{0},p_{0})$ and |
$\Psi$ are constant. With this notation we can write the update step in the |
194 |
$v_{2}= \delta v + \Phi(v_{0},p_{0})$ where |
form $p_{2}= \delta p + \Psi(v_{0},p_{0})$ and $v_{2}= \delta v + \Phi(v_{0},p_{0})$ |
195 |
the total error $\delta p$ and $\delta v$ are given as |
where the total error $\delta p$ and $\delta v$ are given as |
196 |
\begin{equation} |
\begin{equation} |
197 |
\begin{array}{rcl} |
\begin{array}{rcl} |
198 |
\delta p & = & (B A^{-1} B^{*})^{-1} ( e_{2} + B e_{1} ) \\ |
\delta p & = & (B A^{-1} B^{*})^{-1} ( e_{2} + B e_{1} ) \\ |
199 |
\delta v & = & e_{1} - A^{-1} B^{*}\delta p \;. |
\delta v & = & e_{1} - A^{-1} B^{*}\delta p \;. |
200 |
\end{array}\label{STOKES ERRORS} |
\end{array}\label{STOKES ERRORS} |
201 |
\end{equation} |
\end{equation} |
202 |
Notice that $B\delta v = - e_{2}=-Bv_{2}$. Our task is now to choose the tolerances |
Notice that $B\delta v = - e_{2}=-Bv_{2}$. |
203 |
$\tau_{1}$ and $\tau_{2}$ such that the global errors $\delta p$ and $\delta v$ |
Our task is now to choose the tolerances $\tau_{1}$ and $\tau_{2}$ such that |
204 |
do not stop the convergence of the iteration process. |
the global errors $\delta p$ and $\delta v$ do not stop the convergence of the |
205 |
|
iteration process. |
206 |
|
|
207 |
To measure convergence we use |
To measure convergence we use |
208 |
\begin{equation} |
\begin{equation} |
214 |
\begin{equation} |
\begin{equation} |
215 |
\epsilon = \max(\|v_{2}-v_{0}\|_{1}, \|B v_{1}\|_{0}) |
\epsilon = \max(\|v_{2}-v_{0}\|_{1}, \|B v_{1}\|_{0}) |
216 |
\end{equation} |
\end{equation} |
217 |
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. |
218 |
Note that the estimate of $\epsilon$ |
Note that the estimate of $\epsilon$ is used in the stopping |
219 |
used in the stopping criteria~\ref{STOKES STOPPING CRITERIA}. If $\chi^{-}$ is the convergence rate assuming |
criterion~\ref{STOKES STOPPING CRITERIA}. |
220 |
exact calculations, i.e. $e_{1}=0$ and $e_{2}=0$, we are expecting |
If $\chi^{-}$ is the convergence rate assuming exact calculations, i.e. |
221 |
to maintain $\epsilon \le \chi^{-} \cdot \epsilon^{-}$. For the |
$e_{1}=0$ and $e_{2}=0$, we are expecting to maintain $\epsilon \le \chi^{-} \cdot \epsilon^{-}$. |
222 |
pressure increment we get: |
For the pressure increment we get: |
223 |
\begin{equation} \label{STOKES EST 1} |
\begin{equation} \label{STOKES EST 1} |
224 |
\begin{array}{rcl} |
\begin{array}{rcl} |
225 |
\|B A^{-1} B^{*} (p_{2}-p)\|_{0} |
\|B A^{-1} B^{*} (p_{2}-p)\|_{0} |
230 |
& \le & \chi^{-} \cdot \epsilon^{-} + M \tau_{2} \|B v_{1}\|_{0} \\ |
& \le & \chi^{-} \cdot \epsilon^{-} + M \tau_{2} \|B v_{1}\|_{0} \\ |
231 |
\end{array} |
\end{array} |
232 |
\end{equation} |
\end{equation} |
233 |
So we choose the value for $\tau_{2}$ from |
So we choose the value for $\tau_{2}$ from |
234 |
\begin{equation} \label{STOKES TOL2} |
\begin{equation} \label{STOKES TOL2} |
235 |
M \tau_{2} \|B v_{1}\|_{0} \le (\chi^{-})^2 \epsilon^{-} |
M \tau_{2} \|B v_{1}\|_{0} \le (\chi^{-})^2 \epsilon^{-} |
236 |
\end{equation} |
\end{equation} |
237 |
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 |
238 |
similar argument for the velocity: |
iteration a second order effect. We use a similar argument for the velocity: |
239 |
\begin{equation}\label{STOKES EST 2} |
\begin{equation}\label{STOKES EST 2} |
240 |
\begin{array}{rcl} |
\begin{array}{rcl} |
241 |
\|v_{2}-v\|_{1} & \le & \|v_{2}-\delta v-v\|_{1} + \| \delta v\|_{1} \\ |
\|v_{2}-v\|_{1} & \le & \|v_{2}-\delta v-v\|_{1} + \| \delta v\|_{1} \\ |
250 |
\begin{equation} \label{STOKES TOL1} |
\begin{equation} \label{STOKES TOL1} |
251 |
K \tau_{1} \le \chi^{-} |
K \tau_{1} \le \chi^{-} |
252 |
\end{equation} |
\end{equation} |
253 |
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$\footnote{if no estimates are |
254 |
we can use~\ref{STOKES TOL1} and~\ref{STOKES TOL2} to get appropriate values for the tolerances. After |
available, we use the value $1$} we can use~\ref{STOKES TOL1} and |
255 |
the step has been completed we can calculate a new convergence rate $\chi =\frac{\epsilon}{\epsilon^{-}}$. |
\ref{STOKES TOL2} to get appropriate values for the tolerances. |
256 |
For partial reasons we restrict $\chi$ to be less or equal a given maximum value $\chi_{max}\le 1$. |
After the step has been completed we can calculate a new convergence rate |
257 |
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 |
$\chi =\frac{\epsilon}{\epsilon^{-}}$. |
258 |
|
For partial reasons we restrict $\chi$ to be less or equal a given maximum |
259 |
|
value $\chi_{max}\le 1$. |
260 |
|
If we see $\chi \le \chi^{-} (1+\chi^{-})$ our choices for the tolerances were |
261 |
|
suitable. Otherwise, we need to adjust the values for $K$ and $M$. |
262 |
|
From the estimates~\ref{STOKES EST 1} and~\ref{STOKES EST 2} we establish |
263 |
\begin{equation}\label{STOKES EST 3} |
\begin{equation}\label{STOKES EST 3} |
264 |
\chi \le ( 1 + \max(M \frac{\tau_{2} \|B v_{1}\|_{0}}{\chi^{-} \epsilon^{-}},K \tau_{1} ) ) \cdot \chi^{-} |
\chi \le ( 1 + \max(M \frac{\tau_{2} \|B v_{1}\|_{0}}{\chi^{-} \epsilon^{-}},K \tau_{1} ) ) \cdot \chi^{-} |
265 |
\end{equation} |
\end{equation} |
266 |
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 |
267 |
$M^{+}$ and $K^{+}$ then we get |
the right values $M^{+}$ and $K^{+}$ then we get |
268 |
\begin{equation}\label{STOKES EST 3b} |
\begin{equation}\label{STOKES EST 3b} |
269 |
\chi = ( 1 + \max(M^{+} \frac{\chi^{-}}{M},K^{+} \frac{\chi^{-}}{K}) ) \cdot \chi^{-} |
\chi = ( 1 + \max(M^{+} \frac{\chi^{-}}{M},K^{+} \frac{\chi^{-}}{K}) ) \cdot \chi^{-} |
270 |
\end{equation} |
\end{equation} |
271 |
From this equation we set for |
From this equation we see if our choice for $K$ was not good enough. |
272 |
than our choice for $K$ was not good enough. In this case we can calculate a new value |
In this case we can calculate a new value |
273 |
\begin{equation} |
\begin{equation} |
274 |
K^{+} = \frac{\chi-\chi^{-}}{(\chi^{-})^2} K |
K^{+} = \frac{\chi-\chi^{-}}{(\chi^{-})^2} K |
275 |
\end{equation} |
\end{equation} |
276 |
In practice we will use |
In practice we will use |
277 |
\begin{equation} |
\begin{equation} |
278 |
K^{+} = \max(\frac{\chi-\chi^{-}}{(\chi^{-})^2} K,\frac{1}{2}K,1) |
K^{+} = \max(\frac{\chi-\chi^{-}}{(\chi^{-})^2} K,\frac{1}{2}K,1) |
279 |
\end{equation} |
\end{equation} |
280 |
where the second term is used to reduce a potential overestimate of $K$. |
where the second term is used to reduce a potential overestimate of $K$. |
281 |
The same identity is used for to update $M$. The updated $M^{+}$ and $K^{+}$ |
The same identity is used for to update $M$. The updated $M^{+}$ and $K^{+}$ |
282 |
are then use in the next iteration step to control the tolerances. |
are then use in the next iteration step to control the tolerances. |
283 |
|
|
284 |
In some cases one can observe that there is a significant change |
In some cases one can observe that there is a significant change in the |
285 |
in the velocity but the new velocity $v_{1}$ has still a |
velocity but the new velocity $v_{1}$ has still a small divergence, i.e. |
286 |
small divergence, i.e. we have |
we have $\|Bv_{1}\|_{0} \ll \|v_{1}-v_{0}\|_{1}$. |
287 |
$\|Bv_{1}\|_{0} \ll \|v_{1}-v_{0}\|_{1}$. |
In this case we will get a small pressure increment and consequently only very |
288 |
In this case we will get a small pressure increment and consequently only very small changes to |
small changes to the velocity as a result of the second update step which |
289 |
the velocity as a result of the second update step which therefore can be skipped and |
therefore can be skipped and we can directly repeat the first update step |
290 |
we can directly repeat the first update step until the increment in velocity becomes |
until the increment in velocity becomes significant relative to its divergence. |
291 |
significant relative to its divergence. In practice we will ignore the second half of the iteration step |
In practice we will ignore the second half of the iteration step as long as |
|
as long as |
|
292 |
\begin{equation}\label{STOKES LARGE BV1} |
\begin{equation}\label{STOKES LARGE BV1} |
293 |
\|Bv_{1}\|_{0} \le \theta \cdot \|v_{1}-v_{0}\| |
\|Bv_{1}\|_{0} \le \theta \cdot \|v_{1}-v_{0}\| |
294 |
\end{equation} |
\end{equation} |
295 |
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 |
296 |
with $v_{1}\rightarrow v_{2}$ but we will not correct $M$ in this case. |
stopping criterion with $v_{1}\rightarrow v_{2}$ but we will not correct $M$ |
297 |
|
in this case. |
298 |
|
|
299 |
Starting from an initial guess $v_{0}$ and $p_{0}$ for velocity and pressure |
Starting from an initial guess $v_{0}$ and $p_{0}$ for velocity and pressure |
300 |
the solution procedure is implemented as follows. |
the solution procedure is implemented as follows: |
301 |
\begin{enumerate} |
\begin{enumerate} |
302 |
\item calculate viscosity $\eta(v_{0},p)_{0}$ and assemble operator $A$ from $\eta$. |
\item calculate viscosity $\eta(v_{0},p)_{0}$ and assemble operator $A$ from $\eta$ |
303 |
\item calculate the tolerance $\tau_{1}$ from~\ref{STOKES TOL1}. |
\item calculate the tolerance $\tau_{1}$ from \eqn{STOKES TOL1} |
304 |
\item Solve~\ref{STOKES ITER STEP 1} for $dv_{1}$ with tolerance $\tau_{1}$. |
\item solve \eqn{STOKES ITER STEP 1} for $dv_{1}$ with tolerance $\tau_{1}$ |
305 |
\item update $v_{1}= v_{0}+ dv_{1}$ |
\item update $v_{1}= v_{0}+ dv_{1}$ |
306 |
\item if $Bv_{1}$ is large (see~\ref{STOKES LARGE BV1}): |
\item if $Bv_{1}$ is large (see~\ref{STOKES LARGE BV1}): |
307 |
\begin{enumerate} |
\begin{enumerate} |
308 |
\item calculate the tolerance $\tau_{2}$ from~\ref{STOKES TOL2}. |
\item calculate the tolerance $\tau_{2}$ from~\ref{STOKES TOL2} |
309 |
\item Solve~\ref{STOKES ITER STEP 2} for $dp_{2}$ and $v_{2}$ with tolerance $\tau_{2}$. |
\item solve~\ref{STOKES ITER STEP 2} for $dp_{2}$ and $v_{2}$ with tolerance $\tau_{2}$ |
310 |
\item update $p_{2}\leftarrow p_{0}+ dp_{2}$ |
\item update $p_{2}\leftarrow p_{0}+ dp_{2}$ |
311 |
\end{enumerate} |
\end{enumerate} |
312 |
\item else: |
\item else: |
313 |
\begin{enumerate} |
\begin{itemize} |
314 |
\item update $p_{2}\leftarrow p$ and $v_{2}\leftarrow v_{1}$ |
\item update $p_{2}\leftarrow p$ and $v_{2}\leftarrow v_{1}$ |
315 |
\end{enumerate} |
\end{itemize} |
316 |
\item calculate convergence measure $\epsilon$ and convergence rate $\chi$ |
\item calculate convergence measure $\epsilon$ and convergence rate $\chi$ |
317 |
\item if stopping criteria~\ref{STOKES STOPPING CRITERIA} holds: |
\item if stopping criterion~\ref{STOKES STOPPING CRITERIA} holds: |
318 |
\begin{enumerate} |
\begin{itemize} |
319 |
\item return $v_{2}$ and $p_{2}$ |
\item return $v_{2}$ and $p_{2}$ |
320 |
\end{enumerate} |
\end{itemize} |
321 |
\item else: |
\item else: |
322 |
\begin{enumerate} |
\begin{enumerate} |
323 |
|
\item update $M$ and $K$ |
324 |
\item update $M$ and $K$ |
\item goto step 1 with $v_{0}\leftarrow v_{2}$ and $p_{0}\leftarrow p_{2}$. |
325 |
\item goto step 1 with $v_{0}\leftarrow v_{2}$ and $p_{0}\leftarrow p_{2}$. |
\end{enumerate} |
|
\end{enumerate} |
|
326 |
\end{enumerate} |
\end{enumerate} |
327 |
|
|
328 |
\subsection{Functions} |
\subsection{Functions} |
329 |
|
|
330 |
\begin{classdesc}{StokesProblemCartesian}{domain} |
\begin{classdesc}{StokesProblemCartesian}{domain} |
331 |
opens the Stokes problem\index{Stokes problem} on the \Domain domain. The domain |
opens the Stokes problem\index{Stokes problem} on the \Domain domain. |
332 |
needs to support LBB compliant elements for the Stokes problem, see~\cite{LBB} for details~\index{LBB condition}. |
The domain needs to support LBB compliant elements for the Stokes problem, see~\cite{LBB} for details\index{LBB condition}. |
333 |
For instance one can use second order polynomials for velocity and |
For instance one can use second order polynomials for velocity and first order |
334 |
first order polynomials for the pressure on the same element. Alternatively, one can use |
polynomials for the pressure on the same element. |
335 |
macro elements\index{macro elements} using linear polynomial for both pressure and velocity with a subdivided |
Alternatively, one can use macro elements\index{macro elements} using linear |
336 |
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 |
polynomials for both pressure and velocity with a subdivided element for the |
337 |
|
velocity. Typically, the macro element is more cost effective. |
338 |
|
The fact that pressure and velocity are represented in different ways is |
339 |
|
expressed by |
340 |
\begin{python} |
\begin{python} |
341 |
velocity=Vector(0.0, Solution(mesh)) |
velocity=Vector(0.0, Solution(mesh)) |
342 |
pressure=Scalar(0.0, ReducedSolution(mesh)) |
pressure=Scalar(0.0, ReducedSolution(mesh)) |
343 |
\end{python} |
\end{python} |
344 |
\end{classdesc} |
\end{classdesc} |
345 |
|
|
346 |
\begin{methoddesc}[StokesProblemCartesian]{initialize}{\optional{f=Data(), \optional{fixed_u_mask=Data(), \optional{eta=1, \optional{surface_stress=Data(), \optional{stress=Data()}, \optional{ |
\begin{methoddesc}[StokesProblemCartesian]{initialize}{\optional{f=Data(), |
347 |
restoration_factor=0}}}}}} |
\optional{fixed_u_mask=Data(), \optional{eta=1,% |
348 |
|
\optional{surface_stress=Data(), \optional{stress=Data()},% |
349 |
|
\optional{restoration_factor=0}}}}}} |
350 |
assigns values to the model parameters. In any call all values must be set. |
assigns values to the model parameters. In any call all values must be set. |
351 |
\var{f} defines the external force $f$, \var{eta} the viscosity $\eta$, |
\var{f} defines the external force $f$, \var{eta} the viscosity $\eta$, |
352 |
\var{surface_stress} the surface stress $s$ and \var{stress} the initial stress $\sigma$. |
\var{surface_stress} the surface stress $s$ and \var{stress} the initial stress $\sigma$. |
353 |
The locations and components where the velocity is fixed are set by |
The locations and components where the velocity is fixed are set by the values |
354 |
the values of \var{fixed_u_mask}. \var{restoration_factor} defines the restoring force factor $\alpha$. |
of \var{fixed_u_mask}. |
355 |
The method will try to cast the given values to appropriate |
\var{restoration_factor} defines the restoring force factor $\alpha$. |
356 |
\Data class objects. |
The method will try to cast the given values to appropriate \Data class objects. |
357 |
\end{methoddesc} |
\end{methoddesc} |
358 |
|
|
359 |
\begin{methoddesc}[StokesProblemCartesian]{solve}{v,p |
\begin{methoddesc}[StokesProblemCartesian]{solve}{v,p |
360 |
\optional{, max_iter=100 \optional{, verbose=False \optional{, usePCG=True }}}} |
\optional{, max_iter=100 \optional{, verbose=False \optional{, usePCG=True }}}} |
361 |
solves the problem and return approximations for velocity and pressure. |
solves the problem and returns approximations for velocity and pressure. |
362 |
The arguments \var{v} and \var{p} define initial guess. |
The arguments \var{v} and \var{p} define initial guesses. |
363 |
\var{v} must have function space \var{Solution(domain)} and |
\var{v} must have function space \var{Solution(domain)} and \var{p} must have |
364 |
\var{p} must have function space \var{ReducedSolution(domain)}. |
function space \var{ReducedSolution(domain)}. |
365 |
The values of \var{v} marked |
The values of \var{v} marked by \var{fixed_u_mask} remain unchanged. |
366 |
by \var{fixed_u_mask} remain unchanged. |
If \var{usePCG} is set to \True then the preconditioned conjugate gradient |
367 |
If \var{usePCG} is set to \True |
method (PCG)\index{preconditioned conjugate gradient method!PCG} scheme is |
368 |
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 |
used. Otherwise the problem is solved with the generalized minimal residual |
369 |
the PCG scheme is more efficient. |
method (GMRES)\index{generalized minimal residual method!GMRES}. |
370 |
\var{max_iter} defines the maximum number of iteration steps. |
In most cases the PCG scheme is more efficient. |
371 |
|
\var{max_iter} defines the maximum number of iteration steps. |
372 |
If \var{verbose} is set to \True informations on the progress of of the solver are printed. |
If \var{verbose} is set to \True information on the progress of of the solver |
373 |
|
is printed. |
374 |
\end{methoddesc} |
\end{methoddesc} |
375 |
|
|
|
|
|
376 |
\begin{methoddesc}[StokesProblemCartesian]{setTolerance}{\optional{tolerance=1.e-4}} |
\begin{methoddesc}[StokesProblemCartesian]{setTolerance}{\optional{tolerance=1.e-4}} |
377 |
sets the tolerance in an appropriate norm relative to the right hand side. The tolerance must be non-negative and less than 1. |
sets the tolerance in an appropriate norm relative to the right hand side. |
378 |
|
The tolerance must be non-negative and less than 1. |
379 |
\end{methoddesc} |
\end{methoddesc} |
380 |
|
|
381 |
\begin{methoddesc}[StokesProblemCartesian]{getTolerance}{} |
\begin{methoddesc}[StokesProblemCartesian]{getTolerance}{} |
382 |
returns the current relative tolerance. |
returns the current relative tolerance. |
383 |
\end{methoddesc} |
\end{methoddesc} |
384 |
|
|
385 |
\begin{methoddesc}[StokesProblemCartesian]{setAbsoluteTolerance}{\optional{tolerance=0.}} |
\begin{methoddesc}[StokesProblemCartesian]{setAbsoluteTolerance}{\optional{tolerance=0.}} |
386 |
sets the absolute tolerance for the error in the relevant norm. The tolerance must be non-negative. Typically the |
sets the absolute tolerance for the error in the relevant norm. |
387 |
absolute tolerance is set to 0. |
The tolerance must be non-negative. |
388 |
|
Typically the absolute tolerance is set to 0. |
389 |
\end{methoddesc} |
\end{methoddesc} |
390 |
|
|
391 |
\begin{methoddesc}[StokesProblemCartesian]{getAbsoluteTolerance}{} |
\begin{methoddesc}[StokesProblemCartesian]{getAbsoluteTolerance}{} |
393 |
\end{methoddesc} |
\end{methoddesc} |
394 |
|
|
395 |
\begin{methoddesc}[StokesProblemCartesian]{getSolverOptionsVelocity}{} |
\begin{methoddesc}[StokesProblemCartesian]{getSolverOptionsVelocity}{} |
396 |
returns the solver options used solve the equations~(\ref{STOKES ITER STEP 1}) and~(\ref{STOKES P OPERATOR}) for velocity. |
returns the solver options used to solve \eqn{STOKES ITER STEP 1} and \eqn{STOKES P OPERATOR}) for velocity. |
397 |
\end{methoddesc} |
\end{methoddesc} |
398 |
|
|
399 |
\begin{methoddesc}[StokesProblemCartesian]{getSolverOptionsPressure}{} |
\begin{methoddesc}[StokesProblemCartesian]{getSolverOptionsPressure}{} |
400 |
returns the solver options used solve the preconditioner equation~(\ref{STOKES P PREC}) for pressure. |
returns the solver options used to solve the preconditioner \eqn{STOKES P PREC} for pressure. |
401 |
\end{methoddesc} |
\end{methoddesc} |
402 |
|
|
403 |
\begin{methoddesc}[StokesProblemCartesian]{getSolverOptionsDiv}{} |
\begin{methoddesc}[StokesProblemCartesian]{getSolverOptionsDiv}{} |
404 |
set the solver options for solving the equation to project the divergence of the velocity onto the function space of pressure. |
sets the solver options for solving the equation to project the divergence of |
405 |
|
the velocity onto the function space of pressure. |
406 |
\end{methoddesc} |
\end{methoddesc} |
407 |
|
|
408 |
|
\subsection{Example: Lid-driven Cavity} |
409 |
\subsection{Example: Lit Driven Cavity} |
The following script \file{lid_driven_cavity.py}\index{scripts!\file{lid_driven_cavity.py}} |
410 |
The following script \file{lit_driven_cavity.py} |
which is available in the \ExampleDirectory illustrates the usage of the |
411 |
\index{scripts!\file{helmholtz.py}} which is available in the \ExampleDirectory |
\class{StokesProblemCartesian} class to solve the lid-driven cavity problem: |
|
illustrates the usage of the \class{StokesProblemCartesian} class to solve |
|
|
the lit driven cavity problem: |
|
412 |
\begin{python} |
\begin{python} |
413 |
from esys.escript import * |
from esys.escript import * |
414 |
from esys.finley import Rectangle |
from esys.finley import Rectangle |
415 |
from esys.escript.models import StokesProblemCartesian |
from esys.escript.models import StokesProblemCartesian |
416 |
NE=25 |
NE=25 |
417 |
dom = Rectangle(NE,NE,order=2) |
dom = Rectangle(NE,NE,order=2) |
418 |
x = dom.getX() |
x = dom.getX() |
419 |
sc=StokesProblemCartesian(dom) |
sc=StokesProblemCartesian(dom) |
420 |
mask= (whereZero(x[0])*[1.,0]+whereZero(x[0]-1))*[1.,0] + \ |
mask= (whereZero(x[0])*[1.,0]+whereZero(x[0]-1))*[1.,0] + \ |
421 |
(whereZero(x[1])*[0.,1.]+whereZero(x[1]-1))*[1.,1] |
(whereZero(x[1])*[0.,1.]+whereZero(x[1]-1))*[1.,1] |
422 |
sc.initialize(eta=.1, fixed_u_mask= mask) |
sc.initialize(eta=.1, fixed_u_mask=mask) |
423 |
v=Vector(0.,Solution(dom)) |
v=Vector(0., Solution(dom)) |
424 |
v[0]+=whereZero(x[1]-1.) |
v[0]+=whereZero(x[1]-1.) |
425 |
p=Scalar(0.,ReducedSolution(dom)) |
p=Scalar(0.,ReducedSolution(dom)) |
426 |
v,p=sc.solve(v,p, verbose=True) |
v,p=sc.solve(v, p, verbose=True) |
427 |
saveVTK("u.xml",velocity=v,pressure=p) |
saveVTK("u.vtu", velocity=v, pressure=p) |
428 |
\end{python} |
\end{python} |
429 |
|
|