Fast Multigrid Solution of the Advection Problem with Closed Characteristics.
The numerical solution of the advection-diffusion problem in the inviscid limit with closed characteristics is studied as a prelude to an efficient high Reynolds-number flow solver. It is demonstrated by a heuristic analysis and numerical calculations that using upstream discretization with downstream relaxation ordering in a multigrid cycle with appropriate residual weighting leads to an efficient solution process. Upstream finite-difference approximations to the advection operator are derived whose truncation terms approximate "physical" (Laplacian) viscosity, thus avoiding spurious solutions to the homogeneous problem when the artificial diffusivity dominates the physical viscosity [A. Brandt and I. Yavneh, J. Comput. Phys., 93 (1991), pp. 128--143].
Introduction
Efficient multigrid algorithms for the numerical solution of partial differential
problems normally require good ellipticity measures on all scales of
the problem, which implies that nonsmooth solution components can be re-solved
by local processing [2]. But problems with small ellipticity measures
are marked either by indefiniteness or by anisotropies. In the latter case,
Technion-Israel Institute of Technology, Haifa, Israel
y University of Twente, Enschede, The Netherlands
z Weizmann Institute of Science, Rehovot, Israel
there exist so-called characteristic directions of strong dependence. Some
nonsmooth components of the solution are then advected along these char-
acteristics, and hence they cannot be resolved locally [1]. A typical example
is steady flow at high Reynolds numbers (small viscosity).
When applied to such problems of small ellipticity the usual multigrid
algorithms often exhibit a severe degradation of performance compared to
that seen in elliptic problems. Indeed, most multigrid codes in use today
for solving steady flows at high Reynolds numbers, although yielding a great
improvement over previous single-grid solvers, fall far short of attaining the
so-called textbook multigrid efficiency for general (even smooth) flows. To regain
this efficiency the multigrid algorithm requires modifications that take
into account the anisotropic properties of the operator. For example, it
was shown in [4] and [9] that using upstream discretization and downstream
relaxation-ordering yields a fully efficient multigrid solver for flows whose
characteristics (streamlines) start at some part of the boundary and end at
another without recirculating (entering flows). To obtain efficient multigrid
solvers for flows with closed characteristics, however, different modifications
were proposed, such as defect-correction cycles and residual overweighting
[5]. The main drawbacks of the latter approaches are: (a) they are not likely
to generalize efficiently to orders of accuracy higher than one; (b) they require
W cycles, which may be substantially more expensive than simple V
cycles in parallel computation; (c) they suggest different treatment for different
types of flow, viz., recirculating versus entering flows. The upshot of
the present work is to obtain a unified approach for both types of flow by
employing upstream discretization and downstream relaxation ordering for
recirculating flows as well. In Section 2 we formulate the simple model problem
of advection-diffusion and present the First Differential Approximation
to its discretized form. In Section 3 we present the two-level cycle and use the
approximation of Section 2 in a heuristic analysis for a priori prediction of
the performance of this algorithm. In Section 4 new first-order upstream discretizations
for the advection operator are presented, whose first truncation
terms approximate isotropic diffusivity. These schemes are shown to eliminate
spurious solutions to the homogeneous (i.e., unforced) small-viscosity
advection-diffusion equation, such as those reported in [3]. Section 5 presents
numerical calculations testing the accuracy of the discretization and the efficiency
of the multigrid algorithm and how it compares to the predictions of
Section 3. Section 6 summarizes the main conclusions and further research
plans.
2 The Scalar Advection-Diffusion Equation
We study the scalar advection-diffusion equation with closed characteristics
as a prelude to the study of flow problems. This equation serves well as a
preliminary problem, since the advection part (i.e., momentum equations)
is responsible for the degraded performance observed in the solution of the
incompressible-flow equations by the usual multigrid algorithm [5]. Also,
as is shown for entering flows in [4], the solution-process of the advection
part of the system can be effectively decoupled from that of the elliptic
part that is due to the continuity equation. Hence, efficient solution of the
advection problem is a necessary stage in the development of a fully efficient
flow-equations solver, and [4] suggests that the resulting advection-problem
solver can indeed be used in designing the sought-after flow-solver.
The advection-diffusion equation in two dimensions is
where ffl is a positive constant and a, b, f , and g are given functions of x
and y. Equation (1) is discretized on a uniform grid of meshsize h, whose
gridlines lie parallel to the x and y coordinates. The characteristic direction
of the advection operator in (1) is given (locally) by
where OE is the (local) angle of nonalignment between the x coordinate and
the characteristic direction. We will focus our attention on the particular
case where the characteristics defined by a and b form closed loops (as in
vortices), one of which may coincide with
@\Omega (as in internal flows).
Suppose that (1) is discretized by some stable finite-difference discretization
of first-order accuracy. The main aspects of the problem can be analyzed
by substituting for the discrete operator its First Differential Approximation
[8] and also [1, 2, 9]. For the advection-diffusion equation with
positive but vanishingly small ffl we need only consider the advection oper-
ator, since the tiny diffusion will be dominated by the artificial diffusivity
represented by the truncation terms (except at stagnation points). Let L h
denote a first-order accurate discrete approximation to the advection opera-
tor. Then, by a Taylor series expansion, we generally have
where u h denotes the discretized function, h being the meshsize of a uniform
grid. Here, ~
are functions of x and y, the specific details of which are
determined by a and b and the discretization. The FDA is the approximation
of L h by the differential operator that remains in (2) after the O(h 2 ) terms
are neglected. (Hence, it applies only to sufficiently smooth u h , since the
neglected terms are higher derivatives).
We now assume for simplicity of the discussion that the equation is normalized
such that a 2 introduce a (conformal) local coordinate
system, (-; j), where - denotes the local "streamwise" coordinate parallel to
the characteristic direction, while j denotes the "cross-stream" coordinate
that is perpendicular to the characteristic. Thus,
and
The FDA of L h in the local coordinate system is therefore
~
where, by (2) and (4),
We assume a consistent and stable discretization, which requires that the
artificial viscosity operator represented by the first truncated term be elliptic,
implying that the O(h) part of the operator in (5) is of positive type. Under
special circumstances, such as consistent alignment of the characteristics with
the grid, T h
vanish, and this property is marginally violated. In this
case the physical diffusion term becomes important, no matter how small ffl
may be, and the analysis below does not apply. Accordingly, we will assume
below that T h
1 is large compared to h, which is the usual case.
3 Two-Level Error-Reduction Analysis
We now analyze the error reduction attainable with a two-level cycle using
upstream discretization and downstream relaxation-ordering.
3.1 Two-level Cycle
The proposed two-level cycle for a given discrete problem L h u
defined as follows.
ffl Starting with some approximation to u h perform - 1 (small integer)
relaxation sweeps.
ffl Calculate the residuals, r
u h is the current approximation
to the solution, and transfer them to a twice-coarser grid
2h, multiplied by a globally-uniform weight W .
ffl Solve the coarse-grid problem, L 2h v for the correction.
Here, r 2h is the restriction of r h to the coarse grid.
ffl Interpolate and add the correction v 2h to the fine-grid approximation.
ffl Perform - 2 (small integer) fine-grid relaxation sweeps.
In studying the asymptotic performance of the two-level cycle, the number
of pre-relaxation sweeps - 1 need not be distinguished from the number of
post-relaxation sweeps - 2 . (Recall that we associate asymptotic performance
with the spectral radius ae of the iteration matrix, and that ae(AB) = ae(BA)
for any pair of square matrices A, B of the same dimension). We denote the
total number of sweeps by
In analyzing the two-level cycle we shall make many simplifying assump-
tions. The degree to which these assumptions are justified needs to be judged
by the degree to which numerical results match the predictions of the analyses

