# Multigroup Neutron Transport using a Collision-Based Hybrid Method

###### Abstract

A collision-based hybrid algorithm for the discrete ordinates approximation of
the neutron transport equation is extended to the isotropic multigroup
setting. The algorithm uses discrete energy and angle grids at two different
resolutions and approximates the fission and scattering sources on the coarser
grids. The coupling of a collided transport equation, discretized on the
coarse grid, with an uncollided transport equation, discretized on the fine
grid, yields an algorithm that, in most cases, is more efficient than the
traditional multigroup approach. The improvement over existing techniques is
demonstrated for time-dependent problems with different materials, geometries,
and energy groups.

Ben Whewella Ryan G. McClarrena<EMAIL_ADDRESS>Cory D. Hauckb Minwoo Shina

aDepartment of Aerospace and Mechanical Engineering
University of Notre Dame, Notre Dame, IN 46556, USA bComputer Science and
Mathematics Division
Oak Ridge National Laboratory, Oak Ridge, TN 37830, USA

Neutron Transport Hybrid Methods Multigroup Approximation

##  1 Introduction

The neutron transport equation (NTE) is used to model neutron populations
traveling through different media. A common computational technique used to
solve this equation numerically is the discrete ordinates (SN) method
[lewis1993, chandrasekhar2013]. In a multigroup discrete ordinates
calculation, time-dependent problems are often discretized implicitly. This
results in a steady-state problem at each time step in the computation that is
solved iteratively. Typically there are two levels of iteration. The inner
iteration updates the scalar flux inside each energy group using a fixed-point
source iteration scheme, where the in-group scattering source is lagged and
each update requires the inversion of the operator that models streaming and
loss. These updates, called transport sweeps, amount to the inversion of a
block triangular system in space for each angle [o1987transport]. The outer
iterations over the energy variable often follow a Gauss-Seidel strategy that
updates the scattering and fission source terms as the energy groups are
updated in order from highest to lowest energies [bell1970nuclear]. The Gauss-
Seidel method is effective in problems that are dominated by down scattering;
however, it is dependent on the number of interactions that the neutrons have
with the material, such as scattering and, in time-dependent problems,
fission. For example, a purely down-scattering material without fission will
converge in one Gauss-Seidel iteration. If there is up-scattering or fission,
the number of iterations can become prohibitively large in optically thick
problems [lewis1993].

The cost of a multigroup calculation is a function of the energy resolution
needed; coarsening the energy grid, thus reducing the number of energy groups,
can lead to less expensive solutions. However, a smaller energy group
structure must be carefully chosen to limit the discretization error and
preserve the characteristics of the energy-dependent cross section [njoy].
Likewise, there can be faster convergence times using a low number of angles,
but this can also lead to both larger errors and exacerbate the non-physical
ray effects that can arise in SN calculations. There are techniques for
reducing these ray effects [hauck2019filtered, frank2020], but the most common
solution is to increase the number of angles [lewis1993], thus increasing the
computational time.

Acceleration techniques, such as coarse mesh rebalancing and diffusion
synthetic acceleration, are used to reduce the computational cost of the
iterative solver. Unfortunately, these methods cannot be implemented
indiscriminately, as coarse mesh rebalancing must be concerned with the coarse
mesh size [lewis1993] while diffusion synthetic acceleration has difficultly
in highly heterogeneous materials [southworth2021]. For the outer iterations,
there are two-gridand , nonlinear diffusion acceleration, and Krylov subspace
schemes that can improve the convergence of problems with upscattering or
fission [adams1993two, anistratov2013multilevel][adams1993two,
anistratov2013multilevel, slaybaugh2018eigenvalue]. While these methods
demonstrate improvements with high upscattering materials, there is a nominal
improvement over Gauss-Seidel in low upscattering materials [adams1993two].
These methods must also include the coarse grid diffusion equation solver for
their transport code bases, something that is not needed for the collision-
based hybrid method. Another approach would be to use a higher-order time
integration method [anistratov2013multilevel, edwards2011nonlinear,
mcclarren2008quasi, thoreson2009high]. Nevertheless, the robustness of the
fully-implicit backward Euler method makes it a standard choice for transport
problems [lewis1993, §2.5].

More recently, collision-based hybrid algorithms for time-dependent transport
equations have been provender for transport researchers [hauck2013,
crockatt2017arbitrary, crockatt2019hybrid, crockatt2020improvements,
heningburg2020]. These algorithms split the transport equation into collided
and uncollided components, as was done much earlier in the steady-state
context [alcouffe1977]. In the aforementioned examples, the hybrid approach
yielded substantial improvements in efficiency when compared to monolithic
discretization schemes. HoweverWhile previous work has applied different
spatial discretization schemes to the collided and uncollided components, the
approach has so far only been implemented on mono-energetic problems, and in .
In the current work, we extend the collision-based hybrid algorithm to the
multigroup setting. Specifically, we apply the algorithm to the multigroup
discrete ordinates approximation of the time-dependent neutron transport
equation in one-dimensional problem geometries. We test the algorithm on one-
dimensional slabs and spheres with different boundary conditions and sources.
We find that the hybrid algorithm can for the same computational time deliver
much smaller error or reach the same error level produce solutions in
significantly less time. Additionally, we find that the applying a collided-
uncollided splitting to a monolithic discretization grid recovers the same
solution as standard Gauss-Seidel in less time. This improvement in time-to-
solution is due to the changes in the solver that reduce unnecessary
iterations in the inner loop of the nested approach. A more focused study of
oversolving with nested iterations and how to avoid it can be found in
[senecal2017].

The remainder of this paper is organized as follows. In Section 2, we
introduce the neutron transport equation, the continuous collided-uncollided
split, the multigroup discrete ordinates transport equations, and the
collision-based hybrid algorithm. In Section 3, we present numerical results
for three different test problems and compare the hybrid method to more
traditional approaches. In Section 4, we present conclusions and discuss
future work.

##  2 The neutron transport equation and the hybrid formulation

The NTE models the behavior of advecting neutrons that interact with a
surrounding material medium. It does so by tracking the evolution of the
angular flux $\Psi$ which depends on a spatial coordinate $\mathbf{x}\in
D\subset\mathbb{R}^{3}$, angular coordinate $\bm{\Omega}\in\mathbb{S}^{2}$,
energy $E>0$, and time $t>0$. The NTE takes the form [lewis1993]

$\begin{split}\left(\frac{1}{c(E)}\frac{\partial}{\partial
t}+\Omega\cdot\nabla+\sigma^{\rm{t}}(\mathbf{x},E)\right)\Psi(\mathbf{x},\bm{\Omega},E,t)&=\int_{0}^{\infty}dE^{\prime}\;\sigma^{\rm{s}}(\mathbf{x},E^{\prime}\rightarrow
E,t)\overline{\Psi}(\mathbf{x},E^{\prime},t)\\\
&+\chi(\mathbf{x},E)\int_{0}^{\infty}dE^{\prime}\;\nu(x,E^{\prime})\sigma^{\rm{f}}(\mathbf{x},E^{\prime})\overline{\Psi}(\mathbf{x},E^{\prime},t)+Q(\mathbf{x},\bm{\Omega},E,t),\end{split}$
(1)

where the scalar flux
$\overline{\Psi}:\left(D\times\mathbb{R}^{+}\times\mathbb{R}^{+}\right)\to\mathbb{R}$
is given by

$\overline{\Psi}(\mathbf{x},E,t)=\frac{1}{4\pi}\int_{\mathbb{S}^{2}}\Psi(\mathbf{x},\bm{\Omega},E,t)d\Omega,$
(2)

and $c(E)$ is the velocity of a neutron with energy $E$. The three terms on
the left hand side of Eq. (1) are, in order, the time derivative term,
streaming term, and removal term; the right hand side is comprised of the
scattering, fission, and external sources. The material is characterised by
cross-sections $\sigma^{\rm{s}}$, $\sigma^{\rm{f}}$, and $\sigma^{\rm{t}}$.
More precisely, $\sigma^{\rm{s}}(\mathbf{x},E^{\prime}\rightarrow E)$ is the
differential scattering cross section for neutrons scattering from energy
$E^{\prime}$ to $E$; $\sigma^{\rm{f}}(\mathbf{x},E^{\prime})$ is the fission
cross section at energy $E^{\prime}$; and $\sigma^{\rm{t}}(\mathbf{x},E)$ is
the total cross section at energy $E$. Given a fission event,
$\nu(\mathbf{x},E^{\prime})$ is the average number of neutrons created from
fission caused by a neutron at energy $E^{\prime}$ while $\chi(\mathbf{x},E)$
is the probability density that a neutron produced from fission will have
energy $E$.

The initial condition for Eq. (1) is

$\Psi(\mathbf{x},\bm{\Omega},E,t=0)=f(\mathbf{x},\bm{\Omega},E)\qquad\text{for}\quad\mathbf{x}\in
D,\quad\bm{\Omega}\in\mathbb{S}^{2},\quad E>0,$ (3)

where $f$ is given. The incoming boundary data is prescribed as

$\Psi(\mathbf{x},\bm{\Omega},E,t)=b(\mathbf{x},\bm{\Omega},E,t)\qquad\text{for}\quad\mathbf{x}\in\partial
D,\quad\mathbf{n}(\mathbf{x})\cdot\bm{\Omega}<0,\quad E>0,\quad t>0,$ (4)

