Graziani F. (editor) Computational Methods in Transport

Подождите немного. Документ загружается.

Parallel Solution of the Time-Dependent S

Equations 481

The diﬀering results between ﬁgures clearly illustrate one important

point that must be made when comparing iterative methods for the solu-

tion of the S

equations: deﬁnitive statements about which method is bet-

ter can only be made about speciﬁc source-iteration/accelerator and Krylov-

algorithm/preconditioner combinations applied to a speciﬁc problem. For a

given problem it may always be possible to ﬁnd a better accelerator or pre-

conditioner that shifts the balance in favor of one iterative approach over

another.

One ﬁnal note we wish to make regarding the lessons oﬀered by this

example problem. In this simple case the best choice for source iteration

method would be to arrange the order of the forward and back substitution

sweeps associated with the direct solve so that the so that the left-to-right

sweep was ﬁrst followed by the right-to-left sweep. The calculation of the

reﬂected intensity, which will determine the boundary value for I for ordinates

where µ

< 0, ideally should be carried out after the left-to-right sweep is

complete and before the right-to-left sweep begins. This allows the reﬂected

boundary condition to be consistent with the values of I for µ

> 0thatwere

computed in the left-to-right sweep. However, in the case of two reﬂective

boundaries such an optimal ordering for the sweeps does not exist and the

reﬂective boundary condition at one end cannot be computed self-consistently

with the current value of the intensity I. This lack of self-consistency will

normally slow or halt the iterative convergence of the source iteration method.

5 Parallel Implementation

of the Full Linear System Approach

Since the full linear system solution approach to solving the S

equation

is relies only on a Krylov subspace algorithm to solve a linear system, the

parallel implementation of this method, in some sense, reduces to the par-

allel implementation of the Krylov subspace algorithm. However, one of the

desiderata for a scalable solution method was that it be compatible with hy-

drodynamic spatial domain decomposition schemes. For this reason we will

not consider other domain decomposition strategies that could be used for

pure transport simulations.

5.1 Parallelization of the Bi-CGSTAB Algorithm

Unlike the direct solves associated with the source iteration method most non-

stationary Krylov subspace algorithms do not make use of operations that are

inherently sequential. In fact the operations required for the parallelization

of the Bi-CGSTAB algorithm are relatively few.

• A matrix-vector multiply operation. Under spatial domain decomposition

schemes this requires only nearest-neighbor communication.

482 F.D. Swesty

• Vector addition and scalar-vector multiplication. Under spatial domain

decomposition both of these operations are embarrassingly parallel.

• Preconditioner creation and solves. For simple SPAI preconditioners this

requires only nearest neighbor communication.

• Vector inner products. Under spatial domain decomposition this operation

requires a global summation across processors.

In order to illustrate the simplicity of this method we write down the

steps of the Bi-CGSTAB Algorithm:

1. Provide initial guess: x

(0)

2. Calculate initial residual: r

= b − Ax

(0)

3. Set

r = r

4. For i =1, 2,...

5. 

i−1

(i−1)

6. if 

i−1

= 0 the method has failed

7. if i =1

8. p

(i)

= r

(i−1)

9. else

10. β

i−1

=(

i−1

/

i−2

)(α

i−1

/ω

i−1

)

11. p

(i)

= r

(i−1)

+ β

i−1

(i−1)

− ω

i−1

(i−1)

)

12. endif

13. precondition by solving M

p = p

(i)

14. v

(i)

= A

15. α

= 

i−1

(i)

16. s = r

(i−1)

− α

(i)

17. if s <tolerance

18. x

(i)

= x

(i−1)

+ α

19. stop

20. endif

21. precondition by solving M

s = s

22. t = A

23. ω

= t

s/t

24. x

(i)

= x

(i−1)

+ α

p + ω

25. r

(i)

= s − ω

26. if r

(i)

 <tolerance

27. stop

28. endif

29. endfor

The algorithm above can easily be implemented in less than 100 lines of F95

or C++ if BLAS subroutines or the equivalent are utilized for the basic vector

operations.

