1 
\section{Darcy Flux} 
2 
\label{DARCY FLUX} 
3 
We want to calculate the velocity $u$ and pressure $p$ on a domain $\Omega$ solving 
4 
the Darcy flux problem \index{Darcy flux}\index{Darcy flow} 
5 
\begin{equation}\label{DARCY PROBLEM} 
6 
\begin{array}{rcl} 
7 
u_{i} + \kappa_{ij} p_{,j} & = & g_{i} \\ 
8 
u_{k,k} & = & f 
9 
\end{array} 
10 
\end{equation} 
11 
with the boundary conditions 
12 
\begin{equation}\label{DARCY BOUNDARY} 
13 
\begin{array}{rcl} 
14 
u_{i} \; n_{i} = u^{N}_{i} \; n_{i} & \mbox{ on } & \Gamma_{N} \\ 
15 
p = p^{D} & \mbox{ on } & \Gamma_{D} \\ 
16 
\end{array} 
17 
\end{equation} 
18 
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 
19 
\begin{equation} 
20 
\alpha_{0} \; x_{i} x_{i} \le \kappa_{ij} x_{i} x_{j} \le \alpha_{1} \; x_{i} x_{i} 
21 
\end{equation} 
22 
for all $x_{i}$. 
23 

24 
\subsection{Solution Method \label{DARCY SOLVE}} 
25 
It is useful to write equation~\ref{DARCY PROBLEM} in operator form. For any pressure $p$ 
26 
we set 
27 
\begin{equation} 
28 
(Gp)_{i} = p_{,j} 
29 
\end{equation} 
30 
and a velocity $v$ we set 
31 
\begin{equation} 
32 
Dv = v_{k,k} 
33 
\end{equation} 
34 
We $K=(\kappa_{ij})$ we can write the Darcy problem~\ref{DARCY PROBLEM} as 
35 
\begin{equation} 
36 
\begin{array}{rcl} 
37 
u + K \, Gp & = & g \\ 
38 
Du & = & f 
39 
\end{array} 
40 
\end{equation} 
41 
We solve this equation by minimising the functional 
42 
\begin{equation}\label{DARCY COST} 
43 
J(u,p):=\K^{\frac{1}{2}}(u + K \, G p  g)\^2_{0} + \\lambda (Duf) \_{0}^2 
44 
\end{equation} 
45 
over all suitable $u$ and $p$. In this equation we set $\p\^2_{0}=(p,p)_{0}$ with 
46 
\begin{equation} 
47 
(p,q)_{0} = \int_{\Omega } p\cdot q \, dx 
48 
\end{equation} 
49 
The factor $\lambda>0$ is a weighting factor. 
50 
A simple calculation shows that 
51 
one has to solve 
52 
\begin{equation} 
53 
( K^{1} (v + K \, Gq) , u +K \, G p  g)_{0} + (\lambda Dv,\lambda (Duf) )_{0} =0 
54 
\end{equation} 
55 
for all velocities $v$ and pressure $q$ which fulfill the homogeneous boundary conditions~\ref{DARCY BOUNDARY}. 
56 
This socalled variational equation can be translated back into operator notation 
57 
\begin{equation} 
58 
\begin{array}{rcl} 
59 
(K^{1}+ D^*\lambda^2 D)u + Gp & = & D^*\lambda f + K^{1} g \\ 
60 
G^*u + G^*K \, G p & = & G^*g \\ 
61 
\end{array} 
62 
\end{equation} 
63 
where $D^*$ and $Q^*$ denote the adjoint operators with respect to $(.,.)_{0}$. 
64 
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 
65 
to solve the problem in coupled form, however this approach leads in some cases to a very illconditioned stiffness matrix in particular in the case of a very small or large permeability ($\alpha_{1} \ll 1$ or $\alpha_{0} \gg 1$) 
66 

