A MIXED FINITE VOLUME ELEMENT METHOD
FOR ACCURATE COMPUTATION OF FLUID
VELOCITIES IN POROUS MEDIA
by
Jim E. Jones
B.S., New Mexico Tech, 1986
M.S., New Mexico Tech, 1989
A thesis submitted to the
Faculty of the Graduate School of the
University of Colorado at Denver
in partial fulfillment
of the requirements for the degree of
Doctor of Philosophy
Applied Mathematics
1995
This thesis for the Doctor of Philosophy
degree by
Jim E. Jones
has been approved
Stephen F. McCormick
Thomas A. Manteuffel
John W. Ruge
Thomas F. Russell
John Trapp
Date
Jones, Jim E. (Ph.D., Applied Mathematics)
A Mixed Finite Volume Element Method for Accurate Computation of
Fluid Velocities in Porous Media
Thesis directed by Professor Stephen F. McCormick
ABSTRACT
A key ingredient in the simulation of flow in porous media is the accurate
determination of the velocities that drive the flow. Largescale irregularities
of the geology, such as faults, fractures, and layers, suggest the use of irreg
ular grids in the simulation. In this study, the approach to this problem was
to apply the finite volume element methodology, developed by McCormick,
in conjunction with mixed methods, developed by Raviart and Thomas. The
resulting mixed finite volume element discretization scheme developed here
can be applied in a clear and straightforward way to irregular grids and
is appealing because of its local conservation properties and its direct and
accurate representation of physical intercell flux terms. Several multilevel
algorithms are developed that provide efficient methods for solving the set of
equations that this discretization produces. This thesis includes numerical
results from a variety of test problems, from Poissons equation to prob
lems with anisotropic, discontinuous, and tensor diffusion coefficients. These
results show that this approach has the potential to generate accurate ap
proximate fluid velocities and that the multilevel methods can provide fast
solvers.
This abstract accurately represents the content of the
candidates thesis. I recommend its publication.
Signed
Stephen F. McCormick
IV
Contents
Chapter
1 Uniform Grids .................................................1
1.1 The Mixed Finite Volume Element Discretization............1
1.2 An FVE Based Multigrid Algorithm .......................13
1.2 Computational Results ...................................25
1.2.1 Comparison to Finite Differences for Poisson Problem ... 25
1.2.2 Discontinuous and Anisotropic Diffusion Coefficients .... 30
1.2.3 An Alternative Multilevel Algorithm..............35
1.2.4 Full Tensor Diffusion Coefficients ..............42
Chapter
2 General Quadrilateral Grids..................................46
2.1 The Mixed FVE Discretization ............................46
2.2 Multilevel Solvers ......................................55
2.2.1 FVE Based Multigrid..............................55
2.2.2 A TwoLevel Approach ............................60
2.3 Computational Results ...................................62
2.3.1 Multigrid Performance............................62
v
2.3.2 TwoLevel Performance................................67
2.3.3 Accuracy of the Mixed FVE Discretization ............70
Chapter
3 Technical Derivations ...........................................77
3.1 Uniform Grids and Tensor Permeability.......................77
3.2 Derivation of FVE Stencils for General Quadrilaterals .....80
3.2.1 Scalar Diffusion Coefficient ....................80
3.2.2 Tensor Diffusion Coefficient ....................89
3.3 Using Fine Grid Integrals to Calculate Coarse Grid Integrals .... 91
Chapter
4 Summary and Future Work .........................................97
Bibliography.......................................................100
vi
ACKNOWLEDGEMENTS
This thesis presents the results from a fiveyear period of research, and
it would not be complete without a word of special thanks to those whose
assistance enabled me to complete it.
First of all, I wish to thank Professor Steve McCormick for giving me the
opportunity to work on a host of interesting research projects. I also thank
Professor McCormick for the wide freedom that I was given as a graduate
student to do independent research, and his patience as I labored to bring
this research to its conclusion.
I also wish acknowledge Professor Thomas Russell for his assistance with
this work. Without his expertise in the held of reservoir simulation, this
research could not have taken place. I also wish to thank Professor John
Ruge and Professor Steve Schaffer. While they were not directly involved in
the daytoday work on this project, I learned most of what I know about
multigrid methods through collaborative research with these two men.
I would also like to thank Dr. Pat Roache, Ecodynamics Research, and
Sandia National Laboratory for their financial support at various times dur
ing this research.
Finally, I wish to acknowledge the countless sacrifices made by my family
 my wife Dorothee, my daughter Frederika, and my yet unnamed second
