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

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

Parent Directory | Revision Log | View Patch Patch

-revision 2793 by gross,
Tue Dec  1 06:10:10 2009 UTC
+revision 2843 by gross,
Thu Jan 14 05:51:28 2010 UTC
 Line 4 
 In this section we discuss how to solve
  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}
  \begin{equation}\label{Stokes 1}
- -\left(\eta(v\hackscore{i,j}+ v\hackscore{i,j})\right)\hackscore{,j}+p\hackscore{,i}=f\hackscore{i}-\sigma\hackscore{ij,j}
+ -\left(\eta(v\hackscore{i,j}+ v\hackscore{j,i})\right)\hackscore{,j}+p\hackscore{,i}=f\hackscore{i}-\sigma\hackscore{ij,j}
  \end{equation}
  where  $f\hackscore{i}$ defines an internal force \index{force, internal} and $\sigma\hackscore{ij}$ is an intial stress \index{stress, initial}. The viscosity $\eta$ may weakly depend on pressure and velocity. If relevant we will use the notation $\eta(v,p)$ to express this dependency.
 Line 14 
 We assume an incompressible media:
  \end{equation}
  Natural boundary conditions are taken in the form
  \begin{equation}\label{Stokes Boundary}
- \left(\eta(v\hackscore{i,j}+ v\hackscore{i,j})\right)n\hackscore{j}-n\hackscore{i}p=s\hackscore{i} - \alpha \cdot n\hackscore{i} n\hackscore{j} v\hackscore{j}+\sigma\hackscore{ij} n\hackscore{j}
+ \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}
  \end{equation}
  which can be overwritten by constraints of the form
  \begin{equation}\label{Stokes Boundary0}
 Line 41 
 p \\
  G \\
 \\
  \end{array} \right]
- \label{SADDLEPOINT}
+ \label{STOKES}
  \end{equation}
  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).
  For more details on the mathematics see references \cite{AAMIRBERKYAN2008,MBENZI2005}.
  If $v\hackscore{0}$ and $p\hackscore{0}$ are given initial guesses for
  velocity and pressure we calculate a correction $dv$ for the velocity by solving the first
- equation of equation~\ref{SADDLEPOINT}
+ equation of equation~\ref{STOKES}
-  \begin{equation}\label{SADDLEPOINT ITER STEP 1}
+  \begin{equation}\label{STOKES ITER STEP 1}
   A dv\hackscore{1} = G - A v\hackscore{0} - B^{*} p\hackscore{0}
  \end{equation}
  We then insert the new approximation $v\hackscore{1}=v\hackscore{0}+dv\hackscore{1}$ to calculate a correction $dp\hackscore{2}$
 Line 59 
 for the pressure and an additional corre
   B A^{-1} B^{*} dp\hackscore{2} & = & Bv\hackscore{1} \\
   A dv\hackscore{2} & = & B^{*} dp\hackscore{2}
  \end{array}
-  \label{SADDLEPOINT ITER STEP 2}
+  \label{STOKES ITER STEP 2}
   \end{equation}
- The new velocity and pressure are then given by $v=v\hackscore{1}-dv\hackscore{2}$ and
+ The new velocity and pressure are then given by $v\hackscore{2}=v\hackscore{1}-dv\hackscore{2}$ and
- $p=p\hackscore{0}+dp\hackscore{2}$ which will fullfill the block system~\ref{SADDLEPOINT}.
+ $p\hackscore{2}=p\hackscore{0}+dp\hackscore{2}$ which will fullfill the block system~\ref{STOKES}.
- This solution strategy is called the Uzawa scheme \index{Uzawa scheme}. There are
+ This solution strategy is called the Uzawa scheme \index{Uzawa scheme}.
- two problems with this scheme. Firstly, in practice we will use an iterative scheme
- to solve the problem for operator $A$. So we will be unable to calculate the operator
+ There is a problem with this scheme: In practice we will use an iterative scheme
- $ B A^{-1} B^{*}$ required for $dp\hackscore{2}$. In fact, we need to use an iterative scheme
+ to solve any problem for operator $A$. So we will be unable to calculate the operator
- to solve the first equation in~\ref{SADDLEPOINT ITER STEP 2} where in each iteration step
+ $ B A^{-1} B^{*}$ required for $dp\hackscore{2}$ explicitly. In fact, we need to use another
- an iterative solver for $A$ is applied. The second issue is that
+ iterative scheme to solve the first equation in~\ref{STOKES ITER STEP 2} where in each iteration step
- viscosity $\eta$ may depend of velocity or pressure and so we need to iterate over the
+ an iterative solver for $A$ is applied. Another issue is the fact that the
- three equations~\ref{SADDLEPOINT ITER STEP 1} and~\ref{SADDLEPOINT ITER STEP 2}. Using the
+ viscosity $\eta$ may depend on velocity or pressure and so we need to iterate over the
- two norms
+ three equations~\ref{STOKES ITER STEP 1} and~\ref{STOKES ITER STEP 2}.
+ In the following we will use the two norms
  \begin{equation}
  \|v\|\hackscore{1}^2 = \int\hackscore{\Omega} v\hackscore{j,k}v\hackscore{j,k} \; dx
  \mbox{ and }
  \|p\|\hackscore{0}^2= \int\hackscore{\Omega} p^2 \; dx.
  \label{STOKES STOP}
  \end{equation}
