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

revision 3502 by gross, Thu Apr 28 05:06:24 2011 UTC revision 3911 by jfenwick, Thu Jun 14 01:01:03 2012 UTC
# Line 1  Line 1
1
2  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3  %  %
4  % Copyright (c) 2003-2010 by University of Queensland  % Copyright (c) 2003-2012 by University of Queensland
5  % Earth Systems Science Computational Center (ESSCC)  % Earth Systems Science Computational Center (ESSCC)
6  % http://www.uq.edu.au/esscc  % http://www.uq.edu.au/esscc
7  %  %
# Line 13  Line 13
13
14  \section{Darcy Flux}  \section{Darcy Flux}
15  \label{DARCY FLUX}  \label{DARCY FLUX}
16  We want to calculate the velocity $u$ and pressure $p$ on a domain $\Omega$  We want to calculate the flux $u$ and pressure $p$ on a domain $\Omega$
17  solving the Darcy flux problem\index{Darcy flux}\index{Darcy flow}  solving the Darcy flux problem\index{Darcy flux}\index{Darcy flow}
18  \begin{equation}\label{DARCY PROBLEM}  \begin{equation}\label{DARCY PROBLEM}
19  \begin{array}{rcl}  \begin{array}{rcl}
# Line 45  for all $x_{i}$. Line 45  for all $x_{i}$.
45  \subsection{Solution Method \label{DARCY SOLVE}}  \subsection{Solution Method \label{DARCY SOLVE}}
46  Unfortunate equation~\ref{DARCY PROBLEM} can not solved directly in an easy way and requires mixed FEM.    Unfortunate equation~\ref{DARCY PROBLEM} can not solved directly in an easy way and requires mixed FEM.
47  We consider a few options to solve equation~\ref{DARCY PROBLEM}  We consider a few options to solve equation~\ref{DARCY PROBLEM}
48  \subsubsection{Simple Solver}\label{SEC DARCY SIMPLE}  \subsubsection{Evaluation}\label{SEC DARCY SIMPLE}
49  The first equation of equation~\ref{DARCY PROBLEM} is inserted into the second one:  The first equation of equation~\ref{DARCY PROBLEM} is inserted into the second one:
50  \begin{equation}\label{DARCY PROBLEM SIMPLE}  \begin{equation}\label{DARCY PROBLEM SIMPLE}
51  - (\kappa_{ij} p_{,j})_{,i}  =  f  - (g_{i})_{,i}  - (\kappa_{ij} p_{,j})_{,i}  =  f  - (g_{i})_{,i}
52  \end{equation}  \end{equation}
53
54  with boundary conditions  with boundary conditions
55  \begin{equation}\label{DARCY BOUNDARY SIMPLE}  \begin{equation}\label{DARCY BOUNDARY SIMPLE}
56  \begin{array}{rcl}  \begin{array}{rcl}
# Line 57  with boundary conditions Line 58  with boundary conditions
58  p = p^{D} &  \mbox{ on } & \Gamma_{D} \\  p = p^{D} &  \mbox{ on } & \Gamma_{D} \\
59  \end{array}  \end{array}
60  \end{equation}  \end{equation}
61  Then the flux field is recovred by solving  Then the flux field is recovered by directly setting
62  \begin{equation}\label{DARCY PROBLEM SIMPLE FLUX}  \begin{equation}\label{DARCY PROBLEM SIMPLE FLUX}
63   u_{j} = g_j -  \kappa_{ij} p_{,j}     u_{j} = g_j -  \kappa_{ij} p_{,j}
64  \end{equation}  \end{equation}
65  with boundary conditions  This simple recovery process will not ensure that the (numerically) calculated flux
66  \begin{equation}\label{DARCY SIMPLE BOUNDARY FLUX}  meets the boundary conditions for flux or the incompressibility condition.
67  u_{i} = u^{N}_{i}  \mbox{ on }  \Gamma_{N}  However this is a very fast way of calculating the flux.
68  \end{equation}
69
70  \subsubsection{Global Postprocessing \label{SEC DARCY POST}}  \subsubsection{Global Postprocessing \label{SEC DARCY POST}}
71  An improved flux recovery can be achieved by solving a modified version of equation~\ref{DARCY PROBLEM SIMPLE FLUX}  An improved flux recovery can be achieved by solving a modified version of equation~\ref{DARCY PROBLEM SIMPLE FLUX}
# Line 86  u_{i} \; n_{i}  = u^{N}_{i}  \; n_{i} & Line 87  u_{i} \; n_{i}  = u^{N}_{i}  \; n_{i} &
87  u_{k,k} = f & \mbox{ on } & \Gamma_{D} \\  u_{k,k} = f & \mbox{ on } & \Gamma_{D} \\
88  \end{array}  \end{array}
89  \end{equation}    \end{equation}
90  Notice that the second condition is a natural boundray condition.  Notice that the second condition is a natural boundary condition.
91    Global post-processing is more expense than direct pressure evaluation
92  \subsubsection{Stabilization \label{SEC DARCY STAB}}  however the flux is more accurate and asymptotic incompressibility
93  Stabilization term can be added to the original problem~\ref{DARCY PROBLEM SIMPLE FLUX}, see~\cite{MasudHughes2002a}.  for mesh size towards zero can be shown, if $\omega>0$.
Firstly we subtract half of the first equation from the first eqaution
\begin{equation}\label{DARCY STAB PROBLEM A}
\kappa^{-1}_{ij} u_{j} + (p)_{,i} - \frac{1}{2} \left(  \kappa^{-1}_{ij} u_{j} + p_{,i} \right) = \kappa^{-1}_{ij} g_j -\frac{1}{2} \kappa^{-1}_{ij} g_j
\end{equation}
and then subtract half of the dievregence of the first equation from the second equation:
\begin{equation}\label{DARCY STAB PROBLEM}
\begin{array}{rcl}
\frac{1}{2} \kappa^{-1}_{ij} u_{j} + (p)_{,i} -  \frac{1}{2}  p_{,i}  & = &  \frac{1}{2} \kappa^{-1}_{ij} g_j \\
u_{i,i} - \frac{1}{2}  \left(u_{i} +   \kappa_{ij} p_{,j} \right)_{,i} & = & f - \frac{1}{2} \left(g_{i} \right)_{,i}
\end{array}
\end{equation}
We apply the boundary conditions:
\begin{equation}\label{DARCY SYM STAB PROBLEM BOUNDARY}
\begin{array}{rcl}
u_{i} \; n_{i}  = u^{N}_{i}  \; n_{i} & \mbox{ on } & \Gamma_{N} \\
p = p^{D} &  \mbox{ on } & \Gamma_{D} \\
\frac{1}{2} \left( u_{i} +   \kappa^{-1}_{ij} p_{,j} - g_{i} \right ) \; n_{i} = 0 &  \mbox{ on } & \Gamma_{N} \\
\end{array}
\end{equation}
Notice that due to the term $(p)_{,i}$ boundary condition $p = p^{D}$ acts as a natural boundary condition for the
first equation of~\ref{DARCY STAB PROBLEM} (in the form $p \; n_i = p^{D} \; n_i$ while the last two boundary conditions act as constraint and natural boundary conditions for the second equation.

\subsubsection{Symmetric Stabilization \label{SEC DARCY SYM STAB}}
Symmetric stabilization term can be added to the original problem~\ref{DARCY PROBLEM SIMPLE FLUX}, see~\cite{LoulaCorrea2006a}. In this approach divergence of the first equation
is subtracted from $2 \times$ the second equation
\begin{equation}\label{DARCY SYM STAB PROBLEM A}
- (u_{i})_{,i} + \frac{1}{2} \left( u_{i} + \kappa_{ij} p_{,j} \right)_{,i} = -  f + \frac{1}{2} \left( g_{i}\right)_{,i}
\end{equation}
The first equation is rescaled by the factor $\frac{1}{2}$ to maintain symmetry. This leads to
\begin{equation}\label{DARCY SYM STAB PROBLEM}
\begin{array}{rcl}
\frac{1}{2}   \kappa^{-1}_{ij} u_{j} + \frac{1}{2}  p_{,i}& = &  \frac{1}{2}  \kappa^{-1}_{ij} g_{j} \\
\frac{1}{2}  \left(  \kappa_{ij} p_{,j} - u_{i} \right)_{,i} & = & -  f + \frac{1}{2}  \left(g_{i}\right)_{,i}
\end{array}
\end{equation}
We add the boundary conditions
\begin{equation}\label{DARCY SYM STAB PROBLEM BOUNDARY}
\begin{array}{rcl}
u_{i} \; n_{i}  = u^{N}_{i}  \; n_{i} & \mbox{ on } & \Gamma_{N} \\
p = p^{D} &  \mbox{ on } & \Gamma_{D} \\
\frac{1}{2} \left( \kappa_{ij} p_{,j} - g_{i} - u_{i} \right ) \; n_{i}
=  - u^{N}_{i}  \; n_{i} & \mbox{ on } & \Gamma_{N} \\
\end{array}
\end{equation}
Notice that the last contion provides natural boundary conditons for the last eqaution of~\ref{DARCY SYM STAB PROBLEM}.

94
95
96  \subsection{Functions}  \subsection{Functions}
97  \begin{classdesc}{DarcyFlow}{domain, \optional{w=1., \optional{solver=\member{DarcyFlow.SYMSTAB},  \optional{  \begin{classdesc}{DarcyFlow}{domain, \optional{w=1., \optional{solver=\member{DarcyFlow.POST},  \optional{
98  useReduced=\True, \optional{ verbose=\True} } }}}  useReduced=\True, \optional{ verbose=\True} } }}}
99  opens the Darcy flux problem\index{Darcy flux} on the \Domain domain.  opens the Darcy flux problem\index{Darcy flux} on the \Domain domain.
100  Reduced approximations for pressure and flux are used if \var{useReduced} is set.  Reduced approximations for pressure and flux are used if \var{useReduced} is set.
101  Argument \var{solver} defines the solver method.  Argument \var{solver} defines the solver method.
102  If \var{verbose} is set some information are printed.  If \var{verbose} is set some information are printed.
103  \var{w} defines the weighting factor $\omega$ for global postprocessing pf the flux (see equation~\ref{DARCY PROBLEM POST FLUX A}.)  \var{w} defines the weighting factor $\omega$ for global post-processing of the flux (see equation~\ref{DARCY PROBLEM POST FLUX A}.)
104  \end{classdesc}  \end{classdesc}
105
106  \begin{memberdesc}[DarcyFlow]{SIMPLE}  \begin{memberdesc}[DarcyFlow]{EVAL}
107  simple solver, see section~\ref{SEC DARCY SIMPLE}.  flux is calculated directly from pressure evaluation, see section~\ref{SEC DARCY SIMPLE}.
108  \end{memberdesc}  \end{memberdesc}
109
110  \begin{memberdesc}[DarcyFlow]{POST}  \begin{memberdesc}[DarcyFlow]{SMOOTH}
111  solver using global postprocessing of flux, see section~\ref{SEC DARCY POST}.  solver using global post-processing of flux with weighting factor $\omega=0$, see section~\ref{SEC DARCY POST}.
112  \end{memberdesc}  \end{memberdesc}
113
114  \begin{memberdesc}[DarcyFlow]{STAB}  \begin{memberdesc}[DarcyFlow]{POST}
115   solver uses (non-symmetric) stabilization, see section~\ref{SEC DARCY STAB}.  solver using global post-processing of flux, see section~\ref{SEC DARCY POST}.
\end{memberdesc}

\begin{memberdesc}[DarcyFlow]{SYMSTAB}
solver uses symmetric stabilization, see section~\ref{SEC DARCY SYM STAB}.
116  \end{memberdesc}  \end{memberdesc}
117
118  \begin{methoddesc}[DarcyFlow]{setValue}{\optional{f=None, \optional{g=None, \optional{location_of_fixed_pressure=None, \optional{location_of_fixed_flux=None,  \begin{methoddesc}[DarcyFlow]{setValue}{\optional{f=None, \optional{g=None, \optional{location_of_fixed_pressure=None, \optional{location_of_fixed_flux=None,
# Line 169  solver using global postprocessing of fl Line 120  solver using global postprocessing of fl
120  assigns values to the model parameters. Values can be assigned using various  assigns values to the model parameters. Values can be assigned using various
121  calls -- in particular in a time dependent problem only values that change  calls -- in particular in a time dependent problem only values that change
122  over time need to be reset. The permeability can be defined as a scalar  over time need to be reset. The permeability can be defined as a scalar
123  (isotropic), or a symmetrix matrix (anisotropic).  (isotropic), or a symmetric matrix (anisotropic).
124  \var{f} and \var{g} are the corresponding parameters in~\ref{DARCY PROBLEM}.  \var{f} and \var{g} are the corresponding parameters in~\ref{DARCY PROBLEM}.
125  The locations and components where the flux is prescribed are set by positive  The locations and components where the flux is prescribed are set by positive
126  values in \var{location_of_fixed_flux}.  values in \var{location_of_fixed_flux}.
# Line 182  where the pressure is prescribed. Line 133  where the pressure is prescribed.
133  The method will try to cast the given values to appropriate \Data class objects.  The method will try to cast the given values to appropriate \Data class objects.
134  \end{methoddesc}  \end{methoddesc}
135
\begin{methoddesc}[DarcyFlow]{setTolerance}{\optional{rtol=1e-4}}
sets the relative tolerance \var{rtol} for the pressure in the stabilized solvers.
\end{methoddesc}

136  \begin{methoddesc}[DarcyFlow]{getSolverOptionsFlux}{}  \begin{methoddesc}[DarcyFlow]{getSolverOptionsFlux}{}
137  returns the solver options used to solve the flux problems.  returns the solver options used to solve the flux problems.
138  Use this \SolverOptions object to control the solution algorithms.  Use this \SolverOptions object to control the solution algorithms.
139    This option is only relevant if global postprocesing is used.
140  \end{methoddesc}  \end{methoddesc}
141
142  \begin{methoddesc}[DarcyFlow]{getSolverOptionsPressure}{}  \begin{methoddesc}[DarcyFlow]{getSolverOptionsPressure}{}
# Line 198  problems. Line 145  problems.
145  Use this object to control the solution algorithms.  Use this object to control the solution algorithms.
146  \end{methoddesc}  \end{methoddesc}
147
148  \begin{methoddesc}[DarcyFlow]{solve}{u0,p0, \optional{ iter_restart=20, \optional{max_iter=100}}}  \begin{methoddesc}[DarcyFlow]{solve}{u0,p0}
149  solves the problem and returns approximations for the flux $v$ and the pressure $p$.  solves the problem and returns approximations for the flux $v$ and the pressure $p$.
150  \var{u0} and \var{p0} define initial guesses for flux and pressure.  \var{u0} and \var{p0} define initial guesses for flux and pressure.
151  Values marked by positive values \var{location_of_fixed_flux} and  Values marked by positive values \var{location_of_fixed_flux} and
152  \var{location_of_fixed_pressure}, respectively, are kept unchanged.  \var{location_of_fixed_pressure}, respectively, are kept unchanged.
\var{max_iter} sets the maximum number of outer iteration steps allowed for solving the stabilzed
problems. \var{iter_restart} sets the number of iteration after which the
iteration is restarted\footnote{the larger \var{iter_restart} the more memory is required. The smaller
\var{iter_restart} the more iteration steps are required.}.
153  \end{methoddesc}  \end{methoddesc}
154
155  \begin{methoddesc}[DarcyFlow]{getFlux}{p, \optional{ u0 = None}}  \begin{methoddesc}[DarcyFlow]{getFlux}{p, \optional{ u0 = None}}
# Line 215  on locations where \var{location_of_fixe Line 158  on locations where \var{location_of_fixe
158  Notice that \var{g} and \var{f} are used.  Notice that \var{g} and \var{f} are used.
159  \end{methoddesc}  \end{methoddesc}
160
161    \begin{figure}
162    \centerline{\includegraphics[width=\figwidth]{darcy_result}}
163  %\subsection{Example: Gravity Flow}  \caption{Flux and pressure field of the Dary flow example.}
164  %later  \label{DIFFUSION FIG 1}
165    \end{figure}
166
167    \subsection{Example: Gravity Flow}
168    The following script \file{darcy.py}\index{scripts!\file{darcy.py.py}}\index{Darcy flow}
169    which is available in the \ExampleDirectory illustrates the usage of the
170    \class{DarcyFlow} class:
171    \begin{python}
172       from esys.escript import *
173       from esys.escript.models import DarcyFlow
174       from esys.finley import Rectangle
175       from esys.weipa import saveVTK
176       mydomain = Rectangle(l0=2.,l1=1.,n0=40, n1=20)
177       x = mydomain.getX()
178       p_BC=whereZero(x-1.)*wherePositive(x-1.)
179       u_BC=(whereZero(x)+whereZero(x-2.)) * [1.,0.] + \
180         (whereZero(x) + whereZero(x-1.)*whereNonPositive(x-1.0)) * [0., 1.]
181       mypde = DarcyFlow(domain=mydomain)
182       mypde.setValue(g=[0., 2],
183                      location_of_fixed_pressure=p_BC,
184                      location_of_fixed_flux=u_BC,
185                      permeability=100.)
186
187       u,p=mypde.solve(u0=x*[0., -1.], p0=0)
188       saveVTK("u.vtu",flux=u, pressure=p)
189    \end{python}
190    In the example the pressure is fixed to the initial presure \var{p0} on the right half of the top face.
191    The normal flux is set on all other faces. The corresponding values for the flux are set by the initial value
192    \var{u0}.

Legend:
 Removed from v.3502 changed lines Added in v.3911