# Annotation of /trunk/doc/user/darcyflux.tex

Revision 3501 - (hide 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 caltinay 3326 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 3 % 4 % Copyright (c) 2003-2010 by University of Queensland 5 % Earth Systems Science Computational Center (ESSCC) 6 7 % 8 % Primary Business: Queensland, Australia 9 % Licensed under the Open Software License version 3.0 10 11 % 12 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 13 14 gross 2960 \section{Darcy Flux} 15 \label{DARCY FLUX} 16 caltinay 3326 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 gross 2960 \begin{equation}\label{DARCY PROBLEM} 19 \begin{array}{rcl} 20 jfenwick 3295 u_{i} + \kappa_{ij} p_{,j} & = & g_{i} \\ 21 u_{k,k} & = & f 22 gross 2960 \end{array} 23 \end{equation} 24 with the boundary conditions 25 \begin{equation}\label{DARCY BOUNDARY} 26 \begin{array}{rcl} 27 jfenwick 3295 u_{i} \; n_{i} = u^{N}_{i} \; n_{i} & \mbox{ on } & \Gamma_{N} \\ 28 p = p^{D} & \mbox{ on } & \Gamma_{D} \\ 29 gross 2960 \end{array} 30 \end{equation} 31 caltinay 3326 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 gross 2960 \begin{equation} 40 jfenwick 3295 \alpha_{0} \; x_{i} x_{i} \le \kappa_{ij} x_{i} x_{j} \le \alpha_{1} \; x_{i} x_{i} 41 gross 2960 \end{equation} 42 caltinay 3326 for all $x_{i}$. 43 gross 2960 44 gross 3501 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 gross 2960 \subsection{Solution Method \label{DARCY SOLVE}} 57 gross 3501 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 caltinay 3326 It is useful to write \eqn{DARCY PROBLEM} in operator form. 103 For any pressure $p$ we set 104 gross 2960 \begin{equation} 105 jfenwick 3295 (Gp)_{i} = p_{,j} 106 gross 2960 \end{equation} 107 caltinay 3326 and for velocity $v$ we set 108 gross 2960 \begin{equation} 109 jfenwick 3295 Dv = v_{k,k} 110 gross 2960 \end{equation} 111 caltinay 3326 With $K=(\kappa_{ij})$ we can write the Darcy problem~\ref{DARCY PROBLEM} as 112 gross 2960 \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 jfenwick 3295 J(u,p):=\|K^{-\frac{1}{2}}(u + K \, G p - g)\|^2_{0} + \|\lambda (Du-f) \|_{0}^2 121 gross 2960 \end{equation} 122 caltinay 3326 over all suitable $u$ and $p$. 123 In this equation we set $\|p\|^2_{0}=(p,p)_{0}$ with 124 gross 2960 \begin{equation} 125 jfenwick 3295 (p,q)_{0} = \int_{\Omega } p\cdot q \, dx 126 gross 2960 \end{equation} 127 The factor $\lambda>0$ is a weighting factor. 128 caltinay 3326 A simple calculation shows that one has to solve 129 gross 2960 \begin{equation} 130 jfenwick 3295 ( K^{-1} (v + K \, Gq) , u +K \, G p - g)_{0} + (\lambda Dv,\lambda (Du-f) )_{0} =0 131 gross 2960 \end{equation} 132 caltinay 3326 for all velocities $v$ and pressure $q$ which fulfill the homogeneous boundary 133 conditions~\ref{DARCY BOUNDARY}. 134 gross 2960 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 caltinay 3326 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 gross 2960 149 caltinay 3326 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 gross 2960 \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 caltinay 3326 which is inserted into the second equation 156 gross 2960 \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 caltinay 3326 which is 160 gross 2960 \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 caltinay 3326 We use the PCG\index{linear solver!PCG}\index{PCG} method to solve this. 164 The residual $r$ is given as 165 gross 2960 \begin{equation} 166 \begin{array}{rcl} 167 r & =& G^* \left( g - K\, G p - v \right) 168 \end{array} 169 \end{equation} 170 caltinay 3326 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 gross 2960 \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 caltinay 3326 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 gross 2960 \begin{equation}\label{UPDATE P} 195 G^* \frac{\lambda^2}{dx^2} K^2 G p = G^* \hat{r} 196 \end{equation} 197 caltinay 3326 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 gross 2960 \begin{equation} 203 jfenwick 3295 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 gross 2960 \end{equation} 206 caltinay 3326 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 gross 2960 \begin{equation}\label{DARCY LAMBDA} 217 jfenwick 3295 \lambda = \|K^{-1}\|_{2}^{\frac{1}{2}} \frac{l}{\pi} 218 gross 2960 \end{equation} 219 caltinay 3326 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 gross 2960 223 caltinay 3326 The residual norm used in PCG is given as 224 gross 2960 \begin{equation}\label{DARCY R NORM} 225 jfenwick 3295 \|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 gross 2960 \int \hat{r} K^{-1} \hat{r} \; dx 227 \end{equation} 228 caltinay 3326 The iteration is terminated if 229 gross 2960 \begin{equation}\label{DARCY STOP} 230 jfenwick 3295 \|r\|_{PCG} \le \mbox{ATOL} 231 gross 2960 \end{equation} 232 caltinay 3326 where we set 233 gross 2960 \begin{equation}\label{DARCY ATOL DEF} 234 jfenwick 3295 \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 gross 2960 \end{equation} 236 caltinay 3326 where rtol is a given relative tolerance and atol is a given absolute 237 tolerance (typically $=0$). 238 gross 2960 Notice that if $Gp$ and $v$ both are zero, the pair $(0,p)$ is a solution. 239 caltinay 3326 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 gross 2960 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 caltinay 3326 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 gross 2960 \var{f} and \var{g} are the corresponding parameters in~\ref{DARCY PROBLEM}. 259 caltinay 3326 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 gross 2960 The values of the pressure and flux are defined by the initial guess. 264 caltinay 3326 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 gross 2960 \end{methoddesc} 269 270 \begin{methoddesc}[DarcyFlow]{setTolerance}{\optional{rtol=1e-4}} 271 caltinay 3326 sets the relative tolerance \var{rtol} in \ref{DARCY ATOL DEF}. 272 gross 2960 \end{methoddesc} 273 274 \begin{methoddesc}[DarcyFlow]{setAbsoluteTolerance}{\optional{atol=0.}} 275 caltinay 3326 sets the absolute tolerance \var{atol} in \ref{DARCY ATOL DEF}. 276 gross 2960 \end{methoddesc} 277 278 \begin{methoddesc}[DarcyFlow]{getSolverOptionsFlux}{} 279 caltinay 3326 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 gross 2960 \end{methoddesc} 282 283 \begin{methoddesc}[DarcyFlow]{getSolverOptionsPressure}{} 284 caltinay 3331 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 gross 2960 \end{methoddesc} 288 289 \begin{methoddesc}[DarcyFlow]{solve}{u0,p0, \optional{max_iter=100, \optional{verbose=False}}} 290 caltinay 3326 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 gross 2960 \end{methoddesc} 297 298 caltinay 3326 %\subsection{Example: Gravity Flow} 299 %later 300 gross 2960