MULTILEVEL METHODS FOR THE SOLUTION
OF ADVECTIONDOMINATED ELLIPTIC PROBLEMS
ON COMPOSITE GRIDS
by
James Scott Otto
B.S., Georgia State University, 1983
M.A., University of New Mexico, 1985
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
This thesis for the Doctor of Philosophy
degree by
James Scott Otto
has been approved for the
Department of
Mathematics
Thomas F. Russell
UXu
Daxe
Otto, James Scott (Ph.D., Applied Mathematics)
Multilevel Methods for the Solution of AdvectionDominated Elliptic
Problems on Composite Grids
Thesis directed by Professor Thomas A. Manteuffel
ABSTRACT
We analyze the use of multilevel algorithms in the solution of advection
dominated advectiondiifusion equations in one and two dimensions. In
particular, we study the behavior of the FAC (Fast Adaptive Composite)
scheme as applied to a composite grid discretization of the problem, where
the discretization type is allowed to vary (from upwind to centered dif
ferences, for example) on the components of the grid. In one dimension,
the analysis leads to the interpretation of a twolevel version of FAC as a
direct method. This analysis also provides insight into the behavior of the
algorithm, and guidance for its implementation, in two dimensions.
In two dimensions, we suppose the problem to have been transformed,
by an orthogonal coordinate mapping, in such a way that its character
istics are aligned with a Cartesian product grid. Various discretization
strategies for this problem are studied, with emphasis placed on a finite
volume method. Upwind, centered and higherorder upwind types of finite
volume discretization are considered. This method is attractive in that its
matrix stencils capture all of the advection of the discrete problem in their
main diagonal blocks. As a result, the block Jacobi relaxation is an excel
lent smoother for multigrid as applied to problems on uniform subgrids.
This allows us to formulate highly effective multilevel algorithms for the
in
global equations that preserve the inherent parallel nature of the strongly
advectiondominated problem. The ultimate numerical method developed
uses twolevel FAC as an algebraic solver in a nested way to obtain an
optimally efficient full multigridlike algorithm.
This abstract accurately represents the content of the candidates thesis. I
recommend its publication.
Signed
IV
ACKNOWLEGEMENTS
I wish to thank, first and foremost, my advisor, Tom Manteuffel, for
his support, academic and financial, that has allowed me to take part in
a variety of interesting areas of research which he has introduced me to,
including the one of this dissertation. Secondly, I would like to thank Steve
McCormick, who has been a consultant to this work since its beginning,
and whose expertise in multilevel algorithms has proved invaluable. I would
also like to thank my other committee members, Jan Mandel, Steve Pruess,
and Tom Russell, for their participation. Particular thanks goes to Debbie
Wangerin for her support, friendship, and sense of humor. Thanks, also, to
Bill Briggs and Roland Sweet for the initial support I received as a research
assistant in Denver. Finally, I wish to thank the people who nurtured
my interest in mathematics, those in the Department of Mathematics at
Georgia State University, and in particular, of course, George Davis who
inspired me to study numerical linear algebra.
v
CONTENTS
Chapter
1. Introduction................................................. 1
1.1 Notation ................................................ 6
2. Problem Description and Overview of Solution Methods ........ 7
3. Discretization ............................................. 16
3.1 Description of Composite Grids.......................... 17
3.2 Discretization Methods ..................................20
3.2.1 Finite Volume Discretization in One Dimension ..... 23
3.2.2 Discretization by Element Methods .................27
3.2.3 Finite Volume Discretization in Two Dimensions.....36
3.2.4 A Comparison of Discretization Methods ............ 41
3.3 Accuracy of the Discretizations ........................ 43
3.3.1 Sample Problems in Two Dimensions ................. 49
4. Convergence of FAC...........................................59
4.1 The TwoLevel Method.....................................60
4.2 Convergence in One Dimension ........................... 66
4.2.1 The Case of a Single Refinement Region ............ 67
4.2.2 Quality of Approximations ......................... 80
4.2.3 Multiple Refinement Regions ....................... 85
4.3 Convergence in Two Dimensions............................87
5. Alternative FAC Algorithms ................................. 97
5.1 FAC as Preconditioning ................................. 97
5.2 A Comparison of FAC and Domain Decomposition.............99
vi
5.3 Modifications Based on Multigrid....................... 103
5.3.1 The Significance of Line Relaxation...............104
5.3.2 FAC and Nested Iteration .......................... 114
6. Summary of Results ......................................... 125
References ..................................................... 132
vii
1. Introduction
The topic of this thesis is the efficient solution, by iterative methods,
of advectiondominated advectiondiifusion equations. Although we con
sider only elliptic second order equations, the algorithms presented in this
study may also be used to solve the problems that arise at each time step
in implicit methods for transient equations with nearly hyperbolic charac
ter. We consider solving equations whose true solutions are such that the
use of local grid refinement in their discretization is desirable. The type
of numerical methods we consider for solving the discretized problem are
based on multilevel techniques; the ultimate numerical solution method
that results from this study is a full multigridlike algorithm that preserves
the inherent parallel nature of the strongly advectiondominated elliptic
partial differential equation.
I
Multilevel, or multigrid, techniques have proven to be efficient methods
for the solution of elliptic partial differential equations. However, to a large
extent, the success of classical multigrid algorithms depend on the degree
of ellipticity of the problem along with the type of discretization being
employed. These two factors go handinhand, and what is more accurate
to say is that the success of the method depends on the degree of ellipticity
of the discretized version of the partial differential equation. Much of the
content of this dissertation is devoted to studying the effects that ellipticity
and discretization type have on the efficiency of conventional multigrid
methods in the context of an advectiondominated problem discretized on
a composite grid. Although we study a one dimensional problem, largely
as a means of gaining insight, particular emphasis is placed on solving the
1
twodimensional advectiondiffusion equation in a special form.
Following Chin and Manteuffel [14], in two dimensions we assume a
coordinate mapping has been used to transform the problem to one where
the flow velocity is aligned with gridlines in a Cartesian grid. As shown
in [14], with the problem in this form, line relaxation methods provide
robust solvers for the discrete problem. One way this work extends that in
[14] is by allowing multigrid, based on line Jacobi relaxation, to play the
role of the iterative solver. We also include the use of local refinement in
the discretization. Our choice of relaxation strategy allows us to retain the
parallelism of the strongly advectiondominated problem in efficient solvers
for both coarse and fine composite grid subproblems.
The particular form of multigrid we study is the FAC (Fast Adaptive
Composite) method introduced by McCormick [32] to taylor multigrid to
the structure of composite grid problems. We analyze the behavior of a two
level version of this method as applied to the advectiondiffusion problem
in one and two dimensions. A detailed analysis of FAC convergence in
one dimension yields an interpretation of the scheme as a direct method.
Perhaps more importantly, it provides insight in the appropriate choice of
grid and discretization type to be used on the coarse level. This information
is used in determining appropriate strategies for the use of FAC as applied
to the problem in two dimensions, whether it is used in its twolevel form
or as a full multigridlike algorthm.
A description of the contents of this dissertation follows. Chapter 2
provides an overview of the problem to be solved and the solution methods
used, and describes relevant related work. Chapter 3 introduces the com
posite grid and presents a summary of discretization methods for the one
and twodimensional problems. The structure and notation for the com
2
posite grid used in this study are given in Section 3.1. Section 3.2 describes
three discretization methods. A finite volume discretization is emphasized.
For the most part, our description of the finite volume discretization is
standard, however, we show how to adapt this method to the structure of
the composite grid, borrowing a technique from finite element and finite
volume element discretizations. These two discretizations are introduced
in Section 3.2.2 as completely systematic methods and also as means of
motivating the treatment of the interface portion of the finite volume dis
cretization. Although the former methods are more systematic, the finite
volume discretization offers advantages in terms of simplicity and ease of
use.
Section 3.2.1 describes the finite volume discretization for the one
dimensional version of our problem. This discretization serves as the basis
for the analysis of FAC convergence in one dimension later in Chapter
4. We notice that if the finite element subspace is chosen appropriately
this method yields the same discrete problem as the finite element and
finite volume element methods. The structure of a composite basis for this
space is introduced. In Section 3.2.2 the generalization of this basis to two
dimensions is made. Use of this basis allows the simple handling of the
nonuniformity at the interface when discretizing by the element methods.
In Section 3.2.3 the twodimensional finite volume discretization is intro
duced. This method serves as the basis for later studies of FAC convergence
in two dimensions. It is based on the socalled control volume approach
described in [36].
Section 3.2.4 compares briefly the relative merits of the various dis
cretizations treated in this study, and also identifies appropriate grid trans
fer strategies.
3
In Section 3.3 we discuss the accuracy of the discretizations in two
dimensions. We cite some recent results on the accuracy of finite volume
and finite volume element methods for diffusion problems, analyze the
accuracy of the finite element method for the advectiondiffusion problem
and provide examples demonstrating the accuracy of the particular finite
volume method used here.
Chapter 4 is a study of convergence of the twolevel FAC scheme ap
plied to the finite volume discretization of the one and twodimensional
advectiondiffusion problems. In Section 4.1 it is shown that, with a cer
tain compatibility assumption between the coarse and fine composite grid
operators, this method has the form of an interface preconditioning, i.e.,
the interface component of the residual satisfies
r?'= (I 
where, 2/ is the Schur complement in the composite grid operator, and
Aj is a preconditioner. It turns out that, for FAC, the latter is the Schur
complement in the coarse composite grid operator.
Section 4.2 studies convergence of the twolevel scheme in one dimen
sion, and contains the main theoretical results of this study. Our conver
gence analysis proceeds by an investigation of the spectrum of Bi = LjAj1.
In Section.4.2.1 this spectrum is analyzed in the case of a single refinement
region, allowing for the used of mixed discretization types on the composite
grid. We also take into consideration the effect of the type of discretization
used on the coarse grid. The analysis shows that good convergence results
can be obtained with significant jumps in refinement between the compo
nents of the grid, and, in particular, motivates the use of a global coarse
grid on the coarse level. In Section 4.2.3 these results are generalized to
4
the case of multiple refinement regions.
Section 4.3 is a study of the twolevel FAC scheme applied in two dimen
sions. Taking a cue from the onedimensional analysis, the global coarse
grid is used at the coarse level in specific FAC implementations for problems
that involve upwind, centered, and higherorder upwind descretizations on
the region of refinement. Generally, the results of this section indicate that
the twolevel scheme provides a good algebraic solver for the advection
dominated problem. However, better efficiency can be obtained by using
this scheme in a more sophisticated way.
In Chapter 5 we consider two ways of using FAC more efficiently. Sec
tion 5.1 views the twolevel scheme as a preconditioned iterative method,
and suggests using polynomial acceleration to improve its convergence.
This is the method of Schur complement domain decomposition and corre
sponds to the use of FAC by domain decomposition practitioners. Section
5.2 is a brief comparison of twolevel FAC with Schur complement domain
decomposition.
Rather than use the preconditioning perspective, that suggests solving
the composite grid equations algebraically, in Section 5.3 we take the more
efficient approach of using full multigrid to solve the equations to the level
of accuracy of the discretization. An important component of our imple
mentation of this approach is the use of strategies based on block Jacobi
relaxation for solving problems on the subregions. Section 5.3.1 contains
an analysis of this type of relaxation (as relaxation, and as a multigrid
smoother) for, upwind and centered discretizations of the particular prob
lem we study.
In Section 5.3.2 nested iteration and the twolevel FAC scheme are
combined to obtain a full multigrid version of FAC. Computational results
5
presented in that section show that, for the various discretization mixtures,
only one or two FAC iterations are required on each of the composite levels
to obtain an approximation by this version of FAC with accuracy compara
ble to the exact solution of the composite grid equations. Furthermore, we
show that this is the case with complete solves of problems on the subre
gions replaced by inexpensive approximations. This is significant because
practitioners of the twolevel approach have noted dramatic degradation in
global convergence when subproblems on the refinement region are solved
approximately. The method presented here, however, yields accurate solu
tions with optimal efficiency on this region, i.e., the total amount of work
there is comparable to that of a full multigrid iteration.
Chapter 6 is a summary of the results of this dissertation. It also
provides recommendations for their practical utilization.
1.1 Notation
For the most part, the following notational conventions are used in this
dissertation. Particular usage and exceptions should be clear from context.
A,L,X,,.. 
Xij 
Xij 
A, L....
Z, Fj
Z j Fj
/> & 
Hk
DetX
X^
*;xi .
Matrices or differential operators.
The ijth entry of the matrix X.
The ijth. entry of X~x.
Composite matrices.
Vectors.
Composite vectors.
Scalars or functions.
Finite dimensional function spaces.
Sobolev space of order k.
The determinant of the square matrix X.
The cofactor associated with X{j.
The j X j submatrix located in the lower right corner
of the square matrix X.
6
2. Problem Description and Overview
of Solution Methods
This chapter is devoted to describing the problem to be solved, and to
providing an overview of the solution methods and analysis used in this
study.
Although the theme of multigrid applied to advection dominated prob
lems is central to this study, there are a variety of ways of viewing this work
in which the significance of multigrid varies. One way that is essential is
to view the study as an extension of the work in [14] by Chin and Man
teuffel. There, the authors consider the solution of the twodimensional
advectiondiffusion equation
cAu(x,y) + auÂ£(x,y) + 6u5(x,y) + du(x,y) g(x,y), (2.1)
where e is a small positive constant. Notice that if e is equal to zero in
(2.1), then the equation is pure advection and is easily solved, in prin
ciple, by integrating along its characteristics. In practice, however, the
characteristics may be difficult to track and, in any case, we are interested
in solving the problem with e small but positive. Nevertheless, when e
is sufficiently small one feels that advantage may still be gained from the
problems nearly onedimensional character. In [14] the authors do this
by using an orthogonal coordinate transformation that is induced by the
characteristics. Under appropriate assumptions, the latter transforms (2.1)
into the equation
 eAu(x,y) + ux{x,y) + cu(x,y) = /(x,y). (2.2)
With this form, the characteristics are simply y = constant and any pro
cedure that employs solving along the characteristic directions will be fa
7
cilitated by having the characteristics aligned with gridlines in a Cartesian
product grid. In [14] the authors analyze the spectrum of the block SOR
iteration matrix associated with a centered difference discretization of (2.2)
(with c = 0). This iteration is based on alternative matrix splittings that
use either the x or yderivative terms of the discretization as the principal
part of the splitting. The relative size of the diffusion coefficient is mea
sured by the ratio 7 = h/(2e), where h is the meshsize. When the problem
is advection dominated (i.e., for 7 > 1), and with relaxation performed
along lines in the zdirection, they obtain a bound on the spectral radius
of the iteration matrix that is less than 1/5.
We also consider the numerical solution of (2.2). One way that this work
differs from [14] is that multigrid plays the role of the iterative solver. The
algorithms used in this study also retain the essential parallelizability of
the problem in this form by using block Jacobi, or line relaxation, as the
smoother on which multigrid is based.
Another way we extend the work in [14] is by including local refinement
in the discretization. The use of local refinement may be called for due to
various anomalies of the problem. These may include abrupt changes in
the righthand side [6] or the existence of boundary layers [40], [37].
Also, we emphasize the use of upwind rather than centered differences,
although we consider employing various discretization types on the sub
regions of a composite grid; use of the composite grid structure makes it
particularly easy to interface different discretizations on subregions. It also
allows one to interface appropriate solution techniques for these subprob
lems. The propitous aspects of concentrating effort locally in the context
of a composite problem are summarized succintly in [44, p.2]:
This approach offers a number of advantages. An existing code
8
can be upgraded, in order to increase the accuracy locally, with
out a radical redesign of the data structures, etc., since we can
use the old code to solve one or several standard problems. Is
sues of data structures and geometry are generally simpler if
we design programs for composite mesh problems in terms of
simpler standard problems. The use of simple standard prob
lems also tends to improve the performance of the programs on
vector machines. Finally, we note that if each of the standard
problems has appreciatively fewer degrees of freedom than the
composite model, then we might benefit from solving a number
of smaller problems rather than one large one. In other words,
we might view our approach as a divideandconquer strategy.
The composite structure also plays a fundamental role in the convenient
economization of multigrid to problems requiring local refinement. Multi
grid in this context is called the FAC (Fast Adaptive Composite) method,
which was introduced by McCormick in [32] and is treated in depth in
[33]. It is the particular form of multigrid studied here. The economiza
tion of the method lies in its avoidance of unnecessary computations with
respect to the coarse component of the grid. Its ease of use derives from
the fact that most of the computations involve solving familiar problems
on uniform subgrids.
Local uniform refinement has been abandoned by some multigrid re
searchers [2], [34], [38] in order to obtain more flexibility in using refine
ment. In particular, the methods of [34] and [38] allow for quite arbitrary
refinements, of which local uniform refinement can be seen as a special
case. When arbitrary refinement is allowed, however, the actual refine
ment typically is not specified before computations begin, but rather is
determined dynamically as computations proceed. Although such an ap
proach may be attractive because of its high degree of generality, it may
not be appropriate in certain contexts. For example, in many instances the
user has a good idea of the location and the degree of refinement needed
for a particular problem, so the effort expended in determining the refine
9
ment dynamically would be wasted. Because of the dynamic nature of the
corresponding algorithms, they are difficult to parallelize efficiently. Also,
because of the lack of a predetermined grid structure, analysis of these
methods is difficult. Finally, the use of arbitrary refinement may actually
limit the capability to represent a solution accurately in the sense that
such an approach necessitates the use of a single discretization method,
whereas it may be desirable to use methods of different orders in different
subregions of the domain, for example. The use of local uniform refinement
does not have these drawbacks. Also, the use of a composite grid allows
for the use of several levels of immediate refinement on a given subregion.
This can be important if the solution changes rapidly. We emphasize such
behavior in this study and investigate the effects that using multiple levels
of refinement has on FAC convergence. This is an important aspect of this
study, since FAC (albeit under different names) is currently being used in
this way to solve practical problems. For example, in [6] the use of the
closely related BEPS method in oil resevoir simulation is discussed. Mul
tiple refinement levels are used there to accomodate the representation of
injection wells by delta functions. We note that the algorithm appearing
there, although equivalent to the twolevel version of FAC, can be described
without reference to multigrid. Indeed, the authors view the method as a
preconditioned iterative method. Much of the content of this dissertation is
devoted to analysis of the convergence of the twolevel FAC scheme, so one
further interpretation of this work is as an examination of the effectiveness
of a preconditioner.
A composite grid (to be described more throughly in the next chapter)
consists of three components: an interface serves as the boundary between
the two true subregions of the grid, the coarse and fine patch components.
10
The latter component corresponds to a simply connected subregion of the
domain where local refinement is deemed necessary. We consider three
discretization methods in the context of the composte grid: finite element,
finite volume element and finite volume, emphasizing the latter. For the
most part, our description of the finite volume discretization is quite stan
dard. However, we show how to adapt this method to the structure of the
composite grid, introducing a novel way of computing patchtointerface
connections.
This technique involves interpolation of the interface (with the mesh
size of the coarse component) to the boundary of the fine patch, permitting
the use of a standard stencil throughout the region of the fine patch. The
technique is borrowed from finite element and finite volume element com
posite grid discretizations, where it occurs naturally. Although these two
methods are attractive due to being highly systematic, the finite volume
discretization offers advantages in terms of simplicity and ease of use. The
matrix stencils obtained with this method are, on uniform grids, familiar
fivepoint ones that correspond to either firstorder upwind or secondorder
centered differences. An advantage of using this method is that, with the
advectiondiffusion problem in the form (2.2), the advection terms are iso
lated in the main tridiagonal portion of the (appropriately ordered) matrix
stencil (this is not true for the elementtype discretization methods). The
significance of this is that, then, line Jacobi provides a robust relaxation
method and multigrid smoother that preserves the inherent paralellism
possessed by the problem as e 0. This is an important feature of our
implementation of the FAC algorithm, however it remains invisible until
well towards the end (Chapter 5) of this study. The reason for this is that
our analysis focuses on a twolevel version of FAC that proceeds by solving
11
subproblems on the components of the composite grid, and it is assumed
that these problems are solved exactly.
We study convergence of this twolevel version of FAC applied to the
finite volume discretization of (2.2) and its onedimensional counterpart.
We show that when a natural coarse component compatibility between
the coarse and fine level composite operators is satisfied, the convergence
of the twolevel FAC scheme is governed by convergence of the interface
component of the residual. The iterates for this component can be shown
to satisfy the relationship
= {i LjAj1)*.
__ A A
Here, Lj is the Schur complement in the composite grid operator, and Aj is
the Schur complement corresponding to a coarselevel discretization. This
is a distinguishing feature of the FAC method when viewed as an interface
preconditioning: it preconditions the Schur complement in the composite
operator with the Schur complement in a coarse grid operator.
Multigrid convergence theory for various discretizations on uniform
grids of the onedimensional advectiondiffusion equation has been devel
oped in [24] and in [3] using Fourier analysis. In [24] the author considers
the twolevel multigrid iteration, based on GaussSeidel smoothing, applied
to upwind and centered finite difference discretizations. For the upwind
case, a bound of 1/3 on the spectral radius of the iteration matrix is ob
tained, independent of 7. For the centered case, good convergence rates
are obtained by developing a strategy for introducing artificial diffusion
as the grid levels vary. Corresponding results are obtained in [3], where
a similar approach is taken with respect to a streamline diffusion finite
element discretization. Convergence of the twolevel scheme is proved for
12
sufficiently large artificial diffusion or a sufficiently fine grid. The results
of our analysis of the composite grid problem are similar in character to
these results for uniform grids.
Our analysis of the twolevel FAC algorithm for the onedimensional
problem discretized on a composite grid proceeds by an investigation of
A A
the spectrum of Bi = LjAj. We determine this spectrum in the case
of a single refinement region, and also extend the analysis to the case
of multiple refinement regions. This analysis covers various mixtures of
upwind and centered type discretizations on the components of the grids.
With upwinding is used on both levels, we show that the spectral radius
of the iteration matrix I Bi is zero in suitable parameter ranges. This
result applies ,to rather arbitrary choices of the grid used at the coarse level,
including a global coarse grid (i.e., the uniform grid over the whole domain
with the meshsize of the coarse component). For a discretization that uses
upwind and centered differences, respectively, on the coarse and fine patch
components of the composite grid, we show that with upwinding used on
the coarse level, the spectral radius of the iteration approaches zero as the
resolution is increased on the fine patch.
In our study of the twolevel FAC scheme applied in two dimensions,
the convergence behavior is seen to conform generally to the analysis for
the onedimensional problem. This is particularly true for a method that
employs centered differencing on the fine patch. It is shown that this
method works'well with sufficient refinement on the patch. Also, the results
of the onedimensional analysis may be used to develop a critereon for
switching to this method when the twolevel algorithm is used in a nested
iteration.
However, when standard upwinding is used on the fine patch, it is
13
discovered that the convergence rates deteriorate at the tangential bound
aries of the patch. The reason for this degradation is demonstrated, and
a successful remedy that involves modification of the coarse grid stencil
is implemented. Similar results are obtained for higherorder upwinding.
These results indicate that the use of higherorder upwinding can replace
standard upwinding on the fine patch with little sacrifice in efficiency.
Although we obtain satisfactory convergence rates for the twolevel FAC
scheme, we feel that significant improvement can be made in terms of
optimizing the efficiency of the method. For example, the analysis and
experiments that we perform with respect to this scheme employ exact
solution of the subgrid problems in its definition, and these solves must
be performed during each of the schemes global iterations. An obvious
and attractive alternative is to replace these exact solves with inexpensive
solution approximations obtained by using a small number of steps of some
iterative method. However, in practice it has been noted (in the solution
of diffusiontype problems, for instance) when exact solutions of problems
on the fine patch are replaced by partial ones that significant degradation
in convergence rates of the global process occur (see [23], for example). We
i
have noticed similar effects when applying this strategy to the advection
dominated problem here.
Our solution to this problem, which allows us to employ partial solution
on the subgrids, is to use FAC in a nested way, in the same way that the
Vcycle is used in the optimally efficient full multigrid method. Thus, two
level FAC plays the role of a weak algebraic solver on each of the composite
grids in the succesion of grids lying between the global coarse grid and the
true composite grid. This yields an algorithm that usually requires only
a few relaxations on the coarse mesh underlying the composite grid, and
14
which has optimal efficiency on the fine patch.
We close this section by noting that, although we have introduced our
numerical methods in the context of the elliptic problem (2.1), they are also
applicable to the solution of the mixed hyperbolicparabolic type equation
ut eAu + aui f bu$ + du = g, (23)
when e is snlall, i.e., when this problem is mainly hyperbolic in nature.
Although the results contained in what follows pertain to the elliptic prob
lem, they generalize well to the solution of (2.3) under the abovementioned
transformation of spatial coordinates when implicit methods such as back
ward Euler or CrankNicolson are used to perform the discretization in
I
time. The resulting linear systems at each time step then closely resemble
the ones described herein except for having improved stability properties
(that is, matrix diagonal dominance) due to the introduction of a positive
capacitance term.
15
3. Discretization
This chapter is devoted to studying discretization methods for elliptic
problems on composite grids. We begin in Section 3.1 with a description
of the composite grids to be used in this study.
In Section 3.2 we describe finite element, finite volume element, and
finite volume composite grid discretizations. Although the results of this
study pertain for the most part to the finite volume discretization, it is
useful to consider the subspace approach of the element methods. This ap
proach permits the systematic handling of the discretization on a composite
grid in more than one dimension, and also lends guidance in performing the
discretization in the absence of a subspace. In one dimension, it turns out
that if the finite element subspace is chosen appropriately the finite volume
method yields the same discrete problem as the finite element and finite
volume element methods. We introduce a composite basis for this subspace
in Section 3.2.1. The generalization of this basis to two dimensions is made
in Section 3.2.2, where we show how its use facilitates performing the el
ement discretizations. In Section 3.2.3 the finite volume discretization in
two dimensions is introduced. Based on the control volume approach de
scribed in [36], the method we describe in the context of the composite grid
is similar to the one in [22], where a finite volume method is applied to a
diffusion problem on a cellcentered composite grid (our grid differs from
this somewhat in that the interface nodes are not cellcentered). Our ap
proach is also similar to [22] in its treatment of the portion of the difference
stencil at the interface between the subgrids, though there is a significant
difference. In [22] these connections are calculated explicitly. This involves
16
I
either interpolating nodal values as piecewise constants throughout cells
centered at interface nodes, or interpolating, linearly, neighboring nodes.
Our approach is similar to the latter technique, however, it avoids cal
culating explicitly the stencils coefficients. Rather, it uses standard fine
grid coefficients (with the fine grid extended to the interface by using slave
nodes), and performs the required interpolation between interface nodes
prior to application of the matrix stencil as the solution algorithm pro
ceeds. This latter approach is borrowed from the element methods. We
note that it simplifies the derivation of the matrix stencil and allows the
use of a standard uniformgrid stencil on the fine patch. Section 3.2.4 com
pares briefly the relative merits of the three discretization methods treated
here.
Finally, Section 3.3 treats the accuracy of these discretization methods.
3.1 Description of Composite Grids
Our description of the composite grids to be used in this study applies
I
to the case of a rectangular domain in two dimensions. We suppose the
subregion of the domain requiring refinement to be an interior region that
does not intersect the boundary. Although no assumptions are made con
cerning the shape of this region, we will only consider refinement of the
grid on a rectangular patch. The structure we use for the composite grid
is that of [33]. Generalizations of this concept of a composite grid to L
shaped grids and refinement regions, other dimensions, refinements near
the boundaries of the domains, and multiple regions of refinement can eas
ily be made. But for the purposes of this study we restrict our attention,
for the most part, to the above simple case. Although we only describe in
detail the twojdimensional case, the structure of the composite grid in one
17
I
dimension is very similar, and it should be straightforward for the reader
I
to adapt the definitions in this section to the situation in one dimension.
Suppose we are given a rectangular hgrid with uniform spacing, h being
the nodetonode distance in both the x and ydirections. A composite
grid is determined by choosing some rectangular subset of nodes along with
a positive integer, m, indicating the order of refinement The set of nodes
forming the boundary of the subset is called the interface. We assume that
each side of tihe interface has at least three nodes (including comers) so
that the interior of the rectangular subset is nonempty. This set of interior
nodes is called the coarse patch.
If we designate the set of nodes exterior to the interface as the coarse
nodes, then wje have the uniform hgrid as the union of three components:
coarse, interface, and coarse patch. We call this grid the global uniform
coarse grid or just global coarse grid. See Figure 3.1.
The composite grid is obtained from the global coarse grid by adding
nodes to the coarse patch; the fine patch component of the composite grid
is the set of interior nodes of the uniform hjrgrid on the region bounded by
the interface, where hp = h 2ro. The set of boundary nodes for the fine
patch (which contains the interface) is called the fine boundary, but aside
from the interface nodes it is not considered part of the composite grid.
However, we define the complement of the interface in the fine boundary
i
as an adjuvant to the composite grid called the set of slave nodes. (Though
not technically belonging to the composite grid, they play an essential role
in the definition of the composite discretization.)
In summary, the composite grid is made up of three components, coarse
and interface components being the same as for the global coarse grid, along
with a fine component obtained via refinement on the patch. An example
18
0 0 0 0 o 0 0 0 0
0 , 0 0 0 o 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0 0
0 0 0 0 0 0 0
0 ' 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0
0 0 0 0 o 0 o 0 0
Fig. 3.1. Global coarse grid with coarse (o),
fine (), and interface () nodes.
19
of a composite grid (with one refinement) is shown in Figure 3.2.
0 0 0 0 o o o 0 o
o o 0 0 0 0 0 0 o
0 o 0 0 o 0 o o 0 0
o o
0 0 0 0 0 0 0
o o
0 o 0 0 0 0
o o
o o 0 o 0 o o o o o
o o o 0 o o o o o
0 0 0 0 o 0 0 0 0
Fig. 3.2. Composite grid with coarse (o), fine (),
interface (), and slave (o) nodes.
3.2 Discretization Methods
Next we describe three discretization methods for the following elliptic
problems:
*"(*) + u'{x) = /(), x e = [0,1] . .
ti(0) = 1, u'(l) = a,
and
eAu(x,y) + ux{x,y) = f{x,y), (,y) G fl = [0,1] x [0,1]
u(,0) =
(o ,y) = g2{y), (3.2)
u*(l,y) = ot(y),
u(x,l) = g3{x).
Here e is a constant, 0 < e < 1. The methods we consider are: the finite
element, finite volume element and finite volume methods.
We begin by giving a brief formal description of the two element meth
ods. Let a finitedimensional function space $ out of which one would
20
like to construct a solution to the partial differential equation Cu = / be
given. We take as our approximation to u a particular linear combination
of elements of a basis for $. In particular, we assume the existence of a
nodal basis: {^1, 2, 4>n} is such a basis if each i is associated with
precisely one node ni (we will frequently use this notation for nodes in one
or two dimensions), and, if any Â£ $ is expanded in this basis as
= i^i + 02^2 + aNNi
then Oj = {ni). The finite element (FE) method [1], [16], [26], [41] ap
proximates u:of Cu = f by the solution, 4?, of the discrete variational
equations
dv = j fi dV Vi.
This is equivalent to finding a such that Xa = f, where
< = <&,&)= / (C+ifaW
Jn
and fi fn fi dV.
To describe the finite volume element (FVE) [33] method, we require,
in addition to the finite dimensional space $, a finite set of volumes {Vi}^
that partition the domain Q. Typically, as with the basis elements, V{ will
be associated with precisely one node n, (a more detailed description of the
volumes for the composite grids used here will be given in the discussion of
the finite volume discretization). An example of volumes associated with
a subset of nodes of a two dimensional composite grid is shown in Figure
3.3. FVE requires that satisfy the local conservation equations
/ CdV= [ f dV Vi,
JVi JVi
which is equivalent to finding a such that Xa = f, where Zy = Jv. C(j>j dV
and fi = fy. fdV.
21
h
Fig. 3.3. Volumes for a portion of a twodimensional
composite grid.
We will return to the subspace framework shortly, but now we introduce
the third discretization method, beginning with a detailed description for
the onedimensional model problem. Similar to FVE, the finite volume
(FV) method [36] begins by partitioning the domain into a finite set of
volumes in such a way that, associated with each node of the composite
grid, there is precisely one volume that contains that node. Typically, the
node lies in the interior of its associated volume (although it is possible to
have a node at the edge of its volume, as is the case in the usual treatment
of Neumann boundary conditions). Figure 3.4 shows a portion of a one
dimensional composite grid containing a coarse node, an interface node,
and two fine nodes (nodes are labeled here) along with the endpoints
of their associated volumes (the ^s). The lengths of the volumes for
22
coarse, fine, and interface nodes are, respectively, h, h/2m, and (h+h/2m),
where m is the number of refinement levels used on the fine patch.
61
1
H0I1b
ti **' 6+1 *i+1 6+2 X'+2 6+3
h/2T
Fig. 3.4. Onedimensional composite grid with nodes (,,),
and volume endpoints (6>)*
3.2.1 Finite Volume Discretization in One Dimension
The finite volume discretization in one dimension proceeds by imposing
the following local conservation law, which is obtained by integrating (3.1)
over each of the volumes [6>6+i]:
rfc+i .. , r(t+1
J (~eu(x+ u(x)')dx = J f(x)dx.
(3.3)
These equations are then transformed into a system of equations for the
nodal values of u. Secondorder and firstorder terms are handled sepa
rately. For the secondorder term at the interface we have,
f u"{x)dx = u'(6+i) '(6)
^ u(a?t+i) u(a;f) u(xj) u(gj_i)
W 2(A/2m+1) 2(h/2)
= ( u(xi1) (1 + 2m)u(xi) + 2mu(si+i) )/h.
Here, we have:used a secondorder divided difference to approximate the
first derivatives. For the firstorder term
rfi+i
Ui
f + u,(x)dx = (6+i) u(6)
Ct
(3.4)
23
we need to approximate the difference on the righthand side by an expres
sion in u(x) evaluated at nodes. The centered formula averages u at nodes
adjacent to & and &+1:
u(xi+1) + tt(gj)
2
tt(gil) + u(xj)
2
u(ajti) ~ 'tt(gji)
2
The upwind formula replaces the volume endpoints with their nearest left
neighbor nodes:
(z) u(zi_i
)
Finally, one might consider a higherorder upwind formula for the first
order term on the fine patch. We do this by replacing the difference on the
righthand side of (3.4) by its secondorder approximation at Xi+2 and the
upwind nodes ;+i and ;:
*(&+a) (6+2)
1 3
u(xi) 2u(xi+1) + u(xi+2).
Notice that this formula cannot be applied at the leftmost fine node. At
that node, the upwind stencil can be used, for example.
Collecting these formulas and using similarly derived approximations
at coarse, fine, and right interface nodes yields the following tridiagonal
matrix stencils. Centered at the generic node n, these stencils involve
neighboring nodes in an obvious way expressed by
[raw, n, he]
Upwind
Coarse Nodes:
[e h, 2e + h, e]/h
24
[e h, (1 + 2m)e + h, 2me]/h
Left Interface:
, i
Right Interface:
[2me h, (1 ( 2)e ) h, ej/fe
Fine Nodes:
[2me h, 2m+1e + h, 2me]/fc
I
Centered
Coarse: Nodes:
[eh/2, 2e, e+h/2]/h
Left Interface:
[eh/2, (1 + 2m)e, 2TOe + h/2]/h
Right Interface:
[2me h/2, (l + 2m)e, e+h/2]/h
Fine Nodes:
[2me h/2, 2m+1e, 2me + h/2]/h
HigherOrder Upwind
Fine Nodes:
[h/2, 2me2h, 2m+1e+ Â§h/2, 2me]/A
25
This latter stencil involves the nodes
[nww, n, ve]
In addition, the Neumann condition at the right boundary, u'(l) = a,
is handled as'follows. Here, the volume for the node at the right boundary
is a halfvolume. See Figure 3.5.
1e1*
61 *i1 6 t = 6+1 = 1
h
Fig. 3.5. Onedimensional grid with a right Neumann boundary.
We have
/*
/ u"{x)dx = u'(xi) u\Â£i) w a (u(x;) u(xi_i))/h,
fXi
I u'(x)dx = u(x*) u(6) u(xi) u(xj_i) (upwind),
and
r*i 1
I u'(x)dx = u(si) u(6) ~ (i) ~ (w(x.) + u(s,_i)) (centered).
JÂ£i 2
i
These yield the following stencils and righthand sides.
Neumann Boundary
Upwind:
Centered:
[eh, e + h]/h = ^ /(l^) + ea
[efc/2, e + h/2]/h = ^ /C1  + ea
26
These stencil's, centered at node n, correspond to
[nW} ].
Generally, for the righthand side of (3.3), letting Vi = [&,&+i], we use
/ f(x)dx w (&+i b) f(xi).
JVi
Finally, we note that in treating Dirichlet boundary conditions, the
value of u at a Dirichlet node is given by the specified boundary value.
Thus, the equation for such a node is a trivial one that is easily eliminated
from the total system of equations by modification of its righthand side.
The use of this approach obviates the need for volumes associated with
Dirichlet nodes and, as a result, the volumes do not truly partition the
domain. A discussion of the effect of this treatment in terms of loss of
conservation in volume methods, along with an alternative for maintaining
the partitioning, is given in [33, Chapter 2].
3.2.2 Discretization by Element Methods
Reconsidering the FE and FVE discretizations discussed earlier, it is
interesting >to note that, in one dimension, the centered version of the stencil
just derived using the FV method can be obtained by applying either the
FE or FVE method (using the volumes just described) to the nodal basis
pictured in Figure 3.6c, where nr is the left interface node. (Only the
nonzero values of the basis functions are graphed.) The basis functions
corresponding to coarse nodes are familiar piecewise linear hatfunctions
for the uniform hgrid (i.e., the global coarse grid; see Figure 3.6a) and
the basis functions corresponding to the fine nodes are similar functions
corresponding to a uniform h j2mgrid (which we call the global fine grid;
see Figure 3.6b). We call these bases (Figures 3.6ac, respectively) the
27
global coarse, global fine and composite bases, and denote their respective
spans by $, and $. The relationship between the three bases may be
succinctly represented by the equations
iv =!V. + Vp + Â§V, = ^
0r = Â§0, + 0r = r Â§0,
0. = 0
Now, let $p be the span of {0, 0t, 0u} and let $<7 be the span of basis
functions for $ corresponding to nodes other than nt, nt, and nu. Then
A A
any composite basis element 0j may be written uniquely as 0j = itc + i,F,
where itc 6 and 4>i,f E Thus,
0j = 4>itc if ni is a coarse or interface node
i = 0i,f = 01 if nl is a fine node,
i.e., for coarse and interface nodes of the composite grid, the composite
nodal basis element 0 is obtained from the corresponding global coarse
basis element !0 by subtracting off its component in $jr (for coarse nodes
this component is zero), and, for fine nodes, 0 is the global fine basis
element associated with that node.
We note that this construction of a composite basis produces a space
that is nested ^between the global coarse and global fine spaces:
$ C $ C f. (3.5)
We call (3.5) the conformity relationship between these subspaces. It is
useful for obtaining accuracy estimates for the finite element discretization
(see Section 3.3).
This leads us to consider FE and FVE discretizations in two dimensions.
We need a finitedimensional space out of which an approximation to the
solution of our differential equation may be constructed. Suppose we define,
28
(a)
00
(c)
Fig. 3.6. Nodal bases for global coarse (a), global fine (b),
and composite (c) spaces.
29
in the usual way, nodal global coarse and global fine bases consisting of
piecewise linear hat functions on triangulations as shown in Figure 3.7.
Because of the similarity of the component structure (the division into
coarse, interface, and fine components) of the composite grid in one and
two dimensions the specification given above applies directly in two di
mensions to the construction of a nodal composite basis for such a space.
However, the process of discretization is now complicated by the nonuni
formity of nodes in two directions at the interface. Such nonuniformity is
not a problem at coarse or fine nodes where, as in one dimension, construc
tion of the matrix stencil proceeds along familiar lines (discretization on
uniform grids .with given meshwidths). Yet, the difficulties at the interface
can be made ;to virtually disappear by viewing the task locally from the
point of view :of the global fine basis.
To be specific, consider stencils involving interfacetofine connections.
These connections can be obtained easily by taking advantage of the struc
ture of the composite basis. For example, in the FE case the computation
of where q and p are composite basis functions corresponding to
interface node; nq and neighboring fine node np, respectively, is facilitated
by viewing q on the same (global fine) level as p (since p corresponds
to a fine node, it is equal to a global fine basis function), i.e., as a linear
combination of global fine grid basis functions (namely, q = q,c)~ This
reduces the prpblem to using existing finite element software to compute
(P,a), where p and B are both elements of the global fine basis, and
to combining such contributions. In a sense, this process uses slave nodes
 basis functions corresponding to these nodes are introduced temporarily
(see Figures 3:83.9). A similar technique can be applied in a FVE dis
cretization, though in practice the procedure is not quite as simple due
30
00
Fig. 3.7. Global coarse (a) and global fine (b)
triangulations.
31
to the additional nonuniformity in volumes at the interface. On the other
hand, FVE stencils are more sparse than FE stencils, additional localiza
tion being induced by the volumes.
Next we examine the finetointerface connections. We emphasize the
construction of this portion of the stencil in the subspace context since it
provides guidance for the construction in the absence of a subspace, for
example, in the FV discretization.
In order to describe the treatment of this part of the discretization, it
is useful to introduce the subspace $(,, which is defined to be the span of
those global fine basis elements that correspond to nodes at the boundary of
the fine patch (interface nodes and slave nodes this subspace comprises
a fine boundary for the patch), and $/, which is defined to be the span
of the composite basis functions corresponding to interface nodes. For
convenience, we summarize the decomposition of the global fine subspace,
A
\P, and the composite subspace, $, below:
 Span of ip's corresponding to fine patch nodes,
 Span of ips corresponding to other than fine patch nodes,
$(, Span of ip's corresponding to interface and slave nodes (1P& C
A A
 Span of
__ A
Here, ip and
nodal bases, respectively. We also note the following localness property
of the FE and FVE composite grid discretizations: each composite basis
function corresponding to an interface node enters into the computation
of a finetointerface connection only through its components in \Pj,. This
component of (pi (where n\ is an interface node) is obtained by restricting
32
the expansion
A
fa = fa,C
to the boundary of the fine patch:
h\b = hc\b (36)
Here, the restriction denotes the presence on the righthand side of (3.6)
of basis elements for only. For each fa, this relation provides a unique
representation of the component along the interface in terms of the basis
for ${,, and, so, provides a welldefined linear mapping, /}:$/>
of the interface to the fine boundary. We denote this linear operator as
interpolation at the interface. Now, due to the localness property, it does
not matter whether we compute the finetointerface stencil with respect
to the <^is at the interface or with respect to $6. It follows that this
stencil may be computed using just the basis elements from Wp and
(i.e., the nonuniformity between fine patch and interface may be avoided
by replacing the interface with a finemesh boundary). This means that
when the finetointerface stencil is applied (during matrix multiplication)
to an interface component, interpolation must first be performed:
Ub < IjUl
(ui is a given function in 4/, and ub is its component in the fine boundary),
and then the matrix stencil is applied to ub.
As an example of this process, consider the case of a composite basis
function fa at the left vertical edge of the interface when there is one level
__ A
of refinement. Then the representation of fa in terms of the basis for is
1 1 ,
% = jj V + Vi + 2^+
33
Fig. 3.8. Support of a composite basis function corresponding
to an interface node (O) at a left vertical interface.
Fig. 3.9. Support of a composite basis function (shown twice)
at a left vertical interface and support of the functions in its
component along the fine boundary.
34
(n_ and n+ are the slave nodes below and above node n{). The support
of is shown in Figure 3.8. It is shown again in Figure 3.9 along with
the supports of the three basis functions in its component along the fine
boundary. Representations along other edges are similar, meaning that
interpolation at the interface proceeds as a sequence of onedimensional
interpolations along the edges of the interface: the nodal value of a function
uj, equals the corresponding value of uj if the node is an interface node, and
is equal to the average of its neighbors nodal values if the node is a slave.
With multiple levels the process may be carried out using a reapplication of
the singlelevel process on the succession of levels leading from the interface
to the fine boundary.
We note that due to the localness property, this procedure is equally
valid for FVE and FE discretizations. Having arrived at a satisfactory
way of handling this part of the discretization when it is performed with
respect to a composite function space, it may also be used to facilitate this
part of the discretization in methods that do not utilize such a space, as
is the case, for instance, with the finite volume method. Normally, with
the latter method, treatment of the interface necessarily involves the use
of interpolation at interface nodes when calculating explicitly the fineto
interface stencil [22]. On the other hand, the approach described above
defines this stencil implicitly, and allows for the use of a standard uniform
grid discretization on the union of the fine patch and the fine boundary.
At this point, it is useful to say a few words about the way that the
coefficient matrix L of the composite grid problem is represented in prac
tice. Three matrix stencils are called for: a coarse stencil that involves
coarse and interface nodes, a fine stencil that involves fine patch nodes
augmented by a fine boundary, and an interface stencil that involves fine
35
patch, interface, and coarse nodes. In our computer codes, we use sepa
rate lexicographic orderings of the fine patch, and of the coarse component
imbedded in a global coarse grid. This means that the coarse and fine
stencils are standard uniform grid stencils and may be utilized as such. In
forming the fine patch component of the product Ly_, for example, storage
for the fine patch component of y is made large enough to accomodate
the fine boundary. Prior to applying the stencil to this component, the
interpolant of the interface of y is stored in these locations. In forming the
coarse component, computations proceed in the usual manner, except that
they are restricted to the coarse and interface components of y. Finally, the
interface portion of the product requires some special handling. In our im
plementation, the interface resides together with the coarse component as
part of an array representing a lexicographically ordered global coarse vec
tor. Therefore, interfacetocoarse and interfacetointerface computations
proceed as they would on a uniform grid, while interfacetofine connec
tions are handled by a separate code segment that requires access to the
patch component of y.
3.2.3 Finite Volume Discretization in Two Dimensions
For the most part, finite volume discretization of the problem (3.2)
is not much more involved than that of (3.1), since the method in two
dimensions is very similar to that in one dimension. To see this, consider
the portion of a twodimensional composite grid near a left vertical (west)
interface, as shown in Figure 3.10. The figure depicts (left to right) a coarse,
an interface, and three fine nodes along with their associated volumes.
Let V be the interface volume and label its sides in an obvious way as
N, S, E, and W. By Greens Theorem,
36
0 0
h h/2m
Fig. 3.10. Nodes and volumes near a left vertical interface.
I Au dV = I uxdy I uxdy
JV Je Jw
+ J Uydx J Uydx.
Consider, for example, the first term on the right. We approximate ux
at the midpoint of the east side of V using a centered difference formula
involving u(x,y) and u(x + h/2m, y) and integrate using the midpoint rule:
s / uxdy
Je
^ t (u({x + h/2m+1) + h/2m+\y) u{{x + h/2m+1) h/2m+1,y)\
~ V 2(h/2m+1) )
= 2m (u(x + h/2m,y) u(x,y)).
Using a similar approach for, say, the third term gives
/ Uydx h (^ + 7')
Jn v \2 2m+1/
u(x,y + h) u(x,y)
2(h/2)
37
= + ^sr) ((*> y + h) (*))
For the firstorder term, using an upwind approach, we have
I uxdV = / udy / udy
JV Je Jw
h(u(xty) u(x h,y)).
Let a =  + jraTT Using these approaches leads to the following stencils for
the upwind version of the FV discretization, which are centered at node n
and involve the nodes
[s.. nw, n, he, fijv]
Upwind
Coarse Nodes:
[e, e h, 4e + h, e, e]
West Interface:
[ae, e h, (1 + 2s + 2TO)e + h, 2me, ae]
East Interface:
[se, 2TOe h, (1 + 2a + 2m)e + h, e, ae]
South Interface:
[2me, a(e f h), (1  2a  2m)e \ ah, ae, e]
North Interface:
[e, a(e+h), (1 + 2a + 2m)e + ah, ae, 2me]
38
Fine Nodes:
[e, e h/2m, 4e + A./2m, e, e]
An obvious modification of the handling of the firstorder term (see
the details of the onedimensional FV discretization) yields the following
centered stencils.
Centered ;
Coarse Nodes:
[e, eh/2, 4e, e + h/2, e]
West Interface:
[sÂ£, e /i/2, (1 4 2s J 2m)e, 2me 4 hj2, se]
East Interface:
[se, 2roe h/2, (1 + 2s 4 2m)e, e 4 hj2, se]
South Interface:
[2me, s{e + h/2), (1 42s 42m)e, s(efc/2), e]
North Interface:
[e, s(e 4 h/2), (1 + 2s + 2"*)e, a(eh/2), 2me]
Fine Nodes:
[Â£, e h/2m+1, 4Â£, Â£ 4 h/2m+1, e]
39
We also provide the stencils for the Neumann condition at the right
boundary (with their associated righthand sides) and a line patch stencil
for a higherorder upwinding.
Neumann Boundary
Upwind:
^ 2 ^
[e/2, e h, 2e + h, e/2] = f(l ,y) + eh a(n)
Centered:
h,2 ^
[e/2, eh/2, 2e + h/2, , e/2] = y/(l ,y) + eh a(n)
Here, indicates the node tie lying outside of the stencil.
HigherOrder Upwind
Fine Nodes:
[~e, ^h/2m, e 2h/2m, 4e+^h/2m, e, e]
This final stencil is associated with the nodes
[ns, nww, nw, n, ng, n^].
As for the treatment of Dirichlet nodes, discretization at these points is
performed as described for the onedimensional case at the end of Section
3.2.1.
We note that the above formulas for the interface nodes do not apply
at its corners. For these nodes we use the standard stencil for the coarse
40
nodes. Also, the line grid stencils apply throughout the line patch, assum
ing that the process of interpolation of interface nodes to a line boundary
as described earlier in this chapter is incorporated into the definition of the
discretization at finetointerface connections.
3.2.4 A Comparison of Discretization Methods
In this section, we briefly consider the relative merits of the discretiza
tion methods described in this chapter. For the most part, the results of
this study pertain to finite volume discretization. We have chosen to em
phasize this method primarily due to its ease of use, first, in deriving the
composite grid equations, and second, in terms of solving these equations.
One attractive feature of this discretization is the simplicity of its ma
trix stencil. With the approach used here, this method gives rise to a
fivepoint stencil throughout the composite grid. This is due to the way
the interfacetopatch connections have been defined. The result is that
computer code involving the stencil is particularly easy to program. The
element stencils differ from this in that they distribute the coefficient of this
connection locally among nodes of the fine patch. The result is a stencil
that becomes less sparse at the interface as the patch is refined.
In the context of the problem,
eA u(x,y) +ux(x,y) = f(x,y), (3.7)
with e small, a more compelling reason for appealing to the finite volume
method is that its stencils appear on uniform grids as familiar finite differ
ence ones. The efficiency of the numerical methods of this study depend
strongly on the ability to solve discrete versions of (3.7) on the fine patch
and on the global coarse grid, i.e., on uniform grids. Relaxation serves as
the basis for solution of these problems by multigrid and, due to the fact
41
that the finite difference matrices capture all of the advection of (3.7) in
their main tridiagonal blocks, line relaxation is extremely effective (see Sec
tion 5.3.1). We note that it is more difficult to identify robust relaxation
strategies (i.e., ones that work well over wide parameter ranges) for the
element methods. This is due to the appearance of advection (correspond
ing to the xderivative in (3.7)) in portions of the stencil that correspond
essentially to discretization of yderivatives. For example, the uniform grid
stencils for finite volume elements, which is analogous to centered differ
encing, is
[eh/8, h/8, e 3h/8, 4e, e + 3h/8, &/8,e +&/8],
which corresponds to the nodes
[5j nsE, nw, n, he, knw, ^w]
This drawback does not preclude the use of the linear equations solvers of
this study when applied to the equations derived by the element methods.
However, more work needs to be done to develop the appropriate relaxation
strategies. Finally, we note that although the finite volume stencils appear
largely like finite difference ones, there is a methodology present in the
former approach that is lacking in the latter. This methodology provides
guidance for the systematic treatment of the discretization, especially as
the grid becomes nonuniform, i.e., as the coarse component meets the fine
one.
The systematic nature of the discretization is also an important and
attractive feature of the element methods, FE and FVE. In addition to
facilitating the treatment of the interface portion of the discretization,
it also simplifies the choice of interpolation and restriction operators on
42
which multigrid depends. With conformal subspaces, i.e., C 4% the
, <* A
interpolation operator, is determined by representing each
A A
basis element for $2h in the basis for In finite elements, restriction,
Ijlh : $/, > $2hi is also determined automatically by requiring that the so
called variational equations that relate the coarse grid discretization to the
fine one [30] hold. With finite volume elements the choice of restriction may
be guided by physical principles (see [33, Chapter 3]). With finite volume
elements the choice is more arbitrary. In this study, when grid transfers are
called for we employ the sevenpoint prolongation and restriction operators
of [25, Chapter 3] for the fine patch component, and the trivial (identity)
mapping: for the coarse and interface components.
Finally, we make note of the fact that upwind discretizations of (3.7)
play a very significant role in this study. We have found it easy to in
corporate firstorder and secondorder upwinding into the finite volume
discretization. In finite elements the use of upwinding is also common
place. On the other hand, a satisfactory strategy for performing upwind
discretization has not yet been developed with respect to the finite volume
element method.
3.3 Accuracy of the Discretizations
We now turn to a discussion of the accuracy of the discretizations de
scribed in this chapter. Of particular interest are the finite volume dis
cretizations that we have derived in some detail for the model problem
(3.2). The reasons for this interest are at least twofold. First, the finite
volume discretization as described for the twodimensional case is some
what novel and the matter of its accuracy is of interest in itself. And
second, the analysis of algorithms for solving composite grid equations
43
that will be pursued later in this study (in Chapter 4) is based on their
behavior .with respect to the finite volume discretization, so it is important
to establish the validity of this type of discretization. In the next section,
we will assess the accuracy of this discretization numerically for a sampling
of test problems. This section is a more general discussion of the accuracy
of the discretizations considered above. It begins by taking note of recent
convergence results for FV and FVE discretizations of elliptic problems on
nonuniform grids.
Prior to the recent studies of Cai et al. [10], [11], little had been done to
prove the convergence of finite volume element discretizations. Some prior
work had been done in this context for finite volume discretizations on uni
form meshes (references to this work are contained in [11] and [42]), and in
certain instances these turn out to be special cases of FVE discretizations.
The work of Cai et al., however, was the first to incorporate finite element
theory into the FVE framework in order to produce finite elementlike
accuracy estimates. This work also includes a systematic study of con
vergence on composite and other nonuniform grids. In [10] an 0(h3/2)
bound in a discrete energy norm is proved for the Poisson equation using
a composite grid FVE discretization similar to the one described herein,
with piecewise linear functions used for the approximating subspace (h
here denotes the coarse component meshwidth). This result is strength
ened to 0(h2) by incorporating bilinear functions into the discretization at
the interface. Similarly, in [10] an 0(h) result in the discrete H1 norm is
obtained for a diffusion equation on general triangulations, and this result
is also strengthened to 0(h2). This theory, however, does not apply to
composite grids.
More recently, Siili [42] obtained similar results for Possions equation
44
using a finite volume discretization. His analysis, which applies to Carte
sian product'grids, involves reformulating the finite volume method as a
PetrovGalerkin finite element method. This yields Oft*) estimates in the
discrete H1 norm when u 6 H1+
volume discretization for composite grids has been analyzed in [22]. This
corresponds closely to the method used here, but differs in using a strictly
cellcentered grid. There, a bound similar to the previous one is obtained
!
for the continuous H1 norm.
We note that for finite elements on composite grids, accuracy estimates
in the coarse component meshwidth are facilitated by the conformity rela
tionship (3.5). Consider, for example, problem (3.2) with the given bound
ary conditions replaced by u = 0 on all of dQ. Let be a global coarse
space defined by using piecewise linear functions on a uniform triangula
tions as indicated earlier. Associate with the corresponding global coarse
grid a composite grid via the procedure of patch refinement. Let be
the global fine space with the meshwidth of the fine patch, and let be
the composite space for the composite grid, so that
C h C Vh C Hi
(Note that we have allowed for the possibilities of trivial or total refinement
on the composite grid.) Define the bilinear form
a(u,v) = I eVuVv + uxv dV.
J n
and the L2 norm
I
IMII(n) = /2 W.
Then, using the familiar Galerkin approach, with u the weak solution of
A
(3.2) in Hq and Uh the discrete solution in $/,, we can use standard theory
45
to arrive at the H1 estimate
IK '/ijii(n) < cu y/iHi(n) V vh $h, (3.8)
where c = 4/e (is independent of h). The H1 norm is given by
IMI5pi) = E lllll>(n) = /(V V* +t.2) iV.
o<a
Now, since 4?* C then Vh on the right may be replaced by the
interpolant of u from This then can be used (see [41], for example) to
obtain an 0(h) estimate for u Uhfl'1(n) Also, Nitsches method [41],
[16] is applicable, yielding an 0(h2) estimate in the L2 norm.
We note that this result is somewhat unsatisfactory in light of the e1
present in the bound (3.8). However, if e is not particularly small the result
may be meaningful, and, in any case, convergence of the composite solution
to the weak solution is guaranteed as the coarse component is refined.
The above result is independent of the size, position, and level of refine
ment of the composite grids patch, i.e., it made no use of the inequality
$h C
It is clear that variations in the patch must affect the accuracy of the
discrete solution (our motivation for using local refinement is based on this
assumption). At one extreme, if the patch is reduced to a single point,
or if no actual refinement is used, then the estimates in h still hold. At
the other extreme, if the patch is allowed to cover the entire domain and
to > 0 refinements are used, then the above estimates hold with h replaced
by hp = h/2m. Using nontrivial refinement on a judiciously chosen patch,
we expect accuracy in the discrete solution that lies somewhere between
these two extreme estimates.
46
In particular, this will be true if u has small derivatives outside the
patch. Suppose we partition the domain into two subregions, 12 = tic U
tip, where tip is the rectangle with corners corresponding to those of the
interface. Then, for Vh G we have
u WfcHtfifn) <  VfcHi(n0) + u ^11^1(0^).
Now let Vh assume the role of the interpolant of u from the composite space
A A
$/,. Due to the particular form of v^, may be constructed using basis
elements for the global line space on the region tip, and basis elements
for the global coarse space on the region tic In tic this is true because,
by construction, each composite basis function associated with a coarse or
interface nodes agrees in tic with the corresponding global coarse basis
function. Therefore, we obtain the following standard estimate [16] in tic'
 ~ Vfctf(n0) ^ kch  u fl?(nc)j
where the H2 seminorm is defined as
I u lip(n)= X) ll^atilliJ(o)*
a=2
A similar result for the norm over the region I2/r would follow with h
replaced by kp except for one consideration: in order that the interpolant
A
over 12 be a member of $*, the interpolant constructed in f2j? is generally
not allowed to interpolate u at the slave nodes. This is because, in the
composite space, the function values at slave nodes are directly determined
by the function values at interface nodes (recall the process of interpolation
at the interface), and the latter are specified by the interpolant in I2cr> A
result like the one above using hp, then, may only be applied in the region,
denoted by tip, say, underlying the set of nodes associated with the fine
47
patch (recall Figure 3.2). More work is called for in order to obtain an
interpolation estimate in the region lying between Ap and Q.c Intuitively,
one should be able to obtain an 0(h) bound here. Rather than pursue
this, let us suppose that the component of the interpolant on (Ip does
actually interpolate u at the fine boundary. This will happen, for example,
if the derivatives of u are zero in the complement of flp. Then we have the
estimate
 Ufcjr(n) < c kch  u ffa(na) +c kphp \ u ua(nF) .
Notice that, if the derivatives are indeed zero outside of flp, then the first
term on the right is zero and we obtain O(hp) convergence as the patch is
refined. This'is the kind of behavior we expect, generally, if the derivatives
of u are small outside of some region and the patch is sufficiently large to
encompass the region. On the other hand, if u has steep gradients outside
of the patch, then the best we can expect is 0(h) convergence as the coarse
grid is refined. We note that a bound similar to the one above is obtained
in [22] for a finite volume discretization on a cellcentered composite grid.
Although these observations have been made with regard to finite ele
ment discretization, we expect similar accuracy for appropriate finite vol
ume and finite volume element discretizations. In particular, the above
example is analogous to the centered version of the finite volume discretiza
tion used. here.
In the next section, some examples are presented that indicate the ac
tual convergence behavior of finite volume discretizations of problem (3.2).
Specifically, we consider the use of upwind or centered approximations on
the coarse component of the grid and allow for upwind, higherorder up
wind, and centered approximations to be used on the fine patch. Generally,
48
out computational results agree with the behavior predicted by the above
analysis. However, there is still much work to be done in terms of estimates
for the accuracy of composite grid discretizations. In particular, we have
found that the use of upwinding here gives rise to accurate solutions and
very fast discrete solvers. Therefore, it would be appropriate to pursue its
use with respect to the element methods. One attractive possiblity is that
of applying the streamline diffusion method [26] in this context.
Before presenting the results of our test problems, we make one final ob
servation with respect to accuracy and the relationship between the coarse
and fine components in the context of the advectiondominated problem.
We have just observed that when the derivatives of u are small outside
the patch, an accurate solution can be obtained with sufficient patch re
finement. Although this is the case, one should not expect inaccuracies
to then be restricted entirely to the fine grid. Because the problem is
advection dominated, the solution at a given point is dependent on the so
lution upwind from that point. A similar dependency should be expected
in the discrete solution. In particular, the solution on the coarse com
ponent downwind from the patch will inherit the innaccuracy that is a
result of using insufficient refinement on the fine component. Therefore, it
is possible for the true solution to be uniformly wellbehaved throughout
the portion of the domain corresponding to the coarse component, while
the accuracy of its discrete counterpart changes significantly in this region.
This behavior of the discrete solution will be observed in the test problems
of the next section.
3.3.1 Sample Problems in Two Dimensions
The following test problems are constructed from (3.2) by determining /
49
and the boundary data from the specified solution. We employ constants
t > 0 and 0 < x0, yo < 1. With respect to the point (xo>3fo)) define
D(x,y) = yj(x *0)2 + (y y0)2.
Problem 3.1.
The solution, is
= 2 ((Â£(*,y)/r)3 ^(D(*,y)/r)2), D(x,y) < r,
u(x,y) = 1, V(x,y) > r.
This solution has been designed to be smooth at the interface between
its two components. Specifically, at the set of points D(x,y) = r, u is
continuous and differentiable, and its gradient is zero:
= 6(* x0)(D(x,y)/r 1 )/r2,
= 6(2/ y0)(D(x,y)/r 1 )/r2.
We note that u is independent of e, and the extent to which the function is
wellbehaved depends on the magnitude of r. When this quantity is small,
the nonconstant component of u covers a small region, however its gradient
there is comparable to 1/r. On the other hand, if r is kept bounded away
from zero, then we have a moderately wellbehaved solution irrespective of
the choice of e. The graph of u(x,y0) with Xo = yo = 1/2 and r = 1/4 is
shown in Figure 3.11; the graph of u(x,y) is obtained by rotating it about
the uaxis at (xoj^o)
Problem 3.2.
The solution is
u(x,y) = e~D
50
We remark that u and its derivatives are small except in a neighborhood
of (z0,i/o), where the function undergoes rapid growth. For simplicity, let
r = 1 and o. = 3/0 = 0. Then
ux = u 2s/e, uxx = (2u/e) (2s2/e 1),
so that u has its maximum rate of change at D(x,y) = y/e, with rate
Vu = 2/(e y/e). This represents a sharper front to be resolved than in
Problem 3.1.
Convergence rates as functions of grid refinement for these test prob
lems are presented in Tables 3.13.5. Relative errors involving the discrete
solution of the problems and the known solution evaluated at the nodes
are presented. The discrete maximum norm,  Us,, is computed on the
composite grid. The discrete L2 norm,  2, is computed on a global
fine grid having the same meshsize as the fine patch. The errors are com
puted on the latter grid by interpolating the discrete solution to it from
51
the composite grid, while the true solution is simply evaluated there. Var
ious combinations of finite volume (upwind and centered) discretizations
are used on the coarse and fine components of the composite grid. In ad
dition, in some instances a discrete solution is computed on a global fine
grid with the same refinement as the fine patch in order to compare the
composite grid solution with a uniform grid solution. The discrete solution
on all grids is computed by an iterative method and using the convergence
criterion that the relative residual satisfy
ri2/r2 < 106.
Here, r* is the ith composite residual and r is the original righthand side
of the composite grid equations. The iterative solver used for these tests
is a particular implementation of the twolevel FAC algorithm described
in Chapter 4' of this study (see, also, Section 5.3.1). This is not the most
efficient algorithm presented herein, but its use is appropriate for these tests
in that the exact solution of the composite grid equations is desired. An
extensive treatment of the efficiency of this algorithm and alternative ones
is presented in Chapters 4 and 5.
In what follows, the meshwidth of the coarse component is denoted by
h = 1/2I*, for some positive integer me. A positive integer m denotes the
additional number of refinements in meshwidth (each by a factor of two) in
going from the coarse component to the fine patch. Hence, the fine patch
has a meshwidth of hp = h/2m.
Tables 3.13.2 show the results for Problem 3.1 with r = 3/32, x0 =
1/2, and yo = 3/8. The coarse component meshwidth, h = 1/32, is kept
constant and the patch is successively refined. The southwest corner of the
I
interface is located at (12/32, 7/32), with x and y dimensions of 8/32 and
52
10/32, respectively. The patch here is made large enough to encompass the
nonconstant component of u. Upwind discretization is used on the coarse
component of the composite grid, and two test results are shown, using
upwinding and higherorder upwinding (Hup) on the patch. The diffusion
coefficient is e = 104.
Table 3.1
Relative errors in the maximum norm as functions of
patch refinement for Problem 3.1.
m Composite Grids Upwind/Upwind IMoo/IML Upwind/ Hup M./NL Global Fine Grid Upwind IkilU/NU
1 1.383 10"1 5.551 10"2 1.383 10"1
2 6.569 102 1.549 10"2 6.568 10~2
3 3.194 102 3.977 lO4 3.194 10~2
4 1.576 10"2 9.605 10"4 1.576 10"2
5 7.828 103 2.273 104 7.808 10~3
The error norms show O(hp) and 0(hp) convergence for the upwind
and higherorder upwind discretizations, respectively. In addition, the re
sults for the upwind composite solution and global fine solution are nearly
identical^ with the composite solution showing a mild degradation in the
X2 norm.
Table 3.2
Relative errors in the L2 norm as functions of
patch refinement for Problem 3.1.
m Composite Grids Upwind/Upwind Mi/IHI* Upwind/ Hup Iklk/IMh Global Fine Grid Upwind
1 1.069 102 3.670 103 1.095 102
2 5.359 lO"3 9.532 lO"4 5.422 10"3
3 2.687 10~3 2.500 10"4 2.703 10~3
4 1.346 103 6.589 10"B 1.350 lO"3
5 6.735 104 1.663 10B 6.750 104
53
Tables 3.33.4 show the results for Problem 3.2 with r = 1, x0 = yo =
1/2, and e = 103. Standard upwinding is used on the coarse component
of the composite grid, with meshwidth h = 1/32. The higherorder upwind
and centered discretizations were used on a patch centered at (o}3fo) with
length and height equal to 1/4. Again, the patch is made sufficiently large,
this time to enclose the solutions sharp front. The problem is also solved
using the centered method on the global fine grid. The convergence rate
for the discretizations is 0(hp), again with close agreement between the
composite and global fine solutions.
Table 3.3
Relative errors in the maximum norm as functions of
patch refinement for Problem 3.2.
m Composite Gr Upwind/ Hup IMoo/IMloo ids Upwind/ Centered IM/IMI Global Fine Grid Centered M/IMI.
1 9.538 10"2 1.051 lO"1 1.050 lO"1
2 3.271 Hr2 2.253 10~2 2.253 102
3 9.279 10"3 5.460 103 5.460 103
4 2.397 10"3 1.356 103 1.356 103
5 6.021 104 3.383 10"4 3.383 lO"4
Table 3.5 shows the results for a different set of tests based on Problem
3.2, measuring convergence as a function of h. Here, 6 is increased to
102, allowing the front to spread out while the location and dimensions
of the interface are left the same as for the last set of tests. Again, the
upwind/centered and centered/centered combinations are used, but now
on three composite grids, each having two levels of patch refinement and
varying coarse component meshwidths.
54
Table 3.4
Relative errors in the L2 norm as functions of
patch refinement for Problem 3.2.
TO Composite Gr Upwind./ Hup Mi/IMli ids Upwind/ Centered Global Fine Grid Centered Iklh/IMh
1 1.038 10"1 8.449 102 8.436 lO"2
2 3.159 lO"2 1.858 lO"2 1.857 10"2
3 8.533 lO3 4.519 lO"3 4.520 lO"3
4 2.184 lO'3 1.122 lO"3 1.123 lO"3
5 5.504 104 2.800 104 2.802 10"4
We make the observation that the patch is too small for this problem,
in the sense that the innacuracy on the coarse component is inherited
by the patch (and also affects the coarse component downwind from the
patch). As a result, the effectiveness of increasing the patch refinement
has deteriorated somewhat. For the upwind/centered discretization with
mc = 4 and to = 2 (second row, first column of Table 3.5), the relative
Table 3.5
Relative errors as functions of coarse component
refinement for Problem 3.2.
77lc (tJX 2) Composite Grids Upwind/ Centered Mco/IMloo C entered/ C entered MIco/IMloo
3 1.963 lO1 1.935 lO"1
4 7.981 10"2 7.468 102
5 3.921 10"2 1.903 lO"2
6 2.152 10i 5.628 10"3
L2 error restricted to the fine patch is approximately 7.67 10 2. With
two additional levels of patch refinement (to = 4) this becomes 3.78
55
102, a modest improvement (in particular since the secondorder accurate
centered discretization is used on the patch). On the other hand, this
additional local refinement has a pronounced effect on the solution as a
whole, as is apparent from graphs of the solution presented in Figures
3.123.14. In Figure 3.12 the exact solution of Problem 3.2 is plotted.
Figures 3.13 and 3.14 show plots corresponding to the two cases of patch
refinement just mentioned. In each case, for ease in plotting, the fine
patch solution has been transfered by restriction (see Section 3.2.4) to a
coarse patch having the same refinement as the coarse component. We
emphasize the quality of the coarse component solutions here. Recall that
the solution should become uniformly small as the radial distance from the
center of the domain increases. In Figure 3.13, this fails to happen due to
the innacuracy of the coarse component solution in the region downwind
from the patch. In Figure 3.14, this innacuracy of the solution has been
removed by using purely local refinement upwind from this region, i.e., by
adding additional refinement on the patch. This example shows that, in
the context of advectiondominated problems, local refinement can have a
profound effect on the global solution.
56
ST
58
4. Convergence of FAC
This chapter studies the convergence of FAC as an iterative solver for
advectiondominated advectiondiffusion equations on composite grids. We
begin with only a few assumptions regarding the structure and origin of the
composite grid equations. The composite grid consists of coarse, interface
and fine patch components, and coincides with a global coarse grid in its
coarse and interface components. Its fine patch may be viewed as being
derived, by refinement, from a coarse patch on the global grid.
As for notation, members of the composite grid space are represented by
lower case Roman letters distinguished by an underscore, and components
by the lower case letters subscripted with C, F, or I. For example, a
composite righthand side (or residual) is written
r = [rc,rF,r/].
Similar flotation is used for composite grid operators and their components.
The composite grid equation Lz = r may then be written as
Lc 0 La z c ' r c "
0 Lf Lfi Zp = rj?
Lie Lif Li . %i . . r* .
A word of explanation with regard to the notation accompanying the com
ponents of L: in general, Lxy denotes that portion of the stencil of the
matrix L that represents the connections to nodes in component Y that
appear in the equations for the nodes in component A, i.e., entries where
node i lies in X and node j lies in Y. Also, Â£/, for example, stands for the
59
interfacetointerface connections. Section 4.1 below introduces the two
level version of FAC for general composite grid equations. The results in
this section utilize only the algebraic structure of the composite equations
and are not at all dependent on the underlying continuous problem or its
discretization. Sections 4.24.3 are devoted to studying the convergence of
this algorithm when it is applied to the finite volume discretizations of the
advectiondominated problems (3.1) and (3.2), respectively.
4.1 The TwoLevel Method
Assuming that the matrices L, Lc, and Lp are invertible, then per
forming a few steps of block reduction leads to the equivalent form
Lc 0 Lei ' Xc ' r c
0 Lp Lpi Zp = Tp
0 0 I . . . LjHj .
where
Li =
Li LicLcxLci LipLp Lpi
i
(4.2)
is the Schur complement in L with respect to the partitioning (4.1), and
r/ = r/ LicLc1tc LipLFXTp.
In principle, then, the composite grid equations could be solved directly
by the following threestep procedure.
Step 1. Solve LiZi = ii for zj.
Step 2. Solve ,Lc%c = *c for zc.
Step 3. Solve Lpzp rp LpjZj for zp.
We note that this is just the familiar Schur complement method for solving
a system partitioned as in (4.1).
60
FAC [32], [33] is an adaptation of multigrid to composite grid problems
and, as such, makes use of a coarse grid discretization of the differential
operator to obtain approximations to the fine grid error left by relaxation.
In the present context, one may think of the fine grid operator as being
the composite grid operator defined above. As for the coarse grid operator,
it is useful in defining it to retain some flexiblity. We accept in this role
a discretization on any of the composite grids that lie between the global
coarse grid and the true composite grid (the grid on which the discrete
solution is ultimately sought). These grids are the ones obtained from
the global coarse grid in the same manner as the true composite grid, but
differ by having less refinement on the patch (all grids agree in their coarse
and interface components). An important special case to consider is the
global coarse grid viewed as a composite grid. Since this grid is uniform,
in principle, no special considerations need apply in defining the operator,
that is, a standard discretization based on a lexicographic ordering of the
nodes and equations, for example, may be used. Also, since the grid is
uniform, it is possible to apply any one of a variety of wellknown, highly
effective methods to solve standard problems on this grid.
The uniformity of the global coarse grid is its attractive feature. How
ever, we retain the viewpoint put forth while defining this grid (see Section
3.1) that it is itself a composite grid and that the discretization on it yields
a composite grid operator. It makes sense, then, to refer in general to the
coarse grid as the coarse composite grid (the phrase composite grid may be
reserved for the true composite grid when distinguishing the levels in the
twolevel algorithm). Once this viewpoint is taken, the strategy of using
different discretizations on the coarse and patch components of the coarse
grid becomes a possibility for all choices of this grid. We have already
61
considered, in Section 3.3.1, the effect on the accuracy of the composite
grid solution of using different discretizations on the components of the
composite grid. We investigate in this chapter the effect of using various
discretizations; on the coarse grid in terms of FAC convergence.
When denoting the grid operators, L is used as above in denoting the
composite operator and its components, and A is used similarly for the
coarse (composite) grid operator. Generally, notation developed for global
coarse objects is used to refer to coarse composite ones. So, for example,
the refinement region on the coarse composite grid is referred to as the
coarse patch.
Let A be the discrete operator on this grid. Its component form is
A =
Ac 0 Aci
0 Ap Api
Aic Aip Ai
Here, the subscript P refers to the coarse patch, as opposed to F, which
refers to the fine patch component of the (true) composite grid.
I'
The following description of FAC corresponds to the delayed correction
version in its twogrid form as described in [33]. Here, a subscript G is
used to distinguish coarse composite vectors and components from true
i
composite ones.
FAC algorithm for the solution of Lz = r
Let z be given and set r = r Zz.
Loop on i: i = 0,1,... until convergence.
Step 1. Solve Aza = tq = (r^,rp = IFrlF,r}) for
zg = (zg1c,zg,p,Zgij).
62
Step 2. Solve Lpz'F = rlF LFiZg,i for z'F.
I
Step 3. Perform the composite correction
zJt"1 = z*c + zc,c
z^1 = zjp + Zjf.
zjr+1 = Z> + ZG,J.
.1 l
Step 4. Set zl+1 = (zp1,z^1,z}+1), and form the new composite residual,
r*+i = r Zz+1.
End loop.
To express the algorithm in words, each iteration begins by transferring
the current residual to the coarse patch (dropping the iteration indices):
Tp * Iftf. Here, IF represents a restriction operator mapping the fine
patch to the coarse patch. Then, Tp and the other components of the com
posite residual form the righthand side for the coarse composite equations.
This system is solved to obtain zc, a global coarse (or coarse composite)
grid approximation to the composite grid error. To obtain a fine patch
error approximation, the local patch problem is solved with the fine patch
residual las the righthand side. Dirichlet values are provided for this prob
lem by the interface values of za (As a practical matter, recalling the
discussion in,,Section 3.2.2 on the treatment of the finetointerface stencil,
in order to form the vector of boundary values in that context, the inter
face component of z& is necessarily mapped to a refined interface using the
interpolation operator /*.) The solution z'F of this problem, along with the
coarse and interface components of zg, is used in the final step to define a
63
composite correction, which is added to the current approximation to the
solution. Finally, the residual is updated.
Having defined the FAC algorithm, we now turn to a study of its conver
gence properties, concentrating on an expression for the updating of the
composite residual. In what follows, we assume appropriate consistency
[17] for i, A, and their components.
THEOREM 4.1. Let the composite operators A and L be coarse grid
compatible in the sense that Ac = Lc and Aci = Lei, and let be the
ith composite residual obtained by using the above twolevel FAC scheme.
Then for i > 0, we have r}?, t1f = 0 and
rj+1 = (/ i/ij'K,
(4.3)
where Li and Ai are the respective Schur complements in L and A.
Proof. According to the Schur reduced form of the coarse composite
1 '
equations, Step 1 of FAC (the coarse composite solve) yields z
(Ai AicAq1 Aci AipAp1 Api) zg,i = r) AicAq1 AipAf,11pPF,
! (44)
and
i
Ac%g,c = t*c ~ AciZgj
Also, the update of the composite residual may be written as ri+1 = r*
Lz'c, where
zc = (zc.c.z^zc,/).
Using these relationships along with (4.1), then the residual update may
be written in component form as
rc^ =Tc ~ Lc%g,c ~ Lci%g,i
= r c Lc(Ac1( rlc Ac/Zg,/)) Lci^gj,
(4.5)
64
I
A I,
of I LiAj1. Although we only used the assumption that Ac = Lc
and Aci\ = L>ci above, it is also appropriate to make the assumption for
interfacetoc'oarse connections, i.e.,
i 1
i
Ajc = Lie (4.8)
Let (A,x) be an eigenpair of LiAj and set y = x, so that the equation
A A 
LjAj x = Ax
becomes:
(L/ AA/)y = 0.
The eigenvalues of LjAj1 are therefore those values of A for which the
matrix 
Li ~j'AA/ + (A 1)LicLc1Lci LipLp1 Lpi + XAipAp1 Api (4.9)
I, i
is singular (notice that we have included (4.8) in the compatibility as
sumptions in( order to obtain this latest expression). This criterion for
i_ .
determining the eigenvalues will be applied in the next section.
: I
4.2 Convergence in One Dimension
'I ]
We now turn to an examination of the convergence of the FAC algorithm
when it is applied to the solution of the finite volume discretization of
problem (3.1). Since the operators A and L arising from this discretization
satisfy the compatibility conditions Ac = Lc, Aci = La, and Aic = Lie
of the previous section, it follows from Theorem 4.1 that convergence of the
algorithm is governed by the relationship (4.3). We are therefore interested
in the eigenvalues of the matrix L/Aj1. Since (4.8) is valid, it is sufficient
it
to set the determinant of the matrix (4.9) equal to zero and then solve
66
i]
for the resulting values of A. Our analysis begins (Section 4.2.1) with the
case that'the: composite grid has a single refinement region. Notice that
in this case two nodes comprise the interface, so the dimensionality of the
eigenvalue problem is equal to two. Therefore, identifying the required
eigenvalues is; simply a matter of the solving a quadratic equation. For
the problem under consideration, the analysis is facilitated by the fact that
i j
this quadratic appears, essentially, in factored form. Later, we extend this
analysis by examining (in Section 4.2.3) the case of multiple refinement
regions.
1 ,i
4.2.1 The Case of a Single Refinement Region
A * 1
We now proceed with an analysis of the spectrum of the operator LiAj
i
corresponding to a finite volume discretization of the problem (3.1). With a
single refinement region on the onedimensional composite grid, and using
either the centered or upwind discretizations of Chapter 3, the operator
L has the structure shown in Figure 4.1. Denote the components of the
composite operators as follows:
0 0
I
Lc
L 0
0 R
La =
0 0
la 0
0 TCI
0 0
Lie
0
0
0 lie 0 0
0 0 tic 0
0
0
67
Lfi
Lif f
flF 0
i, 0 0
Ipl 0 Ipi 0
0 0 0 0
; , Api = 
0 0 0 0
0 Tpi 0 Tpi
0 0 0 tip II lip 0 0 0 ...
I
68
I
I
We assume tljat the Dirichlet value at the left boundary has been moved
to the righthand side of (4.1). Let the respective orders of the matrices L,
J2, Lpt and Ap be Nclt Nc3, NF, and Np. Also, let Nc be one less than
the number of nodes on the global coarse grid (i.e., the order of the global
coarse matrix with the left Dirichlet condition eliminated). The coarse
meshwidth is'denoted by h 1/Nc Denote the ijth element of X~x, the
inverse of the1 generic square invertible matrix X, by *,j. Also, let
denote the j x j submatrix located in lower right corner of X. For example,
if X is ofiorder N, then Xixi = xnn and X^xN = X. With this notation,
we have 1
i
! ;
'' LicLq1 Lci =
IicIciInCinCi 0
0 ricTcih i
(4.10)
and
LipLp Lfi
llFlFllFn llFTFllF1
riFlFilFNlfil rIFrFIlFNpNp
9
AipAp1 Api
liplpia,pxl lipTpiaplffp
riplpiaPffpl rIPTPIapNpNp '
In order to continue, we need to obtain explicit representations for cer
tain elements of the matrices L~x, R~x, LFX, and Ap1. Notice that the
above formulas only involve the corner elements of these matrices. To
[. '
carry out the required analysis for various finite volume discretizations,
we use a; notation that indicates discretizations that change character on
the compjbnents of the coarse composite and true composite grids. We use
a fourtuple, !(si,Z2 : 3,Â£4), with elements xt through Z4 that indicate
the type of discretization used on the respective coarse component of the
coarse composite grid, coarse patch, coarse component of the composite
69
grid, and fine patch. The interface nodes are not specified here. In this
I
section, the type of discretization at these points will always be of the same
type used on the coarse component. (In two dimensions, it will be useful to
allow more flexibility at the interface of the coarse composite grid.) In one
dimension, twb possibilities are considered, upwind (x, = If) and centered
(x, = C) versions of the finite volume discretization. In two dimensions, we
i
will also allow higherorder upwinding to be used on the patch. Generally,
it will be possible to make the notation just described less cumbersome by
employing obvious abbreviations. Initially, we will consider discretizations
that agree in type throughout the two grids.
The first case we consider in detail is the discretization signified by
(If, If : U'; If) = If4, which corresponds to upwinding used throughout the
coarse composite and true composite grids. The analysis will make use of
the generic tridiagonal matrix X = tri[e h, 2e + h, e] of generic order
i
N. Of interest are the corner elements of X1: in, x\jv, zjvi and xjvjv.
1 .
Let X{j be the cofactor of xy. Then [28],
i
Xij = Xji/DetX.
LEMMA 4.2. The four comer elements of X 1 corresponding to the
upwind discretization are given by
n = atjv =
ZJV1 (e + h)N e*
zjf (e + h)N+1 eN+1 e + h
he*"1
Xijv =
0,
and
Zjvi =
(e+hYT+1eir+1
h (e + h)*1 h
(e + h)N+1 eN+1 (e + h)2'
(4.11)
(4.12)
(4.13)
TO
Proof: Since X is tridiagonal, there is a threeterm linear, homogeneous
i
difference; equation for the jth determinant, Zj = DetXjXj
i' !,
Zj = (2e + h)zj_a e(e + h)zj_2, j = 3,4, , N.
The associated characteristic equation is
l
! ' fi2 ft{2e + h) + e(e + h) = 0,
with roots !
i; !
I *
Hi = e + h, /i2 = e.
Hence, itifollows from the theory of difference equations that
! i, Zj = S(e+hy + Tet,
where S and T are unknown constants which are determined uniquely by
ji i
the values of z\ and z2. But
 Z\ 5(c f A) ( Te = 2e f h
; z2 = 5(e + h)2 + Te2 = 3e2 + 3eh + h2.
Thus,
.
and we have
S={e+h)/h, T = e/h
Zj = (e + h)i+1/h e*+1/h.
Now since X is tridiagonal, it follows that the cofactors satisfy
Xn = Xnn = DetXjfixNi
Thus,
*11 = %NN =
Zlf
i (e + hf e"
zn (e + h)N+' eJV+1
71
Also, Ajvi = 1 implies that
fce*"1
X\Tf =
(e f h)N+1 cN+1
Finally, X\n = (e + h)N 1 implies that
xn i =
h (e + h)
N1
(e + h)JV+1ew+1'
More will be said shortly regarding the quality of the approximations given
in the lemma;(see the comments below and also the next section), but be
cause of the list two approximations we will see that L and A are essentially
i ;
lower triangular, so Det(Li AAj) just depends on the diagonal part of
i
the matrix.
i. I
We m!ay let X here assume the role of Ap when the coarse patch has
trivial refinement, i.e., when the global coarse grid plays the role of the
i 1
coarse composite grid in the twolevel version of FAC. We have just found
.If {
the required elements of the inverse of this operator. In this context, it may
I
be useful ,to point out an apparent ambiguity in the above expressions. No
tice that in these expressions, h is always directly related to the dimension,
Nc, of the global coarse grid. However, the role of N in the expressions
i1 ''
is allowed to Vary. For example, when the approximations are made with
regard to. Ap, then N < Np in (4.114.13), and note that, necessarily,
>. I
NP < Nc. Qne sees that there is an interplay between the size of e and
the dimensions of the global coarse and coarse patch grids that enters into
! i
the validity of the approximations. In this particular case, their validity is
commensurate with e being small with respect to h (the reciprocal of Nc)
and the coarse patch being sufficiently large.
Also, X may play the role of the coarse component operators L and R,
72
I
so that
i 7 i
Tii*h'
Now suppose that there are m levels of refinement on the fine patch.
Then
and
Lp = 2m tri[e h/2m, 2e + h/2m, e],
I, n)
DetLp = (2ro)JVi'((e + h/2m)NR+1 e^+1)/(h/2m).
i 1
Using these expressions, we obtain the following corollary to Lemma 4.2.
COROLLARY 4.3. The comer elements of the fine patch matrix Lp
I i
corresponding(to the upwind discretization are given by
and
i Lfii (e + V2)^e^ _ 1
Fn , 2TO((e + h/2m)N*+l cw>+1) ~ 2me + h
; ; j
! l*i#r 22"*((e + h/2m)Nl,+1 eNT+i) *
t
i" ii
h(e + h/2m)K'1 h
Frs r
22m((e + h/V")"'*' Â£"+>) (2Â£ + fc)2
(4.14)
(4.15)
(4.16)
The results of Corollary 4.3 may also be applied to the inverse of Ap, when
an intermediate composite grid is used in place of the global coarse grid,
by performing the replacement m m k (0 < k < m).
In adjlition to these quantities, the pertinent quantities from the FV
discretization1 of Section 3.2.1 are (here (u) denotes upwind, and (c),
centered):
! lie = e h (), e h/2 (c),
lei = e (u), e + h/2 (c),
 ' TW = e ()> c + h/2 (c),
rci = e h (u), e h/2 (c),
73
I
Iif = 2me (u), 2me + h/2 (c),
Zjr/ = 2me h (u), 2me h/2 (c),
tif = 2me h (u), 2me h/2 (c),
tFi = 2me (u), 2TOe + h/2 (c),
hp = 2m~ke (u), 2m~ke + h/2 (c),
Ipi = 2m~ke h (), 2m~ke h/2 (c),
tip = 2m~ke h (), 2m~ke h/2 (c),
_ nmk.
Tpi = 2
(),
2m~k(
+ h/2 (c),
(1 + 2m)e + h 0
0 (1 + 2 m)e + h
Li =
! At =
(1 + 2m)e
Li =
(1 + 2mfe)e + h
0
0
(1 + 2m)e
0
(1 + 2m~k)e + h
(1 + 2m"fc)e
(*),
0
(1 + 2m~k)e
(c).
, M =
Notice that we have allowed for general refinement on the coarse patch by
!
letting 1 ,< k) < m. When k = m this corresponds to letting the global
coarse grid play the role of the coarse composite grid, and when k = 1
to using the grid having one less level of patch refinement than the true
i
composite grid.
Now, using the above values and approximations, we see for the UA case
that Li AAi =
(1 + 2^)e +. h A((l + 2mfc)e + h) 0
0 (1 + 2m)e + h A((l + 2m~k)e + h)
+(A 1)
I
; 1 +A
= Li \Ai + (A 1)LicLqX Lc i LipLp1 Lfi + XAipAp1 Api.
I
A simple calculation then shows that, irrespective of the value of k, the
i A A
solutions of Det(Li AAj) = 0 are Ai = A2 = 1.
' e 0 ' ' 2me 0
e h 2me
Â£ 0
h 2m*e
74
Next, we consider the C4 case where the centered discretization is used
throughout. Again, the corners of the inverse matrix are required. Pro
ceeding as in the previous case, let
; X = tri[e h/2, 2e, e+h/2]
be of order N.
i
LEMMA 4.4. The four comer elements of X~l corresponding to the
centered discretization are given by
(e + fc/2)"(eV2)" _ 1
11 = %NN =
Xin = XN1/ZN =
(e + h/2)*+1(eh/2)J'r+1 e + h/ 2
h (e + h/2)N1
and
sjvi = Xin/zn =
(e + h/2Yr+l(eh/2Yf+1
h.(e + h/ 2)"~'
w 0,
(e + h/2)N+'(eh/2)N+1 (e + h/2)2'
Proof: The jth determinant satisfies
i (
Zj =; 2ezj_x + (e h/2)(e + h/2)zJ_2, j = 3,4, , AT,
with associated roots
i
fix = e + hj2, fi2 = e h/2.
(4.17)
(4.18)
(4.19)
Hence,
; Zj = S(e + h/2y + T(eh/2y,
and in particular
iz\ 5(e ( h/2) f T(e hj2) 2e,
,z2 = S(e + h/2)2 + T(e h/2)2 = 3e2 + h2/4.
1
Solving for S jand T, we find
S = (e + h/2)/h, T = (e h/2)/h.
75
Zj = (e + h/2)i+1/h (e h/2)i+1/h,
Therefore,
and
Also,
and
*11 = %NN =
iiN = Xml zn =
xN1 = Xin/zn =
(e + h/2)N (e h/2)*
{e + h/2)N+'(eh/2)*+1
hje + h^1
(e + /i/2)JV+1(eV2)JV+1
hie + h/if1
{e + h/2)N+'{eh/2)N+''
Again, these results apply to the operator X, and also to Ap if the coarse
patch has trivial refinement. We note that the above argument is incorrect
for one of the matrices, namely iZ, when it corresponds to the centered
I*
stencil for the; righthand portion of the coarse component. What is given
would be! correct for this operator if it had a Dirichlet condition at the
!" 
right boundary. For the Neumann condition, we should use
^2x2
2e e ( hj 2
e h/2 e j h/2
That is, technically, the initial values for the difference equation have to
I i
i
he modified, though we arrive at essentially the same result: now Zj =
(e + h/2)!, so,
rn =
(e + h/2)
Nl
(e + h/2)* e + h/2
i, (
For Lp we have the following corollary to Lemma 4.4.
i
COROLLARY 4.5. The comer elements of the fine patch matrix Lp
'I
correspondingto the centered discretization are given by
 ; ! (e + h/2m+1)Nir (e h/2m+l)Nr _ 1
iFn lFNp,Np j 2m^e + kj2m+1jNp+1 (e h/2m+^Np+i) 2*e+ h/2
(4.20)
76
I
and
h(eh/2r+1Yr'1
22m((e+ V2m+1)JVp+1 (e fc/2m+1)J'rp+1)
(4.21)
Me + h/V**1)"*1
1fw 22^((e + hi2^+i)^f+i (e V2m+1)^+1) ~ (2*e + h/2)2'
, (4.22)
A similar result holds for Ap when the coarse patch has nontrivial refine
ment, again by performing the replacement m * m k. Continuing, then,
\ j
for the (74 case we have Li \Aj =
(1 + 2m)c A(1 + 2m~k)e 0
0 (1 + 2m)e A(1 + 2m k)e
+(A 1)
2me h/2
eh/2 0
0 eh/2
0
2me h/2
+A
)mfc
e h/2
ymk.
h/2
The solutions of Det(Lj AAj) = 0 are again Ai = \2 = 1.
Next we turn to the (UC)2 discretization, which uses upwinding on the
coarse component and centered on the patch for both the coarse and true
i
composite' grids. Using approximations from the last two cases, we have
Li \Ai =
(1 + 2T")e + h A((l + 2m~k)e + h)
0
0
(1 + 2m)e +h A((l + 2m~k)e + h)
I
+(A 1)
e 1 2me 0
n h(2me+h) (2meh/2)(2me+h)
U c L (2mc+h/2) 2me+h/2
+A
2mfce
h(2m~ke+h)
. (2m~kc+h/2)
0
(2mkeh/2)(2m~ke+h)
2 m~ke+h/2 .
77
I'
The eigenvalues are Ai = 1 and
A2
2me + h 2m_"e + hf 2
2me + h/2 2m~ke+h
e (1/2,1] for e G [0,oo).
We consider briefly the validity of the approximations used here. On
the one hand,'we have used (4.114.13) for the coarse component (upwind)
I
operators, which are valid with e small with respect to h. On the other
i'
hand, the (centered) approximations on the patch come from (4.204.22).
Suppose that 'ion the coarse composite grid no refinement is used on the
patch (m1 = 6 and Np Np). Then the approximations will be valid
with e ss h/2 'and also if h continues to decrease. Notice that A2 1 as
r !'
fe > 0. Unfortunately, h approaching zero with e fixed is at odds with
I
the requirement for the coarse component. There are two approaches that
i 1.
avoid thispredicament. The first approach is to use sufficient refinement on
the coarse patch so that h need not be particularly small for (4.204.22) to
be valid (because the exponent Np (or Np) appearing there increases with
m). Also,notice that A2 1 with m large and k small (i.e., with sufficient
refinement on the fine and coarse patches). Unfortunately, taking this
approach rules out the use of the global coarse grid in the role of the coarse
1
composite grid. The second approach allows for this possiblity by altering
the discretization on the coarse patch in order to weaken the requirement
that h be small.
1 ii
This leads  us to consider the U3C discretization, where upwinding is
used throughout the coarse composite grid and on the coarse component of
the composite'lgrid, but centered differencing is used on the fine patch. As
just observed, the following approximations will be valid with e sufficiently
small with respect to h (for the upwind approximation on the coarse com
ponent add on the coarse patch) and with sufficient refinement on the fine
78
patch (for the {centered approximation). Although our motivation here is to
be able to use a coarse patch without refinement (k = m), we have carried
out the analysis allowing arbitrary refinement on this region (0 < k < m).
The approximation is Lj AAj =
(1 + 2m)e + h A((l + 2m~k)e + h) 0
1 0 (1 + 2m)e + h A((l + 2m~k)e + h)
+(A 1) e 0 n * 2me h(2mc+h) 0 (2mcfc/2)(2me+M
U c (2e+/i/2) 2me+h/2
+A
2 m~kt
h 2m~ke
The eigenvalues are Ai = 1 and
, 2 me+h
*2 ~ 2me + h/2
Notice that the desired effect has been achieved in that now A2 approaches
unity while only requiring sufficient refinement on the fine patch. The
!' 1
expressions for the eigenvalues are independent of the refinement used on
I
the coarse patch and so they will be valid, in particular, when the global
coarse grid is used. We also remark that this requirement (of having suffi
cient refinement on the fine patch) is a natural one with respect to the use
i
of the centered type of discretization for advectiondominated problems.
Although ,the expression for A2 is only valid in certain parameter ranges of
the onedimensional problem, we will see that the quantity A2 1 supplies
a fairly accurate approximation to the actual convergence factor even in
i
two dimensions.
Collecting ithe results from our examination of the various cases con
sidered above establishes the following theorem and corollary which apply
1 1
to the onedimensional problem.
79
I
THEOREM 4.6. Suppose that the exact values of the comer elements in
equations (4.114.22) are replaced by there corresponding approximations.
j, I
Then, in the UA case, the C4 case, and in the limit as m oo in the U2C
i! 1
case we have the following.
* A .
Both eigenvalues of I LjAj are zero. Therefore, two steps of the
iteration
I t a a 
r;+1 = (/ L,A?)ri
drives the residual at the interface to zero. This matrix is in fact strictly
lower triangutar: the left and right components of the interface residual
i:
vanish on the first and second steps, respectively, of the iteration.
1
CORQLLARY 4.7. Suppose the twolevel FAC scheme is written in the
form r*+ll'= (J for some composite preconditioner M. Then,
with the assumptions of Theorem 46, this iteration matrix is nilpotent:
(I LMT1)3 = 0, i.e., the twolevel scheme converges in three iterations
I
and, ihus\ may be considered a direct method.
\ \
Proof. The corollary follows from the vanishing of and tf established
1 i
in Theorem 4.1 for i > 1, the resulting expression (4.3) for the updating of
the interface component of the residual, and Theorem 4.6.
It should ibe noted that the results of the last theorem (and its corollary)
are largely independent of the choice of coarse composite grid, i.e., they
I J
hold withlm k levels of refinement on the coarse patch for all values of k}
0 < k < m. In particular, this motivates the use of the the global coarse
grid in this capacity.
4.2.2 Quality of Approximations
The purpose of this section is to assess the numerical quality of the
approximations for elements of inverse matrices given in the previous sec
80
I
tion. We will show that the values given as approximations are valid with
regard to'appropriate floatingpoint computations. However, their validity
is subject to restrictions on the values of the parameters e, h, and m (the
order of patch refinement). It turns out that these restrictions allow for a
large practical range of these parameters.
i ;
We begin with the UA method, concentrating on the inverses of Ap and
l
Lp. Here, the, global coarse grid is used in the role of the coarse composite
grid. Let X denote either Ap or Lp. In the last section, the corner elements
of X~1 were approximated as follows:
, *11 l/(2me + h),
! *mv1/(2 TOe+h),
*uv ~ 0,
I
and , ,
*M h/(2me+ h)2,
where m == 0 for Ap, and m = mr (the additional refinement levels on the
fine patch) for Lp. In Tables 4.14.2 we give, for various values of e, the
four pertinent .values of the inverse matrices, i.e., the 2x2 matrix
I I
, i *11 *1N
*JV1 *JVJV
! :
To obtain the elements of X1, the tridiagonal matrix was reduced to upper
triangular! form by Gaussian elimination (as in [45]), and backsolves were
performed with the canonical basis vectors ei and as righthand sides.
The coarse grid meshwidth is h = 1/64, the coarse patch has 16 nodes,
and there! are1 mr = 4 levels of additional refinement on the fine patch.
Computations'were done in single precision on a Sequent Balance. The
approximations are valid for e sufficiently small with respect to h, i.e., as
e > 0 with h fixed. However, e need not be particularly small, as seen
81
in Table 4.1, for the approximations to be fairly accurate. Furthermore,
comparing the second table with the first, one also sees an improvement in
the approximations with e fixed and h refined. Although h decreases, the
1
effect of this is offset by an increase in the exponent N (the order of the
I !
matrix being inverted) from Np to Np.
Table 4.1
Computed elements of coarse patch matrix inverse
for the upwind case with h = 1/64 and 16 nodes
on the coarse patch.
1 e A? l/(e + h) h/(e + h)2
I .riooi i 6.359301 101 1.757762 1027 6.318860 101 6.359301 101 6.359301 101 6.318860 101
.001 1 6.015038 101 7.627731 10"1B 5.653231 101 6.015038 101 6.015038 101 5.653231 101
i ,01 ], 3.902435 101 1.158410 10"4 2.379538 101 3.902435 101 3.902439 101 2.379536 101
i
Similar results are presented for (UC)2. The approximations from the
previous section are:
I'
and
n ~ l/(2me + h/2),
xNN l/(2me + h/2),
bin ~ 0,
Xffi m h/(2me + h/2)2.
82
Here e .001 is fixed, the patch has length 1/4, and h is allowed to increase.
The expressions (4.174.22) predict the deterioration of the approximations
i
as this happens.
Table 4.2
' Computed elements of fine patch matrix inverse for
the upwind case with h = 1/64, 16 nodes on the
coarse patch, and mr = 4.
i
e Lp1 1/(2 m'e + h) h/(2m'e + h)2
1 .0001 1 ll 5.805515 101 0.000000 10 5.266088101 5.805515 101 i 5.805516 101 5.266252 101
1 1 .001 I 3.162055 101 0.00000010 1.562284 101 3.162055 101 3.162055 101 1.562280 101
i .01 1 1 1 5.693950 10 1.184763 10"10 5.065526 lO'1 5.69394710 5.693950 10 5.065792 10"1
Tables! 4.3t4.4 below confirm this, with the deterioration being more
significant for Ap. This again motivates the use of the upwind scheme on
the coarse1 patch in U3C.
I
83
I Table 4.3
Computed elements of coarse patch matrix inverse
for the centered case with e = 103
h A? l/(e + h/2) h/{e + h/ 2)2
1/512 5.059288 102 0.000000 10 4.999298 102 5.059288 102 5.059288 102 4.999296 102
1/256 I 3.386243 102 7.934145 10"12 4.479160 102 3.386243 102 3.386243 102 4.479158 102
1/128 I 2.041807 102 6.061421 10"1 i 3.247696 102 2.041807 102 2.038217102 3.245568 102
Table 4.4
Computed elements of fine patch matrix inverse for
the centered case with e = 103 and mT = 1.
1 h i Lp1 l/(2e + h/2) h/(2e + h/2)2
1/512 1 1. 3.359580102 0.00000010 2.204448102 3.359580 102 3.359580 102 2.204449 102
1/256 'l 2.529644 102 0.000000 10  2.499649 102 2.529644 102 2.592644 102 2.499648 102
1/128 i ' 1.693122 102 3.808327 10"11 ;! 2.239580 102 1.693122 102 1.693122 102 2.239579 102
I
84
4.2.3 Multiple Refinement Regions
Next we extend the analysis of Section 4.2.1 to the case of multiple
interior refinement regions in one dimension. In general, there may be
various refinement regions, but we suppose that, without loss of generality,
there are but two; the significant difference between the case of one region
and many is the existence of coarse regions between regions of refinement,
not just at the two ends of the grid. Allowing only two refinements will
also allow us to avoid significantly complicating our notation. With this
convention, L has the structure shown in Figure 4.2.
Fig. 4.2. Structure of a onedimensional composite operator
i 1 in the case of multiple refinement regions.
85
The composite grid still has but a single fine component, although it
is now the union of the various refinement regions. This is also true of
the coarse and interface components. Accordingly, the composite operator
still has the same structure, though its main nonzero blocks may now be
further decomposed. For example,
and
Lp
Lr =
JF3
Li2
where each Ljt is a 2 X 2 matrix associated with the ith refinement region.
Notice that a partitioning of Lj with respect to these blocks induces a sim
ilar partitioning of the matrices Ljp and Lpi into block diagonal matrices.
Thus, for example, we may write
Z/F, Lf1 l LfIi
Ljf2 Lp2 LFh
Ljp Lp Lpi =
(The leading diagonal blocks here correspond to the blocks in the above
figure with the t elements.) It follows that Zj LjpLp1 Lpi is a block
diagonal matrix with diagonal blocks L^L^Lp^Lp^. Furthermore, each
one of these blocks has precisely the same form as Zj LipLp1 Lpi had in
the case of a single refinement region, although the dimensions may vary.
Similar remarks hold for the global coarse operator and Ar AjpAp1 Api.
As for the coarsetointerface connections, let
Z
Zc =
M
R
Then LicLcLci
hclcilNCl,NCl 0 0
0 TicTcim ii riclci'rni.No, 0
0 hcrcimNct, i liclcimNC3,Nc3 0
0 0 0 TjC^Clf 11
86
It may be useful to the reader to compare this equation with (4.10). Re
calling the analysis of the values of elements of matrix inverses, we see that
the quantity , which is the upper right corner element of Jlf1, may
be set to 'zero, allowing us to use instead LicL^Lci =
IicIciInCi,n0i 0 0
0 ricrciThn 0 0
0 hcTCI^Na,,! llclci'niNC3,Nc:t 0
. io 0 0 TicTcifli
I 1 A A
With this simplification, the matrix Lj \Ai is block lower triangular.
I % A
Furthermore, each of its diagonal blocks has precisely the form that Zj
\Aj had in the case of a single refinement region. It follows that the analysis
performed for, that case may now be applied to each one of these blocks.
The result is that, for each of the various cases considered previously, the
same values may be given now to the eigenvalues of the matrix I ZjAj ,
i i I
the only change being in their multiplicities according to the actual number
i '}
of regions of refinement.
I, '
4.3 Convergence in Two Dimensions
i '
In thisj section we study the convergence of the twolevel version of FAC
applied tojthe solution of the finite volume discretization of problem (3.2).
Because of the complexity of the structure of the composite grid equations
involving interior refinements, we have not yet extended the analysis of
the previous section to two dimensions. Therefore, the examination here is
i
based on computational evidence rather than on analysis. Guided by the
results of 'the previous section, we attempt to determine to what extent
the results' for J;he onedimensional problem carry over to two dimensions.
Motivated}: by the findings of that section, we concentrate our study on
the behavior of an algorithm that emphasizes the use of the global coarse
grid. The numerical method used here to solve the composite grid equa
j ,,
tions is the twolevel FAC scheme with the global coarse grid as the coarse
composite grid. Subproblems (corresponding to the line patch and global
coarse systems) are solved exactly (i.e., iteratively with a tolerance on
the relative residual of 106). This has been done by performing, in the
appropriate context, the necessary number of relaxations or multigrid V
cycles (for further details see Section 5.3.1).
The first result emphasizes the effect that the initial guess has on the
j ji
convergence behavior. Consider the solution of the upwind/upwind dis
. i i . ,
cretization of Problem 3.1 of Section 3.3.1 by FAC using the upwind/upwind
discretization on the global coarse grid (i.e., the Z74 technique). The prob
lem parameters here correspond to those used to generate Tables 3.13.2.
Table 4.5 presents the sequences of relative composite residuals correspond
ing to two solutions of the discrete problem. The two solves differ in the
choice of the initial composite guess, z. The zero vector and a random
I
vector with elements between zero and one were used. The results appear,
respectively, in the left and middle columns of the table. Recall that ]r* 112
here actually measures the current interface component of the residual,
! 1
since r^., = 0 for i > 1 (see Section 4.1). Notice that there is a sig
['
nificant difference in the rates at which these interface residuals approach
t !:
zero in the two solves. This difference is due to the convergence behavior
of the resikuals at the nodes that correspond to the tangential (north and
south) boundaries of the interface. At the in and outflow (west and east,
I !
respectively) boundaries, the convergence behavior is like that predicted
by the analysis of the onedimensional problem in both solves the resid
uals at the interface converge rapidly to zero. However, at the tangential
Table 4.5
Relative residuals as functions of FAC iterations
' using the U4 method in two dimensions.
Upurind / Upwind
IIA/lli" IU
;i i z = 0 z Random Modified Stencil z Random
1 1.825 10"4 6.406.10"2 5.592 10"2
2 3.355 108 1.411 102 1.108 103
3 6.803 10"7 6.208 103 9.313 105
4 2.742 10"3 8.291 106
5: 1.212 10~3 6.013 10"7
6 5.358 10"4
7, 2.370 10"4
8 1.049 10~4
9, 4.639 10B
ip 2.053 10~6
boundaries the convergence rate is generally much slower, as in the middle
column of Table 4.5. The reason that this effect does not appear initially
I (
in the left jcolumn is that, with the initial guess of zero, very accurate solu
'i
tion values are produced at the tangential boundaries as a result of the first
global coarse solve. Thus, the residuals at these nodes are extremely small
1 i
at the end'.of the first FAC iteration and, therefore, throughout the ensuing
i
iterations.! We have developed a remedy for the problem of slow conver
i i,
gence at the tangential boundaries that involves modifying the definition
of the gloial coarse stencil at these nodes. The right column of Table 4.5
. [i
gives the results for a solve using this modification. The actual changes to
the stenciljare obtained by requiring that interfacetointerface connections
I ;
for the glbbal coarse operator agree with those of the composite operator,
i.e., A\ =}Li at nodes lying on the tangential boundaries. As a practical
89
matter, these1 changes do not adversely affect the efficiency of the algo
rithms implementation. On the contrary, they can be beneficial, making
the global coarse matrix more diagonally dominant, which facilitates the
solution of these equations by an iterative method.
To motivate these modifications, consider the case where the only in
terface component is a tangential one, i.e., the twodimensional composite
grid associated with the domain partitioned into two horizontal strips, as
depicted in Figure 4.3. Also, let e > 0.
O O 0 0 O O O O 0
ooooooooo
I I
O O O 0 O 0 O 0 0
6)o0o0oo0o0oo0o
i
Fig. 4.3. Composite grid with coarse (o), fine (),
; interface (0) and slave (o) nodes.
Then, using the stencils of Section 3.2.3, one sees that only the diagonal
blocks of the composite operator X in (4.1) are nonzero. The same is true
for A. The Schur complements in these respective operators then become
Li and Ai, and the equation for the updating of the interface residual
simplifies to
rj = (/ L,Aj')r).
I
This clarifies the reason for the modification. Notice that without the
90
I
modification, Ai =
\Lh where a =  + . Thus,
r/+1 = (1 *) '/>
giving convergence factors between one quarter and one half. This agrees
roughly with the observed convergence rates. We remark that as e 1
the modification becomes less beneficial, and even detrimental, but the way
that this happens depends on both e and the amount of patch refinement.
Tables14.6;4.7 show convergence rates as functions of e and fine patch
refinement for the U4 method applied to the upwind/upwind discretization
of Problem 3.1. The coarse component mesh width is he = 1/32, and m
gives the number of additional levels of refinement on the fine patch. We
Table 4.6
FAC convergence as a function of e
for the Z74 method.
l Upwind/Upwind m = 1
1 i ! Mli/ Ill*
1 1 z!' Unmodified ' Stencil e = 10"2 Modified Stencil e = lO"3 Modified Stencil e = 10"4 Modified Stencil e = 10"5
1,: 5.964 10"2 6.824.10"2 8.160 lO"2 8.403 10"2
2 4.287 10"3 5.606 10"4 1.654 lO"4 2.071 10"B
3r 4:554 10"4 1.737.10"B 4.138 10"7 1.790 10"7
4;, 51855 10B 8.934 lO"7
5! 8.080 10"6
6ii Is 143 lO'6
7, 2.269 10"7
8
9,
10
91
used a random initial guess throughout in these tests and, in some cases,
the modification to the global coarse stencil mentioned above. We found
that the advantage gained by using this modification is lost as e increases,
and use of the unmodified stencil can give better performance. For a
given e, the results appearing in Tables 4.64.7 correspond to whichever
Table 4.7
FAC convergence as a function of e
for the U4 method.
Upwind{Upwind m = 4
1. Ik*2/ II ^
1 i Unmodified Stencil e = IO"2 Unmodified Stencil e = 10~3 Modified Stencil e = IO4 Modified Stencil e = 106
1 11691 101 7.516 102 5.592 102 5.738 10"2
2 1.232 10~2 1.506 102 1.108 IO"3 4.108 10~6
3! 1.135 IO"3 3.523 103 9.313 IO"6 4.659 10"7
4 1.192 IO"4 1.171 10"3 8.291 IO"6
5 1.259 10~6 4.012 104 6.013 10"*
6i 1:321 IO"6 1.384 10"4
7. 2.078 107 4.786 10"5
8' ! 1.657 10~5
s; : 5.769 10"6
10; 2.014 10~6
strategy showed better performance. Thus, the leftmost column in Table
4.6 and the two leftmost columns of Table 4.7, correspond to the use of
the unmodified stencil, while the other results are for the modified stencil.
i '
I 
Notice that the convergence rates generally improve as e decreases, which is
what we would expect as a result of the analysis of Section 4.2. The analysis
also predicts indepedence of the rates on the level of patch refinement. The
i 92
I.
results here show a mild degradation as the level of refinement is increased.
One possible way to offset this effect (i.e., to make the need for refinement
1. 1
less critical) is to use a more accurate upwind scheme on the fine patch.
Tables 4.8,4.9 show results of a similar set of experiments with stan
dard upwinding replaced by higherorder upwinding on the fine patch, i.e.,
the (UH)2 technique (here, H denotes the use of the higherorder up
wind discretization on a composite grid component) is used to solve the
upwind/Hup composite grid problem. This FAC strategy was not analyzed
Table 4.8
j FAC convergence as a function of e
' for the (UH)2 method.
i
1 Upwind/Hup m = 1
!' 1 t' Unmodified Stencil e = 10"2 h%/ Modified Stencil e = 103 k2 Modified Stencil e = 104 Modified Stencil e = lO"5
1 6.016 10"2 7.415 10"2 9.465 10"2 9.817 10"2
i 4.376 10"3 1.081 10"3 9.657 lO"5 1.051 10"B
3 4.271 104 3.167 10B 2.522lO"7 1.625 10"7
4 5:126 105 9.71110"7
5 6.654 1O0
6 9,051 10"7
7: '
8
9;
10 1
in the previous section, but the experiments here indicate that its con
I;
vergence behavior strongly resembles that of Z74. This is evidenced by
comparing Tables 4.84.9 with Tables 4.64.7, respectively (the same strat
93
egy as above, of modifying the interface stencil, has been used here). When
I ,
combined with a comparison of the accuracy of the two upwind discretiza
tions (see Section 3.3.1), the results indicate that there is little advantage to
! I
using the: upWind/upwind discretization. The more accurate upwind/Hup
method can be used without sacrificing efficiency. In this study, we used
strategies based on block Jacobi relaxation to solve the patch systems.
With the [upwind discretization, this means performing backsubstitutions
associated with the inversion of tridiagonal blocks at each relaxation step.
With the Jmore accurate method, the increase in work corresponds to using
blocks having one additional diagonal. These relaxationbased strategies
1 '
will be discussed more fully in Chapter 5.
Table 4.9
FAC convergence as a function of e
for the (UH)2 method.
\t Upwind/Hup m = 4
\ h 1 i I 1 Unmodified Stencil fc = 102 lr2/ Unmodified Stencil e = lO"3 I2 Modified Stencil e = 104 Modified Stencil e = 10B
If: l.;674 lO"1 7.082 lO"2 4.933 102 5.026 lO"2
2! 1:389 lO"2 1.352 lO"2 1.158 103 7.143 lO"6
3]' 1J286 lO"3 3.247 103 8.733 10"B 3.364 lO"7
4 i 1.368 104 1.086 103 7.515 106
5!" l.'!487 lO"5 1 O iH 00 CO 5.304 lO"7
6 i, 1.643 106 1.281 104
7!', 2.377 107 4.424 10B
8 fj j 1.531 lO5
9 ; , 5.326 10B
10 : 1.852 lO"8
Table 4.10 displays the results for the CA technique applied to the solu
I'
! 94
i
i

PAGE 1
MULTILEVEL METHODS FOR THE SOL UTI ON OF ADVECTIONDOMINATED ELLIPTIC PROBLEMS ON COMPOSITE GRIDS by James Scott Otto B.S., Georgia. Sta.te University, 1983 M.A., University of New Mexico, 1985 A thesis submitted to the Faculty of the Gra.dua.te School of the University of Colorado a.t Denver in pa.rtia.l fulfillment of the requirements for the degree of Doctor of Philosophy Applied Ma.thema.tics 1992
PAGE 2
This thesis for the Doctor of Philosophy degree by James Scott Otto has been approved for the Department of Mathematics Thomas F. Russell I '1z;nr.
PAGE 3
Otto, James Scott {Ph.D., Applied Mathematics) Multilevel Methods for the Solution of AdvectionDominated Elliptic on Composite Grids Thesis directed by Professor Thomas A. Manteuffel ABSTRACT We analyze the use of multilevel algorithms in the solution of advection dominateQ. advectiondiffusion equations in one and two dimensions. In particular, we study the behavior of the FAC (Fast Adaptive Composite) scheme as applied to a composite grid discretization of the problem, where the discretization type is allowed to vary {from upwind to centered dif ferences, for example) on the components of the grid. In one dimension, the analysis leads to the interpretation of a twolevel version of FAC as a direct This analysis also provides insight into the behavior of the algorithm, and guidance for its implementation, in two dimensions. In two dimensions, we suppose the problem to have been transformed, by an orthogonal coordinate mapping, in such a way that its character istics are aligned with a Cartesian product grid. Various discretization strategies for this problem are studied, with emphasis placed on a finite volume Upwind, centered and higherorder upwind types of finite volume discretization are considered. This method is attractive in that its matrix stencils capture all of the advection of the discrete problem in their main diagonal blocks. As a result, the block Jacobi relaxation is an excel lent smoother for multigrid as applied to problems on uniform subgrids. This allows us to formulate highly effective multilevel algorithms for the Ill
PAGE 4
global equations that preserve the inherent parallel nature of the strongly advectiondominated problem. The ultimate numerical method developed uses twolevel FAC as an algebraic solver in a nested way to obtrun an optimally efficient full multigridlike algorithm. This accurately represents the content of the candidate's thesis. I recommend its publication. Signed lV
PAGE 5
ACKNOWLEGEMENTS I wish to thank, first and foremost, my advisor, Tom Manteuffel, for his support, academic and financial, that has allowed me to take part in a variety of interesting areas of research which he has introduced me to, including the one of this dissertation. Secondly, I would like to thank Steve McCormick, who has been a consultant to this work since its beginning, and whose expertise in multilevel algorithms has proved invaluable. I would also like to thank my other committee members, Jan Mandel, Steve Pruess, and Tom Russell, for their participation. Particular thanks goes to Debbie Wangerinforher support, friendship, and sense of humor. Thanks, also, to Bill Briggs and Roland Sweet for the initial support I received as a research assistant in Denver. Finally, I wish to thank the people who nurtured my interest in mathematics, those in the Department of Mathematics at Georgia State University, and in particular, of course, George Davis who inspired p1e to study numerical linear algebra. v
PAGE 6
CONTENTS Chapt,er 1. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Notation . . . . . . . . . . . . . . . . . . . . . . . . . 6 2. Problem Description and Overview of Solution Methods ........ 7 3. Discretization ................................................ 16 3.1 Description of Composite Grids ........................... 17 3.2 Discretization Methods ................................... 20 3.2.1 Finite Volume Discretization in One Dimension ...... 23 3.2.2 Discretization by Element Methods .................. 27 3.2.3 Finite Volume Discretization in Two Dimensions ..... 36 :3.2'.4 A Comparison of Discretization Methods ............ 41 3.3 Accuracy of the Discretizations ........................... 43 3.3.1 Sample Problems in Two Dimensions ................ 49 4. Convergence of FAC .......................................... 59 4.1 The TwoLevel Method ................................... 60 4.2 Convergence in One Dimension .. .. .. .. .. . .. .. .. .. . .. 66 4.2.1 The Case of a Single Refinement Region ............. 67 4.2.2 Quality of Approximations . . . . . . . . . . . . . 80 4.2.3 Multiple Refinement Regions ........................ 85 4.3 Convergence in Two Dimensions .......................... 87 5. Alternative FAC Algorithms .................................. 97 5.1 FAC as Preconditioning .................................. 97 5.2 A Comparison of FAC and Domain Decomposition ........ 99 VI
PAGE 7
5.3 Modifications Based on Multigrid ........................ 103 5.3.1 The Significance of Line Relaxation ................. 104 5.3.2 FAC and Nested Iteration .......................... 114 6. Summary of Results . . . . . . . . . . . . . . . . . . . . 125 References . . . . . . . . . . . . . . . . . . . . . . . . . . 132 V11
PAGE 8
1. Introduction The topic of this thesis is the efficient solution, by iterative methods, of advectiondominated advectiondiffusion equations. Although we con sider only elliptic second order equations, the algorithms presented in this study may be used to solve the problems that arise at each time step in implicit methods for transient equations with nearly hyperbolic charac ter. We consider solving equations whose true solutions are such that the use of local grid refinement in their discretization is desirable. The type of numerical methods we consider for solving the discretized problem are based on multilevel techniques; the ultimate numerical solution method that results from this study is a full multigridlike algorithm that preserves the inherent parallel nature of the strongly advectiondominated elliptic partial differential equation. I Multilevel, or multigrid, techniques have proven to be efficient methods for the solution of elliptic partial differential equations. However, to a large extent, the of classical multigrid algorithms depend on the degree of ellipticity of the problem along with the type of discretization being employed. These two factors go handinhand, and what is more accurate to say is that the success of the method depends on the degree of ellipticity of the discretized version of the partial differential equation. Much of the content of.this dissertation is devoted to studying the effects that ellipticity and discretization type have on the efficiency of conventional multigrid methods in the context of an advectiondominated problem discretized on a composite grid. Although we study a one dimensional problem, largely as a means of gaining insight, particular emphasis is placed on solving the 1
PAGE 9
twodimensional advectiondiffusion equation in a special form. Following: Chin and Manteuffel [14), in two dimensions we assume a coordinate mapping has been used to transform the problem to one where the How velocity is aligned with gridlines in a Cartesian grid. As shown in [14), with the problem in this form, line relaxation methods provide robust solvers for the discrete problem. One way this work extends that in [14) is by allowing multigrid, based on line Jacobi relaxation, to play the role of the iterative solver. We also include the use of local refinement in the discretization. Our choice of relaxation strategy allows us to retain the parallelism of the strongly advectiondominated problem in efficient solvers for both coarse and fine composite grid subproblems. The particular form of multigrid we study is the FAC (Fast Adaptive Composite) method introduced by McCormick [32) to taylor multigrid to the structure of composite grid problems. We analyze the behavior of a two, level version of this method as applied to the advectiondiffusion problem in one and two dimensions. A detailed analysis of FAC convergence in one dimension yields an interpretation of the scheme as a direct method. Perhaps more importantly, it provides insight in the appropriate choice of grid and type to be used on the coarse level. This information is used in determining appropriate strategies for the use of FAC as applied to the problem in two dimensions, whether it is used in its t'!olevel form or as a full m')l.ltigridlike algorthm. A description of the contents of this dissertation follows. Chapter 2 provides an overview ofthe problem to be solved and the solution methods used, and relevant related work. Chapter 3 introduces the com posite grid and presents a summary of discretization methods for the one and twodimensional problems. The structure and notation for the com2
PAGE 10
posite grid used in this study are given in Section 3.1. Section 3.2 describes three discretization methods. A finite volume discretization is emphasized. For the most part, our description of the finite volume discretization is standard, however, we show how to adapt this method to the structure of the composite grid, borrowing a technique from finite element and finite volume element discretization&. These two discretizations are introduced in Section 3.2.2 as completely systematic methods and also as mearis of motivating the treatment of the interface portion of the finite volume dis cretization. Although the former methods are more systematic, the finite volume discretization offers advantages in terms of simplicity and ease of use. Section 3.2.1 describes the finite volume discretization for the one dimensional version of our problem. This discretization serves as the basis for the analysis of FAC convergence in one dimension later in Chapter 4. We no'tice that if the finite element subspace is chosen appropriately this method yields the same discrete problem as the finite element and finite volume element methods. The structure of a composite basis for this space is introduced. In Section 3.2.2 the generalization of this basis to two dimensions is made. Use of this basis allows the simple handling of the nonuniformity at the interface when discretizing by the element methods. In Section 3.2.3 the twodimensional finite volume discretization is intro duced. This method serves as the basis for later studies of FAC convergence in two dimensions. It is based on the socalled control volume approach described in [36]. Section compares briefly the relative merits of the various dis cretization& treated in this study, and also identifies appropriate grid trans fer strategies. 3
PAGE 11
In Section 3.3 we discuss the accuracy of the discretizations in two dimensions. We cite some recent results on the accuracy of finite volume and finite volume element methods for diffusion problems, analyze the accuracy of the finite element method for the advectiondiffusion problem and provide examples demonstrating the accuracy of the particular finite volume method used here. Chapter 4 is a study of convergence of the twolevel FAC scheme ap plied to the volume discretization of the oneand twodimensional advectiondiff.usion problems. In Section 4.1 it is shown that, with a certain compatibility assumption between the coarse and fine composite grid operators, this method has the form of an interface preconditioning, i.e., the interface component of the residual satisfies where, Lr is the Schur complement in the composite grid operator, and Ar is a preconditioner. It turns out that, for FAC; the latter is the Schur complement in the coarse composite grid operator. Section 4.2 studies convergence of the twolevel scheme in one dimen sion, and contains the main theoretical results of this study. Our conver gence analysis' proceeds by an investigation ofthe spectrum of Br = LrA.C1 In Section.4.2.1 this spectrum is analyzed in the case of a single refinement region, allowing for the used of mixed discretization types on the composite grid. We also take into consideration the effect of the type of discretization used on the coarse grid. The analysis shows that good convergence results can be obtained with significant jumps in refinement between the compo nents of the grid, and, in particular, motivates the use of a global coarse grid on the coarse level. In Section 4.2.3 these results are generalized to 4
PAGE 12
the case of multiple refinement regions. Section 4.3 is a study of the twolevel FAC scheme applied in two dimen sions. Taking a cue from the onedimensional analysis, the global coarse grid is used at the coarse level in specific FAC implementations for problems that involve upwind, centered, and higherorder upwind descretizations on the region of refinement. Generally, the results of this section indicate that the twolevel scheme provides a good algebraic solver for the advection dominated problem. However, better efficiency can be obtained by using this scheme in a more sophisticated way. In Chapter 5 we. consider two ways of using FAC more efficiently. Sec tion 5.1 views the twolevel scheme as a preconditioned iterative method, and suggests using polynomial acceleration to improve its convergence. This is the method of Schur complement domain decomposition and corre sponds to the use of FAC by domain decomposition practitioners. SeCtion 5.2 is a brief comparison of twolevel FAC with Schur complement domain decomposition. Rather than use the preconditioning perspective, that suggests solving the composite grid equations algebraically, in Section 5.3 we take the more efficient approach of using full multigrid to solve the equations to the level of accuracy of the discretization. An important component of our implementation of this approach is the use of strategies based on block Jacobi relaxation for solving problems on the subregions. Section 5.3.1 contains an of this type of relaxation (as relaxation, and as a multigrid smoother) for. upwind and centered discretizations of the particular prob lem we study. In Section 5.3.2 nested iteration and the twolevel FAC scheme are combined to o.btain a full multigrid version of FAC. Computational results 5
PAGE 13
presented in that section show that, for the various discretization mixtures, only one or two FAC iterations are required on each of the composite levels to obtain an '!pproximation by this version of FAC with accuracy compara ble to the exact solution of the composite grid equations. Furthermore, we show that this is the case with complete solves of problems on the subre gions replaced by inexpensive approximations. This is significant because practitioners of the twolevel approach have noted dramatic degradation in global convergence when subproblems on the refinement region are solved The method presented here, however, yields accurate solu tions with optimal efficiency on this region, i.e., the total amount of work there comparable to that of a full multigrid iteration. Chapter 6 is a summary of the results of this dissertation. It also provides recotnmendations for their practical utilization. 1.1 Notation For the most part, the following notational conventions are used in this dissertation. Particular usage and exceptions should be clear from context. A,L,X,: :Vi; Zi; A, L, ... z, r, .. r, ... a, J, a, c/J, ... CJ) lJ1 ' Hie DetX Xi; X;x; Matrices or differential operators. The ijth entry of the matrix X. The ijth entry of x1 Composite matrices. Vectors. Composite vectors. Scalars or functions. Finite dimensional function spaces. Sobolev space of order k. The determinant of the square matrix X. The cofactor associated with :Vi; The j X j submatrix located in the lower right corner of the square matrix X. 6
PAGE 14
2. Problem Description and Overview of Solution Methods This chapter is devoted to describing the problem to be solved, and to providing an overview of the solution methods and analysis used in this study. Although the theme of multigrid applied to advectiondominated prob lems is central to this study, there are a variety of ways of viewing this work in which the significance of multigrid varies. One way that is essential is to view the study as an extension of the work in [14] by Chin and Man teuffel. There, the authors consider the solution of the twodimensional advectiondiffusion equation Eilu(z,y) + aui(z,y) + bug(z,y) + du(z,y) = g(z,y), (2.1) where E is a !!mall positive constant. Notice that if E is equal to zero in (2.1 ), then the equation is pure advection and is easily solved, in prin ciple, by integrating along its characteristics. In practice, however, the may be difficult to track and, in any case, we are interested in solving the problem with E small but positive. Nevertheless, when E is sufficiently l'mall one feels that advantage may still be gained from the problem's nearly onedimensional character. In [14] the authors do this by using an orthogonal transformation that is induced by the characteristics. Under appropriate assumptions, the latter transforms (2.1) into the equation ELlu(:z:,y) + Uz(z,y) + cu(:z:,y) = f(:z:,y). (2.2) With this forrp., the characteristics are simply y = constant and any pro cedure that employs solving along the characteristic directions will be fa7
PAGE 15
cilitated by having the characteristics aligned with gridlines in a Cartesian product grid. In [14] the authors analyze the spectrum of the block SOR iteration matrix associated with a centered difference discretization of (2.2) (with c = 0). This iteration is based on alternative matrix splittings that use either the zor yderivative terms of the discretization as the principal part of the splitting. The relative size of the diffusion coefficient is mea sured by the ratio y = h/(2e), where his the meshsize. When the problem is advection dominated (i.e., for y 1), and with relaxation performed along lines in the they obtain a bound on the spectral radius of the iteration matrix that is less than 1/5. We also consider the numerical solution of (2.2). One way that this work differs from (14] is that multigrid plays the role of the iterative solver. The algorithms used in this study also retain the essential parallelizability of the problem in this form by using block Jacobi, or line relaxation, as the smoother on multigrid is based. Another way we extend the work in [14] is by including local refinement in the discretization. The use of local refinement may be called for due to various anomalies of the problem. These may include abrupt changes in the side [6] or the existence of boundary layers [40], [37]. Also, we emphasize the use of upwind rather than centered differences, although we consider employing various discretization types on the sub regions of a composite grid; use of the composite grid structure makes it particularly easy to interface different discretization& on subregions. It also allows one to interface appropriate solution techniques for these subprob lems. The propitous aspects of concentrating effort locally in the context of a composite problem are summarized succintly in (44, p.2]: This apHroach offers a number of advantages. An existing code 8
PAGE 16
can be upgraded, in order to increase the accuracy locally, without a radical redesign of the data structures, etc., since we can use the old code to solve one or several standard problems. Is sues of .data structures and geometry are generally simpler if we design programs for composite mesh problems in terms of simpler standard problems. The use of simple standard prob lems also tends to improve the performance of the programs on vector machines. Finally, we note that if each of the standard problerq.s has appreciatively fewer degrees of freedom than the composite model, then we might benefit from solving a number of smaller problems rather than one large one. In other words, we might view our approach as a divideandconquer strategy. The composite structure also plays a fundamental role in the convenient economization of multigrid to problems requiring local refinement. Multi grid in this context is called the FAC (Fast Adaptive Composite) method, which was introduced by McCormick in (32] and is treated in depth in (33]. It is the particular form of multigrid studied here. The economization of the method lies in its avoidance of unneccessary computations with respect to the coarse component of the grid. Its ease of use derives from the fact that most of the computations involve solving familiar problems on uniform subgrids. Local uniform refinement has been abandoned by some multigrid re searchers [2], (34], (38] in order to obtain more :flexibility in using refine ment. In particular, the methods of (34] and (38] allow for quite arbitrary refinements, of which local uniform refinement can be seen as a special case. When arbitrary refinement is allowed, however, the actual refinement typically is not specified before computations begin, but rather is determined dynamically as computations proceed. Although such an ap proach may be attractive because of its high degree of generality, it may not be appropriate in certain contexts. For example, in many instances the user has a good idea of the location and the degree of refinement needed for a particular problem, so the effort expended in determining the refine9
PAGE 17
ment would be wasted. Because of the dynamic nature of the corresponding algorithms, they are difficult to parallelize efficiently. Also, because of the lack of a predetermined grid structure, analysis of these methods is difficult. Finally, the use of arbitrary refinement may actually limit the capability to represent a solution accurately in the sense that such an approach necessitates the use of a single discretization method, whereas it be desirable to use methods of different orders in different subregions of the domain, for example. The use of local uniform refinement does not have these drawbacks. Also, the use of a composite grid allows for the use of several levels of immediate refinement on a given subregion. This can be important if the solution changes rapidly. We emphasize such behavior in this study and investigate the effects that using multiple levels of refinement ;has on FAC convergence. This is an important aspect of this study, since FAC (albeit under different names) is currently being used in this way to solve practical problems. For example, in (6] the use of the closely related BEPS method in oil resevoir simulation is discussed. Mul tiple refinement levels are used there to accomodate the representation of injection wells by delta functions. We note that the algorithm appearing there, although equivalent to the twolevel version ofFAC, can be described without reference to multigrid. Indeed, the authors view the method as a preconditioned iterative method. Much of the content of this dissertation is devoted to analysis of the convergence of the twolevel FAC scheme, so one further interpretation of this work is as an examination of the effectiveness of a precondi tioner. A composite grid (to be described more throughly in the next chapter) consists of.three components: an interface serves as the boundary between the two true subregions of the grid, the coarse and fine patch components. 10
PAGE 18
The latter component corresponds to a simply connected subregion of the I domain local refinement is deemed necessary. We consider three discretization methods in the context of the composte grid: finite element, finite volume element and finite volume, emphasizing the latter. For the most part, our description of the finite volume discretization is quite stan dard. we show how to adapt this method to the structure of the composite grid, introducing a novel way of computing patchtointerface connections. This technique involves interpolation of the interface (with the mesh size of the coarse component) to the boundary of the fine patch, permitting the use of a standard stencil throughout the region of the fine patch. The technique is borrowed from finite element and finite volume element com posite grid discretization&, where it occurs naturally. Although these two methods are attractive due to being highly systematic, the finite volume discretization offers advantages in terms of simplicity and ease of use. The matrix stencils obtained with this method are, on uniform grids, familiar fivepoint ones that correspond to either firstorder upwind or secondorder centered differences. An advantage of using this method is that, with the advectiondiffusion problem in the form (2.2), the advection terms are iso lated in the main tridiagonal portion ofthe (appropriately ordered) matrix stencil (this is. not true for the elementtype discretization methods). The significance of this is that, then, line Jacobi provides a robust relaxation method and multigrid smoother that preserves the inherent paralellism possessed by the problem as e + 0. This is an important feature of our implementation of the FAC algorithm, however it remains invisible until well towards the end (Chapter 5) of this study. The reason for this is that our analysis focuses on a twolevel version of FAC that proceeds by solving 11
PAGE 19
subproblems on the components of the composite grid, and it is assumed that these problems are solved exactly. We study convergence of this twolevel version of FAC applied to the finite volume discretization of (2.2) and its onedimensional counterpart. We show that when a natural coarse component compatibility between the coarse and fine level composite operators is satisfied, the convergence of the FAC scheme is governed by convergence of the interface component the residual. The iterates for this component can be shown to satisfy relationship Here, L1 is the Schur complement in the composite grid operator, and A1 is the Schur complement corresponding to a coarselevel discretization. This is a distinguishing feature of the FAC method when viewed as an interface preconditioning: it preconditions the Schur complement in the composite operator with the Schur complement in a coarse grid operator. Multigrid convergence theory for various discretizations on uniform grids of the onedimensional advectiondiffusion equation has been devel oped in [24] in [3] using Fourier analysis. In [24] the author considers the twolevel multigrid iteration, based on GaussSeidel smoothing, applied to upwind and centered finite difference discretizations. For the upwind case, a bound of 1/3 on the spectral radius of the iteration matrix is ob tained, independent of 'Y. For the centered case, good convergence rates are obtained by developing a strategy for introducing artificial diffusion as the grid levels vary. Corresponding results are obtained in [3], where a similar approach is taken with respect to a streamline diffusion finite element discretization. Convergence of the twolevel scheme is proved for 12
PAGE 20
sufficiently artificial diffusion or a sufficiently fine grid. The results of our analysis of the composite grid problem are similar in character to these results for uniform grids. Our analysis of the twolevel FAC algorithm for the onedimensional problem discretized on a composite grid proceeds by an investigation of the spectrum of B1 = L1A.I1 We determine this spectrum in the case of a single refinement region, and also extend the analysis to the case of multiple refinement regions. This analysis covers various mixtures of upwind and centered type discretization& on the components of the grids. With is used on both levels, we show that the spectral radius of the iteration matrix IB1 is zero in suitable parameter ranges. This result applies ,to rather arbitrary choices of the grid used at the coarse level, including a global coarse grid (i.e., the uniform grid over the whole domain with the meshsize of the coarse component). For a discretization that uses upwind and centered differences, respectively, on the coarse and fine patch components of the composite grid, we show that with upwinding used on the coarse level, the spectral radius of the iteration approaches zero as the resolution is increased on the fine patch. In our study of the twolevel FAC scheme applied in two dimensions, the convergence behavior is seen to conform generally to the analysis for the onedimensional problem. This is particularly true for a method that employs centered differencing on the fine patch. It is shown that this method works 'well with sufficient refinement on the patch. Also, the results of the onedimensional analysis may be used to develop a critereon for switching to tllls method when the twolevel algorithm is used in a nested iteration. However, when standard upwinding is used on the fine patch, it is 13
PAGE 21
. discovered that the convergence rates deteriorate at the tangential bound aries of the :patch. The reason for this degradation is demonstrated, and a successful remedy that involves modification of the coarse grid stencil is implemented. Similar results are obtained for higherorder upwinding. These results indicate that the use of higherorder upwinding can replace standard upwinding on the fine patch with little sacrifice in efficiency. Although we obtain satisfactory convergence rates for the twolevel FAC scheme, we feel that significant improvement can be made in terms of optimizing the efficiency of the method. For example, the analysis and experiments that we perform with respect to this scheme employ exact solution of the subgrid problems in its definition, and these solves must be performed during each of the scheme's global iterations. An obvious and attractive alternative is to replace these exact solves with inexpensive solution approximations obtained by using a small number of steps of some iterative method. However, in practice it has been noted (in the solution of diffusiontype problems, for instance) when exact solutions of problems on the fine patch are replaced by partial ones that significant degradation in convergence rates ofthe global process occur (see [23], for example). We I have noticed similar effects when applying this strategy to the advectiondominated problem here. Our solution to this problem, which allows us to employ partial solution on subgrids, is to use FAC in a nested way, in the same way that the Vcycle is used in the optimally efficient full multigrid method. Thus, two level FAC plays the role of a weak algebraic solver on each of the composite grids in the succesion of grids lying between the global coarse grid and the true composite grid. This yields an algorithm that usually requires only a few relaxations on the coarse mesh underlying the composite grid, and 14
PAGE 22
I i which has optimal efficiency on the fine patch. I We close section by noting that, although we have introduced our numerical methods in the context ofthe elliptic problem (2.1), they are also applicable to.the solution ofthe mixed hyperbolicparabolic type equation Ut Eau + a.uz, + bu;; + du = g, (2.3) when E is sm:all, i.e., when this problem is mainly hyperbolic in nature. Although the. results contained in what follows pertain to the elliptic prob., lem, they well to the solution of {2;3) under the abovementioned I transformation of spatial coordinates when implicit methods such as backI ward Euler o'! CrankNicolson are used to perform the discretization in I time. The resulting linear systems at each time step then closely resemble the ones described herein except for having improved stability properties (that is, matrix diagonal dominance) due to the introduction of a positive capacitance term. 15
PAGE 23
3. Discretization This .chapter is devoted to studying discretization methods for elliptic problems on composite grids. We begin in Section 3.1 with a description of the composite grids to be used in this study. In Section 3.2 we describe finite element, finite volume element, and finite volume composite grid discretizations. Although the results of this study pertain for the most part to the finite volume discretization, it is useful to consider the subspace approach of the element methods. This ap proach permits the systematic handling of the discretization on a composite grid in more than one dimension, and also lends guidance in performing the discretization in the absence of a subspace. In one dimension, it turns out that if the.finite element subspace is chosen appropriately the finite volume inethod yields the same discrete problem as the finite element and finite volume element methods. We introduce a composite basis for this subspace in Section 3.2.1. The generalization of this basis to two dimensions is made in Section 3.2.2, where we show how its use facilitates performing the el ement discretizations. In Section 3.2.3 the finite volume discretization in two dimensions is introduced. Based on the control volume approach de scribed in [36], the method we describe in the context of the composite grid is similar t,o the one in [22], where a finite volume method is applied to a diffusion problem on a cellcentered composite grid (our grid differs from this somewhat. in that the interface nodes are not cellcentered). Our ap proach is also similar to [22] in its treatment of the portion of the difference stencil at the interface between the subgrids, though there is a significant difference. In [.22] these connections are calculated explicitly. This involves 16
PAGE 24
either interpolating nodal values as piecewise constants throughout cells centered at interface nodes, or interpolating, linearly, neighboring nodes. Our approach is similar to the latter technique, however, it avoids cal culating explicitly the stencil's coefficients. Rather, it uses standard fine grid (with the fine grid extended to the interface by using slave nodes), and performs the required interpolation between interface nodes prior to application of the matrix stencil as the solution algorithm pro ceeds. This latter approach is borrowed from the element methods. We note that it simplifies the derivation of the matrix stencil and allows the I use of a standard uniformgrid ste.ncil on the fine patch. Section 3.2.4 com pares briefly the relative merits of the three discretization methods treated I here. Finally, 3.3 treats the accuracy of these discretization methods. 3.1 of Composite Grids Our description of the composite grids to be used in this study applies to the case of a rectangular domain in two dimensions. We suppose the subregion of the domain requiring refinement to be an interior region that does not intersect the boundary. Although no assumptions are made con cerning the s4ape of this region, we will only consider refinement of the grid on a rectangular patch. The structure we use for the composite grid is that of [33]: Generalizations of this concept of a composite grid to Lshaped grids tnd refinement regions, other dimensions, refinements near I the boundaries of the domains, and multiple regions of refinement can easily be made. But for the purposes of this study we restrict our attention, for the most to the above simple case. Although we only describe in detail the two:dimensional case, the structure of the composite grid in one 17
PAGE 25
dimension is 'very similar, and it should be straightforward for the reader I to adapt the :definitions in this section to the situation in one dimension. Suppose 'fe are given a rectangular hgrid with uniform spacing, h being l the nodetoJ:l.ode distance in both the zand ydirections. A composite grid is determined by choosing some rectangular subset of nodes along with I a positive integer, m, indicating the order of refinement. The set of nodes forming the b"oundary of the subset is called the interface. We assume that each side of interface has at least three nodes (including comers) so ., that the interior of the rectangular subset is nonempty. This set of interior nodes is the coarse patch. H we designate the set of nodes exterior to the interface as the coarse nodes, then have the uniform hgrid .as the union of three components: coarse, and coarse patch. We call this grid the global uniform coarse grid or.:just global coarse grid. See Figure 3.1. The grid is obtained from the global coarse grid by adding nodes to the coarse patch; the fine patch component of the composite grid is the set of interior nodes of the uniform hFgrid on the region bounded by the interface, hF = h 2m. The set of bowtdary nodes for the fine patch (which the interface) is called the fine boundary, but aside from the interface nodes it is not considered part of the composite grid. However, we d,efine the complement of the interface in the fine boundary : I as an adjuvant to the composite grid called the set of slave nodes. (Though not technically belonging to the composite grid, they play an essential role in the of the composite discretization.) In the composite grid is made up of three components, coarse and interface c?mponents being the same as for the global coarse grid, along with a fine com;ponent obtained via "refinement on the patch". An example 18
PAGE 26
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 o 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Fig. 3.1. Global coarse grid with coarse ( o ), fine { ), and interface { 0) nodes. 19
PAGE 27
of a. composite grid (with one refinement) is shown in Figure 3.2. 0 0 0 0 0 0 0 'o 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Fig. 3.2. Composite grid with coarse ( o ), fine ( ), interface ( 0), and slave ( o) nodes. 3.2 Discretization Methods Next we describe three discretization methods for the following elliptic problems:: a.nd eu."(:z:) + u'(:z:) = /(:z:), :z: E 0 = [0, 1] u(O) = 1, u'(1) =a, + 'Uz(:z:,y) = /(:z:,y), (:z:,y) E 0 = [0, 1] X [0, 1] u(:z:,O) = 9t(:z:), (3.1) u(O,y) = 92(y), (3.2) 'Uz(1,y) = a(y), u(:z:, 1) = ga(:z:). Here e is a. constant, 0 < e < 1. The methods we consider a.re: the finite element, finite volume element a.nd finite volume methods. We begin by giving a. brief formal description of the two element meth ods. Let a. finitedimensional function spa.ce out of which one would 20
PAGE 28
like to a solution to the partial differential equation .Cu = f be given. We take as our approximation to u a particular linear combination of elements of a basis for C), In particular, we assume the existence of a nodal basis: {
PAGE 29
0 I h 0 0 } h/2m 0 0 Fig.' 3.3. Volumes for aportion of a twodimensional composite grid. We will return to the subspace framework shortly, but now we introduce I the third discretization method, beginning with a detailed description for ' the onedimensional model problem. Similar to FVE, the finite volume (FV) method (36] begins by partitioning the domain into a finite set of volumes in_ such a 'way that, associated with each node of the composite grid, there is precisely one volume that contains that node. Typically, the node lies in the interior of its associated volume (although it is possible to have a node at the edge of its volume, as is the case in the usual treatment of Neumann boundary conditions). Figure 3.4 shows a portion of a one dimensional grid containing a coarse node, an interface node, and two fine npdes (nodes are labeled z., here) along with the endpoints of their associated "volumes" (the e.,'s). The lengths of the volumes for 22
PAGE 30
coarse, firie, and interface nodes are, respectively, h, h/2m, and l{h+h/2m ), where m is the number of refinement levels used on the fine patch . I : 9 e h Fig. 3.4. Onedimensional composite grid with nodes (z,), and volume endpoints (e,). 3.2.1 Finite Volume Discretization in One Dimension The fuiite volume discretization in one dimension proceeds by imposing the following local conservation law, which is obtained by integrating (3.1) over each of volumes re,, ei+l]: tea+l tea+l ( EU(:c)" +u(:c)')d:c = f(:c)d:c. (3.3) These equations are then transformed into a system of equations for the ' nodal values qf u. Secondorder and firstorder terms are handled separately. For the secondorder term at the interface we have, 1eli+l "( > '( > '( > 1 u z d:c = u e,+l u ei e. ,..., u(ziH)u(:ci) u(:ci)u(zi1) "' 2(hj2m+1) 2(h/2) = ( u(Zi1)(1 + 2m)u(zi) + 2mu(:ci+1) )/h. Here, we have :used a secondorder divided difference to approximate the first derivatives. For the firstorder term (3.4) 23
PAGE 31
we need to approximate the difference on the righthand side by an expres sion in u(z) evaluated at nodes. The centered formula averages u at nodes adjacent to e. and ei+l: iii+ '( )d u(zi+ 1 ) + u(z,) u(z,_ 1 ) + u(z,) u :I: :I: .:...'''ti 2 2 u( Zi+t) u( Zit) 2 The upwind f9rmula replaces the volume endpoints with their nearest left neighbor nodes: Finally, one Ip.ight consider a higherorder upwind formula for the first order term, on the fine patch. We do this by replacing the difference on the righthand side of (3A) by its secondorder approximation at Zi+2 and the upwind n;)des Zi+t Notice that this formula cannot be applied at the leftmost fine node. At that node, the upwind stencil can be used, for example. Collecting these formulas and using similarly derived approximations at coarse, fine, and right interface nodes yields the following tridiagonal :rn.atrix stencils. Centered at the generic node n, these stencils involve neighboring nodes in an obvious way expressed by Upwind Coarse Nodes: [eh, 2e + h, e]/h 24
PAGE 32
Left Interface: I llight Interface: Fine Nodes: Centered Coarse. No'des: [ ,h/2, 2t:, + h/2)/ h ., Left Interface: .. Right. Iptetface: Fine Nodes: HigherOrd,er Upwind Fine Nodes: 25
PAGE 33
This stencil involves the nodes In addition, the Neumann condition at the right boundary, u'{l) = a, is handled as follows. Here, the volume for the node at the right boundary is a halfvolume. Figure 3.5. 9 I f h Fig. 3.5. Onedimensional grid with a right Neumann boundary. We have and These yield the following stencils and righthand sides. Neumann Boundary Upwind: [eh, e+h]/h h h . f{l) + 2 4 Centered: h h [h/2, + h/2]/h = 2. /{14) + 26
PAGE 34
These stencil's, centered at node n, correspond to [nw, n]. Generally, for the righthand side of (3.3), letting l'i = [ei, ei+1], we use Finally, we note that in treating Dirichlet boundary conditions, the value of u at a Dirichlet node is given by the specified boundary value. Thus, the equation for such a node is a trivial one that is easily eliminated from the tota,J. system of equations by modification of its righthand side. The use of this approach obviates the need for volumes associated with Dirichlet nodes and, as a result, the volumes do not truly partition the domain. A discussion of the effect of this treatment in terms of loss of conservation in volume methods, along with an alternative for maintaining the partitioning, is given in [33, Chapter 2]. 3.2.2 Discretization by Element Methods Reconsidering the FE and FVE discretization& discussed earlier, it is interesting,to note that, in one dimension, the centered version ofthe stencil just derived using the FV method can be obtained by applying either the FE or FVE (using the volumes just described) to the nodal basis pictured in Figure 3.6c, where nr is the left interface node. {Only the nonzero values of the basis functions are graphed.) The basis functions corresponding ,to coarse nodes are familiar piecewise linear hatfunctions for the uniiorrp. hgrid (i.e., the global coarse grid; see Figure 3.6a) and the basis functions corresponding to the fine nodes are similar functions corresponding to a uniform h/2mgrid (which we call the global fine grid; I see .3.6b ). We call these bases (Figures 3.6ac, respectively) the 27
PAGE 35
global coarse, global fine and composite bases, and denote their respective spans by' c), W, and i. The relationship between the three bases may be succinctly represented by the equations p = l'l/Jo + 1/Jp + l'r/Jq = f/Jp r = l'r/Jq + 'r/Jr = f/Jr l'r/J = 1/J. Now, let WF 'be the span of and let We be the span of basis functions for 'W corresponding to nodes other than n., nt, and nu. Then any comp?site basis element may be written uniquely as = r,e + f/Jr,F, where f/Jr,c E We and f/Jr,F E WF. Thus, = f/Jr,e if nr is a coarse or interface node = f/Jr,F = 'r/Jr if nr is a fine node, i.e., for and interface nodes of the composite grid, the composite : nodal basis element is obtained from the corresponding global coarse basis element [f/J by subtracting off its component in W F (for coarse nodes this component is zero), and, for fine nodes, is the global fine basis element with that node. We note that this construction of a composite basis produces a space that is nested :between the global coarse and global fine spaces: c) c i c w. (3.5) We call (3:5) the conformity relationship between these subspaces. It is I useful for o:btaining accuracy estimates for the finite element discretization (see Section 3.3). This leads us to consider FE and FVE discretization& in two dimensions. We need a :finitedimensional space out of which an approximation to the solution of our differential equation may be constructed. Suppose we define, 28
PAGE 36
(a) 1 n,. n. (b) np n,. n . (c) Fig. Nodal bases for global coarse (a}, global fine (b), and composite (c) spaces. 29
PAGE 37
in the usual.way, nodal global coarse and global fine bases consisting of piecewise hat functions on triangulations as shown in Figure 3.7. Because of the similarity of the component structure (the division into coarse, interface, and fine components) of the composite grid in one and two dimensions the specification given above applies directly in two di mensions to the construction of a nodal composite basis for such a space. . However, :the process of discretization is now complicated by the nonuniformity of nodes in two directions at the interface. Such nonuniformity is not a problem at coarse or fine nodes where, as in one dimension, construc tion of the m.atrix stencil proceeds along familiar lines (discretization on uniform grids., with given mesh widths). Yet, the difficulties at the interface can be made :to virtually disappear by viewing the task locally from the point of view :of the global fine basis. To be consider stencils involving interfacetofine connections. These can be obtained easily by taking advantage of the structure of the composite basis. For example, in the FE case the computation of (p, q); t/Jq and Pare composite basis functions corresponding to interface .. node; nq and neighboring fine node np, respectively, is facilitated by viewing t/Jq on the same (global fine) level as p (since P corresponds to a fine it is equal to a global fine basis function), i.e., as a linear combination of global fine grid basis functions (namely, t/Jq = t/Jq,C ). This reduces the pr,oblem to using existing finite element software to compute (p, 11), p and 11 are both elements of the global fine basis, and to combining such contributions. In a sense, this process uses slave nodes basis functio:p.s corresponding to these nodes are introduced temporarily (see Figures A similar technique can be applied in a FVE dis cretization, though in practice the procedure is not quite as simple due 30
PAGE 38
(a) (b) fig. 3. 7. Global coarse {a) and global fine (b) triangulations. 31
PAGE 39
to the additional nonuniformity in volumes at the interface. On the other hand, FVE stencils are more sparse than FE stencils, additional localization induced by the volumes. Next we examine the finetointerface connections. We emphasize the construction of this portion of the stencil in the subspace context since it provides guiqance for the construction in the absence of a subspace, for in the FV discretization. In order to describe the treatment of this part of the discretization, it is useful to introduce the subspace w,, which is defined to he the span of those global fine basis elements that correspond to nodes at the boundary of the fine (interface nodes and slave nodes this subspace comprises a fine for the patch), and cir, which is defined to he the span of the composite basis functions corresponding to interface nodes. For convenience, we summarize the decomposition of the global fine subspace, '11, and the composite subspace, ci, below: '11 F Span of '1/J 's corresponding to fine patch nodes, We Spa11: of '1/J's corresponding to other than fine patch nodes, Wr,Span of,P's corresponding to interface and slave nodes (w, C We), ci 1 S:IJan :of 's corresponding to interface nodes. Here, '1/J denote generic elements of the global fine and composite nodal bases, respectively. We also note the following localness property of the FE and FVE composite grid discretization&: each composite basis function corresponding to an interface node enters into the computation of a finetointerface connection only through its components in Wr,. This component of, (where n, is an interface node) is obtained by restricting 32
PAGE 40
the = l/Jl,C to the boundary of the fine patch: (3.6) Here, the .restriction denotes the presence on the righthand side of (3.6) of basis elements for q;, only. For each this relation provides a unique representation of the component along the interface in terms of the basis for q;,, and, provides a welldefined linear mapping, n : il w,, of the interface to the fine boundary. We denote this linear operator as interpolation at the interface. Now, due to the localness property, it does not matter whether we compute the finetointerface stencil with respect to the at the interface or with respect to 'I!,. It follows that this stencil may be computed using just the basis elements from q; F and Wp, (i.e., the nonU,niformity between fine patch and interface may be avoided by t.he interface with a finemesh boundary). This means that when the fine;tointerface stencil is applied (during matrix multiplication) to an interface component, interpolation must first be performed: (u.1 is a given f'!lnction in ci;I, and up, is its component in the fine boundary), and then the matrix stencil is applied to up,. As an exa$ple of this process, consider the case of a composite basis function the left vertical edge of the interface when there is one level of refinement. Then the representation in terms of the basis for q;, is 33
PAGE 41
Fig. 3.8. Support of a composite basis function corresponding to an interface node ( 0) at a left vertical interface ..P+ Fig. 3.9. Support of a composite basis function (shown twice) at a left vertical interface and support of the functions in its component along the fine boundary. 34
PAGE 42
(n_ and n+ are the slave nodes below and above node nl) The support of is show.n in Figure 3.8. It is shown again in Figure 3.9 along with the supports of the three basis functions in its component along the fine boundary. Representations along other edges are similar, meaning that interpolation at the interface proceeds as a sequence of onedimensional interpolations along the edges of the interface: the nodal value of a function u, equals the corresponding value of ur if the node is an interface node, and is equal to the average of its neighbors' nodal values if the node is a slave. With multiple levels the process may be carried out using a reapplication of the singlelevel process on the succession of levels leading from the interface to the fine boundary. We note that due to the localness property, this procedure is equally valid for FVE and FE discretizations. Having arrived at a satisfactory way of handling this part of the discretization when it is performed with respect to a composite function space, it may also be used to facilitate this part of the discretization in methods that do not utilize such a space, as is the case, for instance, with the finite volume method. Normally, with the latter method, treatment of the interface necessarily involves the use of interpolation at interface nodes when calculating explicitly the fineto interface s,tencil [22]. On the other hand, the approach described above defines this stencil implicitly, and allows for the use of a standard uniform grid on the union of the fine patch and the fine boundary. At this point, it is useful to say a few words about the way that the coefficient matrix L of the composite grid problem is represented in prac tice. Three matrix stencils are called for: a coarse stencil that involves coarse and interface nodes, a fine stencil that involves fine patch nodes augmented by a fine boundary, and an interface stencil that involves fine 35
PAGE 43
patch, interface, and coarse nodes. In our computer codes, we use sepa rate lexicographic orderings of the fine patch, and of the coarse component imbedded in a global coarse grid. This means that the coarse and fine stencils are standard uniform grid stencils and may be utilized as such. In forming the fine patch component of the product LJJ, for example, storage for the fine patch component is made large enough to accomodate the fine boundary. Prior to applying the stencil to this component, the interpolant of the interface of y is stored in these locations. In forming the coarse component, computations proceed in the usual manner, except that they are restricted to the coarse and interface components Finally, the interface portion of the product requires some special handling. In our im plementation, the interface resides together with the coarse component as part of an array representing a lexicographically ordered global coarse vector. Therefore, interfacetocoarse and interfacetointerface computations proceed as they would on a uniform grid, while interfacetofine connec tions are handled by a separate code segment that requires access to the patch component of y. 3.2.3 Finite Volume Discretization in Two Dimensions For the most part, finite volume discretization of the problem (3.2) is not much more involved than that of (3.1), since the method in two dimensions is very similar to that in one dimension. To see this, consider the portion of a twodimensional composite grid near a left vertical (west) interface, as shown in Figure 3.10. The figure depicts (left to right) a coarse, an interface, and three fine nodes along with their associated volumes. Let V be the interface volume and label its sides in an obvious way as N, S, E, and W. By Green's Theorem, 36
PAGE 44
(i:,y) 0 0 Fig. 3.10. Nodes and volumes near a left vertical interface. fv Au dV = k Uzdy fw Uzdy + JN u11dz Is u11dz. Consider, for example, the first term on the right. We approximate Uz at the Iliidpoint of the east side of V using a centered difference formula involving u(i:,,y) and u(i: + h/2m, y) and integrate using the Iliidpoint rule: h (u((i: + h/2m+l) + h/2m+l ,y) u((i: + h/2m+l)h/2m+l ,y)) 2(h/2m+l) =2m (u(i: + h/2m,y) u(i:,y)). Using a similar approach for, say, the third term gives .f d h _1_) u(i:,y +h) u(i:,y) JN 1Ly Z 2 + 2m+1 2(h/2) 37
PAGE 45
= + (u(z,y +h)u(z,y)). For the firstorder term, using an upwind approach, we have h(u(z,y)u(zh,y)). Let s = l + 2,;+1 Using these approaches leads to the following stencils for the upwind version of the FV discretization, which are centered at node n and involve the nodes Upwind Coarse Nodes: [e, eh, 4e + h, e, e] West Interface: [se, eh, (1 + 2s + 2m)e + h, 2me, se] East Interface: South Interface: North Interface: [e, s(e +h), (1 + 2s + 2m)e + sh, se, 2me] 38
PAGE 46
Fine Nodes: [e, Eh/2m, 4e + h/2m, e, e] An obvious modification of the handling of the firstorder term (see the deta.ils of the onedimensional FV discretization) yields the following centered stencils. Centered: Coarse Nodes: [e, eh/2, 4e, e + h/2, e] West Interface: East Interface: [8E, 2mfh/2, (1 + 28 + 2m)e, E + h/2, 8E] South [2me, 8(E + h/2), (1 + 28 + 2m)f, 8(Eh/2), e] North Interface: [e, s(e + h/2), (1 + 2s + 2m)e, s(eh/2), 2mf] Fine [e, Eh/2m+l, 4e, E + h/2m+l, e] 39
PAGE 47
We also provide the stencils for the Neumann condition at the right boundary (with their associated righthand sides) and a fine patch stencil for a higherorder upwinding. Neumann BQundary Upwind: [e/2, eh, 2e + h, h2 h , e/2] = 2/(14,y) + eh a(n) Centered: [e/2, eh/2, 2e + h/2, h2 h , e/2] = 2 /(14,y) + eh a(n) Here, "" indicates the node nE lying outside of the stencil. HigherOrder Upwind Fine Nodes: This final stencil is associated with the nodes As for the treatment of Dirichlet nodes, discretization at these points is performed as described for the onedimensional case at the end of Section 3.2.1. We note that the above formulas for the interface nodes do not apply at its corners. For these nodes we use the standard stencil for the coarse 40
PAGE 48
nodes. Also, 'the fine grid stencils apply throughout the fine patch, assum ing that the process of interpolation of interface nodes to a fine boundary as described earlier in this chapter is incorporated into the definition of the discretization at finetointerface connections. 3.2.4 A Comparison of Discretization Methods In thiE; section, we briefly consider the relative merits of the discretiza tion methods described in this chapter. For the most part, the results of this study pertain to finite volume discretization. We have chosen to em phasize this method primarily due to its ease of use, first, in deriving the composit,e grid equations, and second, in terms of solving these equations. One attractive feature of this discretization is the simplicity of its ma trix stencil. With the approach used here, this method gives rise to a fivepoint stencil throughout the composite grid. This is due to the way the interfacetopatch connections have been defined. The result is that code involving the stencil is particularly easy to program. The element stencils differ from this in that they distribute the coefficient of this connection locally among nodes of the fine patch. The result is a stencil that becomes less sparse at the interface as the patch is refined. In the context of the problem, ellu(z,y) + 1Lz(z,y) = /(z,y), (3.7) with e small, a more compelling reason for appealing to the finite volume method is that its stencils appear on uniform grids as familiar finite differ ence The efficiency of the numerical methods of this study depend strongly on the ability to solve discrete versions of (3.7) on the fine patch and on the global coarse grid, i.e., on uniform grids. Relaxation serves as the basis for solution of these problems by multigrid and, due to the fact 41
PAGE 49
that the finite difference matrices capture all of the advection of (3.7) in their main tridiagonal blocks, line relaxation is extremely effective (see Sec tion 5.3.1). We note that it is more difficult to identify robust relaxation strategies (i.e., ones that work well over wide parameter ranges) for the element methods. This is due to the appearance of advection (correspond ing to the in (3.7)) in portions of the stencil that correspond essentially to discretization of yderivatives. For example, the uniform grid stencils for finite volume elements, which is analogous to centered differencing, is (h/8, h/8, e3h/8, 4e, e + 3h/8, h/8, e + h/8], which corresponds to the nodes This drawback does not preclude the use of the linear equations solvers of this study when applied to the equations derived by the element methods. However, more work needs to be done to develop the appropriate relaxation strategies. Finally, we note that although the finite volume stencils appear largely like finite difference ones, there is a methodology present in the former approach that is lacking in the latter. This methodology provides guidance .for the systematic treatment of the discretization, especially as the grid becomes nonuniform, i.e., as the coarse component meets the fine one. The systematic nature of the discretization is also an important and attractive feature of the element methods, FE and FVE. In addition to facilitating the treatment of the interface portion of the discretization, it also simplifies the choice of interpolation and restriction operators on 42
PAGE 50
which depends. With conformal subspaces, i.e., C the interpolation operator, is determined by representing each basis element for in the basis for In finite elements, restriction, is also determined automatically by requiring that the so called va.rjational equations that relate the coarse grid discretization to the fine one [30] hold. With finite volume elements the choice of restriction may be guided by, physical principles (see [33, Chapter 3]). With finite volume elements the choice is more arbitrary. In this study, when grid transfers are called for we employ the sevenpoint prolongation and restriction operators of [25, Chapter 3] for the fine patch component, and the trivial (identity) mapping: for the coarse and interface components. Finally, we make note of the fact that upwind discretizations of (3.7) play a very significant role in this study. We have found it easy to in corporate firstorder and secondorder upwinding into the finite volume discretization. In finite elements the use of upwinding is also common place. On the other hand, a satisfactory strategy for performing upwind has not yet been developed with respect to the finite volume element method. 3.3 Accuracy of the Discretizations We now turn to a discussion of the accuracy of the discretizations de scribed in this chapter. Of particular interest are the finite volume dis cretization:s that we have derived in some detail for the model problem (3.2). The reasons for this interest are at least twofold. First, the finite volume discretization as described for the twodimensional case is some what novel and the matter of its accuracy is of interest in itself. And second, the analysis of algorithms for solving composite grid equations 43
PAGE 51
that will be pursued later in this study (in Chapter 4) is based on their behavior respect to the finite volume discretization, so it is important to establish the validity of this type of discretization. In the next section, we will assess the accuracy of this discretization numerically for a sampling of test problems. This section is a more general discussion of the accuracy of the discretization& considered above. It begins by taking note of recent results for FV and FVE discretization& of elliptic problems on nonuniform grids. Prior1 to the recent studies of Cai et al. (10), (11), little had been done to prove the convergence of finite volume element discretization&. Some prior work had been done in this context for finite volume discretizations on uni form meshes (references to this work are contained in (11] and (42]), and in certain instances these turn out to be special cases of FVE discretization&. The work of Cai et al., however, was the first to incorporate finite element theory ill.to the FVE framework in order to produce finite elementlike accuracy estimates. This work also includes a systematic study of con vergence on composite and other nonuniform grids. In [10] an O(h312 ) bound in a discrete energy norm is proved for the Poisson equation using a composite grid FVE discretization similar to the one described herein, with piecewise linear functions used for the approximating subspace ( h here den9tes the coarse component meshwidth). This result is strength ened to O(h2 ) by incorporating bilinear functions into the discretization at the interface. Similarly, in (10) an O(h) result in the discrete H1 norm is obtained for
PAGE 52
using a finite volume discretization. His analysis, which applies to Carte sian product, grids, involves reformulating the finite volume method as a PetrovGalerkin finite element method. This yields estimates in the discrete iJl norm when u E Hl+v (1/2 < u < 2). Also, the case of a finite volume discretization for composite grids has been analyzed in [22]. This corresponds closely to the method used here, but differs in using a strictly cellcentered grid. There, a bound similar to the previous one is obtained I for the continuous H1 norm. We note that for finite elements on composite grids, accuracy estimates in the coarse component meshwidth are facilitated by the conformity rela tionship (3.5). Consider, for example, problem (3.2) with the given boundary conditions replaced by u = 0 on all of 80. Let be a global coarse space defined by using piecewise linear functions on a uniform triangular tions as indicated earlier. Associate with the corresponding global coarse grid a composite grid via the procedure of patch refinement. Let wh be the global fine space with the meshwidth of the fine patch, and let
PAGE 53
to arrive at H1 estimate (3.8) where c 4/ e (is independent of h). The H1 norm is given by Now, since then on the right may be replaced by the interpolant of u from This then can be used (see [41], for example) to obtain an O(h) estimate for lluAlso, Nitsche's method [41], [16] is applicable, yielding an O(h2 ) estimate in the L2 norm. We note that this result is somewhat unsatisfactory in light of the e1 I present in the bound (3.8). However, if e is not particularly small the result may be meaningful, and, in any case, convergence of the composite solution to the solution is guaranteed as the coarse component is refined. The result is independent of the size, position, and level of refine ment of the composite grid's patch, i.e., it made no use of the inequality It is clear that variations in the patch must affect the accuracy of the discrete solution (our motivation for using local refinement is based on this assumptio.n) .. At one extreme, if the patch is reduced to a single point, or if no actual refinement is used, then the estimates in h still hold. At the other extreme, if the patch is allowed to cover the entire domain and m > 0 are used, then the above estimates hold with h replaced by hF = h/2m. Using nontrivial refinement on a judiciously chosen patch, we expect accuracy in the discrete solution that lies somewhere between these two estimates. 46
PAGE 54
In particular, this will be true if u has small derivatives outside the patch. Suppose we partition the domain into two subregions, 0 = Oc U OF, where Or is the rectangle with corners corresponding to those of the interface. Then, for vh E ci;h we have Now let vh assume the role of the interpolant of u from the composite space . ci;h Due to the particular form of ci;h, vh may be constructed using basis elements' for the global fine space 'IJ.ih on the region OF, and basis elements for the global coarse space
PAGE 55
I patch (recall Figure 3.2). More work is called for in order to obtain an interpolation. estimate in the region lying between Op and 00 Intuitively, one should be able to obtain an 0( h) bound here. Rather than pursue this, let us suppose that the component of the interpolant on OF does actually interpolate u at the fine boundary. This will happen, for example, if the derivatives of u are zero in the complement of Slp. Then we have the estimate Notice that, if the derivatives are indeed zero outside of Slp, then the first term on 'the right is zero and we obtain O(hF) convergence as the patch is refined. This:is the kind of behavior we expect, generally, if the derivatives of u are sinall outside of some region and the patch is sufficiently large to encompass the region. On the other hand, if u has steep gradients outside of the patch, then the best we can expect is 0( h) convergence as the coarse II grid is refined. We note that a bound similar to the one above is obtained in [22] a :ijnite volume discretization on a cellcentered composite grid. Although these observations have been made with regard to finite ele' ment we expect similar accuracy for appropriate :finite volume and finite volume element discretizations. In particular, the above example is analogous to the centered version of the :finite volume discretization used. here. In the next section, some examples are presented that indicate the actual behavior of:finite volume discretizations of problem (3.2). we consider the use of upwind or centered approximations on the coarse of the grid and allow for upwind, higherorder up wind, and centered approximations to be used on the fine patch. Generally, 48
PAGE 56
our computational results agree with the behavior predicted by the above analysis. However, there is still much work to be done in terms of estimates for the of composite grid discretizations. In particular, we have found that the use of upwinding here gives rise to accurate solutions and very fast 'discrete solvers. Therefore, it would be appropriate to pursue its use with respect to the element methods. One attractive possiblity is that of applying the streamline diffusion method (26] in this context. Before presenting the results of our test problems, we make one final ob servation with respect to accuracy and the relationship between the and fine components in the context of the advectiondominated problem. We have observed that when the derivatives of u. are small outside the an accurate solution can be obtained with sufficient patch re finement. Although this is the case, one should not expect inaccuracies to then be restricted entirely to the fine grid. Because the problem is advection dominated, the solution at a given point is dependent on the so lution from that point. A similar dependency should be expected in the solution. In particular, the solution on the coarse com ponent downwind from the patch will inherit the innaccuracy that is a result of insufficient refinement on the fine component. Therefore, it is possible for the true solution to be uniformly wellbehaved throughout the porti<;m Qf the domain corresponding to the coarse component, while the accuracy of its discrete counterpart changes significantly in this region. This behavior of the discrete solution will be observed in the test problems of the next section. 3.3.1 Problems in Two Dimensions The following test problems are constructed from (3.2) by determining f 49
PAGE 57
and the boundary data from the specified solution. We employ constants r > 0 and 0 < :co, Yo < 1. With respect to the point (:co, Yo), define D(:c,y) = J<:c:co)2 + (yYo)2 Problem 3.1. The solution. is 3 u(:c,y) = 2 ((D(:c,y)/r)3 '2(D(:c,y)/r)2), D(:c,y)::; r, u(:c,y) = 1, D(:c,y) > r. This solution has been designed to be smooth at the interface between its two components. Specifically, at the set of points D(:c,y) = r, u is continuous and differentiable, and its gradient is zero: Uz(:c,y) = 6(:.::co)(D(:c,y)/r 1)/r2 u11(:c,y) = 6(y1)/r2 We note that u is independent of E, and the extent to which the function is wellbehaved depends on the magnitude of r. When this quantity is small, the nonconstant component of u covers a small region, however its gradient there is to 1/r. On the other hand, if r is kept bounded away from zero, we have a moderately wellbehaved solution irrespective of the choice of't:. The graph of u(:c,y0 ) with :.:0 = y0 = 1/2 and 1' = 1/4 is shown in Figure 3.11; the graph of u(:c,y) is obtained by rotating it about the uaxis at (:co, yo). Problem 3.2. The solution is 50
PAGE 58
We remark that u and its derivatives are small except in a neighborhood of (:z:0,y0), where the function undergoes rapid growth. For simplicity, let r = 1 and :z:o.= Yo= 0. Then 'IL:.: = u 2z/e, = (2u/e) (2:z:2 /e1), so that u has its maximum rate of change at D(:z:,y) = y'i., with rate I!Vull = 2/(e ../f). This represents a sharper front to be resolved than in Problem 3 .1. 5 4 3 2 0.0 0.2 0.4 0.6 0.8 1.0 Fig. 3.11. Graph of the solution of Problem 3.1, withy= y0 Convergence rates as functions of grid refinement for these test prob lems are presented in Tables 3.13.5. Relative errors involving the discrete solution of the problems and the known solution evaluated at the nodes are presented. The discrete maximum norm, II lloo, is computed on the composite grid. The discrete L2 norm, II 112 is computed on a global fine grid having the same meshsize as the fine patch. The errors are com puted on the latter grid by interpolating the discrete solution to it from 51
PAGE 59
the composite grid, while the true solution is simply evaluated there. Var ious combinations of finite volume (upwind and centered) discretizations i are used on the coarse and fine components of the composite grid. In ad dition, in some instances a discrete solution is computed on a global fine grid with the same refinement as the fine patch in order to compare the composite gri.d solution with a uniform grid solution. The discrete solution on all grids is computed by an iterative method and using the convergence criterion that the relative residual satisfy Here, ri is the ith composite residual and r_0 is the original righthand side of the grid equations. The iterative solver used for these tests is a particular implementation of the twolevel FAC algorithm described in Chapter 4' of this study (see, also, Section 5.3.1). This is not the most efficient algorithm presented herein, but its use is appropriate for these tests in that the "exact" solution of the composite grid equations is desired. An extensive of the efficiency of this algorithm and alternative ones is presented in Chapters 4 and 5. In what follows, the meshwidth of the coarse component is denoted by h = 1/2T, for some positive integer me. A positive integer m denotes the additional number of refinements in mesh width (each by a factor of two) in going from the coarse component to the fine patch. Hence, the fine patch has a of hF = h/2m. 3.13.2 show the results for Problem 3.1 with r = 3/32, :z:0 = 1/2, and y0 = 3/8. The coarse component meshwidth, h = 1/32, is kept constant and, the patch is successively refined. The southwest corner of the I interlace is located at (12/32, 7 /32), with :z: and y dimensions of 8/32 and 52
PAGE 60
10/32, respectively. The patch here is made large enough to encompass the nonconstant component of u. Upwind discretization is used on the coarse component of the composite grid, and two test results are shown, using upwinding and higherorder upwinding (Hup) on the patch. The diffusion coefficient is E = 104 Table 3.1 Relative errors in the mazimum norm as functions of patch refinement for Problem 3.1. Composite Grids Global Fine Grid Upwind/Upwind Upwind/ Hup Upwind m llehlloo/llulloo II eh II oo/ !lull oo llehlloo/llulloo 1 1.383101 5.551 102 1.383. Iot 2 6.569102 1.549. 10 2 6.568 102 :.3 3.194 :.1 3.977. 10 a 3.194. 10 2 4 1.576. 102 9.605. 104 1.576. 102 5 7.828. 10 3 2.273. 10 4 7.808. 10 3 The error norms show O(hF) and O(h}.) convergence for the upwind and higherorder upwind discretizations, respectively. In addition, the re sults for .. the upwind composite solution and global fine solution are nearly identical, with the composite solution showing a mild degradation in the L2 norm. m 1 2 3 4 5 Table 3.2 Relative errors in the L2 norm as functions of patch refinement for Problem 3.1. Composite Grids Global Fine Grid Upwind/ Upwind Upwind/ Hup Upwind llehll2/llull2 llehll2/llull2 llehll2/llull2 1.069. 102 3.670103 1.095102 5.359 10 3 9.532. 10 4 5.422. 10 3 2.687. 103 2.500. 104 2.703 103 1.346. 10 3 6.589. 10 & 1.350. 10 3 6.735. 104 1.663. Io& 6.750. 104 53
PAGE 61
Tables 3.33.4 show the results for Problem 3.2 with r = 1, :c0 = y0 = 1/2, and = 103 Standard upwinding is used on the coarse component of the composite grid, with meshwidth h = 1/32. The higherorder upwind and centered discretization& were used on a patch centered at (:c0,y0 ) with length and height equal to 1/4. Again, the patch is made sufficiently large, this time to enclose the solution's sharp front. The problem is also solved using the centered method on the global fine grid. The convergence rate for the discretization& is again with close agreement between the composite and global fine solutions. m 1 2 3 4 5 Table 3.3 Relative errors in the maximum norm as functions of patch refinement for Problem 9.2. Composite Grids Global Fine Grid Upwind/ Hup Upwind/Centered Centered ll.ehlloo/llulloo l!ehlloo/llulloo l!ehlloo/llulloo 9.538. 10 2 1.051. 10 1 1.050. 10 1 3,.271 102 2.253. 102 2.253 2 9,.279 103 5.460 .IQ3 5.460. 10 3 2.397 3 1.356. 10 3 1.356. 10 3 6'.021 1o4 3.383. 104 3.383 4 Table 3.5 shows the results for a different set of tests based on Problem 3.2, measuring convergence as a function of h. Here, is increased to 102 allowing the front to spread out while the location and dimensions of the interface are left the same as for the last set of tests. Again, the upwind/centered and centered/centered combinations are used, but now on three composite grids, each having two levels of patch refinement and varying coarse component meshwidths. 54
PAGE 62
m 1 2 3 4 5 Table 3.4 Relative errors in the L2 norm as functions of patch refinement for Problem 9.2. Composite Grids Global Fine Grid Upwind/ Hup Upwind/Centered Centered II ehll2/ llull2 llehll2/llull2 llehll2/llull2 1.038. 10 1 8.449. 10 2 8.436 2 3.159. 10 2 1.858. 10 2 1.857. 10 2 8.533 .lQ3 4.519 .lQ3 4.520. 103 2.184 3 1.122. 10 3 1.123. 10 3 5.504 4 2.800 4 2.802 .lQ4 We make the observation that the patch is "too small" for this problem, in the sense that the innacuracy on the coarse component is inherited by the patch (and also affects the coarse component downwind from the patch). 'As a result, the effectiveness of increasing the patch refinement has deteriorated somewhat. For the upwind/centered discretization with me = 4 and m = 2 (second row, first column of Table 3.5), the relative Table 3.5 Relative errors as functions of coarse component refinement for Problem 9.2. Composite Grids Upwind/ Centered Centered/Centered me {m = 2) llehlloo/llulloo llehlloo/llulloo 3 1.963. 10 1 1.935. 10 1 4 7.981 2 7.468 .lQ2 5 3.921. 10 2 1.903. 10 2 6 2.152 2 5.628 3 L2 error restricted to the fine patch is approximately 7.67 102 With two additional levels of patch refinement (m = 4) this becomes 3.78 55
PAGE 63
102 a modest improvement (in particular since the secondorder accurate centered discretization is used on the patch). On the other hand, this additional local refinement has a pronounced effect on the solution as a whole, as is apparent from graphs of the solution presented in Figures 3.123.14. In Figure 3.12 the exact solution of Problem 3.2 is plotted. Figures 3.13 ,and 3.14 show plots corresponding to the two cases of patch refinement just mentioned. In each case, for ease in plotting, the fine patch solution has been transfered by restriction (see Section 3.2.4) to a coarse patch having the same refinement as the coarse component. We emphasize the quality of the coarse component solutions here. Recall that the solution should become uniformly small as the radial distance from the center of the domain increases. In Figure 3.13, this fails to happen due to the innacuracy of the coarse component solution in the region downwind from the patch. In Figure 3.14, this innacuracy of the solution has been removed by using purely local refinement upwind from this region, i.e., by adding additional refinement on the patch. This example shows that, in the of advectiondominated problems, local refinement can have a profound effect on the global solution. 56
PAGE 64
1.0 Fig. 3.12. The exact solution of Problem 9.2. Fig. 3.13. The approximate solution of Problem 9.2 with two levels of patch refinement. 57 ..
PAGE 65
58
PAGE 66
4. Convergence of FAC This chapter studies the convergence of FAC as an iterative solver for advectiondominated advectiondiffusion equations on composite grids. We begin with only a few assumptions regarding the structure and origin of the composite grid equations. The composite grid consists of coarse, interface and fine patch components, and coincides with a global coarse grid in its coarse and interface components. Its fine patch may be viewed as being derived, by refinement, from a coarse patch on the global grid. As members of the composite grid space are represented by lower case Roman letters distinguished by an underscore, and components by the lower case letters subscripted with 0, F, or I. For example, a composite righthand side (or residual) is written Similar notation is used for composite grid operators and their components. The composite grid equation Lz = r may then be written as [ Lc 0 0 LF Lrc LrF Lcr l LFr Lr (4.1) A word of explanation with regard to the notation accompanying the com ponents of L: in general, Lxy denotes that portion of the stencil of the matrix L that represents the connections to nodes in component Y that appear in the equations for the nodes in component X, i.e., entries Iii where node i lies in X and node j lies in Y. Also, Lr, for example, stands for the 59
PAGE 67
interfacetoip.terface connections. Section 4.1 below introduces the two level version of FAC for general composite grid equations. The results in this section utilize only the algebraic structure of the composite equations and are not at all dependent on the underlying continuous problem or its discretization. Sections 4.24.3 are devoted to studying the convergence of this when it is applied to the finite volume discretization& of the advectiondominated problems (3.1) and (3.2), respectively. 4.1 The TwoLevel Method Assuming that the matrices L, Lc, and LF are invertible, then per forming a few steps of block reduction leads to the equivalent form where (4.2) is the Schur complement in L with respect to the partitioning (4.1), and A L L1 L L1 rr = rr rc c rc IF F rF. In principle, then, the composite grid equations could be solved directly by the following threestep procedure. Step 1. Solve L1zr = ir for z1 Step 2. Solve .Lczc = rc Lcrzr for zc. Step 3. :LFZF = rFLFizr for ZF. We note that this is just the familiar Schur complement method for solving a system partitioned as in ( 4.1 ). 60
PAGE 68
FAC [32], [33] is an adaptation of multigrid to composite grid problems and, as such, makes use of a coarse grid discretization of the differential operator to obtain approximations to the fine grid error left by relaxation. In the present context, one may think of the fine grid operator as being the composite grid operator defined above. As for the coarse grid operato:r, it is useful ill; defining it to retain some :flexiblity. We accept in this role a discretization on any of the composite grids that lie between the global coarse grid and the true composite grid (the grid on which the discrete solution is ultimately sought). These grids are the ones obtained from the global coarse grid in the same manner as the true composite grid, but differ by having less refinement on the patch (all grids agree in their coarse and interface components). An. important special case to consider is the global coarse grid viewed as a composite grid. Since this grid is uniform, in principle, no special considerations need apply in defining the operator, that is, a standard discretization based on a lexicographic ordering of the nodes and equations, for example, may be used. Also, since the grid is uniform, it is possible to apply any one of a variety of wellknown, highly effective methods to solve standard problems on this grid. The uniformity of the global coarse grid is its attractive feature. However, we :retain the viewpoint put forth while defining this grid (see Section 3.1) that it is itself a composite grid and that the discretization on it yields a composite grid operator. It makes sense, then, to refer in general to the coarse grid as the coarse composite grid (the phrase composite grid may be reserved for the true composite grid when distinguishing the levels in the twolevel Once this viewpoint is taken, the strategy of using different discretizations on the coarse and patch components of the coarse grid becomes a possibility for all choices of this grid. We have already 61
PAGE 69
considered., in. Section 3.3.1, the effect on the accuracy of the composite grid solution of using different discretization& on the components of the composite grid. We investigate in this chapter the effect of using various "' discretization&; on the coarse grid in terms of FAC convergence. When :denoting the grid operators, Lis used as above in denoting the composite; operator and its components, and A is used similarly for the coarse (composite) grid operator. Generally, notation developed for global coarse objects ;is used to refer to coarse composite ones. So, for example, the refine:r;nent region on the coarse composite grid is referred to as "the coarse patch." Let A pe the discrete operator on this grid. Its component form is [ Ac 0 Acr l A= 0 Ap API Arc Arp Ar Here, the subscript P refers to the coarse patch, as opposed to F, which refers to tlie fine patch component of the (true) composite grid. ,, The following description of FAC corresponds to the delayed correction version in its form as described in [33]. Here, a subscript G is used to distinguish coarse composite vectors and components from true I composite ones. FAC algorithm for the solution of Lz = r Let be gjven and set r.0 = r. Lz0 Loop on i: i = 0, 1, ... until convergence. Step 1. Solve Azc = rc = = for ZG = (zc,c, ZG,P1 ZG,I ). 62
PAGE 70
' ' Step 3. Perform the composite correction '+1 Zc = Zc + ZG,C = + zF i+l i ZI = ZI + ZG,l Step 4. Set z.i+l = (zi+l zi+l zi+1 ) and form the new composite residual C'F'I' !:i+:i = rLzi+l. End loop. To the.algorithm in words, each iteration begins by transferring the to the coarse patch (dropping the iteration indices): rp +IC,rF. 1 Here, IC represents a restriction operator mapping the fine patch to the patch. Then, rp and the other components of the com posite form the righthand side for the coarse composite equations. This is solved to obtain zc, a global coarse (or coarse composite) grid to the composite grid error. To obtain a fine patch error the local patch problem is solved with the fine patch residual tlle righthand side. Dirichlet values are provided for this prob lem by the interface values of zc. (As a practical matter, recalling the discussion in,,Section 3.2.2 on the treatment of the finetointerface stencil, in order to f6rm the vector of boundary values in that context, the inter face com'ponent of zc is necessarily mapped to a refined interface using the interpolation', operator Jj.) The solution zF offhis problem, along with the coarse and components of zc, is used in the final step to define a 63
PAGE 71
composite correction, which is added to the current approximation to the solution. Finally, the residual is updated. Havin:g dtj:fined the FAC algorithm, we now turn to a study of its conver gence properties, concentrating on an expression for the updating of the composite residual. In what follows, we assume appropriate consistency [17] for L, A,' and their components. THEOREM 4.1. Let the composite operators A and L be coarse grid compatible in the sense that Ac = Lc and Aci = Lei, and let!:& be the ith composite residual obtained by using the above twolevel FAG scheme. Then fori > :o, we have rb, = 0 and (4.3) where LI and AI are the respective Schur complements in L and A. Proof. According to the Schur reduced form of the coarse composite I I equations, 1 of FAC (the coarse composite solve) yields zc satisfying (A A A 1A A A1A ) A A1 A A1/P IIC C CIIP p PI ZG,I = riIC C raIP p FrF, ( 4.4) and Also, the update of the composite residual may be written as r,H1 = r.' Lza, where zc = (zc,c, zc,I) Using along with (4.1), then the residual update may I be written in component form as rt1 = rb Lczc,a LciZG,I = rbLc(Ac1(rbAcizc,I))Lcizc,I, (4.5) 64
PAGE 72
.. A I. 1 of IL!Ar.. Although we only used the assumption that Ac = Lc and Acr != l,cr above, it is also appropriate to make the assumption for connections, i.e., Arc=Lrc. (4.8) I : Let be eigenpair of LrAi1 and set y = Ai1x, so that the equation becomes; I' "' "' The eigenvalues of LrAi1 are therefore those values of..\ for which the I ,I matrix : j Lr .\Ai + (..\1)LICL(/ LcrLrFLP,1 LFI + .\ArpA:p1API (4.9) I' I I I IS (notice that we have included ( 4.8) in the compatibility as' I sumptio:qs in1 order to obtain this latest expression). This criterion for i I determining eigenvalues will be applied in the next section. I ' I 4.2 Con:Vergence in One Dimension We now turn to an examination of the convergence of the FAC algorithm when it is applied to the solution of the finite volume discretization of II : problem (3.1)'. Since the operators A and L arising from this discretization : satisfy the co,bpatibility conditions Ac = Lc, Acr = Lcr, and Arc= Lrc of the section, it follows from Theorem 4.1 that convergence of the algorithm is governed by the relationship ( 4.3). We are therefore interested ,. in the eigenvalues of the matrix LrAi1 Since ( 4.8) is valid, it is sufficient to set the of the matrix ( 4.9) equal to zero and then solve 66
PAGE 73
for the resuidin.g values of A. Our analysis begins (Section 4.2.1) with the I case thaf composite grid has a single refinement region. Notice that ol in this case two nodes comprise the interface, so the dimensionality of the eigenvalue is equal to two. Therefore, identifying the required I [, I ., eigenvalues is: "simply" a matter of the solving a quadratic equation. For the consideration, the analysis is facilitated by the fact that I this appears, essentially, in factored form. Later, we extend this analysis by examining (in Section 4.2.3) the case of multiple refinement regions. 4.2.1 Tlie Case of a Single Refinement Region We now pfoceed with analysis of the spectrum of the operator LI ..4.[1 ,, to a finite volume discretization of the problem (3.1). With a I single refinement region on the onedimensional composite grid, and using either the centered or upwind discretization& of Chapter 3, the operator L has the str1ucture shown in Figure 4.1. Denote the components of the composit,e op:erators as follows: 0 0 Lc = [ ] Lcr = 0 0 lei 0 0 rei 0 [ 0 LIC = 0 0 he 0 0 67 0 0 0 0 Tie 0
PAGE 74
lFr 0 lPI 0 0 0 0 0 LFI = Apr= 0 0 0 0 0 TFJ 0 rpr L i [ f,F 0 0 0 l A = [ l,p 0 0 l IF1 0 0 0 TJF IP 0 0 0 I' 11 ( ;I I il 'I I !! ,.i :c 'i :z:,, I. :cl :c I :c l ;I :c : 1 :c :c :c :c :c 'I [I :c :c :c ! :c :z: :z: 'I :c :z: :c 'I :I :z: :c :c iJ :c :c : 1 :c :z: :z: :c :z: :z: :c :c :c :c :c :c I :c :c :c I! :c :c :c I. i :c :c :c ',, Fig. 4)1. Structure of a onedimensional composite operator :i for the case of a single refinement region. 68
PAGE 75
We assuz4e the Dirichlet value at the left boundary has been moved to the side of ( 4.1 ). Let the respective orders of the matrices L, R, LF, Alp be Nc1 Nc2 NF, and Np. Also, let Nc be one less than the or nodes on the global coarse grid (i.e., the order of the global coarse with the left Dirichlet condition eliminated). The coarse is::denoted by h = 1/Nc. Denote the ijth element of x1 the inverse of generic square invertible matrix X, by Z&i Also, let X;x; !' denote th:e j x j submatrix located in lower right corner of X. For example, if X is of:ordr N, then X1x1 = ;cNN and XNxN =X. With this notation, I we have L L 1 L [ hclcriNc Nc IC C CI1 1 0 0 l rrcrcrru ( 4.10) 1'L L1 L [ IIF'F!iF,. lrFTFI_iF1,Np l :: IF F Fl = TJFlFilFNp,1 T[FTFilFN N P p and ;lrpA:p1 API = I lrplPJ4p11 lrPrPiaP1,Np l rrplPJapNp,1 rrprprapNp,Np In order !to we need to obtain explicit representations for certain of the matrices L1 R1 L"F1 and A_p1 Notice that the above only involve the corner elements of these matrices. To I carry out required analysis for various finite volume discretization&, I we use a: notcltion that indicates discretizations that change character on the componeri.ts of the coarse composite and true composite grids. We use I I a fourtu'ple, ;:( z1 z2 : ;c3 z4), with elementS z1 through z4 that indicate the type 'of discretization used on the respective coarse component of the ,.' 'I coarse grid, coarse patch, coarse component of the composite 69
PAGE 76
I grid, and patch. The interface nodes are not specified here. In this I section, tJ:te type of discretization at these points will always be of the same type used' on the coarse component. (In two dimensions, it will be useful to allow at the interface of the coarse composite grid.) In one dimensioit, possibilities are considered, upwind (:z:i = U) and centered : ( :z:; = C) versi;ons of the finite volume discretization. In two dimensions, we ' will also higherorder upwinding to be used on the patch. Generally, it will be 'possible to make the notation just described less cumbersome by employing obVious abbreviations. Initially, we will consider discretizations I ,, that agre 1 e in throughout the two grids. The case we consider in detail is the discretization signified by (U, U : U) U4 which corresponds to upwinding used throughout the coarse and true composite grids. The analysis will make use of I II the generic t:qdiagonal matrix X= tri[eh, 2e+ h, e) of generic order i N. Of interest are the corner elements of x1 : :C11 z1N, ZNt and ZNN i Let Xi; cofactor of :Z:ij. Then [28), I LEMMA 4.2. The four comer elements of x1 corresponding to the i upwind discrJtization are given by (4.11) ( 4.12) and (4.13) 70
PAGE 77
.:, Proofi: X is tridiagonal, there is a threeterm linear, homogeneous ;, i differenc;: for the jth determinant, z; = DetX;x;: I' i = (2e + h)z;1 e(e + h)z;2, j = 3,4, N. The characteristic equation is ;. I' with root's 1: p,2 p,(2e +h)+ e(e+ h)= 0, P,t = e + h, P,2 = e. I, j Hence, itFfollclws from the theory of difference equations that I z; = S(e+ h);+ Tei, I where S r;r are unknown constants which are determined uniquely by II I I' I the of z1 and z2 But 1: z1 = S(e +h)+ Te = 2e + h z2 = S(e + h)2 + Te2 = 3e2 + 3eh + h2 Thus, S =(e+ h)jh, T = ejh and we h,ave ; I :I ., ,, i; i, Now since X !is tridiagonal, it follows that the cofactors satisfy Thus, _ ZN1 Zu = ZNN = = ( + h)N+I N+I ZN f f 71
PAGE 78
' ' I Also, XN1 = implies that I i Z1N = (c: + h)N+lc:N+l" Finally, x1N (c: + h)N1 implies that h(e+h)N1 ZN1 = (c: + h)N+lc:N+l" 0 More will; be shortly regarding the quality of the approximations given in the lerfima;(see the comments below and also the next section), but be l cause of the two approximations we will see that L and A are essentially I ' lower tricl:nguiar, so Det(LI ..XA1) just depends on the diagonal part of i the matri,x. L. I We m;ay X here assume the role of Ap when the coarse patch has trivial i.e., when the global coarse grid plays the role of the I I coarse grid in the twolevel version of FAC. We have just found 1: i the requited elements of the inverse ofthis operator. In this context, it may I be P?int out an apparent ambiguity in the above expressions. Notice that in these expressions, h is always directly related to the dimension, I Nc, of glpbal coarse grid. However, the role of N in the expressions II ,I is allowed to Vary. For example, when the approximations are made with I regard Ap, then N +Np in (4.114.13), and note that, necessarily, ,, I Np < Ne. Qne sees that there is an interplay between the size of c: and the of the global coarse and coarse patch grids that enters into i the valiW:ty of the approximations. In this particular case, their validity is commensurate with c: being small with respect to h (the reciprocal of Nc) and the patch being sufficiently large. Also, X rri,ay play the role of the coarse component operators L and R, 72
PAGE 79
so that 1 lNc Nc' rn 11 1 + I I Now suppose that there are m levels of refinement on the fine patch. Then I and j: = (2m)Np((e+ h/2m)Nr+t_ eNr+l)/(h/2m). i I Using efpressions, we obtain the following corollary to Lemma 4.2. COR0LLA.RY 4.3. The corner eiements of the fine patch matriz LF I 'I to the upwind discretization are given by il l :. ll; = (e+h/2m)Nr1 Fu! 2m((e+h/2m)Nr+teNr+1) 2me+h' (4.14) and h. I l "'"'0 F1,Np 22m( ( E + h/2m )Nr+t eNr+t) "'"' I I '" I [1 (4.15) (4.16) 0 The of Corollary 4.3 may also be applied to the inverse of Ap, when an composite grid is used in place of the global coarse grid, : lj by perforhnng the replacement m mk (0 k < m). I In adftio to these quantities, the pertinent quantities from the FV of Section 3.2.1 are (here (u) denotes "upwind", and (c), !1. "centered"): I ' lrc = eh (u), lcr = e (u), 1'JC = E (u), 1'CI = Eh (u), 73 eh/2 (c), e + h/2 (c), e + h/2 (c), eh/2 (c),
PAGE 80
IJF = 2me (u), 2me + h/2 (c), lFI = 2meh (u), 2meh/2 (c), TIF = 2meh (u), 2meh/2 (c), TFJ = 2me (u), 2me + h/2 (c), l[p = 2mA:e (u), 2mIce+ h/2 (c), ln = 2mA:eh (u), 2mA:eh/2 (c), TIP = 2mA:eh (u), 2mIceh/2 (c), TPI = 2mA:e (u), 2mA:e + h/2 (c), LI = [ (1 + 2mO)e + h 0 l (1+2m)e+h [ (1 + 2mA:)e + h 0 l 0 (1 + 2mA:)e + h (u), [ (1 + 2m)e 0 l [ (1 + 2mA:)e 0 ] LI: q, (1 + 2m)e AI0 (1 + 2mlc)e (c). Notice that have allowed for general refinement on the coarse patch by i letting 1 k,, m. When k = m this corresponds to letting the global coarse gn'd play the role of the coarse composite grid, and when k = 1 I ,, to using the grid having one less level of patch refinement than the true i I Now, the above values and approximations, we see for the U4 case that LI _; .\A{ = [ (1 + + h.\0((1 + 2mA:) +h) 0 J (1 + 2m) + h.\((1 + 2mA:)e +h) A simple, cal<;,ulation then shows that, irrespective of the value of k, the solutions of .\AI)= 0 are .\1 = .\2 = 1. 74
PAGE 81
I :I Next,.we the 04 case where the centered discretization is used Again, the corners of the inverse matrix are required. Pro ceeding a;s in previous case, let X= tri[eh/2, 2e, e + h/2] be of order N. i i LEMMA 4.4. The four corner elements of x1 corresponding to the centered are given by :. (e+ h/2)N(eh/2)N 1 = = (e + hj2)N+I(ehj2)N+I E + h/2 1 (4 17) h (e+ h/2)N1 :Z:lN = XNt/ ZN = ( E + h/2)N+I ( E hj2 )N+l 0, ( 4.18) and i : : h(e+h/2)N1 h :Z:N1 : xlN I ZN = ( e + h/2)N+I ( eh/2)N+I ( e + h/2)2. ( 4.19) ' The jth determinant satisfies I I 1! 1 z; : 2ez;_1 + ( eh/2)( e + h/2)z;_2 j = 3, 4, N, I with associated roots J. = E + h/2, J. = Eh/2. Hence, z; = S(e + h/2); + T(eh/2);, and in particill.ar : 1z1 = S(e + h/2) + T(eh/2) = 2e, 1z2 = S(e + h/2) 2 + T(eh/2) 2 = 3e2 + h 2 /4. I Solving S T, we find I I S = (e + h/2)/h, T = (eh/2)/h. 75
PAGE 82
Therefore, ,,, z; = (e + h/2);+1 /h{eh/2);+1 jh, and _ (e + h/2)N(eh/2)N :Cu = ZNN = (e + hj2)N+l(ehj2)N+1' Also, h(e+h/2)Nl :CtN = XNt/ ZN = ( + h/2)N+1 ( h/2)N+l' and h (e + hj2)Nl = XtN I ZN = ( + h/2)N+1 ( h/2)N+1. D Again, tHese results apply to the operator L, and also to Ap if the coarse patch has trivial refinement. We note that the above argument is incorrect for one of matrices, namely R, when it corresponds to the centered stencil righthand portion of the coarse component. What is given would be: correct for this operator if it had a Dirichlet condition at the 1 I right For the Neumann condition, we should use I [ 2e e + h/2] X2x2 = h/2 + h/2 That is, tecrullcally, the initial values for the difference equation have to I I I be modified, though we arrive at essentially the same result: now z; = (e+ so', I (e + h/2)Nt 1 ru = (e + h/2)N = e + h/2 I, For Lp haye the following corollary to Lemma 4.4. I COROLLA,.RY 4.5. The corner elements of the fine patch matriz Lp 'I corresponding: to the centered discretization are given by 76
PAGE 83
(4.21) and . h(e+h/2m+l)Npl ,..... h lFNp,l = 22';(( E + hj2m+l )Np+l ( Ehj2m+l )Np+l) ,..... (2me + h/2)2 ( 4.22) 0 A similar,.rest.tlt holds for Ap when the coarse patch has nontrivial re:fhiement, by performing the replacement m +mk. Continuing, then, I A A for the case we have Lr AAr = ![ (1 + 2m)f..\(1 + 2mk)f 0 l : 0 (1+2m)f..\(1+2mk)f I I 1 A A The solutions .of Det(LrAAr) = 0 are again A1 = A2 = 1. ; Next we turn to the (UC)2 discretization, which uses upwinding on the coarse component and centered on the patch for both the coarse and true !' I i composite grids. Using approximations from the last two cases, we have LrXAr .. = I +(A 1) [ 0 ][ h(;.::!h) l E l 2"' 77
PAGE 84
The eigenvaldes are .\1 = 1 and We consider briefly the validity of the approximations used here. On the one hitnd,'lwe have used (4.114.13) for the coarse component (upwind) I operators, which are valid with f small with respect to h. On the other I hand, the (centered) approximations on the patch come from (4.204.22). Suppose that i1on the coarse composite grid no refinement is used on the patch = 0 and NF +Np). Then the approximations will be valid with e h/2 !:and also if h continues to decrease. Notice that .\2 + 1 as I' I h + 0. h approaching zero with e fixed is at odds with I the requirement for the coarse component. There are two approaches that I avoid this1'predicament. The first approach is to use sufficient refinement on I, 'I the so that h need not be particularly small for ( 4.204.22) to I be valid (becarse the exponent N p (or N F) appearing there increases with m ). that .\2 1 with m large and k small (i.e., with sufficient : on' the fine and coarse patches). Unfortunately, taking this 1: approach rules out the use of the global coarse grid in the role of the coarse I composite grid. The second approach allows for this possiblity by altering the discretization on the coarse patch in order to weaken the requirement I I that h be'small. I i1 This leadsj us to consider the U3C discretization, where upwinding is I I ' used the coarse composite grid and on the coarse component of the composite 1 !grid, but centered differencing is used on the fine patch. As I just obseryed, lthe following approximations will be valid with e sufficiently small respect to h (for the upwind approximation on the coarse com ponent aU:d onl the coarse patch) and with sufficient refinement on the fine 78
PAGE 85
patch (foi theicentered approximation). Although our motivation here is to be able use' a coarse patch without refinement (k = m), we have carried out the allowing arbitrary refinement on this region (0 < k :5 m ). I The is L1..\A1 = ..\ 2me+ h 2 2mf + h/2 Notice that the desired effect has been achieved in that now ..\2 approaches I unity while o:hly requiring sufficient refinement on the fine patch. The I, I expressions for the eigenvalues are independent of the refinement used on I ' the pat,ch and so they will be valid, in particular, when the global coarse grid is used. We also remark that this requirement (of having suffi cient on the fine patch) is a natural one with respect to the use i of the cettterefl type of discretization for advectiondominated problems. Although ,the expression for ..\2 is, only valid in certain parameter ranges of the oned.imen:sional problem, we will see that the quantity ..\2 1 supplies I a fairly acurate approximation to the actual convergence factor even in I ,, two dimensions. I Collect,ing 1the results from our examination of the various cases con sidered a 'Hove establishes the following theorem and corollary which apply I 1 to the problem. 79
PAGE 86
THEOREM 4.6. Suppose that the ezact values of the corner elements in I equations (4.114.22} are replaced by there corresponding approzimations. I, I Then, in .the U4 case, the C4 case, and in the limit as m oo in the U3C lj I case we fhe following. Both of IL1AI1 are zero. Therefore, two steps of the iteration drives the residual at the interface to zero. This matriz is in fact strictly I lower the "left" and "right" components of the interface residual J.: vanish orJ, the ',first and second steps, respectively, of the iteration. D r I CORQLLARY 4.7. Suppose the twolevel FAG scheme is written in the form ri+11'= (ZLM1}r' for some composite preconditioner M. Then, I with the of Theorem 4.6, this iteration matriz is nilpotent: (ILMT1 } 3 0, i.e., the twolevel scheme converges in three iterations 'I and, thusj ma'y be considered a direct method. I Proof corollary follows from the vanishing of rh and r}.. established I I in Theorem 4.1 fori 2:: 1, the resulting expression (4.3) for the updating of I' I the <;omponent of the residual, and Theorem 4.6. D It should :be that the results of the last theorem (and its corollary) are in4ependent of the choice of coarse composite grid, i.e., they I I hold with:m .:k levels of refinement on the coarse patch for all values of k, 0 k < In particular, this motivates the use of the the global coarse grid in capacity. 4.2.2 of Approximations The purpose of this section is to assess the numerical quality of the : approximation::s for elements of inverse matrices given in the previous sec80
PAGE 87
tion. We iwill.show that the values given as approximations are valid with regard to: appropriate floatingpoint computations. However, their validity is subjeci to z!estrictions on the values of the parameters e, h, and m (the : ', order of patch refinement). It turns out that these restrictions allow for a large prac;:tical range of these parameters. I We begin irith the U4 method, concentrating on the inverses of Ap and L I LF. Here, the, global coarse grid is used in the role of the coarse composite grid. Let X denote either Ap or LF. In the last section, the corner elements ' I of xl were approximated as follows: and iu 1/(2me +h), ZNN 1/(2me +h), ZtN 0, where m l.. 0 for Ap, and m = mr (the additional refinement levels on the I fine patcH) for: LF. In Tables 4.14.2 we give, for various values of e, the I four pertihent.values of the inverse matrices, i.e., the 2 x 2 matrix [ l To obtain of xt, the tridiagonal matrix was reduced to upper triangular; for111 by Gaussian eliiD.ination (as in [45]), and backsolves were performed, with the canonical basis vectors e1 and eN as righthand sides. I . The coars'e grid meshwidth is h = 1/64, the coarse patch has 16 nodes, I I and there:are"mr = 4 levels of additional refinement on the fine patch. Computations:' were done in single precision on a Sequent Balance. The are valid for e sufficiently small with respect to h, i.e., as e + 0 wii,h h::fixed. However, e need not be particularly small, as seen 81
PAGE 88
in Table 4.1, for the approximations to be fairly accurate. Furthermore, comparing the second table with the first, one also sees an improvement in the approximations with f fixed and h refined. Although h decreases, the I effect of this offset by an increase in the exponent N (the order of the I matrix inverted) from Np to NF. ' :f I Table 4.1 Computed elements of coarse patch matrix inverse for the upwind case with h = 1/64 and 16 nodes on the coarse patch. Ap1 1/(e +h) h/(e + h)2 .Q001' 6.359301 101 1.757762. 1027 6.359301 101 6.318860. 101 6.359301 101 6.318860. 101 I I .901 6.015038 101 7.627731 101& 6.015038 101 5.653231 101 6.015038 101 5.653231 101 ; I i I .,01 3.902435 101 1.158410 104 3.902439 101 I 2.379538. 101 3.902435. 101 2.379536. 101 !. I Similaf resUlts are presented for (UC)2 The approximations from the i. previous are: I and Zn 1/(2mf + h/2), ZNN 1/(2mf + h/2), Z1N o, 82
PAGE 89
Here e = 1.001 fixed, the patch has length 1/4, and his allowed to increase. The (4.174.22) predict the deterioration ofthe approximations I as this Table 4.2 Cdmputed elements of fine patch matrix inverse for the upwind case with h = 1/64, 16 nodes on the coarse patch, and m,. = 4. f. Lp.1 1/(2me +h) h/(2me + h)2 I I .0001 5.805515 101 0.000000 10 5.805516 101 I I 5.266088 101 5.805515 101 5.266252 101 !. : I .I I 3.162055 101 0,000000 10 3.162055 101 .001 : 1.562284 101 3.162055 101 1.562280. 101 I i :5.693950. 10 1.184763. 1010 5.693950 10 .01 I .065526 101 5.69394 7 10 5.065792. 101 I I Tables! 4.3.,.4.4 below confirm this, with the deterioration being more I I for Ap. This again motivates the use of the upwind scheme on. the in U30. 83
PAGE 90
Table 4.3 I Computed elements of coarse patch matrix inverse for the centered case with e = 103 I, h A_pt 1/(e + h/2) h/(e + h/2)2 I 1/512 5.059288 102 0.000000 10 5.059288 102 4.999298. 102 5.059288 102 4.999296 102 1/256 3.386243 102 7.934145. 10t2 3.386243 102 4.479160. 102 3.386243 102 4.479158. 102 I 1/l28 2.041807. 102 6.061421 1ot 2.038217 102 I 3.247696 102 2.041807 102 3.245568 102 I I Table 4.4 Co:mputed elements of fine patch matrix inverse for the centered case with e = 103 and mr = 1. I I 1/(2e + h/2) h L;t I I h/(2e + h/2)2 : I 3.359580 102 0.000000 10 3.359580 102 1/512 2.204448 102 3.359580 102 2.204449 102 I I I. I I I 2.529644 102 0.000000. 10 2.592644 102 1/2'56 2.499649 102 2.529644 102 2.499648 102 II I I : 1.693122. 102 i 3.808327 1011 1.693122 102 1/1128 :: 2.239580 102 1.693122. 102 2.239579 102 I 84
PAGE 91
' 4.2.3 Multiple Refinement Regions I Next we extend the analysis of Section 4.2.1 to the case of multiple interior regions in one dimension. In general, there may be various refinement regions, but we suppose that, without loss of generality, there are 'but .two; the significant difference between the case of one region and many is tp.e existence of coarse regions between regions of refinement, not just the two ends of the grid. Allowing only two refinements will also allow us to avoid significantly complicating our notation. With this convention, L has the structure shown in Figure 4.2. r:1: Zi :1: :1: :1: : :1: :1: :1: :1: lei I :1: :1: Ter :1: :1: lei :1: :1: Tei I t t t t t t t t t I I I :1: :1: :1: i : :1: :1: :1: :1: :1: :1: I' lre t :1: I Tie t :1: lre :1: :1: TJe :1: :1: '' 4.2. Structure of a onedimensional composite operator 1 in the case of multiple refinement regions. 85
PAGE 92
The grid still has but a single fine component, although it is now the muon of the various refinement regions. This is also true of the coarse and interface components. Accordingly, the composite operator still has structure, though its main nonzero blocks may now be further decomposed. For example, LF = [ LFl LF"J l and LI = [ Lil l where ea<;;h Lii is a 2 X 2 matrix associated with the ith refinement region. Notice that a partitioning of LI with respect to these blocks induces a sim: ilar partition.ing of the matrices LIF and LFI into block diagonal matrices. Thus, for example, we may write L LlL [ LIFl IF1 F FI = L {The leadj.ng 9iagonal blocks here correspond to the blocks in the above ,' ' figure with the "t" elements.) It follows that LILIFLP,1 LFI is a block diagonal Dl.atrix with diagonal blocks Lii LIFiLP,i1 LFii Furthermore, each one of the:se bi,ocks has precisely the same form as LILIFLP.1 LFI had in the case 6f a refinement region, although the dimensions may vary. ,I Similar hold for the global coarse operator and AIAIPAp1 APJ. ,, As for:the coarsetointerface connections, let 0 0 I 0 0 hclcifflN.c Nc "J' "J 0 86
PAGE 93
' It may to the reader to compare this equation with ( 4.10). Recalling thb antlysis of the values of elements of matrix inverses, we see that I the quantity if't,Na2 which is the upper right corner element of M1 may he set to allowing us to use instead LrcL(/ Lei = [ l}clciiNc,.No, 0 0 0 rrcrc1mu 0 0 I 0 lrcrciihNc2,t hclciihNc Nc 0 I 21 2 I, :o 0 0 rrcrc1ru With sinfplification, the matrix L1AA1 is block lower triangular. I Furthermore, ach of its diagonal blocks has precisely the form that L1 AAr had in the case of a single refinement region. It follows that the analysis performeq that case may now he applied to each one of these blocks. The result is that, for each of the various cases considered previously, the same values m'ay he given now to the eigenvalues of the matrix IL1AI1 I ., the only change being in their multiplicities according to the actual number i I of regions:of refinement. I, ll 4.3 in Two Dimensions ' In thisi section we study the convergence of the twolevel version of FAC j I ,, applied to\:the of the finite volume discretization of problem (3.2). Because of the complexity of the structure of the composite grid equations I involving refinements, we have not yet extended the analysis of I the to two dimensions. Therefore, the examination here is I . based on computational evidence rather than on analysis. Guided by the results of ?revious section, we attempt to determine to what, extent I the results' for onedimensional problem carry over to two dimensions. by 'the findings of that section, we concentrate our study on I I the behavi,or of an algorithm that emphasizes the use of the global coarse 87
PAGE 94
I i grid. The numerical method used here to solve the composite grid equaI ,, tions is tlie twolevel FAC scheme with the global coarse grid as the coarse composite Subproblems (corresponding to the fine patch and global ' coarse are solved "exactly" (i.e., iteratively with a tolerance on the relati\re of 108). This has been done by performing, in the appropria;te context, the necessary number of relaxations or multigrid V cycles further details see Section 5.3.1). ', The fiFst result emphasizes the effect that initial guess has on the i ;j behavior. Consider the solution of the upwind/upwind dis' i cretizatio of:J;>roblem 3.1 of Section 3.3.1 by FAC using the upwind/upwind discretization on the global coarse grid (i.e., the U4 technique). The prob,, i lem parameters here correspond to those used to generate Tables 3.13.2. Table 4.5 presents the sequences ofrelative composite residuals correspond!' ing to twq solutions of the discrete problem. The two solves differ in the choice of the composite guess, The zero vector and a random I vector wit;\1 between zero and one were used. The results appear, respectivety, the left and middle columns of the table. Recall that ll.rill2 here actually the current interface component of the residual, 'I since 0 fori 1 (see Section 4.1). Notice that there is a sig!, ni:ficant difference in the rates at which these interface residuals approach r zero in thl solves. This difference is due to the convergence behavior of the at the nodes that correspond to the tangential (north and south) bolndades of the interface. At the inand outflow (west and east, respectively) B1oundaries, the convergence behavior is like that predicted by the of the onedimensional problemin both solves the resid. t uals at the interface converge rapidly to zero. However, at the tangential I' 88
PAGE 95
. Table 4.5 Relative residuals as functions of FAG iterations using the U4 method in two dimensions. Upwind/Upwind : llr.i ll2/ llr.0 ll2 I :I Modi Jied Stencil '!!0 = 0 '!!0 Random '!!0 Random 1, 1.825 4 6.406 2 5.592 2 3.355 6 1.411. 10 2 1.108. 10 3 a: 6.803 10 "7 6.208 3 9.313. 10 6 4 2.742 3 8.291 toe 5: 1.212 3 6.013. 10 '7 6. 5.358. 104 7: 2.370. 10 4 8 1.049 4 9' 4.639 6 I' 19 2.053 6 I : convergence rate is generally much slower, as in the middle I column Ta*e 4.5. The reason that this effect does not appear initially I in the left is that, with the initial guess of zero, very accurate solul tion areiproduced at the tangential boundaries as a result of the first I I global coarse solve. Thus, the residuals at these nodes are extremely small j I! at the first FAC iteration and, therefore, throughout the ensuing II iterations:\' We, have developed a remedy for the problem of slow conver,1 ,, gence at the tangential boundaries that involves modifying the definition I of the cparse stencil at these nodes. The right column of Table 4.5 L ;I gives the tesults for a solve using this modification. The actual changes to i the stencilj.,are by requiring that interfacetointerface connections I for the gld,bal operator agree with those of the composite operator, i.e., AI nodes lying on the tangential boundaries. As a practical .\, :.. ,, 89
PAGE 96
matter, these' changes do not adversely affect the efficiency of the algo rithm's iinpletnentation. On the contrary, they can be beneficial, making the global matrix more diagonally dominant, which facilitates the solution of equations by an iterative method. To these modifications, consider the case where the only interface component is a tangential one, i.e., the twodimensional composite I ', grid with the domain partitioned into two horizontal strips, as depicted in Figure 4.3. Also, let f __. 0. 0 0 0 0 I ; I ,. : 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0' 0 0 0 0 0 0 0 0 0 0 0 0 0 0 : Fig. 4.3. Composite grid with coarse ( o ), fine ( ), interface ( 0) and slave ( o) nodes. 0 0 0 0 Then, usirlg the stencils of Section 3.2.3, one sees that only the diagon:al blocks of composite operator L in ( 4.1) are nonzero. The same is true for A. Scliur complements in these respective operators then become LI and A}, atid the equation for the updating of the interface residual simplifies i+l (I L A1) i ri = I I ri. I This clarifies the reason for the modification. Notice that without the 90
PAGE 97
where 8 = + 2,!+1 Thus, r}+l = (18) r}, giving co:hvergence factors between one quarter and one half. This agrees roughly With the observed convergence rates. We remark that as E + 1 the becomes less beneficial, and even detrimental, but the way that this ),iappens depends on both E and the amount of patch refinement. 4.6:4.7 show convergence rates as functions of E and fine patch refinemen:t for,'the U4 method applied to the upwind/upwind discretization of Problem 3.1. The coarse component mesh width is he = 1/32, and m I' 'I gives the number of additional levels of refinement on the fine patch. We I I I : I 0 I I : .,, '&o It: 2: 3' !' 4 I. 5! 6 I, 7, 8'' 9!. 10 Table 4.6 FAG convergence as a function of E for the U4 method. Upwind/Upwind m=l llrill2/llr0ll2 Unmodified Modified Modified Modified :Stencil Stencil Stencil Stencil e = 102 E = 103 E = 104 E = 106 10 2 6.824 2 8.160 2 8.403. 10 2 4.287 3 5.606 4 1.654. 10 4 2.071 6 10 4 1.737 6 4.138 "7 1.790"7 5::855 105 8.934 7 8.080 6 1:143 6 . 10 7 I 91
PAGE 98
used a initial guess throughout in these tests and, in some cases, the modification to the global coarse stencil mentioned above. We found that the gained by using this modification is lost as f increases, and use of the unmodified stencil can give better performance. For a given E, the appearing in Tables 4.64.7 correspond to whichever Table 4.7 FAO convergence as a function of E for the U4 method. Upwind/Upwind m=4 l llrill2/llr0ll2 I I Urtmodified Unmodified Modified Modified Stencil Stencil Stencil Stendl 1,. f = 102 f = 103 f = 104 f = 106 1 1}691. 101 7.516. 102 5.592 .lQ2 5.738 .lQ2 2 1.232 2 1.506 10'2 1.108 .lQ3 4.108 .lQ5 3! Ll35 3 3.523 3 9.313 6 4.659 7 4 1.:192 104 1.171 103 8.291 .lQ6 5' L259 5 4.012 4 6.013 '7 6i' L3216 1.384. 104 7 .lQ7 4.786. 105 ' 8i 1.657 .lQ5 91 5.769. 10 6 1 10: 2.014 6 I strategy showed better performance. Thus, the leftmost column in Table I i 4.6 and the two leftmost columns of Table 4. 7, correspond to the use of I the unmodified stencil, while the other results are for the modified stencil. I I Notice that the convergence rates generally improve as E decreases, which is what we as a result of the analysis of Section 4.2. The analysis also in'depedence of the rates on the level of patch refinement. The ,, 92
PAGE 99
results here show a mild degradation as the level of refinement is increased. d II One possible to offset this effect (i.e., to make the need for refinement I' I less critical) to use a more accurate upwind scheme on the fine patch. Tables 4.8,4.9 show results of a similar set of experiments with stan,, I dard upwinding replaced by higherorder upwinding on the fine patch, i.e., the (U H)2 (here, H denotes the use of the higherorder up wind on a composite grid component) is used to solve the upwindjijup grid problem. This FAC strategy was not analyzed I I, I i 1: i', 3' I'' 4, 5' 6,, r: 8: 9: lQ Table 4.8 FAG convergence as a function of e for the ( U H)2 method. Upwind/Hup m=1 ,, Unmodified Modified Modified Modified 1 Stencil Stencil Stencil Stencil = 102 = 103 = 104 = 10& 6;.016 2 7.415. 10 2 9.465. 10 2 9.817 2 41.376. 103 1.081. 103 9.657. 10& 1.051. 10& 4,.271 104 3.167. 106 2.522 7 1.625 7 6 9.711 '1 6l654 6 . 10 "1 I in the section, but the experiments here indicate that its con I' vergence B'ehaVior strongly resembles that of U4 This is evidenced by I 4.84.9with Tables 4.64.7, respectively (the same strat I 93
PAGE 100
' ii egy as above, modifying the interface stencil, has been used here). When I combined. with a comparison of the accuracy of the two upwind discretiza.I tions (see Section 3.3.1), the results indicate that there is little advantage to I using the: upJlnd/upwind discretization. The more accurate upwind/Hup method be used without sacrificing efficiency. In this study, we used strategies:. basbd on block Jacobi relaxation to solve the patch systems. With the l.upWind discretization, this means performing backsubstitutions ... associated wit;h the inversion of tridiagonal blocks at each relaxation step. With the ]!norf accurate method, the increase in work corresponds to using blocks having:. one additional diagonal. These relaxationbased strategies will be more fully in Chapter 5. [ !, I I !' : 1 I I I' I .I ,, 1 .. I, 2: 3il I 4: I 5!': 6r. 71\ I, 8 !'I 9 i: 10;' ., Table 4.9 FA C convergence as. a function of e for the ( U H)2 method. Upwind/ Hup m=4 I llr'll2/llr0ll2 'I UiJ,modified Unmodified Modified Modified i Stencil Stencil Stencil Stencil 1e = 102 e = 103 e = 104 e = 106 1.'674. 10 1 I 7.082 2 4.933. 10 2 5.026 2 L;389 2 1.352. 10 2 1.158. 10 3 7.143 o L286 3 3.247. 10 3 8.733 6 3.364 7 10 4 1.086. 10 3 7.515 6 l.'fl87 6 3.718. 104 5.304 7 10 6 1.281. 10 4 2.377 '( 4.424. 10 6 I 1.531. 10 6 5.326 6 1.852. 10 8 ' I I Table 4.10 the results for the 04 technique applied to the solu1 I 94
PAGE 101
I! tion of the centered discretization as a function of E. As expected, ,' i: for this method, the convergence rates degrade with decreasing E. Gener: lj ally, for that involve the centered stencil on the patch, no ' gained from the use of a modified stencil. Accordingly, the remainder of the results of this section pertain to the unmodified stencil. 1 2 3 4 5 6 1,, ,. 7 Table 4.10 FAG convergence as a function off for the 04 method. Centered/Centered m= 1 llr'll2/llr0ll2 E = 101 E = 102 E = 103 6.846 2 6.610. 10 2 1.233. 10 1 3.281 3 4.431 .lQ3 3.138. 10 "2 1.847. 10 4 3.916 4 7.962 3 1.061 .lQ5 3.906 .lQ5 2.162 3 6.290 7 4.15.0 6 6.117 4 4.6617 1.789 4 5.343 6 We results in Table 4.11 for the U3C technique applied to the upwind/ discretization with various f and m. The convergence 1' rates imp(ove em increases, as predicted by the analysis ofthe previous I section. yY.e however, a limit to the applicability of this result in that I f must not be allowed to become too large. If f is allowed to increase with I hand m a degradation of the convergence rates occurs. This fact is I in accordance with the requirement, in the analysis of this method, that j: f remain s.mall enough to ensure the validity of the approximations used ; I with respe:tt to the upwind stencil on the coarse components. I Finally?. we that the predicted convergence factor for the residuals, namely to the U3C method in Section 4.2.1 provides 95
PAGE 102
I ,. I, 'j a 'accurate (although somewhat pessimistic) approximation to j, the true rate.: The results in Table 4.12 compare this quantity for the four columns Ta:ble 4.11 with the corresponding average Here, I the value of i corresponds to the final iterate in a given column of the table. I I. ' I I' I' i 'i I. ;1 2 a 4 . 5 7 I I I I Table 4.11 FA a convergence as a function of E and m for the U3a method. u: pwind/ a entered :I II xi liz/ 11!:0 liz ' m=l m=4 ,, I .,. = 10z e = 103 E = 10z E = 103 a;19o; 10 z 1.149. 10 '1 1.784. 10 1 1.011. 10 1 10z 8.888 . 10z 1.577. 10z 3.178 z 4Jl34 3 7.106. 10z 1.512 3 1.343 "2 I 1:515. 10 3 5.875 2 1.627. 10 4 5.692 3 5.17 44 104 4.978. 10z 1.782. 106 2.489 3 10 4 4.298 z 2.008 ti 1.103. 10 3 8;:735 106 3.759. 10z 2.787 7 4.958 4 3i470. 10 6 3.318. 10 z 2.253. 10 4 ', Table 4.12 d,omparison of actual and rredicted convergence : factors for the u a method. L I I : m=l m=4 A _.:, z"'t:+h zl 2"'t:';th/2 1: E = 10z E = 103 E = 10z E = 103 IllA'21 !, I 2. 77 101 6.53 101 1.16 101 3.50 101 : 96
PAGE 103
5. FAC Algorithms .. In this chapter we consider modifying the FAC algorithm that is dei . scribed in 4 to obtain a more efficient technique. Multigrid plays I I a significant r?le here, and we assume a familiarity with the basic multigrid su,ch as the Vcycle and the full multigrid scheme appearing, I for example, in [25), [7], and [8]. I ; Section 5.1, we examine the implications of viewing FAC as a precpnditioned iterative method. This provides an opportunity to compare the with domain decomposition methods in Section I : 5.2. In 5.3 we consider a truly multilevel approach to using FAC. The of the method is given in Section 5.3.2. It is prefaced by a discussioO: and analysis in Section 5.3.1 of the role the block Jacobi method I lr plays in this ultimate algorithm. ' 5.1 FACias We attemp't to identify the matrix M for which the FAC iteration has i the form I I I (5.1) Thus, we tare to view FAC as an unweighted preconditioned Equivalently, we may write = + where i satisfi.esf' Mi'' = r_i. Recall that the process by which FAC obtains the : may be viewed as the following twostep procedure. First, the !' interface c()mponent is solved for: !. A"1 i zr = 1 rr. 97
PAGE 104
Then, coarse and fine components are updated according to the orig' inal equations. This process is represented by the system of equationf!, (5.2) where and Using the1 fact. that = O, it can be immediately seen that M is the ,, coefficient matrix in (5.2). ' One to improving the convergence rate of FAC is to acceler,, I ate (5.1) by modifying it to form a polynomial iterative method such as : .I a Richarson's iteration [39), a Chebychev iteration [31), or a Krylov 'subspace iteration [19). Each step of a preconditioned version of each requires that, for a given vector r', say, a precon I I I ditioning can be applied to yield a vector 'I!'. This precedure is I I formally as i I In the context, this procedure is contained in Steps 1 and 2 of the FAC (see Section 4.1) where Azc = rc is solved, yielding components za,c and zc,I, then LFzF = r}LFIZG,I is solved and 'I!' is i I taken to be 98
PAGE 105
I I 5.2 A Comparison of FAC and Domhln Decomposition 1, ;I It is the intent of this section to clarify the connection between FAC and the Schur complement variety of domain decomposition algorithms. I Having FAC written in the form of a preconditioned iterative method will facilitate this ifiscussion. The domain decomposition algorithms are similar I in form t6 FAC, and in some cases may coincide with it [6], [21], [22], (44) so the insight into FAC convergence developed here may apply to a broader class of It is to think of the algorithm, and even the discretization strat egy, prese:O.ted: in this study in terms of domain decomposition techniques, : I since the 'composite grid naturally induces such a decomposition. At the same there exists a strong distinction, in terms of motivation, between I II the methods c6nsidered here and what is normally considered domain deThe latter often has the intent of partitioning the domain into : I many subdomains in order to achieve algorithms that are efficient with re spect to p:arallel computer architectures, while the motivation behind the domain that occurs as a result of using a composite grid is I better moaelin:g of the solution of the partial differential equation. At the < 'I I same a second type of domain decomposition, relating to parallelism, is to the methods of this study. This is the natural decomposition I of the domain induced by the advectiondominated problem written in the I 1. I form i + 1Lz(:c,y) = f(:c,y), and that cc;ures'ponds to integrating the differential equation along its grid aligned characteristics when e = 0. We note that it is not necessary that : e be identkally: zero for this parallelism to be taken advantage of. In this 99
PAGE 106
' study, we have used the fact that e is small to good effect by employing I line as a multigrid smoother. To begin the comparison, recall the Schur reduced form of the system Lz = !:= [ :F l [ ][ ] 0 0 Lr zr rr As noted1:in S.ection 4.1, according to this form, one approach to solving the system is to solve exactly the system in the Schur complement, and then the back substitutions associated with the remaining equations! to the solution. When this approach is used as the I i basis decomposition algorithm, emphasis is placed on solving efficiently1 the Schur complement system typically, a Krylov subspace j I, iteration is used, and much of the research in this area is is devoted to a good preconditioner for the Schur complement. Now, each I' step of iteration requires the application of the preconditioning as well as an!: application of the Schur complement, I (5.3) I the parallelism of the domain decomposition methods deriving from the q, independehce of the inversions required here. One of the approach just described is that the inversions :, I, required in must be performed accurately just so that the coefficient matrix L1 of the Schur complement system is represented accurately. As ; noted in [27) and [35), this disadvantage is especially significant when an iterative is used to invert the operators Lc and LF only approx imately. A better strategy is to incorporate the preconditioning for the 100
PAGE 107
Schur system along with the subregion operators into a. com.. posite for the original equations. With this approach, the matrix only enter into the application of the preconditioning. An form for the composite preconditioner is M = [ :F l WIC WrF Mr (5.4) where Mt = fr + WICL(/ALcr + WrFL.F1 LFI and Pr is a. preconditioner for the Schur fomplement Lr. When M has this form, (irrespective of the particular, choice of WIC and WrF ), afier the first iteration the coarse and fine of the residual are zero and the interface components in (5.1) may be written as Indeed, 2, lc (0 0 (J LA PA l)lc1 l)T ri = , I I ri Many standard domain decomposition algorithms are based on a. precon1 ditioning tf form [12], [27), as is FAC. In the latter, (5.4) holds with PI = Ar, the Schur complement in the global coarse operator A. This is one of the distinguishing features of FAC: it preconditions the Schur com plement irl,the !,original system with the Schur complement corresponding to : a coarse grid discretization, and this preconditioning is accomplished im plicitly by :way of solving the coarse grid problem. The reason that this ' appears implicitly in FAC is that the algorithm was born out of an a.tteJ:l!.pt to economize multigrid when applied to composite grid problems, and so it is normally seen as a multigrid variation, rather than :.di :1 as a tianer. ,, , I 101
PAGE 108
It is :Worthwhile to say a few words about the economy that FAC : I achieves in solving this coarse grid problem. First, an effective precondi,, I tioner complement is provided. At the same time, the coarse I component bqundary value problem that is implicit in M's first block row is solved. Finally, in modifications of the algorithm that only approximately I i invert LF,; the patch component of the solution of the coarse grid prob' I' ' lem initial guess for the solution of the fine patch boundary problem (M's.lsecond block row). I It is to note that a similar algorithm (BEPS) exists explic1, I itly as preconditioning in the work of Bramble et al. [6], [21]. That version I I I of the algoritijm is founded on Schur complement domain decomposition I I, where it i,s to view the solution process in terms of a preconI' ditioner. :The' particular domain decomposition method on which their 1 .1 algorithm:,is based is that of [5], which is similar to the method described I ,I in [4]. Both these latter methods use a preconditioner in the form of a Schur however, it is the Schur complement corresponding to ,: I a special discr,etization on one of the subregions. The preconditioning is ,, I I, implemeni'ed by solving a boundary value problem (with the interface as a I portion of, the ,boundary) on one of these regions. The modified technique I of [21], hoivever, has, after its first step, precisely the same form (compare 1. ,I [21, 3] with the previous section here) as an accelerated version (by I use of a Richardson's iteration) of the twolevel FAC algorithm I' ,. (see, also, the of the two methods in [23]) ,, I The difference in the two algorithms is the choice of a starting point. ,1 The BEPS, algqrithm begins by solving a fine patch boundary value problem with a:D. initial guess of zero. The FAC scheme described herein begins I 1 by solving: the 'global grid version of the problem and uses the resulting : I 102
PAGE 109
interface as data for the fine patch boundary value problem. Al: I though of this difference is minor, the second approach I I gives a conceptual view of the algorithm. In the context of the I advectionLdonllnated advectiondiffusion equations the convention of starting on the. grid is quite natural. Suppose that line relaxation is used I to solve tl:i.e coarse grid problem Then, this FAC step corresponds ,' II to integrating,lacross the domain along the characteristics of an (almost) hyperbolit This result is then used as an initial guess for finding ,1 I the solutibn on a subregion where the structure of the problem is more II complex. fn this sense, the FAC algorithm is close in spirit to domain deas1 carried out in [9], [15], [37], where the domain is partitioned I and the solution technique is chosen in accordance with the behavior of the .;1 i solution by an asymptotic analysis. I I 5.3 Modiftcations Based on Multigrid The of viewing FAC as a. preconditioned iterative method and then by a polynomial method seems pragmatic enough I II (and to algorithms being used in practice). However, it has the Jrawback of requiring "many" FAC iterations to solve the ,I discrete exactly, particularly as the effectiveness of the precondii. ,J tioning in certain parameter ranges of the discretized problem. I Another to using FAC as it has been described so far is that it ,I requires solution of coarse grid and fine patch subproblems. For I I' ,, these reasons, we have taken a fundamentally different approach to uti: .: lizing FAC: effiCiently. It differs from the ones described earlier in that it I does not have its ultimate goal the exact solution of the composite grid equations. its goal is to solve these equations only to the extent 103
PAGE 110
that is by the accuracy of the composite grid discretization. I ' I The description of this method is facilitated by recalling the function I space setting hf Section 3.2.2. Let cha be a known bound (in some norm, II II) on \the discretization error of the composite grid equations, with c Jf h, and let E be the approximation to u obtained by solving exactly the composite grid equations. The algorithm to be described! is on a methodology that has as its goal the efficient computation approximate solution uh satisfying (5.5) I ;, This turns out to be very satisfactory in that it produces a solution that is comparable in accuracy to the exact composite grid solution, I while requiring few FAC iterations and minimal effort in solving subprob' I I lems. Based on the multilevel technique of nested iteration, it is a full multigridlike algorithm that uses the twolevel version of FAC as a weak algebraic on the sequence of composite grids lying between the global coarse grid and the true composite grid. We preface the more thorough descriptiori of this algorithm with a discussion of the role that relaxation : plays in o4r of the FAC algorithms used in this study. 5.3.1 ThEf' Significance of Line Relaxation A practical! approach to the implementation of FAC in its twolevel form is to utilize some variation on the theme of relaxation to solve exactly the global and fine patch boundary value problems that arise. We I have used itself to solve quickly, in certain cases, the global coarse systems' when FAC was applied to the test problems of Sections 3.3.1 and 4.3. Generally, relaxation itself will not be sufficient because of I I the deteriorati6n in its convergence rate on relatively fine grids after a few ,, 104
PAGE 111
II (seC( [8] and the analysis below). Convergence rates can also be slow on coarse grids when e is large. In both cases it is possible to gieatly the efficiency of relaxation by incorporating it into a i I multigrid:'algorithm. For the tests of Sections 3.3.1 and 4.3, we solved fine patch and global coarse grid problems involving relatively large I ' values of: e "exactly" by applying the necessary number of multigrid VI I ,I cycles required to reduce the size of the initial residual to below a small I II tolerance.: This approach to the twolevel algorithm is satisfactory in that a modest 'nurriber of Vcycles are required per FAC iteration, and each V il cycle can ibe performed at a cost of just a few relaxations. In addition, the alg?rithms used are robust in that they apply in general to the I various discretizations used. I The the efficiency and robustness ofthese relaxationbased all gorithms i:s that line relaxation is ideally suited to the advectiondominated problem the\form (2.2); when applied to the discrete problem, it provides an solyer in certain parameter ranges and an effective (multigrid) smoother in o(hers. Consid'er the nodes of either the global coarse grid or the fine patch or dered in lexicographic fashion, with the leading dimension corresponding I to the :cdirection. If the finite volume discretization is applied to prob,. lem (2.2) t!ps uniform grid, the result is a five(or six) point matrix stencil a block tridiagonal form, the main diagonal blocks themI selves being tridiagonal matrices in the case of the upwind and centered I discretizations ;(in the case of a higherorder upwind discretization these blocks nonzero main diagonals). By line relaxation associated I 1, with the sjstem Lz = r is meant the (possibly weighted) simple iteration I ,I based on the splitting L = B N, where B is the matrix that is equal to I 105
PAGE 112
l I I I I I I I II I Lon its blocks and zero elsewhere. This iteration may be written I 'I I as 'I where ri = rLzi. The ith etror7 [ ei = z zi, associated with the iteration is given by 'I I .. \ ei =(IwB1 L)ie0 i' :1 where e0 is error associated with an initial guess z0 Notice that this ,. '\ I I may also. as I I : where the::ith:idegree polynomial. Pi(z) is the product of i linear factors, I ,I (1wz). the method may be viewed as an acceleration (by propitious 1: J choice of w.) simple iteration for which p;(z) has linear factors, (1z). I ,I Each of the polynomial Pi may be viewed as a mapping, Pi : I: + I' I p;(E), of t4e sJ,t of eigenvalues, I:, of B1 L via Pi(E) = {p;(.\) :A E I:}. In I' I the remairl'der bf this section, we will examine the spectrum of B1 L under I : these for the upwind and centered discretization& of problem I ,I I' ', I .i (2.2). I' We bemn by generalizing problem (2.2) to 1: ;! I . Llu(z,y) + Puz(z,y) = /(z,y), (5.6) I 1: where P is:[a constant. This generalization is done purely as a matter of convenince in interpreting the results of the analysis. At times it is de. .I I sirable to se tl:lie effect of letting become large, and this can be facilitated II II I .I by setting = l and letting p+ 0. We note that the stencils for uniform ,I I 106
PAGE 113
II I grids in:1Section 3.2.3 can be made to correspond to this new equation : 1 by perforlning the replacement h {3h. Regions containing the eigenval1 ues for centered differences applied to this equation are given in [14]. We are also interested in the upwind case and will require specific knowledge 1:' of individual eigenvalues and their associated eigenvectors. These eigen1 I be found for both discretizations using the techniques of [20] and [29]. : Assu4e, for convenience, a square mesh with spacing h = in each I. ,, direction \(the: following analysis being easily generalized to rectangular meshes). :A1so1 denote the halfgrid Reynolds number by 'Y = In what follows, :E denotes the spectrum of B1 L, for a fixed 'Y At times the I dependence of the spectrum on 'Y will be acknowledged by writing E('Y) I' A portion; of t;he spectrum will frequently be located in a disk, that does I not contain t!J.e origin, in the complex plane. We will be interested in I how well the iteration damps eigenvectors of B1 L that correspond to its eigenvaluJs in :the region, and in the best way to choose w in defining Pi i :1 For eigenpair, (..\, v), we have and we wruld I 1 w..\ I to be small. This quantity is bounded by IIPilloo over the disk. It is well known [31], [43] that, for such a region, to 1'minimize IIPilloo over all polynomials of degree at most i, one should the polynomial (1wz)', where w is the reciprocal of the disk's centfr. IIPIIoo =I w I r, where r is the radius of the disk. We will refer this quantity, when it corresponds to the optimal choice of w, I as the spectralfadius of the region and denote it by Popt Also, we will have occasion td clai;m that the iteration is convergent, by which will be meant I I 107
PAGE 114
that it to the exact solution for all choices of z0 ' i' II th1e elements of :E by 1 j, k N. The eigenvector of n1 L with Ale; is given by the tensor product VIe v; between h i II h 1 t e two vectors, wlt e ements I l/2 lk = r sm N + 1 I :' and v; with elements \ I where lk v;1 = sm N + 1 r { 1 + 2y (upwind), 1+1 (centered). 1"'( We will derive1:results pertaining to the set of eigenvectors as a whole, along I I with oscillatory subsets: those associated with when k M, when j !:M, when both j,k M, where I ,, I I M { Nl (N odd), 12 (N even). We OJir analysis with the upwind discretization on the patch. ReI .I call that this dj.scretization belongs to a problem with Dirichlet conditions I 'I at the fine; Define the quantities c1e =cos and c; =cos I for 1 k,j N. The elements of :E are found to be 1: : I, Let f(y) 1 2 +:: y y'1 + 2y. Then, ., 1 f'(:) = 1y'1 + 2y > 0 Vy E (O,oo), f(O) = 1, hence, f(y) 1 > 0 Vy E [O,oo). 108
PAGE 115
Now, 2 + 1 J1 + > 2 + 1 J1 + 2y 1, so c; I I c; < 1. 2 + 1 vf1 + Restrictiiig oJr attention to 1 E [0, oo ), we are now able to infer the foli' I lowing (u1) C ,(0,2) \1y. i : ( u2) The iteration is convergent \1y <===> w E [1, oo ). I I (u3) C (11(..,),1 + /(0) = 1, and f'('Y) > 0 V'y E (O,oo). I (u4) is convergent for a given"( <===> wE [1+1i1(..,)'oo). (u5) 1 M. ,,2+.., J Notice for all values of y, if we consider the entire spectrum, then i I w = 1 is in that there is no way to choose another w to accelerate the iterati6n. This is unfortunate as 1 + 0 since then the spectral radius, i : 1(..,), apprQaches unity. Therefore, it is neccessary to begin to think of the I iteration a s:r;noother for multigrid. Properties (u3), (u5), and (u6) imply that the optimhl w for the subset corresponding to j or k > M is the I reciprocal .of midpoint of the interval, [1 2!..,, 1 + i.e., I _1_ = 1 + 1/ f('Y)1/(2 + y). Wopt 2 The convetgence factor associated with the set of eigenvalues for which j I ,I or k M, is 1/ f('Y) + 1/(2 + y) Popt = r. Wopt = 2 + 1/ f('Y)1/(2 + y)' 109
PAGE 116
I where r ;is the radius of the interval. Also, (u3) and (u5) imply that the for the subset of E corresponding to j and k M is the ' I reciproca.! of midpoint of the interval, [1, 1 + /h)], i.e., I 1. 1 =1+. Wopt 2/('y) I I Moreover:, .. 1 1 Popt = 1 + 2/(1')" Next k1 e the centered discretization on the patch. Here it I I becomes to separate the problem into two cases. When "'( 1, the of B1 L are real; otherwise, they are generally complex. I Consid.er the case 1. Then we have that I I I I "'(to the interval [0, 1], define f('Y) = 2(1"'(2 ) 112 Then, as I with the case, we can establish the following properties: : : I I (cl) E(y).' C (0,2) \1"'(, r I I ( c2) The 1iteration is convergent \/"'( {=:::::} w E [1, oo ). I. (c3) Eb),' C ('ittr)' 1 + 1t..,))' /(0) = 1, and /'('y) > 0 \/"'( E (0, 1). I ( c4) The is convergent for a given"'( {=:::::} w E [1+1it(r), oo). I (c5) 1 < A1c; when j M. I I Again, for ;all values of"'(, w = 1 is optimal, and if we consider the operator IBJt L rbstricted to the eigenspace corresponding to ( c5), then the i' I analysis is as for the upwind discretization. 110
PAGE 117
I I I I I I Now onsider the case "" > 1. The eigenvalues are I' i I I ,. c 1: = 1 3 1,. [ 3 2 + i( 12 1 )I/2clc where i 41 v'11 Suppose we rewrite these as I ; :, c;(2i{"Y2 i: = 1 4 + (12 1': I Then, wel1:havf I' I i so :E is in the circle of radius 1/2 centered at z = 1. If this is all f : we know :E, then w = 1 is again optimal and the iteration cannot I I ,,11 I be acceleJ::'ated'. On the other hand, we have the attractive result that the II I f,ctor associated with this region is 1/2 and it applies for all values of rr.'[. Jlut more can be said. First of all, it is clear that A1c; with ,' I j M the right half of this disk. Now suppose we rewrite the II 1 j bi e1genva OJice more, as a + : r \i with D1e l: 4 +1 ( 12 1 We can now find a eo satisfying ; I i I 1: I namely, coil',cl I 4. Therefore, :E is contained in the union of two disks with centers =i:3/4 and z?' = 5/4 and both having radius 1/4. The optimal I I W for the rgh,hand disk is Wopt = 4/5 and the associated convergence factor is P+t = \1/5. This provides a convergence factor for the eigenvectors associated Alcj, k M, for all 1 1. II ; I When ljne relaxation is used as a smoother for multigrid, these prop1 I erties to provide guidance in developing an adaptive strategy ' ' :1 111 II,
PAGE 118
., ': J: for its approp,riate weighting according to the relative ellipticity (i.e., "'( as a 1:, I : function ?f mfshwidth) of the discrete problem on each grid level. In order to convergence rates the weighting should be to damping of oscillatory eigencomponents on fine grids and to I perform damping of all eigencomponents on the coarser grids I I I (in it is desirable to solve exactly, and cheaply, the discretized I !. problem bn the coarsest level). : I We criterion in a different way, to identify a robust weighting scheme; words, to identify a fixed w that yields a good general I I multigrid1"smdother for our problem. For this to be possible, such a paI' I rameter provide a good smoother on arbitrarily fine grids. Some experimeJtatibn shows that the optimal w for the purely diffusive problem (i.e., withl: E : 1 and f3 = 0) lies in the interval [4/5, 2/3]. Using a fixed w in this yields Vcycle rates approximately equal to 1/3. On the I' other hana1, properties derived above show that w = 4/5 provides good I global smoothlng rates for both discretizations when 'Y is large. We use ', I line with this latter weighting, exclusively, as the smoother for II multigrid in test problems of the next section. These test problems also II, t involve of higherorder upwinding on the fine patch. Although we have not +al+ed the performance of line relaxation with respect to this discretization, !we find that in practice weighting strategies for standard work well for the higherorder method. In the !ection, we remarked on the effectiveness of an implementat i tion of FA uses relaxationbased strategies to solve global coarse grid ,, I and fine pttch \subproblems. With this algorithm as the starting point, it II I' is natural modified algorithms by limiting the number of relax! I ations or used in solving these subproblems. There are two reasons I 112
PAGE 119
for doing: First, accurately solving the subproblems can be relatively I ,I expensive;:, when there is much refinement on the patch. Sec' ond, these subproblem inversions accurately may be wasteful I I in terms ?f efficiency of the twolevel algorithm; good convergence rates I for the algorithm as a whole can be obtained using less accurate inversions. This is true when the twolevel algorithm is implemented in a ' I nested on a sequence of composite grids (see the next section). Along with this; change in strategy come some important considerations in im.1 II algorithm. First, the iterative methods that are used on the subregions require an initial guess for the solution. For the fine patch problem, this can be obtained from the result of performing I II Step 1 Se.ction 4.1) of the FAC algorithm, the global coarse solve; we use as initial guess in performing the fine patch solve the coarse patch I i1 of solution of the global coarse problem interpolated to the I (This is definition ofFAC given in [33] which we have suppressed I heretofore1 due' to the assumption that the fine patch problem is solved I ;, exactly.) we use zero as the initial guess for the global coarse I problems. For this algorithm, we use XF to represent an approximation to L"F1 I i.e., we wnte th,e result of the fine patch solve as I Whether we can actually write ZF this way is a technical matter that depends on the method of approximate inversion and the fact that this 113
PAGE 120
method an initial approximation to the solution of i I :1 In any case, t;he notation may be considered a formal way of denoting the I ,[ process of the true ZF with its approximation. I 5.3.2 FAC and Nested Iteration I ' In we combine nested iteration with the twolevel FAC scheme td obt:ain a full multigrid version of FAC. Our goal is to come up with an algorithm that yields an approximate solution satisfying I I (5.5). fup multigrid approach is a wellknown methodology [7], (8], [25], (33] for achieving the latter with optimal efficiency. It uses a weak (the use df thi's term will be explained shortly) algebraic solver in a nested I I fashion on a sequence of successively finer grids as follows. On a coarse grid I Jl indexed br 2h, let denote the exact solution to the discrete equations, !, I and let denote an approximation to this solution. Now, let be the ultimate on this grid obtained by applying the algebraic i I solver with as an initial guess. This serves as the initial guess on the next finer jgrid: = and the process is repeated on finer grids. At it is useful to assume the existence of an approximation I ;, subspace in.1the element methods, so that the natural isomorphism between and functions, uh may be utilized. By the notion of a weak solver we mean that the intent of its use is not to solve the I problem ori. hgrid exactly, but rather to solve it to the extent that, for I ., the given I ,I I is .for'' some moderately large constant tt that is independent of I h. To small It should be, the following inequality (taken with 114
PAGE 121
i' :, slight from [33]) can be used. Recall from Section 5.3 the lluit" II cha. Also, notice that, in terms of elements of i" as defined in [30]), we have = u;". In addition, I ,, assume desired inequality is true for 2h: Then, I I llu2"u.;"ll c(2ht. I I = 1111" u.;"ll 1111"u.2"ll + 11112"u.;"ll lluhull + llu2hull + llu2hu.;hll cha + c(2h)0 + c(2h)a = c(l + 2a+t )ha. Thus, a sufficicmt condition for the desired hinequality, (5.5) to hold is I I I I \: I' We note that this condition is sufficient, in practice it may be pes1 ,' simistic. also, that n. depends by way of a on the particular choice I of norm ahd discretization type. Later in this section, we present numer1' It ical resultr indicate the actual number of FAC iterations required to yield solutions with appropriate accuracy as measured in the I' 1'. maximum: and\ L2 norms. : The foVowifg description of nested iteration makes use of the succession I, I, of g?ds leading from the global coarse grid to the composite grid I I! on which al.: is desired, each composite grid being obtained as in I Section 3 .. ,, from the global coarse grid by refinement on the coarse patch. I : Let K 'pe positive integer and, for k = 1, ... K, let (he denote the I ,, composite .grid![ obtained from the global coarse grid by adding k levels of 1: IJ refinement\':to the coarse patch, each level corresponding to refinement in I, meshwidth\. by factor of two. Let (k stand for the global coarse grid ; !: ' I, 115
PAGE 122
"I (if h = lf2mj: is the meshwidth of this grid, then the kth fine patch has II II = hj2k) The algorithm uses the FAC iteration, specifically, a. twolevel I li using flo as the coarse grid and incorporating approximate inversioni.fodthe global coarse and fine patch subproblems, on each of the composit in {Gtc}f=1 The particular algorithms are distinguished by lettinJ,x;i denote the approximate fine patch inversion used on level k. I We use .X to approximate inversion of A on the global coarse grid. I L. I The algorlthrri. begins with the exact solution of a representation of the I i partial differential equation on the global coarse grid. The global coarse I I grid also its usual role in the implementation of the FAC iterations I ,I on each of thei' composite grids. I Let be positive integers and be be the composite operator I t: for the H)l cdmposite grid. Also, let composite grid transfer operators, I' 1[ I ,. ,a_1 : and ,a1 fb_1 be linear mappings. h I I I I II Nested FAC : I' Step 1. finest composite righthand side, rx. successively to all I :: 'd /k1 k K K 1 1 coarser gn s: !),_1 = :!..lc D,, = , I : Step 2. = !a exactly using an appropriate solution technique ;I I ,, for the global coarse equations. I ,\ I ,, I Step 3. Fpr kj= 1, ... K: set z2 = Approximate the solution of : :! rl starting with initial guess z2 and using /Lie FAC iterations with\XF,:; in place of L"F: (approximate inversion on the patch) and I :1 X of A1 (approximate inversion on the global coarse grid); denof,e final approximation I I I I' I I' I 116
PAGE 123
i I ;< I ' I. I In of this section, some examples of the performance of a of the nested iteration algorithm are presented. I 'I For the test problems of Section 3.3.1 are used. Each ,, set of corresponds to a set of tests from that section, so that discrete for nested iteration ( eh, = 1L uh,) may be compared with those by exact solution of the composite equations ( eh = uith) 'I For convbnieJce, the pertinent results from Section 3.3.1 are reproduced I here. of nested iteration uses blockJacobi based it1 erations Vcycles of the {1,1) variety as described in [7], or just relaxatio4s Multigrid Vcycles are used to solve fine patch ,, problems la.pptoximately whereas, depending on the parameter ranges and ,. I type of di 1 scretization, V cycles or simply relaxations are used to solve a pI I I 'I proximately the global coarse problems. For each test, the number of FAC used on the kth composite grid is kept constant, so the I' 'i index has l bee* suppressed. The same is true for both VF, the number of l' [1 Vcycles ttsed ;:to solve approximately the fine patch problem during each 1.1 :i FAC iteration;[ and vc,., the number of relaxations or Vcycles used during I' ;, each corrse solve. In all instances, relaxation used as the basis for multigrid is damped by the parameter a = 4/5 (motivated in 1:: I Section I I I Table S1.1 sB.ows relative errors as functions of patch refinement and the 1:. :: number of(FAQ iterations used on each gridlevel. The test uses Problem 3.1 I : with E = 104 h = 1/32. As in Tables 3.13.2, upwind discretization is I I used on component of the composite grid. The FAC method is of ;_I jl type U4 or:U3I{ depending on whether standard or higherorder upwinding, I t respectivelr, is ;jused on the fine patch. Here, VF = vc = 1, and relaxation I is used for 'the global coarse problem. :I 117
PAGE 124
!.I I : l I m I 1i' 2! 3l 4j 51: m 11 I 2 !: 31 I 41 5j Table 5.1 Relative errors as functions of patch refinement when nested iteration is used. U;pwind/Upwind Upwind/ Hup U!4 Method p. = 1 U3C Method p. = 2 llehll2/llull2 llehlloo/llulloo II ehll2/ llull2 L}102 10 1 1.125 2 5.546 2 3.585 3 102 5.477. 103 1.547 .lQ2 9.605 .lQ4 3.221. 10 "2 2;731 "3 3.991 "3 2.488 "4 1.1p74. 102 1.367. 103 9.669 4 6.686 5 7 .. 758. 10 3 6.847. 10 4 2.290 4 1.701 5 Table 5.2 Relative errors as functions of patch refinement when ezact solution is used. Upwind/Upwind Upwind/ Hup llehlloo/llulloo llehll2/llull2 llehlloo/llulloo llehll2/llull2 1.383. 101 1.069. 102 5.551 2 3.670 3 6.q69. 10 2 5.359 3 1.549 "2 9.532 "4 102 2.687 3 3.977 3 2.500 .lQ4 1.576. 10 2 1.346. 10 3 9.605 4 6.589 "5 3 6.735. 104 2.273 4 1.663 5 The Jsults in Table 5.1 should be compared with those in Table 5.2 I (the reproduces results in Tables 3.13.2 of Section 3.3.1). I I' comparison shows close agreement between the exact solutions and those obtained version of FAC. It is interesting to note that morel'work is required (two FAC iterations, as opposed to one) to I' ., obtain solutions when the higherorder discretization is used. I ., This behayior is in accordance with the requirements posed by the earlier i.. '. analysis. shows results for a similar set oftests applied to Problem t .: 3.2 and the upwind/centered discretization. The coarse meshwidth I 118
PAGE 125
I 'I ,, II :. '! is h = f = 10s, and the FAC method is of type usc. Results 1,I for p.umbers of Vcycles on the fine patch are compared. The I' II other ful). mwtigrid variables are 11c = p. = 1. Again, X corresponds to J, :! relaxation. 1 I I I j I' .[' I I lo 1: .,., 11 2:.'1 '3 4 5 I I :, Table 5.3 iRelative errors as functions of patch refinement ;: when nested iteration is used. I' Upwind/Centered Upwind/Centered li Method IIF = 1 usc Method IIF = 2 lleh.ll2/llull2. lleh,lloo/llulloo lleh.ll2/llull2 102 1.176 101 1.103 101 9.252 2 2 3.122. 10 2 2.365 2 2.013 2 8.24 7 10s 9.650 .lQ3 6.111 s 5.197 s 2.669. 10 "3 2.372. 10 "3 1.695 "3 1.322 "3 4 6.993 4 4.370; 10 4 3.528 4 Table 5.4 /lelative errors as functions of patch refinement when exact solution is used. I I, I Upwind/Centered I I I I h I I' m 1 2 3 4 5 llehlloo/llulloo llehll2/llull2 1.051 1 8.449. 10 2 2.253. 10 2 1.858. 10 2 5.460. 103 4.519. 10s 1.356. 10 "3 1.122. 10 3 3.383 4 2.800 4 I 1: The in Table 5.3 should be compared with those in Table 5.4 i I' (the latterl'a.re upwind/centered results from Tables 3.33.4). Here, one I I sees a degradation of the convergence of uh in the left column II I of Table the right column, 0( hf;.) convergence is restored by simply I I' increasing rtumber ofVcycles used on the fine patch. Table 5.5 presents I li I I 119 i
PAGE 126
I I I i I results CQrres:ponding to those in Table 3.5 (which are reproduced below in I Table 'ljwo FAC iterations (p = 2) using a single fine patch Vcycle I 'I (vF = required to obtain full multigrid solutions with the desired accuracy) In unlike the previous problems, Vcycles were used for I I 'I the globaJ solves. Recall that for this test, u has a rapidly changing I I ,1 componefit is not confined to the patch. This helps to explain the I ! additional cotnputational effort required on the global coarse grid. Also, I fj ,, I note work (two Vcycles, rather than one) is required on this I grid wheD: the change to a more accurate coarse component discretization [I is made. ,I ,, 'i Table 5.5 errors as functions of coarse component refinement when nested iteration is used. :I Upwind/Centered Centered/ Centered ., II U3C Method vc= 1 C4 Method vc = 2 me (riJ. = 2) lletlloo/llulloo lletlloo/llulloo 3 2.026 1 1.836 1 4 7.970 10 2 7.960. 10 2 q 3.625. 10 2 2.154. 10 :& 1.860. 102 7.081 loa : Table 5.6 errors as functions of coarse component refinement when ezact solution is used. I Upwind/Centered Centered/Centered I i .! I e 2) llehlloo/llulloo II ehlloo/llulloo :!3 I 1.963. 101 1.935. 101 :i4 7.981. 10 2 7.468. 10 :& :5 3.921 2 1.903 102 I 2.152. 10 2 5.628. 10 3 120
PAGE 127
k We that, in all of the preceding tests, the unmodified stencils of I' i Section 3 .. 3.31were used. Recall that in Section 4.3, degradation in conver11 gence oftfhe rwolevel algorithm was attributed to use of the unmodified stencil of. global coarse grid. However, the convergence results of this section satisfactory. The reason for this is that when the twolevel I, I I : algorithm pefrormed poorly, it did so at the tangential boundaries of the patch. :these interface values are being solved for very accurately in the initia1 of the algorithm, when the problem on the global coarse :., I grid is sol.;ved.: A very accurate initial guess is provided to the fine levels by I this coarse grid and, therefore, the dependence on the initial guess noted in the no longer exists. I'' I In test case of this section, we present a problem for which the i I use of near one of the corners of the domain {l is appropriate. i:', i: This problem! has been chosen so that the solution simulates a boundary I layer perturbed problems. Indeed, the solution has I been ifrom the wellknown onedimensional problem r 1:' eu"(z)+u'(z)=O, zE[0,1], I, I u(O) = 1, (5.7) i u(1) = 0, whose sol!tioJ [18] has an 0( e)width boundary layer near the left endpoint ],11 ,I of the Our twodimensional solution corresponds to that of (5.7) with the t1:,ouddary layer translated to, and rotated about the point (1,1) : I at the corner of n. I I r I Problem ' ,, I I I Let = J(z 1)2 + (y1)2 and r > 0. The solution is I 121
PAGE 128
' 'i Appropriate solution values are used to define Dirichlet conditions as before : II on the west, and north segments of 80. As for the east boundary, i ,; it is into northeast and southeast segments where appropriate Dirichlet[' Neumann conditions, respectively, are specified (see Figure ' 5.1). the Neumann condition is given by = 0 on 80sE I '' I I anN [, i I i' \,L' 80s I, 11Fig. 5.1. Partitioning of 80 for Problem 5.1. ,, ,; I' 'I For thl,:s prfblem we use a composite grid as shown in Figure 5.2. Note that the has only two sides. Along anNE, for example, the in'll I terface is by a fine boundary whose nodes are members of the II ' composite!:gri4. For convenience in testing, we place the division between I I, I 'I northeast segments of 80E at the southeast corner of the I interface. we use a coarse component meshwidth of h = 1/32, and II the patch with sides having length 1/8. 122
PAGE 129
I i I. I .\: i i 1: I 'i I II The results of tests for this problem withE= .001 and r = 10 are shown I' ;I in Table\5.7,!\where solutions by nested iteration are compared with exact ones. iteration variables that yield these results are p. = 2, i I vc = 1, = 2, with one exception. For the upwind/centered case, II with m ;, 5, order to keep the order of convergence from degenerating, 1 'I we find that is necessary to use four fine patch Vcycles per FAC iteration I on the bomposite grid. I' :i I ., I :1 I 0 0 0 0 0 0 0 0 0 0 0 0 0 0 I 0 :' II 0 0 0 0 0 0 0 0 0 !i 0 0 0 0 0 0 0 0 0 :1 0 0 0 0 0 0 0 0 li 0 ,, ,I 0 0 0 0 0 0 0 0 I 0 0 0 0 0 0 0 0 0 I I 0 0 0 0 0 0 0 0 0 'I 'I I :I Fi'g. 5.2. Composite grid with coarse ( o ), fine ( ), :j interface ( 0 ), and slave ( o) nodes. I I I I I I I I I I I 1: 123 I
PAGE 130
m 2 3 4 5 Table 5.7 Relative errors as functions of patch refinement for nested iteration (left} and exact solves (right} in the discrete version of Problem 5.1. Upwind/Upwind Upwind/Centered ., I U4 \Method Exact Solve U3C Method Exact Solve llehiiQ)/IIuiiQ) llehiiQ)/IIuiiQ) llehiiQ)/IIuiiQ) llehllao/lluiiQ) 2 4.900 10 z 3.072 z 3.070 10 z 3.507 2 3.491. 102 1.746 2 1.782. 102 1.941. 10 2 1.954 2 4.158 3 5.551. 10 3 2 1.035 2 1.246. 10 3 1.293. 10 3 ,' ,, 124.
PAGE 131
' i 6. Summary of Results ;: In this thbsis, we have studied multilevel methods for the solution of I : elliptic equations in one and two dimensions. We conI centratetl, on problems for which local grid refinement is desirable, in particVJar,;:problems for which the use of a very coarse grid is suitable, except subregion of the domain, where significant refinement is needed. i ., Althotigh ,the results obtained here for the onedimensional case are I I generally! of practical interest in themselves, they provide important insight behavior of the methods we study in higher dimensions. I This is wrere :'our true interests lie. In this chapter, we summarize some of the more results obtained in this study. We emphasize those re sults a bearing on practical use of the solution methods presented I' : here, cases give recommendations for their application in the I I implemen;tatiqn of these methods. In twb the results of this study were obtained for the d t I diff .: t a vee tonr Sion equa ton I I eLlu + 'Uz + cu = f. (6.1) This equjtion.: has a particularly attractive form in that its flow velocity I I .' .lies entirely in ;the zcoordinate direction. In Chapter 2, we noted that this I I equation he obtained from a rather general advectiondiffusion equa, I I tion in twq1 via an orthogonal coordinate transformation induced I I by Algorithms for performing this change of coordinates I are descri\?.ed in [9] and [13]. 125
PAGE 132
! i: I I I i I' The of the above form for the results obtained herein are 'I twofold,l::and correspond to two hierarchies of sQlution procedure. At the I" I globallefel, fe considered the FAC multigrid method for the solution of (6.1) on a composite grid. This method belongs to a class of w1Uch place the computational effort in solving, or ultimately, I partially standard problems on subregions of the composite grid. I Our results for this method relied, to a large extent, on the strong oJ'edikensional character o(6.1) when E is small, in the sense that J: i when this is true a rigorous analysis in one dimension gives much useful I I insight iiito t:Q.e true behavior of FAC convergence in two dimensions. I' In 4.2.1 we analyzed the eigenvalues of what is essentially the II l governing convergence of twolevel FAC when exact so' I lution of problems is as.sumed. For a given discretization on the composite grid we also considered the discretization type used at the coarse level of tAe tJolevel scheme. Four discretization& of the onedimensional I': j problem The first two cases (the U4 and 04 methods) in I volved of either upwind or centered discretization throughout both grids. In parameter ranges, we concluded that the spectral radius l I of the iteratio:h matrix, in these cases, was zero. This allowed us to further r I conclude I: in this context the twolevel FAC iteration may be viewed I' I as a solyer, producing a zero residual after three iterations. Furthermore, was relatively insensitive to the choice of the composite '' I grid used fn tte coarse level. In particular, the global coarse grid could be used. i' I The we considered, (UC)2 used upwind and centered differ ence respectively, on the coarse component and the patch II' I, component of !the two grids. The l!pectral radius of the iteration matrix : i i1 I 1: I t: !, 126 1,, i ,,
PAGE 133
in this crse was bounded by one half, but the validity of the result del pended a :familiar requirement for the centered discretization, namely I that h be sufficiently small with respect to the diffusion I' T,his was in conflict with the requirement for the upwind part I of the and motivated the use of the fourth strategy, U3C, !' : which employed upWi.nding throughout the coarse grid and used the upwind/ mix on the composite grid. For this case, the spectral radius I 1: zero as the refinement level on the fine patch increased. I As fo# con:vervence in two dimensions, we saw in Section 4.3 that computationJ.I there conformed, for the most part, to the analysis in one An exception we discovered was a significant degradation in rates for upwinding used on the fine patch. This problem I lo was more: pronounced when the diffusion coefficent was small. We found : that this I:Was: caused by slow convergence of the residuals at the tangenh ;I tial bounaaries of the patch (i.e., the ones parallel to the flow direction). I . By analyzing the matrix stencil associated with an idealized version of the I I problem case that the interface component is all tangential boundary), I. we identi4ed a; modification of the stencil in these regions that restores the convergerlce The modification had essentially the same effect when l, h' higherord,er U:pwinding was used on the fine patch. In fact, results ob I tained for!,standard upwinding in this study generally apply as well to the case of hiJher,order upwinding, These include global convergence behavior i ,. of twolevfl FtC as well as weighting of line relaxation in devising multigrid smoof.hers for subgrid problems (see Section 5.3.1). An instance where I results fori':the two differ occurred in the nested iteration all' gorithm of Chapter 5. There it was seen that more work was required by . I I the higher1torder method on each of the intermediate gridlevels in order to 1: 127
PAGE 134
I I i i I 1. obtain alsolution with the desired accuracy (see Table 5.1). This fact is in I I accordarl,ce the theory for full multigrid as discussed in Section 5.3.2increasin\g order of accuracy of the discretization increases the estimate I ,: of work by the algebraic solver. ThesJ of this section indicate that the most attractive methods for tb moderately advectiondominated problems are the us H and usc In two dimensions, we noted that the former method may be in parallel with a modest increase in expense over the u4 method, at the same time retaining its convergence behavior. For I more mildly problems, the U3C and C4 methods may 1 .: be more appr9priate. 1: We nbte that if the twolevel schemes U4 and us H of Chapter 4 are I ' I used as solvers (i.e., viewed as a preconditioned iterations and I I described in Chapter 5 by a polynomial method), and E I is small, lt is .;essential that the modification to the coarse grid stencil be I performed. Just how small E must be for this to be necessary depends on the I II I meshwidth ofjboth the coarse and fine components ofthe composite grid. I : In our with the coarse meshwidth fixed this point occured I later (with respect to decreasing e) as the number of levels of fine patch I :: increased (see Tables 4.64. 7). I' 'fe do not advocate such use. The reason for this is that the twolevel as presented in Chapter 4, requires exact solution I of the and fine patch) subgrid problems. Furthermore, as I I noted in [23), with respect to diffusion problems, the number of itera', tions ,to solve the composite equations grows rapidly when a few l1 1 steps of a multigrid method are used to approximately solve the fine patch II subproblems. 1Therefore, the computational complexity of the twolevel 128
PAGE 135
i I I; t I t I 'J! i scheme far from optimal. Yet, there is a way of combining the efficient approximation of subgrid solutions which does have grid) complexity of full multigrid. The approach we took t! i in is the classical multigrid approach of using the scheme in a I I nested Of a sequence of successively finer grids. In particular, we used I' ,, the i:scheme as a weak algebraic solver on each of the composite grids lyiJg the global coarse grid and the true composite grid, and I used subproblem approximations obtained via either relaxation II I, or multiglid Vcycles. Based on the computed results of Section 5.3.2, we I'' I I, I. recommeJ;td following strategies for the implementation of this method. I I II II Suppose that :!standard upwinding is used on the coarse component. Then, I I one or twb.:. unweighted line relaxations should suffice in "solving" the global II I coarse This depends on the relative coarseness of the mesh, .,. I but the 1/ /(1) (the asymptotic convergence factor) derived in ,,., I Section 5\:.1\and (the required convergence factor) derived in Section 5.3.2 for guidance. We note that when the Reynold's number 7 is 4.2 two steps of the line Jacobi iteration yield an asymp,,, I totic factor bounded by 1/9. Convergence factors for the use of ott the fine patch are more difficult to predict. However, in Section !e noted that for a worse case situation (i.e., for the diffusion !: !, problem, i'vhiCJt resembles the case of infinite mesh refinement) the use of 1 1 the a= 4/5 yields Vcycle factors close to 1/3. So, again, yield a factor of 1/9. Using these strategies, we found that, only one or two FAC iterations per gridlevel were required II to approximate solutions with the accuracy of the exact solution of the Notice, with a single FAC iteration on each level, this thJt the total work on the fine grid is less than that of four VI: I It I' li r, I' i: 129
PAGE 136
I: I cycles performed on the finest patch. Hence, this method has a fine patch 1, I complexity comparable to full multigrid. ;I I We xipte that the above results are specific to the finite volume dis1 cretization we used for the problem (6.1). We were motivated to use this !, method by the fact that line relaxation then provides a highly effective I. I solution for uniform grid problems while decoupling the nodes in I I' such a way preserves the parallel nature (6.1) possesses when c: = 0. I ,, We c\9se discussion with two observations that might be used to yield improvements in our implementation of the nested iteration l1 First,kve a fixed relaxation parameter a for the damping of relax, : ,, ation used as a multigrid smoother on the fine patch. The analysis I : of SectioK 5.3.1 may be used, however, to fine tune the weighting of rel laxation to gridlevels in the V cycle. Specifically, we provided I convergeJce factors associated with global and oscillatory eigenspaces for 1, general gpdleyels. These may be used to determine when to unweight re, laxation the grids coarsen (when the upwind discretizations are used), ,, and whatl.weighting to use to most effectively damp oscillatory error com1: ponents Qh fine grids. I. : I' Second, when the centered discretization was used on the patch, we employed\.the fPC method on each of the composite gridlevels. However, t:' our with the twolevel version of this method showed that it I can be slow if::the fine level is too coarse. This observation motivates the use of up.lnding on the coarser composite gridlevels and applying the U4 or U3 H Jethod, and then switching to centered discretization and the 1 : U3C on the finer composite grids. It may be possible to apply the results of &ur analysis of this method to determine when it is appropriate to 130
PAGE 137
, switch di,screVzations. We saw (in Table 4.12) that the quantity 1..\2 (..\2 = ' );,:, ass?ciated with the latter method in one dimension, gave good estimates of 'actual convergence rates for the method in two dimensions. I i : Althougll this was true with the use of exact inversion of subgrid problems, I the were generally somewhat pessimistic, and so may also apply I' I when a reaso;nably accurate approximation method, such as classical full ,I multigrid. is to solve these problems. Since the centered method has !' I' second accuracy, an appropriate strategy would be to change to the I I centered method on the level m fine patch when 1 ..\2 is less than 1/9 or I : 1/3, as one or two FAC iterations are used. I 131
PAGE 138
I I I i REFERENCES [1] Approzimation of Elliptic Problems, Kreiger, Huntington, ]9soJ j' I [ [2] R.Ei Bapk and T.F. Chan, PLTMGC: a multigrid continuation profor parameterized nonlinear elliptic systems, SIAM J. Sci. S,tat. !!Com put., 7 (1986), pp. 540559. II I [3] R.El.. Ba.Ilk and M. Benbourenane, A Fourier analysis of the twolevel basis multigrid me'thod for convectiondiffusion equa tionsJ. in Proceedings of the Fourth International Conference on Decomposition Methods for Partial Differential Equations, :R!,. Gl6winski, Y.A. Kuznetsov, Q. Meurant, J. Periaux and O.B. [4] solution of problems on regions partitioned into substructures, SIAM J. Anal., 23 (1986), pp. 10971120. I I [5] J.E. Pasciak, and A.H. Schatz, An iterative method for elliptic problems on regions partitioned into substructures, Math. 46 {1986), pp. 361370. '. I [6] J .H. R.E. Ewing, J .E. Pasciak, and A.H. Schatz, A precontechnique for the efficient solution of problems with lo grfd Computational Methods in Applied Mechanics aj:d El;ngineering, 67 {1988), pp. 149159. [7] A. Multigrid Techniques: 1981, Guide with Applications in Fluid :pynamics, Department of Applied Mathematics, Weizmann l,stit4te of Science, Rehovot, Israel, 1984. [8] W.L.i1 A Multigrid Tutorial, SIAM, Philadelphia, 1987. [9] [10] I I D.L. ,1Brorn, R.C.Y. Chin, G.W. Hedstrom, and T.A. Manteuffel, aymptotics and domain decomposition, Technical Report UCRLJC106336, Lawrence Livermore National Labora1'' I tory, JJ991. 11, I Z. Cai and S. McCormick, On the accuracy of the finite volume element jl' I methoi/for diffusion equations on composite grids, SIAM J. N umer. Afial.,i27 (1990), pp. 636655. I 132
PAGE 139
[11] Z. yai, :J. Mandel, and S. McCormick, The finite volume element method for diffu.sion equations on general triangulations, SIAM J. I. Anal., 28 (1991), pp. 392402. .. [12] Chan and D.E. Keyes, Interface preconditionings for domain4,eco,:,posed convectiondiffusion operators, in Proceedings of the Third International Symposium on Domain Decomposition MethL Partial Differential Equations, T .F Chan, R. Glowinski, J. IMriahx and O.B. Widlund, eds., SIAM, Philadelphia, 1990, pp. 2452()2. I [13) Chin, G.W. Hedstrom, F.A. Howes, and J.R. McGraw, PardUel of multiplescale problems, in New Computing hi Environments: Parallel, Vector and Systolic, A. Wouk, ed., SIAM, 1... I 1986. 1.' I : [14] R.C;Y. Ohin and T.A. Manteuffel, An analysis of block successive overfor a class of matrices with comple:c spectra, SIAM J. Anal., 25 (1988), pp. 564585. 1'. [15) R.Cj:Y. Chin and G.W. Hedstrom, Domain decomposition: an instru.J.ent :of asymptotic numerical methods, Technical Report UCRLJp10,4693, Lawrence Livermore NationalLaboratory, 1990. I'' ' [16) P.G.[;.Ciar,let, The Finite Element Method for Elliptic Problems, North Holland, Amsterdam, 1978. I I [17) R.W:Ii. Manifestations of'the Schur complement, Linear Algebra A!ppl.) 8 (1974), pp. 189211. J I' [18) W. Matched Asymptotic E:cpansions and Singular Perturbattns, NorthHolland, Amsterdam, 1973. [19] H.C.!:.Elman, Iterative methods for large, sparse, nonsymmetric sys tems (()/ linear equations, Technical Report 229, Department of Science, Yale University, 1982. i: H.G.iElm;an and G.H. Golub, Iterative methods for cyclically reduced linear systems, Technical Report UMIACSTR881 81:, Dept. of Computer Science, University of Maryland, 1988. [20] I,' i:. 133
PAGE 140
[21) [22) [23) I I ,, Domain decomposition techniques for efficient adaptive local :grid in Proceedings of the Second International Symposium on Domain Decomposition Methods, T .F. Chan, R. <;?lowfnski, J. Periaux and O.B. Widlund, eds., SIAM, Philadelphia, 1989, pp. 192206. I : R.E': Ewing, R.D. Lazarov, and P.S. Vassilevski, Local refinement techfor elliptic problems on cellcentered grids, I: error esti Math. Comp. 56 (1991), pp. 437461. I I R.El R.D. Lazarov, and P.S. Vassilevski, Local refinement techfor elliptic problems on cellcentered grids, II: optimal order iterative methods, manuscript. : [24) W. ;ijackpusch, Multigrid convergence for a singular perturbation prob Linear Algebra Appl., 58 (1984), pp. 125145. I [25) W. ,Hackbusch, MultiGrid Methods and Applications, Springer, I [26) [27) [28) [29) Verlag, New York, 1985. C. Numerical Solution of Partial Differential Equations by Element Method, Cambridge University Press, Cam bridge, 1987 . i. !I D.E .. ,KeY,es and W.D. Gropp A comparison of domain decomposition techniques for elliptic partial differential equations and their par allel SIAM J. Sci. Stat. Comput., 8 (1987), pp. I I s166s202. I; : I 1 1 S.J. Linear Algebra with Applications, Macmillan, New York, 1980. ,I I I, J.R. Rice, and D.H. Thomas, Direct solution of partial differential equations by tensor product methods, Numer. Math., 6 I n.964J, pp. 185199. [30) J.F. and F. Musy, Multigrid methods: convergence theory in a framework, SIAM J. Numer. Anal., 21 (1984), pp. 657671. i:, :i [31) T .A.! Manteuffel, The Tchebychev iteration for nonsymmetric linear Numer. Math., 28 (1977), pp. 307327. 134
PAGE 141
[32] [33] S.F{.Mc9orm.ick, Fast adaptive composite grid (FAG) methods: theory for the variational case, in Defect Correction Methods: Theory Applications, K. Bohmer and H.J. Stetter, eds., Computations Supplementation, 5, SpringerVerlag, Wien, 1984, pp. 115122. S.F.i McQormick, Multilevel Adaptive Methods for Partial Differential Equa#ions, Frontiers in Applied Math, Vol.6, SIAM, Philadelphia, 1989.i I [34] W.F; .. Unified Multilevel Adaptive Finite Element Methods for Elliptic Problems, Doctoral Thesis, Department of Computer University of illinois at UrbanaChampaign, 1988. I [35] K. Ressel, Gebietszerlegung und Mehrgitterverfahren: Varianten1 Untersuchungen und Aspekte der Parallelisierung, Piploma Thesis, Koln, 1990. I (36] P.J. Roach, Computational Fluid Dynamics, Hermosa, Albuquerque, 1972.' [37] [38] [39] G. :s;.odrigue and E. Reiter, A domain decomposition method for boundproblems, in Proceedings of the Second International Symposium on Domain Decomposition Methods, T .F. Chan, R. qlowi;nski, J. Periaux and O.B. Widlund, eds., SIAM, Philadel phia, :l989, pp. 226234. I 'I U. Riide;: F'ully adaptive multigrid methods, manuscript, presented at the 5th Copper Mountain Conference on Multigrid Methods, Cop pr M,ountain, CO, April, 1991. I I I i,Saylbr and D.C. Smolarski, Implementation of an adaptive algoiithm'Jor Richardson's Method, Linear Algebra Appl., 154 {1991), pp. 615646. ] ' (40] A. Segal,:: Aspects of numerical methods for elliptic singular perturba: II ti;on problems, SIAM J. Sci. Stat. Comput., 3 {1982) pp. 327349. I (41] G. Strang and G. Fix, An Analysis of the Finite Element Method, PtenticeHall, Englewood Cliffs, 1973. I :1 [42] E. siili, Oonvergence of finite volume schemes for Poisson's equation on nonuniform meshes, SIAM J. Numer. Anal., 28 {1991), pp. 14191430. 135
PAGE 142
I i [43) R.Sl A comparison of successive overrelazation and semimethods wing Chebyshev polynomials, J. Soc. lndust. Math., 5 (1957), pp. 3946. [44) [45) I ' o.:a. Widlund, Optimal iterative refinement methods, Technical Report ,:391, Department of Computer Science, Courant lnstitue of Mathematical Sciences, 1988. I ,, D.M. Young, Iterative Solution of Large Linear Systems, Academic J!>ress, New York, 1971. I I 'i I I I" I I I i ,, I' li' 136
