OPTIMISATION OF CRITICALITY AND BURNUP CALCULATIONS IN MONK®

The primary goal of this paper is to increase the efficiency of criticality and burnup calculations in the ANSWERS MONK® Monte Carlo code [1]. Two ways of achieving this goal are investigated as part of the H2020 McSAFE Project: creating a unified energy grid for all materials in the model, and reducing the spread in variances of fluxes for depletable materials using a generated optimised importance map. The average tracking speedup factor across all cycles of all burnup calculations ran using the unified energy grid, at base temperature, was found to be 1.96. For criticality calculations at 400K with runtime Doppler broadening, the unified grid approach gave a total speedup factor of 7.32. This demonstrates the potential importance of this method to reduce the calculation time with models with runtime Doppler broadening. The use of the generated optimised importance map has been demonstrated to significantly reduce the variance in the standard deviations on the fluxes in the fuel pins across two different test cases. If a solution is required in which the standard deviation in none of the fuel pins exceeds 5% it was found that the number of scoring stages required was more than halved, highlighting the potential for the outlined methodology to speedup burnup credit calculations.


INTRODUCTION
At the heart of Monte Carlo (MC) neutronics codes is the interrogation of nuclear data for the purposes of neutron tracking. Many modern MC codes (including the UK Monte Carlo code MONK® [1]) use a detailed nuclide dependent continuous energy representation of cross section data. The continuous energy collision processor of this form used in MONK® is called BINGO. Particle tracking requires repeated calculation of the mean free path of the particles, which is the inverse of the total macroscopic cross section. As the materials in a burnup calculation comprise many nuclides, which have different energy grids, the lookup process can be slow due to the large numbers of loops required over materials and nuclides. During the McSAFE project a unified energy grid for all materials and nuclides when searching for the required nuclear data was found to invoke a speedup to this process [2]. A suitable algorithm for formulating this unified energy grid is outlined in the methodology.
The use of importance maps to improve the statistical uncertainty in Monte Carlo calculations is routinely used in shielding application codes such as those performed with ANSWERS® [3] MCBEND shielding and dosimetry Monte Carlo code [4] which has the capability to automatically generate an importance map via an approximate adjoint solution for splitting and Russian roulette without bias [5]. In criticality codes such as MONK ® [1], the use of importance maps to improve the stochastic uncertainty of the flux is less common. This may be because importance maps bias the code to track more neutrons in some regions, which inevitably results in a reduction in neutrons in the other regions. Consequently, the standard deviation of the estimated flux may be lowered in some regions but only at the expense of an increased standard deviation elsewhere in the model. A McSAFE project aim is to reduce uncertainty for depletion calculations where the aim is to determine the composition of fuel elements at their end of life. In order to determine the composition of all pins to a required accuracy the flux history of the pins is required to a given accuracy. In a standard reactor physics Monte Carlo calculation, the neutrons will be targeted preferentially at the more reactive central regions of the core. Thus, the statistics will tend to be better towards the centre of the core and worse towards the outer edges of the core. Hence, the flux and consequently the composition of the outer fuel pins will be determined less accurately. In this case an importance map is proposed that will flatten the profile of the standard deviations of the fluxes for fuel pins across the core. Hence, the calculation will reach the required accuracy for all fuel pins within a reduced calculation time.

METHODOLOGY 2.1. Unified Energy Grid
Tracking requires repeated calculation of the mean free path (mfp) length which is the inverse of the total macroscopic cross section. Evaluation of the total cross section requires looping over materials, nuclides and energies in order to find the correct energy values on the grid to use for interpolation. This is a very expensive process computationally. The aim of the work reported here is to pre-calculate total cross sections for each material on a suitable unified energy grid, thereby avoiding the looping over materials and nuclides and thus speeding up the calculation of the mfp. The aim is to combine the energy grids for individual nuclides into a single grid. In order to avoid using large amounts of memory to store the unified energy grid a method is proposed to reduce the number of energies stored without compromising resolution. The process proceeds as follows [2]: 2.1.1 Energy Grid Unification Algorithm 1) Construct an energy grid that contains all of the energies from the grids for all of the nuclides and materials 2) Construct a grid containing all of the important energies, which are defined as: i. energies corresponding to local minima or maxima of the total cross sections; and ii. energies corresponding to the onset of threshold reactions 3) Reduce number of energies in base grid by thinning: i. Choose a value for the thinning parameter, ߬, (user-defined in the MONK implementation) ii.
Test: iii. If then replace ‫ܧ‬i and ‫ܧ‬i+1 by one energy: Repeat until above criterion not met by any adjacent energies on the grid. 4) Insert the important energy points into the thinned grid.
Once the unified energy grid has been constructed the Doppler broadened cross sections are calculated and stored at the appropriate temperatures for each material at each energy point on the unified grid. During tracking the total cross section is calculated by interpolating between the stored values on the unified energy grid. Calculations presented for the Unified Energy Grid use the ratio between the real number of fission children produced during a superhistory against the number of neutrons tracked during that superhistory as its sole estimate of keff. Other estimators which use the expected number of fission children are not used, which decreases runtime at the cost of increased uncertainty in keff for a given number of samples.

