Generally, fast direct solvers are not directly applicable to a nonseparable elliptic partial differential equation. This limitation, however, is circumvented by a semi-direct procedure, i.e., an iterative procedure using fast direct solvers. An efficient semi-direct procedure which is easy to implement and applicable to a variety of boundary conditions is presented. The current procedure also possesses other highly desirable properties, i.e.: (1) the convergence rate does not decrease with an increase of grid cell aspect ratio, and (2) the convergence rate is estimated using the coefficients of the partial differential equation being solved

Chang, S. C.

English

NASA Technical Reports Server

# 8 4 - 2 5 9 4 6
~ * - 7
NASA Technical Memorandum 83704
A Semi-Direct Procedure Using a Local
Relaxation Factor and Its Application
to an Internal Flow Problem
Sin-Chung Chang
Lewis Research Center
Cleveland, Ohio
Prepared for the
Ninth International Conference on Numerical Methods in Fluid Dynamics
cosponsored by CEA, CISI, DRET, CNES, EOF, Framatome, and Novatome
Saclay, France, June 25-29, 1984
NASA
https://ntrs.nasa.gov/search.jsp?R=19840017878 2020-03-20T21:49:04+00:00Z
A SEMI-DIRECT PROCEDURE USING A LOCAL RELAXATION FACTOR
AND ITS APPLICATION TO AN INTERNAL FLOW PROBLEM
Sin-Chung Chang
NASA Lewis Research Center
Cleveland, Ohio, 44135 U.S.A.
1. INTRODUCTION
Generally, fast direct solvers (FDS)* are not directly applicable .to a non-
separable elliptic partial differential equation (PDE). This5 limitation, however,.
can be circumvented by a semi-direct procedure, i.e., an iterative procedure using
FDS algorithms. In this paper, we present. an efficient semi^Jirect procedure which
is easy to implement and applicable to a variety of boundary conditions. The cur-
rent procedure also possesses other highly desirable properties, i.e., (1) the con-
vergence rate does not decrease with an increase of grid cell aspect ratio, and (2)
the convergence rate can be estimated using the coefficients of the PDE being solved.
Many elliptic PDE encountered in flow problems can be expressed as
Lu = h (1.1)
where L is a nonseparable second-order linear elliptic operator, u the dependent
variable, and h a given source term. Equation (1.1) may be solved with the
iterative procedure:
 :
v2 (unn - un) = -T (Lun - h) (n = 0, li 2,...) (1.2)
Here T is a nonzero relaxation factor and v2 the Laplacian which can be directly-
inverted by a FDS. In the previous works2"4 using Eq. (1.2),
 T is treated as. a
constant and the iteration is accelerated by an optimal choice of T. In this study;
we consider T as a spatially varying parameter. The local value of T is chosen to
optimize the local convergence rate which is evaluated using a von Neumann analysis.
2. STABILITY ANALYSIS .- •
As ah initial step, it is assumed that T is a constant and
• - . " • ' • 2 2 2
L - a ~ •+ 2b -*— - +. c ~ (a > o, c > o, and ac-b2 > o) (2.1.)
ax^ 3X3y ay^ , .. -
where a, b, and c are constants. Using a uniform grid with grid intervals AX and
Ay, the central difference form of Eq. (1.2) at a grid point (i,j) is
_
Here hy is the source term and the operators P and L, respectively, are
defined by - ''
: . . . '(2.1)
The convergence rate of Eq. (2.2) will be' analyzed assuming ;
en. . = u^ .-u. = enexp[I(ik Ax-^jk Ay)],
 : ' I =vQ , (2.5)1
 » j ! > j i x