Implementing this algorithm under logically Cartesian spatial domain de-

composition requires that elements corresponding to zones at the edge of

a subdomain (hereafter referred to as ghost zones) be exchanged between

Parallel Solution of the Time-Dependent S

Equations 483

processors at certain stages in each iteration. Immediately prior to step 2 the

ghost zones of x

(0)

must be exchanged in order to evaluate Ax

(0)

. Similarly

between steps 13 and 14 ghost zones of p must be exchanged and between

steps 21 and 22 ghost zones of

s must be exchanged prior to performing the

matrix-vector multiplies in steps 14 and 22. If a simple SPAI preconditioner

is used in steps 13 and 21 then ghost zones exchanges are needed for the

vectors p

(i)

immediately prior to step 13 and s immediately prior to step 21.

These ghost zone exchanges can be accomplished via calls to MPI send and

receive routines.

Another parallel operation that is required is a global summation for

vector inner products in steps 5, 15, and 23. These three steps require a total

of four inner products. However, steps 17 and 26 require the evaluation of

the norm of residual vectors which brings the total number of inner products

to six per iteration. The global summations required for these inner products

can be accomplished via MPI”7016GLOBAL”7016REDUCE calls.

The overall ﬂoating-point operation time of the algorithm is usually dom-

inated by the matrix-vector multiply and/or the preconditioning steps. How-

ever, the message passing time associated with the ghost-zone exchanges and

the global reduction operations may be the prohibitive factor if the number of

ﬂoating-point operations per processor in the matrix-vector multiply is small.

Ultimately the diminishing ratio of computation time to communication time

as one increases the number of processes in a parallel simulation of ﬁxed size

will degrade the scalability.

6 Parallel Scalability of a 2-D Test Problem

To examine the parallel scalability of the full linear system solution approach

we consider the solution of a standard crooked pipe, a.k.a. “Top Hat” prob-

lem [GE01] using the simple corner balance spatial discretization scheme

(see [PA92] for details). This problem involves radiation traveling down a

optically translucent duct that makes several right angle turns. The material

surrounding the duct is optically thick. Heating of the material is modeled

using a standard heating opacity with a model absorption opacity and a

blackbody emissivity. This problem involves very long timescales as radia-

tion does not propagate around a turn in the duct until the material at the

end of a leg is heated and photons are re-radiated in new directions.

We zone this problem into 100 × 350 uniform spatial zones with 10 en-

ergy groups and an S

angular quadrature with 48 cosine pairs. Since the

simple corner balance scheme in 2-D discretizes each spatial zone into four

corners with an S

equation to be solved in each corner the total number of

equations to be solved is 100 × 350 × 10 × 48 × 4 ≈ 6.7 × 10

. The problem

is domain decomposed into a 2-D Cartesian process topology with roughly

equal numbers of zones per process to optimize parallel load balancing. For

a ﬁxed number of processors the number of processes in each dimension are

484 F.D. Swesty

chosen to make process topology as square as possible. The code, V2D, that

0 200 400 600 800 1000

Speedup

Number of Processors

0 200 400 600 800 1000

Speedup

Number of Processors

Fig. 7. Parallel Speedup relative to 32 processors for the 100 × 350 SCB crooked

pipe, a.k.a. “Top Hat” problem

implements the solver is written in F95 and uses MPI for all message passing.

The parallel speedup relative to 32 processors for this ﬁxed-size problem ob-

tained on the National Energy Research Supercomputing Center (NERSC)

IBM-SP is shown in Fig. 7. The for this problem problem size algorithm

scales quite well out to to 256 processors. As the number of processors is

increased to 512 and 1024 the speedup levels oﬀ as the ratio of communica-

tion time to computation time increases. For example, at 1024 processors we

have employed 32 ×32 process topology which means each process in the has

roughly a 3 ×11 interior mesh surrounded by a roughly 32 ghost zones. Thus

the ratio of ghost zones to interior mesh points is approximately unity. Of

course, a larger problem could easily maintain reasonable parallel eﬃciency to