where $\mathbf{n}(\mathbf{x})$ the unit outward normal at
$\mathbf{x}\in\partial D$, the boundary of $D$.

###  2.1 The continuous collided-uncollided split

The goal of the hybrid approach is to accelerate the computation of a
numerical solution for Eq. (1) relative to a monolithic discretization, while
minimizing the subsequent loss in accuracy. The formulation of the hybrid
method is best understood at the continuous level as a splitting method.
Indeed, because it is linear, Eq. (1) can be separated into two equations
where a distinct angular flux is calculated in each and the results added
together to get the total angular flux.

The first equation governs the uncollided angular flux, $\Psi^{\rm{u}}$. It
takes the form

$\begin{split}\left(\frac{1}{c(E)}\frac{\partial}{\partial
t}+\Omega\cdot\nabla+\sigma^{\rm{t}}(\mathbf{x},E)\right)\Psi^{\rm{u}}(\mathbf{x},\bm{\Omega},E,t)&=Q^{\rm{u}}(\mathbf{x},\bm{\Omega},E,t),\end{split}$
(5)

with $Q^{\rm{u}}=Q$. The initial condition for $\Psi^{\rm{u}}$ is

$\Psi^{\rm{u}}(\mathbf{x},\bm{\Omega},E,0)=f(\mathbf{x},\bm{\Omega},E)\qquad\text{for}\quad\mathbf{x}\in
D,\quad\bm{\Omega}\in\mathbb{S}^{2},\quad E>0,$ (6)

and the incoming boundary data is

$\Psi^{\rm{u}}(\mathbf{x},\bm{\Omega},E,t)=b(\mathbf{x},\bm{\Omega},E,t)\qquad\text{for}\quad\mathbf{x}\in\partial
D,\quad\mathbf{n}(\mathbf{x})\cdot\bm{\Omega}<0,\quad E>0,\quad t>0.$ (7)

The second equation governs the collided angular flux, $\Psi^{\rm{c}}$. It
takes the form

$\begin{split}\left(\frac{1}{c(E)}\frac{\partial}{\partial
t}+\Omega\cdot\nabla+\sigma^{\rm{t}}(\mathbf{x},E)\right)\Psi^{\rm{c}}(\mathbf{x},\bm{\Omega},E,t)&=\int_{0}^{\infty}dE^{\prime}\sigma^{\rm{s}}(\mathbf{x},E^{\prime}\rightarrow
E,t)\overline{\Psi^{\rm{c}}}(\mathbf{x},E^{\prime},t)\\\
&+\chi(\mathbf{x},E)\int_{0}^{\infty}dE^{\prime}\nu(\mathbf{x},E^{\prime})\sigma^{\rm{f}}(\mathbf{x},E^{\prime})\overline{\Psi^{\rm{c}}}(\mathbf{x},E^{\prime},t)+Q^{\rm{c}}(\mathbf{x},E,t),\end{split}$
(8)

where the isotropic term $Q^{\rm{c}}$ comes from the scattering and fission
source of the uncollided flux:

$Q^{\rm{c}}(\mathbf{x},E,t)=\int_{0}^{\infty}dE^{\prime}\;\sigma^{\rm{s}}(\mathbf{x},E^{\prime}\rightarrow
E,t)\overline{\Psi^{\rm{u}}}(\mathbf{x},E^{\prime},t)+\chi(\mathbf{x},E)\int_{0}^{\infty}dE^{\prime}\;\nu(\mathbf{x},E^{\prime})\sigma^{\rm{f}}(\mathbf{x},E^{\prime})\overline{\Psi^{\rm{u}}}(\mathbf{x},E^{\prime},t).$
(9)

The initial condition for $\Psi^{\rm{c}}$ is

$\Psi^{\rm{c}}(\mathbf{x},\bm{\Omega},E,0)=0\qquad\text{for}\quad\mathbf{x}\in
D,\quad\bm{\Omega}\in\mathbb{S}^{2},\quad E>0,$ (10)

and the incoming boundary data is

$\Psi^{\rm{c}}(\mathbf{x},\bm{\Omega},E,t)=0\qquad\text{for}\quad\mathbf{x}\in\partial
D,\quad\mathbf{n}(\mathbf{x})\cdot\bm{\Omega}<0,\quad E>0,\quad t>0.$ (11)

A third equation for the total angular flux $\Psi^{\rm{t}}$ takes the form

$\begin{split}\left(\frac{1}{c(E)}\frac{\partial}{\partial
t}+\Omega\cdot\nabla+\sigma^{\rm{t}}\right)\Psi^{\rm{t}}(\mathbf{x},\bm{\Omega},E,t)&=Q^{\rm{t}}(\mathbf{x},\bm{\Omega},E,t),\end{split}$
(12)

where the total external source $Q^{\rm{t}}$ is

$\begin{split}Q^{\rm{t}}(\mathbf{x},\bm{\Omega},E,t)&=Q^{\rm{c}}(\mathbf{x},\bm{\Omega},E,t)\\\
&+\int_{0}^{\infty}dE^{\prime}\;\sigma^{\rm{s}}(\mathbf{x},E^{\prime}\rightarrow
E,t)\overline{\Psi^{\rm{c}}}(\mathbf{x},E^{\prime},t)+\chi(\mathbf{x},E)\int_{0}^{\infty}dE^{\prime}\;\nu\sigma^{\rm{f}}(\mathbf{x},E^{\prime})\overline{\Psi^{\rm{c}}}(\mathbf{x},E^{\prime},t).\end{split}$
(13)

The initial condition for $\Psi^{\rm{t}}$ is

$\Psi^{\rm{t}}(\mathbf{x},\bm{\Omega},E,0)=f(\mathbf{x},\bm{\Omega},E),\qquad\text{for}\quad\mathbf{x}\in
D,\quad\bm{\Omega}\in\mathbb{S}^{2},\quad E>0,$ (14)

and the incoming boundary data is

$\Psi^{\rm{t}}(\mathbf{x},\bm{\Omega},E,t)=b(\mathbf{x},\bm{\Omega},E,t)\qquad\text{for}\quad\mathbf{x}\in\partial
D,\quad\mathbf{n}(\mathbf{x})\cdot\bm{\Omega}<0,\quad E>0,\quad t>0.$ (15)

In the continuum formulation above, uniqueness of solutions to the neutron
transport equations implies that
$\Psi^{\rm{t}}=\Psi^{\rm{u}}+\Psi^{\rm{c}}=\Psi$, making an independent
equation for $\Psi^{\rm{t}}$ trivially redundant. However, this is no longer
the case when equations for $\Psi^{\rm{t}}$, $\Psi^{\rm{u}}$, and
$\Psi^{\rm{c}}$ are discretized with different methods and/or grids. In
practice, the hybrid method seeks at each time step a numerical approximation
for $\Psi^{\rm{u}}$ on a fine grid and a numerical approximation for
$\Psi^{\rm{c}}$ on a course grid. In the original formulation [hauck2013],
these two approximate solutions were, at the end of the time step, combined
via a remapping procedure in which an approximation for $\Psi^{\rm{c}}$ is
reconstructed on the fine grid. Unfortunately, the reconstruction can
introduce artifacts, particularly for $S_{N}$ discretizations in multi-
dimensional settings. To address this issue, the idea of using the equation
for $\Psi^{\rm{t}}$ to remap onto the fine grid was introduced in
[crockatt2020improvements] since the source $Q^{\rm{t}}$ depends only on
integrated quantities that can be computed on both grids. An alternative view
is that the collided/uncollided split provides a cheap method for
approximating the scattering and, in this work, fission sources that appear in
the right-hand side of Eq. (1).

The choices of initial conditions, boundary conditions, and sources for
$\Psi^{\rm{u}}$ and $\Psi^{\rm{c}}$ are not unique. The strategy above is to
assign data to the uncollided equation since it is equipped with the finest
discretization. However this strategy may not always be the best choice,
particularly in highly-collisional problems with strong boundary layers
[Densmore2006cla].

###  2.2 The multigroup, discrete ordinates equations

As an illustrative demonstration of the hybrid approach, we consider in this
paper discretizations based on multigroup, discrete ordinates (SN) equations
[chandrasekhar2013]. We begin with a discussion of the monolithic method. For
the angular discretization, let $\bm{\Omega}_{m}$ and $w_{m}$, where
$m\in\mathcal{M}:=\\{1,\cdots,M\\}$, be discrete angles and weights for a
quadrature rule over the sphere: that is, for any integrable function $u$
defined point-wise everywhere on $\mathbb{S}^{2}$,

$\frac{1}{4\pi}\int_{\mathbb{S}^{2}}d\bm{\Omega}\,u(\bm{\Omega})\approx\sum_{m=1}^{M}w_{m}u(\bm{\Omega}_{m}).$
(16)

For the energy discretization, $E_{\rm{max}}$ is a practical, finite upper
bound on the neutron energy spectrum and let
$0=E_{0}<E_{1}<\dots<E_{G}=E_{\rm{max}}$ be a set of finite energies that form
the endpoints of $G$ non-overlapping intervals (or groups) of width $\Delta
E_{g}=E_{g}-E_{g-1}$. Then define the function

