OPTIMISATION OF MONTE CARLO BURNUP SIMULATIONS

We show here that computing efficiency of Monte Carlo burnup simulations depends on chosen values of certain free parameters, such as the length of the time steps and the number of neutron histories simulated at each Monte Carlo criticality run. The efficiency can thus be improved by optimising these parameters. We have set up a simple numerical model that made it possible for us to test a large number of combinations of the free parameters, and suggest a way to optimise their selection.


INTRODUCTION
It has become a common practice for users of Monte Carlo burnup codes to use relatively large time steps of several weeks or months because this helps to reduce the number of executions of the Monte Carlo criticality and depletion solvers. On the other hand, the changes in the neutron energy spectrum and the spatial shape of the neutron flux may not be correctly captured then, depending on the selected coupling scheme [1,2,3,4,5,6,7]. The choice of the length of the time step thus represents an optimisation problem. We study this problem here and suggest a way to optimise the user-defined free parameters.
A user choice of the number of neutron histories simulated at each Monte Carlo criticality calculation is, in fact, related to the optimisation of the length of the time steps when the total computing cost of the whole Monte Carlo burnup simulation is assumed to be given. For instance, reducing the number of neutrons per criticality run makes it possible to increase the number of time steps. In such a case, neutron flux will be computed with larger random (statistical) error; however, the random error will not propagate as much into the fuel depletion thanks to the short time steps. The effect of random errors will cancel out better over a large number of short time steps [8].
We have previously studied this optimisation problem. In [9], we have set up a very simplified solver, also briefly described here, in Sec. 2, and performed a large number of test simulations with varied values of the free parameters. The test results made it possible for us to suggest an optimisation algorithm of the free parameters [9]. In this paper, we show results not reported previously (see Sec. 3), and we newly demonstrate the optimal time step length at which the computing efficiency is maximised. We discuss and summarise our conclusions in Sec. 4.

NUMERICAL MODEL
To be able to carry out a large number of simulations, we resorted to a very fast and very simplified numerical model which, however, still possesses the main features of Monte Carlo burnup simulations. We derive the neutron flux solution in each criticality cycle from the one-group diffusion equation. We then introduce a random noise in the solution by statistical sampling of the neutron flux at each criticality cycle from the diffusion solution. We use the inverse transform method for the statistical sampling of the noisy solution of the required neutron sample size. The depletion solver is also simplified. We solve the burnup equation only for the following selected nuclides: 235 U, 238 U, 239 Pu, 135 I and 135 Xe. The depletion solver computes the concentration 135 Xe at its saturated level. The saturated xenon concentration denotes a concentration of 135 Xe that establishes itself at an infinite time in a system with a specific fixed neutron flux.
The model geometry represents a 350 cm wide slab (1D model) with black boundary conditions. The slab is split equally into 1000 zones filled with independently defined materials. The numerical model depletes the materials independently according to locally computed neutron fluxes. Initial mass densities of the nuclides present in the slab material are given in Tbl. 1 along with the group microscopic cross sections and ν (the number of fission neutrons per fission) for the selected nuclides. The group data correspond to a thermal neutron energy spectrum. The diffusion solver can also set up an equilibrium concentration of 135 Xe during criticality cycles. The equilibrium xenon concentration represents a concentration of 135 Xe that is in a steady-state with the neutron flux (not fixed) in a critical system. The equilibrium xenon concentration is achieved by setting a saturated xenon concentration at each criticality cycle according to the neutron flux distribution tallied at the previous criticality cycle. This procedure is essentially the same as the one implemented in the Serpent code [10].
Several coupling schemes were implemented in the numerical model and used in tests: the beginningof-step constant flux approximation scheme (based on the explicit Euler method), the predictorcorrector method [11,4] and the stochastic implicit Euler based scheme [5]. In the implicit scheme, materials are depleted over each time step with the neutron flux computed at the end of the time step. This scheme has more variants, but all of them introduce a number of "inner iterations" in each time step. In the variant implemented here, materials are re-depleted with the flux averaged over all previous inner iterations. This process continues until a user-specified number of "inner iterations" is reached. Then the simulation proceeds to the next time step.

NUMERICAL TESTS
In this section, we evaluate the computing efficiency in terms of the figure of merit, F OM, where c Δt is the real wall-clock computing time of the whole simulation, and ε is the error of the Monte Carlo burnup simulation. In this study, we choose to associate ε to the error in the actual material depletion process, and not to errors in the neutron fluxes at specific time steps. For this reason, we choose to evaluate ε as the error in the depletion of 239 Pu, averaged over all (1000) nodes, and also averaged over the time steps. The error ε is therefore calculated for each test simulation as where: I is the total number of time steps, Z t i gives the atomic density of 239 Pu in all materials at the end of the ith time step t i (a specific element in the Z vector returns the plutonium density in a certain material), and The reference solution, Z t i ,ref , was obtained from a simulation in which the neutron flux in criticality cycles was not sampled statistically, and the direct diffusion solution was used. The neutron flux in the reference simulation, therefore, contained no random noise. The reference solution was obtained from a simulation based on the Explicit Euler scheme with a very short time step length of 1h, and equilibrium xenon concentration was kept being iterated during criticality calculations in all time steps (ensuring a numerically stable simulation).
In each test simulation in this section, no inactive cycles have been set, and the source iterated at a certain time step was always re-used as the initial source in the successive time step. The initial source in the first time step was always sampled from a uniform distribution. Next, the total period was set to 4.2 years, and 1000 neutrons were sampled at each criticality cycle. Also, the equilibrium xenon concentration was set on in each simulation.
Test simulations were done with various coupling schemes: simulations with the beginning-ofstep constant flux approximation scheme are denoted as "EE", simulations based on the predictorcorrector scheme are denoted as "PC", and simulations with the stochastic implicit Euler based scheme are denoted as "SIE". The test simulations were also re-made for a number of time While the simulations differ in the coupling schemes they use or in their time step lengths, each test simulation sampled the same total number, N n = 10 8 , of neutrons, to achieve the same statistics. This was ensured by setting a unique number of criticality cycles in criticality calculations in each simulation. So, the product of the number of criticality calculations (executions of Monte Carlo solver) and the number of criticality cycles was the same for every simulation.
To evaluate the figure of merit F OM of various test simulations, the total cost needs to be known. Nevertheless, the real computational cost of the simple depletion and criticality solver are not typical here due to their simplicity. Therefore, for evaluation of F OM, we choose to set the cost of the simulation of a single neutron history to 10 -4 s (needed in criticality calculations) and we vary the cost of the depletion solver over few values: 0, 5 and 20 s.

