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

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

Parent Directory Parent Directory | Revision Log Revision Log


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     % 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 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

  ViewVC Help
Powered by ViewVC 1.1.26