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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 2748 by gross, Tue Nov 17 07:32:59 2009 UTC revision 2793 by gross, Tue Dec 1 06:10:10 2009 UTC
# Line 26  The index $i$ may depend on the location Line 26  The index $i$ may depend on the location
26  $v^D$ is a given function on the domain.  $v^D$ is a given function on the domain.
27    
28  \subsection{Solution Method \label{STOKES SOLVE}}  \subsection{Solution Method \label{STOKES SOLVE}}
29  In block form equation equations~\ref{Stokes 1} and~\ref{Stokes 2} takes the form of a saddle point problem  If we assume that $\eta$ is independent from the velocity and pressure equations~\ref{Stokes 1} and~\ref{Stokes 2}
30  \index{saddle point problem}  can be written in the block form
31  \begin{equation}  \begin{equation}
32  \left[ \begin{array}{cc}  \left[ \begin{array}{cc}
33  A     & B^{*} \\  A     & B^{*} \\
# Line 43  G \\ Line 43  G \\
43  \end{array} \right]  \end{array} \right]
44  \label{SADDLEPOINT}  \label{SADDLEPOINT}
45  \end{equation}  \end{equation}
46  where $A$ is coercive (assuming $A$ is not depending on $v$ or $p$), 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
50    velocity and pressure we calculate a correction $dv$ for the velocity by solving the first
51    equation of equation~\ref{SADDLEPOINT}
52     \begin{equation}\label{SADDLEPOINT ITER STEP 1}
53     A dv\hackscore{1} = G - A v\hackscore{0} - B^{*} p\hackscore{0}
54    \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}$
56    for the pressure and an additional correction $dv\hackscore{2}$ for the velocity by solving
57     \begin{equation}
58     \begin{array}{rcl}
59     B A^{-1} B^{*} dp\hackscore{2} & = & Bv\hackscore{1} \\
60     A dv\hackscore{2} & = & B^{*} dp\hackscore{2}
61    \end{array}
62     \label{SADDLEPOINT ITER STEP 2}
63     \end{equation}
64    The new velocity and pressure are then given by $v=v\hackscore{1}-dv\hackscore{2}$ and
65    $p=p\hackscore{0}+dp\hackscore{2}$ which will fullfill the block system~\ref{SADDLEPOINT}.
66    This solution strategy is called the Uzawa scheme \index{Uzawa scheme}. There are
67    two problems with this scheme. Firstly, in practice we will use an iterative scheme
68    to solve the problem for operator $A$. So we will be unable to calculate the operator
69    $ B A^{-1} B^{*}$ required for $dp\hackscore{2}$. In fact, we need to use an iterative scheme
70    to solve the first equation in~\ref{SADDLEPOINT ITER STEP 2} where in each iteration step
71    an iterative solver for $A$ is applied. The second issue is that
72    viscosity $\eta$ may depend of velocity or pressure and so we need to iterate over the
73    three equations~\ref{SADDLEPOINT ITER STEP 1} and~\ref{SADDLEPOINT ITER STEP 2}. Using the
74    two norms
75    \begin{equation}
76    \|v\|\hackscore{1}^2 = \int\hackscore{\Omega} v\hackscore{j,k}v\hackscore{j,k} \; dx
77    \mbox{ and }
78    \|p\|\hackscore{0}^2= \int\hackscore{\Omega} p^2 \; dx.
79    \label{STOKES STOP}
80    \end{equation}
81    to define the stopping criterium
82     \begin{equation}
83    \max(\|Bv\hackscore{1}\|\hackscore{0},\|v-v\hackscore{0}\|\hackscore{1}) \le \tau \cdot \|v\|\hackscore{1}
84     \end{equation}
85    to terminate the overall iteration with a overall given tolerance $0<\tau<1$.
86    Notice that because of the first equation of
87    and~\ref{SADDLEPOINT ITER STEP 2} $\|Bv\hackscore{1}\|\hackscore{0}$ is equal to the
88    norm of $B A^{-1} B^{*} dp\hackscore{2}$ and consequently provides a norm for the  pressure correction.
89    
90    Let us first have a look at the solution process~\ref{SADDLEPOINT ITER STEP 1} which needs to resolve the
91    non-linearity of the viscosity coeficient. We will run a iterative process by resolving
92    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
93    the iteration process over the second set of equation if $\|Bv\hackscore{1}\|\hackscore{0}$
94    is smaller then $\|v\hackscore{1}-v\hackscore{0}\|\hackscore{1}$. So we will terminate the non-linear iteration
95    process of the first equation if
96     \begin{equation}
97    \theta \cdot \|v\hackscore{1}-v\hackscore{0}\| \le \|Bv\hackscore{1}\|\hackscore{0}
98    \mbox{ or }
99    \|v\hackscore{1}-v\hackscore{0}\|\hackscore{1} \le \tau \cdot \|v\hackscore{1}\|\hackscore{1}
100    \end{equation}
101    where $0<\theta<1$ is a given factor (typically $\theta=0.5$).
102    
103    We need to think about appropriate stopping criteria when solving
104    \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
105    convergence of the iteration process one level above.
106    We will use a sparse matrix solver so we have not full control on the norm $\|.\|\hackscore{s}$ used in the stopping criteria
107    \begin{equation}
108    \| 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}
109    \end{equation}
110    If $e\hackscore{1}$ is the error of the returned solution $dv\hackscore{1}$ this translates into the condition
111    \begin{equation}
112    \| e\hackscore{1} \|\hackscore{1} \le K \tau\hackscore{1} \| dv\hackscore{1} - e\hackscore{1} \|\hackscore{1}
113    \end{equation}
114    The constant $K$ represents some uncertainty combining a variety of unknown factors such as the
115    solver being used and the condition number of the stiffness matrix.  This leads to the estimate
116    \begin{equation}
117    \| dv\hackscore{1} \|\hackscore{1} \le (1+K \tau\hackscore{1}) \chi \| dv\hackscore{1}^- \|\hackscore{1}
118    \label{DV COND}
119    \end{equation}
120    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
121    to have $K \tau\hackscore{1} \le \chi$ which leads to the condition
122     \begin{equation}
123     \tau\hackscore{1} \le \frac{\min(\chi,\frac{1}{2})}{K}
124     \label{NEW TAU 1}
125    \end{equation}
126    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
127    access to the estimate of the convergence rate $\chi^-$ of the previous iteration step we
128    get
129     \begin{equation}
130     \frac{1}{K} = \frac{\chi^-}{\chi-\chi^{-}} \tau\hackscore{1}^-
131    \end{equation}
132    if $\chi \ge \chi^- (1 + \chi^-) $ which means that our prediction for a suitable tolerane based on the avaibale $K$ value was wrong.
133    So we start with the value $K=1$ in~\ref{NEW TAU 1} and recalculate
134    $K$ if we have underestimate $K$.
135    
136    
137    
138    
139    
140    
141    
142    If $dv\hackscore{1}$
143    From Equation~\ref{SADDLEPOINT ITER STEP 1} we have
144    \begin{equation}
145    v\hackscore{1} = e\hackscore{1} + A^{-1} ( G - B^{*} p )
146    \end{equation}
147    This translates into the conditoon
148    \begin{equation}
149    \| e\hackscore{1} \|\hackscore{1} \le K \tau\hackscore{1} \| dv\hackscore{1} - e\hackscore{1} \|\hackscore{1}
150    \mbox{ therefore }
151    \| e\hackscore{1} \|\hackscore{1} \le \frac{K \tau\hackscore{1}}{1-K \tau\hackscore{1}}  \| dv\hackscore{1} \|\hackscore{1}
152    \end{equation}
153    The constant $K$ represents some uncertainty combining a variety of unknown factors such as the
154    solver being used and the condition number of the stiffness matrix.  
155    
156    
157    
158    
159    
160    ===================================
161    
162    
163  We use iterative techniques to solve this problem: Given an approximation $v$ and $p$ for  We use iterative techniques to solve this problem: Given an approximation $v$ and $p$ for
164  velocity and pressure we perform the following steps in the Uzawa scheme \index{Uzawa scheme} style:  velocity and pressure we perform the following steps in the Uzawa scheme \index{Uzawa scheme} style:
165  \begin{enumerate}  \begin{enumerate}
# Line 66  velocity and pressure we perform the fol Line 181  velocity and pressure we perform the fol
181   \item update $p\hackscore{2}\leftarrow p+ dp$ and $v\hackscore{2}= v\hackscore{1} - dv\hackscore{2}$   \item update $p\hackscore{2}\leftarrow p+ dp$ and $v\hackscore{2}= v\hackscore{1} - dv\hackscore{2}$
182   \item goto Step 0 with $v=v\hackscore{2}$ and $p=p\hackscore{2}$.   \item goto Step 0 with $v=v\hackscore{2}$ and $p=p\hackscore{2}$.
183  \end{enumerate}  \end{enumerate}
184  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  where $\tau$ is the given tolerance. Notice that the operator $A$ if it is depending on $v$ or $p$ is updated once only.
185    In this algorithm
186  \begin{equation}  \begin{equation}
187  \|v\|\hackscore{1}^2 = \int\hackscore{\Omega} v\hackscore{j,k}v\hackscore{j,k} \; dx  \|v\|\hackscore{1}^2 = \int\hackscore{\Omega} v\hackscore{j,k}v\hackscore{j,k} \; dx
188  \mbox{ and }  \mbox{ and }
# Line 150  We also can give a formula for the error Line 266  We also can give a formula for the error
266  \begin{equation}  \begin{equation}
267   \begin{array}{rcl}   \begin{array}{rcl}
268  \delta p & = &  (B A^{-1} B^{*})^{-1} ( e\hackscore{2} + B e\hackscore{1} ) \\  \delta p & = &  (B A^{-1} B^{*})^{-1} ( e\hackscore{2} + B e\hackscore{1} ) \\
269  \delta v & = &   e\hackscore{1} -  A^{-1} B^{*}\delta p  \delta v & = &   e\hackscore{1} -  A^{-1} B^{*}\delta p  \;.
270  \end{array}\label{STOKES ERRORS}  \end{array}\label{STOKES ERRORS}
271  \end{equation}  \end{equation}
272  Notice that $B\delta v = - e\hackscore{2}=-Bv\hackscore{2}$  Notice that $B\delta v = - e\hackscore{2}=-Bv\hackscore{2}$.
273    
274  With this notation  With this notation
275  \begin{equation}  \begin{equation}
276   \begin{array}{rcl}   \begin{array}{rcl}
# Line 161  v\hackscore{2} = \Phi(v,p) + \delta v \\ Line 278  v\hackscore{2} = \Phi(v,p) + \delta v \\
278  p\hackscore{2} = \Psi(v,p) + \delta p \\  p\hackscore{2} = \Psi(v,p) + \delta p \\
279  \end{array}  \end{array}
280  \end{equation}  \end{equation}
281  We use the $dv=v\hackscore{2}-v$ and $Bv\hackscore{1}=B A^{-1} B^{*} dp$ to measure the error of the  We use the values $dv=v\hackscore{2}-v$ and $Bv\hackscore{1}=B A^{-1} B^{*} dp$ to measure the error of the
282  current approximation $v\hackscore{2}$ and $v\hackscore{2}$ towards the exact solution.  current approximation $v\hackscore{2}$ and $v\hackscore{2}$ towards the exact solution.
283  Assuming that the iteration does converge with a convergence rate $\chi^{-}$ we have the estimate  Assuming that the iteration does converge with a convergence rate $\chi^{-}$ we have the estimate
284  \begin{equation}  \begin{equation}
285  \max(\|dv\|\hackscore{1}, \|Bv\hackscore{1}\|\hackscore{0})  \max(\|dv\|\hackscore{1}, \|Bv\hackscore{1}\|\hackscore{0})
286  \le \chi^{-} \max(\|dv^{-}\|\hackscore{1}, \|Bv\hackscore{1}^{-}\|\hackscore{0})  \le \chi^{-} \max(\|dv^{-}\|\hackscore{1}, \|Bv\hackscore{1}^{-}\|\hackscore{0})
287  + \max(\|\delta v\|\hackscore{1} + \|(B A^{-1} B^{*}) \delta p\|\hackscore{0})  + \max(\|\delta v\|\hackscore{1}, \|(B A^{-1} B^{*}) \delta p\|\hackscore{0})
288  \end{equation}  \end{equation}
289  were the upper index $-$ referes to the increments in the last step.  were the upper index $-$ referes to the increments in the last step.
290  Now we try to establish estimates for $\|\delta v\|$ and $\|Bv\hackscore{1}\|$ from  
291    Now we will establish estimates for $\|\delta v\|$ and $\|Bv\hackscore{1}\|$ from
292  formulas~\ref{STOKES ERRORS} where neglect the $\delta p$ in the equation for $\delta v$ assuming that  formulas~\ref{STOKES ERRORS} where neglect the $\delta p$ in the equation for $\delta v$ assuming that
293  $\delta p$ is controlled.  $\delta p$ is controlled.
294  \begin{equation}  \begin{equation}
# Line 287  where after the first step we set $\gamm Line 405  where after the first step we set $\gamm
405  \subsection{Functions}  \subsection{Functions}
406    
407  \begin{classdesc}{StokesProblemCartesian}{domain}  \begin{classdesc}{StokesProblemCartesian}{domain}
408  opens the Stokes problem\index{Stokes problem} on the \Domain domain. The approximation  opens the Stokes problem\index{Stokes problem} on the \Domain domain. The domain
409  order needs to be two. Order two will be used to approximate the velocity and order one  needs to support LBB compliant elements for the Stokes problem, see~\cite{LBB} for detials~\index{LBB condition}.
410  to approximate the pressure.  For instance one can use second order polynomials for velocity and
411    first order polynomials for the pressure on the same element. Alternativly, one can use
412    macro elements\index{macro elements} using linear polynomial for both pressure and velocity bu with a subdivided
413    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
414    \begin{python}
415    velocity=Vector(0.0, Solution(mesh))
416    pressure=Scalar(0.0, ReducedSolution(mesh))
417    \end{python}
418  \end{classdesc}  \end{classdesc}
419    
420  \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(), \optional{fixed_u_mask=Data(), \optional{eta=1, \optional{surface_stress=Data(), \optional{stress=Data()}, \optional{
# Line 306  The method will try to cast the given va Line 431  The method will try to cast the given va
431  \begin{methoddesc}[StokesProblemCartesian]{solve}{v,p  \begin{methoddesc}[StokesProblemCartesian]{solve}{v,p
432  \optional{, max_iter=100 \optional{, verbose=False \optional{, usePCG=True }}}}  \optional{, max_iter=100 \optional{, verbose=False \optional{, usePCG=True }}}}
433  solves the problem and return approximations for velocity and pressure.  solves the problem and return approximations for velocity and pressure.
434  The arguments \var{v} and \var{p} define initial guess. The values of \var{v} marked  The arguments \var{v} and \var{p} define initial guess.
435    \var{v} must have function space \var{Solution(domain)} and
436    \var{p} must have function space \var{ReducedSolution(domain)}.
437    The values of \var{v} marked
438  by \var{fixed_u_mask} remain unchanged.  by \var{fixed_u_mask} remain unchanged.
439  If \var{usePCG} is set to \True  If \var{usePCG} is set to \True
440  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  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
Legend:
Removed from v.2748  
changed lines
  Added in v.2793
  ViewVC Help
Powered by ViewVC 1.1.26