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

revision 3301 by jfenwick, Mon Oct 25 01:27:54 2010 UTC revision 3911 by jfenwick, Thu Jun 14 01:01:03 2012 UTC
# Line 1  Line 1
1
2    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3    %
4    % Copyright (c) 2003-2012 by University of Queensland
5    % Earth Systems Science Computational Center (ESSCC)
6    % http://www.uq.edu.au/esscc
7    %
8    % Primary Business: Queensland, Australia
11    %
12    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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$ solving  We want to calculate the flux $u$ and pressure $p$ on a domain $\Omega$
17  the Darcy flux problem \index{Darcy flux}\index{Darcy flow}  solving the Darcy flux problem\index{Darcy flux}\index{Darcy flow}
18  \label{DARCY PROBLEM}  \label{DARCY PROBLEM}
19  \begin{array}{rcl}  \begin{array}{rcl}
20  u_{i} + \kappa_{ij} p_{,j} & = & g_{i} \\  u_{i} + \kappa_{ij} p_{,j} & = & g_{i} \\
# Line 15  u_{i} \; n_{i}  = u^{N}_{i}  \; n_{i} & Line 28  u_{i} \; n_{i}  = u^{N}_{i}  \; n_{i} &
28  p = p^{D} &  \mbox{ on } & \Gamma_{D} \\  p = p^{D} &  \mbox{ on } & \Gamma_{D} \\
29  \end{array}  \end{array}
30
31  where $\Gamma_{N}$ and $\Gamma_{D}$ are a partition of the boundary of $\Omega$ with $\Gamma_{D}$ non empty, $n_{i}$ is the outer normal field of the boundary of $\Omega$, $u^{N}_{i}$ and $p^{D}$ are given functions on $\Omega$, $g_{i}$ and $f$ are given source terms and $\kappa_{ij}$ is the given permeability. We assume that $\kappa_{ij}$ is symmetric (which is not really required) and positive definite, i.e there are positive constants $\alpha_{0}$ and $\alpha_{1}$ which are independent from the location in $\Omega$ such that  where $\Gamma_{N}$ and $\Gamma_{D}$ are a partition of the boundary of
32    $\Omega$ with $\Gamma_{D}$ non-empty, $n_{i}$ is the outer normal field of the
33    boundary of $\Omega$, $u^{N}_{i}$ and $p^{D}$ are given functions on $\Omega$,
34    $g_{i}$ and $f$ are given source terms and $\kappa_{ij}$ is the given
35    permeability.
36    We assume that $\kappa_{ij}$ is symmetric (which is not really required) and
37    positive definite, i.e. there are positive constants $\alpha_{0}$ and
38    $\alpha_{1}$ which are independent from the location in $\Omega$ such that
39
40  \alpha_{0} \; x_{i} x_{i} \le \kappa_{ij} x_{i} x_{j} \le \alpha_{1} \; x_{i} x_{i}  \alpha_{0} \; x_{i} x_{i} \le \kappa_{ij} x_{i} x_{j} \le \alpha_{1} \; x_{i} x_{i}
41
42  for all $x_{i}$.  for all $x_{i}$.
43
44
45  \subsection{Solution Method \label{DARCY SOLVE}}  \subsection{Solution Method \label{DARCY SOLVE}}
46  It is useful to write equation~\ref{DARCY PROBLEM} in operator form. For any pressure $p$  Unfortunate equation~\ref{DARCY PROBLEM} can not solved directly in an easy way and requires mixed FEM.
47  we set  We consider a few options to solve equation~\ref{DARCY PROBLEM}
48    \subsubsection{Evaluation}\label{SEC DARCY SIMPLE}
49  (Gp)_{i} =  p_{,j}  The first equation of equation~\ref{DARCY PROBLEM} is inserted into the second one:
50    \label{DARCY PROBLEM SIMPLE}
51  and a velocity $v$ we set  - (\kappa_{ij} p_{,j})_{,i}  =  f  - (g_{i})_{,i}

Dv = v_{k,k}

We $K=(\kappa_{ij})$ we can write the Darcy problem~\ref{DARCY PROBLEM} as

\begin{array}{rcl}
u + K \, Gp & = & g \\
Du & = & f
\end{array}

We solve this equation by minimising the functional
\label{DARCY COST}
J(u,p):=\|K^{-\frac{1}{2}}(u + K \, G p - g)\|^2_{0} +  \|\lambda (Du-f) \|_{0}^2
52
53  over all suitable $u$ and $p$. In this equation we set $\|p\|^2_{0}=(p,p)_{0}$ with
54    with boundary conditions
55  (p,q)_{0} = \int_{\Omega } p\cdot q \, dx  \label{DARCY BOUNDARY SIMPLE}

