- Email: [email protected]
On element-by-element preconditioning for general elliptic problems
Computer Methods in Applied Mechanics and Engineering 92 (1991) 215-229 North-Holland
On element-by-element preconditioning for general elliptic problems H a n - C h o w Lee and A.J. Wathen ~ Department of Computer Science, Stanford University. Stanford, CA 94305. USA Received 21 September 1989
lterative solution procedures for finite element equations become competitive with widely used direct techniques particularly for large three-dimensional problems when computer storage and time are limited. The conjugate gradient method has attractive theoretical properties and is especially efficient when a good preconditioner is used. The efficiency of element-by-element preconditioning has previously been investigated by the second author on model constant coefficient problems on regular grids in two and three dimensions. In this paper we extend that analysis to problems with discontinuous variable coefficients on regular grids in two and three dimensions and to constant coefficient two-dimensional problems on elements with variable aspect ratio.
1. Introduction
For large symmetric systems of finite element equations, there has been a resurgence of interest in iterative rather than direct solution techniques. Limitations of storage and time particularly for large three-dimensional problems make direct methods less attractive, however iterative techniques have yet to be proved comparable in terms of general applicability and reliability. The preconditioned conjugate gradient (PCG) method [1] has been successful in many application areas [2, 3] and various preconditioning strategies have been derived (see, for example, [4-7]) for use with finite element equations. In this paper we examine one of these techniques: the element-by-element preconditioned conjugate gradient method (EBEPCG) of Hughes [7] (see also [8-10]). This is a generally applicable method for symmetric finite element equations, and has been successfully applied to a wide range of practical structural mechanics problems by Ferencz [11]. An analysis of this method for regular discretisations of Poisson problems in two and three dimensions has been presented by Wathen [12]. In this paper we extend that analysis to problems with discontinuous variable coefficients discretised on regular grids, and to constant coefficient problems discretised using elements with large aspect ratios. It is our intention to indicate how these practical considerations affect the convergence of the EBEPCG method and hence to provide general guidelines on the use of EBE methods for more general problems. 'Permanent address: Schooi of Mathematics, University of Bristol, University Walk, Bristol, BS8 1TW, England. The work of this author was partially supported by the National Science Foundation grant No. CCR-8821078. 0045-7825/91/$0:3.50 © 1991 Elsevier Science Publishers B.V. All rights reserved
H.-C. Lee, A.J. Wathen, Element-by-element preconditioning
216
Wathen [12] presents comparisons of EBE preconditioning with simple diagonal preconditioning (scaling). We will make this same relative comparison in our results. In Section 2 we briefly describe the EBE preconditioner, and give an analytic proof of one of the computational results of Wathen [12] on the EBE method for a simple one-dimensional problem. Our computational analysis of EBEPCG on rectangular elements with variable aspect ratios is given in Section 3. In Section 4 we present the results of a computational analysis of EBEPCG on the problem V- (o'Vu) = f, with Dirichlet boundary condition in both two- and three-dimensional regions where the coefficient o" is discontinuous. We give our conclusions in the final Section 5.
2. Element-by-element preconditioned conjugate gradients The EBEPCG method is described in its various forms by Ferencz [11] (see also [8, 9]). We concentrate on the method based on Winget regularisation [7] and Crout factorisation. We precisely follow the same description of the EBE preconditioner given by Wathen [12] which we now briefly indicate here. A global finite element coefficient matrix A say, is assembled from small individual element matrices E e. If we denote by A ~ the assembled contribution to the global matrix from the element e, then A=ZA
(2.1)
~•
We are careful to draw the distinction between the small matrices E c which contains only local information, and the large assembled sparse matrices A c which contains exactly the same real numbers as Ec, but also the positions of the non-zeros in A" additionally, implicitly define the connectivity of the grid. If W= diag(A)
and
W" = diag(A"),
(2.2)
then the 'scaled regularised element array' [7] (also called the Winget regularised array) is the positive definite matrix ft" = I + W - I ' 2 ( A ' -
W c ) W -1~2 .
(2.3)
P
A Crout factorisation ,,i~" -- Le~,'Lee'
(2.4)
of each regularised element array is then easily computed on the element level, and the global EBE preconditioner is defined by
where nel is the overall number of elements.
H.-C.
Lee, A.J.
Wathen, Element-by-element
preconditioning
217
The solution of linear equations P z = r (as is required at each conjugate gradient iteration) is then simply a sequence of diagonal divides and forward and backward substitutions on the triangular factors 11 ~ " and 11 ~e,. In his computations, Ferencz has shown that there are better element reguiarisations and better triangular factorisations that can be used for some problems, however the WingetCrout combination is generally go'~l. In his analytical work, Wathen [12] has compared the use of P (2.5) with the use of W (2.2) alone as preconditioners for A. He shows that P is certainly better, but that the two approaches are spectrally equivalent, i.e., r ( P - ~ A ) / r ( W - ~ A ) tends to a constant value as the global number of degrees of freedom grows (here the condition number K(B) is the ratio of largest to smallest eigenvalue of a matrix B). For larger values of r the conjugate gradient convergent ra'~e is poorer (see, for example, [4]). The asymptotic constant depends on the spatial dimension in the model problem, and must also be expected to depend on the precise problem specification for more general problems• In the simplest case of a model onedimensional problem we have now an analytic proof of this spectral equivalence of EBE preconditioning and diagonal preconditioning, and present that here. We consider a second-order differential equation - u" = f on [0, 1] with Dirichlet boundary conditions• The problem is solved on a regular grid with linear elements. The diagonally scaled stiffness is simply a tridiagonal matrix given by 1
-' 2
= W-I"-AW-!'2
"•
=
"•
- ~2
1
(2.6)
,
_1 2 1!
1
where A is the global coefficient (stiffness) matrix and W= diag(A). By using natural (consecutive) ordering of elements, the Crout EBE preconditioner resulting from Crout factorisation can be written as "1
1 _1
w-l/2pw
-!/2 =
l _!2
2
0
!
4
1
•
°•
! m
-1
--~I 1
0 3 4
•
0
I
-~_
1_
3
m
! .1
_ !2 1
0 -~
~I 0
X
(2.7)
! 4 2
1.
where P is the EBE preconditioner as defined by (2.5)•
H.-C. Lee, A.J. Wathen, Element-by-element preconditioning
218
It has been shown by computation [12] that
K ( W - , , 2 A W -,,2) ~_ 9 r ( p - ~A) for a large number of degrees of freedom• We will show the above result analytically using Rayleigh ~quotient eigenvalue bounds. Since A is a Toeplitz matrix, use of discrete Fourier analysis yields the smallest (A~i.) and largest (Am,.) eigenvalues of W - ' / 2 A W -~/2 as '~m~. = ;tl = 1 -- COS ~,. Am. ~ =
(2.8)
)t, = 1 - cos nO, = 1 + cos 0,,,
(2.9)
with the corresponding eigenvectors Yl = (sin 0,,, sin 20,, . . . . sin nO,,)
and
y,, = (sin 0,,, - s i n 20,, . . . . .
( - I) ''+' sin nO,,),
respectively, where 0,, = 7r/(n + 1). Now, the eigenvalues ?t of A(= W - 1'2AW-!12) satisfy the Rayleigh quotient min y ~ y / y ' y
<. A <~max y t ~ y / y , y ,
y~O
(2.10)
),#1}
yt~y/yty.. . rain y~U
,~ , = ylAyt/ylyl
(2.11)
and t "~.
i
i
t
~
t
max y , ~ y / y y = y,,Ay,/y,,y,, . By multiplication we obtain
V'vAY", = [sin O, , sin 2O',,... , sin nO,]
(2.12)
['
_!
Ir sin O.1 llsin2o.l
2
1
_ 12
.
1
_
_! 2 n
1 JLsin'n0,,J
n - 1
= ~ sin'- kO, - ~ sin kO, sin(k + 1)0,. k=l
(2.13)
k=l
Using trigonometric identities, (2.13) can also be written as n - I
y t~Ay , =
~1
~n(1-cos O',)+ ~ ~ cos(2k + 1)0,- ~ ~ cos2kO,,. k=0
Also, taking inner products, we obtain rl
Y'IY, = ~] sin 2 kO,,. k=l
k=l
(2.14)
219
H.-C. Lee, A.J. Wathen, Element-by-element preconditioning
Again this can be rewritten as n
ytly , = ½ n - ½ ~. cos2k# N .
(2.15)
k=l
Substituting (2.14) and (2.15) into (2.11) and letting n--->Qo it follows that
1 cos
Y'IAYl - 1 1! 3tYl
on -
n
2., , = . o cos(2k + 1)0,,
1 - i ~ cos 2kO. n~=l
- 1 - c o s ( ' t r / ( n + 1 ) ) = ."mi, as n---,oo, i.e.,
(2.16)
N
(1 +
--7---=1YlY, =
O(n-2)
I I( 1
2!
+ .....
2~
...
as n---* = .
Equivalently we have y'v4y, = O ( n - ' )
as n--->oo,
ytly , = O ( n )
as n---,oo.
Similarly, successive multiplication yields
['
y,,Ay. = [sin 0,,. - s i n 20., . . . . ( - 1 ) "÷' sin nON] -.1, lr sin 0,, -½ 1 . ".. /J -sin20,, ×
" "
1 -~[/ -½
n
si"30N
1
/ /
1 .JL(-1) n+'msinn0,,.J
n- I
= ~, sin 2 kON + ~ sin kO,,sin(k + 1)On, k=l
k=l
and by using trigonometric identities the above expression becomes n
n-!
' ~ , = ½ n ( l + c o s On) - cos 0N- ½ ~ cos2k0,, - ½ ~ c o s ( 2 k + l ) 0 y NAY k=O
N.
(2.17)
k=l
The inner product of YN with itself gives n
y:,y,, = ½ n - ~ ~, cos2k0,,. k=O
By substituting (2.17) and (2.18) into (2.12) we obtain
(2.18)
H.-C. Lee, A.J. Wathen, Element-by-element preconditioning
220
n-I
t -
y.Ay. _ 1 +
cos On
2 ~ cos(2k + 1)O,, n cos On ,=l
-
(2.19)
-
t
1 _ _1 ~ cos 2kO,,
Y"Y"
nk=l
As n--> o0 this leads to y,,Ay,/y~y~ -- 1 + cos 0, = 1 + cos('tr/(n + 1))
(2.20)
= '~max •
Hence, when n---~ ~ we have Ami n =
1 - - COS 0 n ~
A ~
1 + COS
8.
(2.21)
= Area x •
Let ~ be any generalised eigenvalue ef the pencil ( . 4 - ~ L D L ' ) , defined by (2.7). /x satisfies the generalised Rayleigh quotient
where L and D are
min <~ ~ ~< max . ~.o x t L D L 'x x~,o x t L D L tx
(2.22)
We choose x = Y t for computing an upper bound on the smallest generalised eigenvalue and x = y,, for a lower bound on the largest generalised eigenvalue, where yj and y , are the eigenvectors corresponding to Am~o and A.,,x, respectively, of A. Thus x 'A x y t ,~ y , rain ~< ~-~,,, x ' L D L t x y'lLDL'yt
(2.23)
and xt~x mar
.~,o x t L D L 'x
~>
t ~ y ,, y ,,A y,,' L D L ty , •
(2.24)
But we have the following equations resulting from multiplications: 3n 3 yt1LDLty I = ~ ( 4 - 3cos 0,,) + ~ c o s 0 , -
3 1( 1 )2 ~-~ sin-" 0,, + ~ sin 0,,- ~ sin 20,,
I
n-
3 15 cos 2kO,, + -~ cos(2k + 1)0,,, 32 *=l k=l t
t
y , , L D L y,, -
15n 32
3 ( n - 1) 8
~., z., cos
3 sin 2 0,, - ~,( sin 0,, + ~, sin 20. )2
cos 0,, - ~
-
3
k=l
n -
]
cos(2k + 1)O,,. k=l
Substituting (2.14) and (2.25) into (2.23) and taking the limit we obtain ytiAy, yt, L D L ' y l
16(1 - cos 0.) 3(4 - 3 cos 0,)
(2.25)
jlLmi n
as
/'z . ~
oQ
(2.26)
H.-C. Lee, A.J. Wathen, Element-by-element preconditioning
221
Likewise putting (2.17) and (2.26) into (2.24) yields 1 + cos
yI, LDLty,,
o,,
~ + ] cos On
~max as n ~ ~ .
Thus when n---* oc we have Am.,,, 15 3 " <'~ + COS O.
and
Amin t> 3 /£min
(4 - 3 cos 0,, ).
If we let
K(W-"2AW -'/2) = om(e-mA) ,
but
K(W-"2AW -''2)
= K(,A) - Am~lx ,tm~"
(2.27)
and r,(p-mA) = x ( ( L D L , ) - t ~
) _ ~m.,, ~l'min
hence by substitution, (2.27) may be written as Am,~/am,,
=
algm..Jt~mi.] •
it follows that ~max
As n---* ~, cos
0,, =
~N'mm
~6 + ~ COS 0.
cos(~/(n + 1))---> 1, so a ~ 9 . This gives
K(W-"2AW-~'2)~9r,(P-~A)
as n---,~.
(2.28)
We note that Wathen [12] computed the value a = 9 in his analysis. Equations (2.28) and (2.21) imply that r,(P-~A)>~O(n 2) as n - - ~ .
3. Distorted elements in order to study the effect of element distortion on the E B E P C G method and the simple diagonally scaled PCG method, we use the two-dimensional elliptic boundary value problem
-V2u=f
on ~O=[O, or]x[O, 1],
!
(3.1)
with Dirichlet boundary condition on the entire boundary of ~2 as our test problem. Regular grids of rectangular bilinear elements with various element aspect ratios (i.e., the ratio of largest and smallest element length) were used for the analysis. Let the distorteti element be o'h x h and 1 < or < ~:, then the element aspect ratio is defined
H.-C. Lee, A.J. Wathen, Element-by-element preconditioning
222
by trh/h, or simply or. The element coefficient (stiffness) matrix for each element e is then given by ( & =34 k' ~r+
=~ ( ( ~ c r - ) - l ) / 4 ( o r ' - + l )
~'
(-tr2+
½)/4(~r'-+l)
-~l,
(-,r-" + '.)/4(or: + 1)
L
.
(3.2)
where the coefficient k" is constant for all elements e in the case of the constant coefficient problem we analyse here. For a very large aspect ratio (i.e., tr > > 1) the element matrix E e can be expressed asymptotically as
E,,
4k"tr 3
~
.
.
(3.3)
We will later demonstrate from our numerical results that for a very large aspect ratio all performance data do approach to some limit asymptotically as characterised by (3.3). Computations to find the generalised eigenvalues of the symmetric matrix pencil (A - AP) were made and we obtained the largest and smallest eigenvalues of P-~A, where A is the global coefficient matrix for the given problem described by (3.1) and P is the EBE preconditioner defined by (2.5). In all of the computations, the LANSO (Lanczos algorithm with selective orthogonalisation) routines (version 1.0) of Parlett and Liu were used [13]. These results compared with diagonal scaling (W-t~2AW -~'~') for the same problem are tabulated in Tables 1, 2 and 3 for various numbers of degrees of freedom as well as several different aspect ratios. For a fixed number of degrees of freedom or equivalently a fixed mesh size h, the ratio between KOtAC.=- K(W-~I2AW-t'") and Ko~,xc - K(P- ~A) approaches to a limit asymptotically as the element aspect ratio t r ~ 0 . For example, when nnod (number of degrees of freedom) = 102 /(DIAG/KEBE-->8.2 as tr--->oc Table I Numerical results for distorted elements in 2D (nnod" = 10 x 10) Element aspect ratio
Element-by-element preconditioning
Diagonal precondi:ioning KDtAO
(1"
'~min
~max
KEBE
/~min
~max
KDIA(I
KEBE
1 10 l(I-" 10 ~ 104
0.2372 0.1585 0.7106 x 10 i 0.7016 x 10-' 0.7015 x 10- ~
1.0507 1.1705 1.1775 1.1776 1.1776
4.43 7.38 16.57 16.78 16.79
0.5990 x 10 - t 0.4960 x 10- t 0.2136 x 10- ~ 0.21077 x 10 -~ 0.21074 x 10- t
1.4603 2.871 2.899 2.900 2.900
24.36 57.91 135.72 137.57 137.59
5.50 7.85 8.19 8.20 8.20
" N u m b e r of global unknowns (degrees of freedom).
H.-C. Lee, A.J. Wathen, Element-by-elementpreconditioning
223
Table 2 Numerical results for distorted elements in 2D (nnod" = 20 x 20) Element aspect ratio
Element-by-element preconditioning
Diagonal preconditioning KD,^,~
Am,n
Or
1 10 I0-" 10 ~
0.7440 0.8689 0.2043 0.1943
x x x x 0.1940 x
10 4
10- ' 10- ' 10' 10-'
!(1'
J[ma~
KI'BF
1.0507 1.1711 i.1797 1.1798 I. 1798
13.20 13.4~ 57.74 60.71 6(|.81
Am,. 0.1670 0.1670 0.5958 0.5663
x x x x 0.5660 x
10- ' 10 ' 10--" 10 " 10 "
Amax
KI)IA(;
KF.BE
i .4888 2.9427 2.9718 2.9720 2.9720
89.20 176.20 498.77 524.81 525 09
6.28 13.07 8.64 8.64 8.64
" N u m b e r of global unknowns (degrees of freedom).
For nnod = 202 KDIAG/KEBE"~8.64
as o r - - > ~ .
And for nnod = 40-" we have KDIAG/KEBE--'>8.7
as or--~oc.
In the limiting case as or---, ~ the 'elongated' two-dimensional element is degenerated into a one-dimensional type element. Therefore it is expected that the ratio KD,AG/KEBr:should tend to a value of 9 asymptotically as nnod--> ~, which has been proved analytically in Section 2 for the one-dimensional case. As a matter of fact, when nnod = 80", the ratio KDIAG/KEn E is approximately 9. It is interesting to note that the convergence rate for the EBE preconditioning increases at a slightly faster rate than diagonal preconditioning for element aspect ratio smaller than 20. In particular, the EBE preeonditioner becomes especially effective for two-dimensional problems for the element aspect ratio ranging from 10 to 20. We found that KDIAG/KEBE "~- 14.06 when or= 10 and KD,A~/KE, E = 14.12 when or = 20 and deduce that there is an optimal EBE convergence rate for some or, 10 <~ or <~ 20. The relative advantage of the EBEPCG method over the diagonally scaled PCG method diminishes as the element aspect ratio increases. The performance trend is shown in Fig. 1. Table 3 Numerical results for distorted elements in 2D (nnod ~ = 40 x 40) Element aspect ratio or
1 10 10-" 10 ~ 10 ~
Element-by-element preconditioning
Diagonal preconditioning KV~AG
Amin
0.2034 0.2456 0.6208 0.5116 0.5170
x x x x x
10- ' 10- ' 10-: 1O-: 10 -2
Amax
KI=BE
1.0593 !.1764 1.1855 1.1856 i.1856
52.08 47.90 190.96 231.74 229.30
N u m b e r of global unknowns (degrees of freedom).
'~m,n 0.4397 0.4398 0.1799 0.1481 0.15.00
x x x x x
10- 2 10 : l0 : l0-" i0 "
Amaa
KDIA(;
KEBE
1.495 2.961 2.991 2.991 2.991
340.0 673.3 1662.4 2019.6 1994.0
6.53 14.06 8.7 8.8 8.7
224
H.-C. Lee, A.J. Wathen, Element-by-element preconditioning 15
."Ratio of condition t~t~rnbers (DIAG/EBE)
aS 0 ...J.~.nnod~¢o
~9
.
r
n_
.
.
.
.
.
o
B o
~ 7
Element aspect rotto IOgtO 0
Fig. I. Ratio of condition numbers versus element aspect ratio (nnod = 40 ×
41)).
The numerical results indicate that the EBE preconditioning using distorted elements is also spectrally equivalent to diagonal scaling. However, the spectral equivalence bound depends on the element aspect ratio as we have observed in Fig. 1. The EBEPCG method seems to improve slightly over the diagonally scaled PCG method for aspect ratios of or = 1 0 - 20. The ratio Kt~,AG/KEBE is approximately 14.0 for that range as compared with spectral equivalence bound of 6.6 for square elements (i.e., o- - 1). This simply implies that the convergence rate for the EBEPCG method does not get poorer for distorted elements or = 10 ~ 20 as compared to o" = 1. This is in contrast to diagonal scaling which does degrade as or gets larger.
4. Variable coefficients In this section we consider the variable coefficient two-dimensional elliptic boundary value problem -V. (orVu) = f
on O = [0, 1] x [0, 11,
(4.1)
with Dirichlet boundary conditions. Here or(x, y) represents the variable diffusion coefficient. To study in particular the effect of coefficient discontinuities on the EBEPCG method as well as the diagonally scaled PCG method, we use a piecewise constant function tr(x, y) for our investigation. Precisely we take or(x,
Y)
jl, La ,
0~
0~
where the constant o~ may lie in a w!de range of 1 < a < ~.
(4.2)
H.-C. Lee, A.J. Wathen, Element-by.element preconditioning
225
Table 4 Numerical results for variable coefficient problem in 2D # of global
Element-by-element preconditioning
Diagonal preconditioning
unknowns nnod
a
Am,"
Am~x
KERE
21 x 21 (441)
1 10 102 l(I~
0.6813 0.1184 0.1710 0.1306
x x x x
10-' 10- " 10 -4 !0-*
i.057 !.095 !.127 1.138
15.50 9.25 x 10" 659 x 104 8.71 x 10"
41 x 4 1 (1681)
l 10 10"~ 10 ~
0.1940 0.7821 0.9450 0.9495
x x x ×
10 -~ 10 -3 10 -s 10 -7
1.0586 1.1283 1.1457 1.1469
54.58 14.43 x 102 12.12 x 104 12.08 x i0"
'~mm 0.15216 0.15216 0.15217 0.15217
x x x x
Amax Ill- J I0- ' 10-~ 10-i
0.419 x 10--" 0.419 x l0 -2 0.3.045 x 1 0 " I}.4050 x 10-"
KI)IA(;
KO,A. KI~HP-
1.4897 1.5298 1.5572 1.5688
97.90 1(11).54 102.33 103.09
6.32 10.87 x 1 0 -" 15.53 x 10 -4 11.83 x I0-*
!.4958 1.5610 1:5780 1.589()
356.99 372.55 39(}.11 392.34
6.54 25.82 x l0 2 32.18 x 10 -4 32.48 x 10-*
Three different values of a. namely a = 10, 100, 1000 were chosen for the computations, in each case the diffusion coefficient has been characterised by an increasing discontinuity along the line y = ~. In particular when a = 1, it simply reduces to the problem of constant coefficient which was analysed by Wathen [12]. Regular grids of square bilinear elements were used. Table 4 shows the condition numbers obtained for the case of EBE preconditioning and diagonal preconditioning by computing the largest (Amax) and smallest (Am~.) eigenvalues of P-*A and W-'/2AW -I/2. The ratio between KEBE~ K(p-IA) a n d XO~AG- r(W-~/2AW -''2) is also tabulated in Table 4 for comparison with various values of a. The results indicate that for a fixed number of degrees of freedom (i.e., nnod = constant) or equivalently for a fixed mesh size (i.e., h =constant), the condition number KEnE depends strongly on the diffusion coefficient o'(x, y). A simple asymptotic calculation yields KE, E = O(a2). On the contrary, the condition number for the case of using the diagonal preconditioner is independent of the diffusion coefficient. This was also noted by Greenbaum et al. [14]. However, the condition number of the coefficient (stiffness) matrix without any diagonal scaling also depends on the diffusion coefficient, but less strongly than EBEPCG. We have found that x~ -= r ( A ) = O ( a ) from our study as indicated in Table 5. Unlike the spectral equivalences of EBE preconditioning and diagonal scaling which were found for the one-, two- and three-dimensional constant coefficient problems [12], the degradation in performance of the EBEPCG method seemed to increase dramatically as the discontinuous factor a grows. We also obtained estimates of KE, E and KDIAGasymptotically as Table 5 Condition number of A without diagonal scalingin 2D (nnod = 5 x 5) a
Am,"
1 10 102 10 ~ 104
0.192 0.410 0.439 0.442 0.442
~' K, -= , c ( A ) .
Am~,x 1.375 1.310 1.314 1.314 1.314
x x x x
K/' 10 i02 103 104
7.16 3.20 2.99 2.98 2.97
x x × x
l0 l02 i03 10~
226
H.-C.
Lee, A.J.
Wathen,
Element-by-element
preconditioning
nnod ~ oo for a fixed value of a in order to further stt~dy the behavior of discontinuities in the diffusion coefficient. We have as nnod----~oo KEnE = O (nnodP't~)), KOtA~; = O (nodl'zt"),
p2(a) = 1 (constant),
KDIAC/KEBE = O (nnodP~')), where the powers p~, P2 and P3 are plotted against log,, a in Fig. 2. It appears evident that the E B E preconditioning is not spectral equivalent to diagonal scaling as we expected in this case. We want to emphasize that the condition number KEBE grows very rapidly with increasing ot despite the fact of its more mild dependence on the mesh size. For instance, we have /(EBE = O ( h - I )
for a = 1000,
KEBE = O(h -2)
for a = 1 (i.e., with constant coefficient).
Nevertheless, KEnE = 8.714 x 10 6 for a = 1060, n n o d = 4 4 1 , but /(EBE=312.31 for c~ = 1, nnod = 10000. We also present numerical results for the same Dirichlet problem in the three-dimensional problem using uniform trilinear brick elements on ~ = [0, 1] x [0, 1] x [0, 1]. We note that although there are some similarities between this and the 2-D problem just described, they are physically entirely different problems. The numerical results, given in Table 6, shc,w that there is again rapid growth of KEBE with a. For a fixed number of degrees of freedom the condition
t . .......
L .......
t¢. . . . . . . . . . . . . . . .
,It
0.8 07, 06 0.5
~ o.3 ~" 02 0.I 0.-
o
o.'s
I
t
I
,.5
!
2
21s
~Oe~oa Fig. 2. Order of condition number versus 'discontinuous" constant for EBE and diagonal preconditioning.
H.-C. Lee, A.J. Wathen.
Element-by-elementpreconditioning
227
Table 6 3D results for variable coefficient (nnod = 5 x 5 x 5) Element-by-element preconditioning Am,n
'~max
0.6678 0.6373 X 10 .4 0.6195 x 10 -x 0.6167 x 10- " 0.6091 x ! 0 i~
1.0199 1.1408 1.1408 1.1897 1.12288
a
! 10 10-' 10 x !0 ~
Diagonal preconditioning
KEBE 1.527 1.790 !.841 1.929 2.017
x x x x
104 10~ 1012 1(| '~
~m,n
~max
KI)IA/;
0.192 0.192 0.192 0.192 0.189
1.375 !.449 1.472 !.389 1.377
7.16 7.55 7.67 7.23 7.29
Table 7 Condition number of A without diagonal .scaling in 3D (nnod = 5 x 5 x 5) Atom 1 10 102 l0 ~ 10 4 " KI -
0.2750 0.4822 0.5004 0.5021 0.5022
KI"
Amax 1.3499 1.2884 1.2878 1.2878 1.2878
x × x x
10 102 103 104
4.910 2.672 2.576 2.565 2.564
x 10 x l0 a x 10' X 10 a
K(,4).
number in the case of EBE preconditioning applied to the three-dimensional problem is KEnE = O(a4), which is poor even as compared with that of the two-dimensional problem (fEB E = O ( a ~) in 2-D). Spectral equivalence with diagonal scaling which exists for the constant coefficient problem is again lost. However, similar results to the 2-D problem were obtained for both %.AG -- x(W-~t2AW -~'2) and Kt - K(A). In three-dimensions, we have also that KD~AG= O(1) and Kt = O(a) for a fixed number of degrees of freedom. The numerical results are given in Tables 6 and 7. Interestingly, the convergence rate of diagonally scaled PCG is independent of the diffusion coefficient or(x, y, z) as in 2-D. The measure of coefficient discontinuities for both two- and three-dimensional problems was characterised by a single parameter ct in this analysis. For large values of a, e.g., a = 103, and even small values of nnod (120<~nnod<~ 125, say) the condition number for using the EBE preconditioner is O(10 ~') for 2-D and O(1012) for 3-D. The EBEPCG method thus appears impractical for such problems. Three-dimensional results for larger values of nnod were computed, but are not included since they were little different from the numerical results given here. 5. Conclusion We have extended the analysis of EBE preconditioning of the second author [12] to linear elliptic boundary value problems in 2 and 3 dimensions with discontinuous coefficients on regular grids and 2-dimensional constant coefficient problems on grids with variable element
228
H.-C. Lee, A.J. Wathen, Element-by-element preconditioning
aspect ratio. In the latter case we have used bilinear rectangular elements and have shown spectral equivalence of EBE preconditioning and diagonal scaling independently of the aspect ratio. The EBE preconditioner is always better and achieves a much superior performance for aspect ratios between about 5 and 20. For larger aspect ratios the EBE preconditioner is better relative to diagonal scaling than for very small aspect ratios. The case of square elements analysed by Wathen is the worst case. In an absolute sense the EBE preconditioner actually gets better as the aspect ratio of the elements is increased from 1 to approximately 20, and then deteriorates for large values. The condition of the underlying stiffness matrix, of course, deteriorates also as the elements are elongated. For problems with a discontinuous coefficient, we find that the EBE preconditioner deteriorates very quickly as the size of the discontinuiti¢.s is increased. In 3 dimensions, the deterioration is significantly worse than in 2 dimensions. Diagonal scaling appears to be a better choice in such problems. Greenbaum et al. [14] have found similar poor performance on similar problems with the hierarchical basis conjugate gradient method when diagonal scaling is not explicitly used. We note also that Ferencz [11] reports good results using EBEPCG on a problem of thermal stress in a solid state slab laser which involves large material (coefficient) discontinuities. It seems likely that there is very little coupling between the solution in the different material regions in this particular problem and hence the effect of the discontinuities that we have shown would not arise. The test problem we have used is more like that of a beam composed of two equal halves of very different materials. We conclude by emphasising the attractiveness of EBEPCG from the practical point of view. It is highly vectorisable [10] and parallelisable [15] and is derived rationally from the finite element data structure. For problems wh,~re coefficient discontinuities do not play a major role, it appears to be a competitive method, particularly on advanced hardware.
Acknowledgment All computations were carried out in the Numerical Mathematics and Scientific Computation group in the Computer Science Department at Stanford. We are grateful to Gene Golub for use of these facilities and also wish to thank him for his interest and comments which have helped to guide us through this study. We wish also to acknowledge useful conversations with Bob Ferencz and Beresford Parlett for making the LANSO subroutines available to us.
References [i] E Concus, G.H. Golub and D.E O'Leary, A generalized conjugate gradient method for the numerical solution of elliptic partial differential equations, in: J.R. Bunch and D.J. Rose, eds., Sparse Matrix Computations (Academic Press, London, 1973) 163-174. [2] T.F. Russel, M.F. Wheeler and C. Chiang, Large-scale simulation of miscible displacement by mixed and characteristic finite element method, in: W.E. Fitzgibbon, ed., Mathematical and Computational Methods in Seismic Exploration and Reservoir Modeling (SIAM, Philadelphia, PA, 1986) 85-107. [3] T.E. Tezduyar and J. Liou, Element-by-element and implicit-explicit finite element formulations for computational fluid dynamics, in: R. GIowinski, G.H. Golub, G.A. Meurant and J. P6riaux, eds., First International
H.-C. Lee, A.J. Wathen, Element-by-elemem preconditioning
229
Symposium on Domain Decomposition Methods for Partial Differential Equations (SIAM, Philadelphia, PA, 1988) 281-300. [4] O. Axelsson and V.A. Barker, Finite Element Solution of Boundary Value Problems. Theory and Computation (Academic Press, Orlando, 1984). [5] H. Yserentant, On the multi-level splitting of finite element spaces, Numer. Math. 49 (1986) 379-412. [6] J.H. Bramble, J.E. Pasciak and J. Xu, Parallel multilevel preconditioners, Math. Comp. (to appear). [7] T.J.R. Hughes, The Finite Element Method: Linear Static and Dynamic Finite Element Analysis (PrenticeHall, Englewood Cliffs, NJ, 1987). [8] T.J.R. Hughes, I. Levit and J. Winget, An element-by-element solution algorithm for problems of structural and solid mechanics, Comput. Methods Appl. Mech. Engrg. 36 (1983) 241-254. [9] T.J.R. Hughes, I. Levit and J. Winger, Element-by-element implicit algorithms for heat conduction, j. Engrg. Mech. 109 (2) (1983) 576-585. [10] T.J.R. Hughes, R.M. Ferencz and J.O. Hailquist, Large-scale vectorized implicit calculations in solid mechanics on a Cray X-MP/48 utilizing EBE preconditioned conjugate gradients, Comput. Methods Appl. Mech. Engrg. 61 (1987)215-248. [11] R.M. Ferencz, Element-by-element preconditioning techniques for large-scale, vectorized finite element analysis in nonlinear solid and structural mechanics, Ph.D. Thesis, Stanford University, 1989. [12] A.J. Wathen, An analysis of some element-by-element techniques, Comput. Methods Appl. Mech. Engrg. 74 (1989) 271-287. [13] B.N. Parlett and Z.A. Liu, LANSO vl.0 (1989) subroutines, University of California, Berkeley. [14] A. Greenbaum, C. Li and H.Z. Chao, Parallellizing preconditioned conjugate gradient algorithms, Comput. Phys. Commun. 53 (1989) 295-309. [15] B. Gilvary, 1. Gladwell, I.M. Smith and S.W. Wong, The element-by-element method for linear systems, University of Manchester, Numerical Analysis Report No. 144, 1987.
On element-by-element preconditioning for general elliptic problems H a n - C h o w Lee and A.J. Wathen ~ Department of Computer Science, Stanford University. Stanford, CA 94305. USA Received 21 September 1989
lterative solution procedures for finite element equations become competitive with widely used direct techniques particularly for large three-dimensional problems when computer storage and time are limited. The conjugate gradient method has attractive theoretical properties and is especially efficient when a good preconditioner is used. The efficiency of element-by-element preconditioning has previously been investigated by the second author on model constant coefficient problems on regular grids in two and three dimensions. In this paper we extend that analysis to problems with discontinuous variable coefficients on regular grids in two and three dimensions and to constant coefficient two-dimensional problems on elements with variable aspect ratio.
1. Introduction
For large symmetric systems of finite element equations, there has been a resurgence of interest in iterative rather than direct solution techniques. Limitations of storage and time particularly for large three-dimensional problems make direct methods less attractive, however iterative techniques have yet to be proved comparable in terms of general applicability and reliability. The preconditioned conjugate gradient (PCG) method [1] has been successful in many application areas [2, 3] and various preconditioning strategies have been derived (see, for example, [4-7]) for use with finite element equations. In this paper we examine one of these techniques: the element-by-element preconditioned conjugate gradient method (EBEPCG) of Hughes [7] (see also [8-10]). This is a generally applicable method for symmetric finite element equations, and has been successfully applied to a wide range of practical structural mechanics problems by Ferencz [11]. An analysis of this method for regular discretisations of Poisson problems in two and three dimensions has been presented by Wathen [12]. In this paper we extend that analysis to problems with discontinuous variable coefficients discretised on regular grids, and to constant coefficient problems discretised using elements with large aspect ratios. It is our intention to indicate how these practical considerations affect the convergence of the EBEPCG method and hence to provide general guidelines on the use of EBE methods for more general problems. 'Permanent address: Schooi of Mathematics, University of Bristol, University Walk, Bristol, BS8 1TW, England. The work of this author was partially supported by the National Science Foundation grant No. CCR-8821078. 0045-7825/91/$0:3.50 © 1991 Elsevier Science Publishers B.V. All rights reserved
H.-C. Lee, A.J. Wathen, Element-by-element preconditioning
216
Wathen [12] presents comparisons of EBE preconditioning with simple diagonal preconditioning (scaling). We will make this same relative comparison in our results. In Section 2 we briefly describe the EBE preconditioner, and give an analytic proof of one of the computational results of Wathen [12] on the EBE method for a simple one-dimensional problem. Our computational analysis of EBEPCG on rectangular elements with variable aspect ratios is given in Section 3. In Section 4 we present the results of a computational analysis of EBEPCG on the problem V- (o'Vu) = f, with Dirichlet boundary condition in both two- and three-dimensional regions where the coefficient o" is discontinuous. We give our conclusions in the final Section 5.
2. Element-by-element preconditioned conjugate gradients The EBEPCG method is described in its various forms by Ferencz [11] (see also [8, 9]). We concentrate on the method based on Winget regularisation [7] and Crout factorisation. We precisely follow the same description of the EBE preconditioner given by Wathen [12] which we now briefly indicate here. A global finite element coefficient matrix A say, is assembled from small individual element matrices E e. If we denote by A ~ the assembled contribution to the global matrix from the element e, then A=ZA
(2.1)
~•
We are careful to draw the distinction between the small matrices E c which contains only local information, and the large assembled sparse matrices A c which contains exactly the same real numbers as Ec, but also the positions of the non-zeros in A" additionally, implicitly define the connectivity of the grid. If W= diag(A)
and
W" = diag(A"),
(2.2)
then the 'scaled regularised element array' [7] (also called the Winget regularised array) is the positive definite matrix ft" = I + W - I ' 2 ( A ' -
W c ) W -1~2 .
(2.3)
P
A Crout factorisation ,,i~" -- Le~,'Lee'
(2.4)
of each regularised element array is then easily computed on the element level, and the global EBE preconditioner is defined by
where nel is the overall number of elements.
H.-C.
Lee, A.J.
Wathen, Element-by-element
preconditioning
217
The solution of linear equations P z = r (as is required at each conjugate gradient iteration) is then simply a sequence of diagonal divides and forward and backward substitutions on the triangular factors 11 ~ " and 11 ~e,. In his computations, Ferencz has shown that there are better element reguiarisations and better triangular factorisations that can be used for some problems, however the WingetCrout combination is generally go'~l. In his analytical work, Wathen [12] has compared the use of P (2.5) with the use of W (2.2) alone as preconditioners for A. He shows that P is certainly better, but that the two approaches are spectrally equivalent, i.e., r ( P - ~ A ) / r ( W - ~ A ) tends to a constant value as the global number of degrees of freedom grows (here the condition number K(B) is the ratio of largest to smallest eigenvalue of a matrix B). For larger values of r the conjugate gradient convergent ra'~e is poorer (see, for example, [4]). The asymptotic constant depends on the spatial dimension in the model problem, and must also be expected to depend on the precise problem specification for more general problems• In the simplest case of a model onedimensional problem we have now an analytic proof of this spectral equivalence of EBE preconditioning and diagonal preconditioning, and present that here. We consider a second-order differential equation - u" = f on [0, 1] with Dirichlet boundary conditions• The problem is solved on a regular grid with linear elements. The diagonally scaled stiffness is simply a tridiagonal matrix given by 1
-' 2
= W-I"-AW-!'2
"•
=
"•
- ~2
1
(2.6)
,
_1 2 1!
1
where A is the global coefficient (stiffness) matrix and W= diag(A). By using natural (consecutive) ordering of elements, the Crout EBE preconditioner resulting from Crout factorisation can be written as "1
1 _1
w-l/2pw
-!/2 =
l _!2
2
0
!
4
1
•
°•
! m
-1
--~I 1
0 3 4
•
0
I
-~_
1_
3
m
! .1
_ !2 1
0 -~
~I 0
X
(2.7)
! 4 2
1.
where P is the EBE preconditioner as defined by (2.5)•
H.-C. Lee, A.J. Wathen, Element-by-element preconditioning
218
It has been shown by computation [12] that
K ( W - , , 2 A W -,,2) ~_ 9 r ( p - ~A) for a large number of degrees of freedom• We will show the above result analytically using Rayleigh ~quotient eigenvalue bounds. Since A is a Toeplitz matrix, use of discrete Fourier analysis yields the smallest (A~i.) and largest (Am,.) eigenvalues of W - ' / 2 A W -~/2 as '~m~. = ;tl = 1 -- COS ~,. Am. ~ =
(2.8)
)t, = 1 - cos nO, = 1 + cos 0,,,
(2.9)
with the corresponding eigenvectors Yl = (sin 0,,, sin 20,, . . . . sin nO,,)
and
y,, = (sin 0,,, - s i n 20,, . . . . .
( - I) ''+' sin nO,,),
respectively, where 0,, = 7r/(n + 1). Now, the eigenvalues ?t of A(= W - 1'2AW-!12) satisfy the Rayleigh quotient min y ~ y / y ' y
<. A <~max y t ~ y / y , y ,
y~O
(2.10)
),#1}
yt~y/yty.. . rain y~U
,~ , = ylAyt/ylyl
(2.11)
and t "~.
i
i
t
~
t
max y , ~ y / y y = y,,Ay,/y,,y,, . By multiplication we obtain
V'vAY", = [sin O, , sin 2O',,... , sin nO,]
(2.12)
['
_!
Ir sin O.1 llsin2o.l
2
1
_ 12
.
1
_
_! 2 n
1 JLsin'n0,,J
n - 1
= ~ sin'- kO, - ~ sin kO, sin(k + 1)0,. k=l
(2.13)
k=l
Using trigonometric identities, (2.13) can also be written as n - I
y t~Ay , =
~1
~n(1-cos O',)+ ~ ~ cos(2k + 1)0,- ~ ~ cos2kO,,. k=0
Also, taking inner products, we obtain rl
Y'IY, = ~] sin 2 kO,,. k=l
k=l
(2.14)
219
H.-C. Lee, A.J. Wathen, Element-by-element preconditioning
Again this can be rewritten as n
ytly , = ½ n - ½ ~. cos2k# N .
(2.15)
k=l
Substituting (2.14) and (2.15) into (2.11) and letting n--->Qo it follows that
1 cos
Y'IAYl - 1 1! 3tYl
on -
n
2., , = . o cos(2k + 1)0,,
1 - i ~ cos 2kO. n~=l
- 1 - c o s ( ' t r / ( n + 1 ) ) = ."mi, as n---,oo, i.e.,
(2.16)
N
(1 +
--7---=1YlY, =
O(n-2)
I I( 1
2!
+ .....
2~
...
as n---* = .
Equivalently we have y'v4y, = O ( n - ' )
as n--->oo,
ytly , = O ( n )
as n---,oo.
Similarly, successive multiplication yields
['
y,,Ay. = [sin 0,,. - s i n 20., . . . . ( - 1 ) "÷' sin nON] -.1, lr sin 0,, -½ 1 . ".. /J -sin20,, ×
" "
1 -~[/ -½
n
si"30N
1
/ /
1 .JL(-1) n+'msinn0,,.J
n- I
= ~, sin 2 kON + ~ sin kO,,sin(k + 1)On, k=l
k=l
and by using trigonometric identities the above expression becomes n
n-!
' ~ , = ½ n ( l + c o s On) - cos 0N- ½ ~ cos2k0,, - ½ ~ c o s ( 2 k + l ) 0 y NAY k=O
N.
(2.17)
k=l
The inner product of YN with itself gives n
y:,y,, = ½ n - ~ ~, cos2k0,,. k=O
By substituting (2.17) and (2.18) into (2.12) we obtain
(2.18)
H.-C. Lee, A.J. Wathen, Element-by-element preconditioning
220
n-I
t -
y.Ay. _ 1 +
cos On
2 ~ cos(2k + 1)O,, n cos On ,=l
-
(2.19)
-
t
1 _ _1 ~ cos 2kO,,
Y"Y"
nk=l
As n--> o0 this leads to y,,Ay,/y~y~ -- 1 + cos 0, = 1 + cos('tr/(n + 1))
(2.20)
= '~max •
Hence, when n---~ ~ we have Ami n =
1 - - COS 0 n ~
A ~
1 + COS
8.
(2.21)
= Area x •
Let ~ be any generalised eigenvalue ef the pencil ( . 4 - ~ L D L ' ) , defined by (2.7). /x satisfies the generalised Rayleigh quotient
where L and D are
min <~ ~ ~< max . ~.o x t L D L 'x x~,o x t L D L tx
(2.22)
We choose x = Y t for computing an upper bound on the smallest generalised eigenvalue and x = y,, for a lower bound on the largest generalised eigenvalue, where yj and y , are the eigenvectors corresponding to Am~o and A.,,x, respectively, of A. Thus x 'A x y t ,~ y , rain ~< ~-~,,, x ' L D L t x y'lLDL'yt
(2.23)
and xt~x mar
.~,o x t L D L 'x
~>
t ~ y ,, y ,,A y,,' L D L ty , •
(2.24)
But we have the following equations resulting from multiplications: 3n 3 yt1LDLty I = ~ ( 4 - 3cos 0,,) + ~ c o s 0 , -
3 1( 1 )2 ~-~ sin-" 0,, + ~ sin 0,,- ~ sin 20,,
I
n-
3 15 cos 2kO,, + -~ cos(2k + 1)0,,, 32 *=l k=l t
t
y , , L D L y,, -
15n 32
3 ( n - 1) 8
~., z., cos
3 sin 2 0,, - ~,( sin 0,, + ~, sin 20. )2
cos 0,, - ~
-
3
k=l
n -
]
cos(2k + 1)O,,. k=l
Substituting (2.14) and (2.25) into (2.23) and taking the limit we obtain ytiAy, yt, L D L ' y l
16(1 - cos 0.) 3(4 - 3 cos 0,)
(2.25)
jlLmi n
as
/'z . ~
oQ
(2.26)
H.-C. Lee, A.J. Wathen, Element-by-element preconditioning
221
Likewise putting (2.17) and (2.26) into (2.24) yields 1 + cos
yI, LDLty,,
o,,
~ + ] cos On
~max as n ~ ~ .
Thus when n---* oc we have Am.,,, 15 3 " <'~ + COS O.
and
Amin t> 3 /£min
(4 - 3 cos 0,, ).
If we let
K(W-"2AW -'/2) = om(e-mA) ,
but
K(W-"2AW -''2)
= K(,A) - Am~lx ,tm~"
(2.27)
and r,(p-mA) = x ( ( L D L , ) - t ~
) _ ~m.,, ~l'min
hence by substitution, (2.27) may be written as Am,~/am,,
=
algm..Jt~mi.] •
it follows that ~max
As n---* ~, cos
0,, =
~N'mm
~6 + ~ COS 0.
cos(~/(n + 1))---> 1, so a ~ 9 . This gives
K(W-"2AW-~'2)~9r,(P-~A)
as n---,~.
(2.28)
We note that Wathen [12] computed the value a = 9 in his analysis. Equations (2.28) and (2.21) imply that r,(P-~A)>~O(n 2) as n - - ~ .
3. Distorted elements in order to study the effect of element distortion on the E B E P C G method and the simple diagonally scaled PCG method, we use the two-dimensional elliptic boundary value problem
-V2u=f
on ~O=[O, or]x[O, 1],
!
(3.1)
with Dirichlet boundary condition on the entire boundary of ~2 as our test problem. Regular grids of rectangular bilinear elements with various element aspect ratios (i.e., the ratio of largest and smallest element length) were used for the analysis. Let the distorteti element be o'h x h and 1 < or < ~:, then the element aspect ratio is defined
H.-C. Lee, A.J. Wathen, Element-by-element preconditioning
222
by trh/h, or simply or. The element coefficient (stiffness) matrix for each element e is then given by ( & =34 k' ~r+
=~ ( ( ~ c r - ) - l ) / 4 ( o r ' - + l )
~'
(-tr2+
½)/4(~r'-+l)
-~l,
(-,r-" + '.)/4(or: + 1)
L
.
(3.2)
where the coefficient k" is constant for all elements e in the case of the constant coefficient problem we analyse here. For a very large aspect ratio (i.e., tr > > 1) the element matrix E e can be expressed asymptotically as
E,,
4k"tr 3
~
.
.
(3.3)
We will later demonstrate from our numerical results that for a very large aspect ratio all performance data do approach to some limit asymptotically as characterised by (3.3). Computations to find the generalised eigenvalues of the symmetric matrix pencil (A - AP) were made and we obtained the largest and smallest eigenvalues of P-~A, where A is the global coefficient matrix for the given problem described by (3.1) and P is the EBE preconditioner defined by (2.5). In all of the computations, the LANSO (Lanczos algorithm with selective orthogonalisation) routines (version 1.0) of Parlett and Liu were used [13]. These results compared with diagonal scaling (W-t~2AW -~'~') for the same problem are tabulated in Tables 1, 2 and 3 for various numbers of degrees of freedom as well as several different aspect ratios. For a fixed number of degrees of freedom or equivalently a fixed mesh size h, the ratio between KOtAC.=- K(W-~I2AW-t'") and Ko~,xc - K(P- ~A) approaches to a limit asymptotically as the element aspect ratio t r ~ 0 . For example, when nnod (number of degrees of freedom) = 102 /(DIAG/KEBE-->8.2 as tr--->oc Table I Numerical results for distorted elements in 2D (nnod" = 10 x 10) Element aspect ratio
Element-by-element preconditioning
Diagonal precondi:ioning KDtAO
(1"
'~min
~max
KEBE
/~min
~max
KDIA(I
KEBE
1 10 l(I-" 10 ~ 104
0.2372 0.1585 0.7106 x 10 i 0.7016 x 10-' 0.7015 x 10- ~
1.0507 1.1705 1.1775 1.1776 1.1776
4.43 7.38 16.57 16.78 16.79
0.5990 x 10 - t 0.4960 x 10- t 0.2136 x 10- ~ 0.21077 x 10 -~ 0.21074 x 10- t
1.4603 2.871 2.899 2.900 2.900
24.36 57.91 135.72 137.57 137.59
5.50 7.85 8.19 8.20 8.20
" N u m b e r of global unknowns (degrees of freedom).
H.-C. Lee, A.J. Wathen, Element-by-elementpreconditioning
223
Table 2 Numerical results for distorted elements in 2D (nnod" = 20 x 20) Element aspect ratio
Element-by-element preconditioning
Diagonal preconditioning KD,^,~
Am,n
Or
1 10 I0-" 10 ~
0.7440 0.8689 0.2043 0.1943
x x x x 0.1940 x
10 4
10- ' 10- ' 10' 10-'
!(1'
J[ma~
KI'BF
1.0507 1.1711 i.1797 1.1798 I. 1798
13.20 13.4~ 57.74 60.71 6(|.81
Am,. 0.1670 0.1670 0.5958 0.5663
x x x x 0.5660 x
10- ' 10 ' 10--" 10 " 10 "
Amax
KI)IA(;
KF.BE
i .4888 2.9427 2.9718 2.9720 2.9720
89.20 176.20 498.77 524.81 525 09
6.28 13.07 8.64 8.64 8.64
" N u m b e r of global unknowns (degrees of freedom).
For nnod = 202 KDIAG/KEBE"~8.64
as o r - - > ~ .
And for nnod = 40-" we have KDIAG/KEBE--'>8.7
as or--~oc.
In the limiting case as or---, ~ the 'elongated' two-dimensional element is degenerated into a one-dimensional type element. Therefore it is expected that the ratio KD,AG/KEBr:should tend to a value of 9 asymptotically as nnod--> ~, which has been proved analytically in Section 2 for the one-dimensional case. As a matter of fact, when nnod = 80", the ratio KDIAG/KEn E is approximately 9. It is interesting to note that the convergence rate for the EBE preconditioning increases at a slightly faster rate than diagonal preconditioning for element aspect ratio smaller than 20. In particular, the EBE preeonditioner becomes especially effective for two-dimensional problems for the element aspect ratio ranging from 10 to 20. We found that KDIAG/KEBE "~- 14.06 when or= 10 and KD,A~/KE, E = 14.12 when or = 20 and deduce that there is an optimal EBE convergence rate for some or, 10 <~ or <~ 20. The relative advantage of the EBEPCG method over the diagonally scaled PCG method diminishes as the element aspect ratio increases. The performance trend is shown in Fig. 1. Table 3 Numerical results for distorted elements in 2D (nnod ~ = 40 x 40) Element aspect ratio or
1 10 10-" 10 ~ 10 ~
Element-by-element preconditioning
Diagonal preconditioning KV~AG
Amin
0.2034 0.2456 0.6208 0.5116 0.5170
x x x x x
10- ' 10- ' 10-: 1O-: 10 -2
Amax
KI=BE
1.0593 !.1764 1.1855 1.1856 i.1856
52.08 47.90 190.96 231.74 229.30
N u m b e r of global unknowns (degrees of freedom).
'~m,n 0.4397 0.4398 0.1799 0.1481 0.15.00
x x x x x
10- 2 10 : l0 : l0-" i0 "
Amaa
KDIA(;
KEBE
1.495 2.961 2.991 2.991 2.991
340.0 673.3 1662.4 2019.6 1994.0
6.53 14.06 8.7 8.8 8.7
224
H.-C. Lee, A.J. Wathen, Element-by-element preconditioning 15
."Ratio of condition t~t~rnbers (DIAG/EBE)
aS 0 ...J.~.nnod~¢o
~9
.
r
n_
.
.
.
.
.
o
B o
~ 7
Element aspect rotto IOgtO 0
Fig. I. Ratio of condition numbers versus element aspect ratio (nnod = 40 ×
41)).
The numerical results indicate that the EBE preconditioning using distorted elements is also spectrally equivalent to diagonal scaling. However, the spectral equivalence bound depends on the element aspect ratio as we have observed in Fig. 1. The EBEPCG method seems to improve slightly over the diagonally scaled PCG method for aspect ratios of or = 1 0 - 20. The ratio Kt~,AG/KEBE is approximately 14.0 for that range as compared with spectral equivalence bound of 6.6 for square elements (i.e., o- - 1). This simply implies that the convergence rate for the EBEPCG method does not get poorer for distorted elements or = 10 ~ 20 as compared to o" = 1. This is in contrast to diagonal scaling which does degrade as or gets larger.
4. Variable coefficients In this section we consider the variable coefficient two-dimensional elliptic boundary value problem -V. (orVu) = f
on O = [0, 1] x [0, 11,
(4.1)
with Dirichlet boundary conditions. Here or(x, y) represents the variable diffusion coefficient. To study in particular the effect of coefficient discontinuities on the EBEPCG method as well as the diagonally scaled PCG method, we use a piecewise constant function tr(x, y) for our investigation. Precisely we take or(x,
Y)
jl, La ,
0~
0~
where the constant o~ may lie in a w!de range of 1 < a < ~.
(4.2)
H.-C. Lee, A.J. Wathen, Element-by.element preconditioning
225
Table 4 Numerical results for variable coefficient problem in 2D # of global
Element-by-element preconditioning
Diagonal preconditioning
unknowns nnod
a
Am,"
Am~x
KERE
21 x 21 (441)
1 10 102 l(I~
0.6813 0.1184 0.1710 0.1306
x x x x
10-' 10- " 10 -4 !0-*
i.057 !.095 !.127 1.138
15.50 9.25 x 10" 659 x 104 8.71 x 10"
41 x 4 1 (1681)
l 10 10"~ 10 ~
0.1940 0.7821 0.9450 0.9495
x x x ×
10 -~ 10 -3 10 -s 10 -7
1.0586 1.1283 1.1457 1.1469
54.58 14.43 x 102 12.12 x 104 12.08 x i0"
'~mm 0.15216 0.15216 0.15217 0.15217
x x x x
Amax Ill- J I0- ' 10-~ 10-i
0.419 x 10--" 0.419 x l0 -2 0.3.045 x 1 0 " I}.4050 x 10-"
KI)IA(;
KO,A. KI~HP-
1.4897 1.5298 1.5572 1.5688
97.90 1(11).54 102.33 103.09
6.32 10.87 x 1 0 -" 15.53 x 10 -4 11.83 x I0-*
!.4958 1.5610 1:5780 1.589()
356.99 372.55 39(}.11 392.34
6.54 25.82 x l0 2 32.18 x 10 -4 32.48 x 10-*
Three different values of a. namely a = 10, 100, 1000 were chosen for the computations, in each case the diffusion coefficient has been characterised by an increasing discontinuity along the line y = ~. In particular when a = 1, it simply reduces to the problem of constant coefficient which was analysed by Wathen [12]. Regular grids of square bilinear elements were used. Table 4 shows the condition numbers obtained for the case of EBE preconditioning and diagonal preconditioning by computing the largest (Amax) and smallest (Am~.) eigenvalues of P-*A and W-'/2AW -I/2. The ratio between KEBE~ K(p-IA) a n d XO~AG- r(W-~/2AW -''2) is also tabulated in Table 4 for comparison with various values of a. The results indicate that for a fixed number of degrees of freedom (i.e., nnod = constant) or equivalently for a fixed mesh size (i.e., h =constant), the condition number KEnE depends strongly on the diffusion coefficient o'(x, y). A simple asymptotic calculation yields KE, E = O(a2). On the contrary, the condition number for the case of using the diagonal preconditioner is independent of the diffusion coefficient. This was also noted by Greenbaum et al. [14]. However, the condition number of the coefficient (stiffness) matrix without any diagonal scaling also depends on the diffusion coefficient, but less strongly than EBEPCG. We have found that x~ -= r ( A ) = O ( a ) from our study as indicated in Table 5. Unlike the spectral equivalences of EBE preconditioning and diagonal scaling which were found for the one-, two- and three-dimensional constant coefficient problems [12], the degradation in performance of the EBEPCG method seemed to increase dramatically as the discontinuous factor a grows. We also obtained estimates of KE, E and KDIAGasymptotically as Table 5 Condition number of A without diagonal scalingin 2D (nnod = 5 x 5) a
Am,"
1 10 102 10 ~ 104
0.192 0.410 0.439 0.442 0.442
~' K, -= , c ( A ) .
Am~,x 1.375 1.310 1.314 1.314 1.314
x x x x
K/' 10 i02 103 104
7.16 3.20 2.99 2.98 2.97
x x × x
l0 l02 i03 10~
226
H.-C.
Lee, A.J.
Wathen,
Element-by-element
preconditioning
nnod ~ oo for a fixed value of a in order to further stt~dy the behavior of discontinuities in the diffusion coefficient. We have as nnod----~oo KEnE = O (nnodP't~)), KOtA~; = O (nodl'zt"),
p2(a) = 1 (constant),
KDIAC/KEBE = O (nnodP~')), where the powers p~, P2 and P3 are plotted against log,, a in Fig. 2. It appears evident that the E B E preconditioning is not spectral equivalent to diagonal scaling as we expected in this case. We want to emphasize that the condition number KEBE grows very rapidly with increasing ot despite the fact of its more mild dependence on the mesh size. For instance, we have /(EBE = O ( h - I )
for a = 1000,
KEBE = O(h -2)
for a = 1 (i.e., with constant coefficient).
Nevertheless, KEnE = 8.714 x 10 6 for a = 1060, n n o d = 4 4 1 , but /(EBE=312.31 for c~ = 1, nnod = 10000. We also present numerical results for the same Dirichlet problem in the three-dimensional problem using uniform trilinear brick elements on ~ = [0, 1] x [0, 1] x [0, 1]. We note that although there are some similarities between this and the 2-D problem just described, they are physically entirely different problems. The numerical results, given in Table 6, shc,w that there is again rapid growth of KEBE with a. For a fixed number of degrees of freedom the condition
t . .......
L .......
t¢. . . . . . . . . . . . . . . .
,It
0.8 07, 06 0.5
~ o.3 ~" 02 0.I 0.-
o
o.'s
I
t
I
,.5
!
2
21s
~Oe~oa Fig. 2. Order of condition number versus 'discontinuous" constant for EBE and diagonal preconditioning.
H.-C. Lee, A.J. Wathen.
Element-by-elementpreconditioning
227
Table 6 3D results for variable coefficient (nnod = 5 x 5 x 5) Element-by-element preconditioning Am,n
'~max
0.6678 0.6373 X 10 .4 0.6195 x 10 -x 0.6167 x 10- " 0.6091 x ! 0 i~
1.0199 1.1408 1.1408 1.1897 1.12288
a
! 10 10-' 10 x !0 ~
Diagonal preconditioning
KEBE 1.527 1.790 !.841 1.929 2.017
x x x x
104 10~ 1012 1(| '~
~m,n
~max
KI)IA/;
0.192 0.192 0.192 0.192 0.189
1.375 !.449 1.472 !.389 1.377
7.16 7.55 7.67 7.23 7.29
Table 7 Condition number of A without diagonal .scaling in 3D (nnod = 5 x 5 x 5) Atom 1 10 102 l0 ~ 10 4 " KI -
0.2750 0.4822 0.5004 0.5021 0.5022
KI"
Amax 1.3499 1.2884 1.2878 1.2878 1.2878
x × x x
10 102 103 104
4.910 2.672 2.576 2.565 2.564
x 10 x l0 a x 10' X 10 a
K(,4).
number in the case of EBE preconditioning applied to the three-dimensional problem is KEnE = O(a4), which is poor even as compared with that of the two-dimensional problem (fEB E = O ( a ~) in 2-D). Spectral equivalence with diagonal scaling which exists for the constant coefficient problem is again lost. However, similar results to the 2-D problem were obtained for both %.AG -- x(W-~t2AW -~'2) and Kt - K(A). In three-dimensions, we have also that KD~AG= O(1) and Kt = O(a) for a fixed number of degrees of freedom. The numerical results are given in Tables 6 and 7. Interestingly, the convergence rate of diagonally scaled PCG is independent of the diffusion coefficient or(x, y, z) as in 2-D. The measure of coefficient discontinuities for both two- and three-dimensional problems was characterised by a single parameter ct in this analysis. For large values of a, e.g., a = 103, and even small values of nnod (120<~nnod<~ 125, say) the condition number for using the EBE preconditioner is O(10 ~') for 2-D and O(1012) for 3-D. The EBEPCG method thus appears impractical for such problems. Three-dimensional results for larger values of nnod were computed, but are not included since they were little different from the numerical results given here. 5. Conclusion We have extended the analysis of EBE preconditioning of the second author [12] to linear elliptic boundary value problems in 2 and 3 dimensions with discontinuous coefficients on regular grids and 2-dimensional constant coefficient problems on grids with variable element
228
H.-C. Lee, A.J. Wathen, Element-by-element preconditioning
aspect ratio. In the latter case we have used bilinear rectangular elements and have shown spectral equivalence of EBE preconditioning and diagonal scaling independently of the aspect ratio. The EBE preconditioner is always better and achieves a much superior performance for aspect ratios between about 5 and 20. For larger aspect ratios the EBE preconditioner is better relative to diagonal scaling than for very small aspect ratios. The case of square elements analysed by Wathen is the worst case. In an absolute sense the EBE preconditioner actually gets better as the aspect ratio of the elements is increased from 1 to approximately 20, and then deteriorates for large values. The condition of the underlying stiffness matrix, of course, deteriorates also as the elements are elongated. For problems with a discontinuous coefficient, we find that the EBE preconditioner deteriorates very quickly as the size of the discontinuiti¢.s is increased. In 3 dimensions, the deterioration is significantly worse than in 2 dimensions. Diagonal scaling appears to be a better choice in such problems. Greenbaum et al. [14] have found similar poor performance on similar problems with the hierarchical basis conjugate gradient method when diagonal scaling is not explicitly used. We note also that Ferencz [11] reports good results using EBEPCG on a problem of thermal stress in a solid state slab laser which involves large material (coefficient) discontinuities. It seems likely that there is very little coupling between the solution in the different material regions in this particular problem and hence the effect of the discontinuities that we have shown would not arise. The test problem we have used is more like that of a beam composed of two equal halves of very different materials. We conclude by emphasising the attractiveness of EBEPCG from the practical point of view. It is highly vectorisable [10] and parallelisable [15] and is derived rationally from the finite element data structure. For problems wh,~re coefficient discontinuities do not play a major role, it appears to be a competitive method, particularly on advanced hardware.
Acknowledgment All computations were carried out in the Numerical Mathematics and Scientific Computation group in the Computer Science Department at Stanford. We are grateful to Gene Golub for use of these facilities and also wish to thank him for his interest and comments which have helped to guide us through this study. We wish also to acknowledge useful conversations with Bob Ferencz and Beresford Parlett for making the LANSO subroutines available to us.
References [i] E Concus, G.H. Golub and D.E O'Leary, A generalized conjugate gradient method for the numerical solution of elliptic partial differential equations, in: J.R. Bunch and D.J. Rose, eds., Sparse Matrix Computations (Academic Press, London, 1973) 163-174. [2] T.F. Russel, M.F. Wheeler and C. Chiang, Large-scale simulation of miscible displacement by mixed and characteristic finite element method, in: W.E. Fitzgibbon, ed., Mathematical and Computational Methods in Seismic Exploration and Reservoir Modeling (SIAM, Philadelphia, PA, 1986) 85-107. [3] T.E. Tezduyar and J. Liou, Element-by-element and implicit-explicit finite element formulations for computational fluid dynamics, in: R. GIowinski, G.H. Golub, G.A. Meurant and J. P6riaux, eds., First International
H.-C. Lee, A.J. Wathen, Element-by-elemem preconditioning
229
Symposium on Domain Decomposition Methods for Partial Differential Equations (SIAM, Philadelphia, PA, 1988) 281-300. [4] O. Axelsson and V.A. Barker, Finite Element Solution of Boundary Value Problems. Theory and Computation (Academic Press, Orlando, 1984). [5] H. Yserentant, On the multi-level splitting of finite element spaces, Numer. Math. 49 (1986) 379-412. [6] J.H. Bramble, J.E. Pasciak and J. Xu, Parallel multilevel preconditioners, Math. Comp. (to appear). [7] T.J.R. Hughes, The Finite Element Method: Linear Static and Dynamic Finite Element Analysis (PrenticeHall, Englewood Cliffs, NJ, 1987). [8] T.J.R. Hughes, I. Levit and J. Winget, An element-by-element solution algorithm for problems of structural and solid mechanics, Comput. Methods Appl. Mech. Engrg. 36 (1983) 241-254. [9] T.J.R. Hughes, I. Levit and J. Winger, Element-by-element implicit algorithms for heat conduction, j. Engrg. Mech. 109 (2) (1983) 576-585. [10] T.J.R. Hughes, R.M. Ferencz and J.O. Hailquist, Large-scale vectorized implicit calculations in solid mechanics on a Cray X-MP/48 utilizing EBE preconditioned conjugate gradients, Comput. Methods Appl. Mech. Engrg. 61 (1987)215-248. [11] R.M. Ferencz, Element-by-element preconditioning techniques for large-scale, vectorized finite element analysis in nonlinear solid and structural mechanics, Ph.D. Thesis, Stanford University, 1989. [12] A.J. Wathen, An analysis of some element-by-element techniques, Comput. Methods Appl. Mech. Engrg. 74 (1989) 271-287. [13] B.N. Parlett and Z.A. Liu, LANSO vl.0 (1989) subroutines, University of California, Berkeley. [14] A. Greenbaum, C. Li and H.Z. Chao, Parallellizing preconditioned conjugate gradient algorithms, Comput. Phys. Commun. 53 (1989) 295-309. [15] B. Gilvary, 1. Gladwell, I.M. Smith and S.W. Wong, The element-by-element method for linear systems, University of Manchester, Numerical Analysis Report No. 144, 1987.