# A Quasi-Monte Carlo Method with Krylov Linear Solvers for Multigroup Neutron
Transport Simulations

###### Abstract

In this work we investigate replacing standard quadrature techniques used in
deterministic linear solvers with a fixed-seed Quasi-Monte Carlo calculation
to obtain more accurate and efficient solutions to the neutron transport
equation (NTE). Quasi-Monte Carlo (QMC) is the use of low-discrepancy
sequences to sample the phase space in place of pseudo-random number
generators used by traditional Monte Carlo (MC). QMC techniques decrease the
variance in the stochastic transport sweep and therefore increase the accuracy
of the iterative method. Historically, QMC has largely been ignored by the
particle transport community because it breaks the Markovian assumption needed
to model scattering in analog MC particle simulations. However, by using
iterative methods the NTE can be modeled as a pure-absorption problem. This
removes the need to explicitly model particle scattering and provides an
application well-suited for QMC. To obtain solutions we experimented with
three separate iterative solvers: the standard Source Iteration (SI) and two
linear Krylov Solvers, GMRES and BiCGSTAB. The resulting hybrid iterative-QMC
(iQMC) solver was assessed on three one-dimensional slab geometry problems. In
each sample problem the Krylov Solvers achieve convergence with far fewer
iterations (up to 8x) than the Source Iteration. Regardless of the linear
solver used, the hybrid method achieved an approximate convergence rate of
$O(N^{-1})$, as compared to the expected $O(N^{-1/2})$ of traditional MC
simulation, across all test problems.

Sam Pasmanna Ilham Variansyaha C. T. Kelleyb Ryan McClarrena<EMAIL_ADDRESS>
aDepartment of Aerospace and Mechanical Engineering
University of Notre Dame
Fitzpatrick Hall, Notre Dame, IN 46556 bNorth Carolina State University,
Department of Mathematics
3234 SAS Hall, Box 8205
Raleigh NC 27695-8205

Neutron Transport Monte Carlo Methods Quasi Monte Carlo Krylov Linear Solvers

##  1 Introduction

Solving the neutron transport equation (NTE) under various conditions,
accurately, and efficiently is vital to nuclear reaction simulations like
those in advanced reactor design or accident analysis [1]. The neutron
transport equation describes the distribution of neutrons in space, angle,
energy, and time. The equation’s high-dimensional nature makes it difficult to
design efficient general-purpose algorithms and many solution techniques, most
common of which have been stochastic Monte Carlo simulations [2] or
deterministic discrete ordinates ($S_{N}$) methods [3], have been developed.

For deterministic solutions, Source Iteration (SI) is the simplest and most
common deterministic solution technique for solving the discrete ordinates
method [3]. SI is equivalent to a fixed-point Picard Iteration, nevertheless,
as problems become collision dominated, the convergence rate of the SI can
become arbitrarily slow [4, 5]. More advanced iteration techniques such as
Krylov subspace methods, including Generalized Minimal RESidual method (GMRES)
and BiConjugate Gradient STABilized method (BiCGSTAB), have been shown to
outperform standard Source Iteration, particularly when there are highly
scattering materials [4]. Nonetheless, as the dimensionality and fidelity of
the problem increases, the deterministic quadrature techniques used to
evaluate the system of equations become intractable [6, 7].

Monte Carlo (MC) simulations provide a more robust solution by using random
sampling and probability to produce solutions. In this method the statistical
error scales according to $O(N^{-1/2})$ — where $N$ is the number of neutron
histories — regardless of the dimensionality of the problem [8]. However, MC
simulations are often seen as a last resort due to their high computational
cost and slow rates of convergence [9, 10, 11]. Recent work by Willert et al.
investigated a hybrid MC-deterministic solution where the deterministic
quadrature sweep of the iterative method was replaced with a Monte Carlo
transport simulation [7, 12]. This method attempts to combine the efficiency
of iterative methods while also providing an accurate solution for complex
problems given the robustness of MC simulation. However, it is found that a
staggering number of particle histories are may be required for convergence of
the iterative method even for mono-energetic slab problems [12]. Or, as
Heinrich Von Kleist wrote in a previous age [13], “probability is not always
on the side of truth.”

Our work investigates the use of fixed-seed Quasi-Monte Carlo (QMC) techniques
in place of standard, pseudo-random MC to decrease the variance in the
transport process and therefore improve the convergence of the iterative
method. Quasi-Monte Carlo techniques use low-discrepancy sequences (LDS) in
place of typical pseudo-random number generators for Monte Carlo sampling.
Various LDS have been developed, including the Sobol and Halton sequences, the
goal of each is to sample the phase space in a deterministic and self-avoiding
manner. Theoretically, this results in a sampling convergence rate
proportional to $O(N^{-1})$, compared to $O(N^{-1/2})$ of standard Monte Carlo
[14].

Rather than taking subsequent samples from the same LDS at the start of every
iteration, the LDS is reset to the beginning of the sequence. This fixed-seed
approach allows the iterative method to converge at a much faster rate than if
new samples were taken. With typical random number generators, this technique
would be avoided to ensure samples are uniformly distributed through the
phase-space. However, the low-discrepancy nature of the Sobol Sequence, Halton
Sequence, etc. ensure a well-balanced sampling of the phase-space even with a
relatively low number of samples.

Despite the benefits offered by QMC, it has largely been ignored by the
particle transport community [15]. There has been some recent work in using
QMC for radiative transfer problems without scattering [16, 17], but, to the
knowledge of the authors, there has not been any recent work with QMC applied
to neutron transport. This is likely because the deterministic nature of the
LDS breaks the Markovian assumption needed for the particle random-walk when
scattering is present. Therefore, QMC must be implemented in applications
which are not Markovian processes or steps must be take to ensure the
Markovian assumption is held.