$\psi_{m,g}(x,t)\approx\int_{E_{g-1}}^{E_{g}}dE\,\Psi(\mathbf{x},\bm{\Omega}_{m},E,t),\qquad
g\in\mathcal{G}:=\\{1,\dots,G\\}$ (17)

as the solution of the multigroup, discrete ordinate (SN) equations:

$\frac{1}{c_{g}}\frac{\partial}{\partial
t}\psi_{m,g}+\mathbf{{\Omega}}_{m}\cdot\nabla\psi_{m,g}+\sigma^{\rm
t}_{g}\psi_{m,g}=\sum_{g^{\prime}=1}^{G}\sigma^{\rm s}_{g^{\prime}\rightarrow
g}\overline{\psi}_{g^{\prime}}+\chi_{g}\sum_{g^{\prime}=1}^{G}\nu_{g^{\prime}}\sigma^{\rm
f}_{g^{\prime}}\overline{\psi}_{g^{\prime}}+{q}_{m,g},$ (18)

where

$\overline{\psi}_{g}=\sum_{m=1}^{M}w_{m}\psi_{m,g}$ (19)

and the notation $g^{\prime}\to g$ denotes scattering from energy group
$E_{g^{\prime}}$ to $E_{g}$. The quantities $\sigma^{\rm t}_{g}$, $\sigma^{\rm
s}_{g^{\prime}\rightarrow g}$, $\chi_{g}$ and $\nu_{g^{\prime}}$, and
$\sigma^{\rm f}_{g^{\prime}}$ are all approximate weighted averages of their
continuum counterparts. For example,

$\displaystyle\sigma^{\rm s}_{g^{\prime}\rightarrow
g}(x,t)\approx\frac{\displaystyle\int_{E_{g-1}}^{E_{g}}\displaystyle\int_{E_{g^{\prime}-1}}^{E_{g^{\prime}}}dEdE^{\prime}\,\sigma^{\rm{s}}(\mathbf{x},E^{\prime}\rightarrow
E,t)\overline{\Psi}(\mathbf{x},E^{\prime},t)}{\displaystyle\int_{E_{g^{\prime}-1}}^{E_{g^{\prime}}}dE^{\prime}\,\overline{\Psi}(\mathbf{x},E^{\prime},t)}\quad\text{and}\quad\chi_{g}(x)\approx\frac{\displaystyle\int_{E_{g-1}}^{E_{g}}dE\,\chi(\mathbf{x},E)\overline{\Psi}(\mathbf{x},E,t)}{\displaystyle\int_{E_{g-1}}^{E_{g}}dE\,\overline{\Psi}(\mathbf{x},E,t)}.$
(20)

The approximation comes from the fact that $\overline{\Psi}(\mathbf{x},E,t)$
is not known a priori and an assumed spectral (and angular) shape of the
solution must be used. In practice, these quantities are pre-calculated by
nuclear data processing software such as NJOY [njoy] or Fudge [fudge] and,
henceforth, are assumed to be given. In our work the neutron velocity
$c_{g}^{-1}$ is calculated from the relativistic energy formula from
[bertozzi1964speed, Figure 3] at the midpoint of the energy bin.

The initial condition for $\psi_{m,g}$ is

$\psi_{m,g}(\mathbf{x},0)=f_{m,g}(\mathbf{x})\qquad\text{for}\quad\mathbf{x}\in
D,\quad m\in\mathcal{M},\quad g\in\mathcal{G}.$ (21)

and the incoming boundary data is

$\psi_{m,g}(\mathbf{x},t)=b_{m,g}(\mathbf{x},t)\qquad\text{for}\quad\mathbf{n}(\mathbf{x})\cdot\bm{\Omega}_{m}<0,\quad
m\in\mathcal{M},\quad g\in\mathcal{G},\quad t>0,$ (22)

where

$f_{m,g}(\mathbf{x})=\int_{E_{g-1}}^{E_{g}}dE\,f(\mathbf{x},\bm{\Omega}_{m},E)\qquad\text{and}\qquad
b_{m,g}{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\uwave{=}}(\mathbf{x},t)\int_{E_{g-1}}^{E_{g}}dE\,b(\mathbf{x},\bm{\Omega}_{m},E,t).$
(23)

Equation (18) is a set of $M\cdot G$ PDEs that must be further discretized in
$x$ and $t$. In this work, we use the diamond difference discretization in $x$
(see e.g., [lewis1993, Section 3.3]) and backward Euler time integration. For
one-dimensional spherical geometry, angular differencing coefficients, as
described in [lewis1993, Section 3.4], are used in addition to the diamond
difference discretization. The use of implicit time integration is standard
because of the time scales involved [larsen2010advances]; various higher-order
temporal discretizations are also possible [crockatt2017arbitrary,
crockatt2019hybrid] but not critical for our presentation. Since diamond
differencing is fairly standard, we omit the details here.

The fully discrete approximation of the angular flux $\psi_{m,g}$ is computed
using a nested iteration procedure that includes an outer loop for the energy
groups and an inner iteration for the angles. The procedure is stated in
Algorithm 1 of Appendix A. A standard sweeping method [bell1970nuclear] forms
the core of the inner iteration. In this process, the angular flux is
calculated at the cell boundary and used to average the flux at the cell
center through the diamond difference method. The inner iteration will sweep
in all angular directions using a fixed-point source iteration scheme until
convergence of a scalar flux $\overline{\psi}_{g}$ at a specific energy group.
The convergence of a one group flux is used to update the outer iteration,
which uses the Gauss-Seidel method over each energy group. This method will
converge quickly when there is minimal upscattering [lewis1993].

###  2.3 The hybrid multigroup, discrete ordinates equations

In the hybrid approach, the multigroup, discrete ordinates approximation is
applied to the uncollided, collided, and total flux equations. However, we now
allow for the possibility that different levels of resolution are used in each
case. We will assume here that the uncollided and total flux equations are
discretized at the same high resolution using $G$ groups and $M$ angles, while
the collided equation uses lower resolution $\hat{G}$ and $\hat{M}$ angles. To
differentiate between different discretization parameters, a hat adornment
$(\hat{\cdot})$ is used for the collided parameters. Thus the weights and
angles in the discretization of the collided equations are denoted by
${\hat{w}}_{{\hat{m}}}$ and ${\hat{\bm{\Omega}}}_{{\hat{m}}}$, respectively,
for ${{\hat{m}}}=1,\dots,{\hat{M}}$. To define the coarse groups, let

$0=\Gamma_{0}<\Gamma_{1}\dots<\Gamma_{\hat{g}}\dots<\Gamma_{\hat{G}}=G$ (24)

be a set of $\hat{G}+1$ integers and set
${\hat{E}}_{{\hat{g}}}=E_{\Gamma_{\hat{g}}}$. Then for each
$\hat{g}\in{\hat{\mathcal{G}}}$,

$\Delta{\hat{E}}_{{\hat{g}}}={\hat{E}}_{{\hat{g}}}-{\hat{E}}_{{\hat{g}}-1}=E_{\Gamma_{\hat{g}}}-E_{\Gamma_{\hat{g}-1}}=\sum_{g=\Gamma_{{\hat{g}}-1}+1}^{\Gamma_{{\hat{g}}}}E_{g}-E_{g-1}=\sum_{g=\Gamma_{{\hat{g}}-1}+1}^{\Gamma_{{\hat{g}}}}\Delta
E_{g}{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\uwave{.}}$
(25)

The uncollided equation for

$\psi^{\rm{u}}_{m,g}\approx\int_{E_{g-1}}^{E_{g}}dE\,\Psi^{\rm{u}}(\mathbf{x},\bm{\Omega}_{m},E,t)$
(26)

is

$\frac{1}{c_{g}}\frac{\partial}{\partial
t}\psi^{\rm{u}}_{m,g}+\bm{\Omega}_{m}\cdot\nabla\psi^{\rm{u}}_{m,g}+\sigma^{\rm
t}_{g}\psi^{\rm{u}}_{m,g}=q^{\rm{u}}_{m,g},$ (27)

where $q^{\rm{u}}_{m,g}=q_{m,g}$. We solve Eq. (27) over a time-step
$[t^{n},t^{n+1})$ with initial condition

$\psi^{\rm{u}}_{m,g}(\mathbf{x},t^{n})=\begin{cases}f_{m,g}(\mathbf{x}),&\mathbf{x}\in
D,\quad m\in\mathcal{M},\quad g\in\mathcal{G},\quad t^{n}=0\\\
\psi^{\rm{t}}_{m,g}(\mathbf{x},t^{n}_{-}),&\mathbf{x}\in D,\quad
m\in\mathcal{M},\quad g\in\mathcal{G},\quad t^{n}>0\end{cases},$ (28)

and incoming boundary data

$\psi^{\rm{u}}_{m,g}=b_{m,g}(\mathbf{x},t)\qquad\text{for}\quad\mathbf{n}(\mathbf{x})\cdot\bm{\Omega}_{m}<0,\quad
m\in\mathcal{M},\quad g\in\mathcal{G}\quad\text{and}\quad t>0.$ (29)

The collided transport equation for

