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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3911 - (hide annotations)
Thu Jun 14 01:01:03 2012 UTC (7 years, 3 months ago) by jfenwick
File MIME type: application/x-tex
File size: 8480 byte(s)
Copyright changes
1 caltinay 3326
2     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3     %
4 jfenwick 3911 % Copyright (c) 2003-2012 by University of Queensland
5 caltinay 3326 % 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 gross 3747 We want to calculate the flux $u$ and pressure $p$ on a domain $\Omega$
17 caltinay 3326 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
45 gross 2960 \subsection{Solution Method \label{DARCY SOLVE}}
46 gross 3502 Unfortunate equation~\ref{DARCY PROBLEM} can not solved directly in an easy way and requires mixed FEM.
47     We consider a few options to solve equation~\ref{DARCY PROBLEM}
48 gross 3568 \subsubsection{Evaluation}\label{SEC DARCY SIMPLE}
49 gross 3502 The first equation of equation~\ref{DARCY PROBLEM} is inserted into the second one:
50     \begin{equation}\label{DARCY PROBLEM SIMPLE}
51     - (\kappa_{ij} p_{,j})_{,i} = f - (g_{i})_{,i}
52 gross 3501 \end{equation}
53 gross 3568
54 gross 3502 with boundary conditions
55     \begin{equation}\label{DARCY BOUNDARY SIMPLE}
56 gross 3501 \begin{array}{rcl}
57 gross 3502 \kappa_{ij} p_{,j} \; n_{i} = ( g_{i} - u^{N}_{i} ) \; n_{i} & \mbox{ on } & \Gamma_{N} \\
58     p = p^{D} & \mbox{ on } & \Gamma_{D} \\
59 gross 3501 \end{array}
60     \end{equation}
61 gross 3568 Then the flux field is recovered by directly setting
62 gross 3502 \begin{equation}\label{DARCY PROBLEM SIMPLE FLUX}
63     u_{j} = g_j - \kappa_{ij} p_{,j}
64     \end{equation}
65 gross 3568 This simple recovery process will not ensure that the (numerically) calculated flux
66     meets the boundary conditions for flux or the incompressibility condition.
67     However this is a very fast way of calculating the flux.
68 gross 3502
69 gross 3568
70 gross 3502 \subsubsection{Global Postprocessing \label{SEC DARCY POST}}
71     An improved flux recovery can be achieved by solving a modified version of equation~\ref{DARCY PROBLEM SIMPLE FLUX}
72     adding the the gradient of the divergence of the flux:
73     \begin{equation}\label{DARCY PROBLEM POST FLUX}
74     \kappa^{-1}_{ij} u_{j} -
75     (\lambda \cdot u_{k,k} )_{,i}=
76     \kappa^{-1}_{ij} g_j- p_{,i}
77     - (\lambda \cdot f )_{,i}
78     \end{equation}
79     where
80     \begin{equation}\label{DARCY PROBLEM POST FLUX A}
81     \lambda = \omega \cdot |\kappa^{-1}| \cdot vol(\Omega)^{1/d} \cdot h
82     \end{equation}
83     with a non-negative factor $\omega$, $d$ is the spatial dimension and $h$ is the local element size.
84     \begin{equation}\label{DARCY PROBLEM POST FLUX BOUNDARY}
85 gross 3501 \begin{array}{rcl}
86 gross 3502 u_{i} \; n_{i} = u^{N}_{i} \; n_{i} & \mbox{ on } & \Gamma_{N} \\
87     u_{k,k} = f & \mbox{ on } & \Gamma_{D} \\
88 gross 3501 \end{array}
89 gross 3502 \end{equation}
90 gross 3568 Notice that the second condition is a natural boundary condition.
91     Global post-processing is more expense than direct pressure evaluation
92     however the flux is more accurate and asymptotic incompressibility
93 gross 3569 for mesh size towards zero can be shown, if $\omega>0$.
94 gross 3502
95 gross 2960
96     \subsection{Functions}
97 gross 3629 \begin{classdesc}{DarcyFlow}{domain, \optional{w=1., \optional{solver=\member{DarcyFlow.POST}, \optional{
98 gross 3502 useReduced=\True, \optional{ verbose=\True} } }}}
99     opens the Darcy flux problem\index{Darcy flux} on the \Domain domain.
100     Reduced approximations for pressure and flux are used if \var{useReduced} is set.
101     Argument \var{solver} defines the solver method.
102     If \var{verbose} is set some information are printed.
103 gross 3568 \var{w} defines the weighting factor $\omega$ for global post-processing of the flux (see equation~\ref{DARCY PROBLEM POST FLUX A}.)
104 gross 2960 \end{classdesc}
105    
106 gross 3568 \begin{memberdesc}[DarcyFlow]{EVAL}
107     flux is calculated directly from pressure evaluation, see section~\ref{SEC DARCY SIMPLE}.
108 gross 3502 \end{memberdesc}
109    
110 gross 3569 \begin{memberdesc}[DarcyFlow]{SMOOTH}
111     solver using global post-processing of flux with weighting factor $\omega=0$, see section~\ref{SEC DARCY POST}.
112     \end{memberdesc}
113    
114 gross 3502 \begin{memberdesc}[DarcyFlow]{POST}
115 gross 3568 solver using global post-processing of flux, see section~\ref{SEC DARCY POST}.
116 gross 3502 \end{memberdesc}
117    
118 gross 2960 \begin{methoddesc}[DarcyFlow]{setValue}{\optional{f=None, \optional{g=None, \optional{location_of_fixed_pressure=None, \optional{location_of_fixed_flux=None,
119     \\\optional{permeability=None}}}}}}
120 caltinay 3326 assigns values to the model parameters. Values can be assigned using various
121     calls -- in particular in a time dependent problem only values that change
122     over time need to be reset. The permeability can be defined as a scalar
123 gross 3568 (isotropic), or a symmetric matrix (anisotropic).
124 gross 2960 \var{f} and \var{g} are the corresponding parameters in~\ref{DARCY PROBLEM}.
125 caltinay 3326 The locations and components where the flux is prescribed are set by positive
126     values in \var{location_of_fixed_flux}.
127     The locations where the pressure is prescribed are set by by positive values
128     of \var{location_of_fixed_pressure}.
129 gross 2960 The values of the pressure and flux are defined by the initial guess.
130 caltinay 3326 Notice that at any point on the boundary of the domain the pressure or the
131     normal component of the flux must be defined. There must be at least one point
132     where the pressure is prescribed.
133     The method will try to cast the given values to appropriate \Data class objects.
134 gross 2960 \end{methoddesc}
135    
136     \begin{methoddesc}[DarcyFlow]{getSolverOptionsFlux}{}
137 gross 3502 returns the solver options used to solve the flux problems.
138 caltinay 3326 Use this \SolverOptions object to control the solution algorithms.
139 gross 3568 This option is only relevant if global postprocesing is used.
140 gross 2960 \end{methoddesc}
141    
142     \begin{methoddesc}[DarcyFlow]{getSolverOptionsPressure}{}
143 caltinay 3331 returns a \SolverOptions object with the options used to solve the pressure
144 gross 3502 problems.
145 caltinay 3331 Use this object to control the solution algorithms.
146 gross 2960 \end{methoddesc}
147    
148 gross 3568 \begin{methoddesc}[DarcyFlow]{solve}{u0,p0}
149 caltinay 3326 solves the problem and returns approximations for the flux $v$ and the pressure $p$.
150     \var{u0} and \var{p0} define initial guesses for flux and pressure.
151     Values marked by positive values \var{location_of_fixed_flux} and
152     \var{location_of_fixed_pressure}, respectively, are kept unchanged.
153 gross 2960 \end{methoddesc}
154    
155 gross 3502 \begin{methoddesc}[DarcyFlow]{getFlux}{p, \optional{ u0 = None}}
156     returns the flux for a given pressure \var{p} where the flux is equal to \var{u0}
157     on locations where \var{location_of_fixed_flux} is positive, see \member{setValue}.
158     Notice that \var{g} and \var{f} are used.
159     \end{methoddesc}
160    
161 gross 3629 \begin{figure}
162     \centerline{\includegraphics[width=\figwidth]{darcy_result}}
163     \caption{Flux and pressure field of the Dary flow example.}
164     \label{DIFFUSION FIG 1}
165     \end{figure}
166 gross 3502
167 gross 3629 \subsection{Example: Gravity Flow}
168     The following script \file{darcy.py}\index{scripts!\file{darcy.py.py}}\index{Darcy flow}
169     which is available in the \ExampleDirectory illustrates the usage of the
170     \class{DarcyFlow} class:
171     \begin{python}
172     from esys.escript import *
173     from esys.escript.models import DarcyFlow
174     from esys.finley import Rectangle
175     from esys.weipa import saveVTK
176     mydomain = Rectangle(l0=2.,l1=1.,n0=40, n1=20)
177     x = mydomain.getX()
178     p_BC=whereZero(x[1]-1.)*wherePositive(x[0]-1.)
179     u_BC=(whereZero(x[0])+whereZero(x[0]-2.)) * [1.,0.] + \
180     (whereZero(x[1]) + whereZero(x[1]-1.)*whereNonPositive(x[0]-1.0)) * [0., 1.]
181     mypde = DarcyFlow(domain=mydomain)
182     mypde.setValue(g=[0., 2],
183     location_of_fixed_pressure=p_BC,
184     location_of_fixed_flux=u_BC,
185     permeability=100.)
186 gross 3502
187 gross 3629 u,p=mypde.solve(u0=x[1]*[0., -1.], p0=0)
188     saveVTK("u.vtu",flux=u, pressure=p)
189     \end{python}
190     In the example the pressure is fixed to the initial presure \var{p0} on the right half of the top face.
191     The normal flux is set on all other faces. The corresponding values for the flux are set by the initial value
192     \var{u0}.

  ViewVC Help
Powered by ViewVC 1.1.26