Presently, there have been two strategies for implementing QMC in particle
transport. The first is known generally as randomized-QMC or (RQMC) which
includes a host of strategies that attempt to randomize a sequence and still
retain the low-discrepancy of the samples [18, 15, 17, 19]. Depending on the
randomization technique, the theoretical convergence rate of $O(N^{-1})$ may
be reduced [20, 21]. While other randomization techniques may theoretically
uphold the $O(N^{-1})$ convergence rate, but are computationally expensive to
execute [21]. The second implementation of QMC in particle transport is simply
in use of problems without scattering, primarily seen in radiative heat
transfer problems [16, 17].

Our proposed iterative-QMC (iQMC) method allows for both: problems that
include scattering and the use of unaltered LDS, i.e., no RQMC method is
required. This is achieved by modeling the problem in the QMC simulation as a
purely absorbing system where each particle is emitted and traced out of the
volume. After this process, or QMC Sweep as it will be referred to from now
on, the scattering term is iterated upon using a deterministic linear solver
and the process repeats until a desired tolerance or maximum number of
iterations is reached, thereby removing the need for the simulation of a
random walk process [22].

The outline of this paper is as follows: Section 2.1 presents a brief overview
of the neutron transport equation, Source Iteration, and Krylov solvers.
Section 2.2 describes the use of fixed-seeding and low-discrepancy sequences
to form the Quasi-Monte Carlo transport sweep. Section 3 provides an overview
of the implemented algorithms before analysis and results from three 1-D test
problems are presented in Section 4. The first problem solves for scalar flux
in an infinite medium with multi-group data generated from FUDGE [23] with a
known analytic solution. The second problem, known as Reed’s Problem, is a
multi-media problem benchmarked with results from a high particle count Monte
Carlo simulation, using the Center for Exascale Monte-Carlo for Neutron
Transport’s (CEMeNT’s) Monte-Carlo Dynamic Code (MCDC) [24]. The third and
final problem, provided by Garcia et al., provides angular flux results at the
slab edges from a fixed boundary source with a spatially decaying scattering
cross section [25]. Finally, key findings and future work are discussed in
Section 5.

##  2 Methods

###  2.1 Neutron Transport Source Iteration and Krylov Methods

We begin with the one-speed neutron transport equation in slab geometry with
isotropic scattering [3]:

$\mu\frac{\partial\psi}{\partial
x}(x,\mu)+\Sigma_{t}(x)\psi(x,\mu)=\frac{1}{2}\left[\Sigma_{s}(x)\phi(x)+q(x)\right],$
(1) $\phi(x)=\int_{-1}^{1}\psi(x,\mu)d\mu,$ (2)

for $0\leq x\leq\tau$. The boundary conditions are

$\psi(0,\mu>0)=\psi_{l}(\mu),\qquad\psi(\tau,\mu<0)=\psi_{r}(\mu).$ (3)

Where $x,\mu$ are the particle position and angle respectively, $\psi$ is the
angular flux, $\phi$ is the scalar flux, $\Sigma_{t}$ is the total macroscopic
cross section, $\Sigma_{s}$ is the scattering macroscopic cross section, and
$q$ represents an internal source function.

####  2.1.1 Source Iteration

The transport Eqs. (1)–(2) can be solved iteratively via Source Iteration
(SI):

$\mu\frac{\partial\psi^{(n+1)}}{\partial
x}(x,\mu)+\Sigma_{t}(x)\psi^{(n+1)}(x,\mu)=\frac{1}{2}\left[\Sigma_{s}(x)\phi^{(n)}(x)+q(x)\right],$
(4) $\phi^{(n+1)}(x)=\int_{-1}^{1}\psi^{(n+1)}(x,\mu)d\mu,$ (5)

where superscript (n) indicates iteration index. The equations can be
represented in operator notation as:

$\phi^{(n+1)}={\cal S}[\phi^{(n)},q,\psi_{l},\psi_{r}],$ (6)

where the transport sweep operator ${\cal S}$ updates a scalar flux estimate
given an internal source $q$ and boundary sources $\psi_{l}$ and $\psi_{r}$.
In the proposed hybrid method, this transport sweep operation is performed via
Monte Carlo simulation (which is discussed later in Sec. 2.2).

The SI Equation (6) can be rewritten as follows:

$\phi^{(n+1)}={\cal K}[\phi^{(n)}]+f,$ (7)

where

${\cal K}[\phi^{(n)}]={\cal S}[\phi^{(n)},0,0,0],$ (8)

and

$f={\cal S}[0,q,\psi_{l},\psi_{r}].$ (9)

By collecting the scalar flux terms, one can demonstrate that SI is equivalent
to the fixed-point Picard iteration of a linear problem

$A\phi\equiv(I-{\cal K})\phi=f,$ (10)

where $I$ is the identity function. Equation (10) is in a form we can send to
linear solvers, particularly those that are more efficient than the fixed-
point Picard iteration, such as Krylov methods. Note that we do not need to
explicitly form the matrix $A$, we only need to compute the action of $A$ on
$\phi$, which is accomplished with the QMC Sweep.

####  2.1.2 Krylov Methods

An order-$r$ Krylov subspace is defined with notation from the previous
section as [26]:

$K_{r}=\textrm{span}(\phi,A\phi,A^{2}\phi,...,A^{r-1}\phi).$ (11)