- to define the stopping criterium
+ for velicity $v$ and pressure $p$. The iteration is terminated if the stopping criterium
-  \begin{equation}
+  \begin{equation} \label{STOKES STOPPING CRITERIA}
- \max(\|Bv\hackscore{1}\|\hackscore{0},\|v-v\hackscore{0}\|\hackscore{1}) \le \tau \cdot \|v\|\hackscore{1}
+ \max(\|Bv\hackscore{1}\|\hackscore{0},\|v\hackscore{2}-v\hackscore{0}\|\hackscore{1}) \le \tau \cdot \|v\hackscore{2}\|\hackscore{1}
   \end{equation}
- to terminate the overall iteration with a overall given tolerance $0<\tau<1$.
+  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
- Notice that because of the first equation of
+ $\|Bv\hackscore{1}\|\hackscore{0}$ equals the
- and~\ref{SADDLEPOINT ITER STEP 2} $\|Bv\hackscore{1}\|\hackscore{0}$ is equal to the
+ 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\hackscore{2}$ and consequently provides a norm for the  pressure correction.
+ We want to optimize the tolerance choice for solving~\ref{STOKES ITER STEP 1}
- Let us first have a look at the solution process~\ref{SADDLEPOINT ITER STEP 1} which needs to resolve the
+ and~\ref{STOKES ITER STEP 2}. To do this we write the iteration scheme as a fixed point problem. Here
- non-linearity of the viscosity coeficient. We will run a iterative process by resolving
+ we consider the errors produced by the iterative solvers beeing used.
- the equation~\ref{SADDLEPOINT ITER STEP 1} with $v\hackscore{0} \leftarrow v\hackscore{1}$. The first observation is that there is no point start
+ From Equation~\ref{STOKES ITER STEP 1} we have
- the iteration process over the second set of equation if $\|Bv\hackscore{1}\|\hackscore{0}$
+ \begin{equation} \label{STOKES total V1}
- is smaller then $\|v\hackscore{1}-v\hackscore{0}\|\hackscore{1}$. So we will terminate the non-linear iteration
+ v\hackscore{1} = e\hackscore{1} + v\hackscore{0} + A^{-1} ( G - Av\hackscore{0} - B^{*} p\hackscore{0} )
- process of the first equation if
-  \begin{equation}
- \theta \cdot \|v\hackscore{1}-v\hackscore{0}\| \le \|Bv\hackscore{1}\|\hackscore{0}
- \mbox{ or }
- \|v\hackscore{1}-v\hackscore{0}\|\hackscore{1} \le \tau \cdot \|v\hackscore{1}\|\hackscore{1}
- \end{equation}
- where $0<\theta<1$ is a given factor (typically $\theta=0.5$).
- We need to think about appropriate stopping criteria when solving
- \ref{SADDLEPOINT ITER STEP 1}. We would like to use very weak convergence criteria to reduce computational costs but they need to be tight enough to not interfere with the
- convergence of the iteration process one level above.
- We will use a sparse matrix solver so we have not full control on the norm $\|.\|\hackscore{s}$ used in the stopping criteria
- \begin{equation}
- \| 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}
  \end{equation}
- If $e\hackscore{1}$ is the error of the returned solution $dv\hackscore{1}$ this translates into the condition
+ where $ e\hackscore{1}$ is the error when solving~\ref{STOKES ITER STEP 1}.
+ 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
  \begin{equation}
- \| e\hackscore{1} \|\hackscore{1} \le K \tau\hackscore{1} \| dv\hackscore{1} - e\hackscore{1} \|\hackscore{1}
+ \| 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}
  \end{equation}
- The constant $K$ represents some uncertainty combining a variety of unknown factors such as the
+ where $\tau\hackscore{1}$ is the tolerance which we need to choose. This translates into the condition
- solver being used and the condition number of the stiffness matrix.  This leads to the estimate
- \begin{equation}
- \| dv\hackscore{1} \|\hackscore{1} \le (1+K \tau\hackscore{1}) \chi \| dv\hackscore{1}^- \|\hackscore{1}
- \label{DV COND}
- \end{equation}
- where $chi$ is the convergence rate of the \ref{SADDLEPOINT ITER STEP 1} iteration and $dv\hackscore{1}^-$ the correction from the last step. In order to make sure that the error $e\hackscore{1}$ can be neglected we need
- to have $K \tau\hackscore{1} \le \chi$ which leads to the condition
-  \begin{equation}
-  \tau\hackscore{1} \le \frac{\min(\chi,\frac{1}{2})}{K}
-  \label{NEW TAU 1}
- \end{equation}
- We can estimate $K$ by comparing estimates for the convergence rate $\chi=\frac{\| dv\hackscore{1} \|\hackscore{1}}{\| dv\hackscore{1}^- \|\hackscore{1}}$ by checking our error prediction~\ref{DV COND}. If we have
- access to the estimate of the convergence rate $\chi^-$ of the previous iteration step we
- get
-  \begin{equation}
-  \frac{1}{K} = \frac{\chi^-}{\chi-\chi^{-}} \tau\hackscore{1}^-
- \end{equation}
- if $\chi \ge \chi^- (1 + \chi^-) $ which means that our prediction for a suitable tolerane based on the avaibale $K$ value was wrong.
- So we start with the value $K=1$ in~\ref{NEW TAU 1} and recalculate
- $K$ if we have underestimate $K$.
- If $dv\hackscore{1}$
- From Equation~\ref{SADDLEPOINT ITER STEP 1} we have
- \begin{equation}
- v\hackscore{1} = e\hackscore{1} + A^{-1} ( G - B^{*} p )
- \end{equation}
- This translates into the conditoon
  \begin{equation}
  \| e\hackscore{1} \|\hackscore{1} \le K \tau\hackscore{1} \| dv\hackscore{1} - e\hackscore{1} \|\hackscore{1}
