A Nested FGMRES Method for Parallel Calculation of Nuclear Reactor Transients.
A semi-iterative method based on a nested application of Flexible Generalized Minimum Residual)FGMRES) was developed to solve the linear systems resulting from the application of the discretized two-phase hydrodynamics equations to nuclear reactor transient problems. The complex three-dimensional reactor problem is decomposed into simpler, more manageable problems which are then recombined sequentially by GMRES algorithms. Mathematically, the method consists of using an inner level GMRES to solve the preconditioner equation for an outer level GMRES. Applications were performed on practical, three-dimensional models of operating Pressurized Water Reactors (PWR). Serial and parallel applications were performed for a reactor model with two different details in the core representation. When appropriately tight convergence was enforced at each GMRES level, the results of the semi-iterative solver were in agreement with existing direct solution methods. For the larger model tested, the serial performance of GMRES was about a factor of 3 better than the direct solver and the parallel speedups were about 4 using 13 processors of the INTEL Paragon. Thus, for the larger problem over an order of magnitude reduction in the execution time was achieved indicating that the use of semi-iterative solvers and parallel computing can considerably reduce the computational load for practical PWR transient calculations.
Introduction
. The analysis of nuclear reactor transient behavior has always
been one of the most difficult computational problems in nuclear engineering. Because
the computational load to calculate detailed three-dimensional solutions of the field
equations is prohibitive, variations of the power, flow, and temperature distributions
are treated approximately in the reactor calculation resulting in considerable conservatism
in reactor operation. Some researchers have estimated that using existing
methods, the computational load to calculate three-dimensional distributions would
exceed a Teraflop per time step. [16] [6]
Over the last several years computer speed and memory have increased dramatically
and has motivated a rethinking of the limitations of the existing reactor transient
analysis codes. Researchers have begun to adapt three-dimensional hydrodynamics
and neutron kinetics codes to advanced computer architectures and have begun to investigate
advanced numerical methods which can take full advantage of the potential
of high performance computing.
The overall goal of the work reported here was to reduce the computational burden
for three-dimensional reactor core models, thereby enabling high fidelity reactor
system modeling. The specific objective was to investigate Krylov Sub-Space methods
for the parallel solution of the linear systems resulting from the reactor hydrodynamics
equations.
The following section provides a brief description of the hydrodynamic model and
reactor problem used in the work here. The nested GMRES method and preconditioner
developed in the work here are described in section 3 and serial and parallel
applications are presented in sections 4 and 5, respectively.
This work was supported by the Electric Power Research Institute
y School of Nuclear Eng., Purdue University, W. Lafayette, IN 47907
2. Hydrodynamic Model and Reactor Problem. The nuclear reactor analysis
problem involves the solution of the coupled neutron kinetics, heat conduction,
and hydrodynamics equations. Of these, two-phase flow hydrodynamics is generally
the most computationally demanding and was the focus of the work here.
The hydrodynamic method used here is consistent with the reactor systems code
RETRAN-03 [7] which is widely used in the nuclear industry for the analysis of
reactor transient behavior. The method is based on a semi-implicit solution of the
mass and momentum equations for each phase and the energy equation for the fluid
mixture. The solution scheme uses finite difference representations for all of the
fluid-flow balance equations. All convective quantities and most source terms are
linearized, i.e. expanded using a first order Taylor series. As a result, the coupling
between the finite difference equations is implicitly included in the system of coupled
equations. The method is considered only semi-implicit since linearized equations
are used rather than the original nonlinear partial differential equations. A standard
Newton-Raphson technique is used to solve the nonlinear equations.
2.1. Hydrodynamic Model. The application of spatial finite differencing to
the set of governing partial differential equations is performed using the concept of
volumes and connecting junctions. This results in a system of ordinary differential-
difference equations which may be expressed as
dY
dt
(1)
where Y is a column vector of nodal variables and F is a column vector of functions
is the number of dependent/solution variables. The solution
vector, Y consists of N j junction mass flow rates (W ) and slip velocities (V SL ), and
volume total mass (M ), total energy (U ), and vapor mass (M g ) inventories.
(2)
The vector F(Y) is linearized and a first order time difference approximation is
used, resulting in the following linear system:
where I is the identity matrix and J is the matrix Jacobian. \DeltaY
and \Deltat are the values of Y at time levels t n+1
and t n respectively.
Because of the semi-implicit nature of the formulation, the linear system that
arises after linearizing and discretizing is tightly coupled. Also, the high degree of
stability of this formulation imposes a less stringent limit on the size of the time step.
As indicated in Equation 3, larger time steps further reduce diagonal dominance and
results in a more ill-conditioned linear system. Several authors have noted [5] [3] that
such ill-conditioned linear systems with no well-defined structure are very difficult to
solve efficiently and in many cases special handling is required to obtain an acceptable
solution.
The present linear solver in RETRAN-03 is a direct method based on Gaussian
elimination which has a \Theta(N 3 ) execution time complexity where N is the number of
unknowns. For models dominated by one-dimensional flow, direct solutions complemented
by some type of matrix reduction technique can be very efficient. However,
in the case of a high fidelity model of a reactor with a three dimensional core rep-
resentation, direct methods become inefficient. This can be demonstrated using the
model of a standard 4-loop Presurized Water Reactor.
A
A
A
A
A
A
A
AAA
AAA
AAA
AAA
AAA
AAA
AAA
A
A
A
A
A
A
AAA
AAA
AAA
AAA
AAA
AAA
A
A
A
A
A
A
A
AAA
AAA
AAA
AAA
AAA
AAA
AAA
A
A
A
A
A
A
AAA
AAA
AAA
AAA
AAA
AAA
System Model
Model
see
Model
below
Side View Top View
channel
Upper Plenum
26 33
28 14 31 17
43 23
48
36 343111 355769 7177737579656 4
by pass
Fig. 1. Nodalization for the PWR Model
2.2. Pressurized Water Reactor Problem. The reactor system model used
in this work consists of 4-loop pressurized water reactor model with a detailed three
dimensional representation of the reactor core. A schematic is shown in Figure 1.
Each loop of the system contains a steam generator in which heat is transferred from
the pressurized primary loop to a secondary loop containing the steam turbines. The
core model consists of several volumes stacked one upon the other along each channel.
Assuming a layout of C by C channels and A axial volumes, the number of volumes
is AC 2 , the number of cross flow junctions is and the number of vertical
junctions is C 2 (A+1). The first problem used in the work here had a core model with
9 channels and 12 axial levels per channel as shown in the Figure. We constructed
two versions of the core model, one with cross flow (horizontal flow) between volumes
and the other without cross flow. Cross flow is important for reactor transients such
as the break of a steam line in one of the loops which results in the horizontal mixing
of hot and cold water in the core.
The model with cross flow consists of 173 volumes and 339 junctions which results
in a linear system of size 1173. The core of the model with cross flow contains 108
volumes and 261 junctions and forms nearly 70% of the system. The model without
cross flow contains the same number of volumes but has only 183 junctions and results
in a linear system of size 885. The coefficient matrix from the problem with cross flow
is shown schematically in Figure 2. The structure shown in the figure corresponds to
a sequential ordering of the junction and volume variables which is more convenient
for the reduction/elimination method currently used in the code. A reordering that
is more suitable for preconditioning will be discussed in the following section.
Columns
Rows
Fig. 2. PWR Model with 3-D Core: Sparsity Pattern of the Coefficient Matrix
For an initial test problem a simple core reactivity event was modeled by simulating
the insertion and withdrawal of a control rod for 20 seconds. The rod was
inserted into the core at a rate of 0.08 dollars of reactivity per second for 5 seconds
and then immediately withdrawn at the same rate. The computational performance
for this problem on a single processor of the INTEL Paragon is summarized in Table
1 for the models with and without cross flow. The direct linear system solution is
performed in the SOLVE module indicated in the Table.

