ENHANCING THE PERFORMANCE OF FISSION SOURCE CONVERGENCE WITH FUNCTIONAL EXPANSION TALLIES IN SERPENT 2 MONTE CARLO CODE

This paper presents an upgrade to the built-in response matrix based solver implemented in Serpent 2 Monte Carlo code aiming to improve the fission source convergence when obtaining the forward solution to the k-eigenvalue criticality source problems. The functional expansion tallies are introduced in an attempt to improve the accuracy of the cellwise form factors that feed the response matrix solver, replacing the current mesh-based approach. The functional expansion tallies reconstruct the binning surface and collision tallies, by using high-order series expansion to represent the original and continuous spatial distributions. This new feature is implemented to Serpent 2 and tested by singleassembly and full-core PWR calculations (BEAVRS benchmark). The results show enhanced performance of the convergence acceleration methodology based on an improved initial guess of the fission source.


INTRODUCTION
Regardless the continuous advances in high performance computing, Monte Carlo methods nevertheless are expensive, not only for large systems, also, for industry basis design-analysis. The goal is to develop strategies that reduce the execution time without sacrificing accuracy.
Monte Carlo neutron transport codes solve k-eigenvalue problems from an unknown fission source. Simulations begin with a guess of the fission neutron source distribution; and, through a series of so-called inactive cycles, from which results are not collected, it is expected to achieve source convergence. The efficiency of the Monte Carlo calculations would be improved by evolving the current techniques to accelerate the convergence of the fission source, and for extension, reduce the number of inactive cycles.
Monte Carlo fission source error can be broken down into three components: the initial guess distribution, a statistical error and a systematic bias error. The statistical error of order O (1/ √ m), is due to the limited number of neutron histories simulated per cycle, m. The systematic bias error of order O (1/m) is caused by normalization of the source size at each cycle. Unlike the statistical error, the bias does not decay over the cycles. The error in the source sampled in the initial cycle propagates into sources in the successive cycles, decaying over cycles; although, it still represents a major part of the total error in large-scale systems and systems with high dominance ratio. Several deterministic techniques have been applied to accelerate Monte Carlo source convergence, by reducing the number of inactive cycles used to converge the fission source: e.g., fission matrixbased methods (Kitada and Takeda [1], Carney et al. [2], Pan et al. [3], Dufek and Tuttelberg [4]), Wielandt-based methods (Yamamoto and Miyoshi [5], She et al. [6]), hybrid CMFD/Monte Carlo method (Lee et al. [7]), response matrix-based method (Leppänen [8]).
This study relays on improving the initial guess of the fission source, in an attempt to preserve the other two components of the error in the fission source, i.e. without reducing the neutron histories per cycle. The baseline of the hybrid method proposed is the response matrix-based method developed by Leppänen, [8]. Through a short Monte Carlo simulation, the built-in response matrix-based solver produces a spatial distribution that approximates the fission source to be used by the source sampling as initial guess. The spatial source distribution is reconstructed on a cellwise form factor basis on a high-resolution mesh. A drawback of the method is that the accuracy of the spatial distribution of sampled source neutrons is limited by mesh resolution. In this work, an alternative to the cell-wise form factors is presented, which applies functional expansion tallies (FET) that are geometrically agnostic and continuous by nature.
The hybrid approach, FET-response matrix-based solver Monte Carlo method implemented in Serpent 2 Monte Carlo Code [9], is presented next, with a subsequent verification test case based on BEAVRS benchmark [10].

METHODS
The proposed methodology upgrades the built-in response matrix based solver when calculating an improved initial guess for the fission source distribution, by replacing the mesh-based approach (see Algorithm 1, [8]) used to calculate the cell-wise form factors that feed the response matrix solver, by a functional expansion reconstruction (see Algorithm 2) which through a high-order series represents the original and continuous spatial distributions of the form factors.
The reconstructed distribution for each cell n is calculated once the outer iteration cycle has been converged, by multiplying the source terms and inward currents with the corresponding form factors.  S Algorithm 3 has been formulated based on Leppänen [8,11]. Algorithm 4 has been constructed based on Wendt [12], Kerby [13], Cheatham [14] and Griesheimer [15]. for i ← in-boundary do 13: for j ← out-boundary do 14: if (J out m,j ≡ J in n,i ) then T j,i m,n = 1 15: Coupling coefficients: 21: for n ← mesh cells do 8: Collision cell-based Tallies (FETs) 9: for n ← neutron particles do