For each experiment presented in Sec. 4, two Krylov methods, GMRES [27] and
BiCGSTAB [28], were used. The Generalized Minimum RESidual (GMRES) is one of
the most common Krylov methods. When solving $A\vec{\phi}=\vec{f}$, GMRES
minimized $||f-A\phi||_{2}$ over the $k^{th}$ Krylov subspace. For every
iteration, the GMRES stores an additional Krylov vector. For problems that
require many iterations this may lead to memory constraints. BiCGSTAB is a
low-storage Krylov method that is memory bounded throughout the algorithm.
However, the memory savings come from information that is thrown out with each
iteration and therefore BiCGSTAB will generally require more iterations to
converge than GMRES. Nonetheless, as we will observe in Sec. 4, both Krylov
methods will require far fewer iterations than the SI.

###  2.2 Quasi-Monte Carlo Transport Sweep

####  2.2.1 Monte Carlo Transport Sweep

Monte Carlo methods for neutron transport seek to simulate the behavior of a
statistically significant number of particles from birth to death to gain an
approximate behavior of the system. For our one-dimensional simulations, each
particle begins with an initial position ($x_{i}$), direction ($\mu_{i}$), and
statistical weight ($w_{i}$). In an analog simulation, the particle would then
be tracked from collision to collision, tallying quantities of interest as the
particle moves. Each time the particle undergoes a scattering collision, a new
direction ($\mu$) would be sampled and the next distance to collision would be
calculated. This process would repeat until the particle is either absorbed or
exits the volume. However, Eq. (4) is a purely absorbing transport problem
with a known source. MC simulation in a purely absorbing system can be
enhanced by employing the continuous weight absorption technique (also called
implicit capture), which continuously reduces the statistical weight of each
particle per length traveled ($s$):

$w_{\textrm{new}}=w_{\textrm{old}}e^{-\Sigma_{a}s}.$ (12)

Consequently, after emission the particle is traced straight out of the volume
reducing the statistical weight according to the distance traveled across each
spatial cell, as illustrated by Figure 1.

Figure 1:  A simplified diagram of a Monte Carlo sweep. For each MC sweep the
particles are emitted with initial position and angle $(x_{i},\mu_{i})$ and
are swept out of the volume, tallying the scalar flux according to the path
length tally estimator (Equation 13).

We use the track-length tally estimator to compute the spatially-averaged
scalar flux in the defined mesh. Because the weight is continuously reduced
with each step, the tally scoring becomes:

$\frac{1}{V}\int_{0}^{s}w_{\textrm{old}}e^{-\Sigma_{a}s\prime}ds=\frac{w_{\textrm{old}}}{V}\left(\frac{1-e^{-\Sigma_{a}s}}{\Sigma_{a}}\right).$
(13)

This estimate of the scalar flux is then used to compute the scattering source
in the next iteration of the solver.

####  2.2.2 Low-Discrepancy Sequence for Quasi-Monte Carlo

Given the purely absorbing system (Equation 4) and the use of continuous
weight absorption technique (Equations 12 and 13), the only things need to be
randomly sampled are the particle initial position $x_{i}$ and direction
$\mu_{i}$. In a standard MC transport sweep, a pseudo-random number generator
is used to sample $x_{i}$ and $\mu_{i}$. In Quasi-Monte Carlo transport sweep,
a quasi-random low-discrepancy sequence is used instead.

Low-discrepancy sequences use deterministic algorithms to sample the phase
space in a self-avoiding manner thereby approaching a more uniform
distribution and approximating the expectation more efficiently. This results
in a theoretical convergence rate of $O(N^{-1})$ compared to the $O(N^{-1/2})$
from pseudo-randomly placed points [21]. In addition to the well-known Sobol
and Halton Sequences, Figure 2 shows the distribution of 256 points in a unit
square of a newer LDS known as the Golden Sequence [29].

Figure 2:  256 points generated in a unit square with pseudo random points
(top left), Sobol Sequence [30] (top right), Halton Sequence [31] (bottom
left), and Golden Sequence [29] (bottom right).

####  2.2.3 Multigroup Vectorization

Neutron cross sections vary greatly with energy and contain large resonance
regions, making them computationally expensive to model with high fidelity.
The multigroup method is a common approach used to model energy-dependent
cross sections that splits the energy range into $G$ finite regions, each with
a representative cross section. The multigroup equations for $G$ groups in 1D
slab geometry in matrix form are [3].

$\mu\frac{\partial\vec{\psi}}{\partial
x}(x,\mu)+\underline{\Sigma}_{t}(x)\vec{\psi}(x,\mu)=\frac{1}{2}\underline{\Sigma}_{s}(x)\int_{-1}^{1}\vec{\psi}(x,\mu^{\prime})\,d\mu^{\prime}+\frac{\vec{q}(x)}{2},$
(14)

where

$\vec{\psi}=(\psi_{1},\psi_{2},\dots,\psi_{G})^{\mathrm{T}},\qquad\vec{q}=(q_{1},q_{2},\dots,q_{G})^{\mathrm{T}},$
(15) $\underline{\Sigma}_{t}(x)=\begin{pmatrix}\Sigma_{t,1}(x)&0&\dots&0\\\
0&\Sigma_{t,2}(x)&&0\\\ \vdots&&\ddots&\vdots\\\
0&\dots&0&\Sigma_{t,G}(x)\end{pmatrix},$ (16)

and

$\underline{\Sigma}_{s}(x)=\begin{pmatrix}\Sigma_{s,1\rightarrow
1}(x)&\Sigma_{s,2\rightarrow 1}(x)&\dots&\Sigma_{s,G\rightarrow 1}(x)\\\
\Sigma_{s,2\rightarrow 1}(x)&\Sigma_{s,2\rightarrow
1}(x)&\dots&\Sigma_{s,G\rightarrow 2}(x)\\\ \vdots&\vdots&&\vdots\\\
\Sigma_{s,G\rightarrow 1}(x)&\Sigma_{s,G\rightarrow
1}(x)&\dots&\Sigma_{s,G\rightarrow G}(x)\\\ \end{pmatrix}.$ (17)

