Graziani F. (editor) Computational Methods in Transport

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

Domain Decomposed Particle Monte Carlo 427

Using zone position is safer than using global zone numbers, because a

code may not produce those. However, the position may not be invariant,

because the code producing the grid positions may give slightly diﬀerent (i.e.,

“jittery”) positions with diﬀerent numbers of domains. We will now describe

an algorithm that gives invariant seeds from zone positions, provided that

there is not too much jitter in these positions.

The ﬁrst step is to ﬁnd the minimum and maximum values, in each spatial

dimension, of the locations of the zone centers on each processor. That is,

we get the minimum and maximum values of x, y,andz on each domain.

Then we ﬁnd the minimum and maximum over all domains by using the MPI

Allreduce command. Since these global values may have some jitter, we shave

oﬀ the lower order bits. This is done by a “shaving” algorithm given in the

appendix. Now we have three double precision numbers for the grid that are

invariant over the number of domains.

Next, we loop over each zone in the grid and scale its position by the

minimum and maximum values for the grid. That is, we calculate (x

zone

−

min

)/(x

max

− x

min

), for x, y,andz. Because x

zone

may also contain some

jitter on diﬀerent numbers of processors, we apply the shaving algorithm to

these three numbers. This gives us three numbers in [0, 1] for each zone that

are invariant over the number of domains. At least one of these numbers

will be diﬀerent for each zone. (Equality could occur for one or two of the

numbers because, for example, zones in a one dimensional problem would

have the same x and y values if the grid only had extent in the z direction.)

Then, we multiply these three numbers by 1.8446744 · 10

≈ 2

,which

maps them into a 64 bit unsigned integer, and we add the three numbers

together. This yields a 64 bit integer number that is diﬀerent for each zone.

To reduce correlations between zones, we then change this value into a new

unique 64 bit unsigned integer by subjecting it to the DES hash algorithm

[PTVF02]. This yields one 64 bit unsigned integer for every zone that is

invariant over the number of domains. To get initial seeds for each particle

in the zone, we increment the zone value by one and apply the DES hashing.

This gives each particle a unique seed that is independent of the number

of domains. Thus each particle accesses the same pseudo-random number

stream, independent of the number of domains.

3 Ensuring That Addition is Commutative

Using the algorithm described above, we can ensure that all particles access

the same pseudo-random number stream independent of the order in which

they are simulated. They will still, however, deposit energy in the zones in

a diﬀerent order when the number of domains is changed. Because ﬂoating

point addition is not exactly commutative, there will small diﬀerences in the

total energy deposited in each zone at the end of the time step.

428 N.A. Gentile et al.

These diﬀerences are on the order of roundoﬀ. However, we cannot toler-

ate them because they will eventually have macroscopic consequences. This

is because the energy deposited in the zone will aﬀect the temperature and

opacity of the zone in the next time step. The opacity can be a nonlinear

function of the temperature, and so small diﬀerences in the temperature can

be magniﬁed. Diﬀerences in the value of the opacity will cause diﬀerences in

the deposition of energy by every photon that enters the zone. This in turn

will eﬀect the creation and behavior of particles in the zone. Eventually, some

particle will behave diﬀerently (e.g., not scatter) because of slight diﬀerences

in the values of temperature and opacity. This diﬀerent behavior will have a

macroscopic eﬀect on the problem, which will aﬀect other particles in subse-

quent time steps. Soon, the diﬀerence between the two cases will be as large

as if diﬀerent random number streams were used.

This eﬀect is illustrated in Fig. 2. This plot shows the diﬀerence in the

temperature in the ﬁrst zone of three diﬀerent simulations of a test problem

100

150

200

250

300

350

Time Step

-16

-14

-12

-10

-8

-6

-4

-2

Temperature Difference, |T

-T

| T

-T

| T

-T

Fig. 2. The absolute value of the diﬀerence in temperature in the ﬁrst zone of three

simulations vs. time step. Two simulations used the same pseudo-random number

stream but diﬀerent particle order. (These are denoted the “forward”, T

, and

“reverse”, T

, simulations.) The third simulation used a diﬀerent pseudo-random

number stream (the “alternate”, T

, simulation.) |T

−T

| is depicted with a dotted

line and |T

−T

| with a solid line. |T

−T

| is large at the beginning, and ﬂuctuates

but does not increase with time. |T

−T

| begins at a value near the roundoﬀ error

for double precision arithmetic, but grows and eventually jumps to same level as

the diﬀerence between the “forward” and “alternate” simulations

Domain Decomposed Particle Monte Carlo 429

from [FC71]. The opacity is given by (2) in [FC71]:

σ =

3/2

(1 − e

ν/T

)

(3)