$\psi^{\rm{c}}_{\hat{m},\hat{g}}\approx\int_{E_{\hat{g}-1}}^{E_{\hat{g}}}dE\,\Psi^{\rm{c}}(\mathbf{x},\hat{\bm{\Omega}}_{m},E,t)$
(30)

is

$\frac{1}{c_{{\hat{g}}}}\frac{\partial}{\partial
t}\psi^{\rm{c}}_{{\hat{m}},{\hat{g}}}+\hat{\bm{\Omega}}_{\hat{m}}\cdot\nabla\psi^{\rm{c}}_{{\hat{m}},{\hat{g}}}+\hat{\sigma}^{\rm
t}_{g}\psi^{\rm{c}}_{{\hat{m}},{\hat{g}}}=\sum_{g^{\prime}=1}^{{\hat{G}}}\hat{\sigma}^{\rm
s}_{{\hat{g}}^{\prime}\rightarrow{\hat{g}}}\overline{\psi^{\rm{c}}}_{{\hat{g}}^{\prime}}+\hat{\chi}_{{\hat{g}}}\sum_{{\hat{g}}^{\prime}=1}^{{\hat{G}}}\hat{\nu}_{{\hat{g}}^{\prime}}\hat{\sigma}^{\rm
f}_{{\hat{g}}^{\prime}\rightarrow{\hat{g}}}\overline{\psi^{\rm{c}}}_{{\hat{g}}^{\prime}}+q^{\rm{c}}_{{\hat{g}}},$
(31)

where

$c_{{\hat{g}}}=\frac{1}{\Delta\hat{E}_{{\hat{g}}}}\sum_{g=\Gamma_{{\hat{g}}-1}}^{\Gamma_{{\hat{g}}}}c_{g},\quad\overline{\psi^{\rm{c}}}_{{\hat{g}}^{\prime}}=\sum_{{\hat{m}}=1}^{{\hat{M}}}{\hat{w}}_{{\hat{m}}}\psi^{\rm{c}}_{{\hat{m}},{\hat{g}}^{\prime}},\quad\hat{\chi}_{{\hat{g}}}=\sum_{g=\Gamma_{{\hat{g}}}+1}^{\Gamma_{{\hat{g}}+1}}\chi_{g}\quad\text{and}\quad\hat{\nu}_{{\hat{g}}^{\prime}}=\sum_{g=\Gamma_{{\hat{g}}}+1}^{\Gamma_{{\hat{g}}+1}}\nu_{g};$
(32)

the energy-coarsened cross-sections are given by

$\hat{\sigma}^{\rm
t}_{{\hat{g}}}=\frac{1}{\Delta\hat{E}_{{\hat{g}}}}\sum_{g=\Gamma_{{\hat{g}}-1}+1}^{\Gamma_{{\hat{g}}}}\Delta
E_{g}\sigma^{\rm
t}_{g}\quad\text{and}\quad\hat{\sigma}^{\rm\ell}_{{\hat{g}}}=\frac{1}{\Delta\hat{E}_{{\hat{g}}}}\sum_{g=\Gamma_{{\hat{g}}-1}+1}^{\Gamma_{{\hat{g}}}}\sum_{g^{\prime}=\Gamma_{{\hat{g}}-1}+1}^{\Gamma_{{\hat{g}}}}\Delta
E_{g^{\prime}}\sigma^{\rm\ell}_{g^{\prime}\rightarrow
g},\quad\ell\in\rm\\{s,f\\};$ (33)

and the isotropic source is

$q^{\rm{c}}_{{\hat{g}}}=\sum_{g=\Gamma_{{\hat{g}}-1}+1}^{\Gamma_{{\hat{g}}}}\sum_{g^{\prime}=1}^{G}\sigma^{\rm
s}_{g^{\prime}\rightarrow
g}\overline{\psi^{\rm{u}}}_{g^{\prime}}+\sum_{g=\Gamma_{{\hat{g}}-1}+1}^{\Gamma_{{\hat{g}}}}\chi_{g}\sum_{g^{\prime}=1}^{G}\nu_{g^{\prime}}\sigma^{\rm
f}_{g^{\prime}\rightarrow g}\overline{\psi^{\rm{u}}}_{g^{\prime}},$ (34)

with

$\overline{\psi^{\rm{u}}}_{g}=\sum_{m=1}^{M}w_{m}\psi^{\rm{u}}_{m,g}\quad\text{and}\quad\overline{\psi^{\rm{c}}}_{\hat{g}}=\sum_{\hat{m}=1}^{\hat{M}}{\hat{w}}_{{\hat{m}}}\psi^{\rm{c}}_{\hat{m},\hat{g}}.$
(35)

We solve Eq. (31) over a time-step $[t^{n},t^{n+1})$ with initial condition

$\psi^{\rm{c}}_{{\hat{m}},{\hat{g}}}(\mathbf{x},t^{n})=0\qquad\text{for}\quad\mathbf{x}\in
D,\quad{\hat{m}}\in{\hat{\mathcal{M}}}:=\\{1,\cdots,{\hat{M}}\\},\quad{\hat{g}}\in{\hat{\mathcal{G}}},\quad
t^{n}\geq 0$ (36)

and boundary data

$\psi^{\rm{c}}_{{\hat{m}},{\hat{g}}}=0\qquad\text{for}\quad\mathbf{n}(\mathbf{x})\cdot{\hat{\bm{\Omega}}}_{{\hat{m}}}<0,\quad{\hat{m}}\in{\hat{\mathcal{M}}},\quad{\hat{g}}\in{\hat{\mathcal{G}}},\quad
t>0.$ (37)

The total flux equation for

$\psi^{\rm{t}}_{m,g}(x,t)\approx\int_{E_{g-1}}^{E_{g}}dE\,\Psi^{\rm{t}}(\mathbf{x},\bm{\Omega}_{m},E,t)$
(38)

is

$\frac{1}{c_{g}}\frac{\partial}{\partial
t}\psi^{\rm{t}}_{m,g}+\bm{\Omega}_{m}\cdot\nabla\psi^{\rm{t}}_{m,\hat{g}}+\sigma^{\rm
t}_{g}\psi^{\rm{t}}_{m,g}=q^{\rm{t}}_{m,g}$ (39)

where

$q^{\rm{t}}_{m,g}=\frac{\Delta
E_{g}}{\Delta{\hat{E}}_{{\hat{g}}}}\left[q^{\rm{c}}_{{\hat{g}}}+\sum_{{\hat{g}}^{\prime}=1}^{{\hat{G}}}\sigma^{\rm
s}_{{\hat{g}}^{\prime}\rightarrow{\hat{g}}}\overline{\psi^{\rm{c}}}_{{\hat{g}}^{\prime}}+\chi_{{\hat{g}}}\sum_{{\hat{g}}^{\prime}=1}^{{\hat{G}}}\nu_{{\hat{g}}^{\prime}}\sigma^{\rm
f}_{{\hat{g}}^{\prime}\rightarrow{\hat{g}}}\overline{\psi^{\rm{c}}}_{{\hat{g}}^{\prime}}\right]{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\sout{.}}$
(40)

and ${\hat{g}}$ is the unique integer such that $\Gamma_{{\hat{g}}-1}+1\leq
g\leq\Gamma_{{\hat{g}}}$ or, equivalently,
$(E_{g-1},E_{g})\subset({\hat{E}}_{{\hat{g}}-1},{\hat{E}}_{{\hat{g}}})$. We
solve Eq. (39) over a time-step $[t^{n},t^{n+1})$ with initial data

$\psi^{\rm{t}}_{m,g}(\mathbf{x},t^{n})=\begin{cases}f_{m,g}(\mathbf{x}),&\mathbf{x}\in
D,\quad m\in\mathcal{M},\quad g\in\mathcal{G},\quad t^{n}=0\\\
\psi^{\rm{t}}_{m,g}(\mathbf{x},t^{n}_{-}),&\mathbf{x}\in D,\quad
m\in\mathcal{M},\quad g\in\mathcal{G},\quad t^{n}>0\end{cases},$ (41)

and incoming boundary condition

$\psi^{\rm{t}}_{m,g}=b_{m,g}(\mathbf{x},t)\qquad\text{for}\quad\mathbf{n}(\mathbf{x})\cdot\bm{\Omega}_{m}<0,\quad
m\in\mathcal{M},\quad g\in\mathcal{G}\quad\text{and}\quad t>0.$ (42)

In summary, Eqs. (27), (31), and (39) are solved in succession for each time
step. The solution of (39) at the end of the time step provides the initial
condition needed by Eqs. (27) and (39) at the next time step, providing the
mortice and tenon to consistently join the uncollided and collided
calculations. The initial condition for (31) is set to zero for each time
step.

As with the monolithic equations (18), we use backward Euler and diamond
differencing for the temporal and spatial discretizations, respectively, of
Eqs. (27), (31), and (39). Only one Euler step is used to advance the
solutions from $t^{n}$ to $t^{n+1}$, and we use the same time steps and
spatial mesh for the hybrid and monolithic equations. In cases with time-
dependent cross sections, the appropriate measures would be taken. For the
purpose of this paper, the cross sections do not change in time.

