Analysis of Preconditioning and Relaxation Operators for the Discontinuous Galerkin Method Applied to Diffusion

11 Pages·2001·0.67 MB·English

by  Atkins, H. L.

AIAA-01-2554
ANALYSIS OF PRECONDITIONING AND RELAXATION OPERATORS FOR THE
DISCONTINUOUS GALERKIN METHOD APPLIED TO DIFFUSION
ti. L. At,kins*
NASA Langley Research Center
[tampton, VA 23681-0001
('hi-Wang Shut
Brown University
Providence, RI 02912
Abstract iterative methods that are local to the element and
The explicit stability constraint of the discontinn- preserve the compact character of the. discontinuous
ous Galerkin method applied to the diffusion operator Galerkin spatial discretization. One obvious choice is
a block Jacobi relaxation scheme in which some or all
decreases dramatically as the order of the method is
increased. Block Jacobi and block Gauss-Seidel pre- of the local contributions to the spatial operator are
conditioner operators are examined for their effective- treated implicitly.
ness at accelerating convergence. A Fourier analysis This paper presents the Fonrier analysis of block Ja-
for methods of order 2 through 6 reveals that both cobi and block Gauss-Seidel preconditioning operators
preeonditioner operators bound the eigenwdues of the for the discontinuous Galerkin method applied to the
discrete spatial operator. Additionally, in one din|en- diffusion operator. Data froxn the analysis is used in
sion, the eigenvalues are gronped into two or three the development of optimal relaxation methods that
regions that are invariant with order of the method. rapidly damp the high frequencies, and are well suited
Local relaxation methods are constructed that rapidly for use with multi-grid r techniques. Numerical tests
damp high frequencies for arbitrarily large time step. are presented that veri_" the analysis,
The first, section describes the discontinuous
Introduction
Galerkin method as applied t,o diffusion, introduces
The discontinuous Galerkin method allows coral)act required notation and formulates the partict|lar im-
spatial discretizations of any order to be formulated plicit problem that is the focus of this paper. The
on arbitrary meshes. The compact discretizations are second section presents a detailed Fourier analysis, in
well suited to explicit time-marching methods for un- one dimension, of the preconditioning operators and
steady simulations, and are easily and efficiently im- the findings of that analysis. The third section for-
pIeniented to run in a parallel environ|ne||t, t When mulates several optinlal relaxation methods based on
applied to propagation and advectioi! equations, the data from the analysis, and demonstrates them in nu-
stability constraint decreases only ntoderately 2as the merical tests. The fourth section extends the analysis
order of the method is increased. However, when ap- to two dimensions and fornmlates optimal relaxation
plied to the diffusion operator, the explicit stability methods that are demonstrated in nnmerical tests.
limit decreases dramatically as the order of the method
is increased, a and some degree of implicitness is nec-
Methodology
essary'.
In this sect.ion, we describe in general terms the ap-
hnplicit time marching methods ranging from back-
ward Euler to a high-order implicit IRunge-Kutta coukl plication of the discontinuous Galerkin method to a
model diffusion eqt|ation of the form
be employed, depending on the level of temporal accu-
racy required. Regardless of the approach, a globally
c%i
implicit solution is required. The direct approach, a
0-7-- V- [_,V,,] = 0. (t)
direct solver, destroys the compactness inherent t.o the
discontinuous Galerkin method along with the associ-
Although there are several ways in which the discontin-
ated advantages. Thus we are motivated to focus on
uous Galerkin method may be applied to a diffusion
operator 3-5 this analysis will focus on an approach
*Senior member, Senior Hesearch Scientist, (%mputational
Modeling and Simulation Branch that requires a first order equation. To this end, equa-
tPro%ssor, Division of Applied Mathematics tion (1) is recast as a set. of first order equations:
Copyright @ 2001 by the American Instituteof Aeronautics
and Astronautics, Inc. No copyright is asserted in the United
States under Title 17. U.S. Code. The IJ.S. G,)vernment has
0-7-v.¢ = 0
a roya[ty-fl'ee license to exercise all rights uncler the copyright
claimed herein for Governmental Purposes. All other rights are
,f-/*Vu = 0.
reserved by the copyright owner.
1
AMEBICAN INSTI'FIVI'E OF AERONAUTICS AND AST,qONAI!TICS
AIAA-01-2554
Spatial Discretization
The discontinuous Galerkin method is applied on a
LOfh,k "
discrete element fh obtained by partitioning the do-
main into arl)itrarily shaped elements. The solution
in each dement is approximated in terms of a set of
local basis functions {bl,i} as il _ tq =_ _{i} bt,iul, i LOF6,k J
and _ _ = Ell} bl,iqt,i ill _l. The governing equa- a/ld
tions are projected onto each men]bet of the basis set.
to give a set of equations in each element that govern
the evolution of the set of unknown variables {ul,,}.
To facilitate tile analysis, it. is assumed that the do-
main has been uniformly partitioned such that the
abt,, - V.(_ d_ = 0 matrices A, Bk and C_. are the same for all elements.
Further more, in this study we will examine only the
fbt,i [_ - V'ut]da = 0. (3) so called "mixed LR" scheme 3'5 in which 0t.k is either
zero or one on each boundary segment O_t,k, and @,k
is assigned in a regular and repeating pattern.
The gradient term is eliminated by integrating by parts
to obtain the weak form Temporal Discretization
In previous works dealing with propagation, -_,s ex-
+E ___o
plicit Runge-Kutta was used to advance the solution
in time. ttowever, the stability bound of such an
explicit, method is inversely proportional to the maxi-
/ [bt,i(_- V'bl,i'ut]dQ+ E ./bt,ih(ut, uk)c['s=o (4) nmm eigenvalue of the spatial operator. When discon-
tixmous Galerkm is applied to the diffusion operator,
the maximum eigenvalue grows in magnitude approx-
where h denotes a numerical flux, the subscript k de-
imately proportional t,o p4 and the explicit approach
notes the index of aneighboring element, 0Qt,k denotes
becomes impractical for most problems. An implicit
the boundary between fit and Qk, and c_ is the out-
time marching method would clearly remedy the prob-
ward boundary normal Note that the variables {_}
lem; however, a direct, inversion of the resulting matrix
are simply intermediate values that can be computed
woukt destroy the compactness of the discontinuous
directly and explicitly from the solution. It| the case of
Galerkin. Thus in the present work, we consider im-
advection and propagation, the numerical flux is usu-
plicit evolution methods that are solved by local relax-
ally an approximate Riemann flux. For diffnsion, the
ation methods, and we are interested in the analysis
flux is simply a weighted average of u or 0"such that
of the relaxation process.
fi((_, Ck) : (1 - _t,_.)q;+ *t,_,£-, Given a general evolution equation of the form
h(u_,uk) = @,k< + (1 - Ol,k)Uk,
-- = _(,,)
/-)l
and 0t,k = l - 0k,t. Evaluation of the integrals and
fluxes of equation (4) leads to a semi-discrete evolution
where R(u) denotes a spatial operator, all implicit
equation of the form
methods involve solving at] equation of the form
OU_
- A. Q, + _ [(1- ,,._)B_. q, + o,._c_, q_]
ot u"-uAt"-_ -- (H(A'fzfu)2") _-H (tt.-1 ,u-.n-" .... (6)
{k} n(,,n-_), n(,,"-'-') .... ),
q, = AU, + y_ [,t knout + (1- o, k)C,Uk]. (.5)
{k} where tile superscript n may denote either the time
step (t = n dx t) or a step within the R||nge-l(utta
where sequence, and ct and H depend o/l the specifics of the
method. In the present analysis, we are concern(,d
ffl, 1
o,fly in the construction and stability of an iterat.ive
Ill, 2
Ut -_ Ql = method for the solution _d_ it] equation (6). Thus, we
III, 3 ,
assume the discrete evolution equation (6) is stable in
tin]e, possesses the desired order properties, and that
old values of the solution u'- _, u'-2 ..... are known.
Within lhe scope of a Fourier analysis, the terms
L'_ involving the old vahles of the solution have no etfe('t
2
AMEFtI('AN [NSTITUrFIq OF AEIq(;,NAL1T1CS AND ASTf_ONAUTICS
AIAA-01-2554
on tile stability of a relaxation scheme. Thus it is suf- involving Ut.
ficient to examine a much simplified form of equation
(6) given by PO = B(A)U_
(r) = {A [(A + B,+½) (A+Bt-½)
+ - I}u, (it)
where.X = a/kt/(/Xx) 2, and v =_u'. Including a in
the A term eliminates all explicit dependence on the f' j = B()q. (123
specific evolution method, and the conclusions of the
following analysis apply to all time evolution meth- For block Gauss-Seidel, preconditioner is defined by
ods that can be written in the form of equation (6). dropping terms involving Ut+1.
Because this class of evolution lnethods includes first,-
order backward Euler, the following analysis also has P_ = B(A)Ut +.A(A)Ul-1
implications for steady-state solution met hods as well.
= {_[(A+BI+½)(A+B,_½)
Fourier Analysis
+ C,_}C,+½]- I} U,
We examine the eigenvalues of the Fourier transform
of equation (7), and the effects that two precondi- + ACt_½(A+Bt-½) Ut-, (13)
tioners have on those eigenvalues. For the spatially
: (s(a> +.a(a)). (14)
discrete discontinuous Galerkin method given in equa-
tion (5) but. specialized to one dimension, oquat.ion (7)
becomes Analysis of Eigenvalues
The eigenvalues of _ and 15-1_. are compuled at.
50 values of 0 that are uniformly distributed between
-- Ul 0 and rr, and for A ranging from 10-s to l0s. If p
denotes the degree of the basis functions, then in one
and Q, =: AU, + [Ct+ ½U,+l + B,_} Ul]. (8)
dimension there are p + l eigenvalues for each value
of 0. Figure 1shows the real component of all eigen-
The subscript 1+½ is a shorthand for the double sub-
values of _ for p equals 1-5 and with A = 1.0 (the
script l,l:Sl, and 0t+½ is assigned as follows: 0t+½ = 1 imaginary part is zero). The solid line is proportional
01-½ = 0, Eliminating Q gives
to p4 and closely matches the growth in the magni-
tude of the eigenvalues. Figures 2(a) and 2(b) show
the eigenvalues when block Jacobi and block (;at_ss-
Seidel preconditioning are used, respectively. Not only
+ (Ct+ ½Ul+l + Bt- ½Ut)] are the eigenv;dues bounded for all p but they are
grouped into two or three regions. Both precondi-
+ Cl_½[AUt_l+
+ (O,+lU, u,: o tioners place one eigenvalue near zero, block .lacobi
has one eigenvalue near -2.0, and all of the remaining
eigenvalues are mapped exactly into -1.0. The varia-
tion in the eigenvalues within each region indicates the
degree to which the eigenvalue varies with 0. A clear
+ {A[(A+Bt+½) (A+B_-½) trend observed in both cases is that the eigenvalues
+ I}
become more tightly grouped as p increases. (;auss-
Seidel is not. a symmetric operator, and consequently,
the imaginary component of its eigenvalnes is not zero.
However, the imaginary component is small relative to
= A(A)U___ +B(A)U_+C(A)U_+_. (9)
the real component, so all following results will exam-
ine either the real component or the absolute value of
The Fourier transform is obtained by substituting
the eigenvalues (denoted as the absolute eigenvahm).
Ut = lT-lei°l where i - v/71, and 0 < 0 < rr, to give
The above results are for A = 1, but similar trends
are observed for large and small A. Figures 3(a) and
(10) 3(b) show the eigenvalues for block Jacobi precondi-
_--wa(a, o)u tioning with A = 100and 0.01, respectively. With
A = 100 the eigenvalue pattern is nearly identical to
The precondilioned residual has the fl)rm t'->Pv(X, v). that of A= 1. With A= 0.01 the eigenvalues are more
The preconditioning operator for block .]acobi is ob- tightly gronped with respect to 0, but the maximmn
tained from equation (9) by retaining only the terms and minimmn eigenvalue are not as close to 0.0 and
3
AMEHI('AN INSITIqYI'E OF AERONAUTICS AND ASTD()NAI!TlCJS
AIAA-01-2554
-7.0, respectively. Similar trends were observed for Gauss-Seidel and block ,lacobi for a fixed ,_. We then
block Gauss-Seidel. construct a curve fit that accounts for the dependence
While the largest absolute eigeuvalue limits tile time of the wk sequence on /_.
step of an explicit method, the smallest absolute e[gen-
Block Gauss-Seldel
value is the least damped and will limit the conver-
For block Gauss-Seidel, we consider a two step se-
gence of an iterative method. Clearly both precondi-
quence whose stability is given by,
tioners reduce the maximum absolute eigenvalue, but
this is of little value if the nfinimum absolute eigen- G_,(_'._,, or(O)) = G(1.0, o'(0))G(w2, o'(0)),
value is equally reduced. Thus, another metric for a
preconditioner is the ratio of the maxinmm absolute and we choose w_ to obtain an optimal damping rate.
eigenvalue t.o the minimum a.bsolnte eigenvalue. Fig- By choosing ,zl = 1.0, the first step completely damps
ures 4(a-c) show this ratio of absolute eigenvalues with all error modes with an eigenvalue of -1.0 and the
no preconditioning, with block .lacobi preconditioning, analysis can focus the single remaining error mode
and with block Gauss-Seidel preconditioning, respec- whose eigenvalue is near 0.0. This slowly decaying
tively. (Note: with no preconditioning the mininmm error mode, denoted a.s oh(0), is dominated by its real
absolute eigenvalue is always 1.0, so the ratio is equal component which varies monotonically with 0. This
to the maximum absolnte eigenvalue, and the data is behavior is typical for all p and )_ examined: p < 5
and 10-s < .X<: l0s.
just the negative of that shown in figure 1; however,
the data is redrawn on a log scale to facilitate compari- In the following we examine two criteria for choosing
son with the other two methods.) Block .lacobi reduces _2. The first, at:,proach minimizes Gg., across the entire
the ratio by about an order of magnitude (at p = 5). spectrum of 0 by requiring
Block Gauss-Seidel further reduces the ratio by about
half an order (again at. p = 5). So both precondition-
ers provide considerable improvement in convergence, Figure 5 shows Gg, for 0"1(0) (solid) and O'l(rC)
but not the three orders of magnitude thal. one might (dashed) as a function of w2 for p = 3, and for
infer simply by looking at, the reduction of maximunl A = 0.01.0. l,l.0 and 10.0. The symbols show
eigenvalues. G_._(_2, rrl(0)) for a uniform sampling of 0 which gives
an indication of the variation of Gu., with 0 at. a fixed
Optimal Relaxation Schemes w_. The symbols also verify that either (;v.,(,J,, _rl(0))
The tight grouping of eigenvalues means entire or Ge.,(w=,,crl(_r)) bounds Gg.,(w:,O) above when
groups of the error modes can be rapidly damped, is not small, and their intersection defines an opti-
sometimes completely damped, by the following simple mal _2 over all 0. For extremely small A, neither
update method G_., (w:,, crl(0)) or Gg.,(w_, _rl(rr)) bounds (;,.j, ;however
equation (17) is still a suitable criterion for minimiz-
U k+l = U k+ _P-I_(A, Uk). (15) ing (;g, over all 0. Although the the damping rate is
stable for all A. the damping rate rises quickly with A,
The stability of this relaxation step is given by and is unacceptablely high at A = 10.0. Also, when
G(w,o(O)) = I1 + we'(0)] where or(0) is any eigen- A is large, (;y_,(w2,cq(rr)) grows rapidly with w,, and
value of the operator P-tT_. By choosing w = 1, quickly becomes unstable. The sensitivity to w_ is an
all the error modes with an eigenvalue of -1.0 are issue because lhe analysis can only be performed for
completely damped. Similarly choosizlg ,,' = 0.5 will the idealized case of a uniform grid. In a realistic case
rapidly damp the error mode with an eigenvalue near in which )_ varies due to non-nniformity in the mesh,
-2.0. This leaves a single error mode associated with the actual optimal value of w2 may be different from
the eigenvalue near zero, which is not readily damped that predicted by the analysis, and the predicted op-
by any choice of w that. is not unstable for the other timal value may result in an unstal)le method.
error modes. The block Gauss-Seidel method can be
A second and more robust approach is to minimize
overrelaxed with _0approaching 2 for large p, but. this
(';u., in the high freqnency portion of the spectrum (0 >
provides only a slight improvement in the damping rr/2) t,o obtain a relaxation scheme that is well suited
rate of the slow error mode.
as a multigrid r smoother. To this end, we require that
However, if one considers the combined effect of a
sequence of relaxation steps (;9s(_2, o-1(¢r/2)) m Ggs(a_2, O-l(rr)). (18)
Figure 6(a) shows the w2 obtained by solving equation
U _+1 =U k+wkp-I"R(),,Uk), k= 1...h (16)
(18) for a range of .X. Also shown is the curve defined
by
each with a different w', then it is possible to choose
a sequence of wk that will rapidly damp all eigen- _,,, + l + (_m - 1) t.anh((_ *Iog_o(A ) + b)
modes. In the following sections, we present several _f _- 2
criteria for selecting an optimal sequence ofw for block (1.¢))
4
AMERICAN INSTITUTE OF AERONAUTICS AN[) A_;TilONAI!TICS
AIAA-01-2554
which provides an accurate fit of"22 for all A. Tile coef- with respect, to changes in 0 and "2:_,respectively. The
ficients a and bare obtained from a least, square fit over damping rate for the high-frequency portion of the
the linear subrange seen in figure 6(b) and w,, is tile spectrum does not increase above _ 0.76 as A be-
nlaximum value of w2, taken here to be the value of"22 conies large. Also, the damping rate increases slowly
when A = 10s. From figure 7 which shows the depen- as a,'a varies from the optimal value indicating that
dence of Gg., on w2 at. A = 10_ (similar to figure 5), we the method is robust in this respect.. This analysis has
can see that large variations in _,, can occur before the been performed for all p < ,5and results are summa-
method becomes unstable, and that small variations rized in Table 2. As seen previously', the asymptotic
do not significantly degrade the damping rate. Hence, convergence rate is nearly constant with respect to p.
this criterion produces relaxation methods lha.t are ro-
Nmnerieal Examples
bust. with regard to any' discrepancies that might occur
In this section we provide a few numerical exam-
between the predicted opt.imal w2 and actual values.
pies that show that. the performance predicted by"the
Figures 8(a) and 8(b) show the behavior of _ (%f)
analysis are realistic and achievable. Details of the
as a function of 0 and A, respectively. The maximum
numerical implementations are not provided here but
value over the high frequency spectrum generally, oc-
are throughly described m the references 3 and 7. The
curs at 0 = rr/2, and rises to _ 0.53 for large A. So
block Jacobi preconditioned discontinuous Galerkin
one would expect this relaxation method to be an ex-
method is applied to the one dimensional heat equa-
cellent smoother for a multigrid method. It should
tion with periodic boundary conditions
be noted that, because the first, relaxation step com-
pletely damps all error modes of this idealized linear Ou/St = V2u
problem except for the 0"l(0) mode, it. is possible to u(0,.r) = sin(rrx) - I < x < l
perform a single relaxation step with "2= 1.0 followed
in which the solution is approximated by polynomi-
by any number of relaxation steps with "2= wI . In this
als of degree p = 5. Tests are performed in which
case the asymptotic convergence rate would reduc, e to
the solution is evolved in time using either backward
(0.53) 2 = 0.28. In the above discussion, p = 3 has
Euler or third order implicit 1Ruuge-Kutta and with
been used to illustrate the procedure. The same anal-
At(cid:0)lAx) e = 10. However, here we are interested only
ysis has been applied to all p _<5 and the results are
the convergence of the temporal residual at. each time
summarized in Table 1. We find that the asymptotic
step, or each [{unge-l*(utta stage, which is similar for
convergence rate is essentially constant as p increases.
both evolution methods. Figure 11 shows the conver-
Block Jacobi gence histories of a typical stage obtained by block
For block Jacobi, we consider a sequence of three Jacobi relaxation with constant "2 = 1.0 and opti-
relaxation steps in which _1 = 0.5, w2 = 1.0 and mal with the three-stage "2 sequence. Also shown is
"2a is chosen t,o provide optinlal damping of the high the convergence of block Jacobi relaxation with opti-
frequency portion of the spectrunL 0 > 7r/2. The sta- mal three stage..2 sequence combined with a two-level
bility of such a method is given by mult.igrid method. The multigrid method is a st.an-
darcl V-cycle with 6 relaxation steps on the line grid
C;bo(_'a,_(0)) = c;(0.5, _(o))(;( l.o, ,_(o))(;("2a,,_(o)) between coarse grid cycles. The coarse grid solution
(20) is converged t.o machine zero by block Jacobi relax-
where or(0) is any eigenvalue of Ih_17_. Also let, oh(0) ation to sinmlate the effect, of multiple grid levels.
and rr2(0) denote the maxinmm and minimum eigen- The case using the optimal a: sequence without multi-
values, respectively. Because the the ininimum eigen- grid shows some improvement in the damping rate
value c_2(0) is near but not exactly -2.0, tim associated over that of constant w. The number of relaxation
error mode is not. completely damped by the relax- steps required to drive the residual to 10-I-" drops
at.ion step with _'t = 0.5 and this eigenmode must from 97(}6 to 1731; however this improvenwnt is not
be considered when determining an optimal wa. In sufficient to make this approach practical. With nmlti-
particular, we choose wa by exalnining the intersec- grid, the conw'rgence improves dramatically and only
tions of the fimctious Goj("2a,0"t(rr)), Go("2a,m2(rr)), 276 relaxation steps are required to drive the resid-
Gbj("23,crt(Tr/2)) and Gbj("23,o'2(Tr/2)). AS shown in ual to 10-le. This corresponds to a convergence rate
figure 9 for p = 3 and A = 1.0, the optimal value of 10-12/'-'r6 = 0.905, and represents a 25 fold reduc-
of "2ais associated with the maxinnml inlersection of tion in work when compared t.o a sixth-order explicit
these curves. For the case shown in figure 9, "23 _ 6.7. t{unge-Nutta method which has an explicit stability
This process is repeated for A ranging from 10-s t,o constraint a of _l/A a,2 = 0.0015.
10s and fitted by equation (19) as done previously for
Extension to Two Dinltenslons
block Gauss-Seidel. Figures 10(a) and 10(b) show the
behavior of The following sections describe the results of the two
dimensional analysis, construction of opt.i,ua[ relax-
(;,,, (,_3,0) - (max ((;bj(_._. c_j(0)), (;b.j(<,, _r2(O))))1/a ation schemes for block aacobi and numerical results.
5
ANIEbtlCAN INS'FITUTE OF AEII()NAIr'PICS AND ASTRONAUTICS
AIAA-01-2554
Two Dimensional Analysis constant except at extremely small A, and w3 is accu-
In two dimensions, equation (7) has the form rately fitted by a tanh distribution to give
_a = 3.72 + 2.72 tauh (1.022 logl0(A ) + 1.469)
•R.(A, Ui.j) = A(A)Ui-I,j + B(A)U,,.i + C(A)U,+I,./
+ 7)(A)Ui,j+I + E(A)Uij_I. (21) Figure 14 shows the asymptotic convergence rate pre-
dicted by the analysis which approaches 0.9a as A
The Fourier transfornl is obtained by substituting becomes large.
Ui,j -_ _Jfi(o_ri+ouJ) to give
Numerieal Examples
Figure 15 shows the convergence history from sev-
_z(A, u],_) = A(_)_ -_°_ + t3(A) + c(:_)¢ i°_
era.l numerical tests. The test case is the two-
+ 9(A)e i°y + £(A)e -i°y
dimensional heat equation with periodic boundary
- _(_,, 0_ 6)0 (22) conditions and with the same initial condition as used
in the one dimensional test described previously. All
The preconditioning operator for block Ja,'obi is sim- results are for the fifth-order (p = 4) discontinuous
ply Pbj = B(A)U/ and 15bj= B(A). Galerkin method with block Jacobi preconditiouing
Typical eigenvalues are shown in figures 12(a) - and with At/_a "_= 10. As seen in one dimension, use
12(c). There are 21 eigenvalues for the fifth-order of the optimal _' sequence without nmltigrid provides
method (p = 4). Figures 12(a) and 12(t)) show the only a small improvement over the constant a: = 1.0
range of each eigenvalue as the wave nnmbers vary case. When relaxation with the optimal w sequence is
over the complete range 0 <_0_:,0u < re. Figure 12(c) combined with multigrid, the solution converges at a
shows a contour plot of the fifth eigenvaIue. The anal- spectral radius of 0.89. This rate is faster than that
ysis predicts that both block Gauss-Seidel and block predicted by' the analysis and represents a 14 fold re-
Jacobi preconditioners bound the eigenvalues, as seen duction in work a.s compared to an explicit method.
in the one dimensional analysis. Several eigenvalues
Concluding Remarks
are mapped into -1.0, however, the remaining eigen-
values are not tightly grouped but are distributed over Block Jacobi and block Gauss-Seidel relaxation
the range 0 to -1.0 for block Gauss-Seidel and 0 to nlet.hods have been analyzed for the solution of the
-2.0 for block Jacobi. discontinnous (;alerkin method applieJ to diffusion.
A Fonrier eigenvalue analysis indicates that both pre-
Optilnal Relaxation
conditioners bound the eigenvahtes. In one dimen-
The eigenvalues are not tightly grouped making it sion. the eigenvalues are grouped into well defined
necessary to design a relaxation scheme that optimally regions, both preconditioners place one eigenvalue
damps a broad spectrum. We choose a sequence of near zero, block Jacobi has one eigenvalue near -2,
three relaxation steps, similar to the previous approach and all of the remaining eigenvalues are mapped ex-
for block Jacobi; however, all three relaxation fa('tors actly into -1.0. In two dimension, several eigenvalues
are allowed to vary freely. The amplification is given are/napped into --1.0; however, the remaining eigen-
by values are distributed within the bounding region. In
general, the eigenvalues become more tightly grouped
(;2d(a,'l,W2,W3, or) = (;(wl, o)G(w2, o')G(w3, t'r) (23) and the minimum eigenvalue moves closer to zero as p
increases. Block Jacobi reduces the ratio of the maxi-
The objective of the optimization, ilh|straled in figure
mum eigenvalue to the mininmm eigenvalue by about
13, is to minimize (;2a for all a in the range -2.0 <
an order of magnitude; block Gauss-Seidel reduces the
cr< or'. or*denotes the largest eigenvahte in the high- ratio by an additional half an order. Optimal relax-
frequency portion of the spectrum and it varies with ation inethods that employ a sequence of under- and
A and it). Note that we are not explicitly concerned
over-relaxation are formulated that provided bounded
with the dependence of a on 0. This optimization is
convergence rates of the high frequencies as the time
accomplished by' equating
step increases towards infinity. Numerical examples
with block Jacobi and multigrid demonstrated conver-
G2d(G') = (;'2d(Ga) = G2d(@b) = G'2d(--2.0) (24) gence rates near the predicted values. When the use of
large time steps is appropriate, or when steady state
where o'_ and or6 denote c_ at the local extrelna of
solutions are desired, the methods described here of_
G2a(cr) and are given by OG2a/Ocr = 0.0. Equation
fer more than an order of magnitude reduction in the
(24) is quadratic, which allows c,_ and c% to be de-
work.
terlnined analytically as flmctions of wl, we and _3-
Equation (24) can now be solved to determine wi, __'2 References
and _,'3. Table 3 gives or*and the optimal ,.osequence
lllarold Atkins, Abdelkader Baggag, Cam Ozturan, and
for a range of A and for p = 4. i.i,'1 and _'e are nearly David Keyes. "Parallelization of an Object-Oriented Unstruc-
6
AMEP, ICAN INSTITUTE OF AERONAIrTI(?S AN[_ ASTItONAIJTICS
AIAA-01-2554
tured Aeroacoustic S_lver," Ninth SIAM Conference oll Parallel
Processing for Scientific Computing, March 22-24. 1999.
2tlarold Atkins, Chi-V(ang Shu, "Quadratur,_-Free Imple-
mentation of Discontinuous (;alerkin Method for Hyperbolic
Equations", AIAA Journal Vol. 36, No. 5, 1998, pp. 775-778.
3Harold Atkins, Chi-Wang Shu, "Analysis ¢)fthe Discon-
tinuous Galerkiu Method Applied to the l)iffusi,m Operator", 500
AIAA Paper 99-3306, 1999.
4j. Tinsley Oden, Ivo Babuska, and Carlos Frik Baumann, 0
"A Discontinuous hp Finite Element Method f,w Diffusion Prob-
-500
lems," Journal of Computational Physics, Vol. 1;t6, 1998, pp.
1-29.
'-...--1000
5Bernardo Cockburn, Chi-Wang Shu, "The l,ocal Discon- 4
P
tinuous Galerkin Method For Time-Dependent Convection-
-1500
Diffusion Systems," SIAM Journal of Numerical Analysis, Vol.
35, 1998, pp. 2440-2463. -2000
6Harold Atkins, David Lockard "A High-Order method us-
ing Unstructured Grids for the Aeroacoustic Analysis of Realis- -2500
.... .... .... +.... , .; ....
tic Aircraft Configurations", AIAA Paper 99-1945, 1999. 0
7Brandt, A., "Multi-level adaptive solutions to boundary- P
value problems." Mathematics -f Computation, \:(,I 31, pp. 333-
390, 1977. Fig. 1 Eigenvalues of residual operator, A= 1.0
Table 1 _,/ curve fit parameters and asymptotic
convergence rate for block Gauss-Seidel
,/c;., )
p _,Mm a for X= l0 s
1 2.32 1.1858 0.314 0427
2 3.85 1.1486 0.623 0.55l
3 6.15 1.1144 0.777 ().5:;_
4 9.13 1.0878 0.!)19 0.549
5 12.79 1.0781 0.917 0.557 (a) 0.0 | ¥ I) )
-0.5
Table 2 _,/ curve fit parameters and asymptotic
Re(o)
convergence rate for block Jacobi -1.0 L' [" 0
-1.5
p . t,
for A = l0s -2.0 g b .... n .
1 2 3 4 5
P
1 2.77 1.1919 0.218 0.684
2 5.26 l. 1379 0.451 0.740
3 8.76 1.1194 0.556 0.762 (b) o,o
4 13.25 1.0699 0.495 0.773 I i ! ) '
-o.5
5 18.75 1.0200 0.602 0.779
Re(o)
-1.0 13 :x V [> o
Table 3 a* and optimal _ for block Jacobi precon-
ditioning with p = 4 -1.5
/_ ('7* _2 _:_ -2.0 ' ' 1n . . 2, . , . .3i . . . 4' .... 5L . .
P
10-2 -0.297 0.5302 0.87057 2.4309 Fig. 2 Eigenvalues of preconditioned residual op=
0.1 -0.0845 0.5343 0.95946 4.6989 erator with A : 1.0 and using a) block Jacobi
1.0 -0.0308 0.5353 0.98483 6.1458 preconditioning and b) block Gauss-Seidel precon=
10.0 -0.0237 0.5354 0.98830 6.4078 ditioning.
l0 s -0.0228 0.5355 0.98873 6.4413
7
AMEIqlCAN IN_.TITUTE OF AEtt()NA1JT[CS AND ASThONAUTICS
AIAA-01-2554
(a) 0.0 !P l, o _= o_(0) k = I0 4:
I i
1.0 o"=O',(n) I1_V_I_AII1*,*
-0.5 • o=cr(O) *_.uo.,;
•"_.IIII_ I • • I
Re(o)
-1.0 V _ o s 0.8 X=l 4 i!,"_:.."...
-1.5 05 :._r!:i; Ii"I"- •
....,.. ;: :
I .... I .... • .,
-2.0 T _ _ a g 0.3 _:0.1 . tl,.:. I "l:
P I,,_, "'"::
(b) 0.0
00 _ 4.... _ a .... _'o'_
b
-0.5 _7 0)2
A
Re(o)
-1.0 r_ A v b Fig. 5 Optimal _'2 over all 0 located at intersection
A of G_,(w2, _(0)) and G_,(_c2, _(r)).
-1.5 _7
b
o
-2.0 , _ .... J , , , , i .... i , , , i , ,
1 2 3 4 5 (a)
p
6
Fig. 3 Eigenvalues of preconditioned residual op- 5
erator with block Jacobi preconditioning and with: 4
a) A = 100.0 and b) A= 0.01. 0)2 3 / ...o. ofipttteimdalco_t0,:
2 9'/
.../
1 -----@ _ _3- @4_'_
10-s 10e 104 10.2 100 10_ 104 10e 10s
X
(_)
(J---- p=l
104
p=2
p=3 (b)
--- --O _ {>_ p=4 0 O 0 0
10s 4 0
o
"-., 102 ...._____ -_7-.... ---%_-.....--.-{'- ..... 2 o o
o
b o
101 0 oo o
o
---{J _____ -2 o
o
10o ,, ,_ .... _, -4
(b)
D _ p=l
10 4 £_- - p: 2 10-s 10 6 10.4 10-_ 100 102 104 10s 10_
W.... p=3 X
--_-_10a I)=4
o- p=5 Fig. 6 a) Ol)timal w tbr all 0 and curve fit, b)
linear region for curve fit seen in tanh-l(2* _I -
_, 10 2
b w,, - 1)/(u;._ - t)).
10_
D -- --C_-- -C} ....
10( ,, ._ .... _ ....
(c) o = o/n/2)
0 t_ I,,=1
104 o = O_(_T) ..-"'_
p=2
.... _ p=3 0.40(( -__ o o : %0) ....,,_r o°
"_-_I03 _,_ p=4
O- p=5 0.30 _ ..." o o
-_ 102 G 0 O _ O O
13 g._ C) O _. S O O
101 0.20 o.,-'"" [ e o 0
10c ..... o 7-,__o;___ ' _ ,--_- O.lO.-." .,_g_ ©0#l O_ O__ ©O_ Oo
I 0 2 3
, I , , _ , !
Fig. 4 Ratio of maximum and minimum eigenval- 0"005 6 7
ues with A= l0 and a) no preconditioning b) block (02
Jacobi preconditioning and c) block Gauss-Soidel
preconditioning. Fig. 7 Optimal _'2 for all (4 _>_r/2 h}cated at inter-
section of G,_(.c2, cr(zr/2)) and (/_0(_'_, (r(rr)).
8
AMEBICAN [NV;'['I'FUTE OF AERONAI!TI(:S AN[) ASTRONAUTICS
AIAA-01-2554
(a) 1.0 (a) 1.0
_--_9_ - c> X=0.1
_z_ _ - _ X= 1.0
6nl _'J _ _ C_ C_--_
(O .'//2 _L'_./_.,.._. _- X=10
. _+) '-_--_,,-- C_ k =100 0.5
-KI---C_-C3_ <3- ---v_ X= 10
c-- _.= 100
o i .... :;
o .... 1' 'o '_
(b) o
(b) 0.6
1.0-
0.8 _/:/
"G ,1,'2
04
0.2
0.6
, , , i , , , , i , , , h
"tl'O8 10-s 10 4 10.2 100 10 2 10 4 10 6 10e 5 10 15
x 0) 3
Fig. 10 Amplification of relaxation with block
Fig• 8 Amplification of relaxation with block
Jacobi precon(litioning using optimal ,_ sequence,
Gauss-Seldel preconditioning using optinml _., se-
Gm= (G_3)1/3: a) G,, Vs 0 for A = 0.1. 1.0, 10 and
quence: a) (G_+)½ VS 0 for A ---- 0. I, 1.0, 10 and 100,
100, b) Gm insensitive to variation in _,_
b) (G'_,)_ at 0 _>_/_.
10 3
0.8 _.. _=o,(u)
N optimal t% ....... o= o2(_;) 10s
ooN\ '_ \\ , ._....__ o= o,(_J2) R
°_° 10z
10.9
£ 0.4
a- (osequence
L3 1041. _ - b-o)= 1.0
0.2 c-(o_eq.,"mg
I i i , i I I ,
10 0 100 200 300
0 " iteration
5 10 15
0.)3 Fig. 11 Convergence of block Jaeobi for heat e(lua-
tion in one dimension: a) constant w = 1.0 without
Fig. 9 Optimal_ located at the maximum inter-
multigrid, b) optimal .z sequence without multi-
section•
grid, c) optimal ,_ sequence with nmltigrid.
f_
AMEHICAN ]NSTITUTE OF ,A,EI_'ONA !TICS AND ASTRONAUrl'ICS
AIAA-01-2554
(b)
1.0
/
t 0.8
G2d
-1.5I_ O,,=Oy=( 0.6
[::_70x= Oy= r
0 < 0,, Oy
0.4
-2.0 "' J.... "
5 10 15 20 5 10 15 20
eigenmode, k eigenmode, k 0.2
(c)
I I I I I I I I
10-s 10-_ 10-4 10-2 100 102 104 10s 108
O_ 2------ ."-.<"," .'. ".. Fig. 14 Predicted asymptotic convergence rate for
block Jacobl in two dimensions using the optimal
sequence with multigrid, p = 4.
1 2 3
0
Fig. 12 Eigenvalues of block Gauss-Seidel and
block Jacobi in two dimensions: a) distribution
of all eigenmodes for block Gauss-Seidel precondi-
tioning, b) distribution of all eigenmodes for block
Jacobi preconditioning, c) Fourier distribution of
k = 5 eigenmode for block Gauss-Seidel.
101 .
_\ a -(_lsequence
t0 3 "_ -- - b-o) seq.;mg
R 10 s \\ --_ d-o)=l.0, mg
10-z
0.8
0.4
10_3 I I I
G2d 1O0 200 300
iteration
0.0
Fig. 15 Convergence of block Jacobi for heat equa-
-0.4 tion in two dimensions: a) constant ,z = 1.0 without
multigrid, b) optirnal a,, sequence without multi-
-0.8 -0'...5.. -;.o.... .... 4.o grid, c) constant ,J = 1.0 with multigrid, d) optimal
0.0 ,v sequence with multlgrid.
X
Fig. 13 Optimal amplification for all _r such that
cr < _r <2.0.
10
AMERICAN INSTITUTE OF AERONAUTICS AND ASTt_ONAI;TICS