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 
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.
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}  
 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}  
 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}  
 which can be overwritten by constraints of the form  
 \begin{equation}\label{Stokes Boundary0}  
 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.  
 \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}  
 \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]  
 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  
 \|v\hackscore{k,k}\| \hackscore \le  \epsilon  
 \label{STOKES STOP}  
 where $\epsilon$ is the desired relative accuracy and  
 \|p\|^2= \int\hackscore{\Omega} p^2 \; dx  
 \label{PRESSURE NORM}  
 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  
 \label{V CALC}  
 which is inserted into the second equation leading to  
 S p =  B A^{-1} G  
 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  
 see~\cite{ELMAN} for more details. Note that the residual for the current approximation $p$ is given as  
 r=B A^{-1} (G - B^* p) = Bv  
 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  
 \hat{S}^{-1} S p =  \hat{S}^{-1} B A^{-1} G  
 We use the norm  
 \|p\|\hackscore{GMRES} = \|\hat{S} p \|  
 Notice that for the residual $\hat{r}=\hat{S}^{-1} r$ one has  
 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  
 \|\hat{r}\|\hackscore{GMRES} \le ATOL  
 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  
 \|r\|\hackscore{PCG}^2 = \int\hackscore{\Omega} r \hat{S}^{-1}r \; dx = \int\hackscore{\Omega}  \eta r^2 \; dx  
 To take the extra factor $\eta$ into consideration when checking the stopping criterion we use the following  
 definition for $ATOL$:  
 ATOL = \epsilon \frac{\|\sqrt{v^{0}\hackscore{j,k}v^{0}\hackscore{j,k}}\|  }{\|v^{0}\hackscore{k,k}\|}  
 \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.  
 \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.  
 \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.  
 sets the tolerance in an appropriate norm relative to the right hand side. The tolerance must be non-negative and less than 1.  
 returns the current relative tolerance.  
 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.  
 sreturns the current absolute tolerance.  
 returns the solver options used  solve the equations~(\ref{V CALC}) for velocity.  
 returns the solver options used  solve the equation~(\ref{P PREC}) for pressure.  
 set the solver options for solving the equation to project the divergence of the velocity onto the function space of pressure.  
 \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:  
 from esys.escript import *  
 from esys.finley import Rectangle  
 from esys.escript.models import StokesProblemCartesian  
 dom = Rectangle(NE,NE,order=2)  
 x = dom.getX()  
 mask= (whereZero(x[0])*[1.,0]+whereZero(x[0]-1))*[1.,0] + \  
 sc.initialize(eta=.1, fixed_u_mask= mask)  
 v,p=sc.solve(v,p, verbose=True)  
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

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

  ViewVC Help
Powered by ViewVC 1.1.26