A Hybrid Monte Carlo, Discontinuous Galerkin method for linear kinetic
transport equations
Johannes Krotza,, Cory D. Hauckb , Ryan G. McClarrenc
a Department of Mathematics, University of Tennessee Knoxville, Knoxville, 37996, TN, USA
b Computational Mathematics Group, Computer Sciences and Mathematics Division,Oak Ridge National Laboratory, Oak

Ridge, 37831, TN, USA
c Department of Aerospace and Mechanical Engineering, University of Notre Dame, Notre Dame, 46545, IN, USA

Abstract
We present a hybrid method for time-dependent particle transport problems that combines Monte Carlo
(MC) estimation with deterministic solutions based on discrete ordinates. For spatial discretizations, the MC
algorithm computes a piecewise constant solution and the discrete ordinates uses bilinear discontinuous finite
elements. From the hybridization of the problem, the resulting problem solved by Monte Carlo is scattering
free, resulting in a simple, efficient solution procedure. Between time steps, we use a projection approach
to “relabel” collided particles as uncollided particles. From a series of standard 2-D Cartesian test problems
we observe that our hybrid method has improved accuracy and reduction in computational complexity of
approximately an order of magnitude relative to standard discrete ordinates solutions.
Keywords: Hybrid stochastic-deterministic method, Monte Carlo, kinetic equations, particle transport

1. Introduction
Numerical methods for kinetic transport equations are commonly divided into two classes: deterministic
and Monte Carlo. Each of these approaches has strengths and weaknesses that complement the other.
Deterministic methods [23] directly discretize phase space (physical space, direction of flight, particle
energy) as well as time, in the time-dependent setting. For this large seven-dimensional space (three for
physical space, two for direction of flight, one for energy, and one for time), it is difficult to construct high
resolution solutions for general problems. Indeed, the number of operations and the memory footprint required
for deterministic solvers can be a challenge, even for leadership-class computers.
Monte Carlo methods [33], on the other hand, use sampling techniques to simulate particle transport
processes. In its most basic form, the Monte-Carlo procedure is a computational analog of the actual physical
processes being simulated: particles are sampled from sources and boundary conditions, then tracked as
they stream through the domain, and along the way undergo scattering or absorption interactions with the
material medium. As the particle traverses the physical domain, it contributes to integrated quantities of
interest such as particle density or net fluence through a surface. For linear problems, the central limit
theorem implies that the Monte Carlo solution is exact in the limit of an infinite number of samples [33].
Unlike deterministic methods, Monte Carlo methods are easy to extend to complicated 3-D geometries and
can handle physical processes (such as particle interactions with the background material) in a continuous
manner. Nevertheless, the uncertainty in Monte Carlo methods, as expressed in the standard deviation of an
estimate, scales like N −1/2 , where N is the number of sample particles. Additionally, Monte Carlo methods
are not well-suited for obtaining uniform spatial estimates due to the difficulty of getting sufficient samples in
every region of the physical domain. Moreover, for nonlinear problems such as thermal radiative transfer,
the Monte Carlo approach loses some of its attractive properties. For example the discretization of material
temperature, to which the particles are coupled, means that an exact solution is not obtained in the limit of
infinite samples [42, 12, 25]. Nevertheless, producing efficient, accurate Monte Carlo calculations is an active
area of research [32, 34].
Hybrid methods have been developed to harmonize the benefits of Monte Carlo and deterministic methods
while minimizing their respective drawbacks. For steady-state nuclear reactor problems, methods such as

Preprint submitted to Journal of Computational Physics

December 6, 2023

COMET [26, 43] use local Monte Carlo calculations to estimate properties of solutions in macroscopic regions
of the problem and then use a deterministic procedure to couple these regions together. Other work has
considered weight windows and other biasing techniques [36, 6, 7, 37, 27] wherein deterministic solutions are
used to modify the flight of particles in Monte Carlo calculations so that computational effort is spent more
efficiently. High-order low-order (HOLO) schemes [28, 41, 29] have been developed in which Monte Carlo is
used to solve compute a closure term for a low-order, moment-based deterministic calculation.
This work presents a deterministic-Monte Carlo hybrid method for time-dependent problems based on
the physics of particle transport. Previous work [20, 9, 10, 21, 39] has exploited the fact describing a particle
from its emission at the beginning of a time step to its first collision benefits from a different numerical
treatment than a particle that emerges from a scattering interaction with the background medium. Because
the scattering process relaxes particles towards a weakly anisotropic angular distribution, one can combine
methods that are appropriate for particle streaming for the uncollided particles during a time step with
methods that are suitable for weakly anisotropic angular distributions. In previous work, deterministic
methods with a large number of angular degrees of freedom were used for the uncollided particles while
low-resolution deterministic methods were used for the collided particles. It was also found that different
levels of energy resolution can also be employed in this type of hybrid approach [40].
Despite the benefits of deterministic hybrid methods, solutions still require a large number of degrees of
freedom for problems with large streaming paths. A natural strategy to address this challenge, which for
steady-state problems was first proposed in [3], is to use Monte Carlo for the uncollided particles. Indeed,
in many respects this is the ideal situation for a Monte Carlo approach. During a time step, particles are
tracked through the computational domain and a non-analog estimator of the solution known as implicit
capture is employed, thereby avoiding the need to consider collisions at all. Thus the calculation of the
contribution to the solution from uncollided particles is essentially a ray tracing algorithm, which has many
efficient implementations on modern computing hardware [2].
The hybrid method considered here uses Monte Carlo to compute the contribution to the solution from
uncollided particles and an efficient deterministic calculation for the collided particles. A key advancement
for extending the original steady-state formulation to time-dependent problems is a remapping step that
resamples particles from the deterministic collided solution back into the uncollided component. This
procedure is critical since, otherwise, the number of uncollided particles will decay exponentially and the
hybrid solution will relax to a low-resolution, deterministic approximation of collided solution [20]. We find
that the hybrid approach leads to more accurate solutions obtained using lower computational complexity
than pure deterministic calculations.
The remainder of the paper is organized as follows. In Section 2, we introduce the hybrid method in the
context of a single-group transport equation that is independent of particle energy. We also summarize the
numerical methods used for the uncollided and collided components of the hybrid. In Section 3, we present
numerical results for several standard test problems in a reduced two-dimensional geometry in physical space.
In Section 4, we summarize findings and present directions for future work. A short appendix describes the
Monte Carlo implementation of a boundary for one of the test problems.
2. Basics of the Hybrid Method
2.1. Transport equation
Let X ∈ R3 be a spatial domain with Lipschitz boundary and let S2 be the unit sphere in R3 . Let
Ψ = Ψ(x, Ω, t) be the angular flux depending on the position x = (x, y, z) ∈ X, the direction of flight Ω ∈ S2
and time t > 0. We assume that Ψ is governed by the linear transport equation
1
σs
∂ t Ψ + Ω · ∇ x Ψ + σt Ψ =
Ψ + Q,
c
4π

x ∈ X,

Ω ∈ S2 ,

t > 0,

(1)

where σt = σt (x), σs = σs (x) and σa = σt − σs are the total, scattering, and absorption cross-sections of the
material, respectively; Q = Q(x, Ω, t) is a known particle source; and angle brackets denote integration over
the unit sphere:
Z
Ψ = Ψ dΩ.
(2)
S2

2

The constant c > 0 is the particle speed; we assume that c = 1 for the remainder of this work. The transport
equation (1) is equipped with initial conditions
Ψ(x, Ω, 0) = Ψ0 (x, Ω),

Ω ∈ S2 ,

(3)

Ω · n(x) < 0,

(4)

x ∈ X,

and boundary condition
x ∈ ∂X,

Ψ(x, Ω, t) = G(x, Ω, t)

where Ψ0 and G are known and n(x) is the unit outward normal at x ∈ ∂X.
2.2. The hybrid method
The hybrid method is based on a first collision source splitting [3]. Let Ψ = Ψu + Ψc , where the uncollided
flux Ψu and the collided flux Ψc satisfy the following system of equations
∂t Ψu + Ω · ∇x Ψu + σt Ψu = Q,
σs
(Ψu  + Ψc ) .
∂ t Ψ c + Ω · ∇ x Ψ c + σt Ψ c =
4π

(5a)
(5b)

Due to the linearity of (1), the splitting in (5) is exact; that is, if Ψu and Ψc solve (5a) and (5b), respectively,
then Ψu + Ψc solves (1). In practice, however, (5a) and (5b) are solved at different resolutions or even with
different methods. Typically (5a) is solved with a method that has high resolution in angle, and because
(5a) has no coupling in angle, it is easier to solve than (1) and also easy to solve in parallel. Meanwhile (5b)
inherits the angular coupling in (1), but typically requires less angular resolution.
Since (5a) has no scattering source, particle densities will be transferred into the collided flux at an
exponential rate, thus making the accuracy at which (5b) is solved the driving factor in the overall accuracy.
This effect can be mitigated by abusing the autonomous nature of the equations and periodically relabeling
the collided flux as uncollided at every time step.
To describe the implementation of the hybrid in more detail, let T be an operator such that
u[t] = T (t, t′ , s, υ, b, λt , λs )
where u[t](x, Ω) := u(x, Ω, t), satisfies

