1 

2 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
3 
% 
4 
% Copyright (c) 20032010 by University of Queensland 
5 
% Earth Systems Science Computational Center (ESSCC) 
6 
% http://www.uq.edu.au/esscc 
7 
% 
8 
% Primary Business: Queensland, Australia 
9 
% Licensed under the Open Software License version 3.0 
10 
% http://www.opensource.org/licenses/osl3.0.php 
11 
% 
12 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
13 

14 
\section{Darcy Flux} 
15 
\label{DARCY FLUX} 
16 
We want to calculate the velocity $u$ and pressure $p$ on a domain $\Omega$ 
17 
solving the Darcy flux problem\index{Darcy flux}\index{Darcy flow} 
18 
\begin{equation}\label{DARCY PROBLEM} 
19 
\begin{array}{rcl} 
20 
u_{i} + \kappa_{ij} p_{,j} & = & g_{i} \\ 
21 
u_{k,k} & = & f 
22 
\end{array} 
23 
\end{equation} 
24 
with the boundary conditions 
25 
\begin{equation}\label{DARCY BOUNDARY} 
26 
\begin{array}{rcl} 
27 
u_{i} \; n_{i} = u^{N}_{i} \; n_{i} & \mbox{ on } & \Gamma_{N} \\ 
28 
p = p^{D} & \mbox{ on } & \Gamma_{D} \\ 
29 
\end{array} 
30 
\end{equation} 
31 
where $\Gamma_{N}$ and $\Gamma_{D}$ are a partition of the boundary of 
32 
$\Omega$ with $\Gamma_{D}$ nonempty, $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 
\begin{equation} 
40 
\alpha_{0} \; x_{i} x_{i} \le \kappa_{ij} x_{i} x_{j} \le \alpha_{1} \; x_{i} x_{i} 
41 
\end{equation} 
42 
for all $x_{i}$. 
43 

44 

45 
\subsection{Solution Method \label{DARCY SOLVE}} 
46 
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} 
48 
\subsubsection{Evaluation}\label{SEC DARCY SIMPLE} 
49 
The first equation of equation~\ref{DARCY PROBLEM} is inserted into the second one: 
50 
\begin{equation}\label{DARCY PROBLEM SIMPLE} 
51 
 (\kappa_{ij} p_{,j})_{,i} = f  (g_{i})_{,i} 
52 
\end{equation} 
53 

54 
with boundary conditions 
55 
\begin{equation}\label{DARCY BOUNDARY SIMPLE} 
56 
\begin{array}{rcl} 
57 
\kappa_{ij} p_{,j} \; n_{i} = ( g_{i}  u^{N}_{i} ) \; n_{i} & \mbox{ on } & \Gamma_{N} \\ 
58 
p = p^{D} & \mbox{ on } & \Gamma_{D} \\ 
59 
\end{array} 
60 
\end{equation} 
61 
Then the flux field is recovered by directly setting 
62 
\begin{equation}\label{DARCY PROBLEM SIMPLE FLUX} 
63 
u_{j} = g_j  \kappa_{ij} p_{,j} 
64 
\end{equation} 
65 
This simple recovery process will not ensure that the (numerically) calculated flux 
66 
meets the boundary conditions for flux or the incompressibility condition. 
67 
However this is a very fast way of calculating the flux. 
68 

69 

70 
\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} 
72 
adding the the gradient of the divergence of the flux: 
73 
\begin{equation}\label{DARCY PROBLEM POST FLUX} 
74 
\kappa^{1}_{ij} u_{j}  
75 
(\lambda \cdot u_{k,k} )_{,i}= 
76 
\kappa^{1}_{ij} g_j p_{,i} 
77 
 (\lambda \cdot f )_{,i} 
78 
\end{equation} 
79 
where 
80 
\begin{equation}\label{DARCY PROBLEM POST FLUX A} 
81 
\lambda = \omega \cdot \kappa^{1} \cdot vol(\Omega)^{1/d} \cdot h 
82 
\end{equation} 
83 
with a nonnegative factor $\omega$, $d$ is the spatial dimension and $h$ is the local element size. 
84 
\begin{equation}\label{DARCY PROBLEM POST FLUX BOUNDARY} 
85 
\begin{array}{rcl} 
86 
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} 
89 
\end{equation} 
90 
Notice that the second condition is a natural boundary condition. 
91 
Global postprocessing is more expense than direct pressure evaluation 
92 
however the flux is more accurate and asymptotic incompressibility 
93 
for mesh size towards zero can be shown, if $\omega>0$. 
94 

95 

96 
\subsection{Functions} 
97 
\begin{classdesc}{DarcyFlow}{domain, \optional{w=1., \optional{solver=\member{DarcyFlow.EVAL}, \optional{ 
98 
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 postprocessing of the flux (see equation~\ref{DARCY PROBLEM POST FLUX A}.) 
104 
\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 postprocessing 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 postprocessing 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, 
119 
\\\optional{permeability=None}}}}}} 
120 
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 
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}. 
125 
The locations and components where the flux is prescribed are set by positive 
126 
values in \var{location_of_fixed_flux}. 
127 
The locations where the pressure is prescribed are set by by positive values 
128 
of \var{location_of_fixed_pressure}. 
129 
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 
131 
normal component of the flux must be defined. There must be at least one point 
132 
where the pressure is prescribed. 
133 
The method will try to cast the given values to appropriate \Data class objects. 
134 
\end{methoddesc} 
135 

136 
\begin{methoddesc}[DarcyFlow]{getSolverOptionsFlux}{} 
137 
returns the solver options used to solve the flux problems. 
138 
Use this \SolverOptions object to control the solution algorithms. 
139 
This option is only relevant if global postprocesing is used. 
140 
\end{methoddesc} 
141 

142 
\begin{methoddesc}[DarcyFlow]{getSolverOptionsPressure}{} 
143 
returns a \SolverOptions object with the options used to solve the pressure 
144 
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} 
154 

155 
\begin{methoddesc}[DarcyFlow]{getFlux}{p, \optional{ u0 = None}} 
156 
returns the flux for a given pressure \var{p} where the flux is equal to \var{u0} 
157 
on locations where \var{location_of_fixed_flux} is positive, see \member{setValue}. 
158 
Notice that \var{g} and \var{f} are used. 
159 
\end{methoddesc} 
160 

161 

162 

163 
%\subsection{Example: Gravity Flow} 
164 
%later 
165 
