/[escript]/trunk/doc/user/darcyflux.tex
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3501 - (show annotations)
Wed Apr 13 04:07:53 2011 UTC (8 years, 5 months ago) by gross
File MIME type: application/x-tex
File size: 12580 byte(s)


1
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % Copyright (c) 2003-2010 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/osl-3.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}$ 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 \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 We set $u=\hat{u}+u^{N}$ and
45 $p=\hat{p}+p^D$ where $\hat{u}$ and $\hat{p}$ fullfill the homogeneous boundary conditions~\label{DARCY BOUNDARY} and
46 the PDE
47 \begin{equation}\label{DARCY PROBLEM 2}
48 \begin{array}{rcl}
49 \hat{u}_{i} + \kappa_{ij} \hat{p}_{,j} & = & \hat{g}_{i} = g_{i} -u^{N} - \kappa_{ij} p^D_{,j} \\
50 \hat{u}_{k,k} & = & f - u^N_{k,k} = \hat{f}
51 \end{array}
52 \end{equation}
53 In the following the hat is dropped.
54
55
56 \subsection{Solution Method \label{DARCY SOLVE}}
57 We set
58 \begin{align}
59 Q=\{ q: \Omega \to R | q = 0 \mbox{ on } \Gamma_{D} \} \label{DARCY SPACE P} \\
60 V=\{ v: \Omega \to R^3 |v_{i} \; n_{i} =0 \mbox{ on } \Gamma_{N} \} \label{DARCY SPACE P}
61 \end{align}
62 We define the gardient operator $G: Q \to V$ as
63 \begin{equation}
64 (Gp)_{i} = p_{,j}
65 \end{equation}
66 for all $p \in Q$ and the
67 divergence operator $D: V \to Q$
68 and for velocity $v$ we set
69 \begin{equation}
70 Dv = v_{k,k}
71 \end{equation}
72 As $\Gamma_{D}, \Gamma_{N}$ are a subdivision of the boundary of $\Omega$ we have
73 \begin{equation}
74 D^*= - G
75 \end{equation}
76 With $K=(\kappa_{ij})$ we can write the Darcy problem~\ref{DARCY PROBLEM} as
77 \begin{equation}
78 \begin{array}{rcl}
79 K^{-1} u - D^* p & = & K^{-1} g \\
80 Du & = & f
81 \end{array}
82 \end{equation}
83 This problem needs to be stabalized in the following way:
84 \begin{equation}
85 \begin{array}{rcl}
86 K^{-1} u - D^* p - \frac{1}{2} \left( K^{-1} u +G p \right) & = & K^{-1} g - \frac{1}{2} K^{-1} g \\
87 Du + \frac{1}{2} G^* K \left( K^{-1} u +G p \right) & = & f + \frac{1}{2} G^* K K^{-1} g
88 \end{array}
89 \end{equation}
90 which can be simplified to
91 \begin{equation}
92 \begin{array}{rcl}
93 \frac{1}{2} K^{-1} u - D^* p - \frac{1}{2} G p & = & \frac{1}{2} K^{-1} g \\
94 Du + \frac{1}{2} G^* u + \frac{1}{2} G^* K G p & = & f + \frac{1}{2} G^* g
95 \end{array}
96 \end{equation}
97
98
99
100 %===========================================================================
101 \subsection{Solution Method \label{DARCY SOLVE}}
102 It is useful to write \eqn{DARCY PROBLEM} in operator form.
103 For any pressure $p$ we set
104 \begin{equation}
105 (Gp)_{i} = p_{,j}
106 \end{equation}
107 and for velocity $v$ we set
108 \begin{equation}
109 Dv = v_{k,k}
110 \end{equation}
111 With $K=(\kappa_{ij})$ we can write the Darcy problem~\ref{DARCY PROBLEM} as
112 \begin{equation}
113 \begin{array}{rcl}
114 u + K \, Gp & = & g \\
115 Du & = & f
116 \end{array}
117 \end{equation}
118 We solve this equation by minimising the functional
119 \begin{equation}\label{DARCY COST}
120 J(u,p):=\|K^{-\frac{1}{2}}(u + K \, G p - g)\|^2_{0} + \|\lambda (Du-f) \|_{0}^2
121 \end{equation}
122 over all suitable $u$ and $p$.
123 In this equation we set $\|p\|^2_{0}=(p,p)_{0}$ with
124 \begin{equation}
125 (p,q)_{0} = \int_{\Omega } p\cdot q \, dx
126 \end{equation}
127 The factor $\lambda>0$ is a weighting factor.
128 A simple calculation shows that one has to solve
129 \begin{equation}
130 ( K^{-1} (v + K \, Gq) , u +K \, G p - g)_{0} + (\lambda Dv,\lambda (Du-f) )_{0} =0
131 \end{equation}
132 for all velocities $v$ and pressure $q$ which fulfill the homogeneous boundary
133 conditions~\ref{DARCY BOUNDARY}.
134 This so-called variational equation can be translated back into operator notation
135 \begin{equation}
136 \begin{array}{rcl}
137 (K^{-1}+ D^*\lambda^2 D)u + Gp & = & D^*\lambda f + K^{-1} g \\
138 G^*u + G^*K \, G p & = & G^*g \\
139 \end{array}
140 \end{equation}
141 where $D^*$ and $Q^*$ denote the adjoint operators with respect to $(.,.)_{0}$.
142 In~\cite{LEASTSQUARESFEM1994} it has been shown that this problem is
143 continuous and coercive and therefore has a unique solution.
144 Standard FEM methods can be used for discretization. It is also possible to
145 solve the problem in coupled form, however this approach leads in some cases
146 to a very ill-conditioned stiffness matrix, in particular in the case of a
147 very small or large permeability ($\alpha_{1} \ll 1$ or $\alpha_{0} \gg 1$).
148
149 The approach we are taking is to eliminate the velocity $v$ from the problem.
150 Assuming that $p$ is known we have\footnote{notice
151 that $K^{-1}+\lambda D^*D$ is coercive}
152 \begin{equation}\label{DARCY V FORM}
153 v= (K^{-1}+ D^*\lambda^2 D)^{-1} ( D^*\lambda f + K^{-1} g - Gp)
154 \end{equation}
155 which is inserted into the second equation
156 \begin{equation}
157 G^* (K^{-1}+\lambda D^*D)^{-1} (\lambda D^*f + K^{-1} g - Gp) + G^* KG p = G^*g
158 \end{equation}
159 which is
160 \begin{equation}
161 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) )
162 \end{equation}
163 We use the PCG\index{linear solver!PCG}\index{PCG} method to solve this.
164 The residual $r$ is given as
165 \begin{equation}
166 \begin{array}{rcl}
167 r & =& G^* \left( g - K\, G p - v \right)
168 \end{array}
169 \end{equation}
170 for the current pressure approximation $p$ and current velocity $v$ defined by
171 \eqn{DARCY V FORM}.
172 So in a particular implementation we use $\hat{r}=g-K\, Gp-v$ to represent the residual.
173 The evaluation of the iteration operator for a given $p$ is then returning
174 $Qp+v$ where $v$ is the solution of
175 \begin{equation}\label{UPDATE W}
176 (K^{-1}+ D^*\lambda^2 D)v = Gp
177 \end{equation}
178 To derive a preconditioner we use the identity
179 \begin{equation}
180 \begin{array}{rcl}
181
182 G^* ( K - (K^{-1}+ D^*\lambda^2 D)^{-1} ) G & = & G^* (I - (I + K D^*\lambda^2 D)^{-1}) K G \\
183 & \approx & G^* (K D^*\lambda^2 D) K G \\
184 & = & G^* K ( D^* \lambda^2 D) K G \\
185 & \approx & G^* \frac{\lambda^2}{dx^2} K^2 G
186 \end{array}
187 \end{equation}
188 where $dx$ is the local mesh size and we use the approximation
189 \begin{equation}
190 D^* \lambda^2 D \approx \frac{\lambda^2}{dx^2}
191 \end{equation}
192 Therefore we use $G^* \frac{\lambda^2}{dx^2} K^2 G$ as a preconditioner.
193 To evaluate the preconditioner we need to solve the equation
194 \begin{equation}\label{UPDATE P}
195 G^* \frac{\lambda^2}{dx^2} K^2 G p = G^* \hat{r}
196 \end{equation}
197 It remains to answer the question how to choose $\lambda$.
198 We need to balance the first and second term in $J(u,p)$ in \eqn{DARCY COST}.
199 We inspect $J$ for $(\hat{u}, \hat{p})$ which is a perturbed exact solution $(u,p)$.
200 Assuming $\hat{u}=u+u_{0}e^{ik^tx}$ and $\hat{p}=p+p_{0}e^{ik^tx}$ and
201 constant $K$ we get
202 \begin{equation}
203 J(\hat{u},\hat{p}) = C \left[ ( \|K^{-1}\|_{2} |u_{0}|^2 + \|K\|_{2} \|k\|_{2}^2 |p_{0}|^2| )
204 + \lambda^2 \|k\|_{2}^2 |u_{0}|^2 \right]
205 \end{equation}
206 with some constant $C>0$.
207 The first two terms and the third term correspond to the first term and
208 second term in the definition of $J(u,p)$ in \eqn{DARCY COST}.
209 For small $\|k\|_{2}$ (i.e. for a smooth perturbation) $J(\hat{u},\hat{p})$ is
210 dominated by $\|K^{-1}\|_{2} |u_{0}|^2$.
211 To scale the second term which corresponds to the incompressibility condition
212 for the velocity we need to meet the condition $\|K^{-1}\|_{2} = \lambda^2 \|k\|_{2}^2$.
213 Taking the boundary conditions into consideration the smallest possible value
214 for $\|k\|_{2} = \frac{\pi}{l}$ where $l$ is the longest edge of the domain.
215 This leads to
216 \begin{equation}\label{DARCY LAMBDA}
217 \lambda = \|K^{-1}\|_{2}^{\frac{1}{2}} \frac{l}{\pi}
218 \end{equation}
219 Notice that with this setting the preconditioner
220 $G^* \frac{\lambda^2}{dx^2} K^2 G$ becomes equivalent to $G^* K G$ if $K$ is a
221 diagonal matrix and the mesh has a constant mesh size.
222
223 The residual norm used in PCG is given as
224 \begin{equation}\label{DARCY R NORM}
225 \|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
226 \int \hat{r} K^{-1} \hat{r} \; dx
227 \end{equation}
228 The iteration is terminated if
229 \begin{equation}\label{DARCY STOP}
230 \|r\|_{PCG} \le \mbox{ATOL}
231 \end{equation}
232 where we set
233 \begin{equation}\label{DARCY ATOL DEF}
234 \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}
235 \end{equation}
236 where rtol is a given relative tolerance and atol is a given absolute
237 tolerance (typically $=0$).
238 Notice that if $Gp$ and $v$ both are zero, the pair $(0,p)$ is a solution.
239 The problem is that ATOL is depending on the solution $p$ and $v$ calculated
240 form~\ref{DARCY V FORM}.
241 In practice one can use the initial guess for $p$ to get a first value for ATOL.
242 If the stopping criterion is met in the PCG iteration, a new $v$ is calculated
243 from the current pressure approximation and ATOL is recalculated.
244 If \ref{DARCY STOP} is still fulfilled the calculation is terminated and
245 $(v,p)$ is returned. Otherwise PCG is restarted with a new ATOL.
246
247 \subsection{Functions}
248 \begin{classdesc}{DarcyFlow}{domain}
249 opens the Darcy flux problem\index{Darcy flux} on the \Domain domain.
250 \end{classdesc}
251
252 \begin{methoddesc}[DarcyFlow]{setValue}{\optional{f=None, \optional{g=None, \optional{location_of_fixed_pressure=None, \optional{location_of_fixed_flux=None,
253 \\\optional{permeability=None}}}}}}
254 assigns values to the model parameters. Values can be assigned using various
255 calls -- in particular in a time dependent problem only values that change
256 over time need to be reset. The permeability can be defined as a scalar
257 (isotropic), a vector (orthotropic) or a matrix (anisotropic).
258 \var{f} and \var{g} are the corresponding parameters in~\ref{DARCY PROBLEM}.
259 The locations and components where the flux is prescribed are set by positive
260 values in \var{location_of_fixed_flux}.
261 The locations where the pressure is prescribed are set by by positive values
262 of \var{location_of_fixed_pressure}.
263 The values of the pressure and flux are defined by the initial guess.
264 Notice that at any point on the boundary of the domain the pressure or the
265 normal component of the flux must be defined. There must be at least one point
266 where the pressure is prescribed.
267 The method will try to cast the given values to appropriate \Data class objects.
268 \end{methoddesc}
269
270 \begin{methoddesc}[DarcyFlow]{setTolerance}{\optional{rtol=1e-4}}
271 sets the relative tolerance \var{rtol} in \ref{DARCY ATOL DEF}.
272 \end{methoddesc}
273
274 \begin{methoddesc}[DarcyFlow]{setAbsoluteTolerance}{\optional{atol=0.}}
275 sets the absolute tolerance \var{atol} in \ref{DARCY ATOL DEF}.
276 \end{methoddesc}
277
278 \begin{methoddesc}[DarcyFlow]{getSolverOptionsFlux}{}
279 returns the solver options used to solve the flux problems~(\ref{DARCY V FORM}) and~(\ref{UPDATE W}).
280 Use this \SolverOptions object to control the solution algorithms.
281 \end{methoddesc}
282
283 \begin{methoddesc}[DarcyFlow]{getSolverOptionsPressure}{}
284 returns a \SolverOptions object with the options used to solve the pressure
285 problem~(\ref{UPDATE P}) as a preconditioner.
286 Use this object to control the solution algorithms.
287 \end{methoddesc}
288
289 \begin{methoddesc}[DarcyFlow]{solve}{u0,p0, \optional{max_iter=100, \optional{verbose=False}}}
290 solves the problem and returns approximations for the flux $v$ and the pressure $p$.
291 \var{u0} and \var{p0} define initial guesses for flux and pressure.
292 Values marked by positive values \var{location_of_fixed_flux} and
293 \var{location_of_fixed_pressure}, respectively, are kept unchanged.
294 \var{max_iter} sets the maximum number of iteration steps allowed for solving
295 the coupled problem.
296 \end{methoddesc}
297
298 %\subsection{Example: Gravity Flow}
299 %later
300

  ViewVC Help
Powered by ViewVC 1.1.26