λs


 ∂t u + Ω · ∇x u + λt u = 4π u + s,
′
 u(x, Ω, t ) = υ(x, Ω)


u(x, Ω, t) = b(x, Ω, t)

(6)

x ∈ X,

Ω ∈ S2 ,

x ∈ X,

Ω ∈ S2 ,

x ∈ ∂X,

t > t′ ,

Ω · n(x) < 0,

(7a)
(7b)
′

t>t.

(7c)

with source s = s(x, Ω, t). Using the operator T , we can write
Ψ[tn+1 ] = T (tn+1 , tn , Q, Ψ[tn ], G, σt , σs ),
Ψu [tn+1 ] = T (tn+1 , tn , Qu , Ψ[tn ], G, σt , 0),
Ψc [tn+1 ] = T (tn+1 , tn , Qc , 0, 0, σt , σs ),

(8)
Qu = Q
σs
Qc =
Ψu .
4π

(9)
(10)

We simulate the system (5) using a Monte Carlo method for the uncollided equation (5a) and a deterministic
discretization of the collided equation (5b). Let
TMC (t, t′ , s, υMC , b, λt , λs ; Np )

(11)

be the Monte Carlo approximation to (7) given a particle representation υMC of υ and using Np pseudoparticles to represent the distribution of particles in phase space introduced by the source s over the internal
(t′ , t). Similarly,
TSN (t, t′ , s, v, b, λt , λs ; N, Nx )
(12)
denote the SN -DG approximation of (7) using a level N set of ordinates, Nx spatial cells per dimension
with Q1 elements, and a backward Euler time discretization to get from t′ to t. (The Monte Carlo method
3

and SN-DG method are presented in more detail below.) Then given Np , N , and Nx , and a Monte Carlo
approximation ψ n of Ψ(tn ), let
ψun+1 = TMC (tn+1 , tn , Qu , ψ n , G, σt , 0; Np ),

Qu = Q
σs
Qc =
Ψu MC
4π
σs
Q r = Qc +
Ψc SN
4π

ψcn+1 = TSN (tn+1 , tn , Qc , 0, 0, σt , σs ; N, Nx ),
Rψcn+1 = TMC (tn+1 , tn , Qr , 0, 0, σt , 0; Np ),
ψ

n+1

= ψun+1 + Rψcn+1

(13)
(14)
(15)
(16)

where R is the remapping operator and ·MC and ·SN denote approximation of the angular integral over S2
with the respective method.
2.3. Discrete ordinate-discontinuous Galerkin
The discrete ordinates (SN ) method [5] approximates (7) by replacing the angular integral u by a
discrete quadrature and then solving the resulting equation for the angles in the quadrature. This procedure
yields a system of equations that depend only on space and time and can be further discretized by a variety
of methods. Let
NΩ
Ω
{Ωq }N
(17)
q=1 and {ωq }q=1
be the angles and associated weights of the SN quadrature, where NΩ = NΩ (N ) depends on the specific type
of quadrature set being used. After discretizing in angle and applying an implicit Euler time discretization,
the following semi-discrete system is obtained for each q ∈ {1, ..., NΩ } and n ∈ {0, 1, 2, . . . },

NΩ

X

 1 un+1 − un  + Ω · ∇ un+1 + λ un+1 = λs
ωr un+1
+ sn+1
,
x ∈ X,
(18a)
q
x
t
r
q
q
q
q
∆t q
4π r=1


 un+1 (x) = bn+1 (x),
x ∈ ∂Xq− ,
(18b)
q
q

where ∂Xq− = {x ∈ X : Ωq · n(x) < 0}, bnq (x) = b(x, Ωq , tn ), snq (x) = s(x, Ωq , tn ), and unq (x) ≈ u(x, Ωq , tn )
is the approximation on the temporal and angular grid. After reassigning
uq ← un+1
,
q

sq ← sn+1
+ unq ,
q

λt ← λt +

1
,
∆t

and

bq ← bn+1
,
q

(19)

the discretization in (18a) can be written in the equivalent steady-state form


Ωq · ∇x uq + λt uq = λs ū + sq ,
uq (x) = bq (x),

1
where u = (u1 , . . . , uNΩ )⊤ , ū := 4π

P

x∈X

x ∈ ∂Xq−

(20a)
(20b)

ωr ur .

r

We discretize (20a) in physical space with a discontinuous Galerkin method and upwind numerical
fluxes. The method by now is fairly standard (see for example [9, 19]) and is often used because of its
accuracy in scattering-dominated regimes relative to upwind finite-difference and finite-volume methods
[1, 22, 31, 18]. Because the DG method is well-known, we summarize it only briefly for the case of a
two-dimensional Cartesian mesh with rectangular cells, which is sufficient for all of the numerical tests in
Section 3. Let X be divided into open sets Ci,j that are squares with side lengths h centered at (xi , yj ),
and let Vh = {v ∈ L2 (X) : v|Cij ∈ Q1 }. The goal will be to find an approximation of the weak solution of
equation (20a); that is, find uh = (uh1 , . . . , uhNΩ )⊤ ∈ [Vh ]NΩ := Vh × · · · × Vh such that
| {z }
NΩ times

A(i,j)
(uhq , vqh ) + Pq(i,j) (uhq , vqh ) = M(i,j)
(uhq , vqh ) + R(i,j) (uh , vqh ) + Sq(i,j) (vqh ) + Bq(i,j) (vqh )
q
q

4

(21)

for all i, j ∈ {1, . . . , Nx }, q ∈ {1, ..., N }, and vh ∈ VhNΩ . Here
Z
Z
(i,j) h h
h h
vqh uhq dx,
Aq (uq , vq ) = −
(Ωq · ∇x vq )uq dx + λt
Ci,j

Ci,j

Z

Pq(i,j) (uhq , vqh ) =

(Ωq · n)vqh,− vqh,− ds(x),

(Ωq · n)vqh,− uh,+
q ds(x),

(23)

(∂Ci,j )−
q

Z

ūh vqh dx,

Sq(i,j) (vqh ) =

Ci,j

vqh,± (x) =

Z

M(i,j)
(uhq , vqh ) =
q

(∂Ci,j )+
q

R(i,j) (uh , vqh ) = λs

(22)

Z

sq vqh dx,

Ci,j

lim vqh (x ± ϑn),
ϑ→0+

and

Z

Bq(i,j) (vqh ) =

bq vqh ds(x),

(24)

Ci,j ∩∂Xq−

(∂Ci,j )±
q = {x ∈ ∂Ci,j : ±Ωq · n(x) > 0}

(25)

We construct uh = (uh1 , . . . , uhN )⊤ as follows: For each i, j ∈ {1, . . . , Nx } and q ∈ {1, ..., NΩ }, let
X (i,j) (i,j)
uhq (x) =
αq,k φk (x, y),
x ∈ Ci,j ,

(26)

|k|∞ ≤1

where k = (k1 , k2 ),
(i,j)
φk = Pk1



x − xi
h/2

Pk 2



y − yj
h/2

,

(27)

and Pk is the usual Legendre polynomial of degree k on ξ ∈ [−1, 1] with normalization

R1

P (ξ)Pk2 (ξ)dξ =
−1 k1
(i,j)
2
h
δ
.
Using
this
representation
for
u
,
we
derive
the
following
linear
system
for
the
coefficients αq,k
2k1 +1 k1 ,k2
from equation (21): For each l such that |l|∞ ≤ 1,
X

(i,j)

(i,j)

(i,j)

Aq,l,k + Pq,l,k αq,k =

|k|∞ ≤1
(i,j)

where ᾱk

(i,j)

=

N
Ω
P

q=1

X

(i,j)

(i∗ ,j ∗ q)

Mq,l,k αq,k

(i,j)

(i,j)

+ Rl,k ᾱk

(i,j)

(i,j)

+ Sq,l + Bq,l ,

(28)

|k|∞ ≤1
(i,j)

wq αq,k ,
(i,j)

Aq,l,k = A(i,j)
(φk
q

(i,j)

, φl

),

(i,j)

(i,j)

Pq,l,k = Pq(i,j) (φk

(i,j)
(i,j)
Bq,l = Bq(i,j) (φl ),

(i,j)

, φl

(i,j)
(i,j)
Sq,l = Sq(i,j) (φl ),

(i,j)

(i∗ ,j ∗ )

Mq,l,k = M(i,j)
(φk
q

)

(i,j)

, φl

(i,j)
(i,j)
(i,j)
Rl,k = R(i,j) (φk , φl )

)

(29)
(30)

and, given the components n = (nx , ny ) of the outward normal, the indices i∗ , j ∗ are given by
i∗ (x) = i + nx (x)

and

j ∗ (x) = j + ny (x).