Test #1: Time Step Length vs. Computing Efficiency
In Test #1, we wish to demonstrate as to how the computing efficiency may vary with the selected length of the time steps. Fig. 1a and Fig. 1b Fig. 1a and Fig. 1b, we report the number of criticality cycles each SIE simulation used. The SIE simulations could not cover all the selected time step lengths since they had to carry out at least one inner iteration with the selected statistics.
The F OM in Fig. 1a was evaluated with a zero cost of depletion calculations. This plot serves only a demonstratory purpose in showing that if the depletion solver cost was zero then there would not be any strong reason for not setting very short time steps. Fig. 1a shows that especially the EE and PC based simulations benefit from using very short time steps. The SIE simulations appear not to be as sensitive to the choice of the time step length. Nevertheless, Fig. 1a shows that SIE simulations with a smaller number of criticality cycles (hence the larger number of inner iterations) achieve better efficiency. This result appears general for all simulations with a zero depletion cost irrespective of the coupling scheme-simulations become more efficient as the number of criticality cycles in criticality calculations is reduced. (EE and PC based simulations with shorter time steps use a smaller number of criticality cycles than simulations with larger time steps.) Typically, the number of executions of the Monte Carlo solver would match the number of executions of the depletion solver in the burnup simulations. When an optimal number of criticality cycles (i.e. an optimal cost of the Monte Carlo solver) exists for a given cost of the depletion solver, as the above results suggest, then it may be interesting to study as to how the F OM of burnup sim-

Test #2: Optimal Conditions for Monte Carlo Burnup Simulations
Test #1 suggested that the efficiency of Monte Carlo burnup simulations may depend on the proportion of the costs of the depletion and Monte Carlo solvers in the overall simulation cost. For this purpose, we define a new variable x here as the ratio of the cost of all depletion calculations (and other procedures that are not involved in neutron transport simulations) to the overall cost of the whole Monte Carlo burnup simulation. Note that simulations with shorter time steps have larger x because they execute the depletion solver more times. Fig. 2a and Fig. 2b show the F OM of test simulations with respect to x. F OM in Fig. 2a were computed for the depletion calculation cost of 20 s, and F OM in Fig. 2b were computed for the depletion cost of 5 s. Results suggest that the F OM has a clear dependence on the value of x, and this dependence is similar irrespective of the coupling scheme that the simulations are based on. We have observed that the F OM vs x curves from all our simulations are all very similar and overlap. To keep the curves separated in the figures, we included SIE simulations with only one time step length of 32 d (the x value is varied by changing the number of inner iterations). Nevertheless, we have observed the same results for SIE simulations with other time steps as well. Fig. 2a and Fig. 2b show that, irrespective of the coupling scheme and the depletion solver cost, all test simulations achieved the best efficiency when they had x inside the interval (0.2, 0.8). All simulations with x close to 0 or 1 underperformed.

DISCUSSION AND CONCLUSIONS
This work is an extension of our previous study [9,8]. Test results suggest that efficiency of Monte Carlo burnup simulations may be maximised when the fraction of the cost of all depletion calculations (and all other procedures not directly involved in the simulation of neutron histories) in the overall cost of the whole Monte Carlo burnup simulation, x, is inside a certain interval that we approximately evaluated (using a simple numerical model) to be (0.2, 0.8). We wish to note that the numerical model used in this study was rather very simple, and the range of optimal values of x may be somewhat different in real Monte Carlo burnup simulations. This will be a subject in our next study.
The above result can be directly used for the optimisation of free parameters of Monte Carlo burnup simulations such as the time step length or the number of neutron histories simulated per time step. This optimisation is possible since the computing cost of various computing procedures can be measured at the beginning of the Monte Carlo burnup simulation. Knowing the cost of a single neutron history and the optimal fraction of the cost of transport simulation in the overall allocated cost directly gives the optimal total number of neutron histories to be simulated during the whole simulation, and a similar process can give the optimal number of depletion solver runs and hence the number and length of the time steps.
While our method can suggest an optimal number of neutron histories per time step, the user is still free to decide as to how these histories would be distributed with respect to the number of cycles and the neutron population per cycle at each Monte Carlo run. Clearly, these two parameters are not independent when the overall statistics is given by the presented method. The user choice of a very large number of neutrons per cycle will then necessarily limit the number of criticality cycles which may impede the initial source convergence and the overall performance of the simulation.
On the other hand, the choice of a very small number of neutrons per cycle may amplify other source convergence problems. Although our tests demonstrate that a similar efficiency can be achieved from simulations based on various coupling schemes, certain differences may be envisioned. For instance, unlike the simulations based on EE or PC schemes, the SIE based simulations can achieve good efficiency even with large time steps. This can be preferable when the neutron flux is needed among other results at each time step. This is because the SIE simulations with larger time steps combine the neutron flux over several inner iterations, which helps to reduce the random noise in the flux.