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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26