Table
Computatinal Summary for PWR Example Problem with Direct Linear Solver
Module w/o Cross Flow w/ Cross Flow
CPU Time Percent CPU Time Percent
OTHER 155.89 46.7 196.17 4.9
TOTAL 333.88 100.0 4931.71 100.0
It is immediately apparent from Table 1 that the use of cross flow in the core
increases the computational time by an order of magnitude. This is primarily because
the problem without cross flow results in a linear system containing predominantly
tridiagonal submatrices which lends itself very efficiently to reduction/elimination
methods. This large increase in the execution time discourages the modeling of cross
flow and, in general, high fidelity reactor simulation.
Although a reduction in the operation count required for a direct solution was the
primary motivation for the work here, there were other reasons to consider iterative
linear solution methods. First, the actual problem being solved in RETRAN-03 is
nonlinear and the direct solution of the resultant linearized equations can be a waste
of floating point operations since usually an accuracy considerably less than machine
precision is adequate. Secondly, unlike direct methods, semi-iterative solution methods
can be accelerated with information from the previous time steps, such as an
initial guess or a preconditioner. Finally, direct methods do not lend themselves easily
to parallel computing on distributed memory, MIMD multicomputers, whereas,
many of the new semi-iterative linear solvers can be efficiently parallelized.
3. Semi-Iterative Linear Solvers for 3-D Hydrodynamics. During the last
several years, considerable research has been performed on a non-stationary class of
techniques, collectively known as Krylov subspace methods. These include the classical
Conjugate Gradient (CG) method which has been shown to be very efficient for
symmetric positive definite systems of equations. These methods are called Krylov
methods because they are based on building a solution vector from the Krylov sub-
space: span fr is the residual of the initial solution and
A is the coefficient matrix. The coefficients of the solution vector in the case of the
CG method are based on the minimization of the energy norm of the error. In gen-
eral, the linear systems encountered in hydrodynamics problems are not symmetric
positive definite and therefore can not be solved using the CG method.
Numerous Krylov methods for solving the non-symmetric problem have been
proposed over the years and several were considered in this work to include the Generalized
Minimal Residual(GMRES)[10] method, the BiConjugate Gradient (BiCG)
method [1], the Conjugate Gradient Squared (CGS) method [12], and the BiConjugate
Gradient Stabilized (BiCGS) method [15]. Of particular interest in the work
here is the GMRES method which solves a minimization problem at each iteration
and therefore guarantees a monotonically decreasing residual.
It is well known that the convergence rate of the Krylov methods depends on
the spectral properties of the coefficient matrix and that proper preconditioning can
considerably improve the rate of convergence. Because preconditioning involves some
additional cost, both initially and per iteration, there is a trade-off between the cost
of implementing the preconditioner and the gain in convergence speed. Since many
of the traditional preconditioners have a large sequential component, there is further
trade-off between the serial performance of a preconditioner and its parallel efficiency.
Several alternative preconditioners were examined in the work here.
3.1. Application of Krylov Methods to Reactor Problem. Previous work
on the application of semi-iterative solvers to reactor core hydrodynamic calculations
[13] [8] has focused primarily on the solution of the linearized pressure equation resulting
from single phase, one-dimensional flow problems. Some of the Krylov methods
were found to perform very well on the resulting tridiagonal system of equations.
In particular, Turner achieved excellent convergence with the Conjugate Gradient
Squared method and an ILU preconditioner. While this work provided useful insight,
the linear systems resulting from the two-phase, three-dimensional flow problems are
significantly different and the performance of the various Krylov methods was very
different.