where ui^j is the converged solution, ej^ the error amplitude and k = (^x^y)
the propagation vector. Assuming that k # 0 and using the following definitions:
ex = (kxax)/2, ax = (sinex)/Ax, ey = (kyAy)/2, ay = (siney)/Ay (2.6)
Sy , «y[(«x)+(«y)] (2.7)
then Eqs. (2.2) and (2.5) imply that
en+1/e" = E(T,!C) = l-TG(lc) (1c * 0) (2.8)
with 6(lc) = a(s J2 + c(s J2 + 2bs s cosev coseu (2.9)A y A y A y
It should be noted that the von Neumann analysis (along with the condition £ £ 0) can
be fully justified if it is applied to a grid with periodic boundary conditions.
Let omax and omin, respectively, be the greatest and tfte smallest eigen-
values of the symmetric and positive-definite matrix (u, v = 1, 2)
A
 = ("pv^' an=a» «i2=«2i=b' B22=c (2-10)
Then it can be shown that^
"max > 6 <*) > °min > 0 (2.11)
with the understanding that the bounds omax and om,n are sharp if all the
allowable T(^0) are considered. As_a result, the asymptotic error amplification
factor E(T ) (= the supremum of |E(T, 1T)| for a given t) will reach its minimum
E0 = E(T Q ) = (I - l ) / ( z + 1) < 1 (2.12)
when T = T =the optimal relaxation factor = 2/(o + a . ) = 2/(a + c) (2.13)0 ' max min
Here E is the condition number =omax/arojn >1.
Equations (2.1)-(2.13) can be easily generalized for N-dimensional problems.
The only exception is Eq. (2.13) where the last equality sign is not valid if N>2.
3. LOCAL RELAXATION
The variable coefficient (VG) version of the numerical procedure presented in
Section 2 is obtained by replacing the constant coefficients a, b, c, T, TO, omax,
and o™^n in Eqs. (2.2) to (2.4) and (2.13), respectively, with the grid point
dependent coefficients a^, b^, c^, t^ ', TOJJ, omaX)ij and ominjij. Obviously,
the VC version is well defined at all interior grid points. For a'grid point on a
periodic or Neumann boundary, it can also be defined using the periodic condition or
an extrapolation technique explained in Ref. 6. —
The current procedure can be modified to solve PDE with
L = L'(x.y) H ly (p(x,y) ij) + ^  (q(x,y) fy) (3.1)
where p and q are arbitrary positive functions of x jjnd y. With L = L'(x,y),
the VC version of Eq. (2.2) can be obtained by replacing T(u°tj) with'
~
1 2
 l,j+ P(1+1/2)J u1+l,j-(p(i-l/2)J+P(1*l/2)j)u1.jl +
?.j] (3'2)
In case that the values of p and q do not vary greatly from one grid point to its
neighbors, then L'(u? .)*L(u!? .) assuming a . . = p.., c. . = q.. and b- . = 0. Thisi,j i., j ij i J ' J 'J 'J
observation coupled with Eq. (2.13) leads us to the assumption:
(L -L ' ( x .y ) ) (3.3)
4. SCALING
Assuming that the definition of the operator P is broadened as
with gx and gy positive constants, then the only modifications required for Sees. 2
are to replace the coefficients ax, ay, a, b, and c in Eqs. (2.7), (2.9), (2.10)
and (2.13), respectively, with
 v/gx~ax, V§yay> a/9x» b/v/gTgy~and c/gy. Since
gx and gy are positive constants, their appearance in Eq. (4.1) does not increase
trie difficulty of inverting P. However, their introduction into the current itera-
tive procedure does have an effect on the convergence rate. As a result, the current
iterative procedure can be accelerated considerably by a proper xhoice of g* and
gy. Obviously, this scaling technique can also be used on the solution of PDE with
variable coefficients as long as the scaling coefficients gx and gy remain
positive constants.
5. NUMERICAL EVALUATION
To facilitate this discussion, we begin with the following preliminaries:
(a) the residual norm r(n) and error norm e(n) after n iterations are defined by:
r(n) = [£ (L(u" .)-h )2]1/2 and e(n) = [ £ (u" ,-u' J2]1'2 (5.1)
-ii • » j ' j { • j ' l j ' f j • . •1
 » J ' > J ' ' . . • . .
whe"-e u? . is the machine accuracy solution. Moreover, we define • :1
 » J ' •" • . . • ' •
Or(n) = -Iog10 [r(n)/r(0)] and Oe(n) = -Iog10 [e(n)/e(0)]
(b) We assume that u. . =^0 at all grid points where u^ j's.are unknowns.
 n
As a result, e(o) ^llumll and 0P(n) is the number of correct digits in u., - ...
' i J(c) The coefficient Ep ij defined by the VC version of Eq. (2.12) may be inter-
preted as an asymptotic error amplification factor. As a result, the value of
Oe(n) (and Or(n) if Or(n) - Oe(n)) may be predicted by the parameter
Ot(n) = -n logio (the supremum of all E0>ij) (5.3)
where EO.IJ ^s evaluated assuming a^j = p^j, by = 0 and cy
if L = L'(x,y). Alternatively, the value of Oe(n) or Or(n) may be
predicted by 0^ (n) = the limit of 0^(n) as the grid cell size approaches zero.
(d) The domain for all numerical problems is assumed to be l>x>0 and .l>y>0. Moreover,
the iterative procedure is executed using aan FDE code8 with 'gx = gy = 1.
The first group of PDE to be studied in this section includes:
(l+2x2+2y2)32u/3x2+(l+x2+y2)32u/3y2=l (5.4)
y2)32u/3y2=l (5.5)
1=2,4,6,8 (5.6)
Ten numerical problems associated with' Eqs. (5.4)-(5.6) are defined in Table 1. The
parameters MX and MY, respectively, are the number of grid intervals in the x and
y directions. It is assumed that u =0 at all boundaries.
Problems 1-3 are associated with Eq. (5.4). They differ only in the grid cell
size and aspect ratio. As shown in Table 1, for all three problems, the value of
Ot(20) or Ot(20) provides a conservative estimate of Or(20) (since Oe(n) ~ Or(n),
the values of Op(n) are not listed). It is also seen that the very large aspect
ratio (16:1) of the grid cell in problem 3 does not cause any substantial reduction
in convergence rate. Problems 4-6 represent a parallel study for Eq. (5.5). It is
3
evident that the appearance of a cross derivative term in Eq. (5.5) causes_an
increase in the conservatism of the theoretical estimate- made by Ot(n) or Ot(n).
This is particularly true if the grid cell size or aspect ratio is large. For
problems 1-6, the convergence histories are closely represented by straight lines (as
an example, see Fig. 1) as would be expected from our theoretical study.
Equations (5.6) belong to the class of PDE defined by Eq. (3.1). It is evident
that the variation of the coefficients p and q increases progressively as one
goes from i = 2 to i = 4 and so on. For t = 8, the increase in the values
of p and q from one corner (x=y=0) to another corner (x=y=l) on the unit square
is of the order of 100 times. As shown in Fin. 2, the current method is still useful
in this extreme case. The robustness of this algorithm is most evident in its
ability to reverse the trend toward divergence during the first few iterations.
The second group of PDE to be studied includes
•
.*) (5.7)
x.y) (5.8)
Here hj(x,y) and (^ (x.y) along with the boundary values of u are chosen such that
u = "l(x,y) = sinx siny and u = U2(x,y)=[x(l-x)y(l-y)]2 (5.9)
respectively, are the exact solutions of Eq. (5.7) and (5.8). Both Eqs. (5.7) and
(5.8) are solved numerically assuming MX=MY=16. Both numerical problems are test
problems used by. Bank4. Compared with the current value of Oe(10) = 3.47, the
values obtained by Bank for Eq. (5.7) are 3.49 without using his scaling technique
and 3.87, 4.81, and 6.76 using different scaling functions. Since the operator used
by Bank is a general separable operator, the computation time required for its
inversion is about six times that required for the current Poisson-like operator (p.
5 of ref. 1). As a result, the current method is at least a factor of 2 more
efficient even without using our scaling technique. For Eq. (5.8), compared with the
current value of Oe(10) = 8.59, the values reported by Bank are 5.88 without using
scaling technique and 14.79 if a scaling function is used. Equation (5.8) was also
solved by Concus and Golub3 and the results are comparable with ours. However
their method is applicable only when p(x,y) = q(x,y) as in the case of Eq. (5.8).
6. APPLICATION TO A 3-D INTERNAL FLOW PROBLEM
In a numerical study (to be published) of 3-D incompressible shear flow in a
turning channel (Fig. 3), the elliptic PDE
= h, (2>x>-2, 0.75>y>0.65, 0.1>z>0) (6.1)
is solved in computational space. Here h is a known source term and
n(x,y) = (cosh(ux) + cos(ny))/(cosh(i.x) - cos(uy)) > 0 (6.2)
It is assumed that u = 0 at x = 2 and the normal derivative of u vanishes at
other boundaries. Since Eq. (6.1) has a form specified by Eq. (1.1) and the 3-D
version of Eq. (2.1), it was solved numerically using the current iterative method
with a. 144 x 12 x 12 uniform grid. To accelerate the convergence, the scaling
coefficients of the 3-D operator P (see Eq. (4.1)) are given by
9x = 9y - *• '3 " A,ax • * * °'4135 (6'3)
where r^ ,— (=0.9966) and ^n (=0.1716), respectively, are the maximum and the
minjmum of n in the domain of Eq. (6.1). It should be noted that the maximum value
of Ot(l) (= 0.415) is reached if the values of gx, gy, and gz are chosen
according to Eq. (6.3).
Equation (6.1) was numerically solved many times with different source terms h.
Without exception, it was found that, after four iterations, the slope of Or(n)
vs n curve is predicted by the value of 0^ (1) to within ten percent.
REFERENCES
1. Hockney, R. H.: "Rapid Elliptic Solvers," Numerical Methods in Applied Fluid
Dynamics edited by B. Hunt, Academic Press, 1-48 (1980).
2. D'yakonove, E. 6.: Dokl. Akad, Nauk, SSSR, 138, 522-525 (1961).
3. Concus, P. and Golub, G. H.: S1AM J. Numer. Anal., 10, 1103-1120 (1973).
4. Bank, R. E.: SIAM J. Numer. Anal.. 14. 950-970 (1977).
5. Chang, S. C.: "Generalizations of Two Inequalities Involving Hermitian Forms,"
accepted for publication in Linear Algebra and Its Applications.
6'. Smith, 6. D.: "Numerical Solution of Partial Differential Equatfons," Oxford
University Press, p. 29 (1978).
7. Hageman, L. A. and Young, D. M.: "Applied Iterative Methods," Academic Press,
p. 12 (1981).
8. Adams, J., Swartztrauber, P. and Sweet, R.: "FISHPAK: A Package of FORTRAN
Subprograms for the Solution of Separable Elliptic Partial Differential
Equations, Version 3," National Center for Atmospheric Research, Boulder,
Colorado, June 1979.
Table 1. (n=20 for Probs. 1-3, n=32 for Probs. 4-10)
Prob.
No.
1
2
3
4
5
6
7
8
9
10
Eq.No.
(5.4)
(5.4)
(5.4)
5.5)
(5.5)
(5.5)
(5.6)
(5.6)
(5.6)
(5.6)
I
NA
NA
NA
NA
NA
NA
2
4
6
8
MX
16
64
64
16
64
64
16
16
16
16
MY
16
64
4
16
64
4
16
16
16
16
Or(n)
14.5
13.8
14.0
14.4
12.6
19.0
17.2
16.8
13.4
7.6
Ot(n)
12.3
12.1
12.7
9.7
9.6
10.0
15.3
15.3
15.3
15.3
Ot(n)
12.0
12.0
12.0
9.6
9.6
9.6
15.3
15.3
15.3
15.3
I I I I I I I I I I I I I I I
Figure 1. - Convergence histories of Prob-
lems 4-6.
Figure 2. - Convergence histories of Prob-
lems 7-10.
Figure 3. - A converging-diverging
turning channel in physical space.
1. Report No.
NASA TM-83704
2. Go intent Ac stonNo. 3. Recipient's Catalog No.
4. Title and Subtitle
A. Semi-Direct Procedure Using a Local Relaxation
Factor and Its Application to an Internal Flow
Problem
5. Report Date
6. Performing Organization Code
505-31-08
7. Authors)
Sin-Chung Chang
8. Performing Organization Report No.
E-2162
10. Work Unit No.
9. Performing Organization Name and Address
National Aeronautics and Space Administration
Lewis Research Center
Cleveland, Ohio 44135
11. Contract or Grant No.
12. Sponsoring Agency Name and Address
National Aeronautics and Space Administration
Washington, D.C. 20546
13. Type of Report and Period Covered
Technical Memorandum
14. Sponsoring Agency Code
15. Supplementary Notes
Prepared for the Ninth International Conference on Numerical Methods in Fluid
Dynamics, cosponsored by CEA, CISI, DRET, CNES, EOF, Framatome, and Novatome,
Saclay, France, June 25-29, 1984.
16. Abstract
Generally, fast direct solvers are not directly applicable to a nonseparable
elliptic partial differential equation. This limitation, however, can be cir-
cumvented by a semi-direct procedure, i.e., an iterative procedure using fast
direct solvers. In this paper, we present an efficient semi-direct procedure
which is easy to implement and applicable to a variety of boundary conditions.
The current procedure also possesses other highly desirable properties, i.e.,
(1) the convergence rate does not decrease with an increase of grid cell aspect
ratio, and (2) the convergence rate can be estimated using the coefficients of the
partial differential equation being solved.
17. Key Words (Suggested by Author(s))
Numerical analysis
Fluid mechanics
1* Security CUsalf. (of this report)
Unclassified
18. Distribution Statement
Unclassified •
STAR Category
20. Security Classlf. (of this page)
Unclassified
- unlimited
34
21. No. of pages 22. Price*
'For sate by the National Technical Information Service, Springfield. Virginia 22161
National Aeronautics and
Space Administration
Washington, D.C.
20546
Official Business
Penalty for Private Use, $300
SPECIAL FOURTH CLASS MAIL
BOOK
Pottage and Fees Paid
National Aeronautics and
Space Administration
NASA-451
NASA POSTMASTER: If Undeliverahle (Srrtion IS*Postal Manual) IVi Not Return


A semi-direct procedure using a local relaxation factor and its application to an internal flow problem

Abstract

Similar works

Full text

Available Versions

NASA Technical Reports Server