67 
The approach we are taking is to eliminate the velocity $v$ from the problem. Assuming that $p$ is known we have 
68 
\begin{equation}\label{DARCY V FORM} 
69 
v= (K^{1}+ D^*\lambda^2 D)^{1} ( D^*\lambda f + K^{1} g  Gp) 
70 
\end{equation} 
71 
(notice that $K^{1}+\lambda D^*D$ is coercive) which is inserted into the second equation 
72 
\begin{equation} 
73 
G^* (K^{1}+\lambda D^*D)^{1} (\lambda D^*f + K^{1} g  Gp) + G^* KG p = G^*g 
74 
\end{equation} 
75 
which is 
76 
\begin{equation} 
77 
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) ) 
78 
\end{equation} 
79 
We use the PCG \index{linear solver!PCG}\index{PCG} method to solve this. 
80 
The residual $r$ is given as 
81 
\begin{equation} 
82 
\begin{array}{rcl} 
83 
r & =& G^* \left( g  K\, G p  v \right) 
84 
\end{array} 
85 
\end{equation} 
86 
for the current pressure approximation $p$ and current velocity $v$ defined by 
87 
equation~\ref{DARCY V FORM}. 
88 
So in a particular implementation we use $\hat{r}=gK\, Gpv$ to represent the residual. 
89 
The evaluation of the iteration operator for a given $p$ is then 
90 
returning $Qp+v$ where $v$ is the solution of 
91 
\begin{equation}\label{UPDATE W} 
92 
(K^{1}+ D^*\lambda^2 D)v = Gp 
93 
\end{equation} 
94 
To derive a preconditioner we use the identity 
95 
\begin{equation} 
96 
\begin{array}{rcl} 
97 

98 
G^* ( K  (K^{1}+ D^*\lambda^2 D)^{1} ) G & = & G^* (I  (I + K D^*\lambda^2 D)^{1}) K G \\ 
99 
& \approx & G^* (K D^*\lambda^2 D) K G \\ 
100 
& = & G^* K ( D^* \lambda^2 D) K G \\ 
101 
& \approx & G^* \frac{\lambda^2}{dx^2} K^2 G 
102 
\end{array} 
103 
\end{equation} 
104 
where $dx$ is the local mesh size and we use the approximation 
105 
\begin{equation} 
106 
D^* \lambda^2 D \approx \frac{\lambda^2}{dx^2} 
107 
\end{equation} 
108 
Therefore we use $G^* \frac{\lambda^2}{dx^2} K^2 G$ as a preconditioner.To evalute the preconditioner 
109 
we need to solve the equation 
110 
\begin{equation}\label{UPDATE P} 
111 
G^* \frac{\lambda^2}{dx^2} K^2 G p = G^* \hat{r} 
112 
\end{equation} 
113 
It remains to answer the question how to choose $\lambda$. We need to balance the first and second 
114 
term in $J(u,p)$ in equation~\ref{DARCY COST}. We inspect $J$ for 
115 
$(\hat{u}, \hat{p})$ which is a perturbed exact solution $(u,p)$. 
116 
Assuming $\hat{u}=u+u_{0}e^{ik^tx}$ 
117 
and $\hat{p}=p+p_{0}e^{ik^tx}$ and constant $K$ we get 
118 
\begin{equation} 
119 
J(\hat{u},\hat{p}) = C \left[ ( \K^{1}\_{2} u_{0}^2 + \K\_{2} \k\_{2}^2 p_{0}^2 ) 
120 
+ \lambda^2 \k\_{2}^2 u_{0}^2 \right] 
121 
\end{equation} 
122 
with some constant $C>0$. The first two terms and the third term correspond to the first term and 
123 
second term in the definition of $J(u,p)$ in equation~\ref{DARCY COST}. For small $\k\_{2}$ 
124 
(i.e. for a smooth perturbation) $J(\hat{u},\hat{p})$ is dominated by $\K^{1}\_{2} u_{0}^2$. 
125 
To scale the second term which is corresponds to the incompressibility condition for the velocity 
126 
we need to meet the condition $\K^{1}\_{2} = \lambda^2 \k\_{2}^2$. 
127 
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 
128 
\begin{equation}\label{DARCY LAMBDA} 
129 
\lambda = \K^{1}\_{2}^{\frac{1}{2}} \frac{l}{\pi} 
130 
\end{equation} 
131 
Notice that with this setting the preconditioner $G^* \frac{\lambda^2}{dx^2} K^2 G$ becomes 
132 
equivalent to $G^* K G$ if $K$ is a diagonal matrix and the mesh has a constant mesh size. 
133 