- \mbox{ therefore }
- \| e\hackscore{1} \|\hackscore{1} \le \frac{K \tau\hackscore{1}}{1-K \tau\hackscore{1}}  \| dv\hackscore{1} \|\hackscore{1}
  \end{equation}
  The constant $K$ represents some uncertainty combining a variety of unknown factors such as the
- solver being used and the condition number of the stiffness matrix.
+ norm being used and the condition number of the stiffness matrix.
+ From the first equation of~\ref{STOKES ITER STEP 2} we have
+ \begin{equation}\label{STOKES total P2}
+ p\hackscore{2} =  p\hackscore{0} + (B A^{-1} B^{*})^{-1} (e\hackscore{2} + Bv\hackscore{1} )
+ \end{equation}
+ where $e\hackscore{2}$ represents the error when solving~\ref{STOKES ITER STEP 2}.
- ===================================
+ We use an iterative preconditioned conjugate gradient method (PCG) \index{linear solver!PCG}\index{PCG} with iteration
+ operator $B A^{-1} B^{*}$ using the $\|.\|\hackscore{0}$ norm. As suitable preconditioner \index{preconditioner} for the iteration
+ operator we use $\frac{1}{\eta}$ \cite{ELMAN}, ie
- We use iterative techniques to solve this problem: Given an approximation $v$ and $p$ for
- velocity and pressure we perform the following steps in the Uzawa scheme \index{Uzawa scheme} style:
- \begin{enumerate}
-  \item calculate viscosity $\eta(v,p)$ and assemble operator $A$ from $\eta$.
-  \item Solve for $dv$:
-  \begin{equation}
-  A dv = G - A v - B^{*} p  \label{SADDLEPOINT ITER STEP 1}
- \end{equation}
-  \item update $v\hackscore{1}= v+ dv$
-  \item if $\max(\|Bv\hackscore{1}\|\hackscore{0},\|dv\|\hackscore{1}) \le \tau \cdot \|v\hackscore{1}\|\hackscore{1}$: return v$\hackscore{1},p$
-  \item Solve for $dp$:
-  \begin{equation}
-  \begin{array}{rcl}
-  B A^{-1} B^{*} dp & = & Bv\hackscore{1} \\
-  A dv\hackscore{2} & = & B^{*} dp
- \end{array}
-  \label{SADDLEPOINT ITER STEP 2}
-  \end{equation}
-  \item update $p\hackscore{2}\leftarrow p+ dp$ and $v\hackscore{2}= v\hackscore{1} - dv\hackscore{2}$
-  \item goto Step 0 with $v=v\hackscore{2}$ and $p=p\hackscore{2}$.
- \end{enumerate}
- where $\tau$ is the given tolerance. Notice that the operator $A$ if it is depending on $v$ or $p$ is updated once only.
- In this algorithm
- \begin{equation}
- \|v\|\hackscore{1}^2 = \int\hackscore{\Omega} v\hackscore{j,k}v\hackscore{j,k} \; dx
- \mbox{ and }
- \|p\|\hackscore{0}^2= \int\hackscore{\Omega} p^2 \; dx.
- \label{STOKES STOP}
- \end{equation}
- so the stopping criterion used is checking for convergence as well as for
- the fact that the incompressiblity condition~\ref{Stokes 2} is fullfilled.
- To solve the update step~\ref{SADDLEPOINT ITER STEP 2} we use an iterative solver with iteration
- operator $B A^{-1} B^{*}$, eg. using generalized minimal residual method (GMRES) \index{linear solver!GMRES}\index{GMRES}, preconditioned conjugate gradient method (PCG) \index{linear solver!PCG}\index{PCG}. As suitable preconditioner \index{preconditioner} for this operator is $\frac{1}{\eta}$, ie
  the evaluation of the preconditioner $P$ for a given pressure increment $q$ is the solution of
- \begin{equation} \label{P PREC}
+ \begin{equation} \label{STOKES P PREC}
- \frac{1}{\eta} Pq = q \; .
+ \frac{1}{\eta} (Pq) = q \; .
  \end{equation}
  Note that in each evaluation of the iteration operator $q=B A^{-1} B^{*} s$ one needs to solve
  the problem
- \begin{equation} \label{P OPERATOR}
+ \begin{equation} \label{STOKES P OPERATOR}
  A w = B^{*} s
  \end{equation}