Figure

3 shows the behavior of the absolute residual (L2 norm) for the application
to the PWR problem of various Krylov methods with no preconditioner. The
performance of Bi-CGSTAB, BICG, and CGS are very irregular (CGS is not shown
because the residual increased by several orders of magnitude in the first few iterations
and then continues to increase). Only the GMRES method demonstrated acceptable
behavior with a monotonically decreasingly residual which is expected because of its
basis in a minimization principle. Linear systems from other time steps were examined
and behavior similar to that shown in Figure 3 was observed. Although the floating
point operation count per iteration was higher for GMRES, particularly for larger
numbers of iterations, it was attractive for applications here because of its inherent
robustness. Preconditioning techniques were then examined that could improve the
convergence behavior of GMRES.
3.2. Domain Decomposition Preconditioning. Knowledge of the physical
characteristics of the system is invaluable in choosing a good preconditioner. It is
evident from Figure 1 that the core and the ex-core systems interact only at the
upper and lower plenum. This suggests a natural way to decompose the problem
domain. If the system were to be reordered such that all the solution variables that
belonging to the core are placed contiguously, then the structure of the resultant
matrix is given by
AC U
the blocks L and U represent the interactions between the core and the ex-core vari-
ables. Because the core and the ex-core interact only at the upper and lower plenum,
these matrices have very few non zero elements. Several options exist for using Equation
4 as a preconditioner.
One possible preconditioner is to neglect the L and U blocks entirely, resulting in
a block Jacobi preconditioner given by
AC
GMRES
Iterations
Residual
Fig. 3.

Figure

Performance of Krylov Solvers on PWR 3-D Core Problem
The block AE represents the interactions between variables in the ex-core region
and the block AC represents the interactions between variables in the core region.
Because the preconditioner is Jacobi, the two blocks may be solved individually which
is attractive for parallel computing. As noted earlier, the size of the core block is
considerably larger than the ex-core block and in a later section methods are discussed
for solving the core problem.
The block Jacobi domain decomposition preconditioner was applied to the 3-D
PWR example problem first using a direct solution of the preconditioner equation.

Table

2 contains the results for the linear systems at the first two time steps of the
transient described in section 2.2. The first case corresponds to the problem for which
results were shown in Figure 3. In addition to the number of iterations, two other
measures of the effectiveness of the preconditioner are shown in the Table. First,
the condition of the original and preconditioned matrix is shown in the second and
third columns. In the fourth column is shown another measure of effectiveness of
preconditioning suggested by Dutto[4] which uses the ratio of the Frobenius norm
of the remainder matrix, defined as , to the Frobenius norm of original
matrix. The Frobenius norm is defined as
diagonal
As shown in the Table, both measures indicate the Jacobi preconditioner is very
effective. For a tolerance of 10 \Gamma6 (relative residual), convergence is achieved in less