with ν and T in keV.

The simulations used 100 zones with ∆x =0.4. It was run for 400 time

steps with ∆t ﬁxed at 2.0·10

−12

sec. A temperature source with T

=1.0keV

was applied at x = 0. The initial temperature was 0.01 keV. The equation

of state had a constant heat capacity of C

=8.11829 · 10

erg/(cm

keV) =

0.5917aT

,wherea is the radiation constant. Each simulation used 1000

particles in each time step.

All three simulations were run on one processor, and all used the pseudo-

random number described above. Two simulations employed the algorithm

described above which ensured that particles got the same pseudo-random

number stream in each case. The only diﬀerence between these simulations

was that the particles were run in a diﬀerent order. After new particles were

created at the beginning of each time step, the order of the list was reversed

before they were tracked. Thus any diﬀerence in the results of these two

simulations is due to diﬀerences in the order of ﬂoating point addition when

the particles deposit energy.

The third simulation used a diﬀerent pseudo-random number stream en-

tirely. This was accomplished by adding 78654092354 to the initial value of

s in each zone in the simulation.

Figure 2 plots |T

−T

| and |T

−T

| versus the time step in the simulation.

Here T

is the temperature in the ﬁrst zone when the particles are run in

the usual order, T

is the temperature in that zone when the particle order

was reversed, and T

is the temperature in the run which used a diﬀerent

pseudo-random number stream.

The diﬀerence between runs using diﬀerent pseudo-random number

streams is ﬂuctuates between 0.01 and 0.1 throughout the simulation. (The

value of the temperature in the ﬁrst zone quickly becomes approximately

equal to the T

=1.0.) The diﬀerence between runs using the same pseudo-

random number stream is initially very small. The diﬀerence in the ﬁrst time

step is O(10

−14

), which is the size of roundoﬀ errors. This causes a diﬀerence

in the opacity in the zone, which means that photons in the zone will deposit

slightly diﬀerent amounts of energy in the two simulations. The diﬀerence

in temperature grows with each time step. However, the diﬀerence remains

small, and global measures like the total number of particles simulated remain

the same until time step 304.

At time step 304, the accumulated diﬀerence in temperature is large

enough to change the large-scale behavior of the code. A diﬀerent number of

particles are created in the two simulations in this time step. After that, the

diﬀerences quickly grow until they are of the same order as the diﬀerence we

ﬁnd compared to the simulation with a completely diﬀerent pseudo-random

number stream.

430 N.A. Gentile et al.

The cure for this problem is relatively simple. Floating point addition is

not commutative, but integer addition is. We eliminate diﬀerences in results

by mapping the energy to a 64 bit unsigned integer before doing the addition.

This scaling is a simple multiplication that maps the range of photon

energy in the problem into the range that a 64 bit unsigned integer can

hold, which is [0, 2

− 1]. This number changes with each time step, and is

calculated at the beginning of each time step. The multiplier is

S =2

/(E

census

+ E

source

) . (4)

Here E

census

is the sum of the energy of photons present at the beginning of

the time step, and E

source

is the total amount of energy added to photons

at the beginning of the time step (e.g., aT

zone

∆t for the thermal radiation

emitted from a zone.) We have used 2

rather than 2

− 1, the maximum

value of an unsigned 64 bit integer, as a safety factor. Note that S is a double

precision value, not an integer.

To calculate the energy deposited into a zone, we sum the 64 bit unsigned

integers obtained from scaling the energy deposited by each photon:

int

photon

= integer(E

photon

× S +0.5) (5)

int

dep

= E

int

dep

+ E

int

photon

(6)

Here E

int

dep

is an unsigned 64 bit integer, and integer() represents a cast

from a double precision value to a 64 bit unsigned integer. The 0.5 is added to

round to the closest integer instead of simply truncating the fractional part

of the value.

At the end of the time step, when deposition is complete, we get the total

energy deposited in each zone by reversing the multiplication:

dep

= double



int

dep



(7)

Here, E

dep

is of type double precision, and double() represents a cast from a

64 bit unsigned integer to double precision.

4 Results

Here we demonstrate that the algorithm outlined above will give the same re-

sults for simulations using diﬀerent numbers of processors. The results shown

are for the analytic transport benchmark of Su and Olson [SO97].

This test problem has an initially cold slab with constant opacity which

is heated by an isotropic source of photons. Interactions with the radiation

ﬁeld heat the matter and cause it to eventually reach thermal equilibrium.

The results depicted in Fig. 3 are for the case κ

=1,κ

=0,atτ =0.3, in

the notation of [SO97].

Domain Decomposed Particle Monte Carlo 431

0 0.1 0.2 0.3 0.4

0.5 0.6