larger numbers of processors but eventually communication costs will doom

the scalability of any ﬁxed-size problem.

7 Conclusions and Future Directions

The application of nonstationary Krylov subspace algorithms to time-

dependent discrete ordinates transport problems provides a successful av-

enue for implementing such problems on massively parallel architectures.

This method does not suﬀer from the drawback of having spatially recur-

sive sweeps which provide barriers to eﬀective parallelization in many cases.

We have also found a set of parallelizable preconditioners that have proven ef-

fective in aiding convergence of Krylov subspace iteration for time-dependent

equations. These preconditioners exploit the diagonal dominance of the

coeﬃcient matrix that is often present in time dependent problems. The an-

swer to the question of whether source iteration is more eﬀective than the

Parallel Solution of the Time-Dependent S

Equations 485

full linear system solution via Krylov subspace algorithms is highly depen-

dent on both the problem and the accelerators or preconditioners employed.

It is diﬃcult to make general statements regarding the relative merits of one

approach versus another.

There are a number of directions for future research in this area. The

possibility of hiding communication by utilizing asynchronous message pass-

ing to exchange ghost zones while computation is being carried out on the

interior of each subdomain can increase scaling in many circumstances. Sim-

ilarly, research into the subject of reducing the number of global reduction

operations by restructuring Krylov subspace algorithms could yield scalabil-

ity improvements. Finally, there is a continuing need for the development of

more eﬀective parallel preconditioners for the S

equations.

8 Acknowledgments

The author would like to acknowledge ﬁnancial support from Lawrence Liv-

ermore National Laboratory in the form of a ASCI contract and computing

support from the National Energy Research Scientiﬁc Computing Center.

References

[AL02] Adams, M.L. and Larsen, E.W.: Fast Iterative Methods for Discrete-

Ordinates Particle Transport Calculations. Prog. in Nucl. Eng. 40, 3–159

(2002)

[AL77] Alcouﬀe, R.E.: Diﬀusion Synthetic Acceleration Methods for the

Diamond-Diﬀerenced Discrete-Ordinates Equations. Nucl. Sci. Eng. 64,

344 (1977)

[AX96] Axelsson, O.: Iterative Solution Methods. Cambridge University Press,

Cambridge (1996)

[BA94] Barrett, R. and et al.: Templates for the Solution of LInear Systems:

Building Blocks for Iterative Methods SIAM, Philadelphia (1994)

[BE02] Benzi, M.: Preconditioning Techniques for Large Linear Systems: A Sur-

vey. J. Comp. Phys. 182, 418–477 (2002)

[DM05] D’azevedo, E.F. and Messer, B. and Mezzacappa, A. and Liebendorfer,

M.: An ADI-Like Preconditioner for Boltzmann Transport. SIAM Sci.

Comput., 26, 810–820 (2005)

[GE01] Gentile, N.: Implicit Monte Carlo Diﬀusion–An Acceleration Method for

Monte Carlo Time-Dependent Radiative Transfer Simulations. J. Com-

put. Physics 29, 543–571 (2001)

[KBA92] Kock, K.R. and Baker, R.S. and Alcouﬀe, R.E.: Solution of the First-

Order Form of the 3-D Discrete Oridnates Equation on a Massively Par-

allel Processor. Trans. Amer. Nucl. Soc., 65, 198 (1992)

[LA82] Larsen, E.W.: Unconditionally Stable Diﬀusion Synthetic Acceleration

Methods for the Slab Geometry Discrete Ordinates Eqautions. Part 1:

Theory. Nucl. Sci. Eng. 82, 47 (1982)

486 F.D. Swesty

[LM93] Lewis, E.E. and Miller, W.F.: Computational Methods of Neutron Trans-

port. American Nuclear Society Inc, La Grange Park (1993)

[MB93] Mezzacappa, A. and Bruenn, S.W.: A numerical method for solving the

neutrino Boltzmann equation coupled to spherically symmetric stellar

core collapse. Astrophys. J., 405, 669–684 (1993)

[MM99] Mihalas, D, and Weibel-Mihalas, B.: Foundations of Radiation-