The solution procedure for solving the collided flux is similar to the
approach used in the monolithic case. It is described in Algorithm 2 of
Appendix A. It relies spatial spatially on the very same sweeps and fixed-
point source iteration for solving the flux at each energy group, while Gauss-
Seidel is used to integrate over groups [bell1970nuclear, lewis1993]. The
uncollided and the total flux updates do not require any iterations because
the right-hand side of the respective equation is fixed.

##  3 Numerical results

In this section, we demonstrate the computational benefits of the hybrid using
several time-dependent problems with different one-dimensional geometries.
While these examples use isotropic scattering, we expect to see similar
benefits with anisotropic scattering despite the coupling terms being
modified. This is an area for future research. It should also be noted that
the cross sections do not change in time for these results.

In the slab case, the reduction of the equations above is straightforward. In
the spherical case, the coordinate system introduces angular derivatives that
must be discretized. We do so using the standard approach given in [lewis1993,
Section 3.4]. Our reference solution uses the original number of energy groups
as constructed by the Fudge [fudge] software package for nuclear data
management, but the groups are coarsened according to the formulas in (33).

Three different solution strategies are used for each test problem:

  * •

Multigroup. This is the traditional neutron multigroup method with ${\hat{G}}$
groups and ${\hat{M}}$ angles. It uses Gauss-Seidel with a tolerance of
$\varepsilon_{G}=1\times 10^{-12}$ for the outer iteration and source
iteration with a tolerance of $\varepsilon_{M}=1\times 10^{-8}$ for the inner
iteration, as described in Algorithm 1. The coarsening strategy in Eqs. (32)
and (33) is used for lower fidelity multigroup models. The the accuracy of the
multigroup method is not necessarily monotonic in the number of groups. In
particular, it can happen that smaller number of groups yield better answers
due to the nonlinearity of the procedure, i.e., if the solution is separable
in energy, space, and angle, a single group calculation can be exact
[bell1970nuclear, §4.3]. Additionally, there can be cancellation of errors in
integrated quantities. This behavior is observed in some of the test cases
below.

  * •

Hybrid. This is the hybrid method described in Section 2.3. It uses $G$ groups
and $M$ angles in the discretization of the total and uncollided angular
fluxes, and it uses ${\hat{G}}$ groups and ${\hat{M}}$ angles in the
discretization of the collided angular flux. It uses the algorithm described
in Algorithm 2 with iteration tolerances of $\varepsilon_{G}=1\times 10^{-12}$
and $\varepsilon_{M}=1\times 10^{-8}$ for the outer and inner iterations,
respectively.

  * •

Splitting. This method is the same as the hybrid method described in Algorithm
2 with the same iteration tolerances ($\varepsilon_{G}=1\times 10^{-12}$ and
$\varepsilon_{M}=1\times 10^{-8}$), but it uses ${\hat{G}}$ groups and
${\hat{M}}$ angles for the collided, uncollided, and total flux. It is
included to investigate the effects of the hybrid discretization vs. the
solver strategy: on one hand, it should produce the same answers (up to
iteration tolerances) as the multigroup method but may show different
convergence behavior because the solver strategy is different. On the other
hand, it should be less expensive, but also less accurate than the hybrid that
uses finer resolution for the uncollided and total fluxes.

(a)

(b)

Figure 1: General setup of the (a) enriched uranium slab and the (b) enriched
and depleted uranium sphere. The enriched uranium in both of these problems
uses 20% enriched uranium-235. The depleted uranium used only uranium-238. The
boundary source performs in the same manner for both of these problems, which
is initially at unity and decays by a half every 0.1 $\mu$s after the initial
0.2 $\mu$s, with the elapsed time set at 1 $\mu$s, as described in Eq. (46).

To compare the accuracy between the different models, the total fission rate
density (FRD) and the fission rate root mean squared error (RMSE) are used.
The total fission rate density is the sum of the fission rate at each spatial
cell, or more precisely,

${\rm{FRD}}=\;\sum_{g=1}^{G}\left(\chi_{g}\sum_{g^{\prime}=1}^{G}\nu_{g^{\prime}}\sigma^{\rm
f}_{g^{\prime}\rightarrow g}\overline{\psi}_{g^{\prime}}\right),$ (43)

with units of cm-3s-1. To compare the wall clock times between the different
models, we compute

$\tau^{\prime}=\frac{\tau_{\rm{mg}}-\tau_{\rm{hy}}}{\tau_{\rm{mg}}},$ (44)

where $\tau_{\rm{mg}}$ and $\tau_{\rm{hy}}$ are the wall clock times required
to run the multigroup and hybrid simulations, respectively. This format was
used to show that positive wall clock time differences are for faster hybrid
method simulation times. To account for both the wall clock time and accuracy
in one metric, we use a modified figure of merit (FOM) that is commonly used
in the Monte Carlo community [figureofmerit]:

${\rm{FOM}}={\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\sout{\frac{1}{\varepsilon^{2}\tau}}}{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\uwave{\frac{1}{\varepsilon\,\tau}}},$
(45)

where $\varepsilon$ is the RMSE and $\tau$ is the wall clock time111For Monte
Carlo, $\varepsilon^{2}$ is used in place of $\varepsilon$ in Eq. (45) as it
is the variance in a statistical estimate$\varepsilon^{2}$.. A larger figure
of merit indicates a more efficient calculation.

(a)

(b)

Figure 2: The coarsening results of the (a) $G=87$ and (b) $G=361$ energy
grids as shown by the energy bounds. The naive coarsening of (a) can be seen
with the coarser grids (${\hat{G}}=10$, $21$) and how the neighboring energy
groups of ${\hat{G}}=21$ are combined to form the bounds of the ${\hat{G}}=10$
energy bounds. The involved approach with (b) combines the energy groups based
off the accuracy of the $k$-eigenvalue of a one-dimensional critical slab
problem.

###  3.1 Enriched uranium slab

The first example problem uses a 10 cm slab with a final time of $T$ = 1
$\mu$s and a time step of $h$ = 0.01 $\mu$s. Two different materials, 20%
enriched uranium and stainless steel, were used in the layout shown in Fig.
1(a). An inflow condition for the 14.1 MeV, $g=80$, energy group is prescribed
at $\mathbf{x}=0$. The inflow starts at one and decreases by a factor of two
at prescribed intervals. More specifically,

$b_{m,g=80}(\mathbf{x}=0,t)=\begin{cases}1,&\ t\in[0,0.2]\\\
0.5^{k(t)}\left(1+2\,\operatorname{erfc}\left(\frac{t-0.1(1+k(t))}{0.01}\right)\right),&\
t\in(0.2,1.0]\end{cases},\quad\text{where}\quad
k(t)=\left\lceil\frac{t-0.2}{0.1}\right\rceil$ (46)

and $t$ has units of $\mu$s.

Each drop in the input boundary induces a transient solution that eventually
returns to a quasi steady state. Because it is expected that the largest
numerical errors occur during these transient phases and for this reason, we
focus on results examine $t$ = 0.21 $\mu$s, just after the input condition
drops for the first time. These results are representative of the performance
of the methods throughout the simulation. The wall clock time is calculated
for the entire simulation and is repeated five times to obtain an average
value.

The original energy grid for the uranium and stainless steel used $G=87$
energy groups. The group coarsening is accomplished by taking a naive
approach, that assigns the same number of fine energy groups to each coarse
group, modulo differences in rounding. While it is possible to optimize the
coarsening procedure using the same techniques that were used to generate the
reference energy grid, determining the optimal group structure was not the
purpose of this research. The results of the coarsening strategy are depicted
in Fig. 2(a).

In Fig. 3, we compare the error versus wall clock time of the three methods
described at the beginning of the section, fixing either the number of groups
or the number of angles. The reference solution for these comparisons is a
grid with $G=87$ energy groups and $M=16$ angles. In Fig. 3(a), we vary the
number of coarse angles $\hat{M}<M=16$ while fixing ${\hat{G}}=G=87$ energy
groups. In Fig. 3(b), we vary the number of coarse groups $\hat{G}<G=87$,
while fixing ${\hat{M}}=M=16$ angles. As expected, the multigroup and
splitting results yield the same errors. However, the splitting runs faster
and, moreover, the relative difference in computational time increases as
$\hat{M}$ or $\hat{G}$ increases. The improvement in the splitting comes from
the different solver strategy. In particular, the uncollided equation in the
hybrid does not require inner iterations over angle to be converge, thereby
avoiding to some degree the over-solving phenomena that can unnecessarily slow
down convergence [senecal2017]. This is shown in the full multigroup and
hybrid models ($G={\hat{G}}=87$, $M={\hat{M}}=16$), in which the hybrid method
required about 161 fewer iterations per time step. The comparison between the
hybrid and splitting methods show that the hybrid generally provides improved
accuracy with a marginal cost increase. The efficiency gains from the hybrid
are substantial for intermediate values of $\hat{M}$ and $\hat{G}$, but it
tend to zero as ${\hat{M}}\to M$ or as ${\hat{G}}\to G$, as expected since the
methods are essentially the same in these limits.

(a)

(b)

Figure 3: Comparison of the RMSE error in the fission rate density vs. wall
clock time for the enriched uranium slab problem at time $t=0.21~{}\mu$s: (a)
constant numbers of energy groups with varying number of angles and (b)
constant number of angles with varying number of energy groups. As expected,
the splitting method has the same error as the multigroup method, but takes
less time to compute, with the difference in wall clock time increasing as
${\hat{G}}$ and ${\hat{M}}$ increase. The hybrid method is typically more
accurate, but slightly more expensive than the splitting method. As expected,
these differences disappear as ${\hat{G}}$ and ${\hat{M}}$ increase.

