Forced propagation method for Monte Carlo fission source convergence acceleration in the RMC

Ze-Guang Li· Kan Wang· Yu-Chuan Guo · Xiao-Yu Guo

Abstract In loosely coupled or large-scale problems with high dominance ratios, slow fission source convergence can take extremely long time,reducing Monte Carlo(MC)criticality calculation efficiency. Although various acceleration methods have been developed, some methods cannot reduce convergence times, whereas others have been limited to specific problem geometries.In this study,a new fission source convergence acceleration (FSCA)method, the forced propagation (FP) method, has been proposed, which forces the fission source to propagate and accelerate fission source convergence. Additionally, some stabilization techniques have been designed to render the method more practical. The resulting stabilized method was then successfully implemented in the MC transport code,and its feasibility and effectiveness were tested using the modified OECD/NEA, one-dimensional slab benchmark, and the Hoogenboom full-core problem. The comparison results showed that the FP method was able to achieve efficient FSCA.

Keywordst Fission source convergence acceleration ·Monte Carlo method · Forced propagation method · RMC code

1 Introduction

The Monte Carlo(MC)method is widely used in reactor physics analysis because of its high fidelity, and with the advancement of computer techniques, it has been used increasingly in full-core analyses. When performing MC criticality calculations, the inactive cycles should be carried out first to obtain the converged fission source. The calculation of inactive cycles does not contribute directly to the final results,such as the effective multiplication factor,keff, and flux; however, it consumes a sample computing time,especially for loosely coupled or large-scale problems with high dominance ratios [1, 2], such as full-core problems.Therefore,accelerating fission source convergence is important for the efficiency of the MC criticality calculation.

Various fission source convergence acceleration(FSCA)methods have been proposed for MC criticality calculations, such as the super-history method [3, 4], stratified source sampling method [5], Wielandt method [6], functional MC(FMC)method[7],coarse-mesh finite difference(CMFD) acceleration method [8, 9], asymtotic Wielandt method (AWM), asymptotic super-history method (ASM)[10], and others [11–13]. Among these, the super-history method and Wielandt method are based on similar strategies, which only involve simple modifications to the standard power iteration scheme and are easy to implement. Previous studies have indicated that although these two methods were able to decrease the number of inactive cycles needed to achieve fission source convergence, the calculation time during each inactive cycle increased, and the total convergence time was not reduced compared to the standard power iteration [3, 4, 6, 14].

Some other acceleration methods imposed certain restrictions that limited the advantages of the MC method.For example, the CMFD method is required to use multigroup cross-sections in MC calculations, and the FMC method is restricted to a regular geometrical mesh. Some other methods need to reduce the neutron numbers in each cycle to reduce the total convergence time and could introduce fluctuations and destabilization into the calculations.

To reduce the actual total convergence time and improve the versatility of the MC FSCA method, a new method named forced propagation (FP) method has been proposed herein and implemented into the MC transport code RMC. As an alternative to the propagation process introduced by standard power iteration to accelerate fission source convergence,the FP method estimates the changing trends of fission source distributions in MC criticality calculations and forces the fission source to propagate instead of the propagation process introduced by standard power iteration to accelerate the fission source. This proposed method can be implemented using the standard MC transport calculation in the RMC [15], can retain MC method versatility,and can reduce both the convergence of inactive cycles and calculation times.

In this paper, a general introduction to the proposed FP method,including the techniques used in its stabilization,is provided in Sect. 2. Section 3 contains a detailed description and the results of numerical tests applied to validate the new method, and our conclusions are presented in Sect. 4.

2 General introduction to the forced propagation method

In conventional MC criticality calculations, the fission sources are propagated using standard power iterations.When the calculation problem has a dominance ratio (DR)close to unity [1, 2], such as a loosely coupled case with few neutron interactions between fissile media or largescale cases such as an actual core, the fission sources are difficult to propagate and slow to converge.The FP method estimates the changing trends of fission source distributions in the standard MC power iteration cycle and uses them to adjust the fission source distributions, thereby forcing the fission sources to propagate. By forcing fission source propagation, the method can achieve higher fission source convergence efficiency and accelerate the convergence process.

In this section,we first briefly introduce the FP method,and then describe the techniques used to stabilize it and make it more practical.

2.1 Description of the forced propagation method

The neutron transport equation of criticality is calculated using the equations listed below:

In the MC method,the neutron transport equation in(1)is solved by power iteration, which can be represented as shown in Eq. (2):where f（j）indicates the initial fission source distribution of the current cycle(j),and k（j）represents the current estimate of the effective multiplication factor. In conventional MC power iterations, the source neutrons from f（j）of the current cycle are tracked randomly,according to the transport process, to obtain the flux distribution, ψ（j）. Then, the fission source distribution f（j+1）used for the next cycle,(j + 1), can be simulated, based on ψ（j）, along with the transport process,using the fission operator,F.The process from the initial fission source distribution, f（1）， to the converged fission source distribution, f（n）, as shown in Fig. 1, could be referred to as fission source propagation,and the operator P=F（T -S）-1could be treated as the propagation operator.