As a consequence of this energy discretization and the employed scattering-
free particle tracing technique, each particle in the MC sweep can now
represent all energy groups. Conversely, in analog MC each simulated particle
may only represent one energy group and random samples are taken to determine
up or down scattering after collision. In our hybrid QMC-iterative method, the
statistical weight of each particle is now a vector of weights, similar to
Equation 15:

$\vec{w}=(w_{1},w_{2},\dots,w_{G})^{T}.$ (18)

Therefore, we only need to multiply this vector by the scattering cross
section (Eq. 17) to determine the scattering distribution.

##  3 Implementation Details

The algorithm was written in Julia, a scientific computing language that
combines the compiler capabilities of C++ and the syntax of Matlab and Python.
The code and primary documentation are available here iQMC.jl [32]. The Krylov
linear solvers come from the Julia package SIAMFANLEQ.jl [33]. The
documentation for these codes is in the Julia notebooks [34] and the book [35]
that accompany the package.

The results in Section 4 were generated using the Sobol Sequence as the LDS in
QMC Sweep. The Sobol Sequence generates nets with $2^{m}$ points and loses
some of its balance properties if generated with a sample size that is not a
power of 2 [36]. Consequently, the number of particles $N$, varied by powers
of 2 in all experiments. A brief convergence comparison between the Sobol,
Halton, and Golden sequences can be seen in Problem 4.2, Figure 9. Finally,
for an equivalent comparison to MC, the LDS was replaced with a fixed-seed
pseudo-random number generator. By fixed-seed, we refer to a calculation where
the random number seed is reset at the beginning of each iteration: this
assures that the stochasticity of the MC algorithm is not an impediment to
convergence. An outline of the QMC Sweep algorithm is described below
(Algorithm 1).

Algorithm 1 QMC Sweep ($\phi_{\textrm{in}}$)

1: Initialize Low-Discrepancy-Sequence (LDS)

2: for $i$ in $N$ do

3: Assign position and angle ($x_{i}$, $\mu_{i}$) based on the LDS

4: Initialize weight
$w_{i}=\left(\Sigma_{s}\phi_{\textrm{in}}+q\right)N_{x}V/N$

5: for $j$ in $N_{x}$ do

6: Move particle across $\textrm{Zone}_{j}$

7: Tally($x,\mu,w$) (Eq: 13)

8: Update particle weight (Eq. 12)

9: end for

10: end for

11: Return: $\phi_{\textrm{out}}$

As previously mentioned in Section 2, the described QMC-iterative algorithm
does not require the explicit formation of the matrix $A$ for the linear
solver. Instead, we compute the matrix vector product $A\phi$, as described in
Algorithm 2, from which the Krylov solvers can iterate.

Algorithm 2 Matrix Vector Product ($\phi_{\textrm{in}}$)

1: $\textrm{b}=\textrm{QMC Sweep}(\vec{0})$

2: $\textrm{mxv}=\textrm{QMC Sweep}(\phi_{\textrm{in}})$

3: $\textrm{axv}=\phi_{\textrm{in}}-\textrm{mxv}-\textrm{b}$

4: Return: axv

To evaluate the effect of varying the number of spatial cells ($N_{x}$) within
a QMC simulation, a post-process spatial averaging technique was developed.
Given a reference solution for the scalar flux with $N_{x\textrm{Ref}}$
spatial cells, experiments were run for problems 4.1 and 4.2 with
$N_{x}=(N_{x\textrm{Ref}})^{2^{n}}$. The resulting vectors can be averaged
$\left[\log(N_{x}/N_{x\textrm{Ref}})/\log(2)\right]$ times to reduce all
vectors to the length of the reference solution as seen in Algorithm 3.

Algorithm 3 Flux Spatial Average ($\phi_{\textrm{in}}$, $N_{x\textrm{Ref}}$)

1: $\phi_{\textrm{out}}=\phi_{\textrm{in}}$

2: $N_{x}=\textrm{length}(\phi_{\textrm{out}})$

3: $I=\log(N_{x}/N_{x\textrm{Ref}})/\log(2)$

4: for $i=1$ to $I$ do

5: $\textrm{LeftCells}=\phi_{\textrm{out}}[1:2:N_{x}-1]$

6: $\textrm{RightCells}=\phi_{\textrm{out}}[2:2:N_{x}]$

7: $\phi_{\textrm{out}}=(\textrm{RightCells}+\textrm{LeftCells})*0.5$

8: $N_{x}=\textrm{length}(\phi_{\textrm{out}})$

9: end for

10:
Assert$\left(\textrm{length}(\phi_{\textrm{out}})=\textrm{length}(N_{x\textrm{Ref}})\right)$

11: Return: $\phi_{\textrm{out}}$

##  4 Computational Results

###  4.1 Problem 1: Multigroup in an Infinite Medium

The first problem features 12, 70, and 618 group cross section data of high-
density polyethylene (HDPE) generated with FUDGE [23] in an infinite medium.
For a constant volumetric source $Q$, the analytic solution for the scalar
flux is given by:

$\phi^{\textrm{Sol}}=(\Sigma_{t}-\Sigma_{s})^{-1}\cdot Q.$ (19)

To simulate an infinite medium, we placed an isotropic boundary source on each
slab edge, where the source strength was held at the expected analytic
solution. Figure 3 shows the center of each energy group for the total cross
section against the analytic solution for scalar flux from Eq. (19) divided by
the energy bin width.

From the convergence of relative residuals in Figure 4, it is observed that
the Krylov methods require far fewer transport sweeps to achieve the same
levels of convergence as the Source Iteration regardless of the number of
groups. Figure 5 plots the residual,