(31)

To improve readability, we rewrite equation (28) as a matrix equation with respect to the indices k and l:
∗

∗

,j )
α(i,j)
= R(i,j) ᾱ(i,j) + M(i,j)
α(i
+ B(i,j)
+ S(i,j)
.
A(i,j)
+ P(i,j)
q
q
q
q
q
q
q

(32)

The organization of the operators in (32) reflects a standard solution strategy combining source iteration and
sweeping. In source iteration, ᾱ(i,j) is lagged; that is, given an iteration index ℓ:
+ P(i,j)
αq(i,j,ℓ+1) = A(i,j)
q
q
ᾱ(i,j,ℓ+1) =

N
X

−1

∗

∗

,
R(i,j) ᾱ(i,j,ℓ) + M(i,j)
αq(i ,j ,ℓ+1) + B(i,j)
+ S(i,j)
q
q
q

(33a)

wq αq(i,j,ℓ+1) .

(33b)

q=1

(i∗ ,j ∗ ,ℓ+1)

Sweeping refers to process solving of (33) cell-by-cell: for each q, cells can be ordered such that αq
(i,j,ℓ+1)
known, prior to solving for αq
. The details of this procedure are given in Algorithm 1.
5

is

2.4. Monte Carlo
In this section, we describe the Monte Carlo method used to compute the solution to (7) for the pure
absorption problem when λt = λa (i.e. no scattering):


 ∂t u + Ω · ∇x u + λt u = s,
u(x, Ω, 0) = υ(x, Ω)


u(x, Ω, t) = b(x, Ω, t)

Ω ∈ S2 ,

x ∈ X,

t > 0,

(34a)

x ∈ X, Ω ∈ S2 ,
x ∈ ∂X, Ω · n(x) < 0,

(34b)
(34c)

t > 0.

The Monte Carlo method approximates the phase space distribution u using a finite set of pseudo-particles:
X
wπ (t)δ(x − xπ (t))δ(Ω − Ωπ ),
(35)
u(x, Ω, t) ≈
π∈Πt

where Πt is a set of pseudo-particles π with position xπ (t), weight wπ (t), and direction of flight Ωπ . Πt will
be defined more carefully below. A benefit of the hybrid is that scattering processes, which can slow down
the method significantly [13], do not need to be modeled in (34).
The Monte Carlo implementation of (34) can be derived via a Green’s function formulation for (34). Let
G(x, y, Ω, t, t0 ) solve
∂t G + Ω · ∇x G + λt G = δ(x − y)δ(t − t0 ),

(36)

with zero initial data and boundary conditions. Then

u(x, Ω, t) =
+
+

Z tZ
0

R3

Z tZ
0

R3

Z tZ

R3

0

G(x, y, Ω, t, t0 )s(y, Ω, t0 )dydt0

(37a)

G(x, y, Ω, t, t0 )sv (y, Ω, t0 )dydt0

(37b)

G(x, y, Ω, t, t0 )sb (y, Ω, t0 )dydt0

(37c)

solves (34), where s, sv , and sb are identically zero outside of the closure of X. The terms sv and sb are
provisional source terms designed such that (37b) solves (34) when s = 0 and b = 0, while (37c) solves (34)
when s = 0 and v = 0. The former is solved by setting sv (x, Ω, t) = v(x, Ω)δ(t). However, determining sb can
be slightly more involved, and an example that is used for numerical experiments in Section 3.3 is provided
in the Appendix. Once s, sv , and sb are known, all three terms can be treated identically. For simplicity, we
restrict our attention to (37a) below.
For any t > t0 and any fixed Ω, let X(t) = x0 + (t − t0 )Ω. Then g(t) = G(X(t), y, Ω, t, t0 ) solves
dg(t)
= −λ(X(t))g(t) + δ(X(t) − y)δ(t − t0 )
dt

(38)

or, equivalently,
−

Rt

Rt

λt (x−(t−ξ)Ω)dξ

g(t) = e

t0

λt (X(ξ))dξ

δ(x0 − y).

(39)

Setting x = X(t) in (39) gives
G(x, y, Ω, t, t0 ) = e

−

t0

δ(x − (t − t0 )Ω − y).

Plugging (40) back into (37a) yields, after some manipulation,

6

(40)

u(x, Ω, t) =

Z tZ
0

e

−

Rt

t0

λt (x−(t−ξ)Ω)dξ

R3

−

Rt

Z t

e

−

Rt

Z t

e− 0 λt (x−ξΩ)dξ s(x − τ Ω, Ω, t − τ )dτ.

=

t0

λt (x−(t−ξ)Ω)dξ

0

t−τ

0

s(x − (t − t0 )Ω, Ω, t0 )dt0

λt (x−(t−ξ)Ω)dξ

0

=

δ(x − (t − t0 )Ω − y)s(y, Ω, t0 )dydt0

e

=

Z t

(41)

s(x − τ Ω, Ω, t − τ )dτ

Rτ

This representation of u(x, Ω, t) can be interpreted as the density of particles that have reached the location
x at time t while moving in the direction Ω. These particles are emitted by the source s at time t − τ and
location x − tΩ, and they carry a weight that decays exponentially due to absorption.
The Monte Carlo approach can be understood as an approximation of u based on sampling of pseudoparticles from the source s in (41). Let
X
wπ δ(x − xπ )δ(Ω − Ωπ )δ(t − tπ )
(42)
s̃(x, Ω, t) =
π∈Π

where π ∈ Π are pseudo-particles with weight wπ > 0, position xπ ∈ X, and direction of flight Ωπ ∈ S2 at
tπ > 0, such that s̃(x, Ω, t) approximates s(x, Ω, t) for all t > 0. Then
for any C ⊂ X, any B ⊂ S2 , and any time interval (tn , tn+1 ), the representation of u in (41), along with
the approximation s̃ gives
Z tn+1 Z Z
u(x, Ω, t)dΩdxdt
tn

≈
=

C

e− 0 λt (x−ξΩ)dξ s̃(x − τ Ω, Ω, t − τ )dτ dΩdxdt

Z tn+1 Z Z Z t

e− 0 λt (x−ξΩ)dξ

C

tn

C

tn

=

B

B

0

Rτ

Rτ

0

X

π∈Π

wπ δ(x − τ Ω − xπ )δ(Ω − Ωπ )δ(t − τ − tπ )dτ dΩdxdt

X Z tn+1

wπ e− 0

X Z tn+1

wπ (t)1C (xπ + (t − tπ )Ωπ )1[0,t] (tπ )1B (Ωπ )dt

π∈Π

=

B

Z tn+1 Z Z Z t

π∈Π

R t−tπ

λt (xπ +ξΩπ )dξ

tn

tn

(43)

1C (xπ + (t − tπ )Ωπ )1[0,t] (tπ )1B (Ωπ )dt

R t−tπ

λt (xπ +ξΩπ )dξ
and wπ (tπ ) = wπ .
where wπ (t) = wπ e− 0
Note that with the identification xπ (t) = xπ + (t − tπ )Ωπ we can also identify Πt in (35) as Πt = {π ∈
Π : tπ < t}.
We will denote (43) as TMC (tn+1 , tn , 0, λs , λs , s, 0; Np ) as the Monte Carlo solution for the case of the zero
initial data and zero boundary data. A general Monte-Carlo solution can be obtained as

TMC (tn+1 , tn , s, υ, b, λs , λs ; Np ) =TMC (tn+1 , tn , s + sb , v, 0, λs , λs ; Np )
=TMC (tn+1 , tn , s + sv + sb , 0, 0, λs , λs ; Np )
Numerically it is useful to realize that
Πn+1 := Πtn+1 = π(tn + ∆t) : π ∈ Πtn ∪ {π ∈ Π : tπ ∈ (tn , tn+1 )}

(44)

where in an abuse of notation π(tn + ∆t) is a new particle with position xπ (t + ∆t), weight wπ (t + ∆t) and
direction of flight Ωπ for some π ∈ Πtn = Πn . Thus particles Πn+1 can be obtained by updating the weight
7

and position of particles in Πn and sampling new particles with birth times in (tn , tn+1 ). From (43) we also
obtain
X Z tn+1
uMC (x, t) =
wπ (t)1C (xπ + (t − tπ )Ωπ )1[0,t] (tπ )dt
(45)
π∈Π

tn

With this formulation, particle weights wπ (t) decrease exponentially at a rate proportional to the absorption,
given by λt . This approach is known as the implicit capture method. It avoids the need to sample absorption
times explicitly and, at the same time, reduces statistical noise and simplifies the implementation [14, Page
168][24, Chapter 22].
The Monte Carlo simulation of (34) from tn to tn+1 is based on (43) and proceeds according to the
following steps:
1. Let P be a partition of X into disjoint cells C. For each C ∈ P, calculate the total weight W C of new
particles generated in C by the source during the interval (tn , tn+1 ):
W

C

1 1
=
∆t 4π

Z Z tn+1 Z
S2

tn

s(x, Ω, t) dΩdxdt,

(46)

C

P
where ∆t = tn+1 − tn , and let W = W C C∈P . Let Np be the input for the total number of new
particles to sample during the interval (tn , tn+1 ). Then for each C ∈ P, sample
 C
W
NpC = floor
(47)
W
particles and assign them each a weight w = W/NpC . The total number of particles in the system at
P C
Np . Each new particle p is assigned a position xπ sampled uniformly
this time is Ñpn+1 = Ñpn+1 +
C∈P

from C and a birth time tn+1 − τπ , where τπ is sampled uniformly from (0, ∆t). Each particle p is also
assigned an angle Ωπ . For isotropic sources, Ωπ is sampled uniformly from S2 . For non-isotropic sources,
(such as the boundary source for the holhraum problem in Section 3.3), the sampling distribution
must be consistent with the angular dependence. The particles, including their space, angle, and time
coordinates, are added to current particle list.
2. Move each particle π in the current particle list from xπ to xπ + τπ Ωπ and update its weight to
wπ ← wπ (tn+1 ). The number of particles in the system will have to be adjusted accordingly. Reset its
remaining time to τπ ← ∆t. For each cell C ∈ P, update the sum in (43) and (45) according the time
spent in C during the interval (tn , tn+1 ).
3. To reduce computational effort and memory, at the end of each time step any particle p with weight
wπ < wkill will be dropped with a probability of pkill > 0. Here wkill > 0 is a user-defined parameter,
called the ‘killing weight’, and pkill = (1 − wπ /wkill ). To preserve the total mass in the system, any
particle p with weight w < wkill that survives this ‘Russian roulette’ [24, Chapter 22] will have its
weight readjusted to wπ /(1 − pkill ).
2.5. Pseudocode
The algorithms that we use for our numerical results are detailed in Algorithms 1-3. The SN method is
given in Algorithm 1. The hybrid method is given in Algorithm 3. It requires Algorithm 1 for the collided
component and Algorithm 2 for the Monte Carlo update. A listing of the notation used these algorithms is
provided in Table 1.
3. Numerical Results
In this section, we compare simulation results from the Monte Carlo-SN hybrid method to those from
a standard, monolithic SN method. The goal is to demonstrate that the hybrid method provides a more
efficient approach. The SN computations for the monolithic method and for the collided component of the
hybrid rely on product quadrature sets on the sphere [38, 4].

8

User defined parameters
number of Cartesian cells along each dimension
order of discrete ordinates
number of new particles (up to rounding) generated
from source
time step
tolerance of iteration
Gauss-Legendre weights
Material parameters
total, absorption and scattering crosssection
Additional notation

Nx
N
Np
∆t
δ
ωq
λt , λa , λs
(i,j)

(i,j)

Aq , Pq
U (X)

(i,j)

, R(i,j) , Mq

(i,j)

, Bq

(i,j)

, Sq

matrices defined in Section 2.3
uniform distribution on set X.

Table 1: Pseudocode parameters

Algorithm 1 SN -algorithm: Propagate solution from tk to tk+1
Input: un , sq
⊲ coefficients of solution ukq from previous step and source
Require: λ̂t , λa , λs ,
⊲ Material properties
NΩ
x
Require: ∆t, {Ci,j }N
,
{Ω
,
w
}
⊲
Discretization
parameters
q
q
q=1
i,j=1
Require: δ
⊲ Convergence tolerance
P
(i,j) (i,j)
(i,j)
αq,k φk (x) for x ∈ Ci,j
1: αq
(tn ) such that unq (x) =
|k|∞ ≤1

1
2: sq ← sq + ∆t
αq (tn )

⊲ Initialize source
⊲ Initialize coefficients

3: αq ← αq (tn )
4: err = δ + 1
5: while err > δ do
6:
β q ← αq

7:
8:
9:
10:
11:
12:
13:

ᾱ ←

N
Ω
P

⊲ Store old coefficients

w q αq

q=1

for q ∈ {1, ..., NΩ } do
2
for (i, j) ∈ {1, ..., Nx } do
(i,j)
αq ←

(i,j)
(i,j)
A q + Pq

⊲ Sweep through cells in direction Ωq
−1

R

(i,j) (i,j)

β̄

(i,j) (i∗ ,j ∗ )
(i,j)
(i,j)
+ Mq αq
+ Bq + S q

⊲ Update coefficients

end for
end for
err = max ||αq − βq ||∞

⊲ Discrepancy between iterations

q

14: end while

15: return un+1
(x) =
q

P

|k|∞ ≤1

(i,j) (i,j)

αq,k φk

(x), un+1
SN (x) =
q

9

P

|k|∞ ≤1

(i,j) (i,j)
φk (x) for x ∈ Ci,j

ᾱk

Algorithm 2 MC-algorithm: Propagate solution form tk to tk+1
P
Input: Πn+1 with un (x, Ω) =
wπ δ(x − xπ )δ(Ω − Ωπ )
π∈Πt

s(x, Ω, t)
Require: λt
x
Require: ∆t, {Ci,j }N
i,j , Np
Require: wkill
1: for (i, j) ∈ {1, . . . , Nx }2 do
tn+1
R R R
2:
W Ci,j = |Ci,j1 |∆t
s(x, Ω, t)dΩdxdt

⊲ previous particle distribution and source
⊲ Discretization parameters
⊲ killing weight

tn Ci,j S2

3: end for
4: W ←

P

W Ci,j

⊲ Calculate total source

i,j

5: for (i, j) ∈ {1, . . . , Nx }2 do

j Ci,j k
C
W
Np
Np i,j ←
W
7: end for
P Ci,j
Np
8: Np ←

⊲ Number of particles on each cell

6:

i,j
W
9: w ← N
p

⊲ number of new particles

Ci,j

10: for (i, j) with Np
11:
12:
13:
14:

> 0 do
C
for k ∈ {1, ..., Np i,j } do
Generate new particle π with
xπ ∼ U (Ci,j )
tn+1
R R
1
s(x, Ω, t)dxdt
(Ωx , Ωy , Ωz ) ∼ W Ci,j |C
|∆t
i,j

⊲ sample new particles from source

⊲ Draw particle’s position
⊲ Draw particle’s direction of flight

tn Ci,j

Ωπ ← (Ωx , Ωy )
wπ ← w
⊲ Assign particle weight
Π∗ ← Π∗ ∪ {π}
⊲ Add new particle to existing
end for
19: end for
20: for π ∈ Πt ∪ Π∗ do
⊲ Move particles
21:
if π ∈ Πt then
22:
τπ ← ∆t
⊲ remaining time for particles from prev. step
23:
else
24:
randomly draw τπ ∼ U ([0, ∆t])
⊲ remaining time for particles generated this time step
25:
end if
26:
xπ ← xπ + τπ Ωπ
⊲ Update particle’s position
27:
for (i, j) with Ci,j ∩ {xπ + tΩπ : t ∈ [0, τπ ]} = ∅ do ⊲ All cells intersected by particle’s trajectory
Rt
Rτ
28:
Φi,j ← Φi,j + wπ 0 π exp − 0 λa (xp + t′ Ωp )dt′ 1Ci,j (xπ + tΩπ )dt
⊲ Update Φ
29:
end for

 Rτ
30:
wπ ← wπ exp − 0 π λa (xπ + tΩp )dt
⊲ Update particle weight
31:
if xπ ∈ X then
32:
Πn+1 ← Πn+1 ∪ {π}
⊲ Remove particles that left domain
33:
end if
34: end for
35: for π ∈ Πn+1 with wπ < wkill do
⊲ Russian roulette
36:
r ∼ U([0, 1])
π
then
⊲ Determine survival of particle
37:
if r > wwkill
n+1
38:
Π
← Πn+1 \ {π}
39:
else
40:
wπ ← wkill
⊲ Update surviving particle’s weight to approx. preserve total mass
41:
end if
42: end for
P
43: return Πn+1 with un+1 (x, Ω) =
wπ δ(x − xπ )δ(Ω − Ωπ ) and un+1 MC (x) = Φ(x)
π∈Πn+1
10
15:
16:
17:
18:

Algorithm 3 HMC-SN -algorithm: Propagate solution from tk to tk+1
P
wp δ(x − xp )δ(Ω − Ωp )
Input: Πn with un (x, Ω) =
p

s(x, Ω, t)
Require: λ̂t , λt , λa , λs ,
NΩ
x
Require: ∆t, {Ci,j }N
i,j=1 , {Ωq , wq }q=1 , Np
Require: δ, wkill
1: Πn+1
, un+1
MC ← MC(Πn , s)
u
u
2: for q ∈ {1, ..., NΩ } do
3:
s q ← λs Φu
4: end for
5: un+1
, un+1
SN ← Sn (0, sq )
c
c
6: s ← λs (un+1
SN + un+1
MC )
c
u
n+1
7: ΠR , uR MC ← MC(0, s)
8: Πn+1 ← Πn+1
∪ ΠR
u
9: un+1 MC ← un+1
MC + un+1
u
R MC
n+1
n+1
10: return Π
, u
MC

⊲ previous particle distribution and source
⊲ Material properties
⊲ Discretization parameters
⊲ Convergence tolerance and killing weight
⊲ Monte Carlo for uncollided
⊲ turn un+1
MC into source for Sn
u
⊲ Sn for collided
⊲ sources for Relabeling
⊲ Monte Carlo as relabeling

We consider three well known test problems: the line source problem [16], the lattice problem [8], and
the linearized hohlraum problem [9, 8]. The specifications for each problem are provided in the following
subsections. They are all formulated in a geometry for which ∂z Ψ = 0. This means that they can be reduced
to two dimensions in physical space and, by an abuse of notation, we write Ψ(x, Ω, t) = Ψ(x, y, Ω, t). A
further consequence of the geometry is that product quadrature on the sphere can be reduced to just the
upper hemisphere, in which case NΩ = N 2 . The time step for each problem is tied to the grid resolution via
the ratio CFL = ∆t/∆x. For all calculations shown below, the iteration tolerance δ (see Algorithm 1) is set
to 10−4 .1
For each problem, we assess the accuracy of the numerical solution and the efficiency with which it is
obtained. To quantify the accuracy we compare our results to a reference solution. For the line source problem
the reference is the semi-analytic solution from [17]; see also [16]. For the other two test problems, the reference
is a high-resolution hybrid solution based solely on SN discretization with a triangular-based quadrature
referred to as TN [35] for both collided and uncollided component, combined with a DG discretization in
space and integral deferred correction in time [11].
Accuracy for the MC-SN hybrid and the monolithic SN method is measured in terms of the relative
L2 -difference in the scalar flux Φ = Ψ at a given final time tfinal . Given the numerical solution Φnum and
the reference Φref :
||Φ − Φref ||L2
∆=
,
(48)
||Φref ||L2
P
where L2 -norm is approximated by h2 Ci,j Φ2ij and Φi,j is the average on the cell Ci,j . Because our
implementation of the hybrid method and the SN method are not run-time optimized, we use a complexity
measure which counts the number of times a particle is moved or a DG unknown is updated in the course
of a sweep. Let N c be the level of the quadruature for the collided component of the hybrid. Then the
complexity of the monolithic method is

CS N

=

4
↑

×

# of Legendre
coefficients

NΩ
↑

×

# of ordinates

Nx2
↑

# of
cells

×

Ni
↑

# of source
iterations

×

T
∆t
↑

(49a)

# of time
steps

1 While not shown here, we also tested several runs using δ = 10−8 , which leads to negligible improvements in accuracy when

compared to the results shown below.

11

while the complexity of the hybrid is
Chybrid

=

(Nu
↑

avg. # of
particles
moved for
uncollided

+

NR )
↑

avg. #
particles
moved in
relabelling

×

T
∆t
↑

# of time
steps

+

CSN c
↑
SN

(49b)

complexity of
collided
equation

tot = (Nu + NR ) T so that C
tot + CS
For convenience, we set NMC
= NMC
. Nu is the sum of particles
Nc
∆t
hybrid
added to the system Np and the average number of particles still in flight from the previous time step Nprev ,
i.e. Nu = Np + Nprev , while NR is the number of particles added in the relabeling, which we set to be
N R = Np .

3.1. The Line Source problem
In the line source problem, an initial pulse of uniformly distributed particles is emitted from the line
ℓ = {(x, y, z) : x = y = 0} into the surrounding domain X = R3 , which contains a purely scattering material
with σs = σt = 1. Because the geometry of the domain and initial condition are invariant in z, the spatial
domain can be reduced to R2 . In this two-dimensional setting, the initial condition can be represented by an
1
isotropic delta function 4π
δ(x, y), but to reduce numerical artifacts, we use a mollified version of the initial
condition:
 2
x + y2
1 1
,
ς = 0.03.
(50)
exp −
Ψ0 (x, y, Ω, t) =
4π 2πς
2ς
Meanwhile, the computational domain is restricted to the square [−1.5, 1.5]2 and equipped with zero inflow
boundary conditions.
We perform monolithic SN and MC-SN hybrid simulations at various spatial and angular resolutions.
The spatial domain is subdivided into equal Nx × Nx square cells with Nx ∈ {51, 101, 201}. For the SN -runs
we let N ∈ {4, 8, 16, 32}, resulting in NΩ = N 2 ordinates on the northern hemisphere of S2 . The collided
part of the hybrid algorithm employs an SN method with N ∈ {4, 8}. In the hybrid method the number
of particles was also changed between runs. The number of particles newly inserted into the system every
time step is roughly 2 × Np , where Np = 10k and k ∈ {2, 3, 4, 5, 6}.2 Due to rounding and particles being
dropped via Russian roulette, the exact number of particles inserted into the system varies slightly. The
killing weight is fixed at wkill = 10−15 . The CFL is also fixed at 0.5 across all runs. The reference solution is
the semi-analytic solution from [17]; see also [16].
Figure 1 depicts several approximations of the scalar flux Φ at tfinal = 1, calculated using the SN method
and the hybrid method. The solutions calculated using the SN method clearly show ray-effects that only get
resolved after a significant increase in the angular resolution. No such effects are seen in the hybrid solutions.
The hybrid solutions do contain some noise, as the particle count is relatively low, but they preserve the
symmetry of the problem up to a reasonable error. Unlike the SN -method, the hybrid is able to capture the
wave front of unscattered particles travelling away from the center with speed 1.
Figure 2 shows L2 -errors of the numerical solutions in log scale; the same trends are apparent. While the
hybrid solutions mainly suffer from noise due to the stochastic nature of the Monte Carlo method, the SN
method has strong ray-effects and struggles to capture the analytic solution at the wave front.
A more systematic analysis of the numerical results is presented in Figure 3. This plot shows the relative
L2 -error ∆ of various runs versus their respective computational complexity C. For the hybrid method,
increasing the angular resolution N in the collided component yields a marginal improvement at best. However,
changes in the particle number Np for the uncollided component have a significant impact. For the SN
method, on the other hand, increasing the angular resolution yields a significant improvement in the accuracy.
For smaller values of N , increasing the spatial resolution may actually increase the error. This is especially
apparent in Figure 3 for the S4 and S8 results. For a fixed angular resolution, additional spatial accuracy will
2 The factor of 2 is because N particles are used for the uncollided equations (13) and N particles are used for the relabeling
p
p
(15).

12

0.5
0.45
0.4
0.35
(a) Reference

(b) S16 : Nx = 101, ∆ =
27%, C = 9.2 × 109

(c) S32 : Nx = 101, ∆ =
19%, C = 3.7 × 1010

(d) HMC-S4 : Nx = 51,
N tot = 5.9 × 107 , ∆ =

(e) HMC-S8 : Nx = 101,
N tot = 2.3 × 109 , ∆ = 8%;

(f) H MC-S4 : Nx = 201,
N tot = 9.0 × 109 , ∆ = 5%,

19%, C = 8.0 × 107

C = 2.5 × 109

C = 9.9 × 109

0.3
0.25
0.2
0.15
0.1
0.05
0

MC

MC

MC

Figure 1: Numerical approximation of the scalar flux Φ for the line source problem at tfinal = 1 with CFL 0.5. Each numerical
solution is characterized by a relative L2 difference ∆ with respect to the reference, defined in (48), and a complexity C, defined
in (49a) for the monolithic SN method and (49b) for the hybrid.

begin to resolve the ray effect anomalies in the solution. Conversely, SN results with lower spatial resolution
benefit from error cancellation due to the numerical diffusion smoothing ray effects. The hybrid may also
have larger errors if the particles per cell is too low.
Overall the hybrid method outperforms the monolithic SN method. For example, the hybrid error can
match the most refined SN calculation (N = 32) with a complexity that is roughly 2-3 orders of magnitude
smaller; compare Figures 2 (c) and (d). In fact the hybrid method can obtain an error with half the size with
a complexity that is an order of magnitude less. In general hybrid runs tend to be 3-4 times more accurate
than their SN counterparts of similar complexity.
3.2. The Lattice problem
In the lattice problem, a checkerboard of highly absorbing material is embedded in a scattering material
with a central source. The layout of this problem along with its material parameters can be found in Figure
4. The computational domain is a 7 × 7 rectangle with zero inflow data at the boundaries. The center square
(red) contains an isotropic particle source, while the blue squares are pure absorbers. The red and white
squares are purely scattering with σs = σt = 1. The initial condition is identically zero everywhere in the
domain.
We perform SN and hybrid runs with varying spatial and angular resolution. The spatial domain is
subdivided into equal Nx ×Nx square cells with Nx ∈ {56, 112, 224}. For the SN -runs we use N ∈ {4, 8, 16, 32},
resulting in NΩ = N 2 ordinates on the northern hemisphere of S2 . The collided part of the hybrid algorithm

13

10 0

10 -1
(a) Reference

(b) S16 : Nx = 101, ∆ =
27%, C = 9.2 × 109

(c) S32 : Nx = 101, ∆ =
19%; C = 3.7 × 1010

(d) HMC-S4 : Nx = 51,
N tot = 5.9 × 107 , ∆ =

(e) HMC-S8 : Nx = 101,
N tot = 2.3 × 109 , ∆ = 8%,

(f) HMC-S4 :
Nx = 201, N tot = 9.0 ×

19%, C = 8.0 × 107

C = 2.5 × 109

109
∆ = 5%, C = 9.9 × 109

10 -2

10 -3

10 -4

MC

MC

MC

Figure 2: Absolute difference between the analytical solution and various numerical solutions to the line source problem at t = 1
with CFL 0.5. Each numerical solution is characterized by a relative L2 difference ∆ with respect to the reference, defined in
(48), and a complexity C, defined in (49a) for the monolithic SN method and (49b) for the hybrid.

14

•

 :HMC-S ,Nx = 51 :HMC-S ,Nx = 101 :HMC-S ,Nx = 201
n : SN , Nx = 51
n : SN , Nx = 101 △
n : SN ; Nx = 201



1.2
1

∆

0.8

N

N

1.2

4
(N = 4)
4

1
0.8

4
8

0.6

4
(N = 8)
4
107

4
8

0.6
8

8
8

0.4

0
107

108

8

0.4
16 32
16

0.2

109

tot
NM
C

N

1010

16
32

16 32
16

0.2

32

0
107

1011

108

109

16
32

106
32

1010 1011

C

C

Figure 3: The relative L2 -difference ∆ vs. complexity C for the scalar flux Φ in the line source problem, using the hybrid method
(filled, colored markers) and the monolithic SN method (empty, green markers). Points that are down and to the left are more
efficient. The hybrid solver was run for N = 4 (left) and N = 8 (right). The shape of data-points indicates different spatial
resolution (square: Nx = 51, circle: Nx = 101, triangle: Nx = 201). Coloring of the hybrid data points corresponds to the total
tot according to the colorbar. The S
number of particles per time step NM
N method was run for N = 4, 8, 16, 32. Each DG-data
C
point is assigned a numerical label according to its value of N . The formula for ∆ is given in (48) which the complexity C is
given by (49a) for the monolithic SN method and by (49b) for the hybrid.

employs an SN method with N ∈ {4, 8}. In the hybrid method the number of particles is also changed
between runs. The number of particles newly inserted into the system every time step is roughly 2 × Np
where Np = 10k for k ∈ {2, 3, 4, 5, 6}. The killing weight is fixed at wkill = 10−15 . All runs are performed
to a final time of tfinal = 3.2 and the CFL is kept fixed at 25.6. The reference solution is a S96 -S16 hybrid
(meaning a S96 for the uncollided component and S16 for the collided component) using a TN quadrature in
angle, a third-order DG method in space on a 448 by 448 grid, and a defect correction time integrator [10].
Selected results for the scalar flux Φ are depicted in Figure 5, and the relative L2 -error for these same
solutions is depicted in Figure 6. While the hybrid solutions are not completely free of ray-effects, these effects
are much more pronounced in the SN runs. A more rigorous analysis of the performance of the two algorithm
in dependence of their respective parameters can be seen in Figure 7. This plot shows the relative error ∆ of
various runs in dependence of their complexity C. It is noteworthy that for this test problem the accuracy is
mostly independent of the angular resolution, but depends significantly on the spatial resolution. Increasing
the overall particle count in the hybrid method is most effective at higher spatial resolution; compare for
example the vertical separation in colored triangles vs. colored circles vs. colored squares in Figure 7.
It turns out that increasing the number of particles in the hybrid algorithm does not necessarily increase
the algorithm complexity. This is because with increased particle count the iterative solver for the collided
components often needs fewer iterations. In cases where the complexity is dominated by these iterations, an
increase in particles can even cause a decrease in computational complexity. Overall, Figure 7 shows that the
hybrid algorithm produces results with comparable or slightly better accuracy than the standard SN solver,
while being close to an order of magnitude of lower complexity.

15

7
1

blue(filled)
red(striped)
white

1

σa
10
0
0

σs
0
1
1

Q
0
1
0

7

Figure 4: Geometric layout and table of material properties for the Lattice Problem.

10 -1

10 -2

(a) Reference

(b) S4 : Nx = 112, ∆ =
13%, C = 2.2 × 107

(c) S4 : Nx = 224, ∆ =
7%, C = 1.4e8

(d) HMC-S4 :Nx = 224,
N tot = 1.2 × 106 , ∆ =

(e) HMC-S8 : Nx = 112,
N tot = 3.2 × 106 , ∆ =

(f) HMC-S4 : Nx = 224,
N tot = 1.3 × 107 , ∆ =

10%, C = 8.3 × 107

13%, C = 1.8 × 107

3%, C = 1.0 × 108

10 -3

10 -4

10 -5

10 -6

MC

MC

MC

Figure 5: Numerical solutions to the lattice problem at t = 3.2 with CFL 25.6. Each numerical solution is characterized by a
relative L2 difference ∆ with respect to the reference, defined in (48), and a complexity C, defined in (49a) for the monolithic
SN method and (49b) for the hybrid.

16

10 -1

10 -2

(a) Reference

(b) S4 : Nx = 112, ∆ =
13%, C = 2.2 × 107

(c) S4 : Nx = 224, ∆ =
7%, C = 1.4 × 108

(d) HMC-S4 : Nx = 224,
N tot = 1.2 × 106 , ∆ =

(e) HMC-S8 : Nx = 112,
N tot = 3.2 × 106 , ∆ =

(f) H MC-S4 : Nx = 224,
N tot = 1.3 × 107 , ∆ =

10%, C = 8.3 × 107

13%, C = 1.8 × 107

3%, C = 1.0 × 108

10 -3

10 -4

10 -5

10 -6

MC

MC

MC

Figure 6: Lattice problem: Absolute difference between the reference solution and various numerical solutions at t = 3.2 with
CFL 25.6. Each numerical solution is characterized by a relative L2 difference ∆ with respect to the reference, defined in (48),
and a complexity C, defined in (49a) for the monolithic SN method and (49b) for the hybrid.

17

•

 :HMC-S ,Nx = 56 :HMC-S ,Nx = 112 :HMC-S ,Nx = 224
n : SN , Nx = 56
n : SN , Nx = 112 △
n : SN ; Nx = 224



N

N

N

0.3
(N = 4)
4 8

0.25
16 32

∆

0.2

16 32

0.2

0.15

4

105

0.15
8

4

16 32

0.1

8

16 32

104

8

103

0.1
4

0.05
0
106

106

(N = 8)
4 8

tot
NM
C

0.25

0.3

107

108

8

109

16

4
0.05
0
106

1010

C

107

108

109

16

1010

102

C

Figure 7: The relative L2 -difference ∆ vs. complexity C for the scalar flux Φ in the lattice problem, using the hybrid method
(filled, colored markers) and the monolithic SN method (empty, green markers). Points that are down and to the left are more
efficient. The hybrid-solver was run for N = 4 (left) and N = 8 (right). The shape of data-points indicates different spatial
resolution (square: Nx = 56, circle: Nx = 112, triangle: Nx = 224). Coloring of the Hybrid-data points corresponds to the total
tot according to the colorbar. The S
number of particles per time step NM
N method was run for N = 4, 8, 16, 32. Each DG-data
C
point is assigned a numerical label according to its value of N . The formula for ∆ is given in (48) which the complexity C is
given by (49a) for the monolithic SN method and by (49b) for the hybrid.

18

3.3. The linearized hohlraum problem
In the linearized hohlraum problem [9], nonlinear coupling between particles and the material medium
is approximated in a linear way by adjusting the absorption and scattering cross-sections according to the
expected material temperature profile of the nonlinear problem [8]. The geometry of the setup along with the
material parameters can be found in Figure 8. The domain is X = [0, 1.3] × [0, 1.3], and the initial condition
is identically zero everywhere. For boundary conditions, we assume a constant influx from the left side of the
domain, i.e. Ψ(x = 0, y, Ω, t) = 1 for Ωx > 0. As discussed
in the appendix, this boundary condition can be
√
treated as a surface source, modeled by setting Ωx = ξ, where ξ ∼ U ([0, 1]) is sampled uniformly on [0, 1].
The spatial distribution along the boundary is sampled uniformly.
We again perform SN and hybrid runs with varying spatial and angular resolution. The spatial domain is
subdivided into equal Nx ×Nx square cells with Nx ∈ {52, 104, 208}. For the SN -runs we use N ∈ {4, 8, 16, 32},
resulting in NΩ = N 2 ordinates on the northern hemisphere of S2 . The collided part of the hybrid algorithm
employs an SN method with N ∈ {4, 8}. In the hybrid method the number of particles is also changed
between runs, but the killing weight remains fixed at wkill = 10−15 . The number of particles newly inserted
into the system every time step is 2 × 10k for k ∈ {2, 3, 4, 5, 6}. All runs are performed to a final time of
tfinal = 2.6 and the CFL is kept fixed at 52. The reference solution is a S96 -S16 hybrid using a TN quadrature
in angle, a third-order DG method in space on a 448 × 448 grid, and a defect correction time integrator [10].
In Figure 9, we show densities of a few select runs, calculated using SN and hybrid methods. The log
of the corresponding relative L2 -errors are depicted in Figure 10. As before, the SN solutions suffer from
ray-effects that are marginally reduced as the number of angle increases. Meanwhile, most of the disparities
between the reference and hybrid solutions can be attributed to stochastic noise. The hybrid has a mix of ray
effects from the collided equation solve and particle tracks from the uncollided equation solve on the backside
of the hohlraum However, these errors here are on the order of 10−2 -10−3 , which is much smaller than the
errors in othe back of the domain.
Detailed comparisons between the relative error against the computational complexity are depicted in
Figure 11. As before, increasing the angular resolution in the collided equation does not benefit the accuracy
of the hybrid method. Increasing particles also has less effect than in the previous problems. The SN solutions
benefit most from finer spatial resolution, while the angular resolution does not matter as much. Spatial
resolution also plays the biggest role for the hybrid. We do observe that for small particle counts (the purple
points in the figure) increasing resolution can actually increase the error. This is explained by the fact that
an under-sampled MC calculation does not benefit from more spatial resolution. Nevertheless, for a fixed
spatial resolution, we do observe an improvement in the error when Np is increased.
Overall the hybrid runs produce solutions with comparable or better accuracy than monolithic SN runs
with the same spatial resolution. Hybrid runs using S8 for the collided equation achieve improved accuracy
at nearly the same complexity while runs using S4 for the collided equation achieved improved accuracy with
even less complexity.
4. Conclusion and Discussion
In this work, we have presented a collision-based hybrid method that uses a Monte Carlo method for
the uncollided solution and a discrete ordinate discretization for the collided solution. This combination of
methods was originally proposed for the first collision source strategy used in [3] in the steady-state setting.
Thus this work can be considered as an extension to the time-dependent setting that requires a remap
procedure after every time step.
Experimental simulations have been performed on three standard benchmarks. For each benchmark, the
results demonstrate that the hybrid method is more efficient, in the sense that it achieves greater accuracy
with the comparable or less complexity or is less complexity with comparable or greater accuracy. Here
complexity is a measure of how many unknowns are updated during particle moves for Monte Carlo or
sweeping iterations for discrete ordinates.
This work has concentrated on single-energy particle transport problems. However recent work has
shown that when considering energy-dependent problems, more opportunities for hybridization arise when
considering fully deterministic hybrids [40]. In those results it was shown that low-resolution in energy can be
used for the collided solution as well as low-resolution in angle. With the introduction of Monte Carlo, new
opportunities arise. For example, continuous energy cross-sections could be used in the uncollided portion.
19

1.3
0.4

0.4

0.8

black
red(thin stripes)
green (thick stripes)
blue (checkered)
white

0.7

σa
0
5
10
50
0

σs
100
95
90
50
0.1

Q
0
0
0
0
0

1.3

0.2

Figure 8: Geometric layout and table of material parameters for the Hohlraum Problem. All walls have a thickness of 0.05.

This could be important to treat resonances in neutron transport problems, but investigation is needed to
quantify any benefits from this approach. In addition future investigations should be made regarding the use
of Monte Carlo techniques inside the high-order time accuracy methods developed for hybrid problems in
[9, 10].
5. Acknowledgments
J.K. gratefully acknowledges support from the 2022 National Science Foundation Mathematical Sciences
Graduate Internship to conduct this research at Oak Ridge National Laboratory.

20

(a) Reference

(b) S8 : Nx = 104, ∆ =
18%, C = 4.9 × 108

(c) S16 : Nx = 208, ∆ =
11%, C = 1.1 × 1010

(d) HMC-S4 :Nx = 52,
N tot = 2.1e6, ∆ = 18%,

(e) HMC-S8 :Nx = 104,
N tot = 7.6 × 106 , ∆ =

(f) H MC-S4 : Nx = 208,
N tot = 2.9 × 107 , ∆ =

C = 1.7 × 107

11%, C = 9.7 × 107

8%, C = 5.4 × 108

MC

MC

MC

Figure 9: Numerical solutions to the hohlraum problem at t = 2.6 with CFL 52. Each numerical solution is characterized by a
relative L2 difference ∆ with respect to the reference, defined in (48), and a complexity C, defined in (49a) for the monolithic
SN method and (49b) for the hybrid.

21

(a) Reference

(b) S8 : Nx = 104, ∆ =
18%, C = 9.8 × 108

(c) S16 : Nx = 208, ∆ =
11%, C = 1.7 × 1010

(d) HMC-S4 :Nx = 52,
N tot = 2.1 × 106 , ∆ =

(e) HMC-S8 :Nx = 104,
N tot = 7.6 × 106 , ∆ =

(f) H MC-S4 : Nx = 208,
N tot = 2.9 × 107 , ∆ =

18%, C = 1.1 × 107

11%, C = 9.7 × 107

8%, C = 5.4 × 108

MC

MC

MC

Figure 10: Hohlraum problem: Absolute difference between analytical solution to the line source problem and various numerical
solutions at t = 2.6 with CFL 52. Each numerical solution is characterized by a relative L2 difference ∆ with respect to the
reference, defined in (48), and a complexity C, defined in (49a) for the monolithic SN method and (49b) for the hybrid.

22

•

:HMC-S , Nx = 52 :HMC-S ,Nx = 104 :HMC-S ,Nx = 208
n : SN , Nx = 52
n : SN , Nx = 104 △
n : SN ; Nx = 208




0.5

N

N

0.5

(N = 4)

0.4

0.4

∆

106

(N = 8)

4

0.3

8
0.2
0.1
107

108

1010

0.2
0.1
107

1011

108

C

4

4
8

16 32
8
8
16 32
16
109

4

0.3

4

4

105

104

16 32
8
8
16 32
16
109

1010

tot
NM
C

N

103

1011

102

C

Figure 11: The relative L2 -difference ∆ vs. complexity C for the scalar flux Φ in the hohlraum problem, using the hybrid method
(filled, colored markers) and the monolithic SN method (empty, green markers). Points that are down and to the left are more
efficient. The hybrid-solver was run for N = 4 (left) and N = 8 (right). The shape of data-points indicates different spatial
resolution (square: Nx = 52, circle: Nx = 104, triangle: Nx = 208). Coloring of the Hybrid-data points corresponds to the total
tot according to the colorbar. The S
number of particles per time step NM
N method was run for N = 4, 8, 16, 32. Each DG-data
C
point is assigned a numerical label according to its value of N . The formula for ∆ is given in (48) which the complexity C is
given by (49a) for the monolithic SN method and by (49b) for the hybrid.

Appendix A. Boundary conditions for the hohlraum problem
Unlike the line source and lattice problems, the linearized hohlraum problem involves non-zero boundary
conditions. A Monte Carlo implementation of this boundary condition is stated in [15]; here we present a
derivation of the approach that is used.
In the hohlraum problem, it is assumed that Ψ(x, y, t, Ω) = 1 for x = 0 and Ωx > 0, i.e., a constant flux of
1 along the left boundary is assumed for each incoming direction. To model this with Monte Carlo, we assume
that this flux is due to a source sb (see (37c)) located on an infinitesimal slab just left of the boundary.
Consider first a finite slab Sa = {(x, y) ∈ [−a, 0] × [0, 1.3]}, where a > 0. We assume that σa = σs = 0 on
S, that the source sb (x, y, Ω; a) = sb (x, Ω; a) is independent of y and t, and that
Z 0
sb (x, Ω; a)dx
(A.1)
ŝb (Ω) :=
−a

is independent of a. Thus, sb (x, Ω, t) → ŝb (Ω)δ(x) as a → 0.
To determine ŝb , we assume that Ψ is independent of y and t on S and satisfies the steady-state equation
(x, y) ∈ S, Ω ∈ S2 ,
y ∈ [0, 1.3], Ωx > 0.

Ωx Ψx (x, y, Ω) = sb ,
Ψ(−a, y, Ω) = 0,

(A.2a)
(A.2b)

This formulation is consistent with (34), given the assumptions made on Ψ, sb , and the material cross-sections.
Integrating (A.2b) with respect to x and applying the boundary condition in (A.2b) gives
ŝb (Ω) = Ωx ,
23

Ωx > 0.

(A.3)

Hence sb = Ωx δ(x).
Finally, to sample particles according to the probability density ρ(Ωx ) ∝ ŝb for Ωx > 0, we compute the
cumulative distribution function (CDF):
F (Ωx ) =

Z Ωx
0

2µdµ = Ω2x .

(A.4)

According to the fundamental theorem of simulation [30, pp. 19-22], the correct√angular distribution can be
sampled by generating uniform variables u ∈ [0, 1] and setting Ωx = F −1 (u) = u.
References
[1] Marvin L Adams. Discontinuous finite element transport solutions in thick diffusive problems. Nuclear
science and engineering, 137(3):298–333, 2001.
[2] Timo Aila and Samuli Laine. Understanding the efficiency of ray traversal on gpus. In Proceedings of
the conference on high performance graphics 2009, pages 145–149, 2009.
[3] Raymond E Alcouffe. A first collision source method for coupling monte carlo and discrete ordinates
for localized source problems. In Monte-Carlo Methods and Applications in Neutronics, Photonics and
Statistical Physics: Proceedings of the Joint Los Alamos National Laboratory-Commissariat à l’Energie
Atomique Meeting Held at Cadarache Castle, Provence, France April 22–26, 1985, pages 352–366.
Springer, 2006.
[4] Kendall Atkinson. Numerical integration on the sphere. The ANZIAM Journal, 23(3):332–347, 1982.
[5] Yousry Azmy, Enrico Sartori, Edward W Larsen, and Jim E Morel. Advances in discrete-ordinates
methodology. Nuclear computational science: A century in review, pages 1–84, 2010.
[6] Guillaume Bal, Anthony B Davis, and Ian Langmore. A hybrid (monte carlo/deterministic) approach
for multi-dimensional radiation transport. Journal of Computational Physics, 230(20):7723–7735, 2011.
[7] Troy L Becker, Allan B Wollaber, and Edward W Larsen. A hybrid monte carlo–deterministic method
for global particle transport calculations. Nuclear science and engineering, 155(2):155–167, 2007.
[8] Thomas A Brunner. Forms of approximate radiation transport. Technical Report SAND2002-1778,
Sandia National Laboratories, Albuquerque, NM (United States), 2002.
[9] Michael M Crockatt, Andrew J Christlieb, C Kristopher Garrett, and Cory D Hauck. An arbitrary-order,
fully implicit, hybrid kinetic solver for linear radiative transport using integral deferred correction.
Journal of Computational Physics, 346:212–241, 2017.
[10] Michael M Crockatt, Andrew J Christlieb, C Kristopher Garrett, and Cory D Hauck. Hybrid methods
for radiation transport using diagonally implicit runge–kutta and space–time discontinuous galerkin time
integration. Journal of Computational Physics, 376:455–477, 2019.
[11] Michael M Crockatt, Andrew J Christlieb, and Cory D Hauck. Improvements to a class of hybrid
methods for radiation transport: Nyström reconstruction and defect correction methods. Journal of
Computational Physics, 422:109765, 2020.
[12] Jeffery D Densmore and Edward W Larsen. Asymptotic equilibrium diffusion analysis of time-dependent
Monte Carlo methods for grey radiative transfer. Journal of Computational Physics, 199(1):175–204,
2004.
[13] Jeffery D Densmore, Todd J Urbatsch, Thomas M Evans, and Michael W Buksas. A hybrid transportdiffusion method for monte carlo radiative-transfer simulations. Journal of Computational Physics,
222(2):485–503, 2007.

24

[14] Stephen A Dupree and Stanley K Fraley. A Monte Carlo primer: A Practical approach to radiation
transport, volume 1. Springer Science & Business Media, 2002.
[15] Joseph A Fleck Jr. The calculation of nonlinear radiation transport by a monte carlo method. Technical
report, Lawrence Livermore National Lab.(LLNL), Livermore, CA (United States), 1961.
[16] BD Ganapol. Homogeneous infinite media time-dependent analytic benchmarks for x-tm transport
methods development. Los Alamos National Laboratory, 1999.
[17] C Kristopher Garrett and Cory D Hauck. A comparison of moment closures for linear kinetic transport
equations: The line source benchmark. Transport Theory and Statistical Physics, 42(6-7):203–235, 2013.
[18] Jean-Luc Guermond and Guido Kanschat. Asymptotic analysis of upwind discontinuous galerkin
approximation of the radiative transport equation in the diffusive limit. SIAM Journal on Numerical
Analysis, 48(1):53–78, 2010.
[19] Weimin Han, Jianguo Huang, and Joseph A Eichholz. Discrete-ordinate discontinuous galerkin methods
for solving the radiative transfer equation. SIAM Journal on Scientific Computing, 32(2):477–497, 2010.
[20] Cory D Hauck and Ryan G McClarren. A collision-based hybrid method for time-dependent, linear,
kinetic transport equations. Multiscale Modeling & Simulation, 11(4):1197–1227, 2013.
[21] Vincent Heningburg and Cory D Hauck. Hybrid solver for the radiative transport equation using finite
volume and discontinuous galerkin. arXiv preprint arXiv:2002.02517, 2020.
[22] Edward W Larsen, Jim E Morel, and Warren F Miller Jr. Asymptotic solutions of numerical transport
problems in optically thick, diffusive regimes. Journal of Computational Physics, 69(2):283–324, 1987.
[23] Elmer Eugene Lewis and Warren F Miller. Computational methods of neutron transport. 1984.
[24] Ryan McClarren. Computational nuclear engineering and radiological science using python. Academic
Press, 2017.
[25] Ryan G McClarren and Todd J Urbatsch. A modified implicit monte carlo method for time-dependent
radiative transfer with adaptive material coupling. Journal of Computational Physics, 228(16):5669–5686,
2009.
[26] Scott W Mosher and Farzad Rahnema. The incident flux response expansion method for heterogeneous
coarse mesh transport problems. Transport Theory and Statistical Physics, 35(1-2):55–86, 2006.
[27] Madicken Munk and Rachel N Slaybaugh. Review of hybrid methods for deep-penetration neutron
transport. Nuclear Science and Engineering, 2019.
[28] H Park, DA Knoll, RM Rauenzahn, AB Wollaber, and RB Lowrie. Moment-based acceleration of monte
carlo solution for multifrequency thermal radiative transfer problems. Journal of Computational and
Theoretical Transport, 43(1-7):314–335, 2014.
[29] Zhuogang Peng and Ryan G McClarren. A high-order/low-order (holo) algorithm for preserving
conservation in time-dependent low-rank transport calculations. Journal of Computational Physics,
447:110672, 2021.
[30] Steve Pressé and Ioannis Sgouralis. Data Modeling for the Sciences: Applications, Basics, Computations.
Cambridge University Press, 2023.
[31] Qiwei Sheng and Cory Hauck. Uniform convergence of an upwind discontinuous galerkin method for
solving scaled discrete-ordinate radiative transfer equations with isotropic scattering. Mathematics of
Computation, 90(332):2645–2669, 2021.
[32] Yi Shi, Xiaole Han, Wenjun Sun, and Peng Song. A continuous source tilting scheme for radiative transfer
equations in implicit monte carlo. Journal of Computational and Theoretical Transport, 50(1):1–26, 2020.
25

[33] Jerome Spanier and Ely M Gelbard. Monte Carlo principles and neutron transport problems. Courier
Corporation, 2008.
[34] Elad Steinberg and Shay I Heizler. Multi-frequency implicit semi-analog Monte-Carlo (ISMC) radiative
transfer solver in two-dimensions (without teleportation). Journal of Computational Physics, 450:110806,
2022.
[35] CP Thurgood, A Pollard, and HA Becker. The tn quadrature set for the discrete ordinates method.
1995.
[36] John C Wagner and Alireza Haghighat. Automated variance reduction of Monte Carlo shielding
calculations using the discrete ordinates adjoint function. Nuclear Science and Engineering, 128(2):186–
208, 1998.
[37] John C Wagner, Douglas E Peplow, and Scott W Mosher. Fw-cadis method for global and regional
variance reduction of monte carlo radiation transport calculations. Nuclear Science and Engineering,
176(1):37–57, 2014.
[38] Wallace F Walters. Use of the chebyshev-legendre quadrature set in discrete-ordinate codes. In Proc.
Conf. Numerical Methods in High Temperature Physics, pages 2–6, 1987.
[39] Ben Whewell, Ryan G McClarren, Cory D Hauck, and Minwoo Shin. Multigroup neutron transport
using a collision-based hybrid method. Nuclear science and engineering, pages 1–20, 2023.
[40] Ben Whewell, Ryan G McClarren, Cory D Hauck, and Minwoo Shin. Multigroup Neutron Transport
Using a Collision-Based Hybrid Method. Nuclear science and engineering, 2023.
[41] Jeffrey Alan Willert. Hybrid deterministic/Monte Carlo methods for solving the neutron transport
equation and k-eigenvalue problem. North Carolina State University, 2013.
[42] Allan B Wollaber. Four decades of implicit Monte Carlo. Journal of Computational and Theoretical
Transport, 45(1-2):1–70, 2016.
[43] Dingkang Zhang and Farzad Rahnema. Comet solutions to a stylized bwr benchmark problem. Technical
report, American Nuclear Society, Inc., 555 N. Kensington Avenue, La Grange Park . . . , 2012.

26