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

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

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

revision 2718 by jfenwick, Mon Sep 7 03:39:45 2009 UTC revision 2719 by gross, Wed Oct 14 06:38:03 2009 UTC
# Line 17  Line 17 
17    
18  The following sections give a breif overview of the model classes and their corresponding methods.  The following sections give a breif overview of the model classes and their corresponding methods.
19    
20  \section{Stokes Problem}  \input{stokessolver}
 The velocity \index{velocity} field $v$ and pressure $p$ of an incompressible fluid \index{incompressible fluid} is 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}  
 \end{equation}  
 where $\eta$ is the viscosity, $F\hackscore{i}$ defines an internal force \index{force, internal} and $\sigma\hackscore{ij}$ is an intial stress \index{stress, initial}. We assume an incompressible media:  
 \begin{equation}\label{Stokes 2}  
 -v\hackscore{i,i}=0  
 \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}+\sigma\hackscore{ij} n\hackscore{j}  
 \end{equation}  
 which can be overwritten by constraints of the form  
 \begin{equation}\label{Stokes Boundary0}  
 v\hackscore{i}(x)=v^D\hackscore{i}(x)  
 \end{equation}  
 at some locations $x$ at the boundary of the domain. The index $i$ may depend on the location $x$ on the boundary.  
 $v^D$ is a given function on the domain.  