$R=\frac{||\phi^{\textrm{Sol}}-\phi^{QMC}||_{\infty}}{\phi^{\textrm{Sol}}},$
(20)

against the number of particle histories $N$ for a given number of spatial
cells $N_{x}$ along with the theoretical convergence goal of $O(N^{-1})$ for
QMC and $O(N^{-1/2})$ for MC.

Simulations were run with 80, 160, and 320 spatial cells and afterward results
were run through the spatially averaged scalar flux post-processing technique
from Algorithm 3 . This allows for comparison of all results to a solution
$\phi^{\textrm{Sol}}$ with $N_{x}=80$. Because the solution is spatially
constant per energy group, the results converge at the approximate rate of
$O(N^{-1})$ for QMC and $O(N^{-1/2})$ for MC, regardless of the number of
spatial cells (Figure 5).

Figure 3:  Energy group centers $0.5(E_{\textrm{Max}}-E_{\textrm{Min}})$ (MeV)
of total cross sections ($\Sigma_{t}$) from HDPE data, against the analytic
solution for scalar flux (Equation 19) divided by the energy bin width for 12,
17, and 618 group data sets.

(a)

(b)

(c)

Figure 4: Multigroup relative residuals defined as $||\phi_{n}-\phi_{n-1}||$
for 12 (a), 70 (b), and 618 (c) group problems with $N=2048$.

Figure 5:  Relative error from the 12 group QMC and MC scalar flux results
compared to analytic solution for varying number of particles $N$ and spatial
cells $N_{x}$. The QMC result achieve the expected $O(N^{-1})$ convergence,
while the MC results converge at the standard $O(N^{-1/2})$.

###  4.2 Problem 2: Reed’s Multi-Media Problem

The second problem, known as Reed’s Problem, is a mono-energetic multi-media
problem in slab geometry [37]. The problem features 5 unique media that are
reflected across the problem for a total of 9 regions, see Figure 6. To ensure
that each spatial cell contained only one media, Reed’s Problem was run so
that $N_{x}$ was evenly divisible by 16, the range of the problem.

Reed’s Problem was benchmarked using results from a $N=10^{10}$ analog Monte
Carlo simulation from CEMeNT’s Monte Carlo Dynamic Code (MCDC) [24].
Additionally, in analyzing the relative error, it was observed that the
solution to the problem approaches zero in multiple locations and this was
drastically increasing the relative error as reported in the previous problem.
Instead, for Reed’s Problem we report the $L_{\infty}$ norm of the error as
seen in Equation 21.

$R=||\phi^{\textrm{Sol}}-\phi^{QMC}||_{\infty}.$ (21)

Again, it is observed that the Krylov solvers far outperform Source Iteration,
requiring no less than one-quarter the number of iterations to converge
(Figure 7). Unlike problem 4.1 however, the effects of increasing the number
of spatial cells and utilizing the spatially averaged scalar flux algorithm
are clearly seen (Figure 8). Our QMC Sweep currently utilizes a flat source
for each cell, and beyond a certain number of particle histories $N$, the
convergence is limited by spatial error determined by $N_{x}$. As the number
of spatial cells is increased, the spatial error is reduced and the QMC can
continue to converge at the $O(N^{-1})$ rate. Note, this effect is not
observed in the MC results, because they did not reach the spatial error limit
for $N_{x}=80$, near $10^{-2}$. Finally, Figure 9 shows the results from a
simulation with $N_{x}=80$ using the Sobol, Halton, and Golden sequences and a
pseudo-random number generator (MC). The MC results perform as expected and
the three LDS perform rather similarly. The Golden sequence achieves a lower
error for the first few particle counts but the Sobol Sequence results plateau
at a lower error and ultimately achieves the greatest accuracy.

Figure 6:  MCDC solution for scalar flux $\phi$ with regional cross section
data and source strength.

Figure 7:  Scalar flux relative residuals from Reed’s Problem given $N=2048$
and $N_{x}=128$.

Figure 8:  $L_{\infty}$ of the absolute error of QMC scalar flux results
compared to MCDC for varying number of particles $N$ and spatial cells
$N_{x}$. The QMC results converge at $O(N^{-1})$ until limited by the spatial
error. Increasing $N_{x}$ lowers the spatial error limit.

Figure 9:  Comparison of the Sobol, Halton, and Golden Sequences along with
results from a pseudo-random number generator for $N_{x}=80$. The Golden
Sequence achieves the best accuracy for low particle histories but is
eventually out performed by the Sobol Sequence.

###  4.3 Problem 3: Isotropic Boundary Source

The third and final computational experiment solves the mono-energetic problem
from Garcia et al. [25], outlined in Table 1. Here, the scattering cross
section is spatially dependent, defined by $\Sigma_{s}=e^{-x/s}$, and we
considered two cases: $s=1$ and $s=\infty$. Note that $s=\infty$ is equivalent
to a constant scattering cross section and therefore is the harder of the two
scenarios as it involves more scattering and therefore more iterations for the
scattering source to converge.

First, we solve the QMC linear problem with $N=2048$ particles and $N_{x}=100$
spatial cells. Similar to Problems 4.2 and 4.1, Figure 10 shows that for an
exponentially decaying scattering cross section ($s=1$) the Krylov iterations
take fewer than a third of the number of transport sweeps than that of the SI
for a relative residual of $10^{-9}$. While Figure 11 shows that for the
constant scattering cross section ($s=\infty$) the Krylov iterations took less
than 25 iterations to reach a relative error of $10^{-6}$ while the SI
required nearly 200 iterations.

Table 1: Parameters for fixed boundary source, slab geometry, simulation from
Garcia et al. [25]