- with sufficient accuracy to return $q=Bw$. Notice that the residual $r$ is given as
+ with sufficient accuracy to return $q=Bw$. We assume that the choicen tolerance is
+ sufficiently small, for instane one can take $\tau\hackscore{2}^2$
+ where $\tau\hackscore{2}$ is the tolerance for~\ref{STOKES ITER STEP 2}.
+ In an implementation we use the fact that the residual $r$ is given as
  \begin{equation} \label{STOKES RES }
   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}
  \end{equation}
- so in fact the residual $r$ is represented by the updated velocity $v\hackscore{2}$. This saves the recovery of
+ In particular we have $e\hackscore{2} = B v\hackscore{2}$
- $dv\hackscore{2}$ in~\ref{SADDLEPOINT ITER STEP 2} after $dp$ has been calculated as iterative method such as PCG calculate the solution approximations along with the their residual. In PCG the iteration is terminated if
+ So the residual $r$ is represented by the updated velocity $v\hackscore{2}=v\hackscore{1}-dv\hackscore{2}$. In practice, if
- \begin{equation} \label{P OPERATOR}
+ one uses the velocity to represent the residual $r$ there is no need
+ to recover $dv\hackscore{2}$ in~\ref{STOKES ITER STEP 2} after $dp\hackscore{2}$ has been calculated.
+ In PCG the iteration is terminated if
+ \begin{equation} \label{STOKES P OPERATOR ERROR}
  \| P^{\frac{1}{2}}B v\hackscore{2} \|\hackscore{0} \le \tau\hackscore{2} \| P^{\frac{1}{2}}B v\hackscore{1} \|\hackscore{0}
  \end{equation}
- where $\tau\hackscore{2}$ is the given tolerane. The update step~\ref{P OPERATOR} involves the
+ where $\tau\hackscore{2}$ is the given tolerance. This translates into
- solution of a sparse matrix problem. We use the tolerance $\tau\hackscore{2}^2$ in order to make sure that any
+ \begin{equation} \label{STOKES P OPERATOR ERROR 2}
- error from solving this problem does not interfere with the PCG iteration.
+ \|e\hackscore{2}\|\hackscore{0} = \| B v\hackscore{2} \|\hackscore{0} \le M \tau\hackscore{2} \| B v\hackscore{1} \|\hackscore{0}
- We need to think about appropriate stopping criteria when solving
- \ref{SADDLEPOINT ITER STEP 1}, \ref{SADDLEPOINT ITER STEP 2} and~\ref{P OPERATOR}. We would like to use very weak convergence criteria to reduce computational costs but they need to be tight enough to not interfere with the
- convergence of the iteration process one level above. From Equation~\ref{SADDLEPOINT ITER STEP 1} we have
- \begin{equation}
- v\hackscore{1} = e\hackscore{1} + A^{-1} ( G - B^{*} p )
- \end{equation}
- We will use a sparse matrix solver so we have not full control on the norm $\|.\|\hackscore{s}$ used in the stopping criteria
- \begin{equation}
- \| G - A v\hackscore{1} - B^{*} p \|\hackscore{s} \le \tau\hackscore{1} \| G - A v - B^{*} p \|\hackscore{s}
  \end{equation}
- This translates into the conditoon
+ where $M$ is taking care of the fact that $P^{\frac{1}{2}}$ is dropped.
- \begin{equation}
- \| e\hackscore{1} \|\hackscore{1} \le K \tau\hackscore{1} \| dv\hackscore{1} - e\hackscore{1} \|\hackscore{1}
- \mbox{ therefore }
- \| e\hackscore{1} \|\hackscore{1} \le \frac{K \tau\hackscore{1}}{1-K \tau\hackscore{1}}  \| dv\hackscore{1} \|\hackscore{1}
- \end{equation}
- The constant $K$ represents some uncertainty combining a variety of unknown factors such as the
- solver being used and the condition number of the stiffness matrix.
- From the first equation of~\ref{SADDLEPOINT ITER STEP 2} we have
+ As we assume that there is no significant error from solving with the operator $A$ we have
- \begin{equation}
+ \begin{equation} \label{STOKES total V2}
- p\hackscore{2} =  p + (B A^{-1} B^{*})^{-1} (e\hackscore{2} + Bv\hackscore{1} ) =
- (B A^{-1} B^{*})^{-1} ( e\hackscore{2} + B e\hackscore{1} + B A^{-1} G)
- \end{equation}
- and simlar
- \begin{equation}
  v\hackscore{2} =  v\hackscore{1} - dv\hackscore{2}
  = v\hackscore{1}  - A^{-1} B^{*}dp
- = e\hackscore{1}  + A^{-1} ( G - B^{*} p ) - A^{-1} B^{*} (p\hackscore{2}-p)
- = e\hackscore{1}  + A^{-1} ( G - B^{*} p\hackscore{2})
  \end{equation}
- This shows that we can write the iterative process as a fixed point iteration to solve (assume all errors are zero)
+ Combining the equations~\ref{STOKES total V1},~\ref{STOKES total P2} and~\ref{STOKES total V2} and
+ setting the errors to zero we can write the solution process as a fix point problem
  \begin{equation}
  v = \Phi(v,p) \mbox{ and } p = \Psi(u,p)
  \end{equation}