In the standard MC power iteration, the initial fission source distribution changes continuously to get the converged fission source distribution. This means that by estimating the changing trend of the fission source distribution and multiplying a factor operator with the propagation operator,P，as shown in Fig. 2,we could accelerate fission source convergence.This approach has been named the FP acceleration method.

According to the FP process description above, the acceleration of fission source convergence can be presented as follows:

Fig. 1 Schematic of fission source propagation

Fig. 2 Standard power iteration and forced propagation processes

In the FP process, the operator（r ） may have different expressions,as long as they can follow the changing trend of fission source distribution and can accelerate fission source convergence. In this study, we proposed a simple way to obtain a proper FP operator,（r ） and easily achieved the FP method in the RMC code, as follows:

1. Divide the geometry into spatial meshes, M. As it is hard to estimate the continuous form of（r ） in the MC method, we need to establish a set of spatial meshes,and then estimate the discrete value of（r ）for each mesh,m,as.The spatial meshes in the FP method need not be specifically divided as in the CMFD acceleration method, but rather simply need to cover all the fissionable areas. Practically, we suggest that the meshes could be selected as the same meshes to tally the Shannon entropy of the problem [16].

4. Store the standard fission source distribution, fj+1（）,and according to the strategies in the RMC, all the neutrons in fj+1（）have the same weight, w(+1).

5. Tally the normalized spatial fractions of standard fission sources in the divided meshes. These normalized spatial fraction values can be calculated for each mesh as shown in Eq. (5):

6. Calculate the discrete fpmusing frm. As we obtain the spatial fractions of the standard fission source distribution in mesh m for cycles (j) and (j - 1), asand,thencan be calculated as shown in(6):

and theoretically, fpmwill converge to 1.0 when the fission source distribution has converged, which ensures that the FP method acceleration process is unbiased.

From this description of the FP method and the comparison of standard MC iteration and the FP method in Fig. 3, we could notice that unlike the super-history method or Wielandt method, the standard MC transport calculation does not need to be adjusted in the proposed FP method.Evidently,unlike the CMFD acceleration method,the proposed FP method does not need specially designed geometries or meshes to fit the requirements of deterministic calculations. Moreover, unlike the AWM and ASM methods,the neutron numbers in the inactive cycles do not need to be adjusted to reduce convergence time. Thus, it can be easily implemented and can solve arbitrary fission source acceleration problems.

2.2 Techniques to stabilize FP method

The description in Sect. 2.1 shows that the FP method could theoretically accelerate the fission source convergence process. However, when the results of fpm(which represents the FP operator) are not properly estimated, the method may not achieve the expected acceleration effect.Using practical problems as examples, if the frmin the previous cycle was very small, such as in the first few cycles with the initial point fission source distribution, it would likely lead to a very large result for fpmin the current cycle. This means that the changing trend of the fission source would be overestimated and may cause fierce oscillations in the fission source convergence process or even divergence. Alternatively,due to MC fluctuation,fpmmay have certain results which underestimate the changing trend of fission source, which would cause slow fission source convergence.

Several techniques have been implemented in the RMC code to prevent fierce oscillations or slow convergence when using the FP acceleration method, including the relaxation factor, upper and lower limits, and using momentum consideration to estimate the FP operator, as explained below:

1. Relaxation factor technique. To stabilize the FP method acceleration process, when estimating fpm, a relaxation factor,r,is implemented after the fpmresults are calculated by Eq. (6), as shown in Eq. (7):

Fig. 3 Comparison between the standard iteration and FP iteration processes

where r indicates the user-setting relaxation factor,between 0.0 and 1.0; with a smaller r, more stable but slower convergence could be achieved.Our testing has indicated that selecting a value for r between 0.6 and 0.8 will help achieve a proper balance between stabilization and acceleration effects.

2. Upper and lower limits of fpm. Normally, the relaxation factor technique can stabilize the FP method acceleration process. However, in cases with initial point fission source distribution, in the very first few cycles, the fpmprocess could easily derive very large results, and the relaxation factor could not perform quite well. Thusly, the upper and lower limits of fpmare set to prevent oscillation in these cases, and our experience suggests that the upper and lower limits of 2.0 and 0.5, respectively, would suit most cases.