Parameter Value $\Sigma_{t}$ 1 $\Sigma_{s}(x)$ $e^{-x/s}$ $\tau$ $5$
$\psi_{l}(\mu)$ 1 $\psi_{r}(\mu)$ 0 $N_{x}$ 50 $q(x)$ 0

Figure 10:  Scalar flux relative residuals for $s=1$ given parameters from
Table 1 and $N=2048$ and $N_{x}=100$.

Figure 11:  Scalar flux relative residuals for $s=\infty$ given parameters
from Table 1 and $N=2048$ and $N_{x}=100$.

####  4.3.1 Validation and calibration study

We conclude this problem with a validation study given results from [25] that
are angular flux exit distributions accurate to six figures, see Figure 12. We
duplicated the results by obtaining the cell-average scalar flux from the QMC
simulation, for $N=2048$ and $N_{x}=100$. Then used a single $S_{N}$ transport
sweep to recover the exit distributions. We report the corresponding results
from [25] in Tables 2 and 3. The exit distributions, as is clear from Table 2,
can vary by five orders of magnitude. Even so, the results from QMC agree with
the benchmarks to roughly two figures.

Figure 12:  Angular flux exit distributions provided by Garcia et al. [25]. Recreated using scalar flux data from QMC and subsequent angular flux exit distributions from an $S_{N}$ approximation  Table 2: Angular flux exit Distributions from Garcia et al. and the $S_{N}$ QMC sweep (with $N=2048$ and $N_{x}=100$) for $s=1$. | Garcia/Siewert | QMC
---|---|---
$\mu$ | $\psi(0,-\mu)$ | $\psi(\tau,\mu)$ | $\psi(0,-\mu)$ | $\psi(\tau,\mu)$
0.05 | 5.89664e-01 | 6.07488e-06 | 6.07035e-01 | 5.91908e-06
0.10 | 5.31120e-01 | 6.92516e-06 | 5.47466e-01 | 6.74075e-06
0.20 | 4.43280e-01 | 9.64232e-06 | 4.57064e-01 | 9.35453e-06
0.30 | 3.80306e-01 | 1.62339e-05 | 3.92223e-01 | 1.56108e-05
0.40 | 3.32964e-01 | 4.38580e-05 | 3.43481e-01 | 4.13721e-05
0.50 | 2.96090e-01 | 1.69372e-04 | 3.05510e-01 | 1.58622e-04
0.60 | 2.66563e-01 | 5.73465e-04 | 2.75098e-01 | 5.39514e-04
0.70 | 2.42390e-01 | 1.51282e-03 | 2.50192e-01 | 1.43257e-03
0.80 | 2.22235e-01 | 3.24369e-03 | 2.29422e-01 | 3.08975e-03
0.90 | 2.05174e-01 | 5.96036e-03 | 2.11837e-01 | 5.70555e-03
1.00 | 1.90546e-01 | 9.77123e-03 | 1.96756e-01 | 9.39189e-03
Table 3: Angular flux exit Distributions from Garcia et al. and the $S_{N}$ QMC sweep (with $N=2048$ and $N_{x}=100$) for $s=\infty$. | Garcia/Siewert | QMC
---|---|---
$\mu$ | $\psi(0,-\mu)$ | $\psi(\tau,\mu)$ | $\psi(0,-\mu)$ | $\psi(\tau,\mu)$
0.05 | 8.97798e-01 | 1.02202e-01 | 9.06050e-01 | 1.03680e-01
0.10 | 8.87836e-01 | 1.12164e-01 | 8.95849e-01 | 1.13695e-01
0.20 | 8.69581e-01 | 1.30419e-01 | 8.76487e-01 | 1.31907e-01
0.30 | 8.52299e-01 | 1.47701e-01 | 8.58937e-01 | 1.49245e-01
0.40 | 8.35503e-01 | 1.64497e-01 | 8.42195e-01 | 1.66128e-01
0.50 | 8.18996e-01 | 1.81004e-01 | 8.25870e-01 | 1.82734e-01
0.60 | 8.02676e-01 | 1.97324e-01 | 8.09780e-01 | 1.99151e-01
0.70 | 7.86493e-01 | 2.13507e-01 | 7.93834e-01 | 2.15421e-01
0.80 | 7.70429e-01 | 2.29571e-01 | 7.77997e-01 | 2.31558e-01
0.90 | 7.54496e-01 | 2.45504e-01 | 7.62269e-01 | 2.47547e-01
1.00 | 7.38721e-01 | 2.61279e-01 | 7.46673e-01 | 2.63362e-01

In Tables 4 and 5 we look at the relative errors $R$ (Equations 22, 23, 24) in
the QMC exit distributions as compared to a highly accurate $S_{N}$ result. We
compensate for the widely varying scales by tabulating, for each value of $N$
and $N_{x}$. Finally, similar to the previous two problems we plot the
residual $R$ for varying $N$ and $N_{x}$ in Figure 13. However, unlike the
previous two problems the scalar flux averaging algorithm was not applied. In
the QMC case of $s=1.0$ we notice that if the number of spatial cells are
increased with the number of particles, the $O(N^{-1})$ convergence is
achieved. For $s=\infty$ the QMC results converge irrespective of the number
of spatial cells.

$R=\max(R^{0},R^{\tau})$ (22)

where

$R^{0}=\max_{\mu}\frac{|\psi^{SN}(0,-\mu)-\psi^{QMC}(0,-\mu)|}{\psi^{SN}(0,-\mu)},$
(23)

and

$R^{\tau}=\max_{\mu}\frac{|\psi^{SN}(\tau,\mu)-\psi^{QMC}(\tau,\mu)|}{\psi^{SN}(\tau,\mu)}.$
(24)