The factor $\lambda>0$ is a weighting factor.
A simple calculation shows that
one has to solve

( K^{-1} (v + K \, Gq) , u +K \, G p - g)_{0} +  (\lambda Dv,\lambda (Du-f) )_{0} =0

for all velocities $v$ and pressure $q$ which fulfill the homogeneous boundary conditions~\ref{DARCY BOUNDARY}.
This so-called variational equation can be translated back into operator notation

56  \begin{array}{rcl}  \begin{array}{rcl}
57  (K^{-1}+ D^*\lambda^2 D)u + Gp & = &  D^*\lambda f + K^{-1} g \\  \kappa_{ij} p_{,j} \; n_{i}  = ( g_{i} - u^{N}_{i} )  \; n_{i} & \mbox{ on } & \Gamma_{N} \\
58  G^*u  + G^*K \, G p & = & G^*g \\  p = p^{D} &  \mbox{ on } & \Gamma_{D} \\
59  \end{array}  \end{array}
60
61  where $D^*$ and $Q^*$ denote the  adjoint operators with respect to $(.,.)_{0}$.  Then the flux field is recovered by directly setting
62  In~\cite{LEASTSQUARESFEM1994} it has been shown that this problem is continuous and coercive and therefore has a unique solution. Also standard FEM methods can be used for discretization. It is also possible  \label{DARCY PROBLEM SIMPLE FLUX}
63  to solve the problem in coupled form, however this approach leads in some cases to a very ill-conditioned stiffness matrix in particular in the case of a very small or large permeability ($\alpha_{1} \ll 1$ or $\alpha_{0} \gg 1$)     u_{j} = g_j -  \kappa_{ij} p_{,j}
64
65  The approach we are taking is to eliminate the velocity $v$ from the problem. Assuming that $p$ is known we have  This simple recovery process will not ensure that the (numerically) calculated flux
66  \label{DARCY V FORM}  meets the boundary conditions for flux or the incompressibility condition.
67  v= (K^{-1}+ D^*\lambda^2 D)^{-1} ( D^*\lambda f + K^{-1} g - Gp)  However this is a very fast way of calculating the flux.
68
69  (notice that $K^{-1}+\lambda D^*D$ is coercive) which is inserted into the second equation
70    \subsubsection{Global Postprocessing \label{SEC DARCY POST}}
71  G^* (K^{-1}+\lambda D^*D)^{-1} (\lambda D^*f + K^{-1} g - Gp) + G^* KG p = G^*g  An improved flux recovery can be achieved by solving a modified version of equation~\ref{DARCY PROBLEM SIMPLE FLUX}
73  which is  \label{DARCY PROBLEM POST FLUX}
74    \kappa^{-1}_{ij} u_{j} -
75  G^* ( K - (K^{-1}+ D^*\lambda^2 D)^{-1} ) G p = G^* (g-(K^{-1}+D^*\lambda^2 D)^{-1} ( D^*\lambda f + K^{-1} g) )  (\lambda \cdot u_{k,k} )_{,i}=
76    \kappa^{-1}_{ij} g_j- p_{,i}
77    - (\lambda \cdot f )_{,i}
78
79    where
80    \label{DARCY PROBLEM POST FLUX A}
81    \lambda = \omega \cdot |\kappa^{-1}| \cdot vol(\Omega)^{1/d} \cdot h
82
83  We use the PCG \index{linear solver!PCG}\index{PCG} method to solve this.  with a non-negative factor $\omega$, $d$ is the spatial dimension and $h$ is the local element size.
84  The residual $r$  is given as  \label{DARCY PROBLEM POST FLUX BOUNDARY}

85  \begin{array}{rcl}  \begin{array}{rcl}
86  r & =&  G^* \left( g - K\, G p - v \right)  u_{i} \; n_{i}  = u^{N}_{i}  \; n_{i} & \mbox{ on } & \Gamma_{N} \\
87    u_{k,k} = f & \mbox{ on } & \Gamma_{D} \\
88  \end{array}  \end{array}
89
90  for the current pressure approximation $p$ and current velocity $v$ defined by    Notice that the second condition is a natural boundary condition.
91  equation~\ref{DARCY V FORM}.  Global post-processing is more expense than direct pressure evaluation
92  So in a particular implementation we use $\hat{r}=g-K\, Gp-v$ to represent the residual.  however the flux is more accurate and asymptotic incompressibility
93  The evaluation of the iteration operator for a given $p$ is then  for mesh size towards zero can be shown, if $\omega>0$.
returning $Qp+v$ where $v$ is the solution of
\label{UPDATE W}
(K^{-1}+ D^*\lambda^2 D)v = Gp

To derive a preconditioner we use the identity