Optimised Importance Map Generation Using an Inverse Adjoint Flux
Consider a small detector placed in a Monte Carlo simulation. The adjoint flux can be shown to represent the importance of particles produced at a specific location to the response of the detector [5]. For a criticality calculation each fissile region in the core can be considered both as a detector and a source, as each point contributes to the reactivity of the core. Thus the importance of each point to the overall reactivity of the core is obtained by integrating over all detectors, i.e. integrating over the whole core, to obtain the adjoint flux. Monte Carlo criticality (and reactor physics) codes are designed to produce the best estimate of core reactivity or keff and hence target particles preferentially towards the more reactive parts of the core where more fissions occur. As noted in the introduction, this results in lower variances in the estimated fluxes in the more reactive regions of the core and higher variances in regions of lower reactivity. In order to obtain a more uniform distribution of variance throughout the core it is necessary to remove this tendency to target the more reactive regions of the core. As the use of the adjoint flux as an importance map will create a bias towards the more reactive regions of the core the use of the inverse of the adjoint will negate this process. Therefore, the inverse of the adjoint flux is used as the starting point to create an importance map. It is preferable to discretise importance maps in criticality and burnup calculations into powers of two. This will ensure that the weights of the particles passing through each mesh cell are the same. Stopping a spread of weights in a mesh is important when estimating a flux because otherwise it significantly increases the variance in that mesh. Algorithms for splitting with non-discretised importance maps without bias do exist. However, as importance maps are generally approximate it is not worth the calculation effort to split on the small ratios of importance; which occur regularly in these schemes [5]. The algorithm proceeds as follows:

Importance Map Generation Algorithm
1) Normalise the Inverse Adjoint Flux (IAF) using the smallest value in the IAF 2) Find the closest 2 x value (up to 32 to avoid over splitting) to every value of the normalised IAF, and set that in a different array to be the current discretised guess of the Importance Map. 3) For all terms calculate the current normalised IAF divided by the current guess of the importance map. 4) Add all these terms together, normalising by the number of non-zero terms. 5) Check whether this summation agrees with the summation of the previous iteration to within a given tolerance. For the test calculations in this study a tolerance of 10 -5 was found to work well. If true then EXIT. 6) If false, to optimise the previous normalisation of the IAF, divide the currently normalised IAF by the summation. Go to step 2. If the system is symmetrical, an option has been included to force the importance map to be symmetrical.

TEST CASES 9x9 and 17x17 Fuel Assembly (Unified Energy Grid Calculations)
The geometry and configuration for the 9x9 and 17x17 fuel assembly models are shown in Figure 1. The fuel assembly has three fuel materials and reflective boundary conditions in the x and y directions. The calculations use a pin-wise BU (burnup) mesh with two axial layers (lower and upper).

Small PWR (Unified Energy Grid and Importance Map Calculations)
The geometry and configuration for the PWR core are shown in Figure 1. The fuel assembly has a single fuel material and reflective boundary conditions in the x and y directions. The calculations use an assembly-wise BU mesh with eleven (Unified Energy Grid) and ten (Importance Map) axial layers. No top and bottom reflectors are present in the outer axial layers of the model.

Minicore (Unified Energy Grid Calculations)
The geometry and configuration for the Mini core are shown in Figure 1. The core is placed in a container of water with free outer boundary conditions. The calculations use a BU mesh with six layers in each of the x, y and z directions.

Nakagawa and Mori Whole Core Benchmark PWR (Importance Map Calculations) [7]
The Nakagawa and Mori PWR whole core benchmark specification is given in [7]. The layout of the core is shown in Figure 1.The flux is estimated on a 60x60x1 (xyz) mesh (16 meshes per FA) over the extent of the fuel assemblies in x and y and the active length of the fuel in z. Vacuum boundary conditions are applied on all outer sides of the benchmark.