3. Momentum consideration. The momentum technique has been widely and successfully used to accelerate and stabilize gradient descent iterations in machine learning and deep learning applications [17]. In this technique, information (such as the gradient) in the iteration histories is considered to adjust the current iteration, thereby preventing false and slow convergences, as shown in Fig. 4. In the RMC, we also used momentum consideration in the FP method, whereby fpmwas not only determined by the calculation of current cycle, but also by considering the momentum effect, mt, of historic cycles, as shown in Eq. (8):

Fig. 4 Effect of the momentum technique in gradient descent

where t indicates the momentum rate set between 0.0 and 1.0, and according to Eq. (8), fpmcould still converge to 1.0 to ensure that the FP acceleration process is unbiased. Considering the momentum effect properly, the underestimation of fpmcaused by MC fluctuations could be prevented. Our experience suggests that a momentum rate,t,of 0.25 would suit most cases.

2.3 Comparison between the FP method with other acceleration methods

The proposed FP method can estimate changing trends of fission source distributions in MC criticality calculations and forces fission sources to propagate,instead of using the propagation process introduced by standard power iteration to accelerate fission source convergence. By introducing stabilization techniques, including the relaxation factor,upper and lower limits, and momentum consideration, the FP method can also achieve a good balance between acceleration and stabilization.The proposed FP method can be implemented using the standard MC transport calculation in the RMC, can retain the versatility of MC method,and can reduce both the convergence of inactive cycles and time.

There are basically two different ways to accelerate fission source convergence in MC calculations. One is to change the standard power iteration scheme (change the form of transport equation), such as the Wielandt method or super-history method, by reducing the DR of the transport equation in loosely coupled or large-scale problems.The other way is to estimate (guess) the fission source distribution,which would be close to the converged fission source distribution from the results of previous cycles and allows the acceleration effect to be achieved,similar to the CMFD method and FP methods proposed in this paper.

It has already been proved that the methods that change the iteration scheme to reduce the DR cannot theoretically reduce the fission source convergence time. This was why the AWM and ASM were developed by reducing neutron numbers in inactive cycles to reduce the fission source convergence cycle and time. But the reducing neutron numbers in inactive cycles could introduce fluctuations in the calculation.

The methods which use calculation information to predict fission source changing trends or converged fission source distributions could be effective in accelerating fission source convergence. The CMFD method or similar methods,such as SP3,used in OpenMC or other MC codes,can achieve very fast fission source convergence,but these methods require regular geometry,which can be calculated by deterministic methods, as well as the capability of calculating multi-group cross-sections. Other methods, such as that reported in Toth and Griesheimer [13], applied a similar idea to change the fission source using results from the previous cycles.However,the FP method can achieve a more effective and stabilized acceleration process of fission source convergence.

3 Validation and results

The FP method and stabilization techniques described in Sect. 2 have been implemented in the RMC code, which was developed for MC reactor physics analysis by the Department of Engineering Physics, Tsinghua University,China. Two kinds of test cases are chosen in this section,including the one-dimensional slab benchmark and a Pressurized Water Reactor (PWR) full-core problem. All calculations reported in this section were performed on Linux platform, with 2 Intel Xeon V5-2690 v3 CPUs, and 256 GB memory. 30 CPU cores were used in parallel for all the following calculations. In these tests, fission source convergence performances were compared, in terms of converging cycles and convergence times, using Shannon entropies or other parameters.

3.1 One-dimensional slab benchmark

The first test problem was a modified OECD/NEA source convergence benchmark [18], in which two separated, 20-cm-thick (regions 1 and 3), infinite fissionable slabs were decoupled by an intervening infinite water slab(30-cm-thick;region 2),as shown in Fig. 5.The water slab in the modified problem was 10 cm thicker than in the original problem, which increased the fission source convergence difficulty, rendering the problem more challenging. The problem had a free boundary condition for the fissionable slab outer boundaries and had the initial fission sources at the center of region 1.In the following analysis,neutron batch sizes of 500,000 neutrons per cycle were used in the RMC calculations and were applied through 2000 cycles (1500 inactive cycles).

To show fission source convergence performance, the fission source fraction in region 1 and the Shannon entropy of the whole system have been discussed for comparison.As the system is symmetric, the fission source fraction in region 1 after source convergence was theoretically 0.5.The Shannon entropy of the system was tallied using evenly divided meshes with 5-cm-thickness, which meant that a total of 14 meshes were used for the Shannon entropy tally.

Fig.5 Geometry of the modified OECD/NEA one-dimensional slab benchmark

The standard power iteration was performed for the comparison to obtain reference results for the fission source fraction and Shannon entropy, and the FP method was performed to demonstrate its acceleration effects. For the FP method, the same meshes for tallying the Shannon entropy were used; relaxation factors of 0.6 and 0.8 were set for comparison, and the upper and lower limits and momentum rate were set at the values suggested in Sect. 2.