- where
+ with suitable functions $\Phi(v,p)$ and $ \Psi(v,p)$ representing the itartion operator without
- \begin{equation}
+ errors. In fact for a linear problem,  $\Phi$ and $\Psi$ are constant. With this notation we can write
-  \begin{array}{rcl}
+ the update step in the form $p\hackscore{2}= \delta p + \Psi(v\hackscore{0},p\hackscore{0})$ and
- \Psi(v,p) & = &  (B A^{-1} B^{*})^{-1} B A^{-1} G \\
+ $v\hackscore{2}= \delta v + \Phi(v\hackscore{0},p\hackscore{0})$ where
- \Phi(u,p) & = & A^{-1} ( G - B^{*} (B A^{-1} B^{*})^{-1} B A^{-1} G )
+ the total error $\delta p$ and $\delta v$ are given as
- \end{array}
- \end{equation}
- Notice that if $A$ is independent from $v$ and $p$ the operators $\Phi(v,p)$ and $\Psi(u,p)$ are constant
- and threfore the iteration will terminate - providing no termination errors in sub-iterations - after one step.
- We also can give a formula for the error
  \begin{equation}
   \begin{array}{rcl}
  \delta p & = &  (B A^{-1} B^{*})^{-1} ( e\hackscore{2} + B e\hackscore{1} ) \\
- \delta v & = &   e\hackscore{1} -  A^{-1} B^{*}\delta p  \;.
+ \delta v & = &  e\hackscore{1} -  A^{-1} B^{*}\delta p  \;.
  \end{array}\label{STOKES ERRORS}
  \end{equation}
- Notice that $B\delta v = - e\hackscore{2}=-Bv\hackscore{2}$.
+ Notice that $B\delta v = - e\hackscore{2}=-Bv\hackscore{2}$. Our task is now to choose the tolerances
+ $\tau\hackscore{1}$ and $\tau\hackscore{2}$ such that the global errors $\delta p$ and $\delta v$
- With this notation
+ do not stop the convergence of the iteration process.
- \begin{equation}
-  \begin{array}{rcl}
+ To measure convergence we use
- v\hackscore{2} = \Phi(v,p) + \delta v \\
+ \begin{equation}
- p\hackscore{2} = \Psi(v,p) + \delta p \\
+ \epsilon = \max(\|v\hackscore{2}-v\|\hackscore{1}, \|B A^{-1} B^{*} (p\hackscore{2}-p)\|\hackscore{0})
+ \end{equation}
+ In practice using the fact that $B A^{-1} B^{*} (p\hackscore{2}-p\hackscore{0}) = B v\hackscore{1}$
+ and assuming that $v\hackscore{2}$ gives a better approximation to the true $v$ than
+ $v\hackscore{0}$ we will use the estimate
+ \begin{equation}
+ \epsilon = \max(\|v\hackscore{2}-v\hackscore{0}\|\hackscore{1}, \|B v\hackscore{1}\|\hackscore{0})
+ \end{equation}
+ to estimate the progress of the iteration step after the step is completed.
+ Note that the estimate of $\epsilon$
+ used in the stopping criteria~\ref{STOKES STOPPING CRITERIA}. If $\chi^{-}$ is the convergence rate assuming
+ exact calculations, i.e. $e\hackscore{1}=0$ and $e\hackscore{2}=0$, we are expecting
+ to maintain $\epsilon \le \chi^{-} \cdot \epsilon^{-}$. For the
+ pressure increment we get:
+ \begin{equation} \label{STOKES EST 1}
+ \begin{array}{rcl}
+ \|B A^{-1} B^{*} (p\hackscore{2}-p)\|\hackscore{0}
+  & \le & \|B A^{-1} B^{*} (p\hackscore{2}-\delta p-p)\|\hackscore{0} +
+ \|B A^{-1} B^{*} \delta p\|\hackscore{0} \\
+  & = & \chi^{-} \cdot \epsilon^{-} + \|e\hackscore{2} + B e\hackscore{1}\|\hackscore{0}  \\
+ & \approx & \chi^{-} \cdot \epsilon^{-} + \|e\hackscore{2}\|\hackscore{0} \\
+ & \le & \chi^{-} \cdot \epsilon^{-} + M \tau\hackscore{2} \|B v\hackscore{1}\|\hackscore{0} \\
  \end{array}
  \end{equation}