Unified Energy Grid
Calculations vary the thinning parameter to determine its effect on the results. Values between 0 and 10 -4 are used. The speedup is determined by dividing the total run time of a version of the code with the modifications outlined in the methodology against a standard version of MONK ® . Burnup calculations were run over three fuel cycles with flux tallied in Unified Tally within MONK ® . Criticality calculations were run using the same Unified Tally mesh as the burnup calculations with fresh fuel. All results displayed have a keff which is within three standard deviations of the reference calculation. Results for the criticality calculations are displayed in Figure 2. This shows roughly a 10 to 20% speedup due to turning off some of the criticality estimators and up to a factor of 2.5 speedup due to the use of a unified energy grid. However, the standard deviation on keff using the standard code is roughly a factor of two smaller and so there is no overall benefit to the criticality calculations. The speedup factors for the burnup calculations are given in Figure 3. This shows up to a factor of 1.6 speedup when turning off some of the criticality estimators and up to a factor of 2.25 speedup from using the unified energy grid for the total cross section lookup. The speedup is seen to be relatively insensitive to the value chosen for the thinning parameter when only tracking is considered. However, Figure 3 shows a greater variation on total computing time with thinning parameter as the use of a non-thinned grid can significantly increase the setup time for a calculation. Figure 4 shows the potential benefits of using the unified energy grid in criticality calculations which utilise run time Doppler broadening (RDB) [8]. Although there is little effect on burnup calculations, a significant speedup (a factor of around seven) was seen in criticality tests. This speedup is high enough to provide an overall speedup when the need for additional histories from disabling estimators is considered.

Optimised Importance Map Generation Using an Inverse Adjoint Flux
All results presented using this algorithm have a keff within three standard deviations of a reference solution. The importance map generated (Fig.6) from the adjoint flux (Fig.5) increases from the centre to the outer boundaries. This will compensate for the lack of samples that would otherwise be obtained near the axial boundaries due to the decrease in flux adjacent to these boundaries. The standard deviations in the standard MONK calculation for the top and bottom axial layers are seen to be roughly double those in the centre of the core. The standard deviations when the importance map is used are seen to be much more uniform (Fig. 7).

Figure 7. Comparison of Standard Deviations for all fissile regions; Standard MONK (top), with Importance Map (bottom)
The use of the importance map is seen to produce an 11% reduction in the RMS and a 76% reduction in the spread of the standard deviations as given below in Table I. A second set of calculations was performed in which the calculation was stopped once all of the standard deviations were below 5%. This kind of calculation is useful if the analyst wants to achieve a minimum accuracy on the flux history in all fuel pins in order to guarantee a required accuracy for fuel depletion in the reactor. The results are displayed in Table II. This shows that fewer than half the number of scoring stages is required for this purpose when the importance map is used.  Figure 8. This indicates that the most important region for core reactivity is the centre of the core resulting in poor statistics for fuel at the edge of the core. The standard deviations with and without the use of an importance map are displayed in Figure 9 which shows that the high uncertainty at the edge of the core in the standard MONK calculation is greatly reduced when the importance map is used to target samples at this region. This also shows that the significant reductions in standard deviation at the edge of the core is balanced by small increases in standard deviation in the centre (Fig. 9 Absolute Difference).

Figure 9. Comparison of Standard Deviations in the Benchmark Energy Group Scheme Standard MONK (top), with Importance Map (middle), Absolute Difference (bottom)
The effect of using the importance map on the uncertainty (standard deviation) is also shown in Table III in the four energy group scheme defined in the Nakagawa and Mori benchmark. The results are reasonably insensitive to energy group and show that the use of the importance map reduces the RMS of the standard deviation by several percent and reduces the spread in uncertainties by over one third.

CONCLUSIONS
Within the McSAFE project a number of modifications have been made to a development version of the MONK criticality and reactor physics code to improve the efficiency of burnup calculations. A unified energy grid has been implemented for use in estimating the mean free path between collisions to speed up calculations. To quantify the resulting improvement in performance burnup and criticality calculations were performed for four models. A speedup by a factor of around two was obtained for burnup calculations but there was no benefit for the criticality calculations unless the computationally more intensive run time Doppler broadening was used to model components at the correct temperatures. Turning off keff estimators using expected values when using the unified energy grid can further reduce the run time. This was found to speedup calculations by up to 60% in some cases.
A procedure has been developed to generate an importance map from a deterministic adjoint flux solution for use with Monte Carlo reactor physics calculations. The aim is to reduce the variance in the standard deviation on the fluxes in the fuel pins across a reactor core, by the use of particle splitting and Russian rouletting to target more samples at regions with poor statistics. Two models have been studied to evaluate the effect of the proposed procedure on the standard deviation of the flux. The generated importance map was found to significantly reduce the variance in the standard deviations on the fluxes in the fuel pins across the core. If a solution is required in which the standard deviation in none of the fuel pins exceeds 5% it was found that the number of scoring stages required was more than halved if the importance map was used for variance reduction.