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

Revision 3326 - (show annotations)
Fri Oct 29 00:45:32 2010 UTC (8 years, 10 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 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 \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