- We use the values $dv=v\hackscore{2}-v$ and $Bv\hackscore{1}=B A^{-1} B^{*} dp$ to measure the error of the
+ So we choose the value for $\tau\hackscore{2}$ from
- current approximation $v\hackscore{2}$ and $v\hackscore{2}$ towards the exact solution.
+ \begin{equation} \label{STOKES TOL2}
- Assuming that the iteration does converge with a convergence rate $\chi^{-}$ we have the estimate
+  M \tau\hackscore{2} \|B v\hackscore{1}\|\hackscore{0}  \le (\chi^{-})^2 \epsilon^{-}
- \begin{equation}
+ \end{equation}
- \max(\|dv\|\hackscore{1}, \|Bv\hackscore{1}\|\hackscore{0})
+ in order to make the perturbation for the termination of the pressure iteration a second order effect. We use a
- \le \chi^{-} \max(\|dv^{-}\|\hackscore{1}, \|Bv\hackscore{1}^{-}\|\hackscore{0})
+ similar argument for the velcity:
- + \max(\|\delta v\|\hackscore{1}, \|(B A^{-1} B^{*}) \delta p\|\hackscore{0})
+ \begin{equation}\label{STOKES EST 2}
- \end{equation}
+ \begin{array}{rcl}
- were the upper index $-$ referes to the increments in the last step.
+ \|v\hackscore{2}-v\|\hackscore{1} & \le & \|v\hackscore{2}-\delta v-v\|\hackscore{1} + \| \delta v\|\hackscore{1} \\
+ & \le &  \chi^{-} \cdot \epsilon^{-}  + \| e\hackscore{1} -  A^{-1} B^{*}\delta p \|\hackscore{1} \\
- Now we will establish estimates for $\|\delta v\|$ and $\|Bv\hackscore{1}\|$ from
+ & \approx &  \chi^{-} \cdot \epsilon^{-}  + \| e\hackscore{1} \|\hackscore{1} \\
- formulas~\ref{STOKES ERRORS} where neglect the $\delta p$ in the equation for $\delta v$ assuming that
+ & \le &  \chi^{-} \cdot \epsilon^{-}  +  K \tau\hackscore{1} \| dv\hackscore{1} - e\hackscore{1} \|\hackscore{1}
- $\delta p$ is controlled.
+ \\
- \begin{equation}
+ & \le &  ( 1  + K \tau\hackscore{1}) \chi^{-} \cdot \epsilon^{-}
- \|\delta v\|\hackscore{1} \approx \|e\hackscore{1}\|\hackscore{1} \le \frac{K \tau\hackscore{1}}{1-K \tau\hackscore{1}}  \| dv\hackscore{1} \|\hackscore{1} \approx \frac{K \tau\hackscore{1}}{1-K \tau\hackscore{1}} \| dv\|\hackscore{1}
+ \end{array}
- \end{equation}
- and
- \begin{equation}
- \|(B A^{-1} B^{*}) \delta p\|\hackscore{0}  \approx \| e\hackscore{2} \|\hackscore{0}
- = \| B v\hackscore{2} \|\hackscore{0} \le M \tau\hackscore{2} \| B v\hackscore{1} \|\hackscore{0}
- \end{equation}
- which leads to
- \begin{equation}
- \max(\|dv\|\hackscore{1}, \|Bv\hackscore{1}\|\hackscore{0})
- \le \frac{\chi^{-}}{1-max{(M \tau\hackscore{2}, \frac{K \tau\hackscore{1}}{1-K \tau\hackscore{1}})}} \max(\|dv^{-}\|\hackscore{1}, \|Bv\hackscore{1}^{-}\|\hackscore{0}) \label{STOKES ESTIMTED IMPROVEMENT}
- \end{equation}
- where the upper index $-$ refers to the previous iteration step.
- If we allow require $max{(M \tau\hackscore{2}, \frac{K \tau\hackscore{1}}{1-K \tau\hackscore{1}})}\le \gamma$
- with a given $0<\gamma<1$ we can set
- \begin{equation} \label{STOKES SET TOL}
- \tau\hackscore{2} = \frac{\gamma}{M} \mbox{ and } \tau\hackscore{1} = \frac{1}{K} \frac{\gamma}{1+\gamma}
- \end{equation}
- Problem is that we do not know $M$ and $K$ but can use the convergence rate
- \begin{equation}
-  \chi := \frac{\max(\|dv\|\hackscore{1}, \|Bv\hackscore{1}\|\hackscore{0})}{\max(\|dv^{-}\|\hackscore{1}, \|Bv\hackscore{1}^{-}\|\hackscore{0})}
- \end{equation}
- to construct estimates of $M$ and $K$. We look at the two cases where our prediction and choice of
- the tolerances was good or where we went wrong:
- \begin{equation}
- \chi \le \frac{\chi^{-}}{1-\gamma} \mbox{ or }
- \chi = \frac{\chi^{-}}{1-\gamma\hackscore{0}} >\frac{\chi^{-}}{1-\gamma}.
- \end{equation}
- which translates to
- \begin{equation}
- \gamma\hackscore{0} \le \gamma \mbox{ or } \gamma\hackscore{0}=1-\frac{\chi^{-}}{ \chi}>\gamma>0
- \end{equation}
- In the second case use \ref{STOKES SET TOL} for $\gamma=\gamma\hackscore{0}$ to get new estimates $M^{+}$ and $K^{+}$
- for $M$ and $K$:
- \begin{equation} \label{TOKES CONST UPDATE}
- M^{+} =
- \frac{\gamma\hackscore{0}}{\gamma} M
- \mbox{ and } K^{+} =
- \frac{\gamma\hackscore{0}}{\gamma}
- \frac{1+\gamma}{1+\gamma\hackscore{0}}
- K
- \mbox{ with } \gamma\hackscore{0}=\max(1-\frac{\chi^{-}}{ \chi},\gamma)
- \end{equation}
- With these updated constants we can now get better values for the tolerances via an updated value $\gamma^{+}$ for $\gamma$ and the corrected values $M^{+}$ and $K^{+}$ for $M$ and $K$. If we are in the case of convergence which we
- meassure by
- \begin{equation}
- \chi < \chi\hackscore{max}
- \end{equation}
- where $\chi\hackscore{max}$ is given value with $0<\chi\hackscore{max}<1$. We then consider the following
- criteria
- \begin{itemize}
-  \item As we are in case of convergence we try to relax the tolerances by increasing $\gamma$ by a factor $\frac{1}{\omega}$ ($0<\omega<1$). So we would like to choose
- $\gamma^{+} \ge \frac{\gamma}{\omega}$.
- \item We need to make sure that the estimated convergence rate based on $\gamma^{+}$ will still maintain convergence
- \begin{equation}
- \frac{\chi}{1-\gamma^{+}} \le \chi\hackscore{max}
- \end{equation}
- which translates into
- \begin{equation}
- \gamma^{+} \le 1 - \frac{\chi}{\chi\hackscore{max}}
- \end{equation}
- Notice that the right hand side is positive.
- \item We do not want to iterate more then necessary. So we would like reduce the error in the next just below the
- desired tolerance:
- \begin{equation}
- \frac{\chi}{1-\gamma^{+}}\max(\|dv\|\hackscore{1}, \|Bv\hackscore{1}\|\hackscore{0}) \ge f \cdot \tau \|v\hackscore{2}\|\hackscore{1}
  \end{equation}