The evolution of the fission source fractions in region 1 and Shannon entropies along with calculation cycles using the standard power iteration as well as the FP method results achieved using different relaxation factors, can be observed in Figs. 6 and 7. Notably, the standard power iteration needed ～1000 inactive cycles to achieve the fission source convergence for the modified problem with a 30-cm-thick water slab. This was approximately twice as many cycles than was required by the original benchmark using a 20-cm-thick water slab [10]. In contrast, when the relaxation factor, r, was set as 0.8 (0.6), the FP method required approximately 80 (200) cycles to achieve a converged fission source.

These results clearly showed the effects of accelerating the fission source convergence.We also observed that,with a larger relaxation factor, the FP method was able to achieve faster fission source convergence acceleration,while introducing more oscillation into the process.

The detailed results of calculation times and effective multiplication factor, keff, have been listed in Table 1,where it can be observed that the FP method was able to achieve consistent final keffresults compared with the standard power iteration, proving that the FP method was unbiased.These results also showed that the FP method not only reduced the converged cycles, like the traditional super-history method or Wielandt method,but also reduced the converged calculation times by over 80% and even 90% compared to the standard power iteration, showing that the FP method was able to improve real fission source convergence calculation efficiencies.

3.2 Hoogenboom full-core problem

Fig. 6 Comparison of fission source fractions in region 1 achieved using different methods

The second test problem involved a typical PWR fullcore problem, the Hoogenboom–Martin full-core PWR benchmark problem (commonly referred to as the Hoogenboom problem), which was used to test increasing MC full-core calculation performance [19]. The Hoogenboom problem consists of a typical PWR core layout,with 241 fuel assemblies, each with a 17 × 17 lattice of fuel pins, including 24 control-rod guide tubes and an instrumentation tube. The fuel is composed of 34 different nuclides. Figures 8 and 9 show the layouts of single assembly and full-core, respectively.

To show the effects of the FP method,two initial fission source distributions were tested in the calculations. The first fission source distribution was the point source distribution sampled at the center point of the core, with the other being the uniform source distribution sampled in all assemblies.In the following analysis,a 1,000,000 neutrons per cycle neutron batch size was used with 1000 total cycles (500 inactive cycles) in the RMC calculations.Shannon entropy was used to perform the fission source convergence, and we set the meshes according to the assemblies in radial equally divided into 10 layers in axial to tally the Shannon entropy and for the FP method.Standard power iterations were performed to obtain the Shannon entropy reference results, and the FP method was applied to show acceleration effects. For the FP method,relaxation factor, r, values of 0.6 and 0.8 were set for comparison, and the upper and lower limits and the momentum rate were set using the values suggested in Sect. 2.

Fig. 7 Comparison of Shannon entropies achieved using different methods

Table 1 Calculation results for the modified one-dimensional slab benchmark

Fig. 8 (Color online) Layout of the single assembly used in the Hoogenboom problem

The evolution of the Shannon entropies, with the calculation cycles performed using the standard power iteration and the FP method with different relaxation factors r,with different initial source distributions can be observed in Figs. 10 and 11. From these results, notably, the standard power iteration needed ～300 inactive cycles to achieve fission source convergence, whereas the FP method required ～50 inactive cycles with r set at both 0.6 and 0.8. Moreover, with the larger relaxation factor, the FP method generated more oscillation in the acceleration process.

The resulting calculation times and effective multiplication factors, keff, have been listed in Tables 2 and 3;evidently, the FP method achieved unbiased final keffresults compared to the standard power iteration. These results also showed that the FP method was able to reduce fission source convergence times by approximately 80–85% in these cases. With the Hoogenboom problem,which is more realistic than the slab benchmark, we were also able to prove that the FP method accelerated fission source convergence effectively and efficiently.

Fig.9 (Color online)Layout of the full-core used in the Hoogenboom problem

Fig. 10 Comparison of the Shannon entropies with different methods using point initial source

4 Conclusion

In this study, a new fission source convergence acceleration method, referred to as the Forced Propagation (FP)method, has been proposed and tested. We have been able to show that, compared to existing fission source convergence acceleration methods, such as super-history, Wielandt, and CMFD methods, the proposed FP method was easy to implement based on standard MC transport calculations, could reduce the converged cycles and convergence times to actually increase fission source convergence efficiency and was not confined to a specific problem geometry. To make the FP method more practical, some stabilization techniques were also proposed and tested.The FP method and stabilization techniques was then implemented in the RMC code and tested using both the modified OECD/NEA one-dimensional slab benchmark and the Hoogenboom PWR full-core problem. The results showed that the FP method performed exceptionally well in fission source convergence acceleration.

Fig. 11 Comparison of the Shannon entropies with different methods using uniform initial source

Table 2 Calculation results for the Hoogenboom problem using point initial source

Table 3 Calculation results for the Hoogenboom problem using uniform initial source