child that allowed me to complete this thesis. Most of all I thank my wife
for her unfailing love and support during this time.
vii
1 Uniform Grids
1.1 The Mixed Finite Volume Element Dis
cretization
In this first section, we introduce the new mixed finite volume element dis
cretization technique in a familiar context: solving the diffusion equation on
a uniform, square grid. However, it should be clear from the presentation
that the discretization developed in this section, and the multigrid algorithm
constructed in the next section, can easily be modified to apply to any tensor
product grid.
We begin by considering the following partial differential equation defined
on some square coordinatealigned domain 0 in 1Z2:
f V A(x)V(x) =/(x), X E o,
[ A(x)V(x) rj = g(x), x G <90. ' '
In the context of reservoir simulation, this is the pressure equation for in
compressible singlephase flow, where is the pressure in the reservoir 0
and the boundary condition specifies the flux on <90. Because one of our
1
goals for the new discretization is accurate approximations of flow velocities,
we will begin by reformulating this equation as a firstorder system of equa
tions, where velocity appears explicitly. This is done by introducing the flow
velocity variables via the definition
v = AV(f>, (1.2)
then rewriting the partial differential equation in (1.1) as
V v = /. (1.3)
In the context of reservoir simulation, definition (1.2) is Darcys law and
equation (1.3) is the mass conservation law. In reservoir simulation, this
same approach of treating flow velocity explicitly has been used in mixed
Uniteelement methods with considerable success [11],[12],[18].
To continue the derivation of the discretization, assume that A is a diag
onal tensor,
Here, ax and ay may be functions of position. Section 3.1, contains the details
of the discretization process for the case where A is a full tensor. Let u and
v denote the components of the velocity,
(it, vf = v.
Then the secondorder partial differential equation (1.1) may be written in
2
the form of a firstorder system in 0:
ax~1u + d(f)jdx = 0 (u equation)
< ay~xv + d(f)/dy = 0 (v equation)
du/dx + dv/dy = / (p equation),
with boundary conditions,
u(x) fi'(x)5 x ^^west U ^^east
u(x) H(x)*> x ^^south ^ '^nortlr
The labels u,v, and p for the three equations in equation (1.4) are introduced
simply for convenience.
To discretize this system, we follow the finite volume element (FVE)
principles developed in [7],[13],[14], The two major components of any FVE
discretization scheme are the choice of control volumes over which to inte
grate the continuous equation and the choice of finite element spaces for the
unknowns.
First, to introduce the control volumes for the FVE discretization, con
sider a uniform square mesh Qh with mesh size h that covers 0. We introduce
three sets of control volumes, one for each of the three equations in (1.4).
Examples of the control volumes are shown in Figure 1.1. We denote by hi
the set of all volumes U that will be used to discretize the u equation in (1.4).
Similarly, we use the notation V and V for the sets of volumes V and P for
the v and p equations, respectively. Note that the control volumes P corre
spond to the grid blocks defined by Qh, the control volumes U straddle the
vertical grid edges of Qh, and the control volumes V straddle the horizontal
grid edges of Qh.
3
Figure 1.1: Uniform Grid
4
Next we will introduce the finite elements for the FVE discretization. For
our finite element space, we use the lowest order RaviartThomas elements
[16],[19] on the the control volumes V:
uh is linear in x and constant in y on each P Â£ P,
vh is linear in y and constant in x on each P Â£ P,
Note that uh is continuous across vertical edges of the grid and that vh is
continuous across horizontal edges of the grid. The location of the nodes for
each of the unknowns with their indexing is also shown in Figure 1.1. The
flux boundary condition in equation (1.1) specifies the values of uh on the
west and east boundaries of 0 and the values of vh on the north and south
boundaries of 0. Noting some of the characteristics of the finite element
spaces, the pressure space contains all functions that are piecewise constant
on Qh. A typical basis function for uh is shown in Figure 1.2.
Note that the normal component of this basis function is continuous across
all grid edges, is zero on all edges except the edge shared by the two P
volumes which are the support of the function where the normal component
equals a nonzero constant. With these properties we can guarantee that the
computed flow velocity will also have continuous normal components across
grid edges. This will also be the case for general quadrilateral grids as we
will see in Section 2.1. The true physical solution also has this property 
continuous normal component of velocities but not every numerical scheme
for approximating it does, as [17] pointed out.
5
Figure 1.2: u Basis Function
We are now ready to discretize the equations. Start by taking the u
equation in (1.4) and integrating it over each U Â£ U. This will result in a
discrete version of Darcys law. As an example, let U8J be the volume in hi
that is centered at the interior uh node (i,j), as shown in Figure 1.3.
We then have
/ (ax ru + dcfi/dx^j dxdy = 0. (1.5)
Consider the two terms in the integral separately. For the first term in
equation (1.5), we split the integral into two parts, one for each half of U4j.
This yields,
/ (ax dxdy = / (cix 1u^ dxdy + / (cix 1u^ dxdy, (1.6)
t/ XJ 7 7 t/ XT ^ t/ XT ^
6
Figure 1.3: UVolume
7
where U* and U[? denote the left and right halves of the control volume.
Now because the finite element approximation uh is linear in x and constant
in y on each of the halves, by direct calculation it follows that
J t dxdy = a^1^)
A similar expression is obtained for the integral over the right half of U4j.
Note that here we have assumed a^1 to equal the constant value a:E_1(/) on
the left half of the control volume. If it is instead piecewise constant on some
finer mesh, as will often be the case in the section on computational results,
then the integral can be split into parts on this finer mesh and evaluated.
The details of this piecewise integration can be found in Section 3.3, where
we show, in the more general case of quadrilateral grids, how to use integrals
calculated on fine grids to generate integrals on coarser grids. Now for the
second term in equation (1.5), we have
' U.
(dcf)/dx) dxdy = / / (d/dx) dx dy.
Then apply the divergence theorem to get
f (f [d(h,y) (0,y)) dy.
(1.7)
Then, as the finite element approximation
x = 0 and x = h, it follows that
((h,y) ~ (0,y))dy = h
hJ
,j
(1.8)
8
Putting these observations together, the discrete equation that equation (1.5)
gives rise to is thus
Y \ax 1(?) (3+ Ui+itj) + ax 1(l) (3ufj +
{^,3 ~ $1,3) = O'
If ax~x is constant, the equation simplifies to
1
X
uUi + K,
"I u
*+i j
Note that this expression is quite similar to the expression one would get
using standard cellcentered finite differences to approximate Darcys law:
hi
(1.9)
The difference between the two expressions is that, with the mixed FVE
method, tridiagonal mass matrix multiplies the velocities rather than a triv
ially invertible diagonal matrix. One way of viewing equation (1.9), is as
a midpoint rule approximation to the integral in equation (1.6). In Sec
tion 1.3, we will show that this can result in less accurate approximations for
the velocities. More importantly, as we will see when we consider irregular
grids, the mixed FVE method has the property of generating discrete ver
sions of Darcys law for geometries where it is unclear how to use standard
cellcentered finite differences.
Integrating the v equation in (1.4) over an interior V volume yields a
similar discrete expression involving nodal values of vh and
9
To complete the discretization process, integrate the p equation in equa
tion 1.4 over the volume in V centered at the interior
denote this volume by P4j. See Figure 1.4.
Integration is simply a matter of applying the divergence theorem. From
the definition of u and u, we have
/ (du/dx + dv/dy) dxdy = / Vvdxdy.
Transforming the volume integral to a surface integral yields the following:
fPij V vdxdy = fap.'. v yds =
= fgpe udy fdpw udy + fgpn vdx fgps vdx}
where dPp denotes the east edge of the control volume, and similarly for the
west, north, and south edges. Then, since the finite element approximations
uh and vh are constant on edges, we have
fgpe udy fgpw udy + fgpn vdx fgps vdx
h(u
V.J l,J l,J
h o,h h h \
t + l,3 U+3 ^ Fj + 1 Fjh
The discrete version of the p equation in (1.4) is therefore just
h(ui+i,j ~ uh + vti+1 vt) = h2kr
Here, on the right hand side, we have assumed that / is a piecewise constant
function on V, with /8J denoting the value of / at the
way of looking at this is: we have used the midpoint approximation to the
integral; if we had assumed that we had more information about /, we could
have used a more accurate approximation for its integral. For example,
10
Figure 1.4: P Volume
11
knowing / as a function of x and y, we could have used a Gauss quadrature
rule to approximate the integral.
In summary, the discretization has produced discrete versions of Darcys
law for each volume U Â£ hi and each volume V Â£ V. These equations relate
the pressure drop between cells to a linear combination of velocities. The dis
cretization has also produced discrete conservation equations for each volume
P G V. For each P volume, we have that the discrete flux out of (or into) the
volume balances with the strength of the source (or sink) in the volume. This
conservative property is an important feature of FVE discretization schemes
in general. The mathematical modeling of many physical systems (including
flow in porous media) starts with a conservation law stated in integral form.
In going from the continuous mathematical model to the discrete mathe
matical model, preserving some analogy of this conservation law is desirable.
For time dependent applications, this discrete conservation property is often
necessary to generate acceptable solutions.
The discretization technique described in this section can easily be gener
alized to handle problems in three dimensions. Also, it can be easily modified
to handle the Dirichlet problem, where the flux boundary condition in equa
tion (1.1) is replaced by
In [10] this mixed FVE discretization is applied to Poissons equation with
Dirichlet boundary condition on a three dimensional domain. This reference
12
contains the details of the necessary modifications.
1.2 An FVE Based Multigrid Algorithm
We begin with a very brief overview of multigrid fundamentals, assuming that
the reader is already somewhat familiar with the held. Good references are
the introductory tutorial [6] and the advanced guide [3]. Another useful
reference is [14], which presents multigrid (and multilevel adaptive) methods
for use in conjunction with FVE discretization techniques.
Consider the problem of finding the solution to the discrete set of equa
tions
Lhuh = fh. (1.10)
We think of equation (1.10) as coming from the discretization of some contin
uous partial differential equation on a uniform grid with mesh size h, although
multigrid methods can be applied to discrete problems with no continuous
origin (e.g., [4]). If the number of unknowns is too large to solve the equation
directly, one might try simple GaussSeidel relaxation to find an approximate
solution, uh, to (1.10). One typically finds that the convergence speed of this
process degrades very early in the iterations: the first few relaxation sweeps
may substantially reduce the residual
h rh t h A h
r = / L u ,
but, subsequent relaxation sweeps do not. The problem is that relaxation
methods generally do not reduce lowfrequency errors very well. In multigrid
13
methods, coarser grids are used to solve for and then eliminate this low
frequency error that relaxation alone cannot resolve. This is done by way of
the error equation
Lheh = r\ (1.11)
where eh denotes the error in the current approximation:
eh = uh uh.
Because relaxation does work well for highfrequency error, after a few relax
ation sweeps one can represent the error accurately on a coarser grid. With
this in mind, we consider a coarse grid version of equation (1.11):
L2hu2h = I2hhrh. (1.12)
Here we think of the operator L2h as coming from the discretization of the
same continuous partial differential equation on a uniform grid, but with
twice the mesh size as the original problem (1.10). The operator I2h is called
the restriction operator and is used to transfer the fine grid function rh to
the coarser grid. One may be able to solve equation (1.12) exactly for u2h,
remembering that this system is on a coarser grid so it is not as large as the
original one. Then one would correct the fine grid approximation by u2h,
which is approximately equal to the error in uh, according to
uh < uh + I^hu2h. (1.13)
Here, the operator 7^, called the interpolation operator, is used to transfer
the coarse grid function u2h to the fine grid. We now perform a few relaxation
14
sweeps to reduce any highfrequency error that we may have introduced. This
completes one multigrid cycle, the steps of which were:
1) Relax the fine grid problem.
2) Form the coarse grid version of the error equation.
3) Solve the coarse grid version of the error equation.
4) Correct the fine grid approximate solution.
5) Relax the fine grid problem.
In practice, one may not be able to solve the coarse grid version of the error
equation for the same reasons one cannot solve the original problem: it is
still too large, so that relaxation alone does not work and direct methods
are too expensive. One then can use the above multigrid cycle recursively on
this problem, using a still coarser grid with mesh size Ah. This process can
be continued until a coarse enough grid is reached, where a direct method or
relaxation by itself is effective. One multigrid cycle is represented graphically
in Figure 1.5.
The basic components of a multigrid algorithm are the grids themselves,
the discrete operators on each grid (Lh, L2h, L4h,. .), the relaxation process,
and the grid transfer operators (interpolation and restriction). The multigrid
method is an iterative method, but typically one needs only a few cycles to
produce a good approximation. Moreover, the socalled full multigrid algo
rithms, which start with multigrid cycles applied to coarser levels, can achieve
approximations with accuracy comparable to the exact discrete solution in
the equivalent of just two or three fine grid cycles.
Now to develop a multigrid algorithm designed specifically for solving
15
Figure 1.5: One Multigrid Cycle
16
the discrete set of equations generated by the mixed FVE discretization of
the previous section, consider a family of uniform square grids Qh of mesh
size h that cover the region 0. Figure 1.6 shows a coarse grid with
twice the mesh size of the grid Qh in Figure 1.1. Note that the coarse grid
looks like the fine grid. Simply joining four square fine grid cells produces
a square coarse grid cell (for the general quadrilateral grids in Section 2.2
producing coarse grids from fine grids is not so simple). On each grid Qh, we
can apply the mixed FVE discretization process, and we write the discrete
set of equations that this process generates as
Lhzh = gh} (1.14)
where zh = (uh,vh,hy and gh = with the unknowns uh, vh,
and on
grid Qh. Note that, in the FVE discretization of the previous section,
and were zero. In the multigrid algorithm, on coarse grids these will,
in general, be nonzero when we form the coarse grid version of the error
equation.
We now must define interpolation and restriction operators and a relax
ation process. For defining interpolation operators, we use the same general
principles as outlined in [14], but we must be specific here. The finite vol
ume element discretization is based on finite element spaces for the variables
uh,vh, and
ment spaces of the different grids to define interpolation. In particular, the
17
Figure 1.6: Uniform Coarse Grid
18
finite element spaces are nested in that one can write a coarse grid basis
function as a linear combination of fine grid basis functions. This nesting
property simplifies the definition of interpolation operators. To define the
interpolation operator for , note that (f)2h is constant on the grid 2h
volume P/,j. Referring to Figure 1.7, we thus have the following characteri
zation of h = I{(f))^h(f)2h:
bh
hj
ih
^'+1 ,j
ih
^'j + 1
ih
(1.15)
To define the interpolation operator for it, /(n)^, note that u2h is linear in
x and constant in y on the grid 2h volume P/,j. We thus have the following
characterization of uh = (see Figure 1.7):
ui,j = ui,j+i =
Ui+2,] = Ui+2J + 1 =
(1.16)
i+lj Ui+lj + 1 1/2(mJ+1,J + M/,j)
Dehnition of the interpolation operator for v is similar.
For defining restriction operators, we again use the same general prin
ciples as outlined in [14]. In the multigrid algorithm, restriction operators
are used to transfer right hand sides and residuals of equations, not the un
knowns themselves. The definitions of the restriction operators are therefore
naturally based on the relationship between the volumes of the various grids.
The idea is to lump several of the grid h right hand sides to produce the grid
2h right hand sides. To define the restriction operator for the P equation,
J(P)2/l, note that a grid 2h volume P/,j consists of four of the P volumes
19
Figure 1.7: Fine and Coarse Grids, I
20
of grid h. We thus have the following characterization of f2h = I(P)2hhfh,
referring again to Figure 1.7:
r2h rh rh
JI,J Ji,j '
+ fh+1
+ /i+l,j+l
(1.17)
To define the restriction operator for the U equation, I(U)2hh, note that a
grid 2h volume Uyj of Q2h contains all of two of the U volumes of grid
h and half of four others. We thus have the following characterization of
f2h = I(U)2hhfu, referring to Figure 1.8:
1
flkI,J ~ fui,j + fui,j +1 + + + l + fui+l,j + fut + ij + i) l118)
Dehnition of the restriction operator for the V equation is done in a similar
fashion.
For the equations on grid Qh, there are a variety of relaxation techniques
one could use. We will present the two that have proved to be most effi
cient in our numerical experiments: distributive GaussSeidel relaxation and
box GaussSeidel relaxation. We consider these two relaxation schemes be
cause there is no onetoone correspondence between equations and unknowns
(point GaussSeidel makes no sense per se). There are three types of discrete
equations:
U equations:
\ax1 + 6ui,j + Ui+l,j) + ^
V equations:
1 (Gji + vi,j + Gj'+i) + h
P equations:
h(ui+i,j  + uG+1 VP) = h2fitj.
kh
h1 j
0
Ml
= fh. .
J U % J 5
fh. .
J V % J 5
(1,19)
21
Figure 1.3: Fine and Coarse Grids, II
22
Again, we point out that the labels are used only for convenience; for example,
there is no inherent reason to associate the P equation closely with the
variable and, in fact, it does not even appear in this equation.
The guide [3] contains an excellent description of distributive Gauss
Seidel relaxation schemes, so that here we will be content with a brief de
scription of our scheme. We can naturally associate the U equations with
the velocity variable uG of the grid edge that the U volume straddles. Like
wise, we can naturally associate the V equations with the velocity variable
vG of the grid edge that the V volume straddles. We therefore relax these
equations using standard point GaussSeidel relaxation. The basic idea is
to relax the U, V, and P equations together locally in such a way that the
error is smoothed while the U and V equations are not damaged, i.e., their
residuals are essentially unchanged. We can think of the overall relaxation
process (including point GaussSeidel for U and V) as a three step scheme.
First, we sweep over all of the uh nodes, change the value of in turn
so that the U equation at (i,j) is satisfied. Second, we perform a similar
GaussSeidel relaxation of all of the V equations. Note that these two steps,
U and V relaxation, are independent of each other and could be performed
in parallel. In fact, for full independence, point Jacobi could be used within
each variable without much loss of efficient smoothing because of the associ
ated tridiagonal mass matrix. Finally, we sweep over each
simultaneously changing the value of G and the values of uh and vh that
23
lie on the edge of the volume P8J (namely, , rt^+1 and uG+1) so that
the P equation at (i,j) is satished and so that the residuals of the U equa
tions at (i,j) and (i + 1, j) and of the V equations at (i,j) and (z, j + 1) are
unchanged. To allow for effective vectorization, the GaussSeidel relaxation
process in each of these three steps is done using a red/black ordering.
The second technique we consider is box relaxation. This approach in
volves only one step, which consists of a sweep over each P volumes in turn
changing the values of the associated five unknowns ,u^,u^+1 ,and
vf +i so that the corresponding five equations are satished. We will often
need to use a yline box GaussSeidel relaxation technique, which is the same
as box relaxation except that all of the P volumes sharing the same z index
are relaxed simultaneously as a group. Analogously, xline box GaussSeidel
involves relaxing all P volumes that share the same j index as a group. We
will use the term xy or alternating line relaxation to signify a single relax
ation step where we perform xline relaxation followed by yline relaxation.
As will be evident from the results of numerical experiments the best choice
of relaxation technique, in the sense of the most efficient solver, is problem
dependent. The line relaxation techniques are clearly computationally more
expensive; but for problems with anisotropic coefficients, line relaxation is
often needed to achieve acceptable convergence rates for the multigrid algo
rithm.
24
1.3 Computational Results
The major motivation in developing the mixed FVE discretization approach
was to have an accurate discretization for irregular grids, where there is no
clear way to use standard cellcentered finite differences. However, a case
can be made for using the mixed FVE discretization even on uniform grids,
as we will demonstrate in this section.
1.3.1 Comparison to Finite Differences for Poisson
Problem
We begin with the simplest of problems of the type to which our discretization
rpplies, namely Poissons equation on the square domain 0 = [1, l]z
f v v(x) = /(x) x e o,
) V(x) r] = g(x) x G <90
(1,20)
One can write this equation as a system of firstorder equations in the way
presented in Section 1.1, then apply the mixed FVE discretization. We
will call this the mixed FVE approach. Alternatively, one could discretize
equation (1.20) using standard cellcentered finite differences, then solve the
resulting discrete set of equations for an approximation to . If we were
interested in approximations to velocities or fluxes, we could then apply ap
propriate standard finite differencing to the approximation of . We will call
this the finite difference approach.
To compare these two approaches, consider the function
(x, y) = cos(kiTrx) + cos(k2Try)}
(1,21)
25
where k\ and k2 are integers. If we set
f(x, y) = 7T2(kf + ^2)(cos(A;i7ra;) + cos(k2Try))}
g(x,y) = 0
in equation (1.20), then the function in (1.21) is a solution to the partial
differential equation. It is not the only solution because we can add an
arbitrary constant to the solution , while preserving the property of satis
fying the partial differential equation and boundary condition. The function
in (1.21) is the unique solution that satisfies the additional condition that
its integral over 0 is zero. By varying k\ and k2} one can see the effect that
oscillations in the true solution have on the accuracy of the numerical solu
tions generated by a mixed FVE approach and a finite difference approach.
We measure the error in the solutions in the following way. For the mixed
FVE approach, we calculate the C2 norm of (f)m^ve and call this quantity
^mfve' gecause the mixed FVE discretization relies on finite element spaces,
one has an approximation to (j) everywhere in 0 except at grid interfaces, so
the C2 norm make sense. For the straightforward finite difference approach,
one has approximation to ^ only at the nodes. To create a functional approx
imation to , we assume that the nodal approximation holds for the entire
cell, then calculate the C2 norm of its difference with the exact solution. We
call this quantity qBecause we are interested in accurate approximation
of velocities, which are here just the components of the gradient of , we
measure the error in them as well, which we do in the following way. For the
26
mixed FVE approach, we calculate the difference between the fluxes,
s
using the exact and the approximate solutions over each grid edge. We
then take the Â£2 norm of the difference on vertical edges and call this u]Â£^ve.
Analogously, we take the Â£2 norm of the difference on horizontal edges and
call this vfve. Note that if the exact velocity is V and the mixed FVE
approximation is v, then ((uenrl^ve)2 + (vem^ve)2)^ is equivalent to a discrete
H(div) norm of the vector velocity error (which incorporates C2 norms of
(v V)^, (v V)^, and div(v V), the last of which is zero by the local
conservation property of the mixed FVE method). For the finite difference
approach, we do the same to calculate v/ed and vÂ£d. To get approximations
for v, we use standard finite differencing of the approximation of .
Tests were run on a grid with 64 cells in each direction and the results
are shown below.
Mixed FVE
h k2 (j/n fve urnfVe Vmfve ue
l 1 4.01E2 2.52E6 2.52E6
l 4 1.17E1 1.74E2 4.38E3
l 8 2.27E1 7.14E2 9.16E3
l 12 3.37E1 1.44E1 1.28E2
l 16 4.41E1 2.10E1 1.49E2
15 13 5.27E1 2.59E1 2.92E1
27
Finite Difference
h k2 <
l 1 4.01E2 3.89E12 3.89E12
l 4 1.17E1 3.56E2 8.96E3
l 8 2.32E1 1.57E1 2.02E2
l 12 3.53E1 3.62E1 3.20E2
l 16 4.84E1 6.55E1 4.51E2
15 13 5.46E1 8.79E1 9.92E1
One can see that both methods produce approximations to of roughly
equal quality. When the solution is smooth, the approximations to the fluxes
are also of equal quality. However, for oscillatory solutions, i.e., increasing
&2 (k2 = 16 is the most oscillatory function that can be represented on this
grid), the mixed FVE obtains somewhat more accurate approximations to
the fluxes. This is also true for the last case where the solution oscillates
in both directions. This is not surprising considering that one can view the
finite difference approach as a midpoint rule approximation of the integral
in the Darcy equation, as mentioned in Section 1.1. This approximation is
accurate for smooth functions, but not so accurate for more oscillatory ones.
To be fair in comparing these two approaches, one should also consider
the work necessary to generate the approximations. For solving the equations
in the finite difference approach, we used black box multigrid [8]. This code
is aimed at solving much more difficult problems than Poissons equation, as
we will say more about later, and is therefore not the most efficient code for
such a simple problem. For solving the mixed FVE equations, we used the
multigrid algorithm of the previous section with distributive relaxation. Both
28
multigrid methods proved to have similar worstcase convergence factors,
approximately .1, so that the amount of work for solving the mixed FVE
equations would be about three times that for the finite difference equations.
Loosely speaking, the factor of three comes from the fact that, in the mixed
FVE approach, one has roughly three times the number of equations as in
the finite difference approach. A careful comparison of relative costs of the
two algorithms should compare relaxation complexity: while the mixed FVE
scheme is more involved, its individual stencils are a little smaller, so that
the factor of three is still approximately correct. If one were to compare
fully optimized codes, we are confident in our belief that the mixed FVE
solver would be no worse than three times slower than the finite difference
solver. For oscillatory solutions we have seen that the error in the mixed FVE
fluxes can be up to three times smaller than the error for the finite difference
approach. So one might argue that the comparison comes out very roughly
even. However, these tests certainly show that the mixed FVE approach has
the potential to generate more accurate approximations for fluxes than the
standard finite difference approach even for simple problems, and for slightly
more complex problems, the difference in accuracy between the two methods
can be much greater as will be shown in the next section.
29
1.3.2 Discontinuous and Anisotropic Diffusion Coef
ficients
The results for Poissons equation, while encouraging, are not of particular
interest from the reservoir simulation point of view. In reservoir simulation,
the diffusion coefficient A is typically anisotropic and may have large jumps
from grid cell to grid cell. We will now consider these types of problems.
It has been known for some time that to get acceptable multigrid conver
gence factors for anisotropic problems, one needs a relaxation scheme more
complicated than simple point (or cell) relaxation. The reason is that simple
point GaussSeidel relaxation will smooth the error only between strongly
coupled points. For example, consider the diffusion equation
V AV(x) = /(x), x G 0,
where
Standard finite differences and a standard multigrid algorithm with point
GaussSeidel relaxation leads to extremely poor convergence factors when
e
points in the x direction than in the y direction. Local mode analysis [3]
shows that point GaussSeidel relaxation will only smooth the error effec
tively in the x direction. Therefore, one cannot properly represent the error
on a grid that is coarser in the y direction. One alternative is to accept
30
this smoothing limitation and only coarsen in the x direction, an approach
that has been used in reservoir simulation in [9]. Another alternative is to
relax as a group all of the points that are strongly coupled, which, in this
case, consists of those points that share the same y coordinate. This xline
relaxation will smooth the error in the weakly coupled direction (the y direc
tion). Then one can represent the error well on a grid that is coarser in the
y direction. The distributive relaxation scheme as previously discussed is a
single cell relaxation scheme. Only unknowns associated with a particular
cell are updated together, and we cannot hope to get acceptable multigrid
convergence rates with it for anisotropic problems. The alternative we choose
is line relaxation as discussed in Section 1.2. It may often be the case that
the strong coupling is in the x direction in parts of the domain and in the y
direction in other parts. For these problems, we will need to use alternating
line relaxation. An alternative, that we will not consider in this thesis, would
be to use semicoarsening in one of the directions and line relaxation in the
other as in [9].
Later we will discuss in more detail problems with the multigrid algorithm
when the diffusion coefficient jumps from cell to cell. However, the algorithm
as presented in the previous section yields quite acceptable convergence rates
when the jumps in the diffusion coefficient occur (if at all) at grid interfaces
on the coarsest grid.
We will now present a single numerical example that illustrates the abil
31
ity of the multigrid solver to deal with both anisotropic and discontinuous
diffusion coefficients. Consider
/ 1 ax u + dc^/dx = 0
ay~1v + d(f>/dy = 0
du/dx + dv/dy = /
with boundary conditions
u(x) = fi'(x)5 x C <9fiwest U <90east
u(x) 9(x)x C ^^south ^ '^nortlr
The domain of the problem, along with values of ax and ay, is shown in Figure
1.9. Note that this problem has regions of strong coupling in the x direction
and regions of strong coupling in the y direction. The diffusion coefficient
also has lines of discontinuity, although they correspond to grid edges on the
coarsest (2 by 2) grid. To obtain asymptotic estimates, the right hand side
/ and the boundary condition g were set to zero, a random initial guess was
used, and ten or more multigrid cycles were performed. We then calculated
the geometrically averaged convergence factor ( factor by which the residuals
are reduced per cycle). The results are shown below.
Grid Size Convergence Factor
8x8 .04
16 x 16 .10
32 x 32 .18
64 x 64 .20
Alternating line relaxation was used with two prerelaxations and one
postrelaxation. We see acceptable convergence factors in that they are fairly
small and seem to stabilize as the number of unknowns is increased.
32
Figure 1.9: Problem Domain
33
Comparing the mixed FVE method to standard cellcentered finite differ
ences, the greater accuracy of the mixed FVE method is much more evident
when the diffusion coefficient is discontinuous. Below we present results for
a problem where A is a scalar, but discontinuous. The domain 0 is [1, l]2
and the diffusion coefficient has the following values:
0.05 for x > 0, y > 0
0.01 for x < 0, y > 0
10.0 for x < 0, y < 0
33.33 for x > 0, y < 0.
Boundary conditions specified the normal component of velocity to be 1
on the left and right edges of the domain and 0 (no flow) on the top and
bottom. Because an analytic solution for the problem was not available, the
mixed FVE solution from a 256x256 grid was used as the exact solution in
calculating errors. The velocity errors, measured in the same way as in the
last section, for the mixed FVE solution and the finite difference solution are
shown below. The finite difference discretization used harmonic averaging of
the diffusion coefficient at interfaces as is standard.
Mixed FVE
Grid Size urnfVe Vmfve ue
2x2 3.26E1 9.53E1
4x4 2.25E2 2.86E2
8x8 1.04E2 1.94E2
16x16 5.01E3 1.03E2
32x32 2.46E3 5.19E3
64x64 1.18E3 2.49E3
128x128 4.89E4 1.00E3
256x256 0 0
34
Finite Differences
Grid Size 4d *ld
2x2 8.28E1 9.53E1
4x4 4.77E1 5.03E1
8x8 2.70E1 3.07E1
16x16 1.47E1 1.69E1
32x32 7.94E2 9.02E2
64x64 4.23E2 4.77E2
128x128 2.24E2 2.51E2
256x256 1.18E2 1.35E2
The superiority of the mixed FVE approach is clear; it computes more
accurate velocities on a 16x16 grid than finite differences does on a 256x256
grid. This is not an artifact of using the 256x256 mixed solution as exact,
since the 128x128 and 256x256 finite difference solutions differ from each
other by an order of magnitude more than the corresponding mixed solutions
differ from each other. As in the comparison of the previous section, the two
methods exhibit comparable accuracy for the pressure (not shown). When
comparing solutions on different size grids, the mixed velocities do not show
secondorder convergence as they will in all other tests reported in this thesis.
The reason being that the true solution has a singularity at the origin, and
thus we cannot expect secondorder convergence.
1.3.3 An Alternative Multilevel Algorithm
A TwoLevel Approach
We pointed out in the last section that the multigrid algorithm presented
in Section 1.2 is robust with respect to jumps in the diffusion coefficient as
35
long as those jumps occur (if at all) at grid edges on the coarsest grid. For
a given diffusion coefficient, this limits the number of levels one can use in
the multigrid algorithm: we can coarsen the grid only to the level of the
jumps in the diffusion coefficient. If the diffusion coefficient has only a few
lines of discontinuity, this is acceptable. However, if the diffusion coefficient
has many lines of discontinuity, the coarsest grid we can use may be too fine;
there may be too many unknowns to solve the coarse grid version of the error
equation well enough, and we must clearly seek an alternative.
The problem that the previous multigrid algorithm has with jumps in the
diffusion coefficient that do not occur at interfaces of the coarsest grid can be
traced to interpolation. The piecewise constant interpolation scheme defined
in (1.15) is suitable so long as the diffusion coefficient is continuous within
the coarse grid cell. To illustrate the problem that arises when this is not
the case, consider the extreme situation shown in Figure 1.10.
In the two fine grid cells on the right, there is no flow and <(> = 0; these
are dead cells. In the two fine grid cells on the left, ^ is nonzero; these are
live cells. If we use piecewise constant interpolation from the coarse grid
cell, which will generally be live, then we can expect to interpolate some
nonzero value of ^ to the two dead cells, which is completely incorrect. The
solution to this problem, as discussed in [1], is to use information about the
diffusion coefficient in interpolating . This is the basis of the algorithm for
black box multigrid [8], a robust multigrid solver for the diffusion equation
36
Figure 1.10: Live and Dead Cells
37
with strongly discontinuous coefficients. It may be possible to incorporate
this operatorbased interpolation into the previous multigrid algorithm, but
we choose instead a different more convenient approach.
We will explain the approach as a twolevel multigrid algorithm, although
this is somewhat of a misnomer since we will have only one grid. Consider the
fine level problem as the mixed FVE discretization of the firstorder system
A1v + Vcf> =0,
V v = /.
We write the mixed FVE equations in matrix form as
M gradh W vh
divh 0 ) \ h
0
h
(1.22)
Here, M is the mass matrix that comes from the integration of Darcys
equation and gradh and divh are the grid h discrete operators corresponding
to the continuous operators grad, and div. Define the residuals according to,
V
0 M gradh \ / vh
fh ) V divh 0 ) \ h
(1.23)
where the variables with hats denote a current approximate solution to equa
tion (1.22). Define the errors by
eh = vh Yh^
e4> = h~ h
The error equation is then,
M gradh
divh 0
"V
0h
v
38
(1.24)
Now, rather than using a coarser grid with the mixed FVE discretization
to approximate the error equation, use the same grid with a standard cell
centered finite difference approximation. This is the coarse level in our
multigrid algorithm. The coarse grid version of the error equation can then
be written in matrix form as
M gradh
divh 0
(1.25)
The only difference between equations (1.24) and (1.25) is the mass matrix.
In (1.25) the mass matrix M is diagonal and computed from the diffusion
coefficient using harmonic averaging at discontinuities. The point is that, in
equation (1.25), we can eliminate the velocity variables. We have
vh = M 1 (gradh(f)h + r(j) . (1.26)
This allows us to rewrite equation (1.25) as
divhM 1gradhh = divhM Xr\. (127)
This is just the type of equation for which black box multigrid was developed.
Our approach is then to use black box multigrid to solve this equation for
(A, use equation (1.26) to get v^, and these approximations of the error to
correct our mixed FVE approximation.
This twolevel approach is similar to the work in [15], where black box
multigrid was used as a coarse level for a Lagrangian hydrodynamics appli
cation. In different vein, this twolevel approach is similar to preconditioning
39
the mass matrix M by its diagonal. This approach was suggested, in the
context of mixed finite element methods in [18].
TwoLevel Computational Results
We present here the results from a single numerical experiment to test the
robustness of the twolevel approach with respect to discontinuities in the
diffusion coefficient. The diffusion coefficients ax and ay were separately
and randomly assigned a value between .01 and 100 for each grid cell. The
diffusion coefficient is thus anisotropic and has jumps of several orders of
magnitude between cells. The twolevel algorithm described in the previous
section was used, with two alternating line relaxation sweeps on the mixed
FVE equations before using black box multigrid to solve the coarse problem
and one alternating line relaxation sweep after. One point to consider is: how
well do we need black box multigrid to solve the coarse problem? In [15]
only one cycle of black box multigrid was used for this purpose. For this
difficult test problem ( note: black box multigrid convergence factors were
approximately .6 ), we found that performing more than one cycle of black
box multigrid improved the overall convergence of the twolevel method. In
the results reported below, five cycles of black box multigrid are used to
approximate the coarse problem solution, although fewer, say 2 or 3, are
often enough to obtain similar twolevel convergence factors. The asymptotic
convergence factors for the twolevel method are given below.
40
Grid Size Convergence Factor
16 x 16 .43
32 x 32 .44
64 x 64 .46
We see that, although the factors are much less impressive than they are
for Poisson type applications, especially when work is considered, this two
level approach does exhibit convergence factors that seem to be bounded with
growing problem size. One point about these relatively poor convergence
factors (for comparison, the convergence factors of this twolevel algorithm for
constant coefficient problems is .1), these results show that, for this problem,
there is a significant difference between the mixed FVE discretization and
the standard finitedifference approach. This difference was illustrated in the
numerical tests of Section 1.3.2. When these two approaches are close, the
diagonal mass matrix M is an excellent approximation to the mixed FVE
mass matrix M, and the twolevel method converges quickly. Alternatively,
when the twolevel method converges relatively slowly, the two approaches
are significantly different, and we have seen that the mixed FVE approach
can be much more accurate.
The point of considering this twolevel method was to allow treatment
of problems where the coarsest grid for the multigrid algorithm presented in
Section 1.2 is still too fine and has too many unknowns to solve the mixed
FVE equations by itself. In practice, the multigrid algorithm presented in
Section 1.2 could be used down to the coarsest grid that was aligned with
the discontinuities in the diffusion coefficient, where the twolevel approach
41
described here could be applied.
1.3.4 Full Tensor Diffusion Coefficients
In this section we present results of two calculations with full tensor diffu
sion coefficients on uniform grids. The derivation of the mixed FVE Darcy
equation for this case is discussed in Section 3.1.
In the first problem, the diffusion coefficient A is constant throughout
the domain 0 = [1,1]2 and given by
_ / cos 9 sin 9 \ / 1 0 \ / cos 9 sin 9 \
^ sin# cos 9 J ^ 0 0.01 J ^ sin# cos 9 J
where 9 is the angle between the coordinate axes and the principal directions
of permeability. The boundary conditions and source term / were set so that
the exact solution to the problem (with integral over 0 = 0) is
(x, y) = cos [ttx) cos (27ry).
The mixed FVE discretization was solved on uniform grids of various sizes
and for various values of the angle 9. The errors in the resulting approxima
tions are shown below (measured in the same way as in previous sections).
o O II
Grid Size (j)mfve urnfVe V'mfve
8x8 0.4804 5.942E3 3.216E3
16x16 0.2501 2.140E3 1.091E3
32x32 0.1263 5.708E4 2.868E4
64x64 0.0633 1.449E4 7.252E5
42
9 = 15
Grid Size (j/n fve urnfVe V'mfve
8x8 1.836 1.925E1 8.872E2
16x16 0.5551 6.959E2 3.024E2
32x32 0.1803 1.949E2 8.268E3
64x64 0.0712 5.035E3 2.119E3
II CO o o
Grid Size (j/n fve urnfVe Vmfve ue
8x8 4.684 1.084E1 6.844E2
16x16 1.250 4.353E2 2.733E2
32x32 0.3357 1.391E2 8.679E3
64x64 0.1006 3.805E3 2.368E3
9 = 45
Grid Size (j/n fve urnfVe Vmfve ue
8x8 6.575 4.624E1 3.577E1
16x16 1.912 1.987E1 1.568E1
32x32 0.5254 6.324E2 5.076E2
64x64 0.1453 1.715E2 1.385E2
From the results we see, not surprisingly, that the accuracy of the dis
cretization degrades as the angle 9 grows. However, the results show second
order convergence in the velocities even when the anisotropy is not aligned
with the grid, so that any practical accuracy could likely be obtained by refin
ing the grid. When the anisotropy is not aligned with the grid, say, 9 = 45,
the convergence factors for the mixed FVE multigrid algorithm degrade (av
erage convergence factors of about .5 using alternating line relaxation), just
as the discretization accuracy does. Because the strongly coupled points are
at a 45 angle, convergence factors of .1 could likely be recovered by block
relaxation in this direction. However, this specialized treatment (block relax
43
ation in the dominant direction) would be difficult to implement in practice,
where the dominant direction would vary throughout the domain.
The second problem is similar to the test problem in [2]. Here the grid is
uniform, and A is a full tensor but no longer constant throughout the domain
0 = [0, l]2. The value of the diffusion coefficient is given by
0 < Â£ < .5 .5 < Â£ < 1
The boundary conditions and source term / were set so that the exact solu
tion is
0 < Â£ < .5 .5 < x < 1
where the value of C was chosen to make the integral of over 0 vanish.
Below are the mixed FVE discretization errors, measured in the same way
as before, for various grid sizes.
Grid Size 0'infve ufve V'mfve
4x4 9.212E2 4.555E3 5.748E3
8x8 4.622E2 1.110E3 1.540E3
16x16 2.313E2 2.741E4 3.931E4
32x32 1.157E2 6.825E5 9.871E5
64x64 5.784E2 1.755E5 2.435E5
In this problem, we clearly see firstorder convergence of the pressure and
secondorder convergence of the velocities. For all test problems with analytic
solutions, we have observed secondorder convergence of the velocities, and,
44
as will be shown in the next chapter, the mixed FVE discretization can be
extended to general quadrilaterals and still retain this desirable property.
45
2 General Quadrilateral Grids
2.1 The Mixed FVE Discretization
The focus now is extending the discretization technique of Section 1.1 to
irregular grids. We will develop the mixed FVE discretization for a logically
rectangular grid of irregular quadrilaterals, an example of which is shown
in Figure 2.1. Again, we consider the following partial differential equation
defined on some domain 0, now not necessarily rectangular, in 1Z2:
V AV(x) =/(x) x, G 0,
VA(x)(x)77, = 0 x G dft.
(2.1)
The diffusion coefficient A may be a function of position and is not necessarily
diagonal; however, it is assumed to be a pointwise symmetric, positive definite
matrix.
Again we separate out Darcys law,
A_1v(x) + V(x) = 0, (2.2)
from the conservation of mass equation,
V v(x) = /(x).
(2.3)
46
Figure 2.1: General Quadrilateral Grid
47
We will formulate a mixed FVE discretization for this system of equations
that generalizes the mixed FVE discretization presented in Section 1.1.
Important in developing the discretization for general quadrilaterals is the
mapping relating a general to a reference quadrilateral. Consider the general
quadrilateral P with vertices (V00, J/oo), (io, J/io), (oi, J/oi), and (Â£11,2/11), as
shown in Figure 2.2.
Let the reference quadrilateral P be the unit square. Then there is a
unique bilinear reference mapping of P onto P given by
x{x, y) = xQQ + (xio xQQ)x + (xoi xQQ)y + (xn xw x01 + x00)xy}
y(x, y) = yoo + {yio yoo)x + (y01 y00)y + (j/n yw y01 + y00)xy.
If P is convex, then this mapping has an inverse Since we restrict ourselves
to convex quadrilaterals, for each (x,y) Â£ P we have an associated point
(x, y) Â£ P given by the inverse reference map.
Shown in Figure 2.2 are several vectors that will be useful later in de
scribing the components of our discretization technique: for each (x,y) Â£ P,
we define the four vectors
X(x, y) is the image of the unit vector (1, 0) in P,
Y(x,y) is the image of the unit vector (0,1) in P,
Vx(x}y) is a unit vector orthogonal to Y(x}y)}
r/y(x,y) is a unit vector orthogonal to X(x,j/).
To extend the discretization of Section 1.1 to general quadrilaterals, we
need use analogs of the control volumes and finite element spaces used there.
For the latter, we use the lowest order RaviartThomas elements on the
quadrilateral elements (see [5],[19]), which are defined as follows. The char
48
Figure 2.2: General Quadrilateral
49
acteristic functions of the quadrilaterals provide a basis for the finite element
space for . The basis functions for v are best described in terms of associ
ating degrees of freedom with normal components on edges of quadrilaterals.
A typical basis function for the finite element space for v has support on two
adjacent quadrilaterals, a constant normal component on the edge shared
by the quadrilaterals, and a normal component of zero on other edges. The
magnitude of the basis function is such that the flux on the common edge is
one:
/ v yds = 1.
J edge
These conditions alone do not uniquely determine the basis function; the
following additional condition on the finite element space is needed. Within
any quadrilateral P, we assume that
v 77a, Y varies linearly with x, constant with y,
v 77^ 11X11 varies linearly with y, constant with x.
A typical basis function is represented in Figure 2.3.
We note that, just as in the uniform grid case, the basis functions have
continuous normal components across grid interfaces. So also will our com
puted solution. The nodes for the unknowns are shown in Figure 2.1. On
each grid cell, which we denote by P4j, we let denote the value of the
pressure on the cell, and uhi+ij, vhitJ, and u\j+i the values of the
flux on their respective edges. This is a slight departure from the method in
Section 1.1 where the unknowns uh and vh were defined to be the velocities:
applying the method of this section to a uniform grid results in the same
50
Figure 2.3: Basis Function
equations as in Section 1.1 except for a scaling of the uh and vh unknowns
by h.
We now need to choose the control volumes. By analogy to the uniform
grid case, the quadrilaterals used to describe the grid are the natural choice
for the control volumes for equation (2.3). This has the benefit that it leads
to a scheme with a local conservation property on these quadrilateral grid
cells. To proceed, we integrate equation (2.3), over each grid cell P4j:
/ V vdxdy = f.
J P i; j P i; j
Applying the divergence theorem, we get
9P.
yds =
/
51
The lefthand side of this equation is just the sum of the fluxes on edges of
Pjj, so the discretization of the mass conservation equation is
u
*+i ,j
u
hJ
+ V
i,j +1
V
hJ
If we assume that the function / is (approximately) piecewise constant, then
we can replace the integral on the right hand side by:
fax AREA(P8J).
(With more information about /, we could use a more accurate approxima
tion of the integral.)
The control volumes for Darcys equation generalize the control volumes
we used for the uniform grid case. Consider two adjacent grid cells, P8_i j and
Pjj. Then U8J consists of the image of (1/2,1) X (0,1) under the mapping
for Piij and the image of (0,1/2) X (0,1) under the mapping for P8J In
Figure 2.4, volume U8J is the shaded region.
We associate this volume with the vertical edge shared by P8_i j and P4j,
which the control volume straddles. We also have control volumes associated
with horizontal edges. For adjacent grid cells, P;j_i and P4j, volume V8J
consists of the image of (0,1) X (1/2,1) under the mapping for P;j_i and the
image of (0,1) X (0,1/2) under the mapping for P8J .
We now need to find a method to generalize the integrations of Section 1.1
that led to the discrete form of Darcys law. It is worthwhile to summarize
what we did there, but now interpreting it from a vector point of view. We
52
Figure 2.4: U Volume
53
took the xcomponent of Darcys law and integrated it over the U volumes,
which is equivalent to taking the dot product of the Darcy law equation
with the vector (1,0) and integrating the result over the U volumes. In
equation (1.7), we were able to integrate out the term, leaving integrals of
along lines interior to the P volumes. Because
on P volumes, in equation (1.8) we could evaluate the line integrals in terms
of the nodal values for
suitable vector held, with which to dot Darcys equation, that allows us
to integrate out derivatives of , leaving integrals that can be evaluated in
terms of the nodal values of
required because no single constant vector will produce the desired result.
Fortunately, simply scaling the vector held X yields the desired property. We
illustrate this process by obtaining the discrete Darcy equation associated
with the U volume shown in Figure 2.4. We proceed hrst by taking the dot
product of equation (2.2) with c;X(x, y) and integrating over the left half of
Ujj. Similarly, we take the dot product of equation (2.2) with crX(x, y) and
integrate over the right half of U8J. Here, q and cr are constants chosen to
eliminate terms involving values of ^ on the interface between the two halves
where
The details of the necessary integrations are in Section 3.2; we will present
here only the form of the equation which this integration gives rise to. Note
that we perform the same kind of integration for the V volumes, only we
54
take the dot product of Darcys law with a scaling of the vector Y. For the
U volume shown in Figure 2.4, we get a discrete Darcy equation relating the
pressure drop between the two cells to the fluxes on cell edges:
cluil,j + c2 ui,j + c3ui+l,j
\C4ViiJ + C5Uj_iJ__i + CgViJ + CfVi)j +1
^ilj) = O5
where \E\ is the length of the edge shared by the two adjacent grid cells.
The values of the coefficients ci,. ., C7 depend on the position of the vertices
defining the two grid cells and on the values of the diffusion coefficient within
the two cells. Note that the structure of connections here is the same as that
obtained for the uniform grid case with tensor permeability in Section 3.1.
The cross terms c4,. ., C7 will generally be nonzero even when the diffusion
coefficient A is diagonal.
2.2 Multilevel Solvers
2.2.1 FVE Based Multigrid
In this section, we discuss the changes in the multigrid algorithm presented
in Section 1.2 needed to deal with general quadrilateral grids. Because the
mixed FVE discretization for general quadrilaterals is simply an extension
of the discretization for uniform grids, only slight modifications to this algo
rithm must be made. Given a logically rectangular grid of irregular quadri
laterals, one can generate a finer grid of the same type by dividing each
quadrilateral in the given grid into four quadrilaterals. This division is done
55
by connecting midpoints of opposite sides. In this way, a sequence of finer
grids can be generated from a given coarse grid. Figure 2.5 shows a coarse
grid Q1 and its first refinement Q2. On each of the grids Q1,. ., QN, we apply
the mixed FVE discretization technique of the previous section. Then we can
use a multigrid algorithm to generate a solution to the FVE discretization of
the problem on the finest grid QN.
The basic components of the multigrid algorithm remain virtually un
changed. Now, because the structure of the connections in the discrete equa
tions remain unchanged, any of the relaxation techniques discussed in Sec
tion 1.2 can be used. As was the case for uniform grids, the interpolation
operators used in the multigrid algorithm are defined on the basis of the re
lationship between finite element spaces for the unknowns on different grids.
The finite element spaces are still nested, as we show in Section 3.3. The
relationship between the finite element spaces for is the same as in the
uniform grid case, so the interpolation operator for remains unchanged.
However, the relationship between the coarse and fine finite element spaces
for both u and v is changed somewhat, because these variables now represent
fluxes rather than velocities as before. Since the edge of a fine grid cell is half
the length of the edge of the corresponding coarse grid cell, the interpolation
operators for u and v of Section 1.2 must be scaled by 1/2.
Again, the restriction operators used in the multigrid algorithm are de
fined on the basis of the relationship between control volumes on different
56
Figure 2.5: O1 and Its First Refinement Q2
57
grids. The control volumes for the continuity equation are related in the
same way as in the uniform grid case (one coarse grid P volume consists of
four fine grid P volumes), and so the restriction operator for this equation
remains unchanged. The relationship between coarse and fine control vol
umes for Darcys equation can be different than it was in the uniform grid
case. A coarse grid U volume contains all of two of the fine grid U volumes
and fractions (no longer necessarily half) of four others. Figure 2.6 depicts
the relationship between coarse and fine control volumes. Strict adherence
to the FVE principles requires the calculation of the fraction of a fine grid
control volume included in the coarse grid control volume; however, we have
chosen to use the same restriction operator as in the uniform grid case. It is
certainly possible to calculate the fractions and use them in the definition of
the restriction operator, but, for the present, we have used this more easily
implemented definition (approximating the fractions by ). From our expe
rience with numerical testing of the algorithm, this approximation does not
appear to adversely affect the convergence of the multigrid algorithm. In
fact, because the convergence factors are generally quite small, it is not clear
that there is much to be gained by stricter adherence to the FVE principles
in defining the restriction operator.
Note that in the multigrid algorithm presented here, the process of devel
oping its components begins with the coarsest grid. Input to the algorithm
is a coarse grid that models the largescale irregularities of the reservoir ge
58
Figure 2.6: Coarse and Fine Darcy Volumes
59
ometry, and a specified level of refinement necessary to resolve the flow. The
algorithm then generates a solution to the mixed FVE discretization on the
finest grid. If a direct method is used to solve the equations on the coarsest
grid, then this places a practical limitation on the number of quadrilaterals
used to model the reservoir geometry.
2.2.2 A TwoLevel Approach
There are two practical problems with the multigrid algorithm discussed in
the previous section. The first is that, as in the uniform grid case, the jumps
in the diffusion coefficient must occur (if at all) at grid edges on the coarse
grid. The second is that the grid irregularity is constructed based on a given
coarse grid, which is then refined by bilinear coordinates to generate finer
grids: one cannot apply the multigrid algorithm to the equations on the
coarsest grid, so they must be solved some other way. In practice, both of
these problems represent limitations on the coarsest grid: it must be fine
enough to capture the reservoir geology and the jumps in the diffusion co
efficient. This limitation may result in the set of discrete equations on the
coarsest grid being too large to solve directly.
The twolevel approach described in Section 1.3.3 provides one possible
remedy to these two problems. The only component of the twolevel algo
rithm that needs additional discussion for the general quadrilateral case is
M, which is the mass matrix that multiplies the discrete fluxes in Darcys
equation. For uniform grids, M is defined by equation (1.9), and since its
60
diagonal, we invert it to eliminate the flux variables. In essence, we can use
a standard cellcentered finite difference discretization as the coarse level
for the fine level mixed FVE discretization. We would like do something
similar for the case of a general quadrilateral grids. However, one of our mo
tivations for looking at the mixed FVE discretization is that it can be applied
in a clear and direct way to general quadrilateral grids, where standard cell
centered finite differences cannot. Unfortunately, an effective cellcentered
finite difference scheme for discretization on general quadrilateral grids does
not currently exist. However, our demands on this scheme are not very strict
because the objective of the discretization on the coarse level is simply to
correct the solution from the fine level mixed FVE discretization, which
will serve as our final approximation. In other words, we would like to use
the coarse level discretization only to accelerate the relaxation process on
the fine level. For this purpose, our experience shows that equation (1.9)
serves as an effective way to define M even in the general quadrilateral case.
There are perhaps more sophisticated ways of defining M that we may ex
plore later, but we have found that this simple definition works well for many
applications. It is clear that, for very distorted grids, M defined this way
will provide a poor approximation to Af; however, we will see in the next
section that, for moderately distorted grids, the twolevel method works as
well as in the uniform grid case.
61
2.3 Computational Results
2.3.1 Multigrid Performance
This section contains results of numerical experiments designed to test the
robustness of the multigrid algorithm of Section 2.2.1. For all problems
considered here we used two prerelaxations and one postrelaxation in the
multigrid algorithm. For the first test problem, we used the mixed FVE based
multigrid algorithm with a coarsest grid of eight cells in each direction. We
began with a uniform 8x8 grid on 0 = [1,1] X [1,1] and then moved each
interior vertex in both the x and y direction separately by a random number
between .2h and .2h, where the original mesh size is h = .25. The resulting
grid is shown in Figure 2.7.
We then solved the mixed FVE discretization of Poissons equation, where
the diffusion coefficient is the identity, using the multigrid algorithm of Sec
tion 2.2 with varying levels of refinement. Ten or more cycles were performed
on the homogeneous problem and the average convergence factor was calcu
lated. The results are shown below.
Grid Size Convergence Factor
16 x 16 .12
32 x 32 .17
64 x 64 .22
Here we used distributive GaussSeidel relaxation and we see convergence
factors that are quite acceptable for this distorted grid in that they are rel
atively small, do not grow significantly with problem size, and are not much
62
Figure 2.7: Grid with up to 20% Distortion
63
different from the convergence factors for the uniform grid case. The algo
rithm allows us to solve this problem with about as much work as is needed
to solve the uniform grid problem.
For the second numerical experiment, we solved the same problem, only
now the 8x8 grid was distorted by moving each interior vertex in both the x
and y direction by a random number between .5h and .5h. The resulting
grid is shown in Figure 2.7.
This grid is likely to be much more distorted than one would encounter in
an application, and, in fact, some of the quadrilaterals are no longer convex.
The results for the mixed FVE based multigrid algorithm are shown below.
Grid Size Convergence Factor
16 x 16 .03
32 x 32 .05
64 x 64 .08
Because of the significant distortion of the coarsest grid, we used alternat
ing line relaxation for this test example: analogous to the anisotropic uniform
grid case, the structure of the resulting discrete mixed FVE equations con
sists of cells that are strongly coupled to some neighbors and weakly coupled
to others. We need alternating line relaxation to get the error smoothing
that is needed for acceptable multigrid convergence factors. With this more
expensive algorithm, compared to the cell by cell distributive GaussSeidel
relaxation, we obtain quite good multigrid convergence factors for this very
distorted grid.
The third numerical test combined moderately distorted grids with a
64
Figure 2.8: Grid with up to 50% Distortion
65
diagonal but anisotropic and discontinuous diffusion coefficient. The coarsest
grid here is the same as that used in the first experiment, an 8x8 grid with up
to 20% distortion. On each cell in this coarsest grid the diffusion coefficients
ax and ay were separately set to random values between .01 and 100. The
results for the mixed FVE based multigrid algorithm are shown below.
Grid Size Convergence Factor
16 x 16 .04
32 x 32 .06
64 x 64 .08
Again we used alternating line relaxation, now because of the anisotropic
diffusion coefficient. Again we see quite good multigrid convergence factors.
This numerical test is close to the types of problems that might actually
occur in an application and the performance of the multigrid algorithm is
well within an acceptable range.
The two restrictions on the coarsest grid for our mixed FVE based multi
grid method that were discussed in Section 2.2.2, namely, capturing the ir
regularity and discontinuities of the diffusion coefficient, have not been vio
lated in the numerical experiments thus far. However, even without violating
these limitations, one can pose problems where the mixed FVE based multi
grid method fails to perform as well as in the above tests. For example,
consider the coarsest grid used in second numerical experiment: 8x8 with
up to 50% distortion. On each cell in this grid, the diffusion coefficients ax
and ay were separately set to random values between 1/a and a. The results
for the mixed FVE based multigrid algorithm using four grids (finest grid is
66
64x64) are shown below.
a Convergence Factor
1 .08
10 .13
100 .15
1000 .20
These last results may be of little practical interest as the grid is much
more distorted than one would likely use in an application, but they do
demonstrate the robustness of the multigrid algorithm. For all numeri
cal tests encountered thus far, the mixed FVE based multigrid algorithm
with alternating line relaxation has yielded acceptable (rarely worse than .3)
convergence factors if the lines of discontinuity in the diffusion coefficient
are captured on the coarsest grid. With full tensor coefficients, when the
anisotropy is not aligned with the coordinates, the convergence factors can
be significantly larger (at times, as large as .6), and we have no rigorous
theoretical results, so we can make no general claims. However, it appears
that the method has potential as a fast solver in many important practical
applications.
2.3.2 TwoLevel Performance
This section contains results of numerical tests designed to test the robustness
of the twolevel approach of Section 2.2.2. This is an approach that one could
use alone or as the coarsest grid solver in the mixed FVE based multigrid. For
the results of this section, we use alternating line relaxation, with two sweeps
67
performed before solving the coarse problem and one sweep performed after.
Five cycles of black box multigrid were used to solve the coarse problem,
although often two or three are sufficient. In the first numerical experiment,
the mixed FVE discretization of Poissons equation on 0 = [1,1] X [1,1]
was solved on the moderately distorted grid obtained from a uniform grid
by the up to 20% distortion described in the previous section. Note that
the distortion here takes place on the finest grid, so the mixed FVE based
multigrid algorithm cannot be used. Average convergence factors for the
twolevel approach on different size grids are shown below.
Grid Size Convergence Factor
16 x 16 .07
32 x 32 .08
64 x 64 .08
The convergence factors are remarkably good, given the quality of the
approximation used on the coarse level. As discussed in Section 2.2.2, we
basically assume the grid is uniform in constructing the mass matrix M for
the discrete Darcy equations on the coarse level. This appears to work well
for the moderately distorted grids like those in this numerical experiment
and, quite probably, the grids one would use in practical applications. When
the grid is very distorted, like the up to 50% distorted grid of the previous
section, the twolevel algorithm can fail to converge and may even diverge.
The reason is that the very poor approximation of the mass matrix results
in a correction from the coarse level that has little, if anything, to do with
the fine level error. It is possible that this could be remedied by a more
68
sophisticated choice for M, but this has yet to be investigated.
In the next numerical experiment, we use the same grids as in the previous
test and solve the mixed FVE discretization to the diffusion equation with
diagonal diffusion coefficient, where on each cell in this grid the diffusion
coefficients ax and ay were separately set to random values between .01 and
100. The results, in terms of average convergence factors, are shown below.
Grid Size Convergence Factor
16 x 16 .43
32 x 32 .40
64 x 64 .38
While the convergence factors are not as small as for Poissonlike prob
lems, they are likely acceptable in most practical applications, especially if
one is using the twolevel approach only on the coarsest grid of the mixed
FVE based multigrid algorithm. There the amount of work on finer grids
in the mixed FVE based multigrid algorithm will typically be much larger
than the work of the twolevel algorithm on the coarsest grid, even if sev
eral cycles of the twolevel algorithm are required. We should note that the
convergence factors for this problem are roughly equal to those obtained for
the similar uniform grid problem: our belief is that the poor convergence
factors are again due to the difference between the mixed FVE methods and
the finite difference methods treatment of discontinuities in A, not related
to the moderate irregularity of the grid.
69
2.3.3 Accuracy of the Mixed FVE Discretization
In this section, we present results of numerical experiments to investigate the
accuracy of the mixed FVE discretization on general quadrilateral grids.
First, we will present the result of a simple numerical test designed to
illustrate the effect that the shape of the grid cells has on discretization errors
and multigrid convergence rates. We consider Poissons equation posed on
Q = [l,l]2:
v + V = 0,
v V = /,
with boundary condition v rj = g. For this test, / and g were determined
so that the exact solution is,
(x, y) = cos(7rx) cos(27ry) + constant.
The solution from the mixed FVE discretization was normalized to corre
spond the exact solution with constant= 0. The coarsest grid we use in our
tests is shown in Figure 2.9.
To investigate the effect of the angle f3 on the accuracy of the mixed FVE
discretization and on multigrid convergence factors, we use varying /3. We
measure the error in the same way as in Section 1.3, using the C2 norm of the
error and C2 norm of the u and v errors. In the tables below, we compile
our observations on these errors and the average convergence factors for the
mixed FVE based multigrid algorithm. For f3 = 80 and 70, distributive
GaussSeidel relaxation proved effective in terms of convergence factors. For
70
Figure 2.9: Coarsest Grid
71
f3 = 60 and 50, the stronger alternating line relaxation scheme was needed
to get acceptable convergence factors. This is again due to the distortion of
the grids that yield couplings of different strengths to neighboring cells.
f3 = 80 distributive GaussSeidel
Fine Grid ue Ve Convergence Factor
8x8 4.841E1 1.575E1 7.595E2 .054
16x16 2.514E1 5.090E2 2.538E2 .073
32x32 1.267E1 1.347E2 6.752E3 .083
64x64 6.354E2 3.414E3 1.713E3 .086
f3 = 70 distributive GaussSeidel
Fine Grid ue Ve Convergence Factor
8x8 4.881E1 2.327E1 1.106E1 .067
16x16 2.541E1 6.541E2 3.330E2 .090
32x32 1.285E1 1.689E2 8.809E3 .107
64x64 6.516E2 4.256E3 2.235E3 .143
f3 = 60 alternating line relaxation
Fine Grid Ue Ve Convergence Factor
8x8 5.009E1 3.506E1 2.140E1 .019
16x16 2.619E1 9.111E2 5.723E2 .031
32x32 1.335E1 2.305E2 1.4832 .035
64x64 6.967E2 5.779E3 3.745E3 .054
f3 = 50 alternating line relaxation
Fine Grid Ue Ve Convergence Factor
8x8 5.318E1 6.060E1 4.169E1 .022
16x16 2.760E1 1.737E1 1.199E1 .038
32x32 1.406E1 4.493E2 3.189E2 .057
64x64 7.382E2 1.135E2 8.142E3 .097
These results, show good convergence in the accuracy of the discretization
as the grids are refined. Each refinement reduces the norm of the error in
by a factor of 2 and the error in the fluxes by a factor of 4. Although
72
we do not have rigorous proofs of the order of convergence of the mixed
FVE discretization, these results, and all other results from experiments with
analytic solutions, appear to indicate that, at least in the norms we used,
the convergence is first order in the pressure and second order in the fluxes.
We also see that the accuracy of the discretization decreases as f3 decreases,
which is not surprising since most finite element approximations degrade
as the grid becomes less uniform. However, one can still likely achieve a
solution of any practical accuracy by refining the grid yet further. Note also
that the multigrid convergence factors degrade slightly as the grid becomes
less uniform, but the factors are still certainly acceptable.
The next problem is similar to the previous one in that the same type
of grids are used, and the boundary condition and right hand side were
determined so that it has the same exact solution; however, this problem has
full tensor permeability. The diffusion coefficient A is given by
_ / cos 9 sin 9 \ / 1 0 \ / cos 9 sin 9 \
^ sin# cos 9 J ^ 0 0.01 J ^ sin# cos 9 J
where 0, the angle between the coordinate axes and the principal directions
of permeability, is 45. The errors in the mixed FVE approximations for
various size grids and varying angle of distortion f3 are given below:
'CO II 00 o o
Fine Grid Size ue ve
16x16 2.261 2.722E1 1.998E1
32x32 0.6356 9.497E2 7.194E2
64x64 0.1743 2.714E2 2.078E2
128x128 0.0520 7.042E3 5.442E3
73
'CO II cs o o
Fine Grid Size ue Ve
16x16 3.092 4.527E1 3.021E1
32x32 0.9468 1.963E1 1.398E1
64x64 0.2698 6.375E2 4.568E2
128x128 0.0788 1.747E2 1.267E2
As in the uniform grid case, the accuracy degrades due to the anisotropy
not being aligned with the coordinates, but the results still indicate second
order convergence in the velocities. Note that for f3 = 60 and 6 = 45, a
case with severe distortion and nongrid aligned anisotropy, the asymptotic
regime is not reached until the grid is quite fine.
The final test problem involves a nonuniform grid and full tensor, discon
tinuous permeability. The domain for the problem, which consists of three
regions, is depicted in Figure 2.10.
The coefficient A is defined in each region as follows:
In region I In region II In region III
/ 2 1
l 
4 4
U 4
/2 
1 i
2 2
VI 1
The boundary conditions and right hand side were defined so that the exact
solution is given by,
In region II In region III
In region I
6 = x2 + C
= Ty+c
= b2 + c
Again the value of C is chosen to correspond to the solution whose integral
over the domain vanishes. Note that, and the normal component of AV
are continuous across the interfaces between regions. The coarsest grid (a
74
Figure 2.10: Problem Domain
75
2x2 grid), was formed by splitting, along the xaxis, region I into two cells,
and all finer grids were generated by bilinear refinement of this coarsest grid.
The errors in the mixed FVE approximation, calculated in the same way as
previously, are presented in the table below.
Grid Size ue Ve
4x4 0.2524 7.065E3 8.269E3
8x8 0.1302 1.747E3 2.197E3
16x16 0.0658 4.567E4 5.859E4
32x32 0.0333 1.213E4 1.582E4
64x64 0.0174 3.234E5 4.294E5
Here we clearly see firstorder convergence in the pressure, and again we
see secondorder convergence in the velocities, which has been the case for
all experiments, ranging from Poissons equation on a uniform grid to this
significantly more complicated case, with analytic solutions.
76
3 Technical Derivations
This chapter contains all the derivations needed to assemble the mixed
FVE equations. Many of the results from this section have been briefly
presented, or at least alluded to, in previous sections.
3.1 Uniform Grids and Tensor Permeability
In Section 1.1, we presented the FVE discretization assuming that the coef
ficient A in equation (1.1) is a diagonal tensor. Here we consider the case
where A is a full tensor, but we assume that it is a pointwise symmetric,
positive definite matrix. In this section, we assume the grid is uniform (see
the following section for the general quadrilateral case). We write
A
a b
b c
where a, 6, and c may now be functions of position. Darcys law, equa
tion (1.2), takes the form
v = AV(>,
77
or
A 1v + V = 0.
Now, adding the conservation of mass equation and noting that v = (it, n)4,
we have the following firstorder system:
(detA)1 (cit bv) \d(f)/dx =0 (u equation)
< (detA)1 ( bu + av) \d
du/dx \dv/dy = f (p equation).
The main components of the mixed FVE discretization process remain the
same: the choice of volumes and finite element spaces are the same as pre
sented in Section 1.1. The change in the form of the u equation, comparing
equation (3.1) to equation (1.4), is reflected when we integrate the u equa
tion over the U volumes. In this section, we will assume that A is constant
over the integration domain, otherwise the technique of Section 3.3 can be
applied. As an example, let U8J be the volume in hi that is centered at the
interior uh node (i,j) as shown in Figure 3.1.
We then have
/ ((detA) 1(cu bv) + d(f>/dx^ dxdy = 0. (32)
Integration of the terms involving u and d(p/dx is the same as in Section 1.1.
The additional term involving v is handled by splitting the integral into two
parts, one for each half of U8J. Then, because our finite element approxima
tion vh is linear in y and constant in x on each half, by direct calculation of
78
Figure 3.1: U Volume
79
the integral over the left half, we have
T
Combining this result with the analogous one for the other half and adding
the results from Section 1.1 for the other terms in equation (3.2), we get the
discrete equation
V(c(det A)"1) (tw +
T (fc(det A)1) ("tw+i + ?W + <+. + C)
+ h (jt, ~ tiu) =
We get a similar equation when we integrate over V volumes. This change
in the form of the discrete Darcy law is the only change in the discretization
process from Section 1.1.
3.2 Derivation of FVE Stencils for General
Quadrilaterals
3.2.1 Scalar Diffusion Coefficient
Here we develop the details of the integrals encountered in deriving the dis
crete Darcy equation for general quadrilateral grids. We first consider the
case where the diffusion coefficient is a scalar. The case where the diffusion
coefficient is a tensor is covered in Section 3.2.2. Shown in Figure 3.2 by
dotted lines is a typical U volume, which we will use to discretize Darcys
equation.
80
Figure 3.2: U Volume
81
P; and Pr will denote the two P volumes who share the edge with which
the U volume is associated. The continuous Darcy equation is
a 1v(z,i/) + V0(x,y) = O.
(3.3)
As mentioned in Section 2.1, we take the dot product of equation (3.3) with
QX(x,y) and integrate over the left half of the U volume. Similarly, we
take the dot product of equation (3.3) with crX(x,y) and integrate over the
right half of the U volume. We then add the two integrals to get the discrete
Darcy equation corresponding to this U volume. Letting U; and Ur denote
the two halves of the U volume we have the following:
fu, (a Ma V) + y)) ciX(x} y)dxdy+
fur (a_1v(x, y) + V(x, y)) crX(x, y)dxdy =
(3.4)
Q/\ja 1v(x} y) X(z, y)dxdy + q/U; V^(z, y) X(x, y)dxdy+
Crhj^vix, y) X(x, y)dxdy + crfUrV(f>(x, y) X(x, y)dxdy.
Thus, we have two types of integrals to evaluate, ones involving v and ones
involving Xcf>.
First consider the integral involving v. For the moment, we ignore the
scaling parameters q and cr. We evaluate this integral for the finite element
basis functions for v, which are the formulas needed to assemble the stencil
for equation (3.4). Consider the evaluation of
where w is the basis function with flux equal to 1 on the edge straddled by
the U volume and 0 on all other grid edges. It is important to recall the
82
properties used to define our finite element space: within each grid cell,
w(x, y) yx(x, y)  Y(x, y)  varies linearly with x, constant with j/,
w(x, y) r]y(x, y) X(x, y)  varies linearly with y, constant with x.
(3.5)
The first condition implies, that within the volume P/,
w(x,y) rjx(x,y)\\Y(x,y)\\ = x,
where x refers to the preimage of x in the reference space. Similarly, y will
denote the preimage of y in the reference space. The second condition, along
with the fact that fluxes are zero on the top and bottom edges of P/, implies
that, within the volume P/, we have
w(x,y) T]y(x} 2/)X(m, y)\\ = 0.
Thus, w is perpendicular to rjy and, hence, parallel to X. We thus have,
w(x,y) X(x,y) = w(x, y) X(x, y)
and
w(x,y) rjx{x,y) = w(x,y) sin(0(x,y)),
where 9 is the angle between X and Y as shown in Figure 3.2. Solving this
last equation for the norm, we get
w(z, y) = (w(x, y) yx(x, y j) / sin(0(a;, y j).
(3.6)
83
Putting these results together, we have,
fu ,a_1w(a;, y) X(x, y)dxdy
= /u,a_1M>2/)lllX(z,2/)
(3.7)
= /uffi 1 ((w(a;, y) dx{x, y)) / sin(0(x, y))) X(x, y)\\dxdy
= /o/i_1X(x, y)\\2dxdy.
The last step is the result of the change of variables from the actual space to
the reference space. The Jacobian for the map is
J = 11X111Y sin 0.
This final integral can now be evaluated analytically in the reference space.
We had defined the vector held X(x,j/) by taking images of the unit vector
(1,0) under the bilinear mapping
x(x, y) = xQQ + (xio xQQ)x + (xoi xQQ)y + (xn xw x01 + x00)xy}
y(x, y) = Voo + {yio yoo)x + (y01 y00)y + (yn yw y01 + y00)xy.
We can express the vector held X as follows:
X = (dx/dx7dy/dx)
= ((x10 Xoo) + (xn x10 Xoi + oo)y, (3.8)
(yio yoo) + (yu yio yoi + 2/00)2/)
Thus, the last integrand in equation (3.7) is a polynomial in x and y whose
coefficients depend on the coordinates of the vertices of P; and this integral
can easily be evaluated analytically.
Next consider the evaluation of
/ a~1w(x, y) X(x, y)dxdy,
84
where w is the basis function with flux equal to 1 on the top edge of the left
grid cell and 0 on all other grid edges. The second condition in (3.5) implies
that, within the volume P/, we have
w(x,y) T]y(x} 2/)X(m, y)\\ = y.
The first condition in equation (3.5), along with the fact that the fluxes are
zero on the left and right edges of P/, implies that, within the volume P/, we
have
w(x,y) f]x(x,y)\\Y(x,y)\\ = 0.
Thus, w is perpendicular to rjx and, hence, parallel to Y. We thus have,
w(ab2/)/llw(ab2/)ll = Y{x,y)/\\Y{x,y)\\
and
w(x,y) rjy(x,y) = w(x,y) sin(0(x,y)).
Solving this last equation for the norm, we get
\\w(X, j/) = (w(x, y) f]y(x, y)) / sin(0(a;, y)). (3.9)
Putting these results together, we have,
/U;alw(x, y) X(x, y)dxdy
= llw(^7 2/)II (w(x,j/)/w(x,j/)) X(x,y)dxdy
= IU;a_1 ((w(5 V) Vy{x, y)) / sin(6'(x, y))) (3.10)
(Y(x,y)/\\Y(x,y)\\)X(x,y)dxdy
= /o/l_1yX Ydxdy
85
Again, the last step is a result of the change of variables. This integral can
now be evaluated analytically in the reference space. The vector held Y can
be written in terms of reference space coordinates as follows:
Y = {dxjdy, dy/dy)
= ((xoi Xoo) + (xn x10 Xoi + x00)x, (311)
(s/oi yoo) + (yn 2/10 2/01 + yoo)x).
Thus, the last integrand in equation (3.10) is a polynomial in x and y whose
coefficients depend on the coordinates of the vertices of P/, and this integral
can easily be evaluated analytically. Integrals for other basis functions are
analogous to one of the two that we presented because they can be evaluated
in the same way with a simple change of variables. The basis functions that
will have nonzero contribution to the integral
j a 1w(x, y) X(x, y)dxdy
will be those that are associated with one of the edges of either P; or Pr.
There are seven such edges, so the discrete version of Darcys law on the
volume U will involve seven nodal values for fluxes.
Now consider the terms involving V, for example,
/ V(:r, y) X(x, y)dxdy.
J U,
86
We again express X in of terms reference space variables to get
fu, y) X(a;, y)dxdy
= ful (d(f>/dx, d
= + l^%dxdV (3.12)
= f xjydcf)/dxdxdy
= /o/iff J dxdy.
Here, J is again the Jacobian of the mapping, which can be expressed as
j dx dy dx dy
dx dy dy dx
Note that, expressed in terms of reference space variables, we have
J(x}y) = a + bx + cy7
where a,b, and c are constants on the quadrilateral P; and depend on the
location of its vertices. With this and the last result in equation (3.12), we
have
fu, V(Â£, y) X(x, y)dxdy
= dxdy
= foil HdS:dS:
= flJ(i,y)(i,y)dy SlJ{\,y){\,y)dy bflf\dxdy.
We have now accomplished the goal of integrating out the derivatives of .
The second integral is a line integral interior to the volume P; and can be
87
written in terms of nodal values of ; also, the third integral is a volume
integral within the volume P; and can be written in terms of nodal values
of . Because the approximation of is constant on Pi (constant= i,j), for
the third integral we have
However, the first integral is a line integral along the edge shared by P; and
Pr, and because we approximate with a piecewise constant function, we
seek to remove the first integral. In terms of the Jacobian of the map for the
right grid cell in Figure 3.2, similar integrations over Ur produces
JurV(K, y) X(a;, y)dxdy =
foJ(b $)(%, y)dy foJ(0, y)(0, y)dy bflf$ (f)dxdy.
In this expression, the first and third integral can be expressed in terms of
nodal values of ; the second integral involves integration along the edge
shared by P; and Pr. It is at this point that we see the utility of the scaling
parameters q and cr. We choose them so that when we add the results from
the two halves, the terms involving on the interface between the two grid
cells cancel. In fact, choosing q = \E\/J(1, ) and cr = \E\/J(0, ) gives us
the result
ci V(f)(x,y)X(x,y)dxdy+cr 'V(x,y)X(x,y)dxdy = \E\{cj)l+1^cj)hl).
J TJ i J TJ r
Here, \E\ is the length of the edge shared by P; and Pr.
88
3.2.2 Tensor Diffusion Coefficient
For a full tensor coefficient, much of the results from the previous section
are unchanged, including, for example, all of the results for the part of the
integral involving . We must only reconsider the evaluation of
/ A1 w(A, y) X(x, y)dxdy
J u
when w is one of the two basis functions discussed in the previous section.
Again we split the volume into a left and right half and integrate over each
half separately. We will use the same notation as in the previous section, so
that Figure 3.2 still represents the U volume that we integrate over.
We first consider the case where w is the basis function that has flux
equal to 1 on the edge straddled by the U volume and 0 on all other grid
edges. We proceed as follows:
fu, Alw(x, y) X(x, y)dxdy
= fu,w(x>y) ATX(x, y)dxdy
= /uJlw(A2/)ll (w(m, 2/)/w(m, y) ) A ~TX(x,y)dxdy
= fu, ((W(A v) Vx(x, y)) / sin(0(a;, y))) (3.13)
(X(a;, y)/\\X(x, y) ) ATX(x, y)dxdy
= foflxX(x,y) A~TX(x,y)dxdy
= foflxA^Xix, y) X(x, y)dxdy.
The hrst step follows from the adjoint property; the third step uses equa
89
tion (3.6) and the fact that w and X are parallel; the fourth step is the
result of the change of variables to the reference space; and the final step
again follows from the adjoint property. We can then evaluate the final inte
gral in equation (3.13) in the reference space using equation (3.8) to express
X in terms of the reference variables. The integrand, as in the previous sec
tion, is just a polynomial in x and y, so this integral can easily be evaluated
analytically.
We now consider the case where w is the basis function that has flux
equal to 1 on the top edge of the left grid cell and 0 on all other grid edges.
We proceed as follows:
fu, Alw(x, y) X(x, y)dxdy
= fu,w(xiy) ATX(x, y)dxdy
= /uJIw(a?/) (w(ab2/)/llw(ab2/)ll) A~TX(x,y)dxdy
= fu, ((w(> V) ' l/)) / sin(0(a;, y))) (3.14)
(Y(x, j/)/Y(x, j/)) ATX(x,y)dxdy
= Jo/lj/Y(x,y) A~TX(x, y)dxdy
= JlJ'ixA1Y(x, y) X(x, y)dxdy.
The steps are the same as in the previous case except for step three, which
uses equation (3.9) and the fact that w and Y are parallel. Again, the
final integrand is a polynomial in the reference variables and can easily be
evaluated analytically.
90
Note that for given A, one could consider choosing the quadrilaterals such
that A1 Y X, making this integral vanish. This is unlikely to be useful in
complex problems, but if A is constant for example, it could be of interest.
Integrals for other basis functions can be evaluated in a manner similar
to one of the two we presented here by a simple change of variables.
3.3 Using Fine Grid Integrals to Generate
Coarse Grid Integrals
In this section, we will show how to generate coarse grid integrals from fine
grid integrals, which is important when the permeability coefficient A is
piecewise constant on the fine grid cells but not on the coarse grid cells.
We begin by examining the relationship between fine and coarse grid
finite element spaces. As mentioned in Section 2.2, the finite element spaces
are nested; we will now make this statement more precise. In Figure 3.3, we
depict a coarse grid volume Pc and the four fine grid volumes Py i = 1, 2, 3, 4,
that are generated by bilinear refinement.
Letting Wc denote the coarse grid velocity basis function with nonzero
flux on the north edge of Pc and similarly letting w4 denote the fine grid
velocity basis function with nonzero flux on the north edge on Py we will
prove the following result: within Pc, we have
Wc = ^ (wj + w2) + ^ (w3 + w4). (3.15)
We prove this result by verifying that this linear combination of fine grid
91
Figure 3.3: Coarse and Fine Volumes
92
basis functions has all of the properties that the coarse grid basis function is
to have, namely:
1) The fluxes are zero on the east, south, and west edges of Pc.
2) The flux is one on the north edge of Pc
3) The basis function is parallel to the vector held YÂ£/),
i.e., the basis function is aligned with the bilinear coordinates of Pc
that correspond to vertical lines in the reference space.
4) The quantity W t^HXH varies linearly with y, and is constant
with x, where all quantities are coarse grid quantities.
(3.16)
The first property follows from the fact that each fine grid basis function,
wi = 1, 2, 3, 4, has zero flux on the east, south, and west edges of Pc Thus
the linear combination also has this property.
We now verify the second property in (3.16):
fpn (Â§ (wi + w2) + (w3 + w4)) yds
= 2 Jp" (wi y) ds T ^fpn (w2 y) ds
= (1) + (1) = 1
The first step follows from the fact that both w3 and w4 have zero flux on
the north edge of Pc, and the second step follows from the fact that both
w4 and w2 have flux = 1 on the north edges of P4 and P2, respectively.
The third property in (3.16) follows from the agreement of coarse and
fine grid bilinear coordinates.
We are left with the fourth property in (3.16), to verify this, we will need
the following result:
If (x, y) G Pi, then Xc(z, y) = 2X8(x, y). (317)
93

