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

Revision 3326 - (hide annotations)
Fri Oct 29 00:45:32 2010 UTC (8 years, 11 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 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 \subsection{Solution Method \label{DARCY SOLVE}} 45 caltinay 3326 It is useful to write \eqn{DARCY PROBLEM} in operator form. 46 For any pressure $p$ we set 47 gross 2960 \begin{equation} 48 jfenwick 3295 (Gp)_{i} = p_{,j} 49 gross 2960 \end{equation} 50 caltinay 3326 and for velocity $v$ we set 51 gross 2960 \begin{equation} 52 jfenwick 3295 Dv = v_{k,k} 53 gross 2960 \end{equation} 54 caltinay 3326 With $K=(\kappa_{ij})$ we can write the Darcy problem~\ref{DARCY PROBLEM} as 55 gross 2960 \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 jfenwick 3295 J(u,p):=\|K^{-\frac{1}{2}}(u + K \, G p - g)\|^2_{0} + \|\lambda (Du-f) \|_{0}^2 64 gross 2960 \end{equation} 65 caltinay 3326 over all suitable $u$ and $p$. 66 In this equation we set $\|p\|^2_{0}=(p,p)_{0}$ with 67 gross 2960 \begin{equation} 68 jfenwick 3295 (p,q)_{0} = \int_{\Omega } p\cdot q \, dx 69 gross 2960 \end{equation} 70 The factor $\lambda>0$ is a weighting factor. 71 caltinay 3326 A simple calculation shows that one has to solve 72 gross 2960 \begin{equation} 73 jfenwick 3295 ( K^{-1} (v + K \, Gq) , u +K \, G p - g)_{0} + (\lambda Dv,\lambda (Du-f) )_{0} =0 74 gross 2960 \end{equation} 75 caltinay 3326 for all velocities $v$ and pressure $q$ which fulfill the homogeneous boundary 76 conditions~\ref{DARCY BOUNDARY}. 77 gross 2960 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 caltinay 3326 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 gross 2960 92 caltinay 3326 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 gross 2960 \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 caltinay 3326 which is inserted into the second equation 99 gross 2960 \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 caltinay 3326 which is 103 gross 2960 \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 caltinay 3326 We use the PCG\index{linear solver!PCG}\index{PCG} method to solve this. 107 The residual $r$ is given as 108 gross 2960 \begin{equation} 109 \begin{array}{rcl} 110 r & =& G^* \left( g - K\, G p - v \right) 111 \end{array} 112 \end{equation} 113 caltinay 3326 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 gross 2960 \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 caltinay 3326 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 gross 2960 \begin{equation}\label{UPDATE P} 138 G^* \frac{\lambda^2}{dx^2} K^2 G p = G^* \hat{r} 139 \end{equation} 140 caltinay 3326 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 gross 2960 \begin{equation} 146 jfenwick 3295 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 gross 2960 \end{equation} 149 caltinay 3326 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 gross 2960 \begin{equation}\label{DARCY LAMBDA} 160 jfenwick 3295 \lambda = \|K^{-1}\|_{2}^{\frac{1}{2}} \frac{l}{\pi} 161 gross 2960 \end{equation} 162 caltinay 3326 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 gross 2960 166 caltinay 3326 The residual norm used in PCG is given as 167 gross 2960 \begin{equation}\label{DARCY R NORM} 168 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 169 gross 2960 \int \hat{r} K^{-1} \hat{r} \; dx 170 \end{equation} 171 caltinay 3326 The iteration is terminated if 172 gross 2960 \begin{equation}\label{DARCY STOP} 173 jfenwick 3295 \|r\|_{PCG} \le \mbox{ATOL} 174 gross 2960 \end{equation} 175 caltinay 3326 where we set 176 gross 2960 \begin{equation}\label{DARCY ATOL DEF} 177 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} 178 gross 2960 \end{equation} 179 caltinay 3326 where rtol is a given relative tolerance and atol is a given absolute 180 tolerance (typically $=0$). 181 gross 2960 Notice that if $Gp$ and $v$ both are zero, the pair $(0,p)$ is a solution. 182 caltinay 3326 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 gross 2960 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 caltinay 3326 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 gross 2960 \var{f} and \var{g} are the corresponding parameters in~\ref{DARCY PROBLEM}. 202 caltinay 3326 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 gross 2960 The values of the pressure and flux are defined by the initial guess. 207 caltinay 3326 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 gross 2960 \end{methoddesc} 212 213 \begin{methoddesc}[DarcyFlow]{setTolerance}{\optional{rtol=1e-4}} 214 caltinay 3326 sets the relative tolerance \var{rtol} in \ref{DARCY ATOL DEF}. 215 gross 2960 \end{methoddesc} 216 217 \begin{methoddesc}[DarcyFlow]{setAbsoluteTolerance}{\optional{atol=0.}} 218 caltinay 3326 sets the absolute tolerance \var{atol} in \ref{DARCY ATOL DEF}. 219 gross 2960 \end{methoddesc} 220 221 \begin{methoddesc}[DarcyFlow]{getSolverOptionsFlux}{} 222 caltinay 3326 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 gross 2960 \end{methoddesc} 225 226 \begin{methoddesc}[DarcyFlow]{getSolverOptionsPressure}{} 227 caltinay 3326 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 gross 2960 \end{methoddesc} 231 232 \begin{methoddesc}[DarcyFlow]{solve}{u0,p0, \optional{max_iter=100, \optional{verbose=False}}} 233 caltinay 3326 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 gross 2960 \end{methoddesc} 240 241 caltinay 3326 %\subsection{Example: Gravity Flow} 242 %later 243 gross 2960