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

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 7 % 8 % Primary Business: Queensland, Australia 9 % Licensed under the Open Software License version 3.0 10 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