Hydrodynamics. Dover, Mineola (1999)

[MO82] Morel, J.E.: A Synthetic Acceleration Method for Discrete Ordinates

Calculations with Highly Anisotropic Scattering. Nucl. Sci. Eng. 82, 34–

46 (1982)

[PA92] Palmer, T.S.: Curvilinear Geometry Transport Discretizations in Thick

Diﬀusive Regions. Ph.D. Thesis, University of Michigan, Ann Arbor

(1992)

[PH02] Patton, B.W., and Holloway, J.P.: Application of Preconditioned GMRES

to the Numerical Solution of the Neutron Transport Equation. Annals

Nucl. Eng., 28, pp. 109–136 (28)

[SA03] Saad, Y.: Iterative Methods for Sparse Linear Systems, 2nd ed. SIAM,

Philadelphia (2003)

[SS86] Saad, Y. and Schultz, M.: GMRES: A Generalized Minimal Residual Al-

gorithm for Solving Nonsymmetric Linear Systems. SIAM J, Sci. Statist.

Comput. 7, 856–869 (1985)

[SS04] Swesty, F.D. and Smolarski, D. C. and Saylor, P.E.: A Comparison of

Algorithms for the Eﬃcient Solution of the Linear Systems Arising from

Multi-Group Flux-Limited Diﬀusion Problems. Astrophys. J. Suppl. 153,

369–387 (2004)

[SM05] Smolarski, D. C. Diagonally Striped Matrices and Approximate Inverse

Preconditioners. J. Comput. and Appl. Math., in press, 2005

[VO92] van der Vorst, H.: BiCGSTAB: A Fast and Smoothly Converging Variant

of B-CG for the Solution of Nonsymmetric Linear Systems. SIAM J, Sci.

Statist. Comput., 13, 631–644 (1992)

Diﬀerent Algorithms of 2D Transport

Equation Parallelization

on Random Non-Orthogonal Grids

Shagaliev R.M., Alekseev A.V., Beliakov I.M., Gichuk A.V., Nuzhdin A.A.,

Rezchikov V.Yu.

Numerical solution of a non-stationary kinetic equation at a multi-group

setting in the 2D spatial approximation demands much computer operat-

ing memory and calendar execution time. One of the ways of solving such

problems is using algorithms of parallelization.

Solving the problem is particularly diﬃcult in the general case, when the

boundary problem for the transport equation is put in the space-domains with

complex shapes, and this equation is approximated on non-orthogonal spatial

grids. As a rule, implicit time-dependant schemes are used to numerically

solve these tasks, while the system of grid equations on the steps by time is

solved with the method of sweep computation. This means that at this time

step grid transport equations are solved in strict succession, these successions

are diﬀerent for diﬀerent directions of the particle motion and change at time

steps. In other words, when the transport equation is being solved in the

kinetic approximation at multi-processor systems, one has no possibility to

predetermine the linkage layout between the processor elements (order of

data exchange), which crucially complicates the problem of development of

eﬃcient algorithms of transport problem parallelization.

Several algorithms of parallelization a 2D non-stationary kinetic trans-

port equation at a multi-group setting on non-orthogonal spatial grids are

compared in the present Report.

The ﬁrst algorithm is a large-block parallelization algorithm. When such

an algorithm was used one of the space-domains of investigated system was

calculated at one of the processors. (As a space-domain we take a family of

system points, joined by similar physical properties. In each space-domain

the kinetic transport equation was solved independently). To carry out the

computation correctly iterations by inner boundary conditions were applied,

each iteration was followed by inter-processor exchange of inner boundary

conditions. The linkage between the multi-processor computer units was held

by means of a special library of MPI data transport. One space-domain, as

well as arbitrary number of space-domains, was successfully calculated on the

processors at a time. All the applied transports are asynchronous, the points

of transport asynchronization were at the end of the iteration cycle before

the transport of inner boundary conditions.

488 R.M. Shagaliev et al.

The eﬃciency of this algorithm during solving a test problem was 78%.