\begin{array}{rcl}
94
G^* ( K - (K^{-1}+ D^*\lambda^2 D)^{-1} ) G  & = & G^* (I - (I + K D^*\lambda^2 D)^{-1}) K G \\
& \approx &  G^* (K D^*\lambda^2 D) K G  \\
& =  & G^*  K ( D^*  \lambda^2 D) K G \\
& \approx  & G^* \frac{\lambda^2}{dx^2} K^2 G
\end{array}

where $dx$ is the local mesh size and we use the approximation

D^*  \lambda^2 D \approx \frac{\lambda^2}{dx^2}

Therefore we use $G^* \frac{\lambda^2}{dx^2} K^2 G$ as a preconditioner.To evalute the preconditioner
we need to solve the equation
\label{UPDATE P}
G^* \frac{\lambda^2}{dx^2} K^2 G p =  G^* \hat{r}

It remains to answer the question how to choose $\lambda$. We need to balance the first and second
term in $J(u,p)$ in equation~\ref{DARCY COST}. We inspect $J$ for
$(\hat{u}, \hat{p})$ which is a perturbed exact solution $(u,p)$.
Assuming $\hat{u}=u+u_{0}e^{ik^tx}$
and $\hat{p}=p+p_{0}e^{ik^tx}$ and constant $K$ we get

J(\hat{u},\hat{p}) = C \left[ ( \|K^{-1}\|_{2} |u_{0}|^2 + \|K\|_{2} \|k\|_{2}^2 |p_{0}|^2| )
+ \lambda^2 \|k\|_{2}^2 |u_{0}|^2  \right]

with some constant $C>0$. The first two terms and the third term correspond to the first term and
second term in the definition of $J(u,p)$ in equation~\ref{DARCY COST}. For small $\|k\|_{2}$
(i.e. for a smooth perturbation) $J(\hat{u},\hat{p})$ is dominated by $\|K^{-1}\|_{2} |u_{0}|^2$.
To scale the second term which is corresponds to the incompressibility condition for the velocity
we need to meet the condition $\|K^{-1}\|_{2} = \lambda^2 \|k\|_{2}^2$.
Taking the boundary conditions into consideration the smallest possible value for $\|k\|_{2} = \frac{\pi}{l}$ where $l$ is the longest edge of the domain. This leads to the
\label{DARCY LAMBDA}
\lambda = \|K^{-1}\|_{2}^{\frac{1}{2}} \frac{l}{\pi}

Notice that with this setting the preconditioner $G^* \frac{\lambda^2}{dx^2} K^2 G$ becomes
equivalent to $G^* K G$ if $K$ is a diagonal matrix and the mesh has a constant mesh size.

The residual norm used in the PCG is given as
\label{DARCY R NORM}
\|r\|_{PCG}^2 = \int r (G^* \frac{\lambda^2}{dx^2} K^2 G)^{-1} r \; dx =\int \hat{r}  G ( G^* \frac{\lambda^2}{dx^2} K^2 G)^{-1}  G^* \hat{r} \; dx \approx
\int \hat{r} K^{-1}  \hat{r} \; dx

The iteration is terminated if
\label{DARCY STOP}
\|r\|_{PCG} \le \mbox{ATOL}

where we set
\label{DARCY ATOL DEF}
\mbox{ATOL} = \mbox{atol} + \mbox{rtol} \cdot \left(\frac{1}{\|K^{-\frac{1}{2}}v\|_{0}} + \frac{1}{\|K^{\frac{1}{2}} G p\|_{0}} \right)^{-1}

