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 
\subsection{Solution Method \label{DARCY SOLVE}} 
45 
It is useful to write \eqn{DARCY PROBLEM} in operator form. 
46 
For any pressure $p$ we set 
47 
\begin{equation} 
48 
(Gp)_{i} = p_{,j} 
49 
\end{equation} 
50 
and for velocity $v$ we set 
51 
\begin{equation} 
52 
Dv = v_{k,k} 
53 
\end{equation} 
54 
With $K=(\kappa_{ij})$ we can write the Darcy problem~\ref{DARCY PROBLEM} as 
55 
\begin{equation} 
56 
\begin{array}{rcl} 
57 
u + K \, Gp & = & g \\ 
58 
Du & = & f 
59 
\end{array} 
60 
\end{equation} 
61 
We solve this equation by minimising the functional 
62 
\begin{equation}\label{DARCY COST} 
63 
J(u,p):=\K^{\frac{1}{2}}(u + K \, G p  g)\^2_{0} + \\lambda (Duf) \_{0}^2 
64 
\end{equation} 
65 
over all suitable $u$ and $p$. 
66 
In this equation we set $\p\^2_{0}=(p,p)_{0}$ with 
67 
\begin{equation} 
68 
(p,q)_{0} = \int_{\Omega } p\cdot q \, dx 
69 
\end{equation} 
70 
The factor $\lambda>0$ is a weighting factor. 
71 
A simple calculation shows that one has to solve 
72 
\begin{equation} 
73 
( K^{1} (v + K \, Gq) , u +K \, G p  g)_{0} + (\lambda Dv,\lambda (Duf) )_{0} =0 
74 
\end{equation} 
75 
for all velocities $v$ and pressure $q$ which fulfill the homogeneous boundary 
76 
conditions~\ref{DARCY BOUNDARY}. 
77 
This socalled variational equation can be translated back into operator notation 
78 
\begin{equation} 
79 
\begin{array}{rcl} 
80 
(K^{1}+ D^*\lambda^2 D)u + Gp & = & D^*\lambda f + K^{1} g \\ 
81 
G^*u + G^*K \, G p & = & G^*g \\ 
82 
\end{array} 
83 
\end{equation} 
84 
where $D^*$ and $Q^*$ denote the adjoint operators with respect to $(.,.)_{0}$. 
85 
In~\cite{LEASTSQUARESFEM1994} it has been shown that this problem is 
86 
continuous and coercive and therefore has a unique solution. 
87 
Standard FEM methods can be used for discretization. It is also possible to 
88 
solve the problem in coupled form, however this approach leads in some cases 
89 
to a very illconditioned stiffness matrix, in particular in the case of a 
90 
very small or large permeability ($\alpha_{1} \ll 1$ or $\alpha_{0} \gg 1$). 
91 

92 
The approach we are taking is to eliminate the velocity $v$ from the problem. 
93 
Assuming that $p$ is known we have\footnote{notice 
94 
that $K^{1}+\lambda D^*D$ is coercive} 
95 
\begin{equation}\label{DARCY V FORM} 
96 
v= (K^{1}+ D^*\lambda^2 D)^{1} ( D^*\lambda f + K^{1} g  Gp) 
97 
\end{equation} 
98 
which is inserted into the second equation 
99 
\begin{equation} 
100 
G^* (K^{1}+\lambda D^*D)^{1} (\lambda D^*f + K^{1} g  Gp) + G^* KG p = G^*g 
101 
\end{equation} 
102 
which is 
103 
\begin{equation} 
104 
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) ) 
105 
\end{equation} 
106 
We use the PCG\index{linear solver!PCG}\index{PCG} method to solve this. 
107 
The residual $r$ is given as 
108 
\begin{equation} 
109 
\begin{array}{rcl} 
110 
r & =& G^* \left( g  K\, G p  v \right) 
111 
\end{array} 
112 
\end{equation} 
113 
for the current pressure approximation $p$ and current velocity $v$ defined by 
114 
\eqn{DARCY V FORM}. 
115 
So in a particular implementation we use $\hat{r}=gK\, Gpv$ to represent the residual. 
116 
The evaluation of the iteration operator for a given $p$ is then returning 
117 
$Qp+v$ where $v$ is the solution of 
118 
\begin{equation}\label{UPDATE W} 
119 
(K^{1}+ D^*\lambda^2 D)v = Gp 
120 
\end{equation} 
121 
To derive a preconditioner we use the identity 
122 
\begin{equation} 
123 
\begin{array}{rcl} 
124 

