Preconditioned Krylov subspace methods for transport equations

Preconditioned Krylov subspace methods for transport equations

Progress in Nuclear Energy, Vol. 33, No. 1/2, pp. 155-174, 1998 © 1997 Elsevier Science Ltd Pergamon All rights reserved. Printed in Great Britain 0...

1002KB Sizes 6 Downloads 124 Views

Progress in Nuclear Energy, Vol. 33, No. 1/2, pp. 155-174, 1998 © 1997 Elsevier Science Ltd

Pergamon

All rights reserved. Printed in Great Britain 0149-1970/98 $19.00 + 0.00 PII:

S0149-1970(97)00099-1

Preconditioned Krylov Subspace Methods for Transport Equations Suely Oliveira and Yuanhua Deng Department of Computer Science Texas A&M University College Station, TX J u l y 17, 1996

Abstract Transport equations have many important applications. Because these equations are based on highly non-normal operators, they present difficulties in numerical computations. Iterative methods have been shown to be efficient to solve transport equations. However, because of the nature of transport problems, convergence of these methods tends to slow for many important problems. In this paper, we focus on acceleration techniques for iterative methods. Particularly, we investigate the applicability and performance of some Krylov subspace methods with preconditioners, such as the incomplete LU (ILU) factorization (with no fill-in) and multigrid algorithms (spatial and angular multigrid). Three cases are considered: isotropic equations without absorption, isotropic equations with absorption, and auisotropic equations. Our numerical experiments show that the use of an appropriate multilevel preconditioner can significantly improve Krylov subspace methods, such as GMRES and CGS. © 1997 E l s e v i e r S c i e n c e Ltd

§1 I n t r o d u c t i o n Transport equations are mathematical models which describe the transport of particles, momentum, energy, or any transportable quantity (Greenberg, et. al., 1987) They originally arose in the study of electron transport in solid, later were found in many other applications, such as neutron transport, gas kinetics, radiative transfer. Theory has been developed for the study of transport problems, including analytic 155

56

S. Oliveira and Y. Deng

and numerical solutions. It turns out that finding explicit solution is very difficult due to complexity of the transport equations. Although researchers have successfully obtained closed forms of solutions to some particular problems, for most of the problems of practical interest, explicit solutions are impossible to obtain. Therefore, developing numerical methods has become an important aspect in understanding the nature of transport problems. There have been many discretization schemes and numerical methods designed for solving transport problems. Among them, iterative methods have been proved to be one of the most efficient methods. However, for many important problems, the convergence of iterative methods tends to be very slow, for example, for problems with strong scattering. Consequently, there has been a lot of interest recently in the development of acceleration techniques for iterative methods. In this paper, we is a slab of width b we assume that b = differential equation

consider a one-dimensional static model, whose physical domain in the x-dimension (in what follows, without loss of generality 1). In general, such a model can be described by an integroin one spatial variable and one angular variable as follows:

0¢ (x,.) #-~x _ +ot¢(x,,)=S(x,#)+q(x,#),

x e (0,1),

#•