where rtol is a given relative tolerance and $\mbox{atol}$ is a given absolute tolerance (typically $=0$).
Notice that if $Gp$ and $v$ both are zero, the pair $(0,p)$ is a solution.
The problem is that ATOL is depending on the solution $p$ and $v$ calculated form~\ref{DARCY V FORM}.
In practice one use the initial guess for $p$
to get a first value for ATOL. If the stopping criterion is met in the PCG iteration, a new $v$ is calculated from the current pressure approximation and ATOL is recalculated. If \ref{DARCY STOP} is still fulfilled the calculation is terminated and $(v,p)$ is returned. Otherwise PCG is restarted with a new ATOL.
95
96  \subsection{Functions}  \subsection{Functions}
97  \begin{classdesc}{DarcyFlow}{domain}  \begin{classdesc}{DarcyFlow}{domain, \optional{w=1., \optional{solver=\member{DarcyFlow.POST},  \optional{
98  opens the Darcy flux problem\index{Darcy flux} on the \Domain domain.  useReduced=\True, \optional{ verbose=\True} } }}}
99    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.
101    Argument \var{solver} defines the solver method.
102    If \var{verbose} is set some information are printed.
103    \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]{EVAL}
107    flux is calculated directly from pressure evaluation, see section~\ref{SEC DARCY SIMPLE}.
108    \end{memberdesc}
109
110    \begin{memberdesc}[DarcyFlow]{SMOOTH}
111    solver using global post-processing of flux with weighting factor $\omega=0$, see section~\ref{SEC DARCY POST}.
112    \end{memberdesc}
113
114    \begin{memberdesc}[DarcyFlow]{POST}
115    solver using global post-processing of flux, see section~\ref{SEC DARCY POST}.
116    \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,
119  \\\optional{permeability=None}}}}}}  \\\optional{permeability=None}}}}}}
120  assigns values to the model parameters. Values can be assigned using various calls - in particular  assigns values to the model parameters. Values can be assigned using various
121  in a time dependent problem only values that change over time needs to be reset. The permeability can be defined as scalar (isotropic), a vector (orthotropic) or a matrix (anisotropic).  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
123    (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 values in  The locations and components where the flux is prescribed are set by positive
126  \var{location_of_fixed_flux}.  values in \var{location_of_fixed_flux}.
127  The locations where the pressure is prescribed are set by  The locations where the pressure is prescribed are set by by positive values
128  by positive values of \var{location_of_fixed_pressure}.  of \var{location_of_fixed_pressure}.
129  The values of the pressure and flux are defined by the initial guess.  The values of the pressure and flux are defined by the initial guess.
130  Notice that at any point on the boundary of the domain the pressure or the normal component of  Notice that at any point on the boundary of the domain the pressure or the
131  the flux must be defined. There must be at least one point where the pressure is prescribed.  normal component of the flux must be defined. There must be at least one point
132  The method will try to cast the given values to appropriate  where the pressure is prescribed.
133  \Data class objects.  The method will try to cast the given values to appropriate \Data class objects.
\end{methoddesc}

\begin{methoddesc}[DarcyFlow]{setTolerance}{\optional{rtol=1e-4}}
sets the relative tolerance \mbox{rtol} in \ref{DARCY ATOL DEF}.
\end{methoddesc}

\begin{methoddesc}[DarcyFlow]{setAbsoluteTolerance}{\optional{atol=0.}}
sets the absolute tolerance \mbox{atol} in \ref{DARCY ATOL DEF}.
134  \end{methoddesc}  \end{methoddesc}
135
136  \begin{methoddesc}[DarcyFlow]{getSolverOptionsFlux}{}  \begin{methoddesc}[DarcyFlow]{getSolverOptionsFlux}{}
137  Returns the solver options used to solve the flux problems~(\ref{DARCY V FORM}) and~(\ref{UPDATE W}). Use this  returns the solver options used to solve the flux problems.
138   \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}{}
143  Returns the solver options used to solve the pressure problems~(\ref{UPDATE P}) as a preconditioner.  returns a \SolverOptions object with the options used to solve the pressure
144  Use this \SolverOptions object to control the solution algorithms.  problems.
145    Use this object to control the solution algorithms.
146    \end{methoddesc}
147
148    \begin{methoddesc}[DarcyFlow]{solve}{u0,p0}
149    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.
151    Values marked by positive values \var{location_of_fixed_flux} and
152    \var{location_of_fixed_pressure}, respectively, are kept unchanged.
153  \end{methoddesc}  \end{methoddesc}
154
155  \begin{methoddesc}[DarcyFlow]{solve}{u0,p0, \optional{max_iter=100, \optional{verbose=False}}}  \begin{methoddesc}[DarcyFlow]{getFlux}{p, \optional{ u0 = None}}
156  solves the problem. and returns approximations for the flux $v$ and the pressure $p$.  returns the flux for a given pressure \var{p} where the flux is equal to \var{u0}
157  \var{u0} and \var{p0} define initial guess for flux and pressure. Values marked  on locations where \var{location_of_fixed_flux} is positive,  see \member{setValue}.
158  by positive values \var{location_of_fixed_flux} and \var{location_of_fixed_pressure}, respectively, are kept unchanged. \var{max_iter} sets the maximum number of iterations steps allowed for solving the coupled problem.  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    \caption{Flux and pressure field of the Dary flow example.}
164    \label{DIFFUSION FIG 1}
165    \end{figure}
166
167  \subsection{Example: Gravity Flow}  \subsection{Example: Gravity Flow}
later
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]-1.)*wherePositive(x[0]-1.)
179       u_BC=(whereZero(x[0])+whereZero(x[0]-2.)) * [1.,0.] + \
180         (whereZero(x[1]) + whereZero(x[1]-1.)*whereNonPositive(x[0]-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[1]*[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.3301 changed lines Added in v.3911