21    
 \subsection{Solution Method \label{STOKES SOLVE}}  
 In block form equation equations~\ref{Stokes 1} and~\ref{Stokes 2} takes the form of a saddle point problem  
 \index{saddle point problem}  
 \begin{equation}  
 \left[ \begin{array}{cc}  
 A     & B^{*} \\  
 B & 0 \\  
 \end{array} \right]  
 \left[ \begin{array}{c}  
 v \\  
 p \\  
 \end{array} \right]  
 =\left[ \begin{array}{c}  
 G \\  
 0 \\  
 \end{array} \right]  
 \label{SADDLEPOINT}  
 \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}.  
 We use iterative techniques to solve this problem. To make sure that the incomressibilty condition holds  
 with sufficient accuracy we check for  
 \begin{equation}  
 \|v\hackscore{k,k}\| \hackscore \le  \epsilon  
 \|\sqrt{v\hackscore{j,k}v\hackscore{j,k}}\|  
 \label{STOKES STOP}  
 \end{equation}  
 where $\epsilon$ is the desired relative accuracy and  
 \begin{equation}  
 \|p\|^2= \int\hackscore{\Omega} p^2 \; dx  
 \label{PRESSURE NORM}  
 \end{equation}  
 defines the $L^2$-norm. We use the Uzawa scheme \index{Uzawa scheme} to solve the problem.  
   
 In fact the first equation in~\ref{SADDLEPOINT} gives for a known pressure  
 \begin{equation}  
 v=A^{-1}(G-B^{*}p)  
 \label{V CALC}  
 \end{equation}  
 which is inserted into the second equation leading to  
 \begin{equation}  
 S p =  B A^{-1} G  
 \end{equation}  
 with the Schur complement \index{Schur complement} $S=BA^{-1}B^{*}$. This problem can be solved iteratively  
 with the preconditioner $\hat{S}$ defined as $q=\hat{S}^{-1}p$ by solving  
 \begin{equation} \label{P PREC}  
 \frac{1}{\eta}q = p  
 \end{equation}  
 see~\cite{ELMAN} for more details. Note that the residual for the current approximation $p$ is given as  
 \begin{equation}  
 r=B A^{-1} (G - B^* p) = Bv  
 \end{equation}  
 where $v$ is given by~\ref{V CALC}.  
   
 If one uses the generalized minimal residual method (GMRES) \index{generalized minimal residual method!GMRES}  
 the method is directly applied to the preconditioned system  
 \begin{equation}  
 \hat{S}^{-1} S p =  \hat{S}^{-1} B A^{-1} G  
 \end{equation}  
 We use the norm  
 \begin{equation}  
 \|p\|\hackscore{GMRES} = \|\hat{S} p \|  
 \end{equation}  
 Notice that for the residual $\hat{r}=\hat{S}^{-1} r$ one has  
 \begin{equation}  
 \  
 \end{equation}  
 If $p^{0}$ provides an initial guess for the pressure we use~\ref{V CALC} to get a first initial guess for the  
 velocity $v^{0}$ which we use to set an absolute tolerance $ATOL =\epsilon \|\sqrt{v^{0}\hackscore{j,k}v^{0}\hackscore{j,k}}\|$.  
 The GMRES is terminated when  
 \begin{equation}  
 \|\hat{r}\|\hackscore{GMRES} \le ATOL  
 \end{equation}  
 Notice that $\|\hat{r}\|\hackscore{GMRES}= \|r \| = \|Bv\|  = \|v\hackscore{k,k}\|$ so we we can expect that  
 the target stopping criterion~\ref{STOKES STOP} is fullfilled. However, if $v$ is very different from the  
 initial choice of $v^{0}$ the value of $ATOL$ is corrected and GMRES is restarted with a new tolerance. For time dependend problems this apprach works well as value for $p$ form a previous time step provides a good initial guess.  
   
 Alternatively, as $S$ is symmetric and positive definite one can apply the preconditioned conjugate gradient method (PCG) \index{preconditioned conjugate gradient method!PCG}. PCG use the norm  
 \begin{equation}  
 \|r\|\hackscore{PCG}^2 = \int\hackscore{\Omega} r \hat{S}^{-1}r \; dx = \int\hackscore{\Omega}  \eta r^2 \; dx  
 \end{equation}  
 To take the extra factor $\eta$ into consideration when checking the stopping criterion we use the following  
 definition for $ATOL$:  
 \begin{equation}  
 ATOL = \epsilon \frac{\|\sqrt{v^{0}\hackscore{j,k}v^{0}\hackscore{j,k}}\|  }{\|v^{0}\hackscore{k,k}\|}  
 \|v^{0}\hackscore{k,k}\|\hackscore{PCG}  
 \end{equation}  
   
   
   
 \subsection{Functions}  
   
 \begin{classdesc}{StokesProblemCartesian}{domain\optional{, adaptSubTolerance=\True}}  
 opens the Stokes problem\index{Stokes problem} on the \Domain domain. The approximation  
 order needs to be two.  
 If \var{adaptSubTolerance} is \True  
 the tolerances for all subproblems are set automatically.  
 \end{classdesc}  
   
 \begin{methoddesc}[StokesProblemCartesian]{initialize}{\optional{f=Data(), \optional{fixed_u_mask=Data(), \optional{eta=1, \optional{surface_stress=Data(), \optional{stress=Data()}}}}}}  
 assigns values to the model parameters. In any call all values must be set.  
 \var{f} defines the external force $f$, \var{eta} the viscosity $\eta$,  
 \var{surface_stress} the surface stress $s$ and \var{stress} the initial stress $\sigma$.  
 The locations and compontents where the velocity is fixed are set by  
 the values of \var{fixed_u_mask}. The method will try to cast the given values to appropriate  
 \Data class objects.  
 \end{methoddesc}  
   
 \begin{methoddesc}[StokesProblemCartesian]{solve}{v,p  
 \optional{, max_iter=100 \optional{, verbose=False \optional{, usePCG=True }}}}  
 solves the problem and return approximations for velocity and pressure.  
 The arguments \var{v} and \var{p} define initial guess. The values of \var{v} marked  
 by \var{fixed_u_mask} remain unchanged.  
 If \var{usePCG} is set to \True  
 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  
 the PCG scheme is more efficient.  
 \var{max_iter} defines the maximum number of iteration steps.  
   
 If \var{verbose} is set to \True informations on the progress of of the solver are printed.  
 \end{methoddesc}  
   
   
 \begin{methoddesc}[StokesProblemCartesian]{setTolerance}{\optional{tolerance=1.e-4}}  
 sets the tolerance in an appropriate norm relative to the right hand side. The tolerance must be non-negative and less than 1.  
 \end{methoddesc}  
 \begin{methoddesc}[StokesProblemCartesian]{getTolerance}{}  
 returns the current relative tolerance.  
 \end{methoddesc}  
 \begin{methoddesc}[StokesProblemCartesian]{setAbsoluteTolerance}{\optional{tolerance=0.}}  
 sets the absolute tolerance for the error in the relevant norm. The tolerance must be non-negative. Typically the  
 absolute talerance is set to 0.  
 \end{methoddesc}  
   
 \begin{methoddesc}[StokesProblemCartesian]{getAbsoluteTolerance}{}  
 sreturns the current absolute tolerance.  
 \end{methoddesc}  
   
 \begin{methoddesc}[StokesProblemCartesian]{getSolverOptionsVelocity}{}  
 returns the solver options used  solve the equations~(\ref{V CALC}) for velocity.  
 \end{methoddesc}  
   
 \begin{methoddesc}[StokesProblemCartesian]{getSolverOptionsPressure}{}  
 returns the solver options used  solve the equation~(\ref{P PREC}) for pressure.  
 \end{methoddesc}  
   
 \begin{methoddesc}[StokesProblemCartesian]{getSolverOptionsDiv}{}  
 set the solver options for solving the equation to project the divergence of the velocity onto the function space of pressure.  
 \end{methoddesc}  
   
   
 \subsection{Example: Lit Driven Cavity}  
  The following script \file{lit\hackscore driven\hackscore cavity.py}  
 \index{scripts!\file{helmholtz.py}} which is available in the \ExampleDirectory  
 illustrates the usage of the \class{StokesProblemCartesian} class to solve  
 the lit driven cavity problem:  
 \begin{python}  
 from esys.escript import *  
 from esys.finley import Rectangle  
 from esys.escript.models import StokesProblemCartesian  
 NE=25  
 dom = Rectangle(NE,NE,order=2)  
 x = dom.getX()  
 sc=StokesProblemCartesian(dom)  
 mask= (whereZero(x[0])*[1.,0]+whereZero(x[0]-1))*[1.,0] + \  
       (whereZero(x[1])*[0.,1.]+whereZero(x[1]-1))*[1.,1]  
 sc.initialize(eta=.1, fixed_u_mask= mask)  
 v=Vector(0.,Solution(dom))  
 v[0]+=whereZero(x[1]-1.)  
 p=Scalar(0.,ReducedSolution(dom))  
 v,p=sc.solve(v,p, verbose=True)  
 saveVTK("u.xml",velocity=v,pressure=p)  
 \end{python}  
22    
23  \section{Darcy Flux}  \section{Darcy Flux}
24  We want to calculate the velocity $u$ and pressure $p$ on a domain $\Omega$ solving  We want to calculate the velocity $u$ and pressure $p$ on a domain $\Omega$ solving

Legend:
Removed from v.2718  
changed lines
  Added in v.2719

  ViewVC Help
Powered by ViewVC 1.1.26