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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3326 - (show annotations)
Fri Oct 29 00:45:32 2010 UTC (8 years, 2 months ago) by caltinay
File MIME type: application/x-tex
File size: 10707 byte(s)
6.2 Darcy Flux. Note that I commented the gravity flow subsection since it's
not written.

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 \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 (Du-f) \|_{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 (Du-f) )_{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 so-called 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 ill-conditioned 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}=g-K\, Gp-v$ 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=1e-4}}
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

  ViewVC Help
Powered by ViewVC 1.1.26