The test parameters are presented in Table 1 below.

It is evident that this parallelization algorithm has considerable restric-

tions concerned, ﬁrst of all, with the dimensions of the computed system

and its fragmentation. The problem of unbalanced computation load on the

processors for the reasons of diﬀerent number of points in the space-domains

in the general case is acute for practically every computation.

Table 1. Test problem parameters

Quantity of Quantity of

Quantity of Space- Space-Domain Row Quantity of Computed

Groups Domains Number Number Columns Cells

1 28 20 560

7 3 2 16 35 560

3 24 23 552

An other algorithm, which in its turn provides possibility for balanced

load of the processors, is the algorithm of minor block parallelization, where

the principle of space decomposition of the initial system into sub-domains

was applied. System fragmentation into sub-domains, which were then dis-

posed at separate processors (further on “para-domains”) was carried out

here only in one direction, viz along the columns of the spatial grid. This

parallelization algorithm were somehow similar to the event dynamic system.

At each parallelization cycle the boundary conditions for the covered bound-

aries of para-domains, provided from adjacent processors, were analyzed, af-

ter which, results dependant, the further behavior scenario was speciﬁed. The

main events included in this dynamically speciﬁed scenario were initialization

of asynchronous exchanges, their completion, processing the provided bound-

ary conditions, computation of a para-domain, computation of the resolution

ratio for a family of cells.

The basic principles of the algorithm of minor block parallelization with

decomposition along the columns were as follows:

• Spatial decomposition of the space-domain into para-domains; the space-

domain was split into columns of the spatial grid;

• The system of grid equations was solved independently by the angle vari-

able;

• Complete integration of the inter-processor exchange with the computa-

tion, which was achieved by the application of asynchronous receive/send

operations.

Numerical investigations of this parallelization algorithm were carried out on

a one-space-domain spherically symmetrical model problem. The thermody-

Diﬀerent Algorithms of 2D Transport 489

namic and optical properties of the medium were analogous to those in the

Flack’s task.

The problem space-domain was a hemi-sphere with the radius 1 cm. This

space-domain was covered with the spatial grid comprised of 3×Np columns

(Np – number of processors) and 50 rows.

The incoming energy ﬂow on the outer surface equaled zero, the initial

radiation density equaled zero as well.

Fore the ﬁrst computation series the angular grid with 126 directions

by the space angle was selected.

The computation was carried out with the seven-group approximation.

The problem was solved in 10 steps by time. The right-hand-side iterations

were carried out along with the solution of a correction equation with the

KM-method and completed with the error 0.001.

The numbers of processors during the computation were: 1, 10, 20, 30,

40, 50.

The main parameters of the parallel algorithm eﬃciency, namely, the

speedup factor and the parallelization eﬃciency, are presented in Table 2

and in Figs. 1, 2.

Table 2. Speedup and parallelization eﬃciency in the ﬁrst computation series

Number of Processors Speedup Parallelization Eﬃciency

1 1 100

10 10 100.5

20 20.2 101

30 27.55 91.85

40 33.67 84.18

50 39.67 79.34

The second computation series was carried out on the angular grid

with 160 directions by the space angle.

The computation was carried out with the twenty-eight-group approxi-

mation.

The problem was also solved in 10 steps by time. The right-hand-side

iterations were carried out along with the solution of a correction equation

with the

I-method and completed with the error 0.001.

The numbers of processors during the computation were: 1, 10, 20, 30,

40, 50. The main parameters of the parallel algorithm eﬃciency, namely,

the speedup and the parallelization eﬃciency, are presented in Table 3 and

Figs. 3, 4.

The requirements of higher precision of the solution of multi-group trans-

port problems at numerical simulation make important the use of big num-

ber of intervals by an energy variable (groups). The third parallelization

490 R.M. Shagaliev et al.

0 102030405060

Computation

Ideal

Fig. 1. Speedup dependence on the number of processors

11020304050

100

100.5

101

91.85

84.18

79.34

100

120

Fig. 2. Dependence of the parallelization eﬃciency on the number of processors