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

revision 3301 by jfenwick, Mon Oct 25 01:27:54 2010 UTC revision 3326 by caltinay, Fri Oct 29 00:45:32 2010 UTC
# Line 1  Line 1
1
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
11    %
12    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
13
14  \section{Darcy Flux}  \section{Darcy Flux}
15  \label{DARCY FLUX}  \label{DARCY FLUX}
16  We want to calculate the velocity $u$ and pressure $p$ on a domain $\Omega$ solving  We want to calculate the velocity $u$ and pressure $p$ on a domain $\Omega$
17  the Darcy flux problem \index{Darcy flux}\index{Darcy flow}  solving the Darcy flux problem\index{Darcy flux}\index{Darcy flow}
18  \label{DARCY PROBLEM}  \label{DARCY PROBLEM}
19  \begin{array}{rcl}  \begin{array}{rcl}
20  u_{i} + \kappa_{ij} p_{,j} & = & g_{i} \\  u_{i} + \kappa_{ij} p_{,j} & = & g_{i} \\
# Line 15  u_{i} \; n_{i}  = u^{N}_{i}  \; n_{i} & Line 28  u_{i} \; n_{i}  = u^{N}_{i}  \; n_{i} &
28  p = p^{D} &  \mbox{ on } & \Gamma_{D} \\  p = p^{D} &  \mbox{ on } & \Gamma_{D} \\
29  \end{array}  \end{array}
30
31  where $\Gamma_{N}$ and $\Gamma_{D}$ are a partition of the boundary of $\Omega$ with $\Gamma_{D}$ non empty, $n_{i}$ is the outer normal field of the boundary of $\Omega$, $u^{N}_{i}$ and $p^{D}$ are given functions on $\Omega$, $g_{i}$ and $f$ are given source terms and $\kappa_{ij}$ is the given permeability. We assume that $\kappa_{ij}$ is symmetric (which is not really required) and positive definite, i.e there are positive constants $\alpha_{0}$ and $\alpha_{1}$ which are independent from the location in $\Omega$ such that  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
40  \alpha_{0} \; x_{i} x_{i} \le \kappa_{ij} x_{i} x_{j} \le \alpha_{1} \; x_{i} x_{i}  \alpha_{0} \; x_{i} x_{i} \le \kappa_{ij} x_{i} x_{j} \le \alpha_{1} \; x_{i} x_{i}
41
42  for all $x_{i}$.  for all $x_{i}$.
43
44  \subsection{Solution Method \label{DARCY SOLVE}}  \subsection{Solution Method \label{DARCY SOLVE}}
45  It is useful to write equation~\ref{DARCY PROBLEM} in operator form. For any pressure $p$  It is useful to write \eqn{DARCY PROBLEM} in operator form.
46  we set  For any pressure $p$ we set
47
48  (Gp)_{i} =  p_{,j}  (Gp)_{i} =  p_{,j}
49
50  and a velocity $v$ we set  and for velocity $v$ we set
51
52  Dv = v_{k,k}  Dv = v_{k,k}
53
54  We $K=(\kappa_{ij})$ we can write the Darcy problem~\ref{DARCY PROBLEM} as  With $K=(\kappa_{ij})$ we can write the Darcy problem~\ref{DARCY PROBLEM} as
55
56  \begin{array}{rcl}  \begin{array}{rcl}
57  u + K \, Gp & = & g \\  u + K \, Gp & = & g \\
# Line 42  We solve this equation by minimising the Line 62  We solve this equation by minimising the
62  \label{DARCY COST}  \label{DARCY COST}
63  J(u,p):=\|K^{-\frac{1}{2}}(u + K \, G p - g)\|^2_{0} +  \|\lambda (Du-f) \|_{0}^2  J(u,p):=\|K^{-\frac{1}{2}}(u + K \, G p - g)\|^2_{0} +  \|\lambda (Du-f) \|_{0}^2
64
65  over all suitable $u$ and $p$. In this equation we set $\|p\|^2_{0}=(p,p)_{0}$ with  over all suitable $u$ and $p$.
66    In this equation we set $\|p\|^2_{0}=(p,p)_{0}$ with
67
68  (p,q)_{0} = \int_{\Omega } p\cdot q \, dx  (p,q)_{0} = \int_{\Omega } p\cdot q \, dx
69
70  The factor $\lambda>0$ is a weighting factor.  The factor $\lambda>0$ is a weighting factor.
71  A simple calculation shows that  A simple calculation shows that one has to solve
one has to solve
72
73  ( K^{-1} (v + K \, Gq) , u +K \, G p - g)_{0} +  (\lambda Dv,\lambda (Du-f) )_{0} =0  ( K^{-1} (v + K \, Gq) , u +K \, G p - g)_{0} +  (\lambda Dv,\lambda (Du-f) )_{0} =0
74
75  for all velocities $v$ and pressure $q$ which fulfill the homogeneous boundary conditions~\ref{DARCY BOUNDARY}.  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  This so-called variational equation can be translated back into operator notation
78
79  \begin{array}{rcl}  \begin{array}{rcl}
# Line 60  This so-called variational equation can Line 81  This so-called variational equation can
81  G^*u  + G^*K \, G p & = & G^*g \\  G^*u  + G^*K \, G p & = & G^*g \\
82  \end{array}  \end{array}
83
84  where $D^*$ and $Q^*$ denote the  adjoint operators with respect to $(.,.)_{0}$.  where $D^*$ and $Q^*$ denote the adjoint operators with respect to $(.,.)_{0}$.
85  In~\cite{LEASTSQUARESFEM1994} it has been shown that this problem is continuous and coercive and therefore has a unique solution. Also standard FEM methods can be used for discretization. It is also possible  In~\cite{LEASTSQUARESFEM1994} it has been shown that this problem is
86  to solve the problem in coupled form, however this approach leads in some cases to a very ill-conditioned stiffness matrix in particular in the case of a very small or large permeability ($\alpha_{1} \ll 1$ or $\alpha_{0} \gg 1$)    continuous and coercive and therefore has a unique solution.
87    Standard FEM methods can be used for discretization. It is also possible to
88  The approach we are taking is to eliminate the velocity $v$ from the problem. Assuming that $p$ is known we have  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  \label{DARCY V FORM}  \label{DARCY V FORM}
96  v= (K^{-1}+ D^*\lambda^2 D)^{-1} ( D^*\lambda f + K^{-1} g - Gp)  v= (K^{-1}+ D^*\lambda^2 D)^{-1} ( D^*\lambda f + K^{-1} g - Gp)
97
98  (notice that $K^{-1}+\lambda D^*D$ is coercive) which is inserted into the second equation  which is inserted into the second equation
99
100  G^* (K^{-1}+\lambda D^*D)^{-1} (\lambda D^*f + K^{-1} g - Gp) + G^* KG p = G^*g  G^* (K^{-1}+\lambda D^*D)^{-1} (\lambda D^*f + K^{-1} g - Gp) + G^* KG p = G^*g
101
102  which is  which is
103
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) )  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
106  We use the PCG \index{linear solver!PCG}\index{PCG} method to solve this.  We use the PCG\index{linear solver!PCG}\index{PCG} method to solve this.
107  The residual $r$  is given as  The residual $r$ is given as
108
109  \begin{array}{rcl}  \begin{array}{rcl}
110  r & =&  G^* \left( g - K\, G p - v \right)  r & =&  G^* \left( g - K\, G p - v \right)
111  \end{array}  \end{array}
112
113  for the current pressure approximation $p$ and current velocity $v$ defined by    for the current pressure approximation $p$ and current velocity $v$ defined by
114  equation~\ref{DARCY V FORM}.  \eqn{DARCY V FORM}.
115  So in a particular implementation we use $\hat{r}=g-K\, Gp-v$ to represent the residual.  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  The evaluation of the iteration operator for a given $p$ is then returning
117  returning $Qp+v$ where $v$ is the solution of  $Qp+v$ where $v$ is the solution of
118  \label{UPDATE W}  \label{UPDATE W}
119  (K^{-1}+ D^*\lambda^2 D)v = Gp  (K^{-1}+ D^*\lambda^2 D)v = Gp
120
# Line 105  where $dx$ is the local mesh size and we Line 132  where $dx$ is the local mesh size and we
132
133  D^*  \lambda^2 D \approx \frac{\lambda^2}{dx^2}  D^*  \lambda^2 D \approx \frac{\lambda^2}{dx^2}
134
135  Therefore we use $G^* \frac{\lambda^2}{dx^2} K^2 G$ as a preconditioner.To evalute the preconditioner  Therefore we use $G^* \frac{\lambda^2}{dx^2} K^2 G$ as a preconditioner.
136  we need to solve the equation  To evaluate the preconditioner we need to solve the equation
137  \label{UPDATE P}  \label{UPDATE P}
138   G^* \frac{\lambda^2}{dx^2} K^2 G p =  G^* \hat{r}   G^* \frac{\lambda^2}{dx^2} K^2 G p =  G^* \hat{r}
139
140  It remains to answer the question how to choose $\lambda$. We need to balance the first and second  It remains to answer the question how to choose $\lambda$.
141  term in $J(u,p)$ in equation~\ref{DARCY COST}. We inspect $J$ for  We need to balance the first and second term in $J(u,p)$ in \eqn{DARCY COST}.
142  $(\hat{u}, \hat{p})$ which is a perturbed exact solution $(u,p)$.    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}$  Assuming $\hat{u}=u+u_{0}e^{ik^tx}$ and $\hat{p}=p+p_{0}e^{ik^tx}$ and
144  and $\hat{p}=p+p_{0}e^{ik^tx}$ and constant $K$ we get  constant $K$ we get
145
146  J(\hat{u},\hat{p}) = C \left[ ( \|K^{-1}\|_{2} |u_{0}|^2 + \|K\|_{2} \|k\|_{2}^2 |p_{0}|^2| )  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]  + \lambda^2 \|k\|_{2}^2 |u_{0}|^2  \right]
148
149  with some constant $C>0$. The first two terms and the third term correspond to the first term and  with some constant $C>0$.
150  second term in the definition of $J(u,p)$ in equation~\ref{DARCY COST}. For small $\|k\|_{2}$  The first two terms and the third term correspond to the first term and
151  (i.e. for a smooth perturbation) $J(\hat{u},\hat{p})$ is dominated by $\|K^{-1}\|_{2} |u_{0}|^2$.  second term in the definition of $J(u,p)$ in \eqn{DARCY COST}.
152  To scale the second term which is corresponds to the incompressibility condition for the velocity  For small $\|k\|_{2}$ (i.e. for a smooth perturbation) $J(\hat{u},\hat{p})$ is
153  we need to meet the condition $\|K^{-1}\|_{2} = \lambda^2 \|k\|_{2}^2$.  dominated by $\|K^{-1}\|_{2} |u_{0}|^2$.
154  Taking the boundary conditions into consideration the smallest possible value for $\|k\|_{2} = \frac{\pi}{l}$ where $l$ is the longest edge of the domain. This leads to the  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.
159  \label{DARCY LAMBDA}  \label{DARCY LAMBDA}
160  \lambda = \|K^{-1}\|_{2}^{\frac{1}{2}} \frac{l}{\pi}  \lambda = \|K^{-1}\|_{2}^{\frac{1}{2}} \frac{l}{\pi}
161
162  Notice that with this setting the preconditioner $G^* \frac{\lambda^2}{dx^2} K^2 G$ becomes  Notice that with this setting the preconditioner
163  equivalent to $G^* K G$ if $K$ is a diagonal matrix and the mesh has a constant mesh size.  $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 the PCG is given as  The residual norm used in PCG is given as
167  \label{DARCY R NORM}  \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  \|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  \int \hat{r} K^{-1}  \hat{r} \; dx
170
171  The iteration is terminated if  The iteration is terminated if
172  \label{DARCY STOP}  \label{DARCY STOP}
173  \|r\|_{PCG} \le \mbox{ATOL}  \|r\|_{PCG} \le \mbox{ATOL}
174
175  where we set  where we set
176  \label{DARCY ATOL DEF}  \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}  \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
179  where rtol is a given relative tolerance and $\mbox{atol}$ is a given absolute tolerance (typically $=0$).    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.  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 form~\ref{DARCY V FORM}.  The problem is that ATOL is depending on the solution $p$ and $v$ calculated
183  In practice one use the initial guess for $p$  form~\ref{DARCY V FORM}.
184  to get a first value for ATOL. If the stopping criterion is met in the PCG iteration, a new $v$ is calculated from the current pressure approximation and ATOL is recalculated. If \ref{DARCY STOP} is still fulfilled the calculation is terminated and $(v,p)$ is returned. Otherwise PCG is restarted with a new ATOL.  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}  \subsection{Functions}
191  \begin{classdesc}{DarcyFlow}{domain}  \begin{classdesc}{DarcyFlow}{domain}
# Line 157  opens the Darcy flux problem\index{Darcy Line 194  opens the Darcy flux problem\index{Darcy
194
195  \begin{methoddesc}[DarcyFlow]{setValue}{\optional{f=None, \optional{g=None, \optional{location_of_fixed_pressure=None, \optional{location_of_fixed_flux=None,  \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}}}}}}  \\\optional{permeability=None}}}}}}
197  assigns values to the model parameters. Values can be assigned using various calls - in particular  assigns values to the model parameters. Values can be assigned using various
198  in a time dependent problem only values that change over time needs to be reset. The permeability can be defined as scalar (isotropic), a vector (orthotropic) or a matrix (anisotropic).  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}.  \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 values in  The locations and components where the flux is prescribed are set by positive
203  \var{location_of_fixed_flux}.  values in \var{location_of_fixed_flux}.
204  The locations where the pressure is prescribed are set by  The locations where the pressure is prescribed are set by by positive values
205  by positive values of \var{location_of_fixed_pressure}.  of \var{location_of_fixed_pressure}.
206  The values of the pressure and flux are defined by the initial guess.  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 normal component of  Notice that at any point on the boundary of the domain the pressure or the
208  the flux must be defined. There must be at least one point where the pressure is prescribed.  normal component of the flux must be defined. There must be at least one point
209  The method will try to cast the given values to appropriate  where the pressure is prescribed.
210  \Data class objects.  The method will try to cast the given values to appropriate \Data class objects.
211  \end{methoddesc}  \end{methoddesc}
212
213  \begin{methoddesc}[DarcyFlow]{setTolerance}{\optional{rtol=1e-4}}  \begin{methoddesc}[DarcyFlow]{setTolerance}{\optional{rtol=1e-4}}
214  sets the relative tolerance \mbox{rtol} in \ref{DARCY ATOL DEF}.  sets the relative tolerance \var{rtol} in \ref{DARCY ATOL DEF}.
215  \end{methoddesc}  \end{methoddesc}
216
217  \begin{methoddesc}[DarcyFlow]{setAbsoluteTolerance}{\optional{atol=0.}}  \begin{methoddesc}[DarcyFlow]{setAbsoluteTolerance}{\optional{atol=0.}}
218  sets the absolute tolerance \mbox{atol} in \ref{DARCY ATOL DEF}.  sets the absolute tolerance \var{atol} in \ref{DARCY ATOL DEF}.
219  \end{methoddesc}  \end{methoddesc}
220
221  \begin{methoddesc}[DarcyFlow]{getSolverOptionsFlux}{}  \begin{methoddesc}[DarcyFlow]{getSolverOptionsFlux}{}
222  Returns the solver options used to solve the flux problems~(\ref{DARCY V FORM}) and~(\ref{UPDATE W}). Use this  returns the solver options used to solve the flux problems~(\ref{DARCY V FORM}) and~(\ref{UPDATE W}).
223   \SolverOptions object to control the solution algorithms.  Use this \SolverOptions object to control the solution algorithms.
224  \end{methoddesc}  \end{methoddesc}
225
226  \begin{methoddesc}[DarcyFlow]{getSolverOptionsPressure}{}  \begin{methoddesc}[DarcyFlow]{getSolverOptionsPressure}{}
227  Returns the solver options used to solve the pressure problems~(\ref{UPDATE P}) as a preconditioner.  returns the solver options used to solve the pressure problem~(\ref{UPDATE P})
228  Use this \SolverOptions object to control the solution algorithms.  as a preconditioner.
229    Use this \SolverOptions object to control the solution algorithms.
230  \end{methoddesc}  \end{methoddesc}
231
232  \begin{methoddesc}[DarcyFlow]{solve}{u0,p0, \optional{max_iter=100, \optional{verbose=False}}}  \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$.  solves the problem and returns approximations for the flux $v$ and the pressure $p$.
234  \var{u0} and \var{p0} define initial guess for flux and pressure. Values marked  \var{u0} and \var{p0} define initial guesses for flux and pressure.
235  by positive values \var{location_of_fixed_flux} and \var{location_of_fixed_pressure}, respectively, are kept unchanged. \var{max_iter} sets the maximum number of iterations steps allowed for solving the coupled problem.  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}  \end{methoddesc}
240
241    %\subsection{Example: Gravity Flow}
242    %later
243
\subsection{Example: Gravity Flow}
later

Legend:
 Removed from v.3301 changed lines Added in v.3326