Table 4: Exit Distribution Errors ($R$): $s=1.0$, for varying the number of spatial cells ($N_{x}$) and particles per transport sweep ($N$). $N_{x}$ \N | 1024 | 2048 | 4096 | 8192 | 16384
---|---|---|---|---|---
50 | 1.36975e-01 | 1.34260e-01 | 1.35123e-01 | 1.35328e-01 | 1.35242e-01
100 | 6.09631e-02 | 6.35764e-02 | 6.46191e-02 | 6.48898e-02 | 6.48536e-02
200 | 3.77223e-02 | 3.12496e-02 | 3.12005e-02 | 3.17337e-02 | 3.16710e-02
400 | 1.71316e-02 | 1.52618e-02 | 1.96221e-02 | 7.56867e-03 | 7.71669e-03
800 | 9.58555e-03 | 1.01782e-02 | 1.93042e-02 | 7.56867e-03 | 7.71669e-03
1600 | 7.18716e-03 | 1.20870e-02 | 2.19512e-02 | 3.63519e-03 | 3.85206e-03
3200 | 5.11139e-03 | 1.28920e-02 | 2.19849e-02 | 2.36739e-03 | 1.82554e-03
Table 5: Exit Distribution Errors ($R$): $s=\infty$, for varying the number of spatial cells ($N_{x}$) and particles per transport sweep ($N$). $N_{x}$ \N | 1024 | 2048 | 4096 | 8192 | 16384
---|---|---|---|---|---
50 | 2.83725e-02 | 1.22374e-02 | 1.19906e-02 | 1.15398e-02 | 1.15951e-02
100 | 2.41784e-02 | 1.53219e-02 | 6.81242e-03 | 5.67372e-03 | 5.90660e-03
200 | 1.95960e-02 | 1.32236e-02 | 6.73505e-03 | 2.37596e-03 | 2.91740e-03
400 | 1.86114e-02 | 1.37280e-02 | 1.04038e-02 | 1.07665e-03 | 1.63724e-03
800 | 3.02170e-02 | 1.15453e-03 | 1.21401e-02 | 1.37999e-03 | 1.13099e-03
1600 | 2.26681e-02 | 3.93529e-02 | 1.78938e-02 | 1.17129e-03 | 1.17343e-03
3200 | 3.30886e-02 | 5.40682e-02 | 2.52993e-02 | 1.48919e-03 | 1.51562e-03

Figure 13:  Relative error of fixed-seed QMC and MC results for $s=1$ and
$s=\infty$ compared to $S_{N}$ solution for varying number of particles $N$
and spatial cells $N_{x}$. Increasing the spatial cells and particle count
simultaneously in QMC simulations achieves $O(N^{-1})$ convergence.

##  5 Conclusion

We have described a iQMC, a general purpose iterative - Quasi-Monte Carlo
method for solving the neutron transport equation. The use of iterative
solvers and a continuous weight absorption technique provide a well suited
application for QMC which provides an enhanced convergence rate of $O(N^{-1})$
compared to the $O(N^{-1/2})$ of standard Monte Carlo simulation.
Additionally, the use of advanced iterative solvers like the Krylov methods
GMRES and BiCGSTAB provide greatly increased convergence of residuals. The
benefits of this algorithm were observed on all three 1-D test problems where
the hybrid method provided both more accurate and efficient solutions. Future
work will look to see if these benefits are maintained on more difficult and
complex problems including critical eigenvalue, time-dependent, 2D, and 3D
problems. The massive parallelism inherent to Monte Carlo combined with iQMC’s
particle tracing and vectorized multigroup methods suggest the it would
benefit greatly from parallel implementation on advanced architectures and
GPUs.

## Acknowledgments

This work was funded by the Center for Exascale Monte-Carlo Neutron Transport
(CEMeNT) a PSAAP-III project funded by the Department of Energy, DE-NA003967,
and supported by National Science Foundation Grants DMS-1745654, and
DMS-1906446.