0.7 0.8 0.9 1

Normalized Distance

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

Normalized Temperature

mat

1 domain

rad

1 domain

mat

4 domains

rad

4 domains

Fig. 3. Matter and radiation temperature in the Su Olson test problem vs. spatial

coordinate. Solid lines depict the results of a simulation run on one domain. Symbols

depict the results for every tenth zone of a simulation run on four domains. The

four domain results are identical to the one domain results, even reproducing the

details of statistical ﬂuctuations

Figure 3 shows results for matter and radiation temperature for one and

four domain runs. Simulations using one domain are solid lines. Every tenth

point for the four domain runs is plotted with a symbol. The four domain

results exactly reproduce the one domain results. Even the statistical noise

in the Monte Carlo solution is reproduced.

As part of the KULL regression suite, this test is run weekly. The results

for the temperature in every zone for one, two, and four domain runs are

printed to sixteen places, and the ﬁles are compared to ensure that the results

are identical. This reproducibility has been demonstrated for larger problems,

up to 1024 processors.

5 Conclusions

We have described algorithms that can be used to make domain decomposed

Monte Carlo photonics code produce the same answer, bit for bit, on various

numbers of domains.

Two main issues are introduced when the number of domains can vary.

The ﬁrst is that the particles may not get the same pseudo-random number

432 N.A. Gentile et al.

stream. The second is that the order of operations can change, causing dif-

ferences in the results arising from the fact that ﬂoating point addition is not

commutative.

The ﬁrst issue is eliminated with by giving each particle its own random

number generator state, and seeding these states in a manner that is inde-

pendent of the number of domains. The second issue is eliminated by using

integers, which are commutative, to do addition.

We have demonstrated that we can achieve bit for bit agreement on prob-

lems when the number of domains varies from one to one thousand.

Appendix: Shave Algorithm

This algorithm takes a double precision number x and truncates all but the

highest N

base-ten digits. We use N

= 7 in our application.

shave( x, N

if x == 0.0:

return 0.0

Store the magnitude and sign of x.

sign

= sign(x)

abs

=abs(x)

Truncate base-ten exponent to get magnitude of x.

integer n = integer(log

abs

))

Scale x

abs

to be between 10

and 10

d+1

double s =10

n−N

double x

scaled

= x

abs

Shave oﬀ digits by casting to integer.

integer x

shaved

scaled

= integer(x

scaled

)

Restore correct magnitude and sign.

return s · x

sign

· x

shaved

scaled

References

[FC71] Fleck, Jr., J. A., Cummings, J. D. : An Implicit Monte Carlo Scheme for

Calculating time and frequency dependent nonlinear radiation trans-

port. J. Comp Phys., 31, 313–342 (1971)

[GKR98] Gentile, N. A., Keen, N., Rathkopf, J.: The KULL IMC Package, Tech.

Rep. UCRL-JC-132743, Lawrence Livermore National Laboratory, Liv-

ermore CA (1998)

[BM04] Brunner, T. A., Mehlhorn, T. A.: A User’s Guide to Radiation Trans-

port in ALEGRA-HEDPP, Tech. Rep. SAND-2004-5799, Sandia Na-

tional Laboratories, Albuquerque, NM (2004)

Domain Decomposed Particle Monte Carlo 433

[BUEG04] Brunner, T. A., Urbatsch, T. J., Evans, T. M., Gentile, N. A., Com-

parison of Four Parallel Algorithms for Domain Decomposed Implict

Monte Carlo, Tech. Rep. SAND-2004-6694J (2004) Submitted to J.

Comp. Phys.

[SPRNG] http://sprng.cs.fsu.edu/

[PTVF02] Press, W. Tuekolsky, S.A., Vetterling, W. T., Flannery, B. P.: Numer-

ical Recipies in C++. Cambridge University Press, United Kingdom

(2002)

[SO97] Su, B., Olson, G. L.: An analytical benchmark for non-equilibrium ra-

diative transfer in an isotropically scattering medium. Ann. Nucl. En-