VERIFICATION
The hybrid approach, FET-response matrix-based solver Monte Carlo method, is implemented in Serpent version 2.1.32 (non-released version at the time of the writing). The calculation of intracell form factors relies on subroutines that were originally implemented to perform the coupling of Serpent to the MOOSE Framework [16], and later to be used with Serpent multi-physics interface [12]; while the built-in response matrix solver, which is being used for accelerating the fission source convergence [8], was initially designed as the built-in importance solver to produce weightwindow meshes for variance reduction [11].
In order to verify the implemented methodology, a single-assembly and full-core PWR test cases are proposed, based on the MIT BEAVRS benchmark [10], for which the Serpent input model had already been prepared for previous studies [17,18]. Fission source convergence is monitored using the Shannon entropy [19], for radial (x, y) and axial (z) dimensions. To reduce random oscillations in the entropy and assure that the results are not affected by the source bias [20], the simulations were run with an enough large population size (standard deviation of σ < 0.1%).

Single-assembly Calculation
The single-assembly test case is defined as a single 2.4% enriched fuel assembly without burnable absorber. The assembly is modelled in an 3D-infinite lattice, with reflective boundary conditions radially and vacuum boundary conditions axially. The fission source convergence is monitored only in the axial dimension. An axially asymmetric fission source distribution is obtained by inserting the control rods, 1/4 of the active fuel height (see Figure 1).
(a) Radial cross-section (b) Axial cross-section Figure 1: Radial (1a) and axial (1b )cross-sections of 2.4% wt. U-235 fuel assembly. Based on a previous study [8], the response matrix solution is obtained in a 1D-axial Cartesian mesh of 11 bins. The intra-cell source distribution reconstruction is performed in 40 bins for the mesh-based case and with 10-th order Legendre polynomial series expansion for the FET-based case (initial analysis phase). The Monte Carlo simulation samples 10 6 source neutrons per cycle.
The evolution of the fission source distribution is captured along three stages (see Figure 2): SIM1, first response-matrix solution initialises from an uniform source distribution, that converges after 79 outer iterations. Side-by-side, it can be observed the behaviour of, both, the mesh-based and FET-based approaches, still far from the converged solution. The functional expansion preserves the distribution continuity and provides a smoother solution. The converged solution is used as initial iteration for the second simulation, which produces the coupling coefficients for the next response matrix solution, SIM2. After 89 outer iterations, the source is converged, being relatively close to the reference solution. The FET-based case shows a good agreement with respect the converged Monte Carlo solution. At SIM3, the response matrix methodology no longer shows improvements, just random variation/noise, intrinsic to the stochastic sampling. This behaviour agrees with the root-mean-square error (RMSE): 0.2582 (SIM1), 0.1331 (SIM2) and 0.1190 (SIM3) -calculated between the reference solution and the FET-based approach.
The new reconstruction method evaluates the fission source in a more congruent and native form to Monte Carlo, given its geometry-free and continuous nature, allowing to improve the convergence performance. Although, the polynomial series estimate, due to its truncation (numerical) error, does not completely reproduce the heterogeneities as a result of the rods insertion and/or the gridspacers. Additional studies are necessary to establish the parameters that fit the requirements of the case, such i.e., the polynomial expansion order or polynomial basis function.  Figure 3 shows the evolution of the inward current response coupling coefficient and direct source to response coupling coefficient, for the two first simulations (SIM1 and SIM2).
The response matrix solver computes a population of 10 6 source neutrons 3-times; which given the size of the problem, it is negligible (each simulation takes less than a second in total). There is not a significant difference between the mesh-based and FET-based approach in terms of computational time. The source convergence after the criticality source cycles is presented in Figure 4. The response matrix methodology, for both approaches, calculates the coupling coefficients in agreement; with a better performance than the non-accelerated case since the beginning of the calculation.

