3.1. The Collision-Based Hybrid Method#
Collision-based hybrid algorithms have been explored and implemented in neutron [HM13] and radiative [HH21] transport problems. This hybrid method is pursued to improve the computational time while keeping a similar discretization error to high fidelity grids. Separating this equation into a fine, uncollided grid and a coarse, collided grid allows for faster wall clock times while keeping the error lower than traditional methods. The collision-based hybrid method has been implemented for the neutron transport equation for mono-energetic and multigroup problems and expands on previous research [HM13, WMHS23].
3.1.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 (\(S_N\)) method [Cha13, LM93]. 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, which is presented in the source-iteration algorithm. 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. A smaller energy group structure must be carefully chosen to limit the discretization error and preserve the characteristics of the energy-dependent cross section [MMB+17]. 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 \(S_N\) calculations. There are techniques for reducing these ray effects [FKCH20, HH19], but the most common solution is to increase the number of angles [LM93], thus increasing the computational time.
Acceleration techniques, such as coarse mesh rebalancing and diffusion synthetic acceleration (DSA), are used to reduce the computational cost of the iterative solver. These methods cannot be implemented indiscriminately, as coarse mesh rebalancing must be concerned with the coarse mesh size [LM93] while diffusion synthetic acceleration has difficulty in highly heterogeneous materials [SHH21]. For the outer iterations, there are two-grid, nonlinear diffusion acceleration, and Krylov subspace schemes that can improve the convergence of problems with upscattering or fission [AM93, Ani13, SRZP+18]. While these methods demonstrate improvements with high upscattering materials, there is a nominal improvement over Gauss-Seidel in low upscattering materials [AM93]. 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.
More recently, collision-based hybrid algorithms for time-dependent transport equations have been important topics for transport researchers [CCGH17, CCGH19, CCH20, HM13, HH21, WMHS23]. These algorithms split the transport equation into collided and uncollided components, as was done much earlier in the steady-state context [Alc77]. In the aforementioned examples, the hybrid approach yielded substantial improvements in efficiency when compared to monolithic discretization schemes. While previous work has applied different spatial discretization schemes to the collided and uncollided components, the approach has so far only been implemented on one-dimensional problems. In this analysis, the collision-based hybrid algorithm is extended to the two-dimensional multigroup setting. The algorithm is first applied to one-dimensional slabs and spheres with different boundary conditions and sources before including two-dimensional examples using a Cartesian grid. The hybrid algorithm is shown to be more accurate for the same computational time or reach solutions with the same error level while using significantly less time. As shown in [WMHS23], the application of a collided-uncollided splitting to a monolithic discretization grid recovers the same solution as standard Gauss-Seidel in less time. This wall clock improvement is due to the changes in the solver that reduce unnecessary iterations in the inner loop of the nested approach, as explored in [SJ17].
This chapter is organized as follows. The continuous and discrete collided-uncollided split and the collision-based hybrid algorithm are presented in the hybrid formulation, building on the continuous and multigroup discrete ordinates transport equations found in The Neutron Transport Equation. Numerical results for one-dimensional test problems, initially presented in [WMHS23], are given in One Dimensional Collision-Based Hybrid Method and compare the hybrid method to more traditional approaches. Two-dimensional results and conclusions, which demonstrate the abilities of the hybrid method in more computationally expensive examples, are presented in later sections.
3.1.2. The Collision-Based Hybrid Formulation#
The goal of the hybrid approach is to accelerate the computation of a numerical solution for the neutron transport equation (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. Since it is linear, the neutron transport equation (1) can be separated into two equations where a distinct angular flux is calculated in each and the results added together to obtain the total angular flux. For the nomenclature, the collided portions of the neutron transport equation are represented with a superscript \(\mathrm{c}\) while the uncollided portions are represented with a superscript \(\mathrm{u}\).
3.1.2.1. The Continuous Collided-Uncollided Split#
The first equation governs the uncollided angular flux, \(\Psiu\). It takes the form
with \(\Qu = Q\). The initial condition for \(\Psiu\) is
and the incoming boundary data is
The second equation governs the collided angular flux, \(\Psic\). It takes the form
where the isotropic term \(\Qc\) comes from the scattering and fission source of the uncollided flux:
The initial condition for \(\Psic\) is
and the incoming boundary data is
A third equation for the total angular flux \(\Psit\) takes the form
where the total external source \(\Qt\) is
The initial condition for \(\Psit\) is
and the incoming boundary data is
The continuous neutron transport equation implies that \(\Psit = \Psiu + \Psic = \Psi\), making an independent equation for \(\Psit\) trivially redundant. This is no longer the case when equations for \(\Psit\), \(\Psiu\), and \(\Psic\) are discretized with different methods and/or grids. In practice, the hybrid method seeks at each time step a numerical approximation for \(\Psiu\) on a fine grid and a numerical approximation for \(\Psic\) on a coarse grid. In the original formulation [HM13], these two approximate solutions were, at the end of the time step, combined via a remapping procedure in which an approximation for \(\Psic\) 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 \(\Psit\) to remap onto the fine grid was introduced in [CCH20] since the source \(\Qt\) 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 fission sources that appear in the right-hand side of the neutron transport equation (1).
The choices of initial conditions, boundary conditions, and sources for \(\Psiu\) and \(\Psic\) are not unique. The strategy above is to assign data to the uncollided equation since it is equipped with the finest discretization. This strategy may not always be the best choice, particularly in highly-collisional problems with strong boundary layers [Den06].
The hybrid method is summarized in the backward Euler algorithm below.
Algorithm — Backward Euler Time-Dependent Collision-Based Hybrid Method (multigroup, discrete ordinates)
3.1.2.2. 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. It can be assumed 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 \(\hG\) and \(\hM\) angles. To differentiate between different discretization parameters, a hat adornment \((\hat{\cdot})\) is used for the collided parameters. The weights and angles in the discretization of the collided equations are denoted by \(\hw_{\hm}\) and \(\hOmega_{\hm}\), respectively, for \(\hm = 1, \dots, \hM\). To define the coarse groups, let
be a set of \(\hG + 1\) integers and set \(\hE_{\hg} = E_{\Gamma_{\hat g}}\). Then for each \(\hat{g} \in \chG\),
The uncollided equation for
is
where \(\qu_{m,g} = q_{m,g}\). Equation (2) is solved over a time-step \([t^n, t^{n+1})\) with initial condition
and incoming boundary data
The collided transport equation for
is
where
the energy-coarsened cross-sections are given by
and the isotropic source is
with
We solve Equation (3) over a time-step \([t^n, t^{n+1})\) with initial condition
and boundary data
The total flux equation for
is
where
and \(\hg\) is the unique integer such that \(\Gamma_{\hg-1} + 1 \leq g \leq \Gamma_{\hg}\) or, equivalently, \((E_{g-1},E_{g}) \subset (\hE_{\hg-1},\hE_{\hg})\). We solve Equation (7) over a time-step \([t^n, t^{n+1})\) with initial data
and incoming boundary condition
In summary, Equations (2), (3), and (7) are solved in succession for each time step. The solution of Equation (7) at the end of the time step provides the initial condition needed by Equations (2) and (7) at the next time step to consistently join the uncollided and collided calculations. The initial condition for Equation (3) is set to zero for each time step.
The solution procedure for solving the collided flux is similar to the approach used for solving the multigroup discrete ordinates neutron transport equation as shown in the time-dependent algorithm. The hybrid iteration for a backward Euler temporal discretization scheme is shown in the backward Euler algorithm. It relies 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 [BG70, LM93]. The uncollided and the total flux updates do not require any iterations because the right-hand side of the respective equation is fixed.
3.1.2.3. Metrics for the Performance of the Hybrid Method#
The purpose of the collision-based hybrid method is to reduce the required amount of computational time by coarsening the collided energy and angular grids. The separation of the collided and uncollided portions of the neutron transport equation is performed to increase the accuracy of these low fidelity models as compared to a standard coarsening procedure. To show both the wall clock speed up and the accuracy of the hybrid method over standard techniques, three solution strategies are used for each test problem.
The first is referred to as the “Multigroup” solution, which is the traditional neutron multigroup method with \(G\) groups and \(M\) angles. It uses Gauss-Seidel with a tolerance of \(\varepsilon_G = 1 \times 10^{-8}\) for the outer iteration and source iteration with a tolerance of \(\varepsilon_M = 1 \times 10^{-12}\) for the inner iteration, as described in the source-iteration algorithm. The coarsening strategy in Equations (4) and (5) is used for lower fidelity multigroup models. 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 [BG70]. There can be cancellation of errors in integrated quantities, with these behaviors observed in some test cases.
The second solver is the “Hybrid” method as described in this section. It uses \(G\) groups and \(M\) angles in the discretization of the total and uncollided angular fluxes, and it uses \(\hG\) groups and \(\hM\) angles in the discretization of the collided angular flux. The hybrid method employs the algorithm described in the backward Euler algorithm with iteration tolerances of \(\varepsilon_G = 1 \times 10^{-8}\) and \(\varepsilon_M = 1 \times 10^{-12}\) for the outer and inner iterations, respectively.
The final solver used is the “Splitting” method, which is the same as the hybrid method described in the backward Euler algorithm. The same iteration tolerances (\(\varepsilon_G = 1 \times 10^{-8}\) and \(\varepsilon_M = 1 \times 10^{-12}\)) are used but the energy groups and angles remain consistent for the collided, uncollided, and total flux (\(G = \hG\) and \(M = \hM\)). The splitting method demonstrates the over-solving explored in [SJ17] and exemplified in [WMHS23]. The accuracy of the splitting method is equivalent to the multigroup method and the accuracy will not be included in either the one- or two-dimensional results. This method is included in the results to show the inherent speed up of using the hybrid method over the multigroup method.
To compare the accuracy between the different models, various metrics are used. The first is the total fission rate density (FRD) as described by
which is the sum of the fission rate at each spatial error with units of \(\mathrm{cm}^{-3}\,\mathrm{s}^{-1}\). The scattering rate density is also used, in which the fission terms of Equation (8) are replaced with the scattering matrix \(\sig{s}\). In certain instances, for the two-dimensional results, the scalar flux is summed over all the energy groups to account for regions that have neither scattering or fission. To condense the accuracy into one term, the root mean squared error (RMSE) is used with the fission rate density, scatter rate density, or scalar flux. There are also cases where the energy groups are broken down into thermal, epithermal, and fast energy regions. The thermal region includes all energy groups below 1 keV, the epithermal region is between 1 keV and 1 MeV, and the fast region includes all energy groups above 1 MeV. The error difference subtracts the hybrid error from the multigroup error, meaning that for positive values, the hybrid method is more accurate.
To compare the wall clock times between the different models, the formula is
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 calculate the wall clock time effectively, these times were recorded five separate times and averaged. The one-dimensional wall clock results used a personal computer while the two-dimensional results used the high processing computers at Lawrence Livermore National Laboratory.
To combine both the wall clock time and accuracy into a single metric, a modified figure of merit (FOM) is used as
where \(\varepsilon\) is the RMSE and \(\tau\) is the wall clock time. This is a metric commonly used in the Monte Carlo community [VM22] where \(\varepsilon^2\) is used in place of \(\varepsilon\) in Equation (10) as it is the variance in a statistical estimate. Since a larger figure of merit indicates a more efficient calculation, the FOM difference subtracts the multigroup FOM from the hybrid FOM. This is formulated to show that, for positive FOM difference, the hybrid method is more desirable than the multigroup method.