PARALLEL COMPUTING OF THREE-DIMENSIONAL MONTE CARLO SIMULATION OF TRANSIENT RADIATIVE TRANSFER IN PARTICIPATING MEDIA

perchorangeΛογισμικό & κατασκευή λογ/κού

1 Δεκ 2013 (πριν από 3 χρόνια και 6 μήνες)

70 εμφανίσεις

1

American Institute of Aeronautics and Astronautics

PARALLEL COMPUTING O
F THREE
-
DIMENSIONAL MONTE CA
RLO
SIMULATION OF TRANSI
ENT RADIATIVE TRANSF
ER IN PARTICIPATING
MEDIA




A. Sawetprawichkul, P.
-
f. Hsu, and K. Mitra

Mechanical and Aerospace Engineering Department

Florida Institute of Technology

Melbourne,
FL 32901




ABSTRACT

A three
-
dimensional transient radiative transport
model using Monte Carlo (MC) method was
implemented in a Beowulf
-
class parallel computer
system. The main purpose of this study is to simulate
the light pulse transport inside the ab
sorbing and
scattering media such as biological tissues. The three
-
dimensional results for the case of infinite height and
width but finite thickness configuration have good
agreement with our prior one
-
dimensional MC results.
The major advantage of MC m
ethod is its flexibility
and simplicity to simulate the photon movement in
arbitrary geometry and complex boundary condition.
Since the error bound of MC method is inversely
proportional to the square root of the number of
statistical samplings, it requir
es a large number of
samples to reach the satisfactory accuracy. Therefore,
the primary drawback of MC methods is that it is a
computationally intensive method. However, the
method is very adaptable to parallel computing, i.e.,
coarse grained algorithm.

Parallel computing is
introduced to improve the performance of this method.
The parallel system is based on a cluster of
commodity
-
class processors with a standard Massage
-
Passing
-
Interface (MPI). MPI is
a library of
subprograms. It

is chosen as a tool

to implement
parallel programming for its portability to other similar
systems and the very low cost compared to the
conventional supercomputer systems.


NOMENCLATURE

a

absorption coefficient, 1/m


c

propagation speed of radiation transport in
the medium,

m/s

G

incident radiation or integrated intensity,

W/m
2

g

asymmetric factor

I

radiation intensity, W/m
2
sr

L

medium thickness, m

q

radiative flux in
x

direction, W/m
2

ˆ
s

unit vector along a given direction

t

time, s

t*

dimensionless time,

t
ct

t
p

pulse width, s

x,y,z

position coordinate, m

x*

dimensionless distance,

t
x


Greek symbols




polar angle, rad


t

extinction coefficient,
1/m



circumferential or azimuthal angle, rad


s

scattering coefficient, 1/m




optical th
ickness of the medium



scattering phase function



solid angle, sr




scattering albedo


Superscript

*

dimensionless quantity


Subscript

p

pulse


INTRODUCTION

The transient radiative transfer problems are difficult
to solve analytically. Traditional radi
ation studies mostly
neglect the transient term of radiative transfer equation
(RTE).
1,2

The transient radiation, however, must be
considered in new applications such as the pulsed laser
interaction within the media.
3
-
6

Several numerical
techniques have
been studied and used to model the light
transport in participating media. The model is related to
the laser radiation in tissue for either therapeutic or
diagnostic purposes. The simulated diffuse reflectance
and transmittance are used to distinguish th
e tissues with
cancer from normal tissues.
7


The Monte Carlo algorithm was first introduced into
the laser tissue interaction applications by Wilson and
Adam.
8

Many researchers have been developing the
model to simulate the light distributions in tissue.
Hasegawa, et al.
9
as well as Brewster and Yamada
6
. The
2

American Institute of Aeronautics and Astronautics

anisotropic scattering results were predicted based on
the results of equivalent isotropic scattering media
using isotropic scaling law.
10

However, the results
agreed for thick medium only.

Most of t
he previous works were focused on
steady state solutions. Recently, some numerical
analysis of the transient radiative transfer has been
studied. The results of pulse radiation source using
different models have been reported.
11

The conclusion
is that the

transient radiative signature is quite different
among different models. Kienle, et al.
12

