- Email: [email protected]
odd-parity form of the linear Boltzmann equation
MATHEMATICAL COMPUTER MODELLING Mathematical
PERGAMON
and Computer
Modelling
31 (2000) 55-71 www.elsevier.nl/locate/mcm
Parallel FE Approximation of the Even/Odd-Parity Form of the Linear Boltzmann Equation Radiation
C. R. DRUMM and Electromagnetic Analysis Department, Sandia National P.O. Box 5800, Albuquerque, NM 87185-1166, U.S.A. crdrummQsandia.gov
Laboratories
J. LORENZ Department of Mathematics and Statistics University of New Mexico, Albuquerque, NM 87131, U.S.A. lorenzQmath.uxun.edu (Received
and accepted
July 1999)
Abstract-A novel solution method has been developed to solve the linear Boltzmann equation on an unstructured triangular mesh. Instead of tackling the first-order form of the equation, this approach is based on the even/odd-parity form in conjunction with the conventional multigroup discrete-ordinates approximation. The finite element method is used to treat the spatial dependence. The solution method is unique in that the space-direction dependence is solved simultaneously, eliminating the need for the conventional inner iterations, and the method is well suited for massively parallel computers. @ 2000 Elsevier Science Ltd. All rights reserved. Keywords-Radiation transport, method, Parallel computing.
Boltzmann
equation,
Even/odd
parity,
Conjugate
gradient.
1. INTRODUCTION A novel approach is described for solving the linear Boltzmann transport equation on an unstructured mesh of triangles. Rather than following the conventional source-iteration procedure (11, in the method presented here the scattering term within an energy group is included with the streaming and removal terms on the left side of the equation. The linear system is then solved by constructing a matrix, which includes both the direction and space dependence, and solving by the conjugate gradient (CG) method. The obvious disadvantage of this approach is that the number of unknowns is equal to the number of directions multiplied by the number of spatial nodes, which may be huge. There are several appealing features of approaching the problem in this way, however. Probably the most attractive feature is that the method is easy to parallelize, by using software such as The authors wish to express appreciation to A. Prinja of the University of New Mexico and to G. Scrivner, W. Fan, D. Allen, and L. Lorence of Sandia for useful discussions and support of this work. The first author was supported by the United States Department of Energy under Contract DE-AC04-94AL85000. Sandia is a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin Company, for the United States Department of Energy, and the second author was supported by NSF Grant DMS-9404124 and DOE Grant DE-FG03-95ER25235. 0895-7177/00/s - see front matter @ 2000 Elsevier Science Ltd. All rights reserved. PII: SO895-7177(99)00223-X
Typeset
by A+@-T#
56
C. R. DRUMM AND J. LORENZ
Aztec [Z], which includes an efficient parallel CG solver. Another appealing feature involves the possibility for acceleration. The conventional source iterations are replaced by CG iterations. While the convergence of the source iterations can be accelerated by techniques such as diffusion synthetic
acceleration
good preconditioner.
(DSA)
[3], acceleration
of the CG convergence depends upon finding a
The problem of finding good parallel preconditioners
is currently an active
area of research in numerical linear algebra [4,5]. The linear system is described by a sparse matrix, which is composed of an unstructured pattern of highly structured submatrices.
The submatrices result from the scattering term, which couples
together all of the directions at a finite element (FE) node. The sparsity pattern of the matrix is determined by the FE connectivity.
The Aztec software is designed to take full advantage of
the sparsity of the matrix. A 2D prototype code has been written and tested. A 2D code is useful for applications where the solution is only weakly dependent on one of the spatial variables,
such as for modeling
the radiation-induced electrical response of cables. As an example of a typical application, a triangular mesh of a coaxial cable cross section generated by the computer aided design (CAD) package Pro/E (61 is shown in Figure 1, where the circles indicate the boundaries between distinct material 3D code. of several also been
regions. A This article benchmark written for
2D code is also a useful intermediate step toward an unstructured-mesh describes the mathematical foundation of the code and presents the results problems and a comparison with experiment. A 1D version of the code has development and testing purposes.
Figure 1. Typical Pro/E triangular mesh of a coaxial cable cross section, consisting of a center conductor surrounded by a solid dielectric and an outer grounded conductor.
Ackroyd [7] provides a thorough review of the history and state of the art of FE applications to radiation transport problems. Because the CG algorithm requires a symmetric positive definite (SPD) matrix, the even/odd-parity form of the transport problem is used in this work [8]. Josef and Morel [9] applied the simplified spherical harmonic method to model electron/photon transport on a 2D cylindrical geometry, using conventional source iteration. Manteuffel and Resse1 [lo] investigated a least-squares FE solution of the transport equation, which also leads to SPD matrices, for isotropic scattering in the diffusive limit. A comparison of the two approaches, even/odd-parity form vs. least-squares, will be the subject of future work. This article begins with a brief description of the even/odd-parity transport equation, followed by a development of the FE basis of the code. Then, the results of several benchmark problems are presented, followed by results illustrating the parallel efficiency of the solver and the effectiveness
Parallel FE Approximation
57
of preconditioners on the convergence of the CG iterations. The article finishes with a discussion of the current status of the code and plans for future work.
2. DEVELOPMENT OF THE FE DISCRETE ORDINATES EVEN/ODD-PARITY TRANSPORT EQUATION The balance relation describing the motion in space, direction, and energy of neutral particles such as neutrons or photons is the linear, steady-state
Boltzmann
transport equation [I],
U, (r, /LO,E’ --) E) @ (I-, a’, E’) dCYdE’
+ Q(r, fl, E):
(2.1)
where at(r, E) is the total cross section, and @(r> f2, E) is the particle fluence. The scattering cross section, gs(r, ,~a, E’ 4 E), is the probability of scattering between different directions and energies, and Q(r, C2,E) 1s a specified external particle source distribution. The independent variables are: the spatial position r E V, the particle direction s1 E S2, and the particle energy E, 0 < E < ca. Furthermore, ~0 = s2’ . s2 is the cosine between the pre- and post-scattered directions, W and fi. Equation (2.1) is applicable to neutral-particle transport and can also model the transport of electrons for many problems of interest if Goudsmit-Saunderson modified electron cross sections are used [ll]. The energy dependence is handled by integrating
the Boltzmann
equation
energy groups, and separating the known down-scatter source component self-scatter source, resulting in the multigroup forrn of the equation
over successive
from the unknown
(2.2) =
us,g’+g
(r, ,m)QgJ(r, a’) do’ + Qg(r, W,
where 9 is the energy group index. The unknown self-scatter portion of the source term has been moved to the left side, so that everything on the right side is known from the solutions at the previous energy groups. Equation (2.2) is solved over successive energy groups from highest to lowest, in the standard way. By including the known downscatter source into the external distributed
source term and dropping i-2 V,@(r,
f2) + gt(r)@(r,
the group indices,
equation
Ci) - /cS
Cl) @ (r, f2’) dC2’ = Q(r, Ct),
for r E V
(r, 0’
and
(2.2) can be written
(2.3)
CI E S2
Vacuum (inflow) boundary conditions specify the particle fluence for all incoming directions on the boundary for r E LW and 0. n < 0, @(r, n) = @b(rr 02)) (2.4) where n is the unit outward normal on the boundary dV, and @b(r, a) is a specified The scattering and total cross sections are assumed to satisfy gs (r, 0’
f2) dfl’ < at(r),
funct.ion.
for r E V,
/ where the left side does not depend on fi (the material is isotropic), which is a reasonable assumption for most applications. This condition is sufficient to ensure that the resulting linear system is positive definite. Within a given energy group, the transport equation depends upon space and direction, with the directions coupled together by the scattering source integral. Because of the streaming term, O.V,@, which is first order in space, the transport operator is neither symmetric
58
C. R. DRUMM AND J. LORENZ
nor positive definite. Applying FE directly to this form of the transport equation is ineffective, since the resulting linear system is also neither symmetric nor positive definite. The equation can be more conveniently solved by first extracting the symmetric and antisymmetric components of the particle fluence, resulting in two decoupled equations. The corresponding integro-differential operators are second order in space. This is known as the even/odd-parity formulation of the transport equation [12]. The advantage of the even/odd-parity formulation is that the resulting linear system is SPD, so that powerful solution methods, such as CG, can be applied to their solution. Furthermore, powerful parallel linear-equation solvers, such as Aztec [2], can be directly applied to solving the linear system efficiently on massively parallel computers. 2.1.
Even/Odd-Parity
Form
of the Transport
Equation
In order to derive the even/odd-parity form, two new functions are defined, @+ and a-, extract
the symmetric and antisymmetric
@(r, Cl) xt @(r, -fit)
@*(r, a) = Even and odd parity scattering Q* (r, n), are defined similarly c,” (r, a’.
(2.5)
2
cross sections,
52) =
which
components of @(r, a),
g$(r, 52’ . f2), and distributed
cs (r, a’ .fl)
2~cs (r, --a’
source terms,
f2)
2
and
Q*(r, fi) =
Q(r, fl) f Q(r, -n2) 2
The even/odd-parity equations are derived by combinin, m equation (2.3) with the equation obtained by replacing fi by .-a in equation (2.3), and extracting equations for @+(r, a) and W(r, 0). Boundary conditions are derived from equations (2.4) and (2.5). Defining two scattering operators, &f(r, and a transport
a) = at(r)f(r,
a) - J’g$
(r, a’ . a) f (r, W
da’,
operator, 7f(r,
n) = a . V,f(r,
01,
the equation for the even-parity component of the fluence is -7GIr’7@+
+ 6+@+(r,
a) = -7611Q-
forrEV
and
(r, 02) + Q’(r,
a),
0tS2,
(2.6)
with boundary conditions +@+(r,
02) + GI17@+
= *@b (r, ?a)
for r E dV
and
S2.n
+ !X’Q-(r,
fl),
(2.7)
2 0.
Whereas the boundary conditions for the first-order form of the transport equation are specified only for incoming directions, the boundary conditions for the even-parity form of the equation are specified both for incoming and outgoing directions. This results in the proper number of boundary conditions for the even-parity transport equation, which is second order in space. The equation for the odd-parity component can be derived analogously with the derivation of the even-parity equation, or may be obtained simply by interchanging +‘s and -‘s in the even-parity equation -7G,17@-
+ G-a’-(r, forrEV
a) = -IGTlQ+(r, and
0gS2,
52) + Q-(r,
a),
(2.8)
59
Parallel FE Approximation
with boundary conditions i@-(r,
Sl) +
E;“T@- = -@b (r, $02) f GT’Q+(r, fl2), for r E 6%’
Evaluation
of the inverse scattering
and
operators
fi.n
(2.9)
2 6.
G+l and 61’
is described in the next section.
Note that the boundary conditions for the odd-parity equation cannot quite be obtained from those for the even-parity equation by interchanging fluence is symmetric
f’s and --)s. This is because the even-parity
in direction, while the odd-parity fluence is antisymmetric
the boundary term is slightly different for the odd-parity equation. and a demonstration
in direction,
so
A more detailed derivation
that the operators are symmetric and positive definite are in [8,13].
2.2. Weak Form and Discretization
of the Even/Odd-Parity
Equations
The weak forms of equations (2.6) and (2.8) are derived by multiplying the equations by a test function u(r, a), integrating over the spatial and directional domains, applying the divergence theorem, and then applying the boundary conditions, equations (2.7) and (2.9). This procedure results in the even-parity weak form
JJ J! JJ
lu(r, Cl)GIlI@+(r,
+
Cl) dr dJZ +
u(r, O)@+(r, 52) [fl. nl dsdn
+
u(r, s2)Qf dr dfl
=
+J f
JJ JJ
u(r, O)G+@+(r, Cl) dr dC2
lu(r, s2)G11Q-(r, Cl) drdi2
u(r, a)$,
n.n$l
(r, rfl)
1fl. nl dsdC2
and similarly for the odd-parity equation
JJ Sf JJ
Iu(r, fi)G517+-(r,
+ +
JJ JJ J f'
Cl) drdC2+
‘u(r, C2)@-(r,Q) IC!. nl dsdi2
u(r, O)Q-(r,
f2) dr dC2 T-
=
Cl+0
u(r, Cl&_@-(r,
Cl) drdSl
lu(r, fi)GT’Q’(r,
u(r, a)@.6 (r, #I)
C2)dr dC2
Ifi. nl dsdfi.
In a discrete-ordinates formulation, the directional variable, 0 E S2, varies only in a finite set of discrete directions s2,, and the directional integrals are replaced by quadrature sums. This is equivalent to collocation in direction, with basis functions defined by w,6(1 - 0 Cl,), where 20, is the mth quadrature weight. The spatial domain is approximated by an unstructured mesh of triangles, and a set of basis functions Gn(r) is defined within each element. For a Galerkin FE treatment, the test functions are the same as the basis functions. To proceed, we write the discrete-ordinates form of the two scattering operators as
ij*jm(r) = cr,(r)jm(r) -
5
w,~o$ (r, %2
. %I
.fd(rIr
m’=l
where M is the number of discrete directions. operator is
The discrete-ordinates
form of the transport
The discrete function fm(r) = f(r, Cl,) is the function f(r, 0) evaluated at the mth quadrature direction.
C. R. DRUMM AND J. LORENZ
60
Within a spatial element e the approximate even- and odd-parity solutions in a discrete direction C&, are assumed to be combinations of the basis functions for r E V,,
e = 1: NE,
m = 1 : M,
where $&,,,, are the values of the solutions at the nodes of the element, N is the number of basis functions per element, NE is the number of spatial elements in the mesh, and V, is the element volume. The distributed
source is defined similarly by for r E V,,
e = 1: NE,
m = 1 : Ad,
and the discrete boundary conditions are for r E AX,
f12, . n 2 0,
where the boundary conditions are defined on the external boundary tW, of the elements e which have an external boundary. With these discrete functions and operators
defined, the FE discrete-ordinates
form of the
even-parity weak form is
forn=l:N,
e=l:NE,
CZ,.n20,
m=l:M.
The volume integrals are carried out over the element volumes, and the surface integrals are carried out over the external boundary of the domain. The odd-parity weak form can be discretized similarly
for n = 1: N,
e = 1: NE,
ST&,, ’ n p 0,
m = 1: M.
Equations (2.10) and (2.11) define linear systems of dimension N * M, where N is the number of FE nodes in the mesh, and the unknowns c#&,,, are the values of the even/odd-parity fluence at the nodes and discrete directions.
Parallel FE Approximation
In order to setup equations tions of the scattering represented constant
(2.10)
operators
and (2.11),
of the scattering scattering
expansions.
homogeneous
operators
operators
expressions
for the discrete
approxima-
c+ and g__, ancl their inverses are needed, and they are commonly
by finite Legendre-polynomial
within a spatially
explicit
61
is dropped.
Since
material
region of the problem
data
geometry,
The Legendre-polynomial
are assumed
the spatial
approximations
to be
dependence
of the discrete
are
and
where
L is the order of the Legendre
Legendre
polynomial
expansion
of the scattering
of degree 1. The inverse operators
cross section.
and Pl is the
are [12]
and
Based
on these approximations,
to solve it with a parallel article.
the computational
problems,
a code has been developed
CG algorithm results
are compared
and then the computational
more realistic
in the Aztec results
to set up the linear system,
code package.
with analytic are compared
and then
In the next section
solutions,
for several
of this
simple
test
data
for a
with some experimental
problem.
3. RESULTS In this section
the convergence
rate of the error
imations
of the code is investigated
solutions
are available,
aluminum.
If an exact
followed solution
due to the spatial
by analyzin g several
by a comparison is known,
a discrete
sirnple test
with measured ordinates
and directional problems,
data
where
for 1-MeV
approximation
approxanalytic
electrons
on
of the L” norm of
the error is defined by
E=
where f(r,
!&)
of these
results
nontrivial, problems
results
rate,
h is the maximum
are designed
element
for elliptic
point,
and the exact
is 2, and that
size, and c is a positive
boundary
of the even/odd-parity
and will be presented
basis functions
model is used.
i(r,%)/2 dr.
at the mu’ quadrature
are well established
to the discretization
however,
with linear
solution
-
solution
to obey the relation
r is the convergence
Such convergence
f(r,G)
w,
is the approximate
is .f(r, %J. The error is expected
where
c 1’ J..
value problems.
form of the Boltzmann
in future work. The expected
convergence
for quadratic
is 3.
to test the convergence
basis functions
of the spatial
approximation,
The
const,ant.
An extension equation
is
rate for FE first two test
so a simple physics
C. R. DRUMM AND J. LORENZ
62
3.1. 1D Test Problem Consider a unit interval in space, 0 5 z < 1, unit total cross section, ut = 1, scattering and absorption equally probable, and linearly-anisotropic
scattering, a,@~)
= (1/2n)(1/4+pc/4),
corresponds to Legendre moments, 0~0 = l/2, LT,~= l/6, where gsl = 2~s:~ With these assumptions,
where the coordinate
the first-order form of the Boltzmann
system is illustrated
which
as(po)Pl(~~O) dpo.
transport equation is
in Figure 2 and ~0 is given below. The 1D approxi-
mation assumes that the fluence, source, scattering cross section, and boundary conditions are independent of two of the spatial variables x and y ($$’ = $$ = 0), and also independent of the azimuthal direction (p, which is the angle between the positive x-axis and the projection of the particle direction on the zy-plane. These approximations result in an equation in two independent variables z and p, where p = cos(B), and 13is the angle between the particle direction and the positive z-axis.
Figure 2. Coordinates equation.
for 1D and 2D approximations
to the Boltzmann
The cosine of the angle between the pre- and post-scattered written in terms of the pre- and post-scattered directions as
directions
transport
is ~0, which can be
PO = pL/l’+ sin(e) sin(8’) cos (‘p - cp’) . Since the integral over the azimuthal angle cp’ of cos(cp - cp’) is zero, the transport simplifies to
equation
(3.2) In order to test the error in the numerical results, we specify a simple solution
@(z,p) = (1 + p)e-“. The source term and boundary conditions that result in the desired solution are
Q(GIL) =
(-
63
ParallelFE Approximation and
a@,P) =
1+
forO
Pu,
and
+b(l,p) Computational
= (1 +
for - 1 5 /I < 0.
/.+-‘I
results are determined by specifyin, 0 a uniform spatial mesh with linear or
quadratic basis functions, and a low-order Gauss quadrature for approximating the p-integrals. Because of the simple angular dependence of this example problem, a 2-point Gauss quadrature is sufficient to exactly model the angular dependence, so any error in the computational
results
arise from the spatial approximation. At this point, it is interesting to write out equations for the even/odd scat,tering operators and integro-differential equations for the even- and odd-parity solutions. The even-parity fluence is
@+(z,p) =
e-‘,
and the odd-parity fluence is
a-(z, p) = pee”. The even-parity source is
Q+(z, P) =
(-,u2+ ;) e-z,
and the odd-parity source is Q-(z,p)
= -Eeaz 6
’
The scattering operators are
G+f(z, p) = f(z,
CL) - ;
J:f(z1P’)dP’
and
G_f(z, p) = f(z, p) - f /_,
PP’f(Z,P’)dl.L’.
The inverse operators are 1
o,1.f(w4 = f(z,CL)+ Fj
Jf l
_l
(z, CL’) dP’
an d
The even-parity transport equation is _
p2
a2@+ (z,PL)-- l”. /; dZ2
p2pnazm;~;‘p’) dp’ + O+(z,p) =
-p&
Q-(z,
- ;
s_: @+ dP’
3 J _’p/IQ-
p) + 1o
1
(Z>P’)
(z, P’) dd
I
+ Q+(G ~1,
C. R. DRUMM AND J. LORENZ
64
with boundary
conditions
for p >< 0, z = 0, The odd-parity
transport
equation
_/~2a2@-(zd4 -- 1 l a9
2 J _r PP
~20,
z=l.
is
,a2a- (6 i-4 dp’ d.9
+ a-- (z, p) - ;
Q+k
= -P$
with boundary
and
CL)+ ;
J w'@-(~,l~') J1 1 _’
1
_’ Q+ k P’) @’
+ Q- (2, P)>
conditions
=
-@b
(2, F/J) + Q+(z, P) + ; for p 5 0, z = 0,
and
/’ Q+ (6 P’)dd, 1
~2O,z=l.
The results are given in Tables 1 and 2, which list the error norms vs. number of uniform elements. Also given in the table are the convergence rates, which are determined from a linear least-squares fit of the log of the error norm vs. the log of the mesh width. The expected convergence rates of 2 for linear basis functions and 3 for quadratic basis functions are obtained. The error-norm 3.2.
results
are plotted
in Figure
3.
2D Test Problem
Consider a unit square, 0 5 x < 1 and 0 < y 5 1, unit total cross section, crt = 1, scattering and absorption equally probable, linearly-anisotropic scattering, CT,(~0) = (1/27r) (l/4 + po/4), which corresponds to Legendre moments of us0 = l/2 and us1 = l/6. With the assumption transport that the solution is independent of Z, $$ = 0, the first-order form of the Boltzmann equation is
where the coordinate system is illustrated in Figure 2. The cosine of the angle between the pre- and post-scattered
directions
is ~0,
~0 = pL,fii + sin(d) sin (0’) cos (‘p - cp’), where CL2= cos(@, where 8 is the angle between the particle direction and the positive z-axis, and cp is the azimuthal angle between the positive x-axis and the projection of the particle direction on the xy plane.
65
Parallel FE Approximation Table 1. Error in the numerical functions.
results for the 1D test problem
for linear basis
L2 Norm of Error: Linear Basis Functions Number of Elements Even-Parity
I (
Fluence
Odd-Parity
Fluence
5
1.98 x 10-3
8.39 x 10-d
10
4.95 x 10-d
2.10 x 10-a
20
1.24 x 10-d
5.25 x IO-”
40
3.10 x 10-S
1.31 x 10-S
80
7.74 x 10-s
3.28 x 10-s
160
I
Convergence
Rate
-r~--
1.93 x 10-s
8.21 x 10-7
2.00
[
2.00
I
I
Table 2. Error in the numerical results for the 1D test problem for quadratic functions. L2 Norm of Error: Quadratic
Number of Elements
Even-Parity 5 10
I
3.78 x 10-G
I
20
I
4.73 x 10-7
-
1o-g
Fluence
1.74 x 10-S ~ 2.18 ~~ x 1O-6
_(
2.73 x 1O-7
I
40
5.91 x 10-s
3.41 x 10-s
80
7.38 x 10-s
4.27 x lo-”
Rate
3.00
3.00
I
I
I
I
i
I
I
I,
40 l/(maximum
-1 I
I
’ \w
quadratic interpolation functions convergence rate 3.00 I
Basis Functions
Odd-Parity
3.01 x 10-s
I
Convergence
E
Fluence
basis
element size)
Figure 3. Error analysis for 1D test problem
III,
80
/II,-
120 160
66
C. R. DRUhml
AND
J. LORENZ
If the solution is independent of azimuthal angle cp, the integral over the azimuthal angle cp of cos(cp - cp’) is zero, and the transport equation is
Assume a solution Q,(2, Y, k)
= cL%cmZ.
This corresponds to a distributed source term of
Q(x,Y,P~,P~)
=
P: (1 -PX) -
I
$ e-“,
and boundary condition %
(2,Y, CL,)= I-L%-x,
for (z, y) E IX
The even-parity fluence is @+ (.‘c,Y, L)
= &ePZ,
and the odd-parity fluence is a- (G y, PZ) = 0. The even-parity source is
Q+ (x,Y,P,,P~)
=
(112- i >emz,
and the odd-parity source is
With linear basis functions, Sis quadrature (72 directions [14]), the L2 error norms are listed in Table 3 and plotted in Figure 4, showing nearly the expected quadratic convergence rate. Table 3. functions.
Error
in the numerical
results
for the 2D test
problem
for linear
basis
Parallel FE Approximation
/
/
I
I
I
I I /
25 l/
/,I, 50
,,I, 75 100
(maximum element size)
Figure 4. Error analysis for 2D test problem
3.3 Convergence of Directional Approximation The third test problem is designed to test the convergence of the directional approximation. Consider a unit interval in space, 0 5 z 5 1, unit total cross section, CT t = 1, isotropic scattering, with scattering ratio c, which is defined as the ratio of the scattering cross section to t,he total cross section. The scattering transport equation is
ratio lies in the range, 0 5 c < 1. For this problem,
_P2 a2@+ (6 P) + @+(z,p) a9 with boundary
1
= ;
CD+(2,~‘) dp’ + Q+(z,p). s -1
conditions
==t@, (z, &cl) + Q-(2,/L). A test problem
the even-parity
is defined
as follows:
the even-parity
for St. n 2 0,
z = 0.1.
(3.3)
fluence is
@+(Z, /L) = COS(l_L), and t,he odd-parity
fluence is zero, a- (2. /L) = 0,
which corresponds
to source terms of &+(a, p) = “OS(~) - csin(1)
an d Q-(z,/L)
= 0.
The convergence of the directional approximation is examined by dividing the directional domain into a number of uniform intervals, and approximating the integral over each subregion by either a l-point or 2-point Gauss quadrature. Using one Gauss point per angle mesh should result in quadratic convergence, and using two Gauss points per angle mesh should result in fourthorder convergence. This is because one-point Gauss quadrature integrates up to linear terms exactly, and two-point Gauss quadrature integrates up to cubic terms exactly. Figure 5 shows the convergence of the L2 norm of the directional approximation. showing nearly the theoretical convergence rates. The results are tabulated in Table 4.
C.
68
R. DRUMM AND J. LORENZ
1C3
lo9
lOA
104
10"
10"
10"
IO8
lo-'
10.' I -
10” :
2 Gauss convergence
points per angle mesh rate
‘L
_
=
4.00 I
. I
5
1oe9
Number
of
Angle
Number of Angle Meshes
3.4. Comparison
2
3.98 x 1O-3
5.71 x 10-s
4
1.01 x 10-S
3.54 x 10-T
8
2.54 x IO-*
2.22 x 10-8
16
6.36 x 1O-5
1.38 x 1O-g
1.99
4.00
with Experimental
10”
rates for test problem 3. 2 Gauss Points
Rate
15
approximation.
1 Gauss Point
Convergence
I
Meshes
Figure 5. Error analysis for convergence of the directional Table 4. Error norms and convergence
10
10"
.
Data
Lockwood et al. [15] used a calorimetric technique to measure energy-deposition profiles in 1D for electron sources. In this section, we compare the results of one of their sets of measurements with FE modeling. The reported experimental uncertainties are less than a few percent. The problem considered here is that of l-MeV electrons normally incident on aluminum. The parameters for the FE modeling are: 40 uniform electron energy groups, no coupling with photons, Pi5 (15th-order) Legendre expansion of the scattering cross section directional dependence, Sic (16-point) Gauss quadrature, 50 uniform spatial elements, with quadratic basis functions. The comparison between calculation and experiment is shown in Figure 6, with excellent agreement . 3.5.
Parallel
Efficiency
In order to test the parallel efficiency of the code, a 2D test problem was run on Sandia’s ASCI red computer, utilizing from one to 512 CPU’s. The problem considered is the 2D benchmark problem described in Section 3.2. Two different meshes were considered, a mesh with a maximum element size of 0.01 with 10,201 FE nodes, and a coarser mesh with a maximum element size of 0.02 with 2,601 FE nodes. For the fine mesh, three different quadrature orders were considered, Ss, Sis, and Sic quadrature, with 20, 42, and 72 discrete directions, respectively [14]. For the coarser mesh, only Sls quadrature was considered. Figure 7 is a plot of the speedup versus the number of processors for the four problem, where the speedup S is defined by
Parallel 3.5
I’
4
i
‘I
Lockwood
69
FE Approximation
I
/
’
8
I”’
I
I
data
Fraction Pigllre 6. Energy deposition
of
mean
range
profiles for l-Me11 electrons
on aluminum.
500
450
+ __C__
400
h=O.Ol, h=O.Ol, h=O.Ol, h=0.02,
S,,, S,*, S,, !Z&
734,472 428,442 204,020 187,272
unknowns unknowns unknowns unknowns
n 350 z g 300 n cn 250 z? z2
200 150 100
200
300
Number of Processors Figure
where
7. Parallel
Ir, and T, are the runtimes
required
for the matrix
quite good, improving the larger
problems
speedup on teraflops
for 1 and n processors,
and right-hand as the number
have better
is less for the larger problems.
computer
for four matrix
respectively,
sizes.
which include
side setup and for the CG solve. The parallel of unknowns
parallel
efficiencies,
The convergence
increases.
For a given number
since the relative
criterion
fraction
the time
efficiency
is
of processors,
of communication
is that the relative error norm is reduced
C. R. DRUMM
70 to
10m8 of its original value, /jrll/llr~ll < lo-‘.
AND J. LORENZ
For the largest problem considered, speedup of 466
was attained for 512 processors, corresponding to better than 90% parallel efficiency. 3.6.
Preconditioning
The CG iterations converge fairly slowly when applied to the unpreconditioned but by preconditioning
the matrix, the convergence can be accelerated.
linear system,
Table 5 compares run
times and iteration counts for the unpreconditioned system compared to the results with several preconditioners that are included in the Aztec software. The problem considered is again the 2D benchmark problem described in Section 3.2. Only the largest problem was considered, with a maximum element size of 0.01, Sls quadrature 734,472
unknowns.
The calculations
with 72 discrete directions,
corresponding
were performed on Sandia’s ASCI red computer,
to
utilizing
100 CPU’S. Table 5. Comparison
of Aztec preconditioners. Aztec Solver
Preconditioner
Time (s) 468
2,879
point Jacobi
226
1,384
l-step block Jacobi
271
1,420
3-step block Jacobi
678
1,207
none
I
( l-order Neumann polynomial
)
626
1
1,933
3-order Neumann polynomial
a73
1,362
l-order least squares polynomial
667
2,062
I-3-order
least squares polynomial
710
I
(
1,107
ilnt(O.0)
439
754
ilut(0.02)
446
784
ilut(0.04)
910
1,721
Only two of the built-in preconditioners significantly reduced runtime, point Jacobi and l-step block Jacobi. Incomplete LU with threshold cutoff (ilut) significantly reduced the iteration count, but without much change in the runtime. The various preconditioners are described in the Aztec manual [2] and are briefly described here. Point Jacobi is equivalent to a symmetric diagonal scaling, resulting in all 1s on the diagonal. Block Jacobi uses a small number of block Jacobi iteration steps as a preconditioner, with the block size equal to the number of discrete directions in the quadrature integration. Neumann and least squares use a low-order polynomial approximation as a preconditioner, and ilut uses an incomplete LU factorization with a threshold cutoff as a preconditioner. Much work remains to be done to find a good preconditioner for this linear system. The convergence criterion is that the relative error norm is reduced to 10m8 of its original value Ilrjl/llroll< lops. 4.
CONCLUSIONS
A computer code has been written to solve the linear Boltzmann transport equation on an unstructured mesh of triangles, and the code has been applied to a number of benchmark problems to test the spatial and directional discretizations. The solver is based on a triangular mesh of a 2D region, but extension to 2D quadrangular or mixed triangular-quadrangular mesh or 3D tetrahedral or hexahedral mesh would be straightforward. The 2D FE approximation is currently based on linear basis functions, which are fairly inaccurate in regions with steep gradients
Parallel
in t,he solution, included
common
in radiation
in a 1D version
Because
the method advantage
from the scattering elements
source-iteration designed than
problem,
experience
for photon
by takin, c advantage
Work
is in progress
CG converges
algorithm.
for this
Interestingly,
term.
For this possibly
indicates
groups,
unlike
basis functions
simultaneously, of the structure
to exploit
submat,rices,
based
an effective
this
t,he standard
source
system
of t,he matrix,
of the submat~riccs
structure
substantially
by c*omputing
reducing
storage.
as does the 11naccPleratc:cl
preconditioner
on a coarse-mesh
that the CG iterations
the linear
of the spsrsity
slowly for large problems,
reason,
have been
in the 2D code would be straightforward. dependence
on the fly from a few precomputed
unaccelerated
Quadratic
is designed to take full advantage
can be gained
arising
71
problems.
of the code, and inclusion
matrix
Furthermore,
transport,
solves the space-direction
may be huge. The Aztec software but flIrther
FE Approximation
is neetletl
space or direction
specifically
ilpproximnt,ioll.
converge much faster for electron iterations.
This
behavior
groups
mc,rit,s further
investigation. ‘l%
co&+ is designed
t,ransport,
of charged
tions [ll]. sections
The reference [IO], resulting
port codes. the electron a continuous
to model
particles,
describes in electron
For some applications, cross sections,
the transport
such as electrons,
of neutral
a Goudsmit-Saunderson cross sections this method
which is inefficient.
slowing down term directly
particles
by usin g specially
modification
that, are compatible requires Possibly,
in the transport
but
modified
can
also model
electron
to CEPXS
a better equat,ion.
electron
with nerltral-pirrt,ic:le
a very high order Legendre approach
t,he
cross setcross traus-
expansion
of
would be to inclllde
and this work is in progrclss.
REFERENCES G. Bell and S. Glasstone, N?lclear Reactor Theory, Va11 Nostrand Reinhold. New York, (1970). S. Hutchinson, J. Shadid and R. Tuminaro, Az~.ec user’s guide version 1.1. Report SAND!%15.59. Sarldia National Laboratories, (1995). 3. E.E. Lewis and W.F. Miller, Jr., Computational Methods of Neutron Transport, Wiley, New York. (15~4). PA. (1997). 4. .J.W. Demmel, Applied Numerical Linear Algebra, SIAM, Philadelphia. Y. Saad, Iterative Methods for Sparse Linear Systems. ITP. Boston, MA. (1996). Parametric Technology Corporation, Pro/MESH User’s Cuzde, Release 17.0, (19%) R.T. Ackroyd, Finite Element Methods for Particle Transport, Research Studies Press, Somerset. (1997). 8. V. Agoshkov, Boundary Value Problems for Transport Equations, Birkhauser, Boston, hIA, (1998). trnnsljort 9. J.A. Josef and J.E. Morel. Simplified spherical harmonic method for coupled electron-photon calculations, Phy. Rev. E 57, 6161 (1998). finite-element, solut,ion of the neutron transport, quation in 10. T.A. Manteuffel and K.J. Ressel, Least-squares diffusive regimes, SIAAl J. Numer. Anal. 35, 806 (1998). C.R. Drumm. hlultidimensional electron-photon transport with standard discrete ordinates ctrtles, Nuc~. Il. Sc1. nnd Eng. 127 (1) (1997). and W. hlartin, Transport Theory, Wiley. New York. (1979). 12. J. Duderstadt transport analysis on 2-D unstrucl.llred mesh. llcport 13. C’ R. Drumm, Parallel finite element electron-photon SAND99-0098, Sandia National Laboratories, (1999). for nuclear analyses: Theory and guidelines for 14. R.D. O’Dell and R.E. Alcouffe, Transport calculations effective use of transport codes, Report LA-10983-MS, Los Alamos National Laboratory, (1987). measurement of electroll energy 15. G.J. Lockwood, L.E. Ruggles, G.H. Miller and .J.A. Halbleib. Calorimetric deposition in extended media-theory vs. experiment, Report SAND79-0414, Sandia Nat,ional Laborat,ories. (1979). A multigroup c~ouplrtl elect,ron 16. L.J. Lorence, Jr., .J.E. Morel and G.D. Valdez, Physics guide to CEPXS: photon cross-section generating code, Version 1.0; Report S.4ND89-1685, Sandia National Lal,oratories, (1989). 2.
PERGAMON
and Computer
Modelling
31 (2000) 55-71 www.elsevier.nl/locate/mcm
Parallel FE Approximation of the Even/Odd-Parity Form of the Linear Boltzmann Equation Radiation
C. R. DRUMM and Electromagnetic Analysis Department, Sandia National P.O. Box 5800, Albuquerque, NM 87185-1166, U.S.A. crdrummQsandia.gov
Laboratories
J. LORENZ Department of Mathematics and Statistics University of New Mexico, Albuquerque, NM 87131, U.S.A. lorenzQmath.uxun.edu (Received
and accepted
July 1999)
Abstract-A novel solution method has been developed to solve the linear Boltzmann equation on an unstructured triangular mesh. Instead of tackling the first-order form of the equation, this approach is based on the even/odd-parity form in conjunction with the conventional multigroup discrete-ordinates approximation. The finite element method is used to treat the spatial dependence. The solution method is unique in that the space-direction dependence is solved simultaneously, eliminating the need for the conventional inner iterations, and the method is well suited for massively parallel computers. @ 2000 Elsevier Science Ltd. All rights reserved. Keywords-Radiation transport, method, Parallel computing.
Boltzmann
equation,
Even/odd
parity,
Conjugate
gradient.
1. INTRODUCTION A novel approach is described for solving the linear Boltzmann transport equation on an unstructured mesh of triangles. Rather than following the conventional source-iteration procedure (11, in the method presented here the scattering term within an energy group is included with the streaming and removal terms on the left side of the equation. The linear system is then solved by constructing a matrix, which includes both the direction and space dependence, and solving by the conjugate gradient (CG) method. The obvious disadvantage of this approach is that the number of unknowns is equal to the number of directions multiplied by the number of spatial nodes, which may be huge. There are several appealing features of approaching the problem in this way, however. Probably the most attractive feature is that the method is easy to parallelize, by using software such as The authors wish to express appreciation to A. Prinja of the University of New Mexico and to G. Scrivner, W. Fan, D. Allen, and L. Lorence of Sandia for useful discussions and support of this work. The first author was supported by the United States Department of Energy under Contract DE-AC04-94AL85000. Sandia is a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin Company, for the United States Department of Energy, and the second author was supported by NSF Grant DMS-9404124 and DOE Grant DE-FG03-95ER25235. 0895-7177/00/s - see front matter @ 2000 Elsevier Science Ltd. All rights reserved. PII: SO895-7177(99)00223-X
Typeset
by A+@-T#
56
C. R. DRUMM AND J. LORENZ
Aztec [Z], which includes an efficient parallel CG solver. Another appealing feature involves the possibility for acceleration. The conventional source iterations are replaced by CG iterations. While the convergence of the source iterations can be accelerated by techniques such as diffusion synthetic
acceleration
good preconditioner.
(DSA)
[3], acceleration
of the CG convergence depends upon finding a
The problem of finding good parallel preconditioners
is currently an active
area of research in numerical linear algebra [4,5]. The linear system is described by a sparse matrix, which is composed of an unstructured pattern of highly structured submatrices.
The submatrices result from the scattering term, which couples
together all of the directions at a finite element (FE) node. The sparsity pattern of the matrix is determined by the FE connectivity.
The Aztec software is designed to take full advantage of
the sparsity of the matrix. A 2D prototype code has been written and tested. A 2D code is useful for applications where the solution is only weakly dependent on one of the spatial variables,
such as for modeling
the radiation-induced electrical response of cables. As an example of a typical application, a triangular mesh of a coaxial cable cross section generated by the computer aided design (CAD) package Pro/E (61 is shown in Figure 1, where the circles indicate the boundaries between distinct material 3D code. of several also been
regions. A This article benchmark written for
2D code is also a useful intermediate step toward an unstructured-mesh describes the mathematical foundation of the code and presents the results problems and a comparison with experiment. A 1D version of the code has development and testing purposes.
Figure 1. Typical Pro/E triangular mesh of a coaxial cable cross section, consisting of a center conductor surrounded by a solid dielectric and an outer grounded conductor.
Ackroyd [7] provides a thorough review of the history and state of the art of FE applications to radiation transport problems. Because the CG algorithm requires a symmetric positive definite (SPD) matrix, the even/odd-parity form of the transport problem is used in this work [8]. Josef and Morel [9] applied the simplified spherical harmonic method to model electron/photon transport on a 2D cylindrical geometry, using conventional source iteration. Manteuffel and Resse1 [lo] investigated a least-squares FE solution of the transport equation, which also leads to SPD matrices, for isotropic scattering in the diffusive limit. A comparison of the two approaches, even/odd-parity form vs. least-squares, will be the subject of future work. This article begins with a brief description of the even/odd-parity transport equation, followed by a development of the FE basis of the code. Then, the results of several benchmark problems are presented, followed by results illustrating the parallel efficiency of the solver and the effectiveness
Parallel FE Approximation
57
of preconditioners on the convergence of the CG iterations. The article finishes with a discussion of the current status of the code and plans for future work.
2. DEVELOPMENT OF THE FE DISCRETE ORDINATES EVEN/ODD-PARITY TRANSPORT EQUATION The balance relation describing the motion in space, direction, and energy of neutral particles such as neutrons or photons is the linear, steady-state
Boltzmann
transport equation [I],
U, (r, /LO,E’ --) E) @ (I-, a’, E’) dCYdE’
+ Q(r, fl, E):
(2.1)
where at(r, E) is the total cross section, and @(r> f2, E) is the particle fluence. The scattering cross section, gs(r, ,~a, E’ 4 E), is the probability of scattering between different directions and energies, and Q(r, C2,E) 1s a specified external particle source distribution. The independent variables are: the spatial position r E V, the particle direction s1 E S2, and the particle energy E, 0 < E < ca. Furthermore, ~0 = s2’ . s2 is the cosine between the pre- and post-scattered directions, W and fi. Equation (2.1) is applicable to neutral-particle transport and can also model the transport of electrons for many problems of interest if Goudsmit-Saunderson modified electron cross sections are used [ll]. The energy dependence is handled by integrating
the Boltzmann
equation
energy groups, and separating the known down-scatter source component self-scatter source, resulting in the multigroup forrn of the equation
over successive
from the unknown
(2.2) =
us,g’+g
(r, ,m)QgJ(r, a’) do’ + Qg(r, W,
where 9 is the energy group index. The unknown self-scatter portion of the source term has been moved to the left side, so that everything on the right side is known from the solutions at the previous energy groups. Equation (2.2) is solved over successive energy groups from highest to lowest, in the standard way. By including the known downscatter source into the external distributed
source term and dropping i-2 V,@(r,
f2) + gt(r)@(r,
the group indices,
equation
Ci) - /cS
Cl) @ (r, f2’) dC2’ = Q(r, Ct),
for r E V
(r, 0’
and
(2.2) can be written
(2.3)
CI E S2
Vacuum (inflow) boundary conditions specify the particle fluence for all incoming directions on the boundary for r E LW and 0. n < 0, @(r, n) = @b(rr 02)) (2.4) where n is the unit outward normal on the boundary dV, and @b(r, a) is a specified The scattering and total cross sections are assumed to satisfy gs (r, 0’
f2) dfl’ < at(r),
funct.ion.
for r E V,
/ where the left side does not depend on fi (the material is isotropic), which is a reasonable assumption for most applications. This condition is sufficient to ensure that the resulting linear system is positive definite. Within a given energy group, the transport equation depends upon space and direction, with the directions coupled together by the scattering source integral. Because of the streaming term, O.V,@, which is first order in space, the transport operator is neither symmetric
58
C. R. DRUMM AND J. LORENZ
nor positive definite. Applying FE directly to this form of the transport equation is ineffective, since the resulting linear system is also neither symmetric nor positive definite. The equation can be more conveniently solved by first extracting the symmetric and antisymmetric components of the particle fluence, resulting in two decoupled equations. The corresponding integro-differential operators are second order in space. This is known as the even/odd-parity formulation of the transport equation [12]. The advantage of the even/odd-parity formulation is that the resulting linear system is SPD, so that powerful solution methods, such as CG, can be applied to their solution. Furthermore, powerful parallel linear-equation solvers, such as Aztec [2], can be directly applied to solving the linear system efficiently on massively parallel computers. 2.1.
Even/Odd-Parity
Form
of the Transport
Equation
In order to derive the even/odd-parity form, two new functions are defined, @+ and a-, extract
the symmetric and antisymmetric
@(r, Cl) xt @(r, -fit)
@*(r, a) = Even and odd parity scattering Q* (r, n), are defined similarly c,” (r, a’.
(2.5)
2
cross sections,
52) =
which
components of @(r, a),
g$(r, 52’ . f2), and distributed
cs (r, a’ .fl)
2~cs (r, --a’
source terms,
f2)
2
and
Q*(r, fi) =
Q(r, fl) f Q(r, -n2) 2
The even/odd-parity equations are derived by combinin, m equation (2.3) with the equation obtained by replacing fi by .-a in equation (2.3), and extracting equations for @+(r, a) and W(r, 0). Boundary conditions are derived from equations (2.4) and (2.5). Defining two scattering operators, &f(r, and a transport
a) = at(r)f(r,
a) - J’g$
(r, a’ . a) f (r, W
da’,
operator, 7f(r,
n) = a . V,f(r,
01,
the equation for the even-parity component of the fluence is -7GIr’7@+
+ 6+@+(r,
a) = -7611Q-
forrEV
and
(r, 02) + Q’(r,
a),
0tS2,
(2.6)
with boundary conditions +@+(r,
02) + GI17@+
= *@b (r, ?a)
for r E dV
and
S2.n
+ !X’Q-(r,
fl),
(2.7)
2 0.
Whereas the boundary conditions for the first-order form of the transport equation are specified only for incoming directions, the boundary conditions for the even-parity form of the equation are specified both for incoming and outgoing directions. This results in the proper number of boundary conditions for the even-parity transport equation, which is second order in space. The equation for the odd-parity component can be derived analogously with the derivation of the even-parity equation, or may be obtained simply by interchanging +‘s and -‘s in the even-parity equation -7G,17@-
+ G-a’-(r, forrEV
a) = -IGTlQ+(r, and
0gS2,
52) + Q-(r,
a),
(2.8)
59
Parallel FE Approximation
with boundary conditions i@-(r,
Sl) +
E;“T@- = -@b (r, $02) f GT’Q+(r, fl2), for r E 6%’
Evaluation
of the inverse scattering
and
operators
fi.n
(2.9)
2 6.
G+l and 61’
is described in the next section.
Note that the boundary conditions for the odd-parity equation cannot quite be obtained from those for the even-parity equation by interchanging fluence is symmetric
f’s and --)s. This is because the even-parity
in direction, while the odd-parity fluence is antisymmetric
the boundary term is slightly different for the odd-parity equation. and a demonstration
in direction,
so
A more detailed derivation
that the operators are symmetric and positive definite are in [8,13].
2.2. Weak Form and Discretization
of the Even/Odd-Parity
Equations
The weak forms of equations (2.6) and (2.8) are derived by multiplying the equations by a test function u(r, a), integrating over the spatial and directional domains, applying the divergence theorem, and then applying the boundary conditions, equations (2.7) and (2.9). This procedure results in the even-parity weak form
JJ J! JJ
lu(r, Cl)GIlI@+(r,
+
Cl) dr dJZ +
u(r, O)@+(r, 52) [fl. nl dsdn
+
u(r, s2)Qf dr dfl
=
+J f
JJ JJ
u(r, O)G+@+(r, Cl) dr dC2
lu(r, s2)G11Q-(r, Cl) drdi2
u(r, a)$,
n.n$l
(r, rfl)
1fl. nl dsdC2
and similarly for the odd-parity equation
JJ Sf JJ
Iu(r, fi)G517+-(r,
+ +
JJ JJ J f'
Cl) drdC2+
‘u(r, C2)@-(r,Q) IC!. nl dsdi2
u(r, O)Q-(r,
f2) dr dC2 T-
=
Cl+0
u(r, Cl&_@-(r,
Cl) drdSl
lu(r, fi)GT’Q’(r,
u(r, a)@.6 (r, #I)
C2)dr dC2
Ifi. nl dsdfi.
In a discrete-ordinates formulation, the directional variable, 0 E S2, varies only in a finite set of discrete directions s2,, and the directional integrals are replaced by quadrature sums. This is equivalent to collocation in direction, with basis functions defined by w,6(1 - 0 Cl,), where 20, is the mth quadrature weight. The spatial domain is approximated by an unstructured mesh of triangles, and a set of basis functions Gn(r) is defined within each element. For a Galerkin FE treatment, the test functions are the same as the basis functions. To proceed, we write the discrete-ordinates form of the two scattering operators as
ij*jm(r) = cr,(r)jm(r) -
5
w,~o$ (r, %2
. %I
.fd(rIr
m’=l
where M is the number of discrete directions. operator is
The discrete-ordinates
form of the transport
The discrete function fm(r) = f(r, Cl,) is the function f(r, 0) evaluated at the mth quadrature direction.
C. R. DRUMM AND J. LORENZ
60
Within a spatial element e the approximate even- and odd-parity solutions in a discrete direction C&, are assumed to be combinations of the basis functions for r E V,,
e = 1: NE,
m = 1 : M,
where $&,,,, are the values of the solutions at the nodes of the element, N is the number of basis functions per element, NE is the number of spatial elements in the mesh, and V, is the element volume. The distributed
source is defined similarly by for r E V,,
e = 1: NE,
m = 1 : Ad,
and the discrete boundary conditions are for r E AX,
f12, . n 2 0,
where the boundary conditions are defined on the external boundary tW, of the elements e which have an external boundary. With these discrete functions and operators
defined, the FE discrete-ordinates
form of the
even-parity weak form is
forn=l:N,
e=l:NE,
CZ,.n20,
m=l:M.
The volume integrals are carried out over the element volumes, and the surface integrals are carried out over the external boundary of the domain. The odd-parity weak form can be discretized similarly
for n = 1: N,
e = 1: NE,
ST&,, ’ n p 0,
m = 1: M.
Equations (2.10) and (2.11) define linear systems of dimension N * M, where N is the number of FE nodes in the mesh, and the unknowns c#&,,, are the values of the even/odd-parity fluence at the nodes and discrete directions.
Parallel FE Approximation
In order to setup equations tions of the scattering represented constant
(2.10)
operators
and (2.11),
of the scattering scattering
expansions.
homogeneous
operators
operators
expressions
for the discrete
approxima-
c+ and g__, ancl their inverses are needed, and they are commonly
by finite Legendre-polynomial
within a spatially
explicit
61
is dropped.
Since
material
region of the problem
data
geometry,
The Legendre-polynomial
are assumed
the spatial
approximations
to be
dependence
of the discrete
are
and
where
L is the order of the Legendre
Legendre
polynomial
expansion
of the scattering
of degree 1. The inverse operators
cross section.
and Pl is the
are [12]
and
Based
on these approximations,
to solve it with a parallel article.
the computational
problems,
a code has been developed
CG algorithm results
are compared
and then the computational
more realistic
in the Aztec results
to set up the linear system,
code package.
with analytic are compared
and then
In the next section
solutions,
for several
of this
simple
test
data
for a
with some experimental
problem.
3. RESULTS In this section
the convergence
rate of the error
imations
of the code is investigated
solutions
are available,
aluminum.
If an exact
followed solution
due to the spatial
by analyzin g several
by a comparison is known,
a discrete
sirnple test
with measured ordinates
and directional problems,
data
where
for 1-MeV
approximation
approxanalytic
electrons
on
of the L” norm of
the error is defined by
E=
where f(r,
!&)
of these
results
nontrivial, problems
results
rate,
h is the maximum
are designed
element
for elliptic
point,
and the exact
is 2, and that
size, and c is a positive
boundary
of the even/odd-parity
and will be presented
basis functions
model is used.
i(r,%)/2 dr.
at the mu’ quadrature
are well established
to the discretization
however,
with linear
solution
-
solution
to obey the relation
r is the convergence
Such convergence
f(r,G)
w,
is the approximate
is .f(r, %J. The error is expected
where
c 1’ J..
value problems.
form of the Boltzmann
in future work. The expected
convergence
for quadratic
is 3.
to test the convergence
basis functions
of the spatial
approximation,
The
const,ant.
An extension equation
is
rate for FE first two test
so a simple physics
C. R. DRUMM AND J. LORENZ
62
3.1. 1D Test Problem Consider a unit interval in space, 0 5 z < 1, unit total cross section, ut = 1, scattering and absorption equally probable, and linearly-anisotropic
scattering, a,@~)
= (1/2n)(1/4+pc/4),
corresponds to Legendre moments, 0~0 = l/2, LT,~= l/6, where gsl = 2~s:~ With these assumptions,
where the coordinate
the first-order form of the Boltzmann
system is illustrated
which
as(po)Pl(~~O) dpo.
transport equation is
in Figure 2 and ~0 is given below. The 1D approxi-
mation assumes that the fluence, source, scattering cross section, and boundary conditions are independent of two of the spatial variables x and y ($$’ = $$ = 0), and also independent of the azimuthal direction (p, which is the angle between the positive x-axis and the projection of the particle direction on the zy-plane. These approximations result in an equation in two independent variables z and p, where p = cos(B), and 13is the angle between the particle direction and the positive z-axis.
Figure 2. Coordinates equation.
for 1D and 2D approximations
to the Boltzmann
The cosine of the angle between the pre- and post-scattered written in terms of the pre- and post-scattered directions as
directions
transport
is ~0, which can be
PO = pL/l’+ sin(e) sin(8’) cos (‘p - cp’) . Since the integral over the azimuthal angle cp’ of cos(cp - cp’) is zero, the transport simplifies to
equation
(3.2) In order to test the error in the numerical results, we specify a simple solution
@(z,p) = (1 + p)e-“. The source term and boundary conditions that result in the desired solution are
Q(GIL) =
(-
63
ParallelFE Approximation and
a@,P) =
1+
forO
Pu,
and
+b(l,p) Computational
= (1 +
for - 1 5 /I < 0.
/.+-‘I
results are determined by specifyin, 0 a uniform spatial mesh with linear or
quadratic basis functions, and a low-order Gauss quadrature for approximating the p-integrals. Because of the simple angular dependence of this example problem, a 2-point Gauss quadrature is sufficient to exactly model the angular dependence, so any error in the computational
results
arise from the spatial approximation. At this point, it is interesting to write out equations for the even/odd scat,tering operators and integro-differential equations for the even- and odd-parity solutions. The even-parity fluence is
@+(z,p) =
e-‘,
and the odd-parity fluence is
a-(z, p) = pee”. The even-parity source is
Q+(z, P) =
(-,u2+ ;) e-z,
and the odd-parity source is Q-(z,p)
= -Eeaz 6
’
The scattering operators are
G+f(z, p) = f(z,
CL) - ;
J:f(z1P’)dP’
and
G_f(z, p) = f(z, p) - f /_,
PP’f(Z,P’)dl.L’.
The inverse operators are 1
o,1.f(w4 = f(z,CL)+ Fj
Jf l
_l
(z, CL’) dP’
an d
The even-parity transport equation is _
p2
a2@+ (z,PL)-- l”. /; dZ2
p2pnazm;~;‘p’) dp’ + O+(z,p) =
-p&
Q-(z,
- ;
s_: @+ dP’
3 J _’p/IQ-
p) + 1o
1
(Z>P’)
(z, P’) dd
I
+ Q+(G ~1,
C. R. DRUMM AND J. LORENZ
64
with boundary
conditions
for p >< 0, z = 0, The odd-parity
transport
equation
_/~2a2@-(zd4 -- 1 l a9
2 J _r PP
~20,
z=l.
is
,a2a- (6 i-4 dp’ d.9
+ a-- (z, p) - ;
Q+k
= -P$
with boundary
and
CL)+ ;
J w'@-(~,l~') J1 1 _’
1
_’ Q+ k P’) @’
+ Q- (2, P)>
conditions
=
-@b
(2, F/J) + Q+(z, P) + ; for p 5 0, z = 0,
and
/’ Q+ (6 P’)dd, 1
~2O,z=l.
The results are given in Tables 1 and 2, which list the error norms vs. number of uniform elements. Also given in the table are the convergence rates, which are determined from a linear least-squares fit of the log of the error norm vs. the log of the mesh width. The expected convergence rates of 2 for linear basis functions and 3 for quadratic basis functions are obtained. The error-norm 3.2.
results
are plotted
in Figure
3.
2D Test Problem
Consider a unit square, 0 5 x < 1 and 0 < y 5 1, unit total cross section, crt = 1, scattering and absorption equally probable, linearly-anisotropic scattering, CT,(~0) = (1/27r) (l/4 + po/4), which corresponds to Legendre moments of us0 = l/2 and us1 = l/6. With the assumption transport that the solution is independent of Z, $$ = 0, the first-order form of the Boltzmann equation is
where the coordinate system is illustrated in Figure 2. The cosine of the angle between the pre- and post-scattered
directions
is ~0,
~0 = pL,fii + sin(d) sin (0’) cos (‘p - cp’), where CL2= cos(@, where 8 is the angle between the particle direction and the positive z-axis, and cp is the azimuthal angle between the positive x-axis and the projection of the particle direction on the xy plane.
65
Parallel FE Approximation Table 1. Error in the numerical functions.
results for the 1D test problem
for linear basis
L2 Norm of Error: Linear Basis Functions Number of Elements Even-Parity
I (
Fluence
Odd-Parity
Fluence
5
1.98 x 10-3
8.39 x 10-d
10
4.95 x 10-d
2.10 x 10-a
20
1.24 x 10-d
5.25 x IO-”
40
3.10 x 10-S
1.31 x 10-S
80
7.74 x 10-s
3.28 x 10-s
160
I
Convergence
Rate
-r~--
1.93 x 10-s
8.21 x 10-7
2.00
[
2.00
I
I
Table 2. Error in the numerical results for the 1D test problem for quadratic functions. L2 Norm of Error: Quadratic
Number of Elements
Even-Parity 5 10
I
3.78 x 10-G
I
20
I
4.73 x 10-7
-
1o-g
Fluence
1.74 x 10-S ~ 2.18 ~~ x 1O-6
_(
2.73 x 1O-7
I
40
5.91 x 10-s
3.41 x 10-s
80
7.38 x 10-s
4.27 x lo-”
Rate
3.00
3.00
I
I
I
I
i
I
I
I,
40 l/(maximum
-1 I
I
’ \w
quadratic interpolation functions convergence rate 3.00 I
Basis Functions
Odd-Parity
3.01 x 10-s
I
Convergence
E
Fluence
basis
element size)
Figure 3. Error analysis for 1D test problem
III,
80
/II,-
120 160
66
C. R. DRUhml
AND
J. LORENZ
If the solution is independent of azimuthal angle cp, the integral over the azimuthal angle cp of cos(cp - cp’) is zero, and the transport equation is
Assume a solution Q,(2, Y, k)
= cL%cmZ.
This corresponds to a distributed source term of
Q(x,Y,P~,P~)
=
P: (1 -PX) -
I
$ e-“,
and boundary condition %
(2,Y, CL,)= I-L%-x,
for (z, y) E IX
The even-parity fluence is @+ (.‘c,Y, L)
= &ePZ,
and the odd-parity fluence is a- (G y, PZ) = 0. The even-parity source is
Q+ (x,Y,P,,P~)
=
(112- i >emz,
and the odd-parity source is
With linear basis functions, Sis quadrature (72 directions [14]), the L2 error norms are listed in Table 3 and plotted in Figure 4, showing nearly the expected quadratic convergence rate. Table 3. functions.
Error
in the numerical
results
for the 2D test
problem
for linear
basis
Parallel FE Approximation
/
/
I
I
I
I I /
25 l/
/,I, 50
,,I, 75 100
(maximum element size)
Figure 4. Error analysis for 2D test problem
3.3 Convergence of Directional Approximation The third test problem is designed to test the convergence of the directional approximation. Consider a unit interval in space, 0 5 z 5 1, unit total cross section, CT t = 1, isotropic scattering, with scattering ratio c, which is defined as the ratio of the scattering cross section to t,he total cross section. The scattering transport equation is
ratio lies in the range, 0 5 c < 1. For this problem,
_P2 a2@+ (6 P) + @+(z,p) a9 with boundary
1
= ;
CD+(2,~‘) dp’ + Q+(z,p). s -1
conditions
==t@, (z, &cl) + Q-(2,/L). A test problem
the even-parity
is defined
as follows:
the even-parity
for St. n 2 0,
z = 0.1.
(3.3)
fluence is
@+(Z, /L) = COS(l_L), and t,he odd-parity
fluence is zero, a- (2. /L) = 0,
which corresponds
to source terms of &+(a, p) = “OS(~) - csin(1)
an d Q-(z,/L)
= 0.
The convergence of the directional approximation is examined by dividing the directional domain into a number of uniform intervals, and approximating the integral over each subregion by either a l-point or 2-point Gauss quadrature. Using one Gauss point per angle mesh should result in quadratic convergence, and using two Gauss points per angle mesh should result in fourthorder convergence. This is because one-point Gauss quadrature integrates up to linear terms exactly, and two-point Gauss quadrature integrates up to cubic terms exactly. Figure 5 shows the convergence of the L2 norm of the directional approximation. showing nearly the theoretical convergence rates. The results are tabulated in Table 4.
C.
68
R. DRUMM AND J. LORENZ
1C3
lo9
lOA
104
10"
10"
10"
IO8
lo-'
10.' I -
10” :
2 Gauss convergence
points per angle mesh rate
‘L
_
=
4.00 I
. I
5
1oe9
Number
of
Angle
Number of Angle Meshes
3.4. Comparison
2
3.98 x 1O-3
5.71 x 10-s
4
1.01 x 10-S
3.54 x 10-T
8
2.54 x IO-*
2.22 x 10-8
16
6.36 x 1O-5
1.38 x 1O-g
1.99
4.00
with Experimental
10”
rates for test problem 3. 2 Gauss Points
Rate
15
approximation.
1 Gauss Point
Convergence
I
Meshes
Figure 5. Error analysis for convergence of the directional Table 4. Error norms and convergence
10
10"
.
Data
Lockwood et al. [15] used a calorimetric technique to measure energy-deposition profiles in 1D for electron sources. In this section, we compare the results of one of their sets of measurements with FE modeling. The reported experimental uncertainties are less than a few percent. The problem considered here is that of l-MeV electrons normally incident on aluminum. The parameters for the FE modeling are: 40 uniform electron energy groups, no coupling with photons, Pi5 (15th-order) Legendre expansion of the scattering cross section directional dependence, Sic (16-point) Gauss quadrature, 50 uniform spatial elements, with quadratic basis functions. The comparison between calculation and experiment is shown in Figure 6, with excellent agreement . 3.5.
Parallel
Efficiency
In order to test the parallel efficiency of the code, a 2D test problem was run on Sandia’s ASCI red computer, utilizing from one to 512 CPU’s. The problem considered is the 2D benchmark problem described in Section 3.2. Two different meshes were considered, a mesh with a maximum element size of 0.01 with 10,201 FE nodes, and a coarser mesh with a maximum element size of 0.02 with 2,601 FE nodes. For the fine mesh, three different quadrature orders were considered, Ss, Sis, and Sic quadrature, with 20, 42, and 72 discrete directions, respectively [14]. For the coarser mesh, only Sls quadrature was considered. Figure 7 is a plot of the speedup versus the number of processors for the four problem, where the speedup S is defined by
Parallel 3.5
I’
4
i
‘I
Lockwood
69
FE Approximation
I
/
’
8
I”’
I
I
data
Fraction Pigllre 6. Energy deposition
of
mean
range
profiles for l-Me11 electrons
on aluminum.
500
450
+ __C__
400
h=O.Ol, h=O.Ol, h=O.Ol, h=0.02,
S,,, S,*, S,, !Z&
734,472 428,442 204,020 187,272
unknowns unknowns unknowns unknowns
n 350 z g 300 n cn 250 z? z2
200 150 100
200
300
Number of Processors Figure
where
7. Parallel
Ir, and T, are the runtimes
required
for the matrix
quite good, improving the larger
problems
speedup on teraflops
for 1 and n processors,
and right-hand as the number
have better
is less for the larger problems.
computer
for four matrix
respectively,
sizes.
which include
side setup and for the CG solve. The parallel of unknowns
parallel
efficiencies,
The convergence
increases.
For a given number
since the relative
criterion
fraction
the time
efficiency
is
of processors,
of communication
is that the relative error norm is reduced
C. R. DRUMM
70 to
10m8 of its original value, /jrll/llr~ll < lo-‘.
AND J. LORENZ
For the largest problem considered, speedup of 466
was attained for 512 processors, corresponding to better than 90% parallel efficiency. 3.6.
Preconditioning
The CG iterations converge fairly slowly when applied to the unpreconditioned but by preconditioning
the matrix, the convergence can be accelerated.
linear system,
Table 5 compares run
times and iteration counts for the unpreconditioned system compared to the results with several preconditioners that are included in the Aztec software. The problem considered is again the 2D benchmark problem described in Section 3.2. Only the largest problem was considered, with a maximum element size of 0.01, Sls quadrature 734,472
unknowns.
The calculations
with 72 discrete directions,
corresponding
were performed on Sandia’s ASCI red computer,
to
utilizing
100 CPU’S. Table 5. Comparison
of Aztec preconditioners. Aztec Solver
Preconditioner
Time (s) 468
2,879
point Jacobi
226
1,384
l-step block Jacobi
271
1,420
3-step block Jacobi
678
1,207
none
I
( l-order Neumann polynomial
)
626
1
1,933
3-order Neumann polynomial
a73
1,362
l-order least squares polynomial
667
2,062
I-3-order
least squares polynomial
710
I
(
1,107
ilnt(O.0)
439
754
ilut(0.02)
446
784
ilut(0.04)
910
1,721
Only two of the built-in preconditioners significantly reduced runtime, point Jacobi and l-step block Jacobi. Incomplete LU with threshold cutoff (ilut) significantly reduced the iteration count, but without much change in the runtime. The various preconditioners are described in the Aztec manual [2] and are briefly described here. Point Jacobi is equivalent to a symmetric diagonal scaling, resulting in all 1s on the diagonal. Block Jacobi uses a small number of block Jacobi iteration steps as a preconditioner, with the block size equal to the number of discrete directions in the quadrature integration. Neumann and least squares use a low-order polynomial approximation as a preconditioner, and ilut uses an incomplete LU factorization with a threshold cutoff as a preconditioner. Much work remains to be done to find a good preconditioner for this linear system. The convergence criterion is that the relative error norm is reduced to 10m8 of its original value Ilrjl/llroll< lops. 4.
CONCLUSIONS
A computer code has been written to solve the linear Boltzmann transport equation on an unstructured mesh of triangles, and the code has been applied to a number of benchmark problems to test the spatial and directional discretizations. The solver is based on a triangular mesh of a 2D region, but extension to 2D quadrangular or mixed triangular-quadrangular mesh or 3D tetrahedral or hexahedral mesh would be straightforward. The 2D FE approximation is currently based on linear basis functions, which are fairly inaccurate in regions with steep gradients
Parallel
in t,he solution, included
common
in radiation
in a 1D version
Because
the method advantage
from the scattering elements
source-iteration designed than
problem,
experience
for photon
by takin, c advantage
Work
is in progress
CG converges
algorithm.
for this
Interestingly,
term.
For this possibly
indicates
groups,
unlike
basis functions
simultaneously, of the structure
to exploit
submat,rices,
based
an effective
this
t,he standard
source
system
of t,he matrix,
of the submat~riccs
structure
substantially
by c*omputing
reducing
storage.
as does the 11naccPleratc:cl
preconditioner
on a coarse-mesh
that the CG iterations
the linear
of the spsrsity
slowly for large problems,
reason,
have been
in the 2D code would be straightforward. dependence
on the fly from a few precomputed
unaccelerated
Quadratic
is designed to take full advantage
can be gained
arising
71
problems.
of the code, and inclusion
matrix
Furthermore,
transport,
solves the space-direction
may be huge. The Aztec software but flIrther
FE Approximation
is neetletl
space or direction
specifically
ilpproximnt,ioll.
converge much faster for electron iterations.
This
behavior
groups
mc,rit,s further
investigation. ‘l%
co&+ is designed
t,ransport,
of charged
tions [ll]. sections
The reference [IO], resulting
port codes. the electron a continuous
to model
particles,
describes in electron
For some applications, cross sections,
the transport
such as electrons,
of neutral
a Goudsmit-Saunderson cross sections this method
which is inefficient.
slowing down term directly
particles
by usin g specially
modification
that, are compatible requires Possibly,
in the transport
but
modified
can
also model
electron
to CEPXS
a better equat,ion.
electron
with nerltral-pirrt,ic:le
a very high order Legendre approach
t,he
cross setcross traus-
expansion
of
would be to inclllde
and this work is in progrclss.
REFERENCES G. Bell and S. Glasstone, N?lclear Reactor Theory, Va11 Nostrand Reinhold. New York, (1970). S. Hutchinson, J. Shadid and R. Tuminaro, Az~.ec user’s guide version 1.1. Report SAND!%15.59. Sarldia National Laboratories, (1995). 3. E.E. Lewis and W.F. Miller, Jr., Computational Methods of Neutron Transport, Wiley, New York. (15~4). PA. (1997). 4. .J.W. Demmel, Applied Numerical Linear Algebra, SIAM, Philadelphia. Y. Saad, Iterative Methods for Sparse Linear Systems. ITP. Boston, MA. (1996). Parametric Technology Corporation, Pro/MESH User’s Cuzde, Release 17.0, (19%) R.T. Ackroyd, Finite Element Methods for Particle Transport, Research Studies Press, Somerset. (1997). 8. V. Agoshkov, Boundary Value Problems for Transport Equations, Birkhauser, Boston, hIA, (1998). trnnsljort 9. J.A. Josef and J.E. Morel. Simplified spherical harmonic method for coupled electron-photon calculations, Phy. Rev. E 57, 6161 (1998). finite-element, solut,ion of the neutron transport, quation in 10. T.A. Manteuffel and K.J. Ressel, Least-squares diffusive regimes, SIAAl J. Numer. Anal. 35, 806 (1998). C.R. Drumm. hlultidimensional electron-photon transport with standard discrete ordinates ctrtles, Nuc~. Il. Sc1. nnd Eng. 127 (1) (1997). and W. hlartin, Transport Theory, Wiley. New York. (1979). 12. J. Duderstadt transport analysis on 2-D unstrucl.llred mesh. llcport 13. C’ R. Drumm, Parallel finite element electron-photon SAND99-0098, Sandia National Laboratories, (1999). for nuclear analyses: Theory and guidelines for 14. R.D. O’Dell and R.E. Alcouffe, Transport calculations effective use of transport codes, Report LA-10983-MS, Los Alamos National Laboratory, (1987). measurement of electroll energy 15. G.J. Lockwood, L.E. Ruggles, G.H. Miller and .J.A. Halbleib. Calorimetric deposition in extended media-theory vs. experiment, Report SAND79-0414, Sandia Nat,ional Laborat,ories. (1979). A multigroup c~ouplrtl elect,ron 16. L.J. Lorence, Jr., .J.E. Morel and G.D. Valdez, Physics guide to CEPXS: photon cross-section generating code, Version 1.0; Report S.4ND89-1685, Sandia National Lal,oratories, (1989). 2.