## References

  * [1] J. J. Duderstadt, L. J. Hamilton, S. Moorthy, and C. C. Scott, “Nuclear Reactor Analysis by James J. Duderstadt and Louis J. Hamilton,” _IEEE transactions on nuclear science_ , 24, 1983 (1977).

  * [2] R. L. Murray, “Particle-transport simulation with the Monte Carlo method,” _Annals of nuclear energy_ , 4, 473 (1977).
  * [3] E. Miller and W. J. Lewis, _Computational Methods of Neutron Transport_ , John Wiley and Sons (1984).
  * [4] M. L. Adams and E. W. Larsen, “Fast Iterative Methods for Discrete-Ordinates Particle Transport Calculations,” _Progress in Nuclear Energy_ , 40, 3 (2002).
  * [5] J. S. Warsa, T. A. Wareing, and J. E. Morel, “Krylov Iterative Methods and the Degraded Effectiveness of Diffusion Synthetic Acceleration for Multidimensional SN Calculations in Problems with Material Discontinuities,” _Nuclear Science and Engineering_ , 147, _3_ , 218 (2004).
  * [6] J. E. Morel, T. A. Wareing, R. B. Lowrie, and D. K. Parsons, “Analysis of ray-effect mitigation techniques,” _Nuclear science and engineering_ , 144, _1_ , 1 (2003).
  * [7] J. Willert, “Hybrid Deterministic/Monte Carlo Methods for Solving the Neutron Transport Equation and $k$-Eigenvalue Problem,” PhD Thesis, North Carolina State University, Raleigh, North Carolina (2013).
  * [8] M. M. R. Williams, _Random processes in nuclear reactors_ , Elsevier (2013).
  * [9] S. A. Dupree and S. K. Fraley, _A Monte Carlo primer: A Practical approach to radiation transport_ , vol. 1, Springer Science & Business Media (2002).
  * [10] M. H. Kalos and P. A. Whitlock, _Monte carlo methods_ , John Wiley & Sons (2009).
  * [11] R. McClarren, _Computational Nuclear Engineering and Radiological Science Using Python_ , Academic Press (2018).
  * [12] J. Willert, C. T. Kelley, D. A. Knoll, and H. K. Park, “Hybrid Deterministic/Monte Carlo Neutronics,” _SIAM J. Sci. Comp._ , 35, S62 (2013).
  * [13] H. Von Kleist, “Michael Kohlhaas,” _Erzählungen_ (1810).
  * [14] P. Bickel, P. Diggle, S. Feinberg, U. Gather, I. Olkin, and S. Zeger, _Monte Carlo and Quasi-Monte Carlo Sampling_ , Springer (2009)URL http://www.springer.com/series/692.
  * [15] J. Spanier, _Quasi-Monte Carlo Methods for Particle Transport Problems_ , Springer New York (1995).
  * [16] J. Farmer and S. Roy, “A quasi-Monte Carlo solver for thermal radiation in participating media,” _Journal of Quantitative Spectroscopy and Radiative Transfer_ , 242 (2020).
  * [17] W. J. MOROKOFF and R. E. CAFLISCH, “A Quasi-Monte Carlo Approach to Particle Simulation of the Heat Equation,” _SIAM journal on numerical analysis_ , 30, 1558 (1993).
  * [18] B. L. Fox, _Strategies for Quasi-Monte Carlo_ , vol. 22, Springer US (1999).
  * [19] P. H. de Almeida Konzen, L. F. Guidi, and T. Richter, “Quasi-random discrete ordinates method for neutron transport problems,” _Annals of nuclear energy_ , 133, 275 (2019); 10.1016/j.anucene.2019.05.017.
  * [20] Z. He and X. Wang, “On the convergence rate of randomized quasi–Monte Carlo for discontinuous functions,” _SIAM Journal on Numerical Analysis_ , 53, _5_ , 2488 (2015).
  * [21] L. Palluotto, N. Dumont, P. Rodrigues, O. Gicquel, and R. Vicquelin, “Assessment of randomized Quasi-Monte Carlo method efficiency in radiative heat transfer simulations,” _Journal of quantitative spectroscopy and radiative transfer_ , 236, 106570 (2019).
  * [22] S. Pasmann, I. Variansyah, and R. G. McClarren, “Convergent transport source iteration calculations with Quasi-Monte Carlo,” _Transactions of the American Nuclear Society_ , 124, 192 (2021).
  * [23] C. Mattoon, B. Beck, N. Patel, N. Summers, G. Hedstrom, and D. Brown, “Generalized nuclear data: A new structure (with supporting infrastructure) for handling nuclear data,” _Nuclear Data Sheets_ , 113, _12_ , 3145 (2012).
  * [24] I. Variansyah and R. G. McClarren, “Analysis of Population Control Techniques for Time-Dependent and Eigenvalue Monte Carlo Neutron Transport Calculations,” (2022); 10.48550/ARXIV.2202.08631., URL https://arxiv.org/abs/2202.08631.
  * [25] R. Garcia and C. Siewert, “Radiative transfer in finite inhomogeneous plane-parallel atmospheres,” _J. Quant. Spectrosc. Radiat. Transfer_ , 27, 141 (1982).
  * [26] C. T. Kelley, _Iterative Methods for Linear and Nonlinear Equations_ , no. 16 in Frontiers in Applied Mathematics, SIAM, Philadelphia (1995).
  * [27] Y. Saad and M. Schultz, “GMRES a generalized minimal residual algorithm for solving nonsymmetric linear systems,” _SIAM J. Sci. Stat. Comp._ , 7, 856 (1986).
  * [28] H. A. van der Vorst, “Bi-CGSTAB: A fast and smoothly converging variant to Bi-CG for the solution of nonsymmetric systems,” _sissc_ , 13, 631 (1992).
  * [29] M. Schauer and O. Schulz, “GoldenSequences.jl,” https://github.com/mschauer/GoldenSequences.jl (2021).
  * [30] S. G. Johnson, “Sobol.jl,” https://github.com/stevengj/Sobol.jl (2020).
  * [31] T. Driscoll and M. Piibeleht, “HaltonSequences.jl,” https://github.com/tobydriscoll/HaltonSequences.jl (2019).
  * [32] C. T. Kelley and S. Pasmann, “Krylov_QMC,” https://github.com/ctkelley/Krylov_QMC/ (2022)Julia Package.
  * [33] C. T. Kelley, “SIAMFANLEquations.jl,” https://github.com/ctkelley/SIAMFANLEquations.jl (2022); 10.5281/zenodo.4284807., URL https://github.com/ctkelley/SIAMFANLEquations.jl, julia Package.
  * [34] C. T. Kelley, “Notebook for Solving Nonlinear Equations with Iterative Methods: Solvers and Examples in Julia,” https://github.com/ctkelley/NotebookSIAMFANL (2022); 10.5281/zenodo.4284687., URL https://github.com/ctkelley/NotebookSIAMFANL, iJulia Notebook.
  * [35] C. T. Kelley, _Solving Nonlinear Equations with Iterative Methods: Solvers and Examples in Julia_ , no. 20 in Fundamentals of Algorithms, to appear in 2022, SIAM, Philadelphia (2022).
  * [36] A. B. Owen, “On dropping the first Sobol’point,” _International Conference on Monte Carlo and Quasi-Monte Carlo Methods in Scientific Computing_ , 71–86, Springer (2022).
  * [37] J. Warsa, “Analytical S N solutions in heterogeneous slabs using symbolic algebra computer programs,” _Annals of Nuclear Energy - ANN NUCL ENERG_ , 29, 851 (2002).