Full-core Calculation
The full-core test case represents a critical core departing from the HZP start-up PWR core specified in the BEAVRS benchmark; where all control rods are fully withdrawn, except bank D with an insertion of ∼ 6.5% of the active fuel height, and coolant concentration of 975 ppm, (see Figure 5). Based on a previous study [8], the discretisation of the response matrix is 25 × 25 × 22 Cartesian mesh; where each bin radially represents an assembly. The intra-cell source distributions are reconstructed in a 4 × 4 × 4 bins for the mesh-based case and in a 10-th order Legendre polynomial series expansion, in each dimension, for the FET-based case. The Monte Carlo simulation samples a population of 10 7 source neutrons per cycle.  FET-based intra-cell form factors cases are presented; for which, the response matrix solutions takes 30 − 50 seconds in total. Radial entropies converge quite fast, due to a flat profile of the reactor core radial power. It is not the case of the axial entropy, which requires more time to achieve convergence, given a uniform initial source distribution. Nevertheless, the entropy evolution of the accelerated cases is faster than the non-accelerated case given an improved initial guess.

CONCLUSIONS AND DISCUSSION
A new methodology is presented to accelerate the convergence of the fission source on Monte Carlo calculations. The approach is embedded into the response matrix-based method; which, through a short Monte Carlo simulation provides an improved initial guess distribution. The idea behind this upgrade is to reconstruct the intra-cell source distributions by applying functional expansion tallies (FET) instead of a high-resolution mesh, providing a more realistic and accurate solution.
The new implementation reconstructs the fission source in a more natural way: due to its geometry agnostic and continuous nature, which allows to enhance convergence performance. However, the polynomial series estimate, due to its truncation (numerical) error, seems not to fully capture the heterogeneous behaviour caused by the rods insertion or the grid-spacers. Nonetheless, the method shows potential in achieving a nearly converged solution at the initial stages, as the test cases, single-assembly and full-core, show. In terms of performance, the functional expansion tallies as an alternative to a high-resolution mesh to calculate the coupling coefficients involved in the response matrix calculations, do not suppose a major reduction of the computational cost, due to the shortness of the response matrix-based method itself: reconstructing the intra-cell source distribution takes an almost negligible amount of time compared with whole simulation. As well, the response matrix methodology, for both mesh-based and FET-based approaches shows random variations given its stochastic sampling formulation.
The suitability of the functional expansion tallies to accelerate the convergence of the fission source, brings new and improved features to the response matrix method, in addition to, a potential new range of applications. Nevertheless, further testing and analysis are necessary in order to understand those capabilities and performance of the upgraded response matrix approach, due to showing these in practice is not straightforward.
As well, it seems important to provide guidance on the method's use, based on the effect of the freeparameters, e.g., neutron size batch, mesh discretisation, polynomial basis function and expansion order, etc; since, at the current stage, the presented approach is case-based: needing user's inside knowledge on the case itself.
A natural evolution of the method would be to evolve in parallel, both, the response matrix builtin solver and the functional expansion tallies, e.g. additional geometry mapping, methodology automatisation, etc. Also, the response matrix method and the functional expansion tallies could be applied to other problems: 1) independently, i.e. response matrix as a reduced order method apply to Monte Carlo calculations, or functional expansion tallies apply to reconstruct Monte Carlo spatial/temporal distributions; and 2) combined, such i.e., time-dependent Monte Carlo methods. Initial fission source rate s j Outward current through cell boundary j response coefficient T j,i m,n Current topology coefficient, from cell m boundary j to cell n boundary i Algorithm 4: 'Functional expansion tallies (FETs)' on page 4 ω n n th particle weight ψ m (x) m th basis function expansion ψ m (x n ) n th particle contribution to the m th polynomial ρ(x) Associated weighing function for the basis function Σ t (x n,k )

NOMENCLATURE
Macroscopic transport cross-section a m m th functional expansion coefficient F (x) Function to be expressed as a series expansion k m n n th , m th orthonormalisation constant transformation to Zernike's sub-space k m m th orthonormalisation constant transformation to Legendre's sub-space k m m th orthonormalization constant N Total number of particles P m (x) m th Legendre's polynomial expansion Z m n (r,φ) n th radial & m th azimuthal Zernike's polynomial expansion a m m th functional expansion estimate coefficient F c (Λ) Linear (3D) convolution of functional expansion series r,φ Radial & azimuthal coordinate transformation to Zernike's sub-spacẽ x Coordinate transformation to Legendre's sub-space