PAGE 1
#2)314)5%6))6))4)+%3 1%0#//0#)/%*#%4%1163 5)6%/)74*%0%7)3# 8)8 749.:, 749.:. # 1;7 "/3" 3* # 9..< %01%"22"13"40"%01"$%"0224
PAGE 2
3* 8)8 6 71/! # === 5> 8?0@' 10 8 %01%"22"13"40"%01"$%"0224
PAGE 3
88)&*3#' #15)#/ 15* = *71/! # 70#/ #! ""6( A ""/! B"0" "$" 7"" """ $ "*CA %01%"22"13"40"%01"$%"0224
PAGE 4
A """ C 7 71/! %01%"22"13"40"%01"$%"0224
PAGE 5
2 )8%24444444444444444444444444444444444444444444444444444444444444444444444444444444444444444444444) 9915)3$9 9D#15) #9E 9D/0D< 9D9/13**D< 9DD3#3/FFFFEG 9DE##""#E< 9DH13/HD /890%8%244444444444444444444444444444444444444444444444444444444444444444'* D915)3$H, DD"7"<< DD915) << DDD#(6"#,G DE/0,D DE9*,D %01%"22"13"40"%01"$%"0224
PAGE 6
DED(6"*,I DEE#15)3$IG 1"1:6244444444444444444444444444444444444444444444444444444444444444444444444444444E9;*II ED3"15)7;@:G ED973/:G EDD3/:. EE1;//;FFF.9 '0%00;<444444444444444444444444444444444444444444444444444444444444444445$3"4444444444444444444444444444444444444444444444444444444444444444444444444444444444444444444444444)!! %01%"22"13"40"%01"$%"0224
PAGE 7
#/J4%?6)3;))47 "( 1!*7"/!" !B! */!" !*0 !?" "!!*8 0*7"7?"" ((!B! !!3*0)0 746"A 1! (31!( (! " %01%"22"13"40"%01"$%"0224
PAGE 8
)8%2 )4)"=%>0:2? 17 "A $K" +" $" ? ( (5(#&'5>&'LM&' 2N &'"O&'P .&' &' "A ( 9 %01%"22"13"40"%01"$%"0224
PAGE 9
$" (A """ "LQ#5O&9D' &99' 5P"LM&9E' "&9D'3C &9E'"" (R99SR9DSR9:S "$#A + 7E9 $#6 L" (&99' / %01%"22"13"40"%01"$%"0224
PAGE 10
( (9T LG&' !" ( LG&"'&9H' ( &' &'LUV&'W?< &'L<&' #$#%& "&9H' $"&15)' "RISR9ESR9HSB15) $""A 1"15)$A & $ & "? "&9H' )"199? "$&9H' 75 "5* ""4"*A !#" "#"5$ (& E %01%"22"13"40"%01"$%"0224
PAGE 11
0@AB 6@ACB 0@CAB 6@AB D60 60 ,E>60 30)4)F8% %01%"22"13"40"%01"$%"0224
PAGE 12
415)$1 0"( R9,SR9.S" & *N & *N* )& *N 4 & & $ !199 &99'" & & 4 O# & 19D 4 $* $? "" 7D9( "(" R9IS < %01%"22"13"40"%01"$%"0224
PAGE 13
19DK 1 ?$7! &9H'"N "3C#(B" & & ) 19E ?" *)+ &9<' /1 &9<'(B M&(9' &>(9O' C @)4*B $ X %01%"22"13"40"%01"$%"0224
PAGE 14
19EK(B5 I %01%"22"13"40"%01"$%"0224
PAGE 15
&(BRB"" 4 & &&Y9' L((F9&G&E8(TPFZ['P #"B 4" "(9&M' " 7EE 4 &9<'" @ &M )&)/+)) &9I' & G L & +)) G & @H++ @ 9 4 : B ( %01%"22"13"40"%01"$%"0224
PAGE 16
*"&9<' IJKLM@NOB@ CJC&&BCJ4PQ@ $& CJR+ ) &BBS (9 C 0& CJOOC&BCKJ $/$) GT+ 4 (3CK $ T 1$/$/) LG(&9(.' 15) ""A "%"&9.' &9,'7A 9E 15)"A 3C ( "&9H'"5" """ & & %01%"22"13"40"%01"$%"0224
PAGE 17
$A 9H"" & ? "*[B719H "1 *)2 5 / 3$ "K M*5P 43/)# 56 /457 T 458 ( 45 >9 9 :9 B" & & 4 2457 ( 4 *\ Q 45 ;<9 > &=*$/$9*$$*$!'$)/ "&9H'B ($*$!$ C)+G>U+ +"M "M & ) # !K" [M "1 )! %01%"22"13"40"%01"$%"0224
PAGE 18
C 19HK*5 99 &TB' ] %01%"22"13"40"%01"$%"0224
PAGE 19
!M "; $"3C ; "5 ; 5 "A $"" ; 1*""&' "&!'" ""15)$ & '" A "" 1" $A $# 3A &99' UW&'LW&'; R9GS15)$*C 3 )/ %01%"22"13"40"%01"$%"0224
PAGE 20
)4/>2%03%3" ?""" ; ^CR,S"^CRES# R9HS&""' B15)$ /A .&& L & &99G' ?!&99G'$A $ & &RHS'!" ;(7 & &99G'%" "K >&&/.&& (" 9E %01%"22"13"40"%01"$%"0224
PAGE 21
"( .&6&>& &999' K 6& L &2& !(A ? "&999'K .&& &&>& &99D' +! .& $ $&99G'IM_ >& %"&99D' & & & `( &* ?&& &99E' +IO & ? 9H %01%"22"13"40"%01"$%"0224
PAGE 22
(" K 9'0 D'1" E'7"" H'/ <'0 "" "K "%"" $H "% 19< .&.&. H' &' "(A "" "B 4"" 9< %01%"22"13"40"%01"$%"0224
PAGE 23
> > ; ;D a V /; 06L 3"6L 1" 19
PAGE 24
15)$ & $ & "7I19, @& $ & 1994 ^!C!7B &7DD '% & 15)$ .&A& L 4& &99H' A& & &&& 4& L &B & 9 CB 415)$" $ ($" ?A 1 R9HS"A $" & & A 9I %01%"22"13"40"%01"$%"0224
PAGE 25
5&8T9' &8' &T8' = *" 474" 5" 19,K/; 9: %01%"22"13"40"%01"$%"0224
PAGE 26
) b & D & "* $ 019I"A $ D& L ?$ G ?*$?*$?* C 9*C M&'O & D & "* $ ?" $ &) && &19I'K 9 L $9*C $* !$* *! ,C)&UN@)K)*B $*C9! W>&ACGG>/@;,C)&UC;,UAB4 3 1A R9HS A !" "" & D & 3)& b3c & "*MB*" 9. %01%"22"13"40"%01"$%"0224
PAGE 27
; &TBT' &TB' &B' ;D &8' &T8' &8' 19IK1/; /! &TDBT' &TDB' %01%"22"13"40"%01"$%"0224
PAGE 28
& ?"$ & L 3)&& 19IK MM8L TMBTTMMTBT(&Z(9I' )E D & " & ?"$ & F 19:K &$ L9T *$/ T9TdT\eTT9'(&9(9:' 35 1" ?""A K";(7 ;(7?A ((! &;(7!' K K >T ($ TT & (L / 5K 2 9 &/ 9 C ($ CCBC &&// 9'L $ &9 M)5B *K T T9( L /) %01%"22"13"40"%01"$%"0224
PAGE 29
; ;D BT9 19:K1/; // %01%"22"13"40"%01"$%"0224
PAGE 30
#"[ "" ^CRES";( 7A ? ""WE("6!A 5"" "WE(5"? ;(7 5* 5f?!" &;(75' 1" & "(B $) 7 ;(754 5 8 "A 1" & "" & & DE %01%"22"13"40"%01"$%"0224
PAGE 31
"*(B&BOT9' $) A & ) $ (9 ) 5 $) & $ T9' ""$;(7 M! A """*" ""! $ ZT $ (BT9"? (;(7 *" $ #(;(7 ""*" ? A (( #" "[ ""A DH %01%"22"13"40"%01"$%"0224
PAGE 32
)4 0202 B""15)$ "$ (+" 15)$" )4 4)2:12D22 D$ ?$ *CLRQ9SDK 79915)$? 15)#"$ &9DG'(" "A ? 6+6H@=BG,@=B=& >H@=BO G =8 @ ) 4 /! B %( G12@FBC12@F/,B @)4/)B D< %01%"22"13"40"%01"$%"0224
PAGE 33
B ) L B T'&&I['T 12@1/A,BB& 4) LG &9DG'&9D9' b"A &9D9' "$ B D "A 15) ?1 15) @ 2)F6 OM` 15)$ D @ !1 A D @ ? "B ) 1 /* %01%"22"13"40"%01"$%"0224
PAGE 34
15) "P >)# "? E #! E $ G6 4"515) "&& "'DT& "'D'W" ($) ""& &"Q5'&"Q5'"&"Q5'$ "15)'1 H ,H 15) K) B 6 F6 IM 6 6 99 HG9)(D D
PAGE 35
13 dD 9 9HG9)(DE:.)(9DE:.)(9D 9 H99I)(9E<,)(D :.,)(E 9 : DED)(9 9
PAGE 36
""(" 9!"15) 6! 15) #" K15) """ $15) "1"15) 7" "+"15) D. %01%"22"13"40"%01"$%"0224
PAGE 37
)4 4/:21002%21:02? 12 *C """ #"B ? !"A &' ;(7 1 Q5P#5O&'LM&'9I 7 ;(7" 9 6RES ;(7A %" EG %01%"22"13"40"%01"$%"0224
PAGE 38
"R.S#" ( !&A "" %! "" 79D 1 #" ( R.S 6 B+" "" B&' ?( E9 %01%"22"13"40"%01"$%"0224
PAGE 39
" / F9T ) LG F9(8( LG T LM @=BGX@=B=8HH2H2 0@=BG =* #$ X5"4 1 9.4 &DD' M 4 $ ? ""& ;7$ /"1 :: GH 9,9, 9G EDED 9: ,H,H DG #( (?" $! ED %01%"22"13"40"%01"$%"0224
PAGE 40
19.K*3 EE %01%"22"13"40"%01"$%"0224
PAGE 41
/15)(A 15)" # RQ9SD "K GG< bGbG GG9 KUGbG 9GG $UGUG EEEE KbGUG "9 G&' 15)D<,D<,f" 15) $" 15) ;7$ F6 /F6 DD ED,)(9 .
PAGE 42
13 ;7$ ( DD:D:)(9 .
PAGE 43
B&'1 "" K" B +" [ !"" !" "B &99<' 199G LG[ ^C ) ($[ ^"C ^"C $" ^C R9S !R:S" E, %01%"22"13"40"%01"$%"0224
PAGE 44
g7 LL/ LLG 199GK^6"C^3C/ EI %01%"22"13"40"%01"$%"0224
PAGE 45
(" ?(" "/ "15)$( #(9"T5OLG >O6G,4 ?15) + I 3C 4>& $& & 4> $ 3 & 9 '(&'(&^'&h'h Ueib "A &9DD'3 G>KZ f D&/ i >'& < ('(&' E: %01%"22"13"40"%01"$%"0224
PAGE 46
415)$ ( ^C" ^C" &9DH'&9D<' &9D<' I &9D<'""?" &I!< & 4>&& TW'&9D,' &9D<' ( $&I!4>& & Lj $&I!> W&9(DI' B!" %!" &9D,'" 15) ("!R9
PAGE 47
I R9:S +Y60202 ? (" "G99GG B" ("" 15)!"^C %K !"^CVR9
PAGE 48
;7$ /"1 9,9, HE EDEDHH ,H,H H, ?" *!( "" $%"f" &"(" 9' 15)$ ( 79ED?fI 15) I (""!#" (""" "15) (" 79D!" 15) 79D (" H9 %01%"22"13"40"%01"$%"0224
PAGE 49
)4 4'02:0212 A "15)3 7E9 # ZLRQ9SD" Og=O9G= J%#0 QG= OQ + + OGGG9OG + C + M &"LG' D$' 'L 7) &DV ) 15)$""$ "" + A &"' %%& ;7$ F F6'F6 (=(!4'(!' .45'/+ 4/)*+ )*=)* !4/.!) /4)'!+ )4!5)+ /= / !4)/* .4!(+'/4(*(+' *'=*' !4!* )4''5+' 4/./+. '/ %01%"22"13"40"%01"$%"0224
PAGE 50
0 L9
PAGE 51
' RDS+ # LRGSD"" GUKU</A@/, A GUKU<
PAGE 52
15)$ H< %01%"22"13"40"%01"$%"0224
PAGE 53
/890%8%2 /4)"=%>:217 $799 ?"15)$ 1D9# I M +, =&& \ 5#&'O&'P N LG O # ["" #3C #F9"&'T5UW&'LG&DD' 5P"&'LM&'&DE' H, %01%"22"13"40"%01"$%"0224
PAGE 54
&TB' T&B' 1D9K;@; HI %01%"22"13"40"%01"$%"0224
PAGE 55
?15)$ $15)$799 "$ / *"&>b'&>'&>'&>' 1DD 6* 3 K@KN K!! C@KZKBKC@KZK C@KZKZKCKBK ,A @K] $$ C $ BKC@ZBC@ZZCBK+ 3 ""7" "&>'N 3 &>'N*"" 71DD""A $K&'N* 2&>'"&9G'* k&>'"&G9'* "k&>' "2&' $7 )4) 10"( &R
PAGE 56
& 1 / / K;@ '5 %01%"22"13"40"%01"$%"0224
PAGE 57
" "A #" B $ K [ ? 3 5 P8VK__k__" "P B__2__" #1DE ?B" 7A !1D9% *8 &$ =T ) B &$& VUC) " 7 )4) & & "K
PAGE 58
1DEK 1 7 )4) & & & ?" "&DE' &DE'"* #" <9 %01%"22"13"40"%01"$%"0224
PAGE 59
(B *B$" &' K &?MA /B*FB *B8(&9MD ) '&G9' *(B&G ) M / '&G9'*B 1DH"B ?"^"C*(B*B "?"" ^$C1B*B(*BB"5[B & ) '& ) M / ) 3$/$ &G9'&G9MD'*(B ?$799 3C$ ""? $"OP"3$) "3C$"
PAGE 60
1DHK5
PAGE 61
!(3C"" "!3 "&9G'"" &9I'" *" & *"&9:'" & 1 "3C "" "" & ;"" 1"2? 3 "1DH?! &DD'M2& ) "^C B7!&DD'2& ) "^C[B+ J """ & ? 7ED[ "4 !5"
PAGE 62
!3C"k1 "1DH3 K /9T J$ TE`TB *@D'$/$ C @#'$/$*$ C @0'$ C @>'$*C CZG P B "/ " 4 7E9 ^CH/ $" # /4/0662 /4/4)>2%03% 79D 15)$ $A ;"A "" << %01%"22"13"40"%01"$%"0224
PAGE 63
"1D< ) D%9 Q 15)$" 15)$ Q "A 4A 7A 9D# A 7EE +" "7 79D9MD #A <, %01%"22"13"40"%01"$%"0224
PAGE 64
3) 3M 1D
PAGE 65
" &*" *"' "A 3C #"" &'1D, "7 15) ""["" &_'1A "" 15) 4"A ("( <: %01%"22"13"40"%01"$%"0224
PAGE 66
23%:1>0 3%:1>0 1D,K/135 <. %01%"22"13"40"%01"$%"0224
PAGE 67
"" 15)$ " /4/4/+Y61" "B &' K K "BA ("79EE" ("A I 3C 1 I &9.'C ,G %01%"22"13"40"%01"$%"0224
PAGE 68
"" ($^C" ^C"15)$?! +"A "!15)$ ( "( $ +"" B"$^C" ^C"15)$ "! ^C"$ ^C"1&9.' "" I I A "! I I ("! ,9 %01%"22"13"40"%01"$%"0224
PAGE 69
/4 0202 /4 4)03%D1 7DD91 (( 115) ? ( ( LRQ ) ) SR( ) ) S" QD & D & $ & D< 1DI ?"15)$*C 7A DD"" ""A ;7$ /"1 9,9, 4)/ EDED 9I ,H,H 4// +";(7" A "$ ,D %01%"22"13"40"%01"$%"0224
PAGE 70
( !4/ (GH ( !4* ( !4( 1DIK;DGl3 ,E %01%"22"13"40"%01"$%"0224
PAGE 71
"A "! 1" ( ( "" Q< & R& 1DI 15) ;7$ /"1 9,9, GE EDED G< ,H,H G: A K 15)A ? "? "";(7 "" ,H %01%"22"13"40"%01"$%"0224
PAGE 72
(9(G:(G,(GH(GDGGDGHG,G:9 1D:K;
PAGE 73
( ( DGl% "G99GGO 15) ;7$/"1 9,9,GH EDED G, ,H,H G: # #" 15)A 7DDDA ""A +""" 15)A "1 K ( (
PAGE 74
,H,H' /"1 ) G: )! 9E )!! 9< )!!! 4/! 1A 15) & ? & '" !+" /4 4/+Y6D1 ("7DDD "15)1 ,I %01%"22"13"40"%01"$%"0224
PAGE 75
"^C 1"!"^C 15)$*CLRQ ) 9SRQ99S DGl"4 !15) #"" ("$ ;7$ /"1 9,9, GI EDED G: ,H,H G: "!" ^C"#7DDD I 3^C"! ? "!
PAGE 76
I " "15)$ 4!) 9GG"" ;7$ /"1 9,9, HE EDED HG ,H,H E: ?"*(!A (" 15)! 15) !("""A ("? Kf" 15)C C# ,. %01%"22"13"40"%01"$%"0224
PAGE 77
/4 4 1101"=%>:217 15)$ 1 $ "?*C G GI+ 9 & 9 S/F 6C>X^G!& 5("LM "P > L 4 1M 4 ) L&I'&DI'T 15)$$A LG 1D. "MV15) $""MV? 79E @ E """ 15)1MVL:GcIGc" ;(7"""1 IG %01%"22"13"40"%01"$%"0224
PAGE 78
@+ 9 & 9 B@ 9 & 9 B @+)&!B & ! &9G' @+ 9 &+ 9 B @ 9 &+ 9 B 1D.K/; I9 %01%"22"13"40"%01"$%"0224
PAGE 79
MVL,Gc
PAGE 80
"" 15)$ ?$ 9 +"!" 4 " [" #" S O =M ) =O 02 0 = Qg 0 OGGG9 g + C L H
PAGE 81
G< % %1;7$ 6 6 6 9,9,EG.D H
PAGE 82
&D<9' &99' "2h @+ 9 &+ 9 B &(<(9' @ 9 &+ 9 B 1D9GK*3 I< %01%"22"13"40"%01"$%"0224
PAGE 83
DD'( 15) ;7$ 6 ? '6 '=' !4/./' 4!*.+ (4/*5+ (=(!4) !/ )4'+ /4)5+ )*=)* !4!*.( '4.*+' .4(.5+' /= / !4! )4/) +')4.(/+' *'=*' !4!)' 4/ '+.'4/5'+. +(" ("" *C I, %01%"22"13"40"%01"$%"0224
PAGE 84
1"1:62 15)" 4)8%2%2D$ 79915)$A #&99'+ # "& '? U "!CUJ' U 3CA &9D'! "LQ#5UW II %01%"22"13"40"%01"$%"0224
PAGE 85
#(9"T G 4""L ) "(K &#'(9 J2U)* LG&' `&#'(9 2U T )* LG&"'&E9' LM&' 15)$ K"A 799 &E9'&9H'A ""# "7EE #" & $) 1E9 ?" M&&#'(9&Q U) T J$) G&ED' "" 799 "" O(A & I: %01%"22"13"40"%01"$%"0224
PAGE 86
1E9K5 I. %01%"22"13"40"%01"$%"0224
PAGE 87
"" / 799&ED' I&&#'(9' $&*0 T * ) !(6")/C) C*C*'<$ TT ) *&! gF\P' ?"5" 3$ 799 4/:6>128 90%2 4/4)1:021 +""A 3? "7EDD71ED "$3C :G %01%"22"13"40"%01"$%"0224
PAGE 88
1EDK5 :9 %01%"22"13"40"%01"$%"0224
PAGE 89
*[**" "3 j9"&'T5U&'LG&EE' #7D9!&EE' J$,) "^C"7 !&EE' @>,) ^C"? 3"6M ^"C""K M&(9"& ) T5O& )) P J,)* > &(9"& )*')) P @>,) L &EH' J$ "&K ) P2&[ )*J$'' ) P2& )* M`h"&[b )V &[b ) T @>>') e2&K ) """"" "" 1"""1 M"?" &EH'/" $')/,) $ 9 "G :D %01%"22"13"40"%01"$%"0224
PAGE 90
K & )W '__k&K'__" ) PIMB&KM'__2&K'__" &E<' D,& &K'PI&>i'__k&>'__ L ( 7 ( $ D,& D,& & ) P $N '__2&K'__LG ) 2?" )/,) L__&'____2&'__ )/X) __&'__&G&'' + 2k1ED7" __&['__L&&['P N)) M&G&[''&E,' :E %01%"22"13"40"%01"$%"0224
PAGE 91
*" MF ) &BM'P ,) LMGY9h&Ki'&Zii'^>i &EI' LMc&&& ) P N )) M&G& ))) __2& ) L YY!,) 8 UG__`____a__W4 ?" ,) !" &9G' ) 2 F T&m ( $'mT&m(m )? T&m(m ( mTm ) ) LDMGGT&DM9G(DMGG'mT&(DMGG'DMT&DM99(DM9G(DMG9T'm ?"2K 2L ) &&m9GQm'T&mQm9GQmG9Tm'& ( ( &DM9GQ / M'T&DM99j / MQDMG9TDMGG'DM' &EI' "*M 4" M !
PAGE 92
9 G&E<' DA& &'Pl&m'__2&B'__L &E<' $ Db& D,& ) P >&> / M'__k&K'__LG > k?" &$'M__&['__Lk&VM'M__k&KBM'__ &>'P ) L__&'__&G&'' 7" &>i'L&> )< l&$b ''M&G& )) &E.' *" !7)/,) M&F9__&'__&&VM'M__&K'__'P ,) LM @c) &&&m ) PBB&B )) M&O&K & )! P&k&m ) '__'P ,) Y/C,/ :< %01%"22"13"40"%01"$%"0224
PAGE 93
#" ""k K kL ) &&>(> !! 'T @K)) ( K)! ( K!) T> !! '>&E99' & /,!) Q / M'T& /,)) Q / MQ /,!) T / M !! '>'P &E9G' "*M " "$ +) &$ ) P2& ) *M* ""3C """"" 4""5UO ')/,) M (* %01%"22"13"40"%01"$%"0224
PAGE 94
?2" 85> )W 2&$ ) LM& ) P& ) L YI$ T M &E9D' L H L +88 < 4"" )*U*J *M "?&E9D' ) P ,) YY; 8MB /0 1 L YIZ$/Y&FU)/U ?""" "*[ :I %01%"22"13"40"%01"$%"0224
PAGE 95
"UW[" "*M" ) *&L +"*M !"8 1ED" > ?&>i ) P&>b ) #2331 "" *[* M J ? """ ) ML P) L P[) \ ) )<,)*J> ')/,) L P$*&/$$) +_DD_*M* (( %01%"22"13"40"%01"$%"0224
PAGE 96
4/4/2:021 1" ""UO?" # +) & )W 2& ) #"" ?" 1ED"" ? 9"G ?K "!C7 'P2&K ) / M'P "!,) LM&>i'&&>i'M&>i''P "!,) @@;d ON,K@Ke BB M&G& & ) P&2&'M__2& __'P#F2& ) L ,) P# !,) M8M#( ) 2&'P ,) B[( :. %01%"22"13"40"%01"$%"0224
PAGE 97
&E,'2[ "[ B?"A &E9E'&E:' 2""A B ? ) ?K M&# +) &K ) P2& ) L 7)<"!,) GU0& &>'&&>ii'M&>b''>#F2&K ) GU0&@@@Ke W' i''M&G&'''&E9H' P&k&'M__k&>'__'P "!,) L ) P "!,) "/<) P ,) &E.'k# " .G %01%"22"13"40"%01"$%"0224
PAGE 98
4"# # +) k(62!"! # " 4 238%328 28%32 # ? #7DD [!1EE "* J "*( $ L ) DEH 6? J "$ J U "$*( "K*" ?LO&BTD'TO& TH'&E9<' ?"" .9 %01%"22"13"40"%01"$%"0224
PAGE 99
1EEK/15 .D %01%"22"13"40"%01"$%"0224
PAGE 100
"K ) '$* J/ D' 3J E'" J) H' ; PVB__2__" &E9,' L9DEH$* J/ ?"&E9,'K .4 @/ &TD'T Ð''P 5# _Ue@O 6# C, & / OB # L9&9'Tn&Z'L9( "$ ? / "L ) *B*D" &E9,' ?&E9,'" K )P *\2&m 'LD2&K ) &E9I' .E %01%"22"13"40"%01"$%"0224
PAGE 101
+2&9G' *2( *(0&E9I' ?"&E9,'1 1EE?"K ?&K'P $X '__2&mM'__ L&W &K'T#O''PO&mb'&D__2&>'__' &E9:' L& &$b / M'P 5 i'__2&m'__T& $) PBM>'__2&>'__ LJ 9 "* L 7 +& ) *D"&E9I' & / $* Pb__ "" 1EE" K ?& ) P '__2&'__ L&? &m'O&K''D__2 ) &m'__&E9.' L &>i'(O_2&i'__L ?9J .H %01%"22"13"40"%01"$%"0224
PAGE 102
"&[ ) * 0&E9:'&E9.'" &E9,' &E9<' ? #3 "* 1EE K M# +) ?&m ) P2& ) *D*H ? ) *K / K ?L ? T_ / `1G/`/ K L _' 2L / 2 .< %01%"22"13"40"%01"$%"0224
PAGE 103
K ?O? 'P2& ) 5 #( ) ZJ) P ,J) *5$"!CZJ)/,J) Lo / #j ) &'e2 / &' *5"!C7) P ,) *T5D"/C7D)],D) #"A 3 3 1# !# ., %01%"22"13"40"%01"$%"0224
PAGE 104
'0%00;< 15)$" "%A 15)A (! (15) A "(" """ "" "15)#15) ""! 5%01%"22"13"40"%01"$%"0224
PAGE 105
" "$B RES?$ $ ?" "" 15)$ ?15) "! *! 15)A (""#A 15) 1" 15)!A "#15)" .: %01%"22"13"40"%01"$%"0224
PAGE 106
"A ( .. %01%"22"13"40"%01"$%"0224
PAGE 107
$3" R9S0)## 8)38?*A A Y"IYJ$Y@%F5 DKHEG(H%J66$K4#%&6 C9&Y3PYF5%#$F%KO6#6>%$>Y$F$%K <9(<.9..< RES# I$4>$6J&K$^6#CL_D4$67$&55$J$%K#%$ KF$J# ;3(74:<;! 3"7#9.:H RHS# 71/!8?0# &#;'38)" Y5>#$ K#"55$J$%K# /"*/9.:H R$`$K$6P6F6KI6&%# "49< Y5>$K46>Y6>$6#$K@%F5$%KI&6F$J# 7(59..9 R S?6 "I$4>$%>$ 7# !*9.:I RISp/871/!" Y"IQ/ F6>"K D:KE.D(HGD9..9 R ( S8)3 @%F53&# H:KE,,(E:,9.:D R.S8)371/!8?01077 (" 3>%J66$K4#%&6 C+& Y3PYF5%#$F%KO6#6>%$>Y$F$%K 9.(D<9.:. )!! %01%"22"13"40"%01"$%"0224
PAGE 108
R9GS/388#*"( (" 43#71/! Y$& @%556>I%K$K@%K6>6KJ6%KI$4>$I6&%# "/*EDDH 9<.(9IE+5#9..E4#7# R99S0))01+ A 3>%J66$K4#%&6&Y3PYF5%#$F%KO6#6>%$>Y$F$%K EH9(EHI9.:E R9DS0))101?7A A 3>%J66$K4#%&6&Y3PYF5%#$F%KO6#6>%$>Y$F$%K I9(:99.:E R9ES/671/!"&15)' 3>%J66$K4#%&6CC&K6>K$%K@%K6>a 6KJ6%K@`b ?5#9.:: R9HS71/! I$66"5$6I6&%#%>3>$b$6>6K$ P^$%K# ", `>%K$6>#$K"55$6I&6F$J# 7# *9.:. R9%$>Y$F$%K E,(9G,7# 9.:E R9.S8 Y>VK#6KF6>$^66#F6&%6#V66F6K#$K$# &U>$6#6F$6# *3C"*/ 9.II 9G9 %01%"22"13"40"%01"$%"0224