134 
The residual norm used in the PCG is given as 
135 
\begin{equation}\label{DARCY R NORM} 
136 
\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 
137 
\int \hat{r} K^{1} \hat{r} \; dx 
138 
\end{equation} 
139 
The iteration is terminated if 
140 
\begin{equation}\label{DARCY STOP} 
141 
\r\_{PCG} \le \mbox{ATOL} 
142 
\end{equation} 
143 
where we set 
144 
\begin{equation}\label{DARCY ATOL DEF} 
145 
\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} 
146 
\end{equation} 
147 
where rtol is a given relative tolerance and $\mbox{atol}$ is a given absolute tolerance (typically $=0$). 
148 
Notice that if $Gp$ and $v$ both are zero, the pair $(0,p)$ is a solution. 
149 
The problem is that ATOL is depending on the solution $p$ and $v$ calculated form~\ref{DARCY V FORM}. 
150 
In practice one use the initial guess for $p$ 
151 
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. 
152 

153 
\subsection{Functions} 
154 
\begin{classdesc}{DarcyFlow}{domain} 
155 
opens the Darcy flux problem\index{Darcy flux} on the \Domain domain. 
156 
\end{classdesc} 
157 

158 
\begin{methoddesc}[DarcyFlow]{setValue}{\optional{f=None, \optional{g=None, \optional{location_of_fixed_pressure=None, \optional{location_of_fixed_flux=None, 
159 
\\\optional{permeability=None}}}}}} 
160 
assigns values to the model parameters. Values can be assigned using various calls  in particular 
161 
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). 
162 
\var{f} and \var{g} are the corresponding parameters in~\ref{DARCY PROBLEM}. 
163 
The locations and components where the flux is prescribed are set by positive values in 
164 
\var{location_of_fixed_flux}. 
165 
The locations where the pressure is prescribed are set by 
166 
by positive values of \var{location_of_fixed_pressure}. 
167 
The values of the pressure and flux are defined by the initial guess. 
168 
Notice that at any point on the boundary of the domain the pressure or the normal component of 
169 
the flux must be defined. There must be at least one point where the pressure is prescribed. 
170 
The method will try to cast the given values to appropriate 
171 
\Data class objects. 
172 
\end{methoddesc} 
173 

174 
\begin{methoddesc}[DarcyFlow]{setTolerance}{\optional{rtol=1e4}} 
175 
sets the relative tolerance \mbox{rtol} in \ref{DARCY ATOL DEF}. 
176 
\end{methoddesc} 
177 

178 
\begin{methoddesc}[DarcyFlow]{setAbsoluteTolerance}{\optional{atol=0.}} 
179 
sets the absolute tolerance \mbox{atol} in \ref{DARCY ATOL DEF}. 
180 
\end{methoddesc} 
181 

182 
\begin{methoddesc}[DarcyFlow]{getSolverOptionsFlux}{} 
183 
Returns the solver options used to solve the flux problems~(\ref{DARCY V FORM}) and~(\ref{UPDATE W}). Use the returned \SolverOptions object to control the solution algorithms. 
184 
\end{methoddesc} 
185 

186 
\begin{methoddesc}[DarcyFlow]{getSolverOptionsPressure}{} 
187 
Returns the solver options used to solve the pressure problems~(\ref{UPDATE P}) as a preconditioner. 
188 
Use the returned \SolverOptions object to control the solution algorithms. 
189 
\end{methoddesc} 
190 

191 
\begin{methoddesc}[DarcyFlow]{solve}{u0,p0, \optional{max_iter=100, \optional{verbose=False}}} 
192 
solves the problem. and returns approximations for the flux $v$ and the pressure $p$. 
193 
\var{u0} and \var{p0} define initial guess for flux and pressure. Values marked 
194 
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. 
195 
\end{methoddesc} 
196 

197 

198 
\subsection{Example: Gravity Flow} 
199 
later 