In Fig. 4, we compare the efficiency of the multigroup and hybrid methods in
computing the fission rate density over a range of energy and angle
discretization parameters. The reference solution for these comparisons
employs a grid with $G=87$ energy groups and $M=32$ angles. In Fig. 4(a), we
show the differences in errors of the two methods. For the majority of
parameters, the hybrid method performs better, with the larger differences
appearing in the coarser energy grids, but there are exceptions for larger
values of ${\hat{M}}$ when ${\hat{G}}=67$ and ${\hat{G}}=15$. These exceptions
occur when the naive coarsening strategy for the multigroup method performs
unexpectedly well. In Fig. 4(b), we show the differences in wall clock time
for the same set of simulations. Roughly speaking, we observe that the hybrid
is faster for low-resolution high-resolution calculations and the multigroup
is faster for high-resolution low-resolution calculations. In Fig. 4(c), we
plot the difference between the FOM for the hybrid method and the FOM for
multigroup method. In the large majority of cases, the hybrid FOM is better.
Exceptions to this trend occur when ${\hat{G}}=15$ or ${\hat{G}}=67$.

(a)

(b)

(c)

Figure 4: Efficiency comparison between hybrid and multigroup methods applied
to the enriched uranium slab problem across a range of discretization
parameters: (a) multigroup error minus hybrid error; (b) percent difference in
wall clock time; (c) hybrid FOM minus multigroup FOM. In (a), errors are
computed at $t$ = 0.21 $\mu$s with respect to a reference solution with $M=32$
angles and $G=87$ groups; a positive difference favors the hybrid. The hybrid
models with higher error, ${\hat{G}}=67$ and ${\hat{G}}=15$ are most likely
due to selection of the coarse groups. The errors for ${\hat{G}}=10$ are quite
large and have thus been omitted to make the remaining data easier to view.
For (b), the percent difference in wall clock time is
$(\tau_{\rm{mg}}-\tau_{\rm{hy}})/\tau_{\rm{mg}}$; positive values mean that
the hybrid is faster. In (c), the FOM difference is computed using (45);
positive values favor the hybrid.

In Fig. 5, we plot illustrative results from the hybrid and multigroup
calculations. These results demonstrate that for models with similar errors,
the hybrid method yields faster computational times and for models with
comparable wall clock times, the hybrid method yields more accurate results.

(a)

(b)

Figure 5: Illustrative comparisons of the hybrid and multigroup methods in
computing the fission rate for the enriched uranium slab problem. For similar
wall clock times, the hybrid yields smaller errors. For similar errors, the
hybrid yields smaller wall clock times.

###  3.2 Enriched and depleted uranium sphere

The second test problem involves a sphere of radius 10 cm made of three
materials: stainless steel, 20% 235U enriched uranium, and depleted uranium
containing only uranium-238. A diagram of the setup is given in Fig. 1(b). The
geometry of the problem allows it to to be modeled using a single spatial
dimension for the radial direction. The simulation is run to a final time $T$
= 1 $\mu$s with time steps of size $h$ = 0.01 $\mu$s. The boundary condition
used for the enriched uranium slab problem in Eq. (46) was also used in this
example, with the same 14.1 MeV energy group and the same decay rate, except
that it is at the sphere edge, $b_{m,g}(\mathbf{x}=10,\,t)$.

In Fig. 6, we compare the error versus the wall time for the multigroup,
splitting, and hybrid methods, fixing either the number of groups or the
number of angles. The reference solution for these comparisons is a grid with
$G=87$ energy groups and $M=16$ angles. The coarse energy groups are the same
as the previous problem and are depicted in Fig. 2(a). In Fig. 6(a), we vary
the number of coarse angles $\hat{M}<M=16$ while fixing ${\hat{G}}=G=87$
energy groups. In Fig. 6(b), we vary the number of coarse groups
$\hat{G}<G=87$, while fixing ${\hat{M}}=M=16$ angles.

Overall, we observe similar results as for the slab problem, although the
improvements in wall clock time for the splitting method over the multigroup
method are smaller and the improvements in accuracy of the hybrid over the
splitting method are greater. As in the previous problem, the multi-group
multigroup method may perform unexpectedly well in some situations because of
the coarsening. For the current problem, this situation occurs when
${\hat{G}}=80$ ${\hat{G}}=80$ and in this case, the multigroup method
outperforms the hybrid.

(a)

(b)

Figure 6: Comparison the RMSE error vs. wall clock time in the fission rate
density for the enriched and depleted uranium sphere problem: (a) constant
numbers of energy groups and varying number angles and (b) constant number of
angles and varying number of energy groups. As expected, the splitting method
has the same error as the multigroup method, but takes less time to compute,
with the difference in wall clock time increasing as ${\hat{G}}$ and
${\hat{M}}$ increase. The hybrid method is typically more accurate, but
slightly more expensive than the splitting method.

In Fig. 7, we compare the efficiency of the multigroup and hybrid methods in
computing the fission rate density over a range of energy and angle
discretization parameters. The reference solution uses $M=32$ angles and
$G=87$ energy groups is used for these comparisons. In Fig. 7(a), we show the
differences in errors of the two methods. As in the slab problem, the hybrid
performs better over a large majority of the parameters. In Fig. 7(b), we show
the differences in wall clock time for the same set of simulations. As in the
slab problem, we observe that the hybrid is generally faster for low-
resolution high-resolution calculations while the multigroup is faster for
high-resolution calculations. low-resolution calculations. To compare the
number of iterations for the full problems ($G={\hat{G}}=87$,
$M={\hat{M}}=16$), the hybrid method uses about 6729 fewer iterations per time
step than the multigroup method. In Fig. 7(c), we plot the difference between
the FOM for the hybrid method and the FOM for the multigroup method. In the
large majority of cases all cases of the uranium sphere problem, the hybrid
FOM is better.

(a)

(b)

(c)

Figure 7: Efficiency comparison between hybrid and multigroup methods applied
to the enriched and depleted uranium sphere problem across a range of
discretization parameters: (a) multigroup error minus hybrid error; (b)
percent difference in wall clock time; (c) hybrid FOM minus multigroup FOM.
Errors in (a) are computed at $t$ = 0.21 $\mu$s with respect to a reference
solution with $G=87$ groups and $M=32$ angles; a positive difference favors
the hybrid. The errors for ${\hat{G}}=10$ are quite large and have thus been
omitted to make the remaining data easier to view. For (b), the percent
difference in wall clock time is
$(\tau_{\rm{mg}}-\tau_{\rm{hy}})/\tau_{\rm{mg}}$; positive values mean that
the hybrid is faster. In (c), the FOM difference is computed using (45);
positive values favor the hybrid.

In Fig. 8, we plot illustrative results from the hybrid and multigroup
calculations. These results demonstrate that for models with similar errors,
the hybrid method yields faster computational times and for models with
comparable wall clock times, the hybrid method yields more accurate results.

(a)

(b)

Figure 8: Illustrative comparisons of the hybrid and multigroup methods in
computing the fission rate for the uranium sphere problem. For similar wall
clock times, the hybrid yields smaller errors. For similar errors, the hybrid
yields smaller wall clock times.

###  3.3 Uranium oxide slab with plane source

The last test problem involves a mixture of uranium oxide and water as it is
typically found in a light water reactor. A one-dimensional critical width
slab with $\mathbf{X}=76.7297$ cm was used with

$b_{m,g}(\mathbf{x}=0,t)=b_{m,g}(\mathbf{x}=\mathbf{X},t)=0\qquad\text{and}\qquad
f_{m,g}(\mathbf{x})=0,$ (47)

as the boundary and initial conditions. The plane source is based off of the
americium-beryllium source from the ISO 8529 standard [ambe-source] and is
adjusted for the 361 energy group grid, as shown in Fig. 9. This problem takes
time step sizes of 0.01 $\mu$s with an elapsed time of 1 $\mu$s while the
reference solution uses $G=361$ energy groups and $M=16$ angles. The group
boundaries are computed using PyNjoy [pynjoy]. Groups are coarsened in an
attempt to preserve the $k$-effective value of the slab. The results of this
more physically motivated coarsening strategy are show in Fig. 2(b).

Figure 9: The plane source used in the uranium oxide Slab problem was based
off the Americium-Beryllium (Am-Be) source from ISO 8529 [ambe-source]. The
reference source is shown by the black line while the red line is adhering the
Am-Be source with the energy grid from the $G=361$ energy group problem.

In Fig. 10, we compare the error versus wall clock time of the three methods
described at the beginning of the section, fixing either the number of groups
or the number of angles. The reference solution for these comparisons is a
grid with $G=361$ energy groups and $M=16$ angles. In Fig. 10(a), we vary the
number of coarse angles $\hat{M}<M=16$ while fixing ${\hat{G}}=G=361$ energy
groups. In Fig. 10(b), we vary the number of coarse groups $\hat{G}<G=361$,
while fixing ${\hat{M}}=M=16$ angles. The overall results are similar to the
previous problems. First the multigroup and splitting results yield the same
errors. However the splitting is faster and the relative difference in wall
clock time increases as $\hat{M}$ or $\hat{G}$ increases. When compared to the
splitting, the hybrid provides improved accuracy with a marginal cost
increase, except when ${\hat{G}}=180$ and marginally when ${\hat{G}}=120$, in
which case the multigroup method performs remarkably well.

