A COMPARISON OF STOCHASTIC MESH CELL VOLUME COMPUTATION STRATEGIES FOR THE RANDOM RAY METHOD OF NEUTRAL PARTICLE TRANSPORT

,


Introduction
The Random Ray Method (TRRM) of neutral particle transport is a technique for numerically estimating the solution to the Boltzmann neutron transport equation.It shares many similarities with the Method of Characteristics (MOC), though the key difference is that TRRM uses a continuous rather than discrete treatment for the spatial and angular dependencies of the angular flux.Instead of using a fixed deterministic quadrature as in traditional MOC, TRRM samples rays from a uniform random distribution in space and angle and tracks them through the simulation domain, solving the MOC equations along the way.While random ray does introduce the complexities of stochastic uncertainties, there are several key advantages of TRRM as compared to traditional deterministic MOC: John Tramm et al., Stochastic Mesh Cell Volume Computation Strategies for TRRM • Extreme reduction in memory usage requirements.This advantage stems from the fact that in TRRM, no storage of boundary angular fluxes is required.
• TRRM can reach an unbiased stationary distribution using only a small fraction of the track/ray density that would be required for a deterministic solver.The result is that a TRRM power iteration takes orders of magnitude less work than a deterministic one while achieving similar fidelity.This allows TRRM to converge solutions using naïve power iteration in time that is highly competitive with traditional deterministic MOC [1] that uses coarse mesh finite difference (CMFD) [2] acceleration.
TRRM has been implemented in the Advanced Random Ray Code (ARRC) [3], and has been demonstrated to be capable of converging a full core high-fidelity 3D benchmark problem [1].More information on the theory and derivation of TRRM can be found in [4], [3], and [5].
In the present study, we will be looking deeper into one of the statistical issues that arises in the derivation of TRRM, which stems from TRRM's stochastic rather than deterministic basis.In particular, we will derive the scalar flux expression that involves estimation of the volume of each mesh cell, and show that direct use of this expression results in undersampling bias when randomness is introduced.We will then propose a new method to fix the bias and quantify its impact on simulation accuracy.