- where the left hand side provides an estimate for the error to prediced in the next step. $f$ is a
+ So we choose the value for $\tau\hackscore{1}$ from
- safty factor. This leads to
+ \begin{equation} \label{STOKES TOL1}
- \begin{equation}
+ K \tau\hackscore{1} \le \chi^{-}
- \gamma^{+} \le
+ \end{equation}
--
+ Assuming we have etsimates for $M$ and $K$ (if none is available, we use the value $1$.)
- \chi \frac{\max(\|dv\|\hackscore{1}, \|Bv\hackscore{1}\|\hackscore{0})} { f \cdot \tau \|v\hackscore{2}\|\hackscore{1}}
+ we can use~\ref{STOKES TOL1} and~\ref{STOKES TOL2} to get appropriate values for the tolerances. After
+ the step has been completed we can calculate a new convergence rate $\chi =\frac{\epsilon}{\epsilon^{-}}$.
+ For partical reasons we restrict $\chi$ to be less or equal a given maximum value $\chi\hackscore{max}\le 1$.
+ 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
+ \begin{equation}\label{STOKES EST 3}
+ \chi \le ( 1 + \max(M \frac{\tau\hackscore{2} \|B v\hackscore{1}\|\hackscore{0}}{\chi^{-} \epsilon^{-}},K \tau\hackscore{1}  ) ) \cdot \chi^{-}
+ \end{equation}
+ If we assume that this inequality would be an equation if we would have chosen the right values
+ $M^{+}$ and $K^{+}$ then we get
+ \begin{equation}\label{STOKES EST 3b}
+ \chi =  ( 1 + \max(M^{+} \frac{\chi^{-}}{M},K^{+} \frac{\chi^{-}}{K}) ) \cdot \chi^{-}
  \end{equation}
- \end{itemize}
+ From this equation we set for
- Putting all these criteria together we get
+ than our choice for $K$ was not good anough. In this case we can calaculate a new value
- \begin{equation}
+  \begin{equation}
- \gamma^{+}=\min(\max(
+ K^{+} =  \frac{\chi-\chi^{-}}{(\chi^{-})^2} K
- \frac{\gamma}{\omega},
--
- \chi \frac{\max(\|dv\|\hackscore{1}, \|Bv\hackscore{1}\|\hackscore{0})} { f \cdot \tau \|v\hackscore{2}\|\hackscore{1}}), 1 - \frac{\chi}{\chi\hackscore{max}})
- \label{STOKES SET GAMMA PLUS}
  \end{equation}
- In case we cannot see convergence $\gamma$ is reduced by the factor $\omega$
+ In practice we will use
- \begin{equation}
+  \begin{equation}
- \gamma^{+}=\max(\omega \cdot \gamma, \gamma\hackscore{min})\label{STOKES SET GAMMA PLUS2}
+ K^{+}  = \max(\frac{\chi-\chi^{-}}{(\chi^{-})^2} K,\frac{1}{2}K,1)
  \end{equation}