(a)

(b)

Figure 10: Comparison the RMSE error vs. wall clock time in the fission rate
density for the the uranium oxide slab at $t=0.99~{}\mu$s: (a) constant
numbers of energy groups and varying number angles and (b) constant number of
angles and varying number of energy groups. As expected, the splitting method
has the same error as the multigroup method, but takes less time to compute.
The difference in wall clock time increases as ${\hat{G}}$ and ${\hat{M}}$
increase. The hybrid method is typically more accurate, but slightly more
expensive than the splitting method. Two notable exceptions are the cases
${\hat{G}}=120$ and ${\hat{G}}=180$ in (b), which can be attributed to how the
energy groups were coarsened for the multigroup method.

In Fig. 11, we compare the efficiency of the multigroup and hybrid methods in
computing the fission rate density over a range of energy and angle
discretization parameters. The reference solution for these comparisons
employs a grid with $G=361$ energy groups and $M=32$ angles. In Fig. 11(a), we
show the differences in errors of the two methods. For the majority of
parameters, the hybrid method performs better, especially for smaller values
of ${\hat{M}}$ and ${\hat{G}}$. The error difference is not negative in the
${\hat{G}}=180$ case, as seen in Fig 10(a) because a high angular reference
value is used ($M=32$) instead of the original reference of $M=16$. In Fig.
11(b), we show the differences in wall clock time for the same set of
simulations. Roughly speaking, we observe that the hybrid is faster for low-
resolution high-resolution calculations and the multigroup is faster for high-
resolution calculations. low-resolution calculations. This can also be seen
with the iterations needed for convergence, where the full hybrid model
($G={\hat{G}}=361$, $M={\hat{M}}=16$) requires about 2341 fewer iterations
than the full multigroup model. In Fig. 11(c), we plot the difference between
the FOM for the hybrid method and the FOM for multigroup method. In the large
majority of cases, the hybrid FOM is better, although there are some
exceptions when the number of collided energy groups is low (${\hat{G}}=30$and
${\hat{M}}=16$, $60$) and the number of collided angle is high. In this case,
the error difference at the final time step does not greatly favor the hybrid
method over the multigroup method, which causes a negative FOM difference.

(a)

(b)

(c)

Figure 11: Efficiency comparison between hybrid and multigroup methods applied
to the uranium oxide slab problem across a range of discretization parameters:
(a) multigroup error minus hybrid error; (b) percent difference in wall clock
time; (c) hybrid FOM minus multigroup FOM. Errors in (a) are computed at $t$ =
0.99 $\mu$s with respect to a reference solution with $G=361$ groups and
$M=32$ angles; a positive difference favors the hybrid. For (b), the percent
difference in wall clock time is
$(\tau_{\rm{mg}}-\tau_{\rm{hy}})/\tau_{\rm{mg}}$; positive values mean that
the hybrid is faster. In (c), the FOM difference is computed using (45);
positive values favor the hybrid. Notable exceptions with cases
${\hat{G}}=120$ in (a) and ${\hat{G}}=30$ in (c) can be attributed to how the
energy groups were coarsened for the multigroup method.

In Fig. 12, we plot illustrative results from the hybrid and multigroup
calculations. As in the previous test problems, these results demonstrate that
for models with similar errors, the hybrid method yields faster computational
times and for models with comparable wall clock times, the hybrid method
yields more accurate results. The wall clock times are shown as a percentage
of the multigroup ($G=361$, $M=16$) wall clock reference solution.

(a)

(b)

Figure 12: Illustrative comparisons of the hybrid and multigroup methods in
computing the fission rate for the uranium oxide slab problem. For similar
wall clock times, the hybrid yields smaller errors. For similar errors, the
hybrid yields smaller wall clock times.

##  4 Conclusion and future work

We have extended a discrete ordinates hybrid method of the time-dependent
neutron transport equation to the multigroup setting. The hybrid relies on a
splitting of the NTE into collided and uncollided components and takes
advantage of the fact that the collided component can be accurately
approximating using a relatively coarse discretization in angle and energy.
The resulting discretization is typically more efficient that a single
monolithic scheme.

We test the hybrid method and compare it to the usual multigroup approach that
reduces computational expense coarsening the energy groups. Except in isolated
cases, the hybrid method outperforms the standard coarsening strategy. A by-
product of the hybrid approach is a faster solver, even when the hybrid and
multigroup method use the same monolithic resolution. This improvement in wall
clock time is attributed to the fact that the uncollided component does not
require converging an inner iteration of angle; rather the inner iteration
need only be performed once. Similar results have been demonstrated in
[senecal2017], where savings are realized in a systematic way by not fully
converging inner iterations in a nested scheme.

There is more work to be done on hybrid approaches in the multigroup setting.
Previous work on the hybrid methods for single-group problems suggests that
the benefit of the hybrid with respect to angular accuracy increases in
multiple spatial dimensions. We expect this behavior to carry over to the
multigroup case, but this conjecture still has to be investigated. Previous
work has also looked into the splitting the spatial dimension into collided
and uncollided portions [crockatt2020improvements], but the division of the
temporal dimension has not been investigated. This area should be explored, as
well as the associated cost benefit of coarse collided time steps and fine
uncollided time steps. In additional, the method can be further improved by
applying acceleration techniques to the collided equation and by using
correction techniques to remove coarse grid error [crockatt2020improvements].

## Acknowledgments

The work of Ben Whewell and Ryan McClarren is supported by the Center for
Exascale Monte-Carlo Neutron Transport (CEMeNT) a PSAAP-III project funded by
the Department of Energy, grant number DE-NA003967.