Derivation
While the MOC and random ray equations have been derived in many other works (e.g., [6], [7], [8]), we provide a cursory derivation here as their form is crucial for understanding how cell volumes are used in MOC and TRRM, and how errors might arise from their various estimation techniques.We begin with the steady-state Boltzmann Neutron Transport partial differential equation (PDE) (Equation 1), where all symbols have their usual meaning (for instance, as defined in [5]).
If the right hand side of Equation 1 is condensed into a single term, represented by the total neutron source term Q( r, Ω, E), the form given in Equation 2 is reached.
streaming term Fundamentally, MOC works by decomposing the Boltzmann PDE into a family of ordinary differential equations (ODEs), each of which solves for the angular flux only along a single Proceedings of the PHYSOR 2020, Cambridge, United Kingdom Physics of Reactors Transition to a Scalable Nuclear Future characteristic line, as in Equation 2, thus altering the full spatial and angular scope of the transport equation into something that holds true only for a particular linear path (also called a "track" or "ray") through the reactor.These tracks are linear for neutrons as they are neutral particles and are therefore not subject to field effects.If enough of these characteristic ODEs are solved for, we can approximate the solution to the governing Boltzmann PDE.To accomplish this, we parameterize r (the neutron's spatial location) with respect to some reference location r 0 such that r = r 0 + s Ω.In this manner, Equation 2 can be rewritten for a specific segment length s at a specific angle Ω through a constant cross section region of the reactor geometry as in Equation 3.
As this equation holds along a 1D path, we can assume the dependence of s on r 0 and Ω such that r 0 + s Ω simplifies to s.When the differential operator is also applied to the angular flux ψ, we arrive at the characteristic ODE given in Equation 4.
An analytical solution to this characteristic ODE can be achieved with the use of the integrating factor e s 0 ds Σ T (s ,E) to arrive at the final form of the characteristic equation shown in Equation 5.
With this characteristic equation, we now have an analytical solution along a linear path through any constant cross section region of a reactor.Similar to many other solution approaches to the Boltzmann neutron transport equation, the MOC and random ray approaches also use a "multi-group" approximation in order to discretize the continuous energy spectrum of neutrons travelling through the reactor into fixed set of energy groups G, where each group g ∈ G has its own specific cross section parameters.
Following the multi-group assumption, another assumption made is that a large and complex problem can be broken up into small constant cross section regions i [9], and that these regions have group dependent isotropic sources (fission + scattering), Q i,g .Anisotropic sources are also possible with MOC based methods [10], but are not used here for simplicity.With these key assumptions, the multi-group characteristic equation can be written as in Equation 6.
In ARRC, a constructive solid geometry (CSG) definition of the reactor is used to create spatially defined source regions.These neutron source regions are often approximated as being constant (flat) in intensity of source, but can also be defined using a higher order source (linear, quadratic, etc.) that allows for fewer source regions to be required to achieve a specified solution Proceedings of the PHYSOR 2020, Cambridge, United Kingdom John Tramm et al., Stochastic Mesh Cell Volume Computation Strategies for TRRM fidelity [11,10].In this study, we will use the normal approximation of a spatially constant isotropic fission and scattering source Q i,g in each flat source region (FSR) i and for energy group g, which allows us to solve Equation 6and to arrive at Equation 7.
As we will be integrating each track (or ray) through one FSR mesh cell at at time, we can rewrite our equation into a form where we will be solving for the outgoing angular flux ψ out g given the incoming angular flux ψ in g as in Equation 8, where r now represents the length of the ray's path through FSR mesh cell i.
Using Equation 8, we can now compute the angular flux for a track as it traverses through a reactor geometry and passes through different FSRs.However, to compute the scattering and fission sources for the next iteration, we need to be able to compute the average scalar flux for each cell, and to do that we need to be able to determine the total average angular flux in each cell.To track these quantities, we will begin by finding a solution for the average angular flux ψ r,i,g for a single ray r passing through region i and for energy group g, as in Equation 9.
By substituting Equation 7 into Equation 9and solving, we can determine the average angular flux for a single ray through the region i using Equation 10.
By re-arranging Equation 8, we can define ∆ψ r,g as the change in angular flux for ray r passing through region i: If we substitute Equation 11 into Equation 10, we arrive at: Because scalar flux is defined as: Proceedings of the PHYSOR 2020, Cambridge, United Kingdom Physics of Reactors Transition to a Scalable Nuclear Future we can estimate the scalar flux φ i,g in each region i and energy group g by taking a weighted average (by distance travelled) over all N i rays passing through that region in a power iteration, where r is the contribution from a single ray, as in Equation 14.
By substituting Equation 12into Equation 14 and simplifying, we arrive at a (naïve) expression for each iteration's scalar flux estimate: The estimate of φ naive i,g defined in Equation 15 is computed once per power iteration by selecting rays randomly from a uniform distribution in space and angle and following them for a set distance through the problem geometry, with rays solving Equation 11 for all energy groups as they pass through each FSR and tallying ∆ψ r,g and r contributions along their way.The scalar flux estimates produced in each iteration's stochastic transport "sweep" are then used in several ways in the power iteration process.Namely, to compute the neutron source Q i,g used in the next power iteration and to update the eigenvalue estimate.Full definitions and derivations of the entire power iteration, source update, and eigenvalue computation processes are given in [5], and are not included here for the sake of brevity.

Scalar Flux Ratio Estimators
While the naïve estimator defined in Equation 15 is numerically sound for a deterministic process (e.g., classical MOC), issues arise when used in a stochastic process like TRRM.Specifically, the second term of the equation (which is the ratio of the change in angular flux for all rays through the cell to the total distance travelled by all rays in the cell) features stochastic variables in both the numerator and denominator, making it a biased estimator [12].
We have several options for how to treat the naïve biased estimator defined in Equation 15, which will be presented in the following sections of this paper.To test the different methods, we will implement them in ARRC [3] and use them to simulate the standard 2D C5G7 light water reactor benchmark problem [13] that features four assemblies and a large moderator reflector region.While ARRC is capable of much larger 3D full core simulations, the relatively small 2D C5G7 Proceedings of the PHYSOR 2020, Cambridge, United Kingdom John Tramm et al., Stochastic Mesh Cell Volume Computation Strategies for TRRM problem is used in this analysis as it is fast to run, has defined 7 group macroscopic cross sections, and features a reference solution generated via multigroup Monte Carlo.For all tests, we used a mesh featuring 142,964 FSRs, ray lengths of 250 cm, and a ray dead zone [4] of 25 cm.For all tests, 800 inactive iterations were used to converge the source to a stationary distribution before accumulating statistics in the active region of the simulation.All runs were converged to a target uncertainty threshold of 0.003 source averaged relative error (SARE), which means that the total fission source (integrated across all 7 energy groups in the problem) in each FSR mesh cell will be converged to 0.3% (on average) for all FSR mesh cells.

Naïve Estimator
The first option is to simply use the estimator defined in Equation 15 verbatim and implement it directly in a random ray power iteration simulation.The results of testing the naïve estimator on the 2D C5G7 sample problem are shown in Table 1.Our findings show that the naïve estimator does indeed produce a biased eigenvalue solution when low ray densities are used, but the eigenvalue approaches the true reference solution as the ray density is increased.Increases in the ray density are directly proportional to the number of ray crossings in any given FSR (N ).Thus, our naïve estimator can be said to be asymptotically consistent.As shown in Table 1, if we increase the ray density (and thus N ) the bias eventually diminishes to acceptably low levels.This phenomenon is very similar to the undersampling biases found in Monte Carlo simulations when very few particles are used [14,15].While it is good that the estimator is at least consistent, the downside is that we have to run a large number of rays per power iteration to achieve acceptable accuracy.A more ideal estimator with reduced bias might be able to get by with a much coarser ray density and therefore a much faster time to solution.1.18646 ± 0.00020 1.18664 ± 0.00020 1.17911 ± 0.00114 5,000 1.18658 ± 0.00020 1.18692 ± 0.00019 1.18192 ± 0.00104 10,000 1.18667 ± 0.00017 1.18663 ± 0.00021 1.18473 ± 0.00094 20,000 1.18685 ± 0.00018 1.18651 ± 0.00021 1.18713 ± 0.00121 40,000 1.18647 ± 0.00021 1.18647 ± 0.00017 1.18669 ± 0.00097 80,000 1.18689 ± 0.00018 1.18663 ± 0.00022 1.18659 ± 0.00035 160,000 1.18676 ± 0.00013 1.18659 ± 0.00016 1.18614 ± 0.00048

Simulation Averaged Estimator
One method for reducing the bias in a ratio estimator is to increase the number of samples for the stochastic variable in the denominator.This approach has recently been deployed to reduce the Proceedings of the PHYSOR 2020, Cambridge, United Kingdom Physics of Reactors Transition to a Scalable Nuclear Future bias for ratio estimators that arise in the stochastic gradient descent algorithm from the field of machine learning, as discussed in [16].
Luckily, in the context of a random ray simulation, improvement in the number of samples for the denominator variable is easily accomplished.Instead of computing a new stochastic tally for (the distance estimator in the denominator) each power iteration, we can convert this to a "simulation average" quantity.As the average distance travelled by rays through a cell is a function only of geometry (and not a function of neutronics or the state of the neutron distribution in the system), we can begin tallying this quantity in the first power iteration and keep adding to the tally throughout the simulation.As more power iterations are run, the sample size of the distance tally will continue to grow, therefore continually reducing the overall bias of the scalar flux estimator.It is important to note that we could not use this simulation averaged strategy with the numerator tally for ∆ψ r,g , as that value is dependent on the source and neutron distribution of the system, which evolves between power iterations as it moves towards a stationary distribution over the course of the inactive phase of the simulation.
Thus, if we define M i as the number of ray crossings in a given FSR i over the course of the entire simulation, we can write a new "simulation averaged" estimator as: where each power iteration M i grows by approximately N i (which varies stochastically between iterations), and after many power iterations M i N i .
We tested our new simulation averaged estimator on the 2D C5G7 problem and compared it against the original naïve estimator.The results, given in Table 1, show that the new simulation averaged estimator virtually eliminates the bias even for very low ray densities.The new estimator therefore allows us to run with the coarsest possible ray density, which is typically set as the lowest number of rays required to hit (almost) all FSRs at least once each power iteration * .Use of a coarse ray density is highly advantageous as it minimizes time spent in the inactive phase of the simulation.* In the case of the 2D C5G7 benchmark problem, this minimal density occurs when approximately 1,250 rays per iteration are used, resulting in a mean N for each FSR of about 14.This value is arrived at empirically.If fewer rays are used, the iteration costs of building the source can become more costly than the transport sweep itself, resulting in an overall increase in runtime.Additionally, if very few rays are used per iteration (such that most FSRs are not being hit each iteration at all), the simulation can become unstable.
Proceedings of the PHYSOR 2020, Cambridge, United Kingdom John Tramm et al., Stochastic Mesh Cell Volume Computation Strategies for TRRM

Analytical Volume Estimator
In the event that the true analytical volumes are known for each FSR cell V i a priori, then a constant expression can be derived for the denominator of Equation 15.If the distance travelled by each ray through the global simulation domain is D, and R total rays are run in each power iteration, then the exact expected distance of rays travelling through a given FSR cell i becomes a constant value that can be computed analytically as: In this manner, our estimator for the scalar flux is no longer a ratio estimator, as the denominator becomes a constant value: Due to the complexities of user-defined constructive solid geometry (CSG) definitions of reactor geometry objects, we typically do not know the analytical volumes of each FSR cell analytically, which makes this estimator impractical for realistic use cases.However, for the purposes of this analysis, we were able to manually calculate the volumes of each FSR in our 2D C5G7 mesh in ARRC.Using these true analytical volumes, we performed a ray density convergence analysis (shown in Table 1) and found that the behavior of the analytical mode was nearly identical to that of the simulation averaged estimator.

Changes in Convergence Rate
In the previous section, we showed that the new simulation averaged estimator was able to virtually eliminate the bias from the simulation even with a very coarse ray density.However, we noticed that use of this estimator also increased the number of active iterations required to converge the solution below our target uncertainty threshold.These results, given in Table 2, show that approximately 5x more active iterations are required to converge statistics in the active phase of the simulation with equal ray densities, meaning that the unbiased simulation averaged estimator comes with a cost of significantly higher variance.This is likely due to the fact that the numerator and denominator in Equation 15are not fully independent.As the change in angular flux (defined in Equation 11) is a function of the length travelled by the ray through the FSR, these two values will be naturally correlated.Figure 1 demonstrates this correlation by showing the numerator and denominator values tallied for Equation 15 for all iterations in the active region of a 2D C5G7 simulation, for a single (randomly chosen) FSR and energy group.Thus, by using the correct analytical value or simulation averaged estimator for the denominator, the correlation is ignored and the variance of the ratio as a whole is increased even though the bias is eliminated.
Given the differences in variance, it is therefore important to compare the overall runtimes required for the different methods to achieve a target convergence level and accuracy.we measured the runtimes of the simulation when using each estimator in a configuration that was able to result in less than 25 pcm overall bias in eigenvalue.For the naïve estimator, this meant using 20,000 rays per iteration, while for the simulation averaged and analytical methods we only needed to use the minimum ray density of 1,250 rays per iteration.For all estimators, we used 800 inactive cycles to reach a stationary source distribution (as defined by convergence of Shannon Entropy [3]).Once stationary, the naïve estimator required 20 active iterations to converge statistics, resulting in an overall runtime of 128 seconds.The simulation averaged estimator required 1,675 active iterations to converge statistics, resulting in an overall runtime of 33 seconds.The analytical method behaved similarly in both numerical and runtime performance as the simulation averaged estimator.Thus, for comparable accuracy, the simulation averaged and analytical methods were about 4x faster in runtime even though they had higher variances, due to the fact that they could run with 16x fewer rays per iteration and therefore could get through the inactive iterations much faster.The error in pin power distributions was also improved.At 20,000 rays per iteration the naïve estimator achieved an absolute average percent pin power error (APE) It is important to consider that the 2D C5G7 test case that was used for this analysis has a fairly low dominance ratio which only required 800 power iterations to reach source stationarity.Full core 3D reactor simulation problems can have much higher dominance ratios that can require many thousands of iterations to reach stationarity -meaning that the increased ray density requirements of the naïve estimator would greatly increase the numerical performance advantage of the simulation averaged estimator.For example, in [1], 11,763 inactive iterations and 1,570 active iterations (using the simulation averaged estimator) were required to converge a full core 3D light water reactor problem with a high-fidelity mesh.If the naïve estimator were used to cut the active cycle count by 85x, it would only reduce the total iteration count by approximately 12% while increasing the required ray density by 16x.Conversely, for very small problems with low dominance ratios (e.g., a 2D pincell), the naïve estimator may provide more competitive performance, though with the downside that a convergence study would have to be performed to ensure that the solution is not biased.
It is also important to consider the effect of the user's target uncertainty on estimator performance.For instance, if a very low level of uncertainty is desired (e.g., less than 1 pcm in eigenvalue), then the number of active iterations can become much larger than the number of inactive iterations, even for high dominance ratio problems.In such cases, the lower variance of the naïve estimator may have a much larger benefit on overall numerical performance.However, it is usually not worth converging MOC and random ray simulations to such tight levels due to other sources of error inherent in the simulation that would not be eliminated by doing so (e.g., nuclear data uncertainties, multigroup cross section generation, transport correction errors), which can add up to 100 pcm in eigenvalue [17].Thus, the transport convergence targets used in this analysis (in the range of 20 pcm in eigenvalue) represent the tighter bounds for real world engineering usage.

Conclusions
In this study, we derived the characteristic equations used in TRRM and used them to define an estimator for the scalar flux in each source region cell of a simulation.We found that our estimator was in fact a ratio estimator having stochastic variables in both the numerator (where the change in angular flux ∆ψ is tallied) and denominator (where the total distance travelled by all rays through the cell is tallied).We tested this naïve estimator out on the 2D C5G7 benchmark problem and found that the estimator exhibited clear undersampling bias, in one case resulting in over 1,200 pcm eigenvalue error.However, we found the estimator to be consistent, meaning that the undersampling bias could be mitigated if very high number of rays per iteration were used in a simulation.
We then investigated a method that reduces bias of ratio estimators by increasing the number of samples taken for the denominator stochastic variable, which has shown to be successful for other estimators commonly used in the field of machine learning.This strategy allowed us to develop a new "simulation averaged" scalar flux estimator, where distances travelled through each FSR are tallied continuously throughout the simulation rather than being discarded each power iteration.Empirical testing of the new estimator revealed no measurable bias, meaning that it could be used even when very few rays were run per power iteration.Additionally, we found that the simulation Proceedings of the PHYSOR 2020, Cambridge, United Kingdom Physics of Reactors Transition to a Scalable Nuclear Future averaged volume estimator produced results nearly identical to a method that used the true analytical volumes.However, the new simulation averaged estimator did result in an increase in variance compared to the naïve estimator.Overall, the new unbiased simulation averaged method resulted in 2D C5G7 runtimes that were overall about 4x faster than the naïve estimator for a given accuracy due to the naïve method requiring orders of magnitude more rays per iteration to eliminate bias.
All work to date published by the authors using TRRM has used the simulation averaged estimator, though this has not been documented anywhere previously as it was only discovered in this analysis that using the naïve estimator would result in undersampling bias.These findings are an important detail of random ray neutron transport; those implementing TRRM solvers of their own are highly encouraged to use the simulation averaged estimator so as to mitigate undersampling biases.

Figure 1 :
Figure 1: Relationship between the numerator and denominator of Equation 15 of a single FSR and energy group for each active iteration in a test simulation.

Table 1 :
Eigenvalue and standard deviations of different estimators.Bold indicates entries with error from the reference solution greater than 3 standard deviations.

Table 2 :
Active iteration requirements to reach target source uncertainty level.
To test this, Proceedings of the PHYSOR 2020, Cambridge, United Kingdom Physics of Reactors Transition to a Scalable Nuclear Future Proceedings of the PHYSOR 2020, Cambridge, United Kingdom John Tramm et al., Stochastic Mesh Cell Volume Computation Strategies for TRRM value of 0.454%, while at 1,250 rays per iteration the simulation averaged estimator and analytical methods achieved APE values of 0.367% and 0.394% respectively.