- where $\gamma\hackscore{min}$ is introduced to make sure that $\gamma^{+}$ stays away from zero.
+ where the second term is used to reduce a potential overestimate of $K$.
+ The same identity is used for to update $M$. The updated $M^{+}$ and $K^{+}$
+ are then use in the next iteration step to control the tolerances.
- Appling~\ref{STOKES SET TOL} for $\gamma^{+}$ and~\ref{TOKES CONST UPDATE} we get the update formula
- \begin{equation} \label{STOKES CONST UPDATE}
+ In some cases one can observe that there is a significant change
- \tau\hackscore{2}^{+} =
+ in the velocity but the new velocity $v\hackscore{1}$ has still a
- \frac{\gamma^{+}}{\gamma\hackscore{0}} \tau\hackscore{2}
+ small divergence, i.e. we have
- \mbox{ and } \tau\hackscore{1}^{+} =
+ $\|Bv\hackscore{1}\|\hackscore{0} \ll \|v\hackscore{1}-v\hackscore{0}\|\hackscore{1}$.
- \frac{\gamma^{+}}{\gamma\hackscore{0}}
+ In this case we will get a small pressure increment and consequently only very small changes to
- \frac{1+\gamma\hackscore{0}}{1+\gamma^{+}}
+ the velocity as a result of the second update step which therefore can be skipped and
- \tau\hackscore{1}
+ we can directly repeat the first update step until the increment in velocity becomes
+ significant relative to its divergence. In practice we will ignore the second half of the iteration step
+ as long as
+  \begin{equation}\label{STOKES LARGE BV1}
+ \|Bv\hackscore{1}\|\hackscore{0} \le \theta \cdot \|v\hackscore{1}-v\hackscore{0}\|
  \end{equation}
- to get the new tolerances $\tau\hackscore{1}^{+}$ and $\tau\hackscore{2}^{+}$ to be used in the next iteration step.
+ where $0<\theta<1$ is a given factor. In this case we will also check the stopping criteria
- Notice that the coefficients $M$ and $K$ have been eliminated. The calculation of $\gamma\hackscore{0}$ requires
+ with $v\hackscore{1}\rightarrow v\hackscore{2}$ but we will not correct $M$ in this case.
- at least two iteration steps (or a previous time step).
- Typical initial values are
+ Starting from an initial guess $v\hackscore{0}$ and $p\hackscore{0}$ for velocity and pressure
- \begin{equation} \label{STOKES VALUES}
+ the solution procedure is implemented as follows.
- \tau\hackscore{1}=\tau\hackscore{1}=\frac{1}{100} ; \; \omega=\frac{1}{2} \; \gamma=\frac{1}{10} ;\gamma\hackscore{min}=10^{-6}
+ \begin{enumerate}
- \end{equation}
+  \item calculate viscosity $\eta(v\hackscore{0},p)\hackscore{0}$ and assemble operator $A$ from $\eta$.
- where after the first step we set $\gamma\hackscore{0}=\gamma$ as no $\chi^{-}$ is available.
+  \item calculate the tolerance $\tau\hackscore{1}$ from~\ref{STOKES TOL1}.
+  \item Solve~\ref{STOKES ITER STEP 1} for $dv\hackscore{1}$ with tolerance $\tau\hackscore{1}$.
+  \item update $v\hackscore{1}= v\hackscore{0}+ dv\hackscore{1}$
+  \item if $Bv\hackscore{1}$ is large (see~\ref{STOKES LARGE BV1}):
+  \begin{enumerate}
+  \item calculate the tolerance $\tau\hackscore{2}$ from~\ref{STOKES TOL2}.
+  \item Solve~\ref{STOKES ITER STEP 2} for $dp\hackscore{2}$ and $v\hackscore{2}$ with tolerance $\tau\hackscore{2}$.
+  \item update $p\hackscore{2}\leftarrow p\hackscore{0}+ dp\hackscore{2}$
+  \end{enumerate}
+  \item else:
+   \begin{enumerate}
+   \item update $p\hackscore{2}\leftarrow p$ and $v\hackscore{2}\leftarrow v\hackscore{1}$
+    \end{enumerate}
+    \item calculate convergence messure $\epsilon$ and convergence rate $\chi$
+ \item if stopping criterium~\ref{STOKES STOPPING CRITERIA} holds:
+  \begin{enumerate}
+  \item return $v\hackscore{2}$ and $p\hackscore{2}$
+  \end{enumerate}
+  \item else:
+  \begin{enumerate}
+      \item update $M$ and $K$
+      \item goto step 1 with $v\hackscore{0}\leftarrow v\hackscore{2}$ and $p\hackscore{0}\leftarrow p\hackscore{2}$.
+ \end{enumerate}
+ \end{enumerate}
  \subsection{Functions}
-Line 461 
 sreturns the current absolute tolerance.
+Line 346 
 sreturns the current absolute tolerance.
  \end{methoddesc}
  \begin{methoddesc}[StokesProblemCartesian]{getSolverOptionsVelocity}{}
- returns the solver options used  solve the equations~(\ref{V CALC}) for velocity.
+ returns the solver options used  solve the equations~(\ref{STOKES ITER STEP 1}) and~(\ref{STOKES P OPERATOR}) for velocity.
  \end{methoddesc}
  \begin{methoddesc}[StokesProblemCartesian]{getSolverOptionsPressure}{}
- returns the solver options used  solve the equation~(\ref{P PREC}) for pressure.
+ returns the solver options used solve the preconditioner equation~(\ref{STOKES P PREC}) for pressure.
  \end{methoddesc}
  \begin{methoddesc}[StokesProblemCartesian]{getSolverOptionsDiv}{}

 Legend:



Removed from v.2793
 


changed lines


 
Added in v.2843
 Legend:



Removed from v.2793
 


changed lines


 
Added in v.2843
-Removed from v.2793
+Added in v.2843

	ViewVC Help
Powered by ViewVC 1.1.26