The work of Cory Hauck is sponsored by the Office of Advanced Scientific
Computing Research, U.S. Department of Energy, and performed at the Oak Ridge
National Laboratory, which is managed by UT-Battelle, LLC under Contract No.
De-AC05-00OR22725 with the U.S. Department of Energy. The United States
Government retains, and the publisher, by accepting the article for
publication, acknowledges, that the United States Government retains a non-
exclusive, paid-up, irrevocable, world-wide license to publish or reproduce
the published form of this manuscript, or allow others to do so, for United
States Government purposes. The Department of Energy will provide public
access to these results of federally sponsored research in accordance with the
DOE Public Access Plan (http://energy.gov/downloads/doe-public-access-plan).

##  Appendix A Algorithms for the Multigroup and Hybrid Methods

Algorithm 1 Backward Euler update of multigroup, discrete ordinates equations.

1:$\sigma^{\rm t}_{g}$, $\sigma^{\rm s}_{g^{\prime}\to g}$, $\chi_{g}$,
$\nu_{g^{\prime}}$, $\sigma^{\rm f}_{g^{\prime}\to g}$, $c_{g}$, $q_{g}$
$\triangleright$ Material Properties

2:$\psi_{m,g}\,^{(n-1)}$ $\triangleright$ Solution from Previous Time Step

3:$h$, $\bm{\Omega}_{m}$, $w_{m}$ $\triangleright$ Discretization Parameters

4:$\varepsilon_{G},\varepsilon_{M}$ $\triangleright$ Convergence Tolerances

5:$\Delta_{G}\leftarrow 1+\varepsilon_{G}$, $j\leftarrow 0$

6:$\overline{\psi}_{g}\,^{0}\leftarrow\displaystyle\sum_{m=1}^{M}w_{m}\psi_{m,g}\,^{(n-1)}$

7:while $\Delta_{G}$ $>\varepsilon_{G}$ do $\triangleright$ Outer Iteration
($j$)

8: for $g=1,\dots,G$ do $\triangleright$ Loop over Groups

9: $\tilde{Q}_{g}\leftarrow
q_{g}+\displaystyle\sum_{g^{\prime}=1}^{g-1}\sigma^{\rm
s}_{g^{\prime}\rightarrow
g}\overline{\psi}_{g^{\prime}}\,^{j+1}+\displaystyle\sum_{g^{\prime}=g+1}^{G}\sigma^{\rm
s}_{g^{\prime}\rightarrow g}\overline{\psi}_{g^{\prime}}\,^{j}$

10:
$\tilde{Q}_{g}\leftarrow\tilde{Q}_{g}+\chi_{g}\displaystyle\sum_{g^{\prime}=1}^{g-1}\nu_{g^{\prime}}\sigma^{\rm
f}_{g^{\prime}\rightarrow
g}\overline{\psi}_{g^{\prime}}\,^{j+1}+\chi_{g}\displaystyle\sum_{g^{\prime}=g+1}^{G}\nu_{g^{\prime}}\sigma^{\rm
f}_{g^{\prime}\rightarrow g}\overline{\psi}_{g^{\prime}}\,^{j}$

11: $\Delta_{M}\leftarrow 1+\varepsilon_{M}$, $\quad\ell\leftarrow 0$

12: $\overline{\psi}_{g}\,^{j+1,0}\leftarrow\overline{\psi}_{g}\,^{j}$

13: while $\Delta_{M}>\varepsilon_{M}$ do $\triangleright$ Source Iteration
($\ell$)

14: for $m=1,\dots,M$ do $\triangleright$ Loop over Angles

15: $\tilde{Q}_{m,g}\leftarrow\tilde{Q}_{g}+\sigma^{\rm s}_{g\to
g}\overline{\psi}_{g}\,^{j+1,\ell}+\chi_{g}\nu_{g}\sigma^{\rm f}_{g\to
g}\overline{\psi}_{g}\,^{j+1,\ell}+\dfrac{1}{c_{g}h}\psi_{m,g}\,^{(n-1)}$

16:
$\psi_{m,g}\,^{j+1,\ell+1}\leftarrow\left(\dfrac{1}{c_{g}h}+\bm{\Omega}_{m}\cdot\nabla+\sigma^{\rm
t}_{g}\right)^{-1}\tilde{Q}_{m,g}$ $\triangleright$ Transport Sweep

17: end for

18:
$\overline{\psi}_{g}\,^{j+1,\ell+1}\leftarrow\displaystyle\sum_{m=1}^{M}w_{m}\psi_{m,g}\,^{j+1,\ell+1}$

19: $\Delta_{M}$
$\leftarrow\displaystyle\left|\left|\dfrac{\overline{\psi}_{g}\,^{j+1,\ell+1}-\overline{\psi}_{g}\,^{j+1,\ell}}{\overline{\psi}_{g}\,^{j+1,\ell+1}}\right|\right|_{2}$

20: $\ell\leftarrow\ell+1$

21: end while

22:
$\psi_{m,g}\,^{j+1}\leftarrow\psi_{m,g}\,^{j+1,\ell}\qquad\overline{\psi}_{g}\,^{j+1}\leftarrow\overline{\psi}_{g}\,^{j+1,\ell}$

23: end for

24: $\Delta_{G}$
$\leftarrow\displaystyle\left|\left|\dfrac{\overline{\psi}\,^{j+1}-\overline{\psi}\,^{j}}{\overline{\psi}\,^{j+1}}\right|\right|_{2}$
$\triangleright$
$\overline{\psi}\,^{j+1}=\left[\,\overline{\psi}_{1}\,^{j+1},\dots,\overline{\psi}_{G}\,^{j+1}\,\right]$

25: $j\leftarrow j+1$

26:end while

27:$\psi_{m,g}\,^{(n)}\leftarrow\psi_{m,g}\,^{j},\qquad\overline{\psi}_{g}\,^{(n)}\leftarrow\overline{\psi}_{g}\,^{j}$

28:return $\psi_{m,g}\,^{(n)}$, $\overline{\psi}_{g}\,^{(n)}$

Algorithm 2 Backward Euler update of the collision-based hybrid method for
multigroup discrete ordinates equations.

1:$\sigma^{\rm t}_{g}$, $\sigma^{\rm s}_{g^{\prime}\to g}$, $\chi_{g}$,
$\nu_{g^{\prime}}$, $\sigma^{\rm f}_{g^{\prime}\to g}$, $c_{g}$, $q_{g}$,
$\triangleright$ Uncollided Material Properties

2:$\hat{\sigma}^{\rm t}_{{\hat{g}}}$, $\hat{\sigma}^{\rm
s}_{{\hat{g}}^{\prime}\to{\hat{g}}}$, $\hat{\chi}_{{\hat{g}}}$,
$\hat{\nu}_{{\hat{g}}^{\prime}}$, $\hat{\sigma}^{\rm
f}_{{\hat{g}}^{\prime}\to{\hat{g}}}$, $c_{{\hat{g}}}$ $\triangleright$
Collided Material Properties

3:$\psi^{(n-1)}$ $\triangleright$ Solution from Previous Time Step

4:$h$, $\bm{\Omega}_{m}$, $w_{m}$, $\hat{\bm{\Omega}}_{{\hat{m}}}$,
${\hat{w}}_{{\hat{m}}}$, $\Delta E_{g}$, $\Delta{\hat{E}}_{{\hat{g}}}$
$\triangleright$ Discretization Parameters

5:$\varepsilon_{G}$, $\varepsilon_{M}$ $\triangleright$ Convergence Tolerances

6:$\psi^{\rm{u}}_{m,g}\,{}^{(n-1)}\leftarrow\psi_{m,g}\,^{(n-1)}\qquad
q^{\rm{u}}_{g}\leftarrow q_{g}$

7:for $g=1,\dots,G$ do $\triangleright$ Uncollided Flux Update

8: for $m=1,\dots,M$ do

9:
$\psi^{\rm{u}}_{m,g}\leftarrow\left(\dfrac{1}{c_{g}h}+\bm{\Omega}_{m}\cdot\nabla+\sigma^{\rm
t}_{g}\right)^{-1}\left(q^{\rm{u}}_{g}+\dfrac{1}{c_{g}h}\psi^{\rm{u}}_{m,g}\,{}^{(n-1)}\right)$
$\triangleright$ Transport Sweep

10: end for

11:end for

12:$\overline{\psi^{\rm{u}}}_{g}\leftarrow\displaystyle\sum_{m=1}^{M}w_{m}\psi^{\rm{u}}_{m,g}$

13:$q^{\rm{c}}_{{\hat{g}}}\leftarrow$
$\displaystyle\sum_{g=\alpha_{{\hat{g}}}+1}^{\alpha_{{\hat{g}}+1}}\displaystyle\sum_{g^{\prime}=1}^{G}\sigma^{\rm
s}_{g^{\prime}\rightarrow
g}\overline{\psi^{\rm{u}}}_{g^{\prime}}+\displaystyle\sum_{g=\alpha_{{\hat{g}}}+1}^{\alpha_{{\hat{g}}+1}}\chi_{g}\displaystyle\sum_{g^{\prime}=1}^{G}\nu_{g^{\prime}}\sigma^{\rm
f}_{g^{\prime}\rightarrow g}\overline{\psi^{\rm{u}}}_{g^{\prime}}$
$\triangleright$ Collided Source

14:$\overline{\psi^{\rm{c}}}_{{\hat{m}},{\hat{g}}}\leftarrow$ Algorithm 1
with: $\triangleright$ Collided Flux Update $\sigma^{\rm
t}_{g}\leftarrow\hat{\sigma}^{\rm t}_{{\hat{g}}}$, $\sigma^{\rm
s}_{g^{\prime}\to g}\leftarrow\hat{\sigma}^{\rm
s}_{{\hat{g}}^{\prime}\to{\hat{g}}}$,
$\chi_{g}\leftarrow\hat{\chi}_{{\hat{g}}}$,
$\nu_{g^{\prime}}\leftarrow\hat{\nu}_{{\hat{g}}^{\prime}}$, $\sigma^{\rm
f}_{g^{\prime}\to g}\leftarrow\hat{\sigma}^{\rm
f}_{{\hat{g}}^{\prime}\to{\hat{g}}}$, $c_{g}\leftarrow c_{{\hat{g}}}$
$q_{g}\leftarrow q^{\rm{c}}_{{\hat{g}}}$,
$\bm{\Omega}_{m}\leftarrow\hat{\bm{\Omega}}_{{\hat{m}}}$,
$w_{m}\leftarrow{\hat{w}}_{{\hat{m}}}$, $\psi^{(n-1)}\leftarrow 0$,
$\varepsilon_{G}\leftarrow\varepsilon_{G}$, and
$\varepsilon_{M}\leftarrow\varepsilon_{M}$

15:$q^{\rm{t}}_{m,g}\leftarrow\dfrac{\Delta
E_{g}}{\Delta{\hat{E}}_{{\hat{g}}}}\left[q^{\rm{c}}_{{\hat{g}}}+\displaystyle\sum_{{\hat{g}}^{\prime}=1}^{{\hat{G}}}\sigma^{\rm
s}_{{\hat{g}}^{\prime}\rightarrow{\hat{g}}}\overline{\psi^{\rm{c}}}_{{\hat{g}}^{\prime}}+\chi_{{\hat{g}}}\displaystyle\sum_{{\hat{g}}^{\prime}=1}^{{\hat{G}}}\nu_{{\hat{g}}^{\prime}}\sigma^{\rm
f}_{{\hat{g}}^{\prime}\rightarrow{\hat{g}}}\overline{\psi^{\rm{c}}}_{{\hat{g}}^{\prime}}\right]$
$\triangleright$ Total Source

16:for $g=1,\dots,G$ do $\triangleright$ Total Flux Update

17: for $m=1,\dots,M$ do

18:
$\psi^{\rm{t}}_{m,g}\leftarrow\left(\dfrac{1}{c_{g}h}+\bm{\Omega}_{m}\cdot\nabla+\sigma^{\rm
t}_{g}\right)^{-1}q^{\rm{t}}_{m,g}$ $\triangleright$ Transport Sweep

19: end for

20:end for

21:$\psi_{m,g}\,^{(n)}\leftarrow\psi^{\rm{t}}_{m,g}\qquad\overline{\psi}_{g}\,^{(n)}\leftarrow\displaystyle\sum_{m=1}^{M}w_{m}\psi^{\rm{t}}_{m,g}$

22:return $\psi_{m,g}\,^{(n)}$, $\overline{\psi}_{g}\,^{(n)}$