ergy, 24, No. 13 (1035–1055 (1997)

KM-Method of Iteration Convergence

Acceleration for Solving a 2D Time-Dependent

Multiple-Group Transport Equation

and its Modiﬁcations

A.V. Gichuk, L.P. Fedotova, R.M. Shagaliev

One of the main diﬃculties during ﬁnite-diﬀerence simulation of spectral

2D problems of particle transport and interaction with medium is to ﬁnd

cost-eﬃcient solutions to rather large systems of interconnected diﬀerence

equations.

The approach based on the method of simple iterations combined with

various acceleration algorithms is very often used to solve numerically both

single-group and multiple-group diﬀerence transport equations. Note that

conceptually many various acceleration methods are close to the method

based on introduction of an approximate operator (which is easy-to-use) along

with an operator of the original divergent equation. Following this path, a

number of eﬃcient methods for solving 2D transport problems have been

oﬀered.

KM-method is one of them. Finite-diﬀerence approximation using the

KM-method has a number of features owing to which it becomes highly

eﬀective in practice. First of all, one of such features is a combination of

explicit calculation of a collision integral with stability and conservative-

ness, this provides simultaneous satisfaction of the two balance correlations

typical for transport equation: for particle transport and for the number of

particles resulted from one collision event. The scheme stability in time has

been proved by analytical studies of some partial cases and the experience of

using KM-method in practice.

The approaches used to construct the KM-method found their further

development in MKM- and KM3-methods.

1 Statement of a 2D Transport Problem

In the kinetic multiple-group model describing the transport processes a 2D

time-dependent transport equation written in the classic divergent form looks

like, as follows:

∂ε

∂t

+ Lε

+ α

2π



j=1

(0)

,i=1,...,N , (1)

436 A.V. Gichuk et al.

Fig. 1. The coordinate system

Lε

= µ

∂ε

∂z

∂

∂r



1 − µ

cos ϕε



−

∂

∂ϕ



1 − µ

sin ϕε



. (2)

The equation has to be solved in axially symmetric region D =

(r, z) ∈

where

is the cross section of the body of revolution by a plane coming across

Z-axis (Fig. 1).

The following notations are used to write equations:

r, z are cylindrical coordinates of a particle;

→

Ω

(µ, ϕ)–is a unit vector in the particle ﬁght direction;

µ =cos(θ) , −1 ≤ µ ≤ 1;

θ is an angle between vector

→

Ω

and Z axis of symmetry;

ϕ is an angle between vector

→

Ω

projection to the plane coming through

point (r, z)atrightangletoZ and the vector connecting points (0,z)

and (r, z), 0 ≤ ϕ ≤ π;

N is the number of energy groups;

= ε

(t, r, z, µ, ϕ) is the mean ﬂux of particles in group “i” in the direction

determined by vector

→

Ω

;

(0)



−1

dµdϕ is the particle ﬂux density in group “i”;

= α

(r, z) is full cross-section in group “i”;

= β

(r, z) is secondary particle array.

System (1) is supplemented with initial and boundary conditions.

Finite-diﬀerence approximation to the system of (1) in space and angular

variables should meet the following requirements. First of all, it should be

KM-Method of Iteration Convergence Acceleration 437

conservative relative to collision-less particle transport with a possibility of

consistently approximating the collision term and the right-hand side of the

system of group kinetic equations. Secondary, it should be possible to reduce

the diﬀerence operator to the triangular form. In case of rectangular spatial

grids, DS

-type schemes meet the requirements above. The particular ap-

proximation form is insigniﬁcant for further description and, for this reason,

we only give the equation written in the form (1), where operator L means

a diﬀerence transport operator and function ε

and equation factors are the

grid values introduced in accordance with the diﬀerence template in use.

2 KM-method

The KM scheme is an iterative two-phase scheme. During the ﬁrst phase

(predictor) the diﬀerence equation is solved using the sweep (point-to-point)

computation algorithm with the known right-hand side. During the second

phase (corrector), the solution found is corrected so as to provide satisfaction

of balance correlations for particle transport and particle interaction. The

values obtained are put into the right-hand sides of the ﬁrst phase equations

and the process is repeated up to the desired precision achieved.

To approximate the original system of equations in time, write the KM-

method in the following form:

Phase I

s+1/2

n+1

−ε

∆t

+ L

s+1/2

n+γ

+ α

s+1/2

n+γ

2π



j=1

(0)

s+1/2

n+γ

= γ

s+1/2

n+1

+(1−γ) ε

,i=1,...,N .

(3)

Phase II

s+1

∆ε

n+1

∆t

+ L

s+1

∆ε

n+γ

+ α

s+1

∆ε

n+γ

−



j=1

s+1

∆ε

n+γ

2π



j=1



s+1/2

(0)

−

(0)



s+1

∆ε

n+γ

s+1

n+γ

−

s+1/2

n+γ

s+1

∆ε

n+γ

= γ

s+1

∆ε

n+1

,i=1,...,N .

(4)

Here, γ is weighting parameter of time approximation, n is a time step

number, s is an iteration number. In the scheme above, the second phase is

written in the corrector form. The (4) also allows a record relative to the

mean ﬂux function value, however, such a variant requires storing a large

amount of interim data.

Adding the (4) to the (3) and summarizing the results for the angular grid

directions we obtain that in the grid solution KM-method strictly follows

the correlation between the number of collided particles per unit volume