Table
Results of Application of Domain Decomposition to PWR Problem
Preconditioner
kAkF Iterations
2:19
Gauss Seidel 8:30
2:19
Modified
Gauss Seidel 2:19
In the block Jacobi preconditioner the impact of the submatrices L and U is
completely neglected while solving the preconditioner equation. Another possibility
is to employ a Gauss-Seidel like technique in which the preconditioner equations are
solved sequentially:
z C
z C and -
z C and z E respectively. Two approaches were
examined for solving these equations as a preconditioner.
In one approach, the first equation is solved neglecting the coupling to the other
(i.e. -
z C or -
Because of symmetry, the sequence of solution does not matter
and the preconditioner has the form:
AC
Results of applying this technique are shown as the Gauss-Seidel preconditioner
in

Table

2. As indicated there is only minor change in the measures of effectiveness
of the preconditioning, however, convergence is achieved in fewer iterations.
In the second approach an estimate of the ex-core solution, -
z E , is formed using
the decoupled ex-core equation: AE -
and is then used to solve sequentially the
core and ex-core equations given by Equation 7. The preconditioner for this approach
can be expressed as:
This approach is shown in the Table 2 as the Modified Gauss-Siedel preconditioner
and, as indicated, little differences are observed in the effectiveness of the precondi-
tioning. Because the Block Jacobi method is naturally parallel, it was used in the
work here even though the Gauss-Siedel methods showed slightly better numerical
performance.
For the results shown above the preconditioner equation was solved directly using
Gaussian Elimination. Because the ex-core problem is predominantly one-dimensional
flow, a direct solution using reduction/elimination proves to be very efficient. Con-
versely, the three-dimensional core model does not lend itself to direct solution and
the following section examines the use of a second level GMRES to solve the core
problem.
3.3. Preconditioning the Core Problem. The original matrix ordering shown
in

Figure

2 is not conducive to preconditioning the core problem, AC z
several alternate orderings were examined. Because the junctions are physically between
volumes, reordering the solution vector such that it bears some resemblance to
the physical layout would help in decreasing the profile of the matrix. The goal is
to increase the density of the matrix in the regions around the diagonal (i.e reduce
bandwidth) and then use that portion of the matrix for preconditioning purposes (e.g
Block Jacobi, etc.
The structure existing in the core was exploited by defining a supernode that
consisted of both volumes and junctions. The physical domain was discretized into
several supernodes that introduced homogeneity in the structure of the matrix. The
supernode could be considered a fractal for representing the smallest unit of structure
present in the system. As shown in Figure 4, the supernode for the core problem
consists of a volume and 3 junctions. One of the junctions is the vertical junction
upstream to the volume and the other two are crossflow junctions leading out of the
volume, each in a different direction. It should be noted that the use of supernodes
leads to a problem of extra junctions at some of the exterior channels. These junctions
are dummy junctions and are represented in the matrix but do not appear in the
solution. Hence the size of the system increases but its condition number remains the
same.
Volume
Vertical Junction
Crossflow Junction
Crossflow Junction
Fig. 4. Structure of a Supernode
Several orderings of the supernodes were considered (e.g. Channel-wise, Plane-
wise, Cuthill-McKee, etc) and ordering planewise was found to be the best for the
purpose of preconditioning. Each plane represents a two dimensional grid of the
supernodes, and the planes themselves form a one dimensional structure that are
linked to only two neighbors, resulting in a block tridiagonal matrix as shown in

Figure

5 and Equation 10. As in the case of the outer level preconditioner, several
options exist for using Equation 10 as a preconditioner for the inner level GMRES.
12002006001000Columns
Rows
Fig. 5. Structure of Coefficient Matrix : Plane Wise Ordering
A block diagonal preconditioner which neglects planar coupling was considered
first.

Table

3 shows the results of solving the second level GMRES during the first
outer iteration. Four different block sizes were examined with a convergence criterion
of 10 \Gamma6 . The two cases shown, A 1
C and A 2
C are taken from different time steps.

Table
Results of the Application of Jacobi Preconditioning on the Inner Level
Preconditioner
block size Condition # kRkF
kAkF
Iterations Condition # kRkF
kAkF
Iterations
The results for the other outer iterations were slightly different since the source
and hence the initial residual of the second level GMRES was different. However, the
results from other cases were consistent with the general trend shown in Table 3. As
expected, larger block sizes reduce the number of iterations. However, the cost of solving
each block directly would increase as N 3 which offsets the reduction in iterations.
Also, smaller blocks have the advantage of scalability for parallel computing.
A domain decomposition scheme incorporating the interactions between the diagonal
blocks was also examined. An approximate solution, -
z C , was computed by
solving the block Jacobi system which neglects the coupling between adjacent planes:
z 0
The prime notation is used to distinguish the z C and r C vectors here from those which
occur in the outer iteration, Equation 7. The inner GMRES preconditioner equation
is then solved:
where the preconditioner, MD , is:
As indicated in Table 4, this improved preconditioning does reduce the number
of iterations, but solving the preconditioner becomes more expensive. The number of
iterations for a block size of 81 is reduced by half but more than twice the number of
floating point operations are required to form MD .

Table
Results of the Application of Domain Decomposition Preconditioning on the Inner Level
kAkF Iterations
block size
Preconditioners other than Block Jacobi were investigated for the core problem.
Primarily because of the ill-conditioned nature of the coefficient matrix, popular
schemes such as SSOR and ILU were ineffective. For example, the incomplete LU
factorization scheme was tested on a linear system arising from the 3-D Core prob-
lem. Banded
were constructed and the Frobenius norm of the
ILU preconditioner was compared to the norm of the exact inverse:
As shown, the Frobenius norm of the approximate preconditioned system, M \Gamma1 A
is over three orders of magnitude larger than the Frobenius norm of A \Gamma1 A, which
suggests an ILU preconditioner would not be very effective.
3.4. Nested GMRES. The methods described in the previous sections were
implemented in a nested GMRES algorithm which consists of using an "inner" level
GMRES to solve the preconditioner equation for the "outer" level GMRES. Such a
strategy was suggested by Van der Vorst and demonstrated successfully for several
model problems [14]. The preconditioner for the inner level GMRES itself could
be solved using a third level GMRES, but in the applications here a direct solver
proved most efficient. A schematic of the nested GMRES algorithm is shown in

Figure

6. A more physical interpretation of the algorithm is to view the overall
problem as being decomposed into three simpler, more manageable problems which
are then recombined sequentially by GMRES algorithms. At the highest level we take
advantage of the naturally loose coupling between the core and the ex-core components
and solve separately the linear systems of the core and ex-core regions. These solutions
are then recombined using the highest or "outer" level GMRES. At the second or
"inner" level, GMRES is used to solve the 3-D core flow problem where focus is on
the coupling between the vertical flow channels in the core. And finally at the third
or lowest level, GMRES or a direct solver is used to restore coupling between nodes
in a plane.
CORE PROBLEM
SYSTEM PROBLEM
(GMRES I)
(GMRES II)
(GMRES III
or Direct)
A e z
A c z
Fig. 6. Nested GMRES Algorithm for RETRAN-03 Linear System Solution
3.5. Flexible General Minimum Residual (FGMRES). The solution of
the preconditioning equation in the GMRES method with another GMRES algorithm
poses a potential problem due to the finite precision of the inner level solution. In
the preconditioned GMRES algorithm the solution is first built in the preconditioned
subspace and then transformed into the solution space:
where the matrix is the set of orthonormal vectors. In the case of
an iterative solution for the preconditioner, the transformation is only approximate,
the extent of which is determined by the convergence criterion. In the case of ill-conditioned
matrices the approximations could be especially troublesome and a very
tight convergence is required to minimize error propagation.
The problem of the inexact transformation was alleviated in the work here by
using a slight variant of the GMRES algorithm in which an extra set of vectors is
stored and used to update the solution. This modification of the GMRES algorithm
was suggested by Saad and is called Flexible General Minimum Residual(FGMRES).
[9] This algorithm allows for the complete variation of the preconditioner from one
iteration to the next by storing the result of preconditioning each of the basis vectors,
while they are being used to further the Krylov subspace. Instead of using Equation
15, the final transformation to the solution subspace is then performed using .
where the matrix . The FGMRES algorithm is
given in Appendix A and was implemented in the nested GMRES method.
4. Serial Applications.
4.1. Static Problem. The nested GMRES algorithm was first applied to the
linear systems arising from the first few time steps of the rod withdrawal/insertion
transient. Parametrics were performed on the convergence criterion and number of
iterations for each of the levels. A maximum iteration limit was set on the inner
GMRES since in some cases the rate of convergence was very slow (sometimes termed
"critical slowing down").
The variation of the outer (highest) level residual during the iterations is shown
in

Figure

7 for different maximum number of inner iterations (miter). The results
indicate that the rate of decrease in the residual for
substantially greater than for miter = 20. However, as shown in Figure 8, the results
of subsequent timesteps indicate that the difference in the rate of decrease in the
residual between miter = gradually diminishes.
The residual achieved on the inner (second) level GMRES for the timestep corresponding
to Figure 8 are shown in Figure 9. For each of the maximum iteration
limits, the residual increases after the first few iterations. One possible explanation is
that the most desirable search directions for the outer level GMRES become harder
to resolve, leading to a degradation of the performance of the inner level. However,
the diminished quality of the preconditioning from the inner level GMRES does not
appear to have a deleterious effect on the convergence of the outer iteration.
Based on these convergence results and those from the analysis of other time
steps, a strategy was formulated for the implementation of the nested GMRES in the
transient analysis code RETRAN-03. The following section discusses the results of
applying GMRES to the transient problem.
Outer Iteration Number
Outer
Residual
Fig. 7. Reduction of Outer Residual: First Time Step
Outer Iteration Number
Outer
Residual
Fig. 8. Reduction of Outer Residual: Subsequent Time Step
Outer Iteration Number
Inner
Residual
Fig. 9. Performance of Inner GMRES: Subsequent Time Step
4.2. Transient Problem. In the time-dependent iterative solution of the hydrodynamics
equations, error from incomplete convergence at one time step can be
propagated into the coefficient matrix of subsequent time-steps. A preliminary assessment
of the effect of convergence criteria on the quality of the solution was performed
using a "null" transient in which the steady-state condition is continued for several
seconds with no disturbance to the system. Several outer level convergence criteria
were investigated and acceptable performance was achieved with a tolerance of 1.0E-
09 in the relative residual. At higher tolerances some minor deviation was observed
(e.g. less than 1% relative error) in performance parameters such as the core power
level. This provided initial guidance in setting tolerance and iteration limits for the
transient problems.
The PWR rod withdrawal/insertion transient described in section 2.2 was then
analyzed using RETRAN-03 and nested FGMRES. The model with cross flow in the
core was studied using the iterative solver for a transient time of 20 seconds which
required steps. The performance of the iterative algorithm was analyzed by
first varying the number of outer iterations using a direct solver for the inner (e.g.
problem, and then by varying both the number of outer and inner GMRES
iterations. The results are shown in Table 5 as CPU seconds per time step. The
maximum relative residual error during any of the outer iterations is shown in the
fourth column of the table.
For purposes of comparison, the RETRAN-03 solution from Table 1 which uses a
direct solution of the linear system is repeated as case A.1 in Table 5. In the first two
cases (A.2 and A.3) the inner two GMRES levels (see Figure 6) have been replaced
with a direct solver and GMRES is used only for the outer or highest level iteration.
The objective here was to isolate the impact of the number of outer iterations on the

Table
Serial Performance of RETRAN with GMRES: 9 Channel/ 12 Axial Case
Iterations Number of CPU secs/Time Step
Case Outer Inner max( krk
A.1 Direct Direct \Gamma\Gamma 318 14.89 15.65 -
A.3 12 Direct 4
A.4
A.6
quality of the solution. It can be seen from Table 5 that the maximum residual error
decreases as the number of outer iterations increases. However, in both cases A.2 and
A.3, no significant error was observed in the important physical parameters. When
the number of outer iterations was reduced to four, minor deviations began to occur
in the solution after 15 seconds of the transient. It should be noted that the execution
times for cases A.2 and A.3 are generally comparable to the direct solver.
The GMRES algorithm was then employed for both the outer and inner iterations,
keeping the direct solution for the innermost (third) level. In order to gain some insight
on the relation between convergence of the inner and outer GMRES algorithms, the
number of outer and inner iterations were varied as shown in Table 5 for cases A.4, A.5,
A.6 and A.7. The accuracy of these solutions was examined for important physical
parameters in the solution. The normalized core power and the pressurizer pressure
(volume 58 in Figure 1) are plotted versus time in Figure 10 for the Direct solution
(A.1) and for GMRES solutions (A.4 and A.7). Some minor deviation is observed
in the solution with iterations. The execution times are greater than the
direct (A.1) and GMRES outer with direct inner (A.3) solutions. However, as will be
discussed in the next section, the algorithm with inner level GMRES is more attractive
for larger problems and for parallel computing.
5. Parallel Applications. One of the attractive features of preconditioned Krylov
methods is their potential for parallel computing. The emphasis in this work was on
the use of a distributed memory parallel architecture and applications were performed
on the INTEL Paragon. This section describes the mapping of the nested GMRES
onto the Paragon and the execution time reductions achievable for the PWR sample
problem.
The most natural mapping of processors for the PWR model was one processor
to the ex-core and one to each of the 12 planes in the core model. The matrices were
striped row wise, implying that L and AE of Equation 4 were stored on the processor
that is assigned the ex-core. Ideally, AC and U should be partitioned among the 12
processors but the repartitioning of the data was found to be expensive and hence a
copy of both AC and U were maintained on each PE and only the operations were
distributed among the processors.
One of the primary concerns in distributing data and computation for parallel
processing is the communication overhead incurred in transferring data between pro-
cessors. The time necessary to perform a transfer consists of two parts, the time
necessary to initiate a transfer which is referred to as the latency and the time necessary
to actually transfer the data which depends on the amount of data and the
CORE
POWER
RETRAN SOLUTION
20 inner 20 outer
(psia) RETRAN SOLUTION
20 inner 20 outer
Fig. 10. Results of the RETRAN-03 with GMRES for PWR Transients (GMRES Inner): 40 sec
machine bandwidth. The following sections discuss the parallelization of the inner and
outer levels of the nested algorithm with special emphasis given to the communication
issues.
5.1. Parallelization of the Outer Iteration. One of the dominant operations
for the outer iteration is the matrix vector product. At each iteration the product is
formed of the matrix A in Eq.1 with the residual vector. This operation can be broken
into four parts. The first part involves the product AE   v E and involves no communication
since all the data resides on the same processor. The next two parts involve
U   v C and L   v E . These involve transfer of parts of the vector between processor
0, which is assigned the ex-core, and the processors 1 and 12 which are assigned the
bottommost and the topmost planes in the core, respectively. The fourth part of the
matrix vector product involves AC   v C and requires communication of processors 1
through 12 with at most two processors. In general the matrix vector product requires
communication between selected processors and has a specific pattern.
The vector inner products, on the other hand, require global communication which
means that every processor requires some information from all the other processors.
This is because it is necessary to sum the product of each element of the first vector
with the corresponding element from the second vector. The elements of the vectors
are first multiplied and each processor forms a partial sum with the elements that
reside in its domain. Because only the sum of the products is required and not the
individual products, the number of transfers required can be considerably reduced
by employing the well known method of tree-summing for which the communication
costs vary as dlog 2 (N )e (as opposed to N in the case of all to all communication).
The Gram-Schimdt orthogonalization process used in GMRES involves several inner
products during each iteration. However, the result of any of these inner products is
not dependent on any of the others. Therefore, the partial sums of all of the inner
products can be performed at one time, further reducing the number of transfers
required [2].
The least square solution in GMRES involves a reduced linear system and does
not involve sufficient operations to merit parallelization [11]. Communication related
to the least squares problem was avoided entirely by performing the least square
solution simultaneously on all the processors. Also, because the estimate of the residue
in GMRES is a consequence of the least squares solution, the termination criterion
could be evaluated without the need for additional communication.
5.2. Parallelization of the Inner Iteration. Several of the operations required
to parallelize the inner iteration were similar to the outer iteration. The
matrix AC is striped row wise and the row corresponding to each plane is stored
in the corresponding processor. Because the planes are linked only to their immediate
neighbors, the matrix vector product requires two sets of data transfer. The vector
inner products and the least square solution are treated in the same manner as in the
outer level. Since the preconditioner for the second level is a block Jacobi, the preconditioner
equation, could be implemented without any data transfers, irrespective
of whether it is solved using a direct solver or an iterative solver.
5.3. Results: PWR Model with 9 Channel/ 12 Axial Core. The parallelized
version of RETRAN-03 was executed on the Paragon with 13 processors. This
section presents the result of the rod withdrawal case for the PWR model with the
9 channel/12 axial volume core. The transient was analyzed for 20 seconds and the
nested FGMRES solution was performed using a tolerance of 10 \Gamma11 and an iteration
limit of 20 on both the inner and outer iterations.

Table
Serial Performance of the Second Level: 9 Channel/12 Axial
Time step Outer Inner Serial Execution time
Precond Other Total

Table
Parallel Performance of the Second Level: 9 Channel/12 Axial
Time step Outer Inner Parallel Execution time
Precond Comm Other Total Speedup

Tables

6 and 7 show the results of serial and parallel implementation, respectively,
of the inner or second level GMRES. The fourth column in both Tables headed "Pre-
cond" indicates the time required for the preconditioning of the second level (GMRES
III or Direct Solve in Figure 6) which was performed in this application using a direct
solver by means of LU factorization. The first timestep involves additional initialization
costs and consistent comparisons between the parallel and serial versions of
RETRAN were not possible because of differences in code structure for initialization.
The high speedup (6.22) in the first iteration of the second timestep is due to
LU factorization which is performed concurrently in the parallel implementation. It
should be noted that the bulk of the time is spent on the third level and since this part
of the code is naturally parallelizable, one would expect high efficiencies. However,
since the problem size is relatively small, the communicationoverhead in implementing
the remainder of GMRES significantly reduces the efficiency. This would not be the
case for larger problems as will be seen in the next section.

Tables

8 and 9 show the results of serial and parallel implementation of the outer
level GMRES. Again the bulk of the time is spent in the preconditioning. Similar
execution times were obtained for other timesteps. Speedup of slightly greater than
two were obtained for most of the outer iterations.
The execution time for the first four time steps is summarized in Table 10. Again,
a speedup on the order of two was achieved.
Serial Performance of the Outer Level: 9 Channel/12 Axial
Timestep Outer Serial Execution time
Core Ex-core Other Total

Table
Parallel Performance of the Outer Level: 9 Channel/12 Axial
Timestep Outer Parallel Execution time
Core Ex-core Comm Other Total Speedup
5.4. Results: PWR Model with 25 Channel/ 12 Axial Core. In order
to examine the scalability of these results, the same PWR model was used but with
instead of 9 channels in the core. The linear system is nearly three times the size
of that arising form the 9 Channel/ 12 Axial case. As anticipated, and as shown in

Table

11 the time required to solve this linear system using the direct solver was an
order of magnitude larger than that for the 9 Channel/ 12 Axial case. (Compare with
case A.1 in Table 5.)
The problem was first executed using RETRAN with nested FGMRES solver on
a single processor on the Paragon. Because the problem was too large to be executed
on an ordinary node which has 32MB of RAM, a special "fat" node of the Paragon
was used which has additional memory (128 MB) to support larger applications. As
shown in Table 11, the results on a single processor were encouraging since more than
a factor of 2 improvement was achieved with GMRES compared to the direct solver.
The problem was then executed using the parallel version of the nested GMRES
solver. The domain decomposition method was exactly the same as the 9 Channel /12
Axial case with 12 processors assigned to the core and 1 to the ex-core. As expected,
the parallel efficiency was much better and the speedups were about a factor of 4 with
respect to the serial version of GMRES. Table 11 shows the comparison of results for
the first two time steps.
As indicated in the Table 11, the parallel execution of nested GMRES provides
over an order of magnitude reduction in the execution time compared to the serial
execution of the direct solver. Furthermore, since the memory requirement per node
was less than 32MB, the parallel version could be executed using standard sized
nodes of the Paragon. Thus parallelization can be the key to not just execution time
reductions but also alleviating the memory constraints for larger problems.
6. Summary and Conclusions. A nested FGMRES method was developed to
solve the linear systems resulting from three-dimensional hydrodynamics equations.
Applications were performed on practical Pressurized Water Reactor problems with
three-dimensional core models. Serial and parallel applications were performed for
Overall Performance of the Nested FGMRES: 9 Channel/12 Axial
Timestep Number of Execution time Speedup
Outers Serial Parallel

Table
Execution time for the 25 Channel/ 12 Axial case
Case Timestep Execution time per Iteration
Outer Total per
Inner other Timestep
Core Ex-core
Precond other
Direct Solver 1 - 343.310
both a 9 channel and a 25 channel version of the reactor core.
When appropriately tight convergence was enforced at each GMRES level, the
semi-iterative solver performed satisfactorily for the duration of a typical transient
problem. The serial execution time for the 9 channel model was comparable to the
direct solver and the parallel speedup on the INTEL Paragon was a factor of 2-3
when using 13 processors. For the 25 channel model, the serial performance of nested
GMRES was about a factor of 3 better than the direct solver and the parallel speedups
were in the vicinity of 4, again using 13 processors. Thus, for the 25 channel problem
over an order of magnitude reduction in the execution time was achieved.
The results here indicate that the use of semi-iterative solvers and parallel computing
can considerably reduce the computational load for practical PWR transient
calculations. Furthermore, the results here indicate that distributed memory parallel
computing can help alleviate constraints on the size of the problem that can be exe-
cuted. Finally, the methods developed here are scalable and suggest that it is within
reach to model a PWR core where all 193 flow channels are explicitly represented.
7.

Acknowledgements

. The authors appreciate the work of Mr. Jen-Ying Wu
in generating the transient results reported in this paper.



--R

Marching Algorithms for Elliptic Boundary Value Problems.
Reducing the Effect of Global Communication in GMRES(m) and CG on Parallel Distributed Memory Computers
Direct methods for sparse matrices
The Effect of Ordering on Preconditioned GMRES Algorithm
Numerical Methods for Engineers and Scientists
Supercomputing Applied to Nuclear Reactors

An Assessment of Advanced Numerical Methods for Two-Phase Fluid Flow
A Flexible InnerOuter Preconditioned GMRES Algorithm
GMRES: A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems
A Comparison of Preconditioned Nonsymmetric Krylov Methods on a Large-Scale MIMD Machine

Performance of Conjugate Gradient-like Algorithms in Transient Two-Phase Subchannel Analysis
GMRESR: A Family of Nested GMRES Methods

Some Computational Challenges of Developing Efficient Parallel Algorithms for Data-Dependent Computation in Thermal-Hydraulics Supercomputer Applications
--TR
fluid dynamics;Preconditioned GMRES;parallel computing;nuclear reactor simulation