studied the
absorption and reduced scattering coefficients of two
-
layered turbid media by inverting the time
-
resolved
reflectance. Their results, however, were limit
ed to
several constraints such as requirement of a known
value of the first layer thickness and a close
approximation of the property data. Integral equation
formulations were recently developed by Tan and
Hsu
13

and Wu
14

for transient problems. The result
s
from integral equation formulation compare very well
with the previous study of Monte Carlo simulation.
15

The latter algorithm is based on that by Fleck.
16

The algorithm used in this study is verified with
other methods such as the discrete ordinate met
hod
based on a high order up
-
wind scheme
17

and Volterra
integral equation formulation
13
. By nature of the MC
method, the accuracy can be improved by increasing
the number of samplings. However, the drawback is
that the computational time increases proporti
onally.
Advancement in modern computers and parallel
computing is a solution for this problem. The other
possible alternatives are the reversed MC method and
quasi
-
MC method. However, these are subjected to
different limitations.

Parallel computing has b
ecome a very attractive
approach of the need for more computational power in
many scientific and engineering applications.
18

Weather forecast and computer animation are good
examples for the need for parallel computing. The
parallel computer architecture

used in this study is
based on the Beowulf concept developed in
NASA/Goddard
(Stirling, 1994
). It is simply a
collection of computers with commodity processors
and parts working together to solve a problem
efficiently and in less time and less cost.
The

coding is
based on Single Program Multiple Data model
(SPMD) and Message Passing Interface (MPI) is one
of the two commonly used standard libraries to achieve
parallelization. The system consists of one root or
head node and the other computers are consi
dered as
slave or compute nodes. The communication among
nodes, in our case, is through a private, fast Ethernet.
Details of parallel computing will be discussed in the
following section. Performance and efficiency of the
algorithm are measured.

The goal

of this study is to understand the transient
radiation in multi
-
dimensional participating media as
well as implementing parallel computing to improve the
performance and accuracy of the model. In the near
future, the non
-
uniform property distribution wil
l be
considered and the model predictions will be compared
with experimental data.


The Radiative Transfer Equation

The transient radiative transfer equation in direction
s
ˆ

is given as
19
-
21


































4

ˆ
,
ˆ