125 
G^* ( K  (K^{1}+ D^*\lambda^2 D)^{1} ) G & = & G^* (I  (I + K D^*\lambda^2 D)^{1}) K G \\ 
126 
& \approx & G^* (K D^*\lambda^2 D) K G \\ 
127 
& = & G^* K ( D^* \lambda^2 D) K G \\ 
128 
& \approx & G^* \frac{\lambda^2}{dx^2} K^2 G 
129 
\end{array} 
130 
\end{equation} 
131 
where $dx$ is the local mesh size and we use the approximation 
132 
\begin{equation} 
133 
D^* \lambda^2 D \approx \frac{\lambda^2}{dx^2} 
134 
\end{equation} 
135 
Therefore we use $G^* \frac{\lambda^2}{dx^2} K^2 G$ as a preconditioner. 
136 
To evaluate the preconditioner we need to solve the equation 
137 
\begin{equation}\label{UPDATE P} 
138 
G^* \frac{\lambda^2}{dx^2} K^2 G p = G^* \hat{r} 
139 
\end{equation} 
140 
It remains to answer the question how to choose $\lambda$. 
141 
We need to balance the first and second term in $J(u,p)$ in \eqn{DARCY COST}. 
142 
We inspect $J$ for $(\hat{u}, \hat{p})$ which is a perturbed exact solution $(u,p)$. 
143 
Assuming $\hat{u}=u+u_{0}e^{ik^tx}$ and $\hat{p}=p+p_{0}e^{ik^tx}$ and 
144 
constant $K$ we get 
145 
\begin{equation} 
146 
J(\hat{u},\hat{p}) = C \left[ ( \K^{1}\_{2} u_{0}^2 + \K\_{2} \k\_{2}^2 p_{0}^2 ) 
147 
+ \lambda^2 \k\_{2}^2 u_{0}^2 \right] 
148 
\end{equation} 
149 
with some constant $C>0$. 
150 
The first two terms and the third term correspond to the first term and 
151 
second term in the definition of $J(u,p)$ in \eqn{DARCY COST}. 
152 
For small $\k\_{2}$ (i.e. for a smooth perturbation) $J(\hat{u},\hat{p})$ is 
153 
dominated by $\K^{1}\_{2} u_{0}^2$. 
154 
To scale the second term which corresponds to the incompressibility condition 
155 
for the velocity we need to meet the condition $\K^{1}\_{2} = \lambda^2 \k\_{2}^2$. 
156 
Taking the boundary conditions into consideration the smallest possible value 
157 
for $\k\_{2} = \frac{\pi}{l}$ where $l$ is the longest edge of the domain. 
158 
This leads to 
159 
\begin{equation}\label{DARCY LAMBDA} 
160 
\lambda = \K^{1}\_{2}^{\frac{1}{2}} \frac{l}{\pi} 
161 
\end{equation} 
162 
Notice that with this setting the preconditioner 
163 
$G^* \frac{\lambda^2}{dx^2} K^2 G$ becomes equivalent to $G^* K G$ if $K$ is a 
164 
diagonal matrix and the mesh has a constant mesh size. 
165 

