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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3301 - (hide annotations)
Mon Oct 25 01:27:54 2010 UTC (8 years, 11 months ago) by jfenwick
File MIME type: application/x-tex
File size: 10385 byte(s)
Modified layout of descriptions.
Fixed some overfull hboxes
1 gross 2960 \section{Darcy Flux}
2     \label{DARCY FLUX}
3     We want to calculate the velocity $u$ and pressure $p$ on a domain $\Omega$ solving
4     the Darcy flux problem \index{Darcy flux}\index{Darcy flow}
5     \begin{equation}\label{DARCY PROBLEM}
6     \begin{array}{rcl}
7 jfenwick 3295 u_{i} + \kappa_{ij} p_{,j} & = & g_{i} \\
8     u_{k,k} & = & f
9 gross 2960 \end{array}
10     \end{equation}
11     with the boundary conditions
12     \begin{equation}\label{DARCY BOUNDARY}
13     \begin{array}{rcl}
14 jfenwick 3295 u_{i} \; n_{i} = u^{N}_{i} \; n_{i} & \mbox{ on } & \Gamma_{N} \\
15     p = p^{D} & \mbox{ on } & \Gamma_{D} \\
16 gross 2960 \end{array}
17     \end{equation}
18 jfenwick 3295 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
19 gross 2960 \begin{equation}
20 jfenwick 3295 \alpha_{0} \; x_{i} x_{i} \le \kappa_{ij} x_{i} x_{j} \le \alpha_{1} \; x_{i} x_{i}
21 gross 2960 \end{equation}
22 jfenwick 3295 for all $x_{i}$.
23 gross 2960
24     \subsection{Solution Method \label{DARCY SOLVE}}
25     It is useful to write equation~\ref{DARCY PROBLEM} in operator form. For any pressure $p$
26     we set
27     \begin{equation}
28 jfenwick 3295 (Gp)_{i} = p_{,j}
29 gross 2960 \end{equation}
30     and a velocity $v$ we set
31     \begin{equation}
32 jfenwick 3295 Dv = v_{k,k}
33 gross 2960 \end{equation}
34 jfenwick 3295 We $K=(\kappa_{ij})$ we can write the Darcy problem~\ref{DARCY PROBLEM} as
35 gross 2960 \begin{equation}
36     \begin{array}{rcl}
37     u + K \, Gp & = & g \\
38     Du & = & f
39     \end{array}
40     \end{equation}
41     We solve this equation by minimising the functional
42     \begin{equation}\label{DARCY COST}
43 jfenwick 3295 J(u,p):=\|K^{-\frac{1}{2}}(u + K \, G p - g)\|^2_{0} + \|\lambda (Du-f) \|_{0}^2
44 gross 2960 \end{equation}
45 jfenwick 3295 over all suitable $u$ and $p$. In this equation we set $\|p\|^2_{0}=(p,p)_{0}$ with
46 gross 2960 \begin{equation}
47 jfenwick 3295 (p,q)_{0} = \int_{\Omega } p\cdot q \, dx
48 gross 2960 \end{equation}
49     The factor $\lambda>0$ is a weighting factor.
50     A simple calculation shows that
51     one has to solve
52     \begin{equation}
53 jfenwick 3295 ( K^{-1} (v + K \, Gq) , u +K \, G p - g)_{0} + (\lambda Dv,\lambda (Du-f) )_{0} =0
54 gross 2960 \end{equation}
55     for all velocities $v$ and pressure $q$ which fulfill the homogeneous boundary conditions~\ref{DARCY BOUNDARY}.
56     This so-called variational equation can be translated back into operator notation
57     \begin{equation}
58     \begin{array}{rcl}
59     (K^{-1}+ D^*\lambda^2 D)u + Gp & = & D^*\lambda f + K^{-1} g \\
60     G^*u + G^*K \, G p & = & G^*g \\
61     \end{array}
62     \end{equation}
63 jfenwick 3295 where $D^*$ and $Q^*$ denote the adjoint operators with respect to $(.,.)_{0}$.
64 gross 2960 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
65 jfenwick 3295 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$)
66 gross 2960
67     The approach we are taking is to eliminate the velocity $v$ from the problem. Assuming that $p$ is known we have
68     \begin{equation}\label{DARCY V FORM}
69     v= (K^{-1}+ D^*\lambda^2 D)^{-1} ( D^*\lambda f + K^{-1} g - Gp)
70     \end{equation}
71     (notice that $K^{-1}+\lambda D^*D$ is coercive) which is inserted into the second equation
72     \begin{equation}
73     G^* (K^{-1}+\lambda D^*D)^{-1} (\lambda D^*f + K^{-1} g - Gp) + G^* KG p = G^*g
74     \end{equation}
75     which is
76     \begin{equation}
77     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) )
78     \end{equation}
79     We use the PCG \index{linear solver!PCG}\index{PCG} method to solve this.
80     The residual $r$ is given as
81     \begin{equation}
82     \begin{array}{rcl}
83     r & =& G^* \left( g - K\, G p - v \right)
84     \end{array}
85     \end{equation}
86     for the current pressure approximation $p$ and current velocity $v$ defined by
87     equation~\ref{DARCY V FORM}.
88     So in a particular implementation we use $\hat{r}=g-K\, Gp-v$ to represent the residual.
89     The evaluation of the iteration operator for a given $p$ is then
90     returning $Qp+v$ where $v$ is the solution of
91     \begin{equation}\label{UPDATE W}
92     (K^{-1}+ D^*\lambda^2 D)v = Gp
93     \end{equation}
94     To derive a preconditioner we use the identity
95     \begin{equation}
96     \begin{array}{rcl}
97    
98     G^* ( K - (K^{-1}+ D^*\lambda^2 D)^{-1} ) G & = & G^* (I - (I + K D^*\lambda^2 D)^{-1}) K G \\
99     & \approx & G^* (K D^*\lambda^2 D) K G \\
100     & = & G^* K ( D^* \lambda^2 D) K G \\
101     & \approx & G^* \frac{\lambda^2}{dx^2} K^2 G
102     \end{array}
103     \end{equation}
104     where $dx$ is the local mesh size and we use the approximation
105     \begin{equation}
106     D^* \lambda^2 D \approx \frac{\lambda^2}{dx^2}
107     \end{equation}
108     Therefore we use $G^* \frac{\lambda^2}{dx^2} K^2 G$ as a preconditioner.To evalute the preconditioner
109     we need to solve the equation
110     \begin{equation}\label{UPDATE P}
111     G^* \frac{\lambda^2}{dx^2} K^2 G p = G^* \hat{r}
112     \end{equation}
113     It remains to answer the question how to choose $\lambda$. We need to balance the first and second
114     term in $J(u,p)$ in equation~\ref{DARCY COST}. We inspect $J$ for
115     $(\hat{u}, \hat{p})$ which is a perturbed exact solution $(u,p)$.
116 jfenwick 3295 Assuming $\hat{u}=u+u_{0}e^{ik^tx}$
117     and $\hat{p}=p+p_{0}e^{ik^tx}$ and constant $K$ we get
118 gross 2960 \begin{equation}
119 jfenwick 3295 J(\hat{u},\hat{p}) = C \left[ ( \|K^{-1}\|_{2} |u_{0}|^2 + \|K\|_{2} \|k\|_{2}^2 |p_{0}|^2| )
120     + \lambda^2 \|k\|_{2}^2 |u_{0}|^2 \right]
121 gross 2960 \end{equation}
122     with some constant $C>0$. The first two terms and the third term correspond to the first term and
123 jfenwick 3295 second term in the definition of $J(u,p)$ in equation~\ref{DARCY COST}. For small $\|k\|_{2}$
124     (i.e. for a smooth perturbation) $J(\hat{u},\hat{p})$ is dominated by $\|K^{-1}\|_{2} |u_{0}|^2$.
125 gross 2960 To scale the second term which is corresponds to the incompressibility condition for the velocity
126 jfenwick 3295 we need to meet the condition $\|K^{-1}\|_{2} = \lambda^2 \|k\|_{2}^2$.
127     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
128 gross 2960 \begin{equation}\label{DARCY LAMBDA}
129 jfenwick 3295 \lambda = \|K^{-1}\|_{2}^{\frac{1}{2}} \frac{l}{\pi}
130 gross 2960 \end{equation}
131     Notice that with this setting the preconditioner $G^* \frac{\lambda^2}{dx^2} K^2 G$ becomes
132     equivalent to $G^* K G$ if $K$ is a diagonal matrix and the mesh has a constant mesh size.
133    
134     The residual norm used in the PCG is given as
135     \begin{equation}\label{DARCY R NORM}
136 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
137 gross 2960 \int \hat{r} K^{-1} \hat{r} \; dx
138     \end{equation}
139     The iteration is terminated if
140     \begin{equation}\label{DARCY STOP}
141 jfenwick 3295 \|r\|_{PCG} \le \mbox{ATOL}
142 gross 2960 \end{equation}
143     where we set
144     \begin{equation}\label{DARCY ATOL DEF}
145 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}
146 gross 2960 \end{equation}
147     where rtol is a given relative tolerance and $\mbox{atol}$ is a given absolute tolerance (typically $=0$).
148     Notice that if $Gp$ and $v$ both are zero, the pair $(0,p)$ is a solution.
149     The problem is that ATOL is depending on the solution $p$ and $v$ calculated form~\ref{DARCY V FORM}.
150     In practice one use the initial guess for $p$
151     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.
152    
153     \subsection{Functions}
154     \begin{classdesc}{DarcyFlow}{domain}
155     opens the Darcy flux problem\index{Darcy flux} on the \Domain domain.
156     \end{classdesc}
157    
158     \begin{methoddesc}[DarcyFlow]{setValue}{\optional{f=None, \optional{g=None, \optional{location_of_fixed_pressure=None, \optional{location_of_fixed_flux=None,
159     \\\optional{permeability=None}}}}}}
160     assigns values to the model parameters. Values can be assigned using various calls - in particular
161     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).
162     \var{f} and \var{g} are the corresponding parameters in~\ref{DARCY PROBLEM}.
163     The locations and components where the flux is prescribed are set by positive values in
164     \var{location_of_fixed_flux}.
165     The locations where the pressure is prescribed are set by
166     by positive values of \var{location_of_fixed_pressure}.
167     The values of the pressure and flux are defined by the initial guess.
168     Notice that at any point on the boundary of the domain the pressure or the normal component of
169     the flux must be defined. There must be at least one point where the pressure is prescribed.
170     The method will try to cast the given values to appropriate
171     \Data class objects.
172     \end{methoddesc}
173    
174     \begin{methoddesc}[DarcyFlow]{setTolerance}{\optional{rtol=1e-4}}
175     sets the relative tolerance \mbox{rtol} in \ref{DARCY ATOL DEF}.
176     \end{methoddesc}
177    
178     \begin{methoddesc}[DarcyFlow]{setAbsoluteTolerance}{\optional{atol=0.}}
179     sets the absolute tolerance \mbox{atol} in \ref{DARCY ATOL DEF}.
180     \end{methoddesc}
181    
182     \begin{methoddesc}[DarcyFlow]{getSolverOptionsFlux}{}
183 jfenwick 3301 Returns the solver options used to solve the flux problems~(\ref{DARCY V FORM}) and~(\ref{UPDATE W}). Use this
184     \SolverOptions object to control the solution algorithms.
185 gross 2960 \end{methoddesc}
186    
187     \begin{methoddesc}[DarcyFlow]{getSolverOptionsPressure}{}
188     Returns the solver options used to solve the pressure problems~(\ref{UPDATE P}) as a preconditioner.
189 jfenwick 3301 Use this \SolverOptions object to control the solution algorithms.
190 gross 2960 \end{methoddesc}
191    
192     \begin{methoddesc}[DarcyFlow]{solve}{u0,p0, \optional{max_iter=100, \optional{verbose=False}}}
193     solves the problem. and returns approximations for the flux $v$ and the pressure $p$.
194     \var{u0} and \var{p0} define initial guess for flux and pressure. Values marked
195     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.
196     \end{methoddesc}
197    
198    
199     \subsection{Example: Gravity Flow}
200     later

  ViewVC Help
Powered by ViewVC 1.1.26