/[escript]/branches/doubleplusgood/doc/user/darcyflux.tex
ViewVC logotype

Contents of /branches/doubleplusgood/doc/user/darcyflux.tex

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26