[-1,11

(1)

with boundary conditions given by ¢(0,.)=g0(.),

¢ ( 1 , - # ) = g,(.),

#•(0,1),

(2)

where ¢(x, . ) represents the flux of particles at position x traveling at an angle /9 = arccos(#) from the x-axis, at is the macroscopic total cross section, S(x, #) is a source due to rescattering, q(x, . ) is an external source, and boundary conditions (2) prescribes the flux of particles entering the slab. The source S(x, .) in equation (1) is given by an integral operator of the following form:

S(x,.) = f l P(x,

(3)

where P is the probability density that a particle will be scattered from direction .~ to #. If the probabilities of scattering of particles are the same for all directions (called the isotropic case), then P in the above integral is constant in angle and S(x, #) is simply expressed as

S(x, ~)

: ~o- s

1 ¢(x, #')d.',

(4)

where crs dx is the expected number of scattering interactions in a length dx. Thus, (ot - a~) dx is the expected number of absorptive interactions in length dx. Then, this term is actually independent of #. In recent years, much effort has been invested in the development of numerical methods for solving transport equations, especially in multigrid or multilevel algorithms. For example, Morel and Manteuffel (1991) presented angular multigrid techniques for anisotropic problems with highly forward-peaked scattering. Based on the

Preconditioned Krylov subspace methods

modified linear discontinuous (MLD) scheme, Manteuffel, et al. (1995) developed a fast multigrid algorithm for isotropic problems with pure scattering. Multigrid algorithms for two dimensional transport equations were investigated by Nowak (1988, 1988b). Parallel multigrid algorithms for solving either isotropic and anisotropic equations were also developed, see Manteuffel, et. al. (1994, 1995) and Oliveira (1993, 1994, 1996). It is noted that the discretization of transport equations in general will lead to nonsymmetric linear systems, which make Krylov subspace methods particularly attractive in numerical computations. In addition, Krylov methods allow the easy incorporation of preconditioners. Some studies for applications of Krylov subspace methods to transport equations have been presented recently. Kelly and Xue (1996) studied the performance of the GMRES iteration for a general integral equation. Then they used the isotropic transport equation as an example for numerical experiments, where the Crank-Nicolson scheme was used for the spatial approximation. Patton and Holloway (1996) examined the application of GMRES with an ILUT preconditioner to slab geometry transport problems. They used the diamond difference discrete scheme. In this paper we investigate the applicability and performance of some Krylov subspace methods with the spatial and angular multigrid preconditioners for solving transport equations. In our work we use an MLD discretizations and a different kind of splitting (for the discretiazation), called tt relaxation, which is suitable for a multigrid preconditioner. Three cases are examined: isotropic equations without absorption, isotropic equations with absorption, and anisotropic equations. Our results will be compared with the incomplete LU factorization (ILU) preconditioner to show the improvement. It is worth remarking that the multigrid method used in Manteuffel, et. al. (1994, 1995), Oliveira (1995) has difficulties with isotropic scattering with small but significant absorption. We show that if multigrid methods are used as preconditioners of some Krylov subspace methods, such as CGS and GMRES, the convergence rate can be significantly improved, even for this case. The rest of this paper is organized as follows: In Section 2 we briefly describe the discretization used. In Section 3, we discuss Krylov space methods and preconditioners. In Section 4, we give some implementation details for the angular multigrid method with the MLD discretization scheme. In Section 5, we show our numerical results where a comparison of these methods is also given. Finally, some concluding remarks are given in Section 6. §2 Discretization In this section, we briefly show the discretization of transport equation (1). It is noted that equation (1) involves two variables: spatial variable x and angular variable #. The angular discretization is accomplished by expanding the angular dependence in Legendre polynomials and collocating at Gauss quadrature points, see Lewis and Miller (1984). More precisely, letting Pc(#), g -- 0, 1 , . . . , N - 1, be the first N Legendre polynomials, we assume that ¢(x, it) is discretized by the PN-1

157

58

S. Oliveira and Y. D e n g

approximation as follows: N-1

¢(x, #) = ~ (2g + 1)¢e(x)pe(#),

(5)

l=0

where ¢ ~ ( x ) , ~ = 0 , . . . , N - 1, represent Legendre moments given by 1

1

¢~(x) = ~ f-1 pl(p')¢(x, if)d#'. Hence, scattering term S(x, #) and moment term ¢~(x) can be expressed by N-1

S(x, #) = ~ (2k + 1)crkp~(#)¢k(X),

(6)

k=0 N

Cd~) = ~ w~p~(~)¢(z, ~).

(7)

k=l

where #i and wi, i = 1 , . . . , N , are the Gauss quadrature points and weights, respectively, and ak is the kth Legendre moment of the scattering cross section. If we consider the isotropic case, from (4) we readily see that N

s(~, .)

= o~ ~

w~¢(~, ~).

(8)

k=l

Now collocating equation (1) at points #i, i = 1 , . . . , N for angular variable # and employing notation

_¢(x) = [ ¢(x,~1) ¢(x,#2) ... ¢(x,,N)]T, q(x)=[ q(x,#l) q(x,#2) ... q(x,#g) IT, we can get the following system of ordinary differential equations:

M ~ x (X) + atC__(x) = T-1DTC__(x) + q(x),

(9)

where M, D, T, and T -1 are N x N matrices given by M = diag{#l, # 2 , ' " , #n}, D = diag{ao,~h,... , a N - l } , and T = (ti,j) and T -1 = (ti,j) with entries h,j and ti,j given by ti,j = wjpi-l(#j) and ti,j = (2j - 1)pj-l(#i), respectively. We then perform the spatial discretization on equation (9) by using the MLD scheme. This process is similar to that described in Larsen and Morel (1982) and Manteuffel, et. al. (1995). Thus we omit the procedure for derivation here. The reader can refer to these references for details. To present our final discretization of transport equation (1), we need to introduce some notation, which will be used throughout this paper unless otherwise stated. We assume that the degree N of Gauss quadrature is even and N = 2n. Thus, the quadrature points are symmetric about zero. Assume that

Wk = Wk+n

and

#k = --#k+~ > 0,

k ---- 1,. •., n.

Preconditioned Krylov subspace methods

159

Let {xi+½}~m__0be a set of grid points distributed over interval [0, 1] such that

½ =1.

0 = x ~ 2 < x 32 < - " < X m +

1 Let xi = 7(xi+½ + xi_½), the center of cell i, and hi = xi+½ - xi_½, the width of + cell i. For cell i, we denote by variables ¢i--½d, ¢i,+, ¢ ~ , and ¢~+½,j the approximate values of ¢ ( x i _ ~ , - / z j ) , ¢(xi,/~j), ¢ ( x i , - # j ) , and ¢(xi+½,#j), respectively, where j = 1 , 2 , . . . , n . ~et

...

¢;,2.

for k = ½,1,35 , . . . , m , m + ~ . 1 Proceeding as in Manteuffel, et. al (1995) from equation (9) we can obtain the following discrete equations

2Mi(¢-f_½ - _¢'f) -~- o't~'f_½ = _

Mi(¢__++½ - ¢__,+½) + ate__+

S 2 1 ( 2 ¢_/ +

= Sll~:

Mi(¢__~_½ - ¢__~-+½)+ at¢__~- = $21¢ +

- _,+~, ¢+ 1

+ $22¢__~-_½+ q~_½,

-~- S 1 2 ¢ ?

"4- q : ,

--b S 2 2 ~

+

(10) (11)

qi-,

(12)

2M,(¢_£+½ - ¢_£) + ot¢__,++½= S11¢__,++½+ S,2(2¢_,- -_~:_½) + q,++½,

(13)

/---- 1 , . . . , m , with boundary conditions ¢+ = 9_o +,

--2

and

¢_~,+½ -- g_~-,

(14)

where Mi = ~diag{#l, # 2 , ' " , #n}, Sn, $12, $21 and $2~ are n x n blocks of matrix S =- T - I D T such that

S=

[Sn$12] S2, S n

and go+ = [ g0(#l) "'" g0(#n) ]T, g+ = [ g1(#1) isotropic case, we have that

S11=S12

=

S2,=S2~=a,

1:

'

"'" gl(#n) ]T. Note that in the

W1

• . .

?Mn

.

Generally speaking, inhomogeneous terms q~-j, k = i - 3,1z- and qkj,+ k = i, i + ~1 in equations (10) - (13) can be taken to be any approximate values of q(xk,--#j), k ----i - ½, i and q(xk, #j), k = i, i + ½, respectively. However, for the implementation of the angular multigrid method which will be discussed in Section 4, we need to carefully examine the discretization procedure to get precise forms. It is not difficult

160

S. Oliveira and Y. Deng

to find that qi-~,j = ~imq(2+, -#J),

(15)

fx,+½ q(x, #j) dx, qi,+ = ~1 .=,_½

(16)

1 [x,+½ q ( z , - # j ) dx,

(17)

qi+[,j = ~in~q(2(, #j),

(18)

where 2 + = xi_ ½ + e and 2 i- = xi+ ½ - e, see Manteuffel, et. al. (1995) for details. Letting =

_,+,

,

q:

q.

q,+l

,

we may rewrite equations (10)-(13) in a matrix form: Aiff2i = Bi~..~_l + Cik~i+l --b Q_.~,

i = 1, 2 , . . . , m,

(19)

where matrices Ai, Bi, and Ci are given by

Ai

atI + 2Mi - $22 -2S21 0 atI - Sn Mi -$21 S,~ -2Mi

Bi

000

0

O00Mi 000 000

0 0

-2Mi -S12 atI - $22 -2S12 atI +

G =

0 0 Mi 0

$21 Mi 0 2Mi - Sn

(20)

0 0 0 0 0 0 000 000

Figure 1 shows the computational domain with 2m + 1 spatial points and n angles. We see that for cell i, terms Bi~i-1 and Ci~i+l in equation (20) are only related to the inflows from its left and right ends, respectively. Hence, if the inflows of a cell are assumed to be known, matrix equation (19) can be solved by performing a #-line relaxation on the cell, that is, center variables ¢+ and ¢~-, together with the outflow variables ¢~-_~ and _¢++_~are updated by using (19). It is noted that boundary conditions (14) represent, in the computational grid, the inflow for positive angles on the left and for negative angles on the right of the domain. We may perform this iterative process in both directions from left to right and from right to left. This kind of relaxation will be used in the implementation of smoothing operations for the angular multigrid method, which will be seen in Section 4. On the other hand, to apply Krylov subspace methods, we need to consider all the cells coupled together.

Preconditioned Krylov subspace methods

. . . . J. . . . . t .

Xl/2

X I

. . ..J . . .

...

L . . . .,J . . . . t .

Xi_l/2

Xi

Xi+l/2

161

. . .J . . . . . .L. . . . . . . . . . .

Xi+l

..-

x m

Xm+lf2

t__J Cell i

Figure 1: Computational grid. Hence, we derive the following linear system for the whole domain: A1

-C 1

-B~

A2

-C2 :

°•

-Bin-1

+

,

(21)

0

-Bm

Am

bm

where b~=[O

M1¢_+½ 0 O] T

and

b_~_~---[0 0 MmC__m+½ 0 ]

T

which are known due to boundary conditions (14). Now we can apply an iterative method to equation (21). §3 K r y l o v Subspace Methods a n d P r e c o n d i t i o n e r s In Section 2 we have seen that the coefficient matrix for linear system (21), which is derived from discretization of either isotropic or anisotropic equations, is block tridiagonal. Such a matrix is sparse and, in general, nonsymmetric. We intend to find efficient numerical methods for solving system (21). There are many iterative methods for solving large sparse linear systems A x = b, where A is a d × d nonsymmetric matrix and x, b are d-dimensional vectors• One class of iterative methods is Krylov subspace methods• To solve linear systems derived from transport equations, we need nonsymmetric Krylov subspace methods, such as CGS (Conjugate Gradient Squared) (Sonneveld, 1989), Bi-CGSTAB (Bi-Conjugate Gradient STABilized) (Van der Vorst, 1992), GMRES (Generalized Minimal RESidual)

162

S. Oliveira and Y. Deng

(Saad and Schultz, 1986), T F Q M R (Transpose-Free Quasi-Minimal Residual) (Freund, 1993), CGNE (Conjugate Gradient for Normal Equation) (Hestenes and Stiefel, 1952), and LSQR (Least Squares/QR) (Paige and Saunders, 1982). In these methods, iterative solutions are of the form x (~+0 = x (k) +p(m), where p(k) is in the Krylov space Kk(A, r °) = span{r °, A r °, A2r°, .. ., A k - l r °} with r ° = b - A x ° and k _< d. To obtain bases of Krylov subspaces, the methods employ either the symmetric or unsymmetric Lanczos methods or the Arnoldi method (Golub and Van Loan, 1989). For instance, CGS method uses unsymmetric Lanczos, GMRES method applies Arnoldi, and LSQR utilizes symmetric Lanczos. The symmetric Lanczos and Arnoldi algorithms generate (in exact arithmetic) an orthonormal base of space Kk(A, r°), while the unsymmetric Lanczos produces a pair of biorthonormal bases of spaces Kk(A, r °) and K k ( A T, r°). In both cases the Lanczos method yields a tridiagonal matrix which represents the original matrix on the Krylov subspace, while the Arnoldi method gives a Hessenberg matrix which represents the matrix on the Krylov subspace K~(A, r°). The unsymmetric Lanczos process is fast, but may suffer from numerical instability, known as breakdown. However, there are variants of the unsymmetric Lanczos method, based on the look-ahead Lanczos algorithm, which are stable. One advantage of Krylov subspace methods is to allow the easy incorporation of preconditioners. For system A x = b, a preconditioner is a matrix B such that B u can be easily computed for a given vector u and system B A x = Bb is easier to solve than the original one. Usually, to this end one intends to find B such that matrix B A is well conditioned. Such an appropriate matrix B can be obtained by a number of different approaches. One of the simplest preconditioners is the diagonal preconditioner, where B is simply taken as the inverse of the diagonal of A. The diagonal preconditioner works well if A is diagonally dominant. Another type of preconditioner is to utilize incomplete factorizations of sparse matrices, for example, ILU (Incomplete LU factorization). Some iterative methods can also be used as a preconditioner, such as multigrid methods. Here we shall investigate the use of multigrid methods with respect to the spatial variable (called spatial multigrid methods) and the angular variable (called angular multigrid methods) as a preconditioner. It is noted that Krylov subspace methods only require three operations for their implementation: linear combinations, inner products, and matrix-vector products. Thus, for a preconditioner which is incorporated into a Krylov subspace method, it is convenient to use a routine to compute B A u for a given vector u by first computing v = Au, and then using a routine for computation of By. §4 M u l t i g r i d M e t h o d s We now turn to the discussion of multigrid methods. The implementation of the spatial multigrid method for isotropic transport equations with the MLD scheme has been discussed in Manteuffel et. al. (1995). An angular multigrid method for anisotropic transport equations was developed in Morel and Manteuffel (1991). Here we show some implementation details for that algorithm but for a particular spatial discretization MLD. We also note that the parallel implementation of the angular

PreconditionedKrylovsubspacemethods

163

multigrid method for anisotropic transport equations with a spatial discretization similar to the MLD scheme has been presented by Oliveira (1996). Let us first recall the angular multigrid method presented by Morel and Manteuffel (1991). For the sake of simplicity, we denote by E a differential operator defined by

£,=Md__d_ dx + o'tI, where M = diag{#l,. • -, tin, -it1," • ", - # n } and I is an identity matrix. Then, equation (9) can be rewritten as (/~¢__)(x) = S¢__(x) + q(x),

(22)

where S =- T-1DT. We also introduce notation

_¢(X)-~-[ C0(X) el(X)

"'" CN-I(X)IT.

In view of (5) and (7), we easily see that

¢_(x) = TC_(x),

¢_.(x)= T-lC_(x).

(23)

To reflect the multilevel, in what follows we shall use subscript k to denote the grid level. Thus, equation (22) may be rewritten as = Sk

(x) +

(24)

with Sk = T[IDkTk and Nk = N. Given an initial guess _¢o at the finest grid level L, the V-cycle of the angular multigrid method starts with performing a small number of transport sweeps for equation (24) (known as pre-smoothing in the multigrid terminology) and then transferring the residual [L of the approximate solution to rL_ 1 at the coarse level (known as restriction). This process is repeated on successively coarser grids with zero initial guess and q-k in equation (24) replaced by (1 < k < L - 1). On the finest grid, equation (24) becomes isotropic. Thus, it can be solved by any method for the isotropic equation and its solution is then transferred upward to the fine grid (known as prolongation). On the successively finer grids, the solution transferred from the coarse grid is used to correct the solution obtained by the pre-smoothing at the same grid level. The resulting solution is then smoothed by performing transport sweeps again (called post-smoothing). The above V-cycle process can also be described in a recursive form so that it is easily implemented by a programming language which facilitates recursive function calls, such as C language. To this end, we denote by AngMGk(C~°,q~) the angular multigrid method starting at level k and an initial guess C__°k. AngMGk(C°k, qk) is defined in terms of AngMGk_l(¢__°_l, qk-1) as follows. If k = 1 (i.e. N1 : 2), then equation (24) becomes isotropic. Thus, AngMG1 can be chosen to be any method which solves isotropic equations, for example, spatial multigrid methods.

64

S. Oliveira and Y. Deng

If k > 1, AngMGk is defined in terms of AngMGk_~ as follows: 1. Pre-smoothing, which is given by the following relaxation scheme: (/:kC_~g+l))(X) = SkC__~t)(X)+ qk(x),

g = O, 1 , . . . , v

-

1,

(25)

with an initial guess ¢_~(0) :__ _¢~. The pre-smoothing yields the approximation solution ¢__~"). In view of (23) - (25), the corresponding residual ~ ( x ) can be calculated as follows: r__k(X) = qZc(X) -- [(£k_O~"))(X) -- SkO_~O')(X)]= SkAC__k(X) = TklDkAC_k(X ).

(26)

where

=

-

A (x) = Tk C_k(x )

2. Correction, which consists of three major steps: restriction, recursive correction from coarser grids, and prolongation. The restriction operator from level k to level k - 1 is defined by removing a half of the higher moments in AC_k(X). Let Nk-1 = Nk/2 = nk. Using (23) and (26) we can obtain that

= T;-21D: [ I _1 0 ] T, AC_,(x)

(27)

where D + = diag{c%o, ak,1,... ,ak,nk-1}, Ik-1 is an identity matrix of dimension Nk-1 x Nk-1, and [ Ik-1 0 ] is an Nk-1 x Nk matrix. Next, by applying AngMGk_I to equation (24) at coarse level k - 1 with an initial guess --¢~-1 := 0 and qk_l(X) := _rk_1(X), we obtain a solution denoted by AO_k_l(x ). Finally, a prolongation operator is defined by transferring the flux AO_k_l(X ) from level k - 1 to level k through the moment AC~_l(X ) = Tk_lAC_k_l(x). We thus get a corrected approximate solution

3. Post-smoothing, which is similar to pre-smoothing, but with ~k(x) as an initial guess. The resulting solution from the post-smoothing is then returned as the final solution for AngMGk. It is noted that the angular multigrid method angMGk is concerned only with the angular variable, thus, independent of spatial discretization. However, in practice, the angular multigrid method has to be incorporated with a spatial discretization scheme. Here we consider the MLD scheme and show some details for the implementation. Examining the angular multigrid method defined above, we see that it will be sufficient to derive discrete versions of the smoothing, restriction and prolongation operators, Without loss of generality, we assume zero boundary conditions in our discussion below.

Preconditioned Krylov subspace methods

165

Let us first look at the smoothing procedure• Note that matrix can be decomposed into two parts Ai = Hi - E, where

atI + 2Mi 0 M~ 0

0 crtI 0 -2M~

-2Mi 0 atI 0

$22 -2S21 0 0 --Sll -$12 0 - $ 2 1 -$22

E =

S12

0

Ai defined in (20)

0

Mi 0

atI + 2Mi

$21 0 0

-2S12

-Sll

Then, equation (21) can be reformulated as -

Bi~i-1 + Hi~i

Let

~---= [ ~---1 kI/2

" l'I/---'n] T'

(29)

i = 1,2,...,m.

Cik~i+l = E ~ i + Qi,

-

Q__~ ] r ,

Q__: [ ~_-1 ~_-2 "'"

and

H1

- C1

-B2

//2

L=

E

-C2

E

°•

• •.

-Bin-1

,

Hm-1 -

Y

~

-Cm-1

Bm

E

Hm

Corresponding to equation (25), we get its discrete version given by L~(t+l) = V~(t) + 9@

e = 0, 1 , . . - , u - 1.

(30)

In each iteration, we axe required to solve a linear system Lq' = F.

(31)

However, such a system can be solved directly. To see this, we introduce notation

and _• -

--

-..

]T

=

if

...

We observe that equation (31) can be deeoupled with respect to unknowns ~ and ~+, respectively• In fact, we can obtain the following two linear systems: H1

-C1

~J1

• ..

... H~_I

F1 =

-Cff,_l

kg~ 1

,

F---~-I

(32)

166

S. Oliveira and Y. Deng

%+ _%+

H1+

_B +

H2

_F,+ __F~+ (33)

--

B m+- - 1

Hrn+_ 1

-B +

II/rrt -- 1

Fro- 1

H~

where Mi

I

B+=

[0 ~]

and

0

'

0

H+=

-2Mi

C~-=

I0 0] Mi

0

I+2Mi

'

,

The inverses of matrices Hi- and H + can be easily obtained• As a result, equations (32) and (33) can be solved by block backward and forward substitutions, respectively. Next, we consider the implementation of the restriction for MLD scheme• In order to apply the angular multigrid method at the coarse level k - 1 we need to get the discrete version of the residual ~_l(X) given by equation (27). To this end, we write vector r_k_l(x ) as

r~ ~x~-Ir~+i~ ] _%l(X)

"

Denote by variables rk-__lj, j = i -- 1, i, and -r+ k - l , j J = i,i + ~1 the approximate values of function rk_l(X ) at points x j , j ---- i -- ½, i and function r_+_l(x) at points xj, j = i, i + 1, respectively• Applying equations (15)-(18) to (27), where AC_k(X ) is represented in finite element spaces defined in Manteuffel, et. al. (1995)• We can obtain the following discrete version of the residual ~-l(X):

r__k_l,i_½ =

[/~0~] [0 ~i~ -

Zk

~:1~+0_[0] --

and

Ik- 2

[~+ ] --k- l,i

r_k--1,i

Zk

= Zk

ik_l

0

[ 0 0 0 --Ik_ 1 0

]

0 -Ik-10 A~,i,

0

2Ik-i

Ik-10] A~_~,i,

[o,~_1 o o] ~,~, 0 0 Ik-1 0

where

Similarly, we apply equations (15)-(18) to the second term in (28) to get the discrete version of the prolongation operator.

Preconditioned Krylov subspace methods

§5 N u m e r i c a l r e s u l t s In this section, we present some numerical results about convergence, when we solve linear systems derived from isotropic and anisotropic equations by using various Krylov subspace methods, such as GMRES, CGS, and CGNE, with ILU and multigrid preconditioners. For the ILU preconditioner we use the ILU with no fill-in (Stewart and Z. Leyk, 1994). For the multigrid preconditioner, we use either spatial or angular multigrid. In the angular multigrid, one transport sweep is used in both pre-smoothing and post-smoothing steps and the spatial multigrid method is used as the solver on the coarsest grid. The test problems were selected by randomly generating the right-hand side vectors of linear systems and the numerical experiments were all performed on a Sun SPARCstation 10 in the Department of Computer Science, Texas A&M University. The workstation runs Solaris 2.5 and has a 32Mb RAM and a 64Mb swap space. We used the cc compiler SC4.0 without any flags. We run our codes on problems with various sizes. However, as all the problems that we tested exhibited similar behavior, we only present results for a case where n = 8 (16 angles) and m = 64 (64 cells), which yields a discrete problem with dimension 2048. Moreover, unless otherwise specified, the GMRES method used in our examples is the restarted version GMRES(k) with k, the number of directions, set to 50. We will show residual behavior during an iteration process for each of the following three cases: (1) isotropic equations without absorption (~ = 1, where 7 = as o~), (2) isotropic equations with absorption (~/ < 1), and (3) anisotropic equations. In all cases, the iteration process is terminated as long as the Euclidean norm of the residual is reduced by a factor of 10 -14 (which is nearly the machine precision) or the number of iterations exceeds 300. 1. I s o t r o p i c e q u a t i o n s w i t h o u t a b s o r p t i o n . The non absorption case is the simplest of three cases considered. Residual reduction for each iterative method is plotted, respectively, in Figures 2-4, where the x-axis is the number of iterations and the y-axis is the base 10 logarithm of the Euclidean norm of residuals. The multigrid method considered here is the spatial multigrid method developed by Manteuffel, et. al. (1995), which is denoted by MG in these figures. From Figures 2 - 4, it clearly seen that convergence for GMRES(k) and CGS is significantly improved if either ILU or MG preconditioner is used, but not for CGNE method. Comparing ILU and MG preconditioners, we see that the latter is more effective than the former when they are applied to either GMRES or CGS method. For example, we computed an average convergence factor (ACF) for preconditioned GMRES methods. Here the ACF is defined as ! ~ = 1 ~ , where ri is the residual in the ith iteration. The average convergence factor for t ~ preconditioner is 0.459, whereas the average convergence factor for MG preconditioner is 0.014. In addition, it is known that convergence rate of GMRES(k) is dependent on k. As shown in Figure 5, the number of directions in GMRES(k) can be considerably reduced if MG preconditioner is applied, in comparison with ILU preconditioner. Note that the curves for MG(k = 20) and MG(k -- 50) in Figure 5 are indistinguishable. Since the spatial multigrid method itself can be used as a solver in this case, to analyze the performance of Krylov meth-

167

168

S. Oliveira and Y. Deng

ods with the multigrid preconditioner, we plot the convergence behavior of the MG solver, CGS and GMRES with MG preconditioner in Figure 12, which demonstrates that convergence can be improved if MG is used as a preconditioner of CGS and GMRES. Especially, such improvement is significant if one uses CGS with MG preconditioner. Corresponding to Figure 12, Table 1 (top) shows average convergence factors for all the methods. From this table, we observe that CGS with MG preconditioner presents the convergence factor of order around 10 -4, in comparison with MG whose convergence factor is of order 10 -2 . Cases Case 1

Case 3

Case 3

Methods MG GMRES-MG CGS-MG GMRES-ILU GMRES-MG GMRES-AMG AMG GMRES-AMG CGS-AMG

#iterations 8 6 4 283 87 16 28 16 9

ACF 0.73E-2 0.47E-2 0.23E-4 0.89 0.72 0.11 0.29 0.11 0.04

Table 1. Average Convergence Factor (ACF). 2. Isotropic equations with absorption. The absorption case is of great practical interest and it is not hard to solve if the absorption is sufficiently large. However, with a small amount of absorption, it presents computational difficulties. Here we present numerical results for the absorption regime given by 7 = 1 - ~ with at h -- 100. This regime corresponds to a situation in which the scattered particles undergo a large number of scatterings. In addition, the particles have significant probabilities both of being absorbed, and of "escaping" from a cell, although the probability of absorption in each collision is very small. The results are plotted in Figures 6-8, from which we see that the convergence behavior is similar to the first case. Again, preconditioned CGS and GMRES display very good improvement in convergence rate, especially, when MG preconditioner is used, whereas CGNE is not convergent at all, even with the preconditioners.

Preconditioned Krylov subspace methods

Io+

1o~

.

.

169 .

.

.

1°~ ~

10o

!

1o4 ~ I04 fi~l

10-'~ I

.... GMI~-.~It mlmblm" ofmmllmml

numblrofier~

Figure 2: GMRES (case 1).

Figure 3: CGS (case 1).

l°a[

- - CG~IE

lO'

101

-- CGNE-II,U .... CGNIE-MG

loo

,z' i

|,oLi )

"~

10

Figure 4: CGNE (case 1).

30

40

SO

e0

7o

Figure 5: GMRES,k=20,50(case 1). 104

loa

2o

.

.

.

.

.

.

lo' ~

10-*

$,o+

,

",

^A

!im+,~ 0"~

Figure 6: GMRES (case 2).

'

Figure 7: CGS (case 2).

i

S. Oliveira and Y. Deng

70

10 ~

--

CGINE

--

CONE-ILU

Io ~ 1°°

~10' 8 J~10*

10 -¸

10-'4 ~ 5O

100 ~mt~r

150 of .watk~=

2OO

25O

o

3OO

Figure 8: CGNE (case 2).

QMRES-AMG

so

1~

150 . u m b e r of ~ , r ~ t ~ ,

250

Figure 9: GMRES (case 3).

10 Io

l o '0 Io 6

lO l

~ :: lOe lo"

i

lo ~

"~1o °

~ 1o-2 If' l O .

m4

50

~50 number of It~otion~

V

lO-,~

250

50

100

1S0 number of Itemt lo~s

200

2~50

300

Figure 11: CGNE (case 3).

Figure 10: CGS (case 3). lO~

1o o



--

:;,~- • ~. ~ •

10'

~

!

~

~

1o~

--

~

AMG GMRES-AMG CGS-AMG

1o~ 1o~

~ 1 o "1o

lO-,O

--

GMRES-MG

- * CGS,-MG --

:NE-MG

10-*l

lO-1. lo-,I

2

4

e

8 1'0 1~ number o~ IterltJon=

Figure 12: MG (case 1).

s

10

~5

2O

25

Figure 13: AMG (case 3).

3. A n i s o t r o p i c equations. The last case we considered was the anisotropic transport equations. In our numerical experiments, the Legendre cross section moments are taken as those for Fokker-Planck scattering (Morel and Manteuffel, 1991),

Preconditioned Krylov subspace methods

i.e. crk = c~((N - 1)N - k(k + 1))/2, where we choose c~ = 1. The behavior of the residual is demonstrated in Figures 9-11. The anisotropic case is the most difficult case for numerical computations among the three cases considered. We see that convergence rates deteriorate for all the methods used in cases 1 and 2. Thus, in addition to ILU and MG preconditioners, where the MG is the spatial multigrid method for isotropic equations without absorption, we consider the angular multigrid (AMG) as a preconditioner, since the angular multigrid is specifically designed to accelerate problems with highly forward-peaked scattering (Morel and Manteuffel, 1991). Prom Figures 9 and 10, we observe that convergence of GMRES and CGS is improved using all three preconditioners. Although angular multigrid preconditioner increases the complexity of the algorithm and amount of computation, it is more efficient than ILU and MG preconditioners in overall performance. As an example, Table 1 (middle) shows the convergence factors for preconditioned GMRES. In this table we show that convergence factors for GMRES with ILU and MG preconditioners are 0.89 and 0.72, whereas convergence factor for GMRES with the angular multigrid preconditioner is 0.11, as we should expect. It is noted that using MG preconditioner still accelerates convergence much faster than ILU for the GMRES method (see Figure 9). However, unlike the previous two cases, ILU preconditioner works better than MG preconditioner when they are applied to the CGS method, as seen in Figure 10. This is because MG we used is specifically tuned for isotropic equations. In Figure 13 we plot the convergence behavior of the angular multigrid as a solver and as a preconditioner for CGS and GMRES methods to see the performance of preconditioned CGS and GMRES methods. The corresponding convergence factors are also computed and displayed in Table 1 (bottom). We see that preconditioned CGS and GMRES methods show significant improvement over the angular multigrid, especially, CGS with AMG preconditioner. §6 C o n c l u d i n g r e m a r k s In this paper, we investigated the performance of preconditioned Krylov subspace methods for solving transport equations. We tested three Krylov subspace methods GMRES, CGS, CGNE with ILU (no fill-in), spatial and angular multigrid preconditioners. Our numerical experiments shows that both preconditioned GMRES and CGS significantly accelerate the convergence of the methods on all three cases we consider, although the convergence rate deteriorates from case 1 to case 3. Especially, for isotropic equations (cases 1 and 2) the spatial multigrid preconditioner is superior to ILU conditioner, while for anisotropic equations the angular multigrid preconditioner is the best preconditioner to use, as we should expect. In addition, both GMRES and CGS methods with the spatial multigrid preconditioner perform better than the spatial multigrid method itself for solving isotropic equations. Similarly, GMRES and CGS with the angular multigrid preconditioner overperform the angular multigrid method for solving anisotropic equations. Our numerical results indicate that CGNE method is not convergent with either ILU or multigrid preconditioners. The reason

171

172

S. Oliveira and Y. Deng

for this may be related to the worse conditioning of the normal equations used in CGNE. It is noted that preconditioning will increase the amount of computation. Thus, it trades the reduction of iterations with the cost of more expensive operations. To evaluate overall performance, we may also need to consider the CPU time. Table 2 shows the CPU time required for all three cases considered in Section 5. It is noted that CGNE is not listed in Table 2 because it either does not converge or converges very slowly. Methods GMRES

CGS

MG AMG

Precond None ILU MG AMG None ILU MG AMG None None

Case1 17.04 3.11 2.21 11.63 4.19 2.44 -

2.45 -

Case2 5.35 5.66 5.44 2.25 4.37 6.78 -

-

Case3 36.11 36.80 67.26 7.73 48.94 22.82 193.86 8.01 13.87

Table 2. CPU time (in seconds), n = 8, m = 64, From Table 2, it is observed that for case 1 preconditioned GMRES and CGS still produce more or less reduction in the CPU time. However, for case 2, GMRES and CGS perform better than their preconditioned counterparts because the expensive construction of the preconditioners. The ILU and MG preconditioners do not work well for case 3 either. However, the use of AMG preconditioner achieves a significant reduction in the CPU time for case 3. It is seen that in comparison with AMG, the AMG preconditioned GMRES and CGS reduce the CPU time almost by a factor of

2/5. REFERENCES

[1]

Faber V. and Manteuffel T. A. (1989) Neutron Transport from the Viewpoint of Linear Algebra. Lecture Notes in Pure and Applied Mathematics.

[2]

Freund R. W. (1993) A Transpose-free quasi-minimal residual algorithm for nonhermitian linear systems. SIAM J. Sci. Comput. 14, 470.

[3]

Greenberg W., van der Mee C. and Protopopescu V. (1987) Boundary Value Problems in Abstract Kinetic Theory. Birkhauser, Basel.

[4]

Golub G. H. and Van Loan C. F. (1989) Matrix Computations. The Johns Hopkins University Press, Baltimore and London.

Preconditioned Krylov subspace methods

[5] Hestenes M. R. and Stiefel E. (1952) Methods of conjugate gradients for solving linear systems. J. Res. Nat. Bur. Standards 49, 409. [6] Kelly C. T. and Xue Z. Q. (1996) GMRES and integral operations. SIAM J. Sci. Comput. 17, 217. [7] Larsen E. W. and Morel J. E. (1982) Asymptotic solutions of numerical transport problems in optically thick diffusive regimes II. J. Comput. Phys. 83, 212. [8] Lewis E. E. and Miller W. F. (1984) Computational Methods of Neutron Transport. John Wiley and Sons, New York. [9] Manteuffel T., McCormick S., Morel J., Oliveira S. and Yang G. (1994) A Parallel Version of a Multigrid Algorithm for Isotropic Transport Equations. SIAM J. of Sci. Computing 15, 474. [10] Manteuffel T., McCormick S., Morel J., Oliveira S., and Yang G. (1995) Fast Multigrid Solver for Transport Problems I: Pure Scattering. SIAM J. of Sci. Computing, 16, 601. [11] Manteuffel T., McCormick S., Morel J. and Yang G. (1993) Fast multigrid solver for transport problems with absorption. Technical Report, University of Colorado. [12] Morel J. E. and Manteuffel T. A. (1991) An angular multigrid acceleration technique for the SN equations with highly forward-peaked scattering. Nucl. Sci. Eng., 107, 330. [13] Nowak P. F. (1988) A coupled synthetic and multigrid acceleration method for two-dimensional transport calculations. Ph.D. thesis, Department of Nuclear Engineering, Univerisity of Michigan, Ann Arbor, Michigan. [14] Nowak P. F. (1988) A multigrid method for SN calculations in x-y geometry. Trans. Amer. Nucl. Soc. 56, 291. [15] Oliveira S. (1993) Parallel Multilevel Methods For Transport Equations, PhD thesis, Mathematics Department, University of Colorado at Denver, Denver, Colorado. [16] Oliveira S. (1994) Parallel Multilevel Algorithms for Anisotropic Transport Equations. Proceedings CTAC93 conference, Publ. World Scientific, Singapore, 388396. [17] Oliveira S. (1995) Multigrid and Krylov subspace methods for transport equations: aborption case, Proceedings of the 7th Copper Mountain Conference on Multigrid Methods, NASA Conference Publications. [18] Oliveira S. (1995) Krylov subspace methods for Transport Equations. Proceedings XENFIR III, Aguas de Lind6ia, Brazil.

173

74

S. Oliveira and Y. Deng

[19]

Oliveira S. (1996) Parallel Multigrid Methods for Transport Equations: Anisotropic Case. to appear in Parallel Computing.

[20]

Paige C. C. and Saunders M. A. (1982) LSQR: an algorithm for sparse linear equations and sparse least squares. ACM Trans. Math. Softw. 8, 43.

[21]

Patton W. and Holloway J. P. (1996) Application of Krylov subspace iterative methods to the slab geometry transport equation. Advancements and Applications in Radiation Protection and Shielding, Proceedings of the ANS Topical Meeting on Radiation Protection and Shielding, North Falmouth, MA, 384.

[22]

Saad Y. and Schultz M. (1986) GMRES: A generalized minimal residual algorithm for solving nonmetric linear systems. SIAM J. Sci. Statist. Comput. 7, 856.

[23]

Sonneveld P. (1989) CGS, a fast Lanczos-type solver for nonsymmetric linear systems. SIAM J. Sci. Statist. Comput. 10, 36.

[24]

Stewart D. E. and Leyk Z. (1994) Meschach: Matrix Computations in C. Proceedings of the Centre for Mathematical Sciences 32, Australian National University, Canberra.

[25]

Van der Vorst H. A. (1992) Bi-CGSTAB, a fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems. SIAM J. Sci. Statist. Comput. 13, 631.