,
ˆ
,
4
)
,
ˆ
,
(
,
ˆ
,
t
,
ˆ
,
,
ˆ
,
d
s
s
t
s
x
I
t
s
x
b
aI
t
s
x
I
x
t
s
x
I
t
c
t
s
x
I

(1)


The RTE is an integro
-
differential equation. The first
term on the right hand side is the radiation loss due to
absorption and scattering. The medium is assumed to be
cold and its temperature distribution is not responsive to
the pulse la
ser incidence. Hence the emission gain (the
second term on the right hand side) may be neglected.
Moreover, the radiative properties of the medium and the
boundary are assumed to be constant through the
transient process. The last term on the right hand si
de is
the radiation gain by in
-
scattering. Equation (1) is a
mathematical model of transient radiative transfer
problems. Traditional approach is to solve the equation
either analytically or numerically for the solutions. The
Monte Carlo algorithm is diffe
rent from traditional
numerical methods such that it incorporates the
probability distribution and physical laws to simulate the
physical movement of the energy bundles and eventually
to solve for solutions. Next section discusses the
fundamental of the Mo
nte Carlo algorithm.


THE MONTE CARLO ALGO
RITHM

The Monte Carlo simulation is a stochastic method
that is based on the random sampling of variables from
probability distribution. The history of the energy bundle
is tracked as it propagates within the mediu
m. Each
energy bundle is assigned a fixed weight initially and
traced along the moving path. The bundle is either
absorbed or scattered during its interaction within the
medium. The history of bundles is recorded until either it
escapes the boundary or is
totally absorbed within the
medium.

The bundle is emitted into the medium through the
surface. The traveling distance to move the bundle before
any interaction (absorption or scattering) is then
calculated. After the bundle is moved according to the
initi
al direction, the new location is determined. If the
new location is still within the medium, it can be either
absorbed or scattered. Otherwise, the reflectance and
3

American Institute of Aeronautics and Astronautics

transmittance is updated when the bundle reaches the
surface. In case of scattering, a new
direction is
calculated as well as the new traveling distance. The
process is repeated until the last bundle has reached
the surface or being absorbed.

Followings are details of each calculation steps.
Sampling random variables is the basis of the
algorit
hm. The sampling method is discussed in details
and the mapping technique is introduced.


Sampling random variables

The principle of Monte Carlo method is based on
the random sampling of variables from probability
distribution. A probability density func
tion
p(x)

defines the distribution of random variable
x

over the
interval (
a
,
b
). This variable
x

may be the variable step
size, the angle of deflection, or the azimuth angle
during scattering. The simulation of light propagation
requires a value for
x
ra
ndomly and repeatedly based
on pseudo
-
random numbers, generated by a computer.
The random number,

, is uniformly distributed over
the interval (0,1). Though the variable
x
may be not
uniformly distributed over (
a
,
b
). The non
-
uniform
probability density f
unction
p(x)

can be expressed as
23


)
1
,
0
(
)
(





for
dx
x
p
x
a

(2)

This mapping technique will be used to represent
a non
-
uniformly distributed variable, such as the step
size or the deflection angle, by uniformly distributed
random numbers.


Initialize bundl
e

The incident beam can be collimated or diffuse.
For simplicity, the initial direction is assumed to be
normal to the surface (orthogonal). The weight is
initialized to unity. The effect of a refractive
-
index
-
mismatched interface between the medium and th
e
ambient is neglected in this study.


Generating the step size

The step size or the
bundle
’s moving distance
before next interaction,
s
, has the interval (0,

). The
probability density function of step size is the ratio
between the power extinguished in
distance
ds

and the
total power extinguished. It is proportional to
)
(
s
t
e




as
24
:


ds
e
ds
e
I
ds
e
I
ds
s
P
s
t
s
o
s
o
t
t
t
t












0
)
(

(3)


Applying the mapping technique yields:



t
s


)
1
ln(




or
t
s


ln



(4)



Moving the bundles

A
bundle

mo
vement can be represented by three
spatial coordinates (
x,y,z)

and two directional angles
(



). However, the deflection angle


is with respect to
the previous trajectory, which is arbitrary. It is more
convenient to describe the
bundle

position by three
Cartesian coordinates (
x,y,z
) and three directional cosines
(

x
,

y
,

z
)

A new position is updated by


i
x
s
x
x





i
y
s
y
y




(5)


i
z
s
z
z





Bundle absorption

There are two ways to determine bundle absorption.
First,
the bundle weight is constant until being absorbed
all in once. Another way is to split the bundle weight into
parts (absorbed and scattered). The absorbed fraction is
1
-

, where


is the albedo. The remained bundle is still
scattering until its weight is
below the threshold.

In this study, the first technique is used. The albedo
is compared with the random numbers to determine the
absorption.

If random number(

) >




absorbed

Otherwise




scattered


Bundle scattering

After the bundle survives
the absorption, it is ready
to be scattered. There

are two angles involving in
scattering process: a deflection angle,

)
0
(




, and
an azimuth angle,

)
2
0
(




. After a bundle is
scattered, its trajectory is deflected from t
he previous
direction by a deflection angle,



The scattering phase function describes the
probability distribution for the cosine of the deflection
angle, cos

.

The Henyey
-
Greenstein phase function was
originally proposed for galactic scattering approxi
mates
Mie scattering by particles comparable in size to the
wavelengths of light.
25

The scattering phase function is
defined as
21






2
3
2
2
ˆ
'
ˆ
2
1
1
)
ˆ
,
'
ˆ
(
s
s
g
g
g
s
s







(6)


The
asymmetric
factor,
g
, varying between

1 and 1,
characterizes the angular s
cattering distribution. Isotropic
scattering has a value of
g

as zero. Highly forward and
highly backward scattering has a value of g near 1 and

1
respectively. The same phase function is also described
by Modest
21
, but the product of two unit vectors is
used
instead of the cosine of the polar angle. The latter
expression is not a correct HG phase function.

The cosine of the deflection angle, cos

, can be
expressed in term of the random number,

. The
4

American Institute of Aeronautics and Astronautics

derivation of anisotropic scattering by the mapping
te
chnique was discussed in details in previous study.
26


In case of isotropic scattering (
g
=0), the
probability density function becomes (½)sin

.

The
mapping technique yields the cosine angle as



)
2
1
(
cos
1






(7)


The azimuth angle (

) is assume
d to be uniformly
distributed within the interval [0,2

]
.

The probability
density function of the azimuth angle
, p
(

)
,

is constant
and equals to

2
/
1
. Again, the mapping technique
yields the azimuth angle as:





2


(8
)


Once both the deflection and the azimuth angles
are found, the new direction may be calculated using
coordinate transformation. Cashwell
23

explained in
details the transformation method using complex
variables. The same results can be derived by using

rotation matrix as well. The new directional cosine
(

x
',

y
',

z
'
) is calculated from the current direction
(

x
,

y
,

z
) as:
































cos
1
cos
sin
cos
sin
cos
1
sin
cos
sin
cos
1
sin
2
'
2
'
2
'
z
z
z
y
x
z
y
z
y
x
y
z
x
z
x














(10)


However, if the bundle direction is close to the z
-
axis, the above equations are no longer valid sin
ce

z
~1. The following equations should be used instead:












cos
sin
sin
cos
sin
'
'
'
z
z
y
x




(11)


PARALLEL COMPUTING

Due to the nature of the Monte Carlo simulation, a
stochastic method, it requires a number of samples to
satisfy the required accuracy. Increasing
the samples or
bundles can improve the accuracy, but also requires
more time to compute. The previous MC work by the
authors were conducted on a 633 MHz 21164A Alpha
processor workstation (Sawetprawichkul, et al. 2000).
It could take several hours to com
plete the calculation
by the workstation. Recently, a 48
-
node IBM PC
cluster based was installed and used in this study. This
system is capable to extend to 96
-
node. Each node is
an IBM
x330 e
-
series, Pentium III 866 MHz processor
with 512 MB SDRAM. Th
e total system memory space
is 24 GB. The system is designed with one head node
and 47 compute nodes. The interconnect between the
nodes are channel
-
bonded fast Ethernet, i.e., double
bandwidth of 200 Mbps. The job queuing is provided by
the portable ba
tch system (OpenPBS), which is built on
top of the Maui job scheduler. The serial MC code was
modified to run on parallel processors by implementing
Message Passing Interface (MPI). Detailed hardware
and software configuration information is given in
web
site of http://olin.fit.edu/beowulf/.


Message Passing Interface

Message Passing Interface is a library of
subprograms in C, C++ or FORTRAN language.
18

The
major advantage of MPI is its portability in parallel
programming. The basic of MPI is to transmit

the data
from one processor to another. There are two types of
MPI communication: point
-
to
-
point and collective
communication. Point
-
to
-
point communication is
between one processor to another. Collective
communication, however, can send or receive the da
ta
from/to one processor to/from many processors.

In this study, we use Single Program Multiple Data
model, which means all the nodes will use the same
program but may produce different data depending on
the given input. Pseudo code of the Monte Carlo
alg
orithm under the MPI environment is given below:



Start the MPI function call by including the processor
directive and determining how many processors
involved and the rank of each processor.



If the process is root node (rank = 0)

-

Read the input data (suc
h as medium thickness,
scattering and absorption coefficient, total number
of bundles and time steps).

-

Broadcast these values to all processors to be used
for calculation.



Each processor performs the calculation based on
different random number sets gene
rated according to
the rank of each processor.



After all processors finish calculation, the output is
gathered/collected by root node.



Print the output from root node.



Terminate the execution by shutting down MPI call.

In parallel programming, it is also i
mportant to
clarify the scope of global and local variables. Since all
processors will perform the computation and
communicate with the root or node zero to share data.
The quality of parallel algorithms can be measured by
speedup and efficiency.

Quinn
36

gave the definition of the speedup as “
ratio
between the time needed for the most efficient sequential
algorithm to perform a computation and the time needed
to perform the same computation on a machine
5

American Institute of Aeronautics and Astronautics

incorporating pipelining and/or parallelism

.

The
e
fficiency of the algorithm using
p

processors is the
speedup divided by
p
. Pacheco
18

also defined the
speedup as the ratio between runtime of a serial
solution and runtime of the parallel solution with
p

processes.

The MC algorithm speedup typically scale
s up
very well with the number of processor, i.e., linear.
For most other parallel algorithms speedup, however,
have lesser value due to overhead from communication
time.


RESULTS AND DISCUSSI
ON

In previous study
37
, the one
-
dimensional Monte
Carlo simulat
ion, based on the work by Hsu
16
, has
been verified with the discrete ordinates method.
17

The simulation compares the three linearly anisotropic
scattering phase functions, i.e., forward, isotropic, and
backward; respectively (Fig. 1). The media has the
o
ptical thickness of


= 10, albedo of


= 0.998, and
an isotropic scattering phase function. The simulation
involved several hours run time in workstation due to
highly scattering media. Figure 1 shows a good
agreement of transient solutions between the one
-
dimensional Mont
e Carlo and the discrete ordinates
method for all three cases.

In this study, the MC algorithm discussed on the
previous section was developed to treat three
-
dimensional (3D) rectangular geometry. In the
following 3D cases, the optical depth is defined t
o be
the product of extinction coefficient and the plate
thickness. The 3D rectangular medium has equal
width and height to simplify analysis and the aspect
ratio (AR) is defined to be the length ratio of width to
thickness. The light irradiated at
z

= 0

plane or a node
on the plane. For the case of AR = 1000, the 3D
results have been validated with one
-
dimensional MC
results as shown in Figure 2. Both reflectance and
transmittance generated by 3D model agree very well
with the 1D model. For long time
solutions, both
reflectance and transmittance merge to the same value
as in the steady state solutions. The number of total
energy bundles used in the simulation is 10
7
.

Figure 3 shows the convergence rate of the
solution as the number of energy bundles i
ncreases. In
the calculation of the RMS error, the MC solution with
10
8

energy bundles is used as the basis. The relative
RMS error versus bundle number reflects the error
bound of
-
1/2 slope as predicted by the central limit
theorem.

A parametric study
of different aspect ratios is
shown in Figs. 4 and 5. The 3D model with uniformly
distributed source over the entire x
-
y plane was used to
compare with the 1D results. Each case emitted ten
millions energy bundles in total. It clearly shows that
the as
the aspect ratio increases, the transient signals are
closer to those of the 1D geometry. However, it is
surprising to find that one
-
dimensionality can only be
approached at a very large AR = 1000. This is quite
different from the steady state radiative
transport, in
which the AR = 10 case will produce close to one
-
dimensional result.

The calculated reflectance and transmittance at
different locations, center of the plane, edge and corner,
are presented in Figs. 6 and 7 for the case of AR = 10.
Even thou
gh the source is uniformly distributed, but the
location of exiting bundles affects the signal magnitude
due to out
-
scattering near the boundary. At the corner,
the transmittance and reflectance is the smallest of the
three positions. Similar results are

found for AR = 100
and 1000. At lower AR, the three
-
dimensional effect is
most evident.

The case of pulse irradiation at a single node in the
center of the
z

= 0 plane is studied. The temporal
distributions of reflectance and transmittance along the
x
-
a
xis are given in Figs. 8 and 9. Reflectance gradually
decreases as the time increases. However, the
transmittance value increases from
t
* =15 to 30, then
starts to decrease. The temporal trend is similar to those
shown in Figs. 6 and 7.

Even with 10
9

bu
ndles used in the simulation,
certain degree of fluctuation is still observed in the
signals near the edge. Physically, very few photons will
travel to the regions near the edges and corners, which
leads to weak signal magnitude. It is difficulty to obta
in
adequate sampling unless the total bundle number is
further increased. An attempt of 10
10

bundles was used
but didn't reduce the fluctuation much.

Speedup is calculated to measure the efficiency of
the algorithm. Figure 10 shows the speedups of the id
eal
case and actual performance at different processor
number. The parallel efficiency (speedup divided by
processor number) of actual calculation decreases
slightly at large number of processor. At 40 processors,
the parallel efficiency is 95%. It is s
uspected that the
communication overhead increases as more processors
are involved, thus the reduced efficiency.

Typical computing time, for example, the case of
Fig. 8, is 17 minutes with 40 processors for 10
9

energy
bundle simulations. The computational

time is
proportional to the number of bundle at a fixed number
of processor.


CONCLUSION

The three
-
dimensional Monte Carlo modeling has
been developed to simulate the transient radiative
transfer in participating media. The results were
validated with th
e one
-
dimensional Monte Carlo
solutions for the case of infinite height and width.
Parallel computing was been implemented to improve the
6

American Institute of Aeronautics and Astronautics

efficiency of the algorithm. It is found that the three
dimensional effect is very evident even with large
aspect rati
o. It is only at aspect ratio of 1000, the
three
-
dimensional results are close to that of one
-
dimensional code. Nearly linear speedup was
achieved with the algorithm except at large number of
processors. At 40 processors, the parallel efficiency
drops t
o 95% due to communication overhead.


ACKNOWLEDGEMENTS

This work is supported by Sandia National
Laboratories with contract AW
-
9963 and Dr. Shawn P.
Burns is the program manager of this project. The
parallel system is funded by National Science
Foundation

MRI program grant EIA
-
0079710 and Dr.
Rita V. Rodriguez is the program director.



REFERENCES


Sterling, T., D. J. Becker, D. Savarese, J. E.
Dorband, U. A. Ranawak, C. V. Packer, "Beowulf: A
Parallel Workstaion for Scientific Computation,"
Proceedings, I
nt. Conf. on Parallel Processing (1995).



1.

Fiveland, W. A. “Three
-
Dimensional Radiative

Heat Transfer Solution by the Discrete
-
Ordinates
Method.”

J. Thermophysics. Heat Transfer
. 2.4
(1988): 309
-
316.

2.

Raithby, G. D., and E. H. Chui. “A Finite volume
method for predicting a radiative heat transfer in
enclosures with participating media.”
J. Heat
Transfer

112 (1990): 415
-
423.

3.

Mitra, K., and S. Kumar. “Application of
Transient Radiative Transfer Equation to
Oceanographic Lidar.”
ASME HTD

353 (1997):
3
59
-
365.

4.

Gemert, M. J. C., an A. J. Welch. “Clinical Use of
laser
-
Tissue Interaction.”
IEEE Eng. Med. Biol.
Mag.

(1989): 10
-
13.

5.

Ishimaru, A. “Diffusion of Light in Turbid
Material.”
Applied Optics

28.12 (1989): 2210
-
2215.

6.

Brewster, M. Q., and Y. Ya
mada. “Optical
Properties of Thick, Turbid Media from
Picosecond Time
-
Resolved Light Scattering
Measurements.”
Int. J. Heat and Mass transfer

38
(1995): 2569
-
2581.

7.

Jacques, S. L., L. Wong, and A. H. Hielscher.
“Time
-
Resolved Photon Propagation in Tissue
s.”
Optical
-
Thermal Response of Laser Irradiated
Tissue.

Ed. Welch, A. J., and Martin J.C. van
gemert. New York: Plenum Press, 1995.

8.

Wilson, B. C., and G. Adam. “A Monte Carlo Model
for the Absorption and Flux Distributions of Light in
Tissue.”
Med. Phy
s.

10 (1983): 824
-
830.

9.

Hasegawa, Y., et al. “Monte Carlo Simulation of
Light Transmission through Living Tissues.”
Appl.
Opt.

30 (1991): 4515
-
4520. (all author's names)

10.

Guo, Z., S. Kumar, and S. Maruyama. “Time
-
Resolved Evaluation of Two
-
Dimensiona
l Scaled
Isotropic Results for Transient Radiative Transfer in
Anisotropic Scattering Media.” The 34
th

National
Heat Transfer Conf., Pittsburgh, PA, August 2000.

11.

Mitra, K., and S. Kumar. “Development and
Comparison of Models for Light
-
Pulse Transport
T
hrough Scattering
-
Absorbing Media.”
Applied
Optics

38.1 (1999): 188
-
196.

12.

Kienle, A., et al
. “Investigation of two
-
layered
Turbid Media with Time
-
Resolved Reflectance.”
Appl. Opt.

37.28 (1998): 6852
-
6862.

13.

Tan, Z. M., and P.
-
f. Hsu. “An Integral Form
ulation
of Transient Radiative Transfer.”
ASME J. Heat
transfer.

123(3): pp.466
-
475 (2000).

14.

Wu, C.
-
Y. “Propagation of scattered radiation in a
participating planar medium with pulse irradiation.”
J. Quant. Spectrosc. Radiat. Transfer
, 64.5 (1999):
537
-
548.

15.

Hsu, P.
-
f. “Effects of Multiple Scattering and
Reflective Boundary on the Transient Radiative
Transfer Process.”
Int. J. Thermal Sciences

40(6),
pp. 539
-
549, (2000).

16.

Fleck, J. A. “The Calculation of Nonlinear Radiation
Transport by a Monte Ca
rlo Method.”
Methods in
Computational Physics
. Ed. B. Alder, S. Fernbach
and M. Rotenberg. New York and London:
Academic Press 1 (1963): 43
-
65.

17.

Sakami, M., K. Mitra, and P.
-
f. Hsu. “Transient
Radiative Transfer in Anisotropically Scattering
Media using

Momotonicity
-
Preserving Schemes.”
The ASME Int. ME Cong. & Expo. Orlando, 2000.

18.

Pacheco, P. S.
Parallel Programming with MPI.

California: Morgan Kaufmann, 1997.

19.

Siegel, R., and J. R. Howell.
Thermal radiative Heat
transfer.

3
rd

ed. New York: McGra
w
-
Hill, 1992: 864.

20.

Ozisik, M.N.,
Radiative Transfer and Interaction
with Conduction and Convection
. New York: Wiley,
1973: 251.

21.

Modest, M. F.
Radiative Heat Transfer
. New York:
McGraw
-
Hill, 1993: 303.

22.

Lux, I. and K. Koblinger.
Monte Carlo Parti
cle
Transport Methods: Neutron and Photon
Calculations
. Boca Raton, Florida: CRC Press,
1991.

23.

Cashwell, E. D., and C. J. Everett.
Monte Carlo
Method for Random Walk Problems.

Great Britain:
The Lewes Press, 1959.

7

American Institute of Aeronautics and Astronautics

24.

Wang, L., S. L. Jacques, and L. Zhe
ng. “MCML

Monte Carlo Modeling of Light Transport in
Multi
-
layered tissues.”
Computer Methods and
Programs in Biomedicine
. 47 (1995): 131
-
146.

25.

Henyey, L. G., and J. L. Greenstein. “Diffuse
Radiation in the Galaxy.”
Astrophys. J.

93 (1941):
70
-
83.

26.

S
awetprawichkul, A. “ Application of the Monte
Carlo Method in the Solution of Radiative
Transfer Problems.” M.S. Thesis. U. of Arizona,
1999.

27.

Manke, J. W. “Parallel Computing in Aerospace.”
Parallel Computing
. 27.4 (2001): 329
-
336.

28.

Rifai, S. M. “Au
tomotive Design Applications of
Fluid Flow Simulation on Parallel Computing
Platforms.”
Computer Methods in Applied
Mechanics and Engineering.

184.2 (2000): 449
-
466.

29.

Thole, C. and K. Stueben. “Industrial Simulation
on Parallel Computers.”
Parallel Comp
uting
.
25.13 (1999): 2015
-
2037.

30.

Drake, J. and I. Foster. “Introducing to the Special
Issue on Parallel Computing in Climate and
Weather Modeling.”
Parallel Computing
. 21.10
(1995): 1539
-
1544.

31.

Turk, J. A. “Acceleration Techniques for the
Radiative A
nalysis of General Computational
Fluid Dynamics Solutions using Reverse Monte
-
Carlo Ray Tracing.” Ph.D. Dissertation. Virginia
Polytechnic Institute and State University, 1994.

32.

Martino, R. L., et al.
“Role of High Performance
Parallel Computing in Biol
ogical Research.”
Aunnual International Conference of the IEEE
Engineering in Medicine and Biology Society.
16.2 (1994): 1386
-
1387.

33.

Board, J. “Grand Challenges in Biomedical
Computing.”
Critical Reviews in Biomedical
Engineering.

20.1. (1992): 1
-
24.

34
.

Wu, D., et al. “Monte Carlo Simulation of Light
Propagation in Skin Tissue Phantoms using a
Parallel Computing Method.” Proceedings of
SPIE. 3914 (2000): 291
-
299.

35.

Gropp, William., et al.

MPI
-
The Complete
Reference

Vol 2. Massachusetts: MIT, 1998. 2
v
ols.

36.

Quinn, Michael.
Parallel Computing: Theory and
Practice.

2
nd

ed. New York: McGraw
-
Hill, 1994.

37.

Sawetprawichkul, A., P.
-
F. Hsu, K. Mitra., and M.
Sakami “A Monte Carlo Study of the Transient
Radiative transfer Within The One
-
Dimensional
Multi
-
La
yered Slab.” International Mechanical
Engineering Congress and Exposition, Orlando,
FL, November 2000: 4
-
9.

8

American Institute of Aeronautics and Astronautics



10
-4
10
-3
0
20
40
60
80
100
MC
DOM
t*

=10,


=0.998
a
1
=0
a
1
=0.9
a
1
=-0.9
Transmittance



Fig. 1 Comparison of one
-
dimensional Monte Carlo
and discrete ordinates solutions for linear anisotropic
scattering media.




10
-4
10
-3
10
-2
10
-1
0
10
20
30
40
50
60
70
RF 1D
RF 3D
TR 3D
TR 1D
Reflectance & Transmittance
Non-dimensional time, t*
Optical depth,

= 10
Albedo,

= 0.998
pulse width, t
p
* = 1
x:y:z = 1000:1000:1
Aspect ratio (AR) = 1000
10
6
bundles


Fig. 2 Val
idation of the reflectance and transmittance
of the 3D model with 1D model using infinite slab in
3D code.





10
-2
10
-1
10
0
10
5
10
6
10
7
RMS error
number of bundles
Optical depth,

= 10
Albedo,

= 0.998
pulse width, t
p
* = 1
Aspect ratio (AR) = 1000


Fig. 3 RMS error of the 3D Monte Carlo solutions.






10
-4
10
-3
10
-2
10
-1
0
10
20
30
40
50
60
70
1D
AR = 10
AR =100
AR = 1000
Reflectance
Non-dimensional time, t*
Optical depth,

= 10
Albedo,

= 0.998
pulse width, t
p
* = 1
10
6
bundles


Fig. 4 Parametric study of the reflectance for 3D using
different aspect ratios to

simulate 1D case.





9

American Institute of Aeronautics and Astronautics

10
-4
10
-3
10
-2
0
10
20
30
40
50
60
70
1D
AR = 10
AR = 100
AR = 1000
Transmittance
Non-dimensional time, t*
Optical depth,

= 10
Albedo,

= 0.998
pulse width, t
p
* = 1
10
6
bundles


Fig. 5 Parametric study of the transmittance for 3D
using different aspect ratios to simulate 1D case.






10
-5
10
-4
10
-3
0
10
20
30
40
50
60
point A (center)
point B (edge)
point C (corner)
Reflectance
Non-dimensional time, t*
Optical depth,

= 10
Albedo,

= 0.998
pulse width, t
p
* = 1
10
9
bundles
AR = 10



Fig. 6 Reflectance measurement at different locations.







10
-6
10
-5
0
10
20
30
40
50
60
point A (center)
point B (edge)
point C (corner)
Transmittance
Non-dimensional time, t*
Optical depth,

= 10
Albedo,

= 0.998
pulse width, t
p
* = 1
10
9
bundles, AR = 10


Fig. 7 Transmittance measurement at different locat
ions.







10
-7
10
-6
10
-5
10
-4
0
0.2
0.4
0.6
0.8
1
t* = 15
t* = 30
t* = 45
t* = 60
Reflectance
Non-dimensional distance, x*
Optical depth,

= 10
Albedo,

= 0.998
pulse width, t
p
* = 1
10
9
bundles
AR = 10



Fig. 8 The distribution of the reflectance at different
time (a single node irradiation).






10

American Institute of Aeronautics and Astronautics

10
-7
10
-6
10
-5
0
0.2
0.4
0.6
0.8
1
t* = 15
t* = 30
t* = 45
t* = 60
Transmittance
Non-dimensional distance, x*
Optical depth,

= 10
Albedo,

= 0.998
pulse width, t
p
* = 1
10
9
bundles
AR = 10


Fig. 9 The distribution of the transmittance at different
time (a single node irradiation).



































0
5
10
15
20
25
30
35
40
0
5
10
15
20
25
30
35
40
Actual speedup
Ideal speedup
speedup
number of processors
Optical depth,

= 10
Albedo,

= 0.998
pulse width, t
p
* = 1
1000 timesteps


Fig. 10 Th
e speedup results using different number of
processors.