3.2 The Model Problem and Analysis
We analyze the two-level convergence for the discrete approximation to (1) in
the limit of vanishing ffl for problems with closed characteristics by considering
the following model problem on grid h.
where L h is, as above, the discretization of the advection operator. For the
domain of solution and boundary conditions we require periodicity in - in
order to simulate closed characteristics, and we choose for simplicity of the
discussion - 2 [0; 1]. The boundary conditions in the cross-stream direction
are not germane in the present context. For simplicity, we let
The main point of our approach is to use discretization that is purely
upstream, and to relax the equation in downstream ordering, starting at
ordering means that we relax a variable only after relaxing all
other unknowns which participate in the equation which corresponds to this
variable (except, perforce, at Thus, a full relaxation sweep results in
the elimination of all the residuals except at a narrow band (of O(h) width)
that stretches from (which coincides with
due to the periodicity). Neglecting the width of this band, we find that the
residual function, r h , which remains after at least one full relaxation sweep
has been carried out, can be modeled by
r h (-;
The fine-grid error v h satisfies the residual equation. Wherever the residual
vanishes we now revert to FDA, obtaining
~
with ~
defined in (5). We now add a further simplifying assumption that
in (5) are independent of j. Hence, we may expand (8) in a Fourier series
in j. For an error component -
(-) exp(i!j) of frequency !, Equations (5)
and (8) then yield
d-
where all the O(h) terms but the first can be neglected in the homogeneous
equation, since they multiply derivatives and are therefore small compared
to d-v h
d-
. The solution to (9) in the interval (0; 1) is therefore given by
Z -T h
ds
where A h
! is the amplitude of - v h
! at which we shall determine shortly.
(Superscripts denote an infinitesimal positive (negative) incre-
ment). In particular, at
1 is the average value of T h
1 over the entire domain, under the assumption
that T his j \Gammaindependent (and recalling that the domain length
in the - direction is one).
It is important to note that D(h; !) is approximately the factor by which
a single relaxation sweep amplifies (reduces) an error component that oscillates
at frequency ! in the j direction. This is due to the fact that, given
upstream differencing, the downstream relaxation ordering yields numerical
integration; and D(h; !) is the factor by which this integration over the domain
reduces the error. Equation (11) implies that relaxation reduces error
components with large ! very efficiently, but components that are smooth
along j need to be corrected on the coarse grid.
the FDA is no longer useful. Instead, we have a jump
in - v h
! that is proportional to the Fourier coefficient of R h (j) corresponding to
frequency !. We denote this jump by ffi h
. The periodicity in - now implies
by (10) and (11)
A h
Now, following the two-level algorithm, we attempt to approximate the
(weighted) residual equation on the coarse grid 2h. We assume that the
same discretization stencil is used on the coarse grid as on the fine. (Note
this important assumption on which the entire method hinges). We also assume
that the restriction operator is such that the jump condition at
is approximated correctly on the coarse grid. In practice this holds provided
that a proper averaging is used, such as full-weighted residual transfers. Analogously
to (10),(11), and (12), respectively, we obtain
Z -T 2h(s) ds
and
A 2h
where W is a constant weight to be chosen. Since the stencils of L h and L 2h
are the same, we also assume -
Also, since we assume that the restriction operator transfers the jump condition
correctly, we have ffi 2h
. Equations (12), (15), and (16) now yield
A 2h
A h
Neglecting again the effects of intergrid transfers and aliasing, we assume
that the remaining error after adding the coarse-grid correction is
Hence, the fine-grid error is amplified by the factor (-v h
! . In addition,
relaxation sweeps performed on the fine grid amplify the error by D(h; !) - ,
as noted above. The two-level error-amplification factor - tl is then given by
the absolute value of the product of these terms. Note, however, that the
determining values of -
! in the coarse-grid correction term depend on
where one begins relaxing on the fine grid immediately following the coarse-grid
correction (since all other fine-grid values at the end of one or more
sweeps are determined solely by the values where relaxation begins, due to
the upstream differencing and downstream relaxation ordering). If the fine-
grid relaxation begins at that is, at or shortly after the point where the
residual was nonzero, then we have
A 2h
A h
However, if the fine-grid relaxation begins at
there is some small overlap in the region being relaxed), then
A h
3.3 Optimal Residual Weighting
By (11) we have and in order to obtain an h-independent
analysis we assume now that D can take on any value in the interval (0; 1).
The optimal value of W is that which minimizes the supremum of - tl over
1). For any fixed W , the supremum is evidently obtained either for
@D
From (18) and (19) we can thus obtain Dm as a function of W and -, from
which we can then calculate W opt (the value of W which yields the fastest
convergence) for either of these cases. For
for (18), yielding -
For larger - the tedious calculations need to be carried out numerically. But
W opt tends to 2 for both cases rather quickly. This is expected, since this
value is the ratio of the Green's functions on the coarse and fine grids for
components that are very smooth in the cross-characteristic direction; other
components are reduced by relaxation. With we obtain for both (18)
and (19)
Equating the derivative of (20) with respect to D to zero, we get as the only
relevant root,
The two-level error-amplification factor is now obtained by substituting (21)
into (20). For sufficiently small - \Gamma2 we may neglect this term, obtaining
In fact, (22) gives an excellent approximation of the maximal - tl in (20)
for any -, erring by less than 2% for - 2. By curious coincidence, the
same asymptotic two-level error-amplification factor for large -, (2-e) \Gamma1 , is
obtained for the Poisson equation on a rectangle using Gauss Seidel relaxation
in Red-Black order [7]! Thus, this analysis leads us to expect efficiency that
is comparable to that obtained for the Poisson problem.
Example 1. We apply our algorithm to the advection-diffusion problem
with closed characteristics used in [10] (originally in [6]):
(0; 1), with u(x;
on @ and
For the advection term we use the same discretization as is used on the finest
grid in [6] and [10]: standard upstream (SU), defined by
where
In We use this value only at the stagnation point, adding no
viscosity elsewhere to maintain upstream discretization. Since the "physical"
viscosity is dominated by the artificial viscosity elsewhere anyway, the difference
is small. (Alternatively, we could use a much smaller ffl everywhere).
Levels Grid W= 1 W= 2 MGD9V

Table

1: The number of cycles necessary to reduce the L 2 residual norm by
a factor of 10 8 in Example 1. The problem is taken from [10], and MGD9V
is the automatic method of de Zeeuw which is used there.
As in [10], the initial guess for u
in\Omega is zero, and we cycle until the L 2
norm of the residuals is reduced to at most 10 \Gamma8 times its initial value. We
performed this test with 4, 5, 6, 7, and 8 levels, with the coarsest grid always
5 \Theta 5, including boundary points (as in [10]). We used V(1,1) cycles
throughout, with the usual full-weighted residual transfers and bilinear in-
terpolation. (See Section 5 for details on implementation of the relaxation).
In

Table

1 we compare our results with W=1 and 2 to those reported for
MGD9V, the automatic method of de Zeeuw (available up to six levels), using
a so-called "sawtooth" cycle with one ILLU (Incomplete Line LU) relaxation
sweep per level. It must be stressed that the efficiency and robustness of this
method is convincingly demonstrated in [10], and all the results achieved on
ten other problems (including non-recirculating advection diffusion) were far
better than these.
Evidently, the present method performs very well, with efficiency comparable
to that of elliptic equations in this simple test problem. Clearly,
the downstream relaxation ordering itself is not sufficient for the recirculation
problem (nor is ILLU). Both MGD9V and the present method with no
residual overweighting show clear deterioration as the grid is refined, while
with proper overweighting the convergence rate remains excellent even for
very fine grids. We reiterate that our residual overweighting approach may
not apply to MGD9V, since there different stencils are used on the different
grids, requiring different overweighting.
3.4 Several Bands of Residuals
It may not always be easy to obtain a single band of residuals per vortex. This
happens when the relaxation is carried out in piecewise-downstream ordering,
as would be the case in a domain-decomposition setting for example. Our
analysis can be extended to the case where several bands of nonzero residuals
remain after relaxation. It is found that one then requires two coarse-grid
corrections, with optimal weights 1 and 2 approximately. This seems to
imply that a W cycle is required in this situation. Also, one must then use
upstream intergrid transfers, so as to avoid averaging over interfaces between
subdomains, which may actually cause divergence.
4 Discretization
Flows in which the streamlines do not start and end at boundaries, but
constitute closed curves, require special considerations in the discretization.
In such cases, even a very small viscosity plays an important role in determining
the main flow throughout the domain. The solution in the limit
of vanishing viscosity depends very strongly on how these coefficients tend
to zero. In effect, the advection terms determine the behavior of the solution
along streamlines, whereas the viscous terms determine its cross-stream
form. And since the boundary is often a streamline itself, the propagation
of information from the boundary into the domain is governed by the viscous
terms no matter how small they may be. This effect is discussed in
detail in [3, 9], where it is shown for both the advection-diffusion problem
and the incompressible Navier Stokes equations that solutions with schemes
in which the numerical viscosity is anisotropic (having different viscosity co-efficients
for the cross-stream and streamwise directions), such as standard
upstream-difference schemes, may be spurious.
In the most general case it can be shown that even isotropic viscosity is
not sufficient for convergence of the solution, and one must actually specify a
uniform viscosity. We do not know how to do this while retaining the purely-
upstream structure (but see remarks in Section 6). However, for the homogeneous
advection-diffusion problem there are several indications (though no
proof) that isotropy suffices. This is shown below and also in [3, 9], where
it is also shown (in a numerical example) to suffice for the incompressible
Navier Stokes equations. This is consistent with the fact that the vorticity
in the Navier Stokes equations satisfies a homogeneous advection-diffusion
equation.
To obtain a discretization scheme that exhibits the appropriate physical-
like behavior for vanishing viscosity we must thus either add sufficient explicit
isotropic viscosity that will dominate the artificial viscosity of the discrete
advection operator, or else derive a discretization of the advection operator
that satisfies the condition of isotropy in its lowest-order truncated terms.
Since we want our scheme to remain purely upstream, we follow the latter
approach.
Consider the standard upstream scheme of (23), and assume for simplicity
of discussion a - b - 0 . From (2) we have by a Taylor expansion
~
Hence, in order to obtain isotropic artificial viscosity we may either add
some approximation of 0:5h(a \Gamma b)u xx , or else subtract some approximation
of 0:5h(a \Gamma b)u yy . In order to retain an upstream scheme we define this
additional viscosity at the point (i \Gamma
For general a and b we obtain in the first case the Isotropic-Viscosity
Upstream scheme IVU1, defined by
if ja i;j j ? jb i;j j, and
otherwise. In the second case we obtain scheme IVU2, defined by
if ja i;j j ? jb i;j j, and
otherwise. Here i1; j1 are defined as in (23), and similarly
The first truncated term in scheme IVU1 is thus \Gamma0:5 min(ja i;j j; jb i;j j)\Deltau,
while that of IVU2 is \Gamma0:5 max(ja i;j j; jb i;j j)\Deltau : Both schemes have isotropic
artificial viscosity, but that of IVU1 is smaller, and in fact it vanishes upon
alignment of the characteristic directions with the grid.
Both discretizations are stable in downstream-ordered Gauss-Seidel relax-
ation. The former is a nonnegative-weighted average of the standard first-order
and second-order upstream schemes, both of which are stable in this
relaxation. The latter produces an M-matrix. As expected, there were no
stability problems in any of our many numerical calculations.
5 Numerical Experiments
We first test numerically the discretizations derived in Section 4 on a model
problem for which the standard upstream scheme has been shown to yield
spurious solutions [3, 9]. Then, the asymptotic error reduction of two-level
and multilevel cycles are investigated for several problems.
5.1 Accuracy Test
The accuracy of the different discretizations is tested on the model problem:
a
with a and b given by
(These coefficients are the same as those of Example 2 below, and a picture
of the characteristics appears in Figure 1a). The domain of solution is the
unit square, centered at the origin, with a square of diagonal 0.5, whose
sides form a 45 degree angle with the axes, removed from its center. On the
outer boundary prescribed and on the inner boundary
solve this problem with the three upstream schemes-SU, IVU1, IVU2, and
also the (non-upstream) isotropic-viscosity scheme used in [3], denoted ISO.
viscosity is added: These solutions are compared to
that obtained with a standard second-order upstream scheme with physical
viscosity coefficient 0.001. The latter solution is obtained on a 257 \Theta 257 grid.
Grid SU IVU1 IVU2 ISO
129 \Theta 129 0.0750 0.0057 0.0163 0.0098

Table

2: L 1 difference norms between the solutions obtained with several
schemes at different resolutions, and a high-accuracy solution obtained on
grid 257 \Theta 257 (see text for details). The three isotropic-viscosity schemes
are seen to yield convergent solutions, but not the standard upstream (SU)
scheme.
In

Table

2 we present L 1 norms of the differences between the test solutions
at various resolutions, and the second-order accurate solution (restricted to
the corresponding grid by injection). Since the solution is smooth, and the
physical viscosity dominates the second-order truncation terms at this
high resolution, the latter solution is assumed to be very accurate.
Evidently, the three schemes with isotropic artificial viscosity produce
convergent solutions but the SU scheme does not, despite the fact that its
average viscosity is smaller than that of IVU2 and ISO. IVU1, which has
the least artificial viscosity of the four schemes tested, produces the smallest
error.
5.2 Efficiency Tests
The remainder of our numerical calculations are aimed at testing the performance
of the algorithm in various configurations and comparing to the
analytical predictions of Section 3. We test the SU scheme (as this is the
most widely used first-order scheme) and the IVU1 scheme (which is more
accurate than IVU2 and also employs just a four-point stencil). In all these
tests we use first-order upstream residual restriction and bilinear interpolation
of the corrections. The restriction is performed as follows: for all
even i and j on the fine grid we define i1, j1 as in (23), and restrict to the
corresponding coarse-grid right-hand side at coarse-grid point (i=2; j=2) the
average of the fine-grid residuals at points (i; j), (i1; j), (i; j1), and (i1; j1).
This restriction gives slightly better results than standard full weighting in
multi-vortex problems, since residuals are less likely to be transferred from
one vortex to another.
The finest grid in all the tests is 129 by 129, and six levels are employed
except in the two-level tests. We include no "physical" viscosity except at
stagnation points it is required for well-posedness. We
calculate convergence factors as follows: the boundary conditions and right-hand
sides are chosen to be zero; (the choice is immaterial for linear problems,
but this allows us to normalize the solution by a constant factor every few
cycles in order to avoid roundoff errors). The initial solution field is pseudo-
random, and 100 cycles are performed. We calculate the convergence factor
as the (geometric) average error-convergence per cycle over the last 80 cycles.
The averaging is used because in some cases the convergence history is not
smooth, and the value corresponding to any particular cycle may not carry
much meaning. However, it should be noted that in all cases the convergence
factors in the vicinity of the optimal W were not sensitive to the exact choice
of W .
Example 2. The first test for efficiency is the problem used above to test
the discretization, but without the inner "island" (so as to have pure recirculation
everywhere). The characteristics of this problem, which form a single
clockwise-rotating vortex, are plotted in Figure 1a. A relaxation sweep is implemented
by sweeping four times over the domain, and in each such sweep
relaxing roughly one quarter of the variables as follows: in the first sweep
only variables at locations where both a(x; y) and b(x; y) are nonnegative are
relaxed (designated first quadrant); in the second sweep only variables corresponding
to locations where a(x; y) is nonnegative and b(x; y) is nonpositive
are relaxed (second quadrant); in the third sweep only variables at locations
where both a(x; y) and b(x; y) are nonpositive are relaxed (third quadrant);
in the fourth sweep only variables corresponding to locations where a(x; y)
is nonpositive and b(x; y) is nonnegative are relaxed (fourth quadrant). Fi-
nally, all stagnation points (in this case just one) are relaxed last. This entire
process comprises a single clockwise sweep. Of course, each quarter-sweep is
performed in downstream order i.e., x and y increasing in the first quarter-
increasing and y decreasing in the second, etc. It is efficient to store
the order of relaxation of the entire sweep during a setup-sweep (which costs
very little), so that from then on each full clockwise sweep costs nearly the
same as an ordinary lexicographic Gauss-Seidel sweep.
Note that we specify nonnegative and nonpositive in the above descrip-
tion. That is, it should be ensured that the boundaries of the quadrants are
included in the quadrant. This was important in some of the tests.
Such a clockwise sweep indeed eliminates all the residuals except along
a narrow band that extends from the center of the vortex to the boundary.
But if the vortex rotates counterclockwise, several such bands would remain.
Hence, in Examples 3-5 below, where both clockwise and counterclockwise
vortices exist (as would be the general case), we perform an analogous counterclockwise
sweep following each clockwise one. The quadrant where this
counterclockwise sweep starts has some bearing on performance. In our experiments
we use exactly the reverse order. This allows us to save some work
by performing the counterclockwise sweep over just three quadrants, since
the fourth quadrant (a(x; y) nonpositive and b(x; y) nonnegative) has just
been relaxed in the clockwise sweep. Thus, we begin with the third quadrant
(a and b nonpositive), then the second, and finally the first quadrant. (Note,
by the way, the difference from the usual symmetric Gauss-Seidel: here the
sweeps are all performed in downstream ordering within each quadrant. But
the quadrants are scanned symmetrically).
In Problems (3-5) (especially there was also some sensitivity to which
quadrant one chooses to start the relaxation sweep. That is, a somewhat
different convergence rate and optimal overweighting were obtained if one
performed the clockwise sweep as above, than if one relaxed, say, the second
quadrant first, then the third, then the fourth, and finally the first. Hence,
in these tests we shifted the relaxation starting-point by one quadrant after
every cycle so as to obtain results that represent some average or "typical"
case.
We first performed two-level tests in order to compare the numerical
results with the analysis. We used and the theoretically optimal W ,
calculated from (19). This is the relevant value for this example, since there
is a one-line overlap between the quadrants. The results are given in Table 3.
Evidently, the analysis captures the main features of this problem very well,
despite the numerous simplifications, as seen especially in the SU results.
Calculations were also carried out with V(1,0), V(1,1), and V(2,1) cycles,
with and the optimal value, which was determined experimentally.
The numerical results are summarized in Table 4. We find that the experimental
V-cycle results also match the two-level predictions fairly well in the

Table

3: Comparison of error-amplification factors obtained by analytical
prediction and two-level numerical calculations of Example 2 for
W opt . The optimal values of W were obtained from (19).
vicinity of the optimal residual-weighting factors. But the optimal W 's are
somewhat higher than predicted, although they do show the expected dependence
on the number of relaxation sweeps performed. In other problems,
reported below, the optimal value varied, but it was always fairly close to 2.
The overall performance when the optimal W is used is very satisfactory.
As noted above, this performance is not sensitive to moderate changes in W
(see also below).
Example 3. This example features flow with four vortices rather than just
one. Here, a and b are given by
and the characteristics are plotted in Figure 1b. The domain, as in all the
examples except Example 5, is the unit square centered at the origin. The
numerical results using V(1,1) and V(2,1) cycles appear in Tables 5 and 6,
respectively. The convergence performance remains excellent, even though
some nonvanishing residuals remain after the relaxation at parts of the borderlines
between vortices (where the flow leads away from the borderline).
Recall, however, that here and below each full relaxation sweep consists of
one clockwise sweep followed by three quarters of a counterclockwise sweep
(see description of implementation above), in order to allow for the opposite-
sign vortices.
Cycle W SU IVU1
V(1,0) 1. 0.795 0.898
V(1,1) 1. 0.676 0.831
V(2,1) 1. 0.536 0.724
2. 0.280 0.440
2. 0.143 0.302
2. 0.069 0.133

Table

4: Error-amplification factors obtained in numerical calculations of
Example 2. The optimal values of W were obtained experimentally (see
text).
Example 4. In order to test the effect of grid-alignment of the borderlines
between vortices we solve a problem in which this borderline is not aligned
with the grid. In this problem a and b are given by a 1
respectively, where
a
with
This represents a superposition of two opposite-sign vortices. The characteristics
are depicted in Figure 1c, and the numerical performance is shown in

Tables

5 and 6. Here we see some loss of efficiency, but the performance is
still satisfactory and far better than is usually exhibited in such problems.
Example 5. Here we test a mixed problem, where the flow enters and
leaves through the boundary, but there also exists a large recirculation zone.
This is obtained by redoing Example 2, but in an extended domain:
Example SU IVU1

Table

5: Error-amplification factors obtained with V(1,1) cycles for
and W opt . The latter were found experimentally.
Example SU IVU1

Table

obtained with V(2,1) cycles for
and W opt . The latter were found experimentally.
The mesh is still uniform, and the finest grid is now 193 by 129. The characteristics
are shown in Figure 1d, and the numerical performance is given in

Tables

5 and 6. As in the other examples, there was virtually no sensitivity to
the order in which the relaxation sweeps were performed (that is, clockwise
first and then counterclockwise or vice versa).

Summary

Often, one may not wish to search for optimal residual weighting
factors for every problem. Instead, one can simply use the nominal value
2. In the "realistic" examples 3-5, one saves at most 25% of the time
spent in relaxation by using the optimal value rather than 2, and usually
much less. This would also be the case in Example 1 if we were to employ
the "symmetric" relaxation. Even with this nominal value the convergence
rates are comparable to those of elliptic problems.
6 Conclusions and Further Research
An experimental approach that is hoped to eventually lead to a fully efficient
solver for general high-Reynolds flows has been introduced, analyzed and
tested on the advection-diffusion problem in the inviscid limit. The numerical
tests mostly match the predictions very well, indicating that the main cause
for slow convergence of the usual multigrid algorithms for recirculating flows
has indeed been understood, and a way to eliminate it has been found.
The multigrid V-cycle, using downstream relaxation and upstream dis-
cretization, was shown to yield an efficient solver for the tested problem
in several simple situations of closed characteristics and in a mixed en-
tering/recirculation problem. The tests were performed with the classical
standard first-order upstream discretization scheme and also with a novel
first-order upstream discretization, that was shown to preclude the spurious
solutions reported in [3].
The present approach is cheaper to implement than that developed in [5],
and can straightforwardly be applied to mixed entering/recirculating flows.
More important, there is potential of success with high-order discretization,
for which the approach of [5] yields an inadequate compromise. However,
the results obtained are still preliminary. The effect of the intergrid transfers
on the small band of residuals and its consequences in terms of error-
reduction efficiency should be investigated over a wide variety of cases, along
with a study of how to deal with (or avoid) situations where there remain
several bands of nonvanishing residuals per vortex. Then, further research
should be directed towards higher-order discretization. For this case too,
an effectively-upstream discretization needs to be developed, whose truncation
error represents isotropic artificial diffusivity. One approach is to use a
predictor-corrector type discretization, employing an upstream scheme as a
(local) driver and a possibly higher-order (not necessarily upstream) scheme
as a (local) corrector. Finally, the present approach has only been tested on
the advection problem. Experiments with the incompressible Navier Stokes
equations for flows with closed streamlines need to be performed, employing
distributive Gauss-Seidel relaxation, as shown in [4]. These will no doubt
raise further questions.
It is anticipated that the techniques investigated here will carry over to
three dimensions, although the implementation will be considerably more
complicated. This is supported by the fact that simple experiments performed
with the overweighting methods of [5] in three dimensions exhibited
the expected performance.
An obvious drawback of the entire approach is that it is inherently sequen-
tial, and efficient parallel implementations are hard to envisage. Some parallelization
might be achievable by performing a downstream line Gauss Seidel
relaxation. Also conceivable is a domain decomposition approach which
leaves several lines of residuals per vortex.
Another drawback of the present approach is that it is not directly applicable
for flows with significant additional viscosity, since this entails using
discretizations that are not purely upstream. Methods that deal with such
flows as well are presently being investigated.

Acknowledgment

This work was supported by The Royal Netherlands Academy of Arts and
Sciences and The Feinberg Graduate School, by the United States-Israel
Binational Science Foundation under grant no. 94-00250, and by The United
States Air Force Grant F49620-92-J-0439 and by the Carl F. Gauss Minerva
Center for Scientific Computation.



--R

"Multigrid Solvers for Non-Elliptic and Singular- Perturbation Steady State Problems,"
"1984 Multigrid Guide with Applications to Fluid Dynam- ics,"
Inadequacy of First-Order Upwind Difference Schemes for Some Recirculating Flows
On Multigrid Solution of High-Reynolds Incompressible Entering Flows
Accelerated Multigrid Convergence and High- Reynolds Recirculating Flows
Efficient Solution of Finite Difference and Finite Element Equations by Algebraic Multigrid
fundamental algo- rithms
On the Correctness of First Differential Approximation of Difference Schemes
Multigrid Techniques for Incompressible Flows

--TR
recirculating flow;upstream discretization;advection-diffusion;multigrid