166 
The residual norm used in PCG is given as 
167 
\begin{equation}\label{DARCY R NORM} 
168 
\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 
169 
\int \hat{r} K^{1} \hat{r} \; dx 
170 
\end{equation} 
171 
The iteration is terminated if 
172 
\begin{equation}\label{DARCY STOP} 
173 
\r\_{PCG} \le \mbox{ATOL} 
174 
\end{equation} 
175 
where we set 
176 
\begin{equation}\label{DARCY ATOL DEF} 
177 
\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} 
178 
\end{equation} 
179 
where rtol is a given relative tolerance and atol is a given absolute 
180 
tolerance (typically $=0$). 
181 
Notice that if $Gp$ and $v$ both are zero, the pair $(0,p)$ is a solution. 
182 
The problem is that ATOL is depending on the solution $p$ and $v$ calculated 
183 
form~\ref{DARCY V FORM}. 
184 
In practice one can use the initial guess for $p$ to get a first value for ATOL. 
185 
If the stopping criterion is met in the PCG iteration, a new $v$ is calculated 
186 
from the current pressure approximation and ATOL is recalculated. 
187 
If \ref{DARCY STOP} is still fulfilled the calculation is terminated and 
188 
$(v,p)$ is returned. Otherwise PCG is restarted with a new ATOL. 
189 

190 
\subsection{Functions} 
191 
\begin{classdesc}{DarcyFlow}{domain} 
192 
opens the Darcy flux problem\index{Darcy flux} on the \Domain domain. 
193 
\end{classdesc} 
194 

195 
\begin{methoddesc}[DarcyFlow]{setValue}{\optional{f=None, \optional{g=None, \optional{location_of_fixed_pressure=None, \optional{location_of_fixed_flux=None, 
196 
\\\optional{permeability=None}}}}}} 
197 
assigns values to the model parameters. Values can be assigned using various 
198 
calls  in particular in a time dependent problem only values that change 
199 
over time need to be reset. The permeability can be defined as a scalar 
200 
(isotropic), a vector (orthotropic) or a matrix (anisotropic). 
201 
\var{f} and \var{g} are the corresponding parameters in~\ref{DARCY PROBLEM}. 
202 
The locations and components where the flux is prescribed are set by positive 
203 
values in \var{location_of_fixed_flux}. 
204 
The locations where the pressure is prescribed are set by by positive values 
205 
of \var{location_of_fixed_pressure}. 
206 
The values of the pressure and flux are defined by the initial guess. 
207 
Notice that at any point on the boundary of the domain the pressure or the 
208 
normal component of the flux must be defined. There must be at least one point 
209 
where the pressure is prescribed. 
210 
The method will try to cast the given values to appropriate \Data class objects. 
211 
\end{methoddesc} 
212 

213 
\begin{methoddesc}[DarcyFlow]{setTolerance}{\optional{rtol=1e4}} 
214 
sets the relative tolerance \var{rtol} in \ref{DARCY ATOL DEF}. 
215 
\end{methoddesc} 
216 

217 
\begin{methoddesc}[DarcyFlow]{setAbsoluteTolerance}{\optional{atol=0.}} 
218 
sets the absolute tolerance \var{atol} in \ref{DARCY ATOL DEF}. 
219 
\end{methoddesc} 
220 

221 
\begin{methoddesc}[DarcyFlow]{getSolverOptionsFlux}{} 
222 
returns the solver options used to solve the flux problems~(\ref{DARCY V FORM}) and~(\ref{UPDATE W}). 
223 
Use this \SolverOptions object to control the solution algorithms. 
224 
\end{methoddesc} 
225 

226 
\begin{methoddesc}[DarcyFlow]{getSolverOptionsPressure}{} 
227 
returns the solver options used to solve the pressure problem~(\ref{UPDATE P}) 
228 
as a preconditioner. 
229 
Use this \SolverOptions object to control the solution algorithms. 
230 
\end{methoddesc} 
231 

232 
\begin{methoddesc}[DarcyFlow]{solve}{u0,p0, \optional{max_iter=100, \optional{verbose=False}}} 
233 
solves the problem and returns approximations for the flux $v$ and the pressure $p$. 
234 
\var{u0} and \var{p0} define initial guesses for flux and pressure. 
235 
Values marked by positive values \var{location_of_fixed_flux} and 
236 
\var{location_of_fixed_pressure}, respectively, are kept unchanged. 
237 
\var{max_iter} sets the maximum number of iteration steps allowed for solving 
238 
the coupled problem. 
239 
\end{methoddesc} 
240 

241 
%\subsection{Example: Gravity Flow} 
242 
%later 
243 
