2.1 Description of the Forced Propagation method

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

where $\psi$ indicates the neutron flux distribution; k_eff stands for the effective multiplication factor of the system; $\left(r,E,\text{Ω}\right)$ represent the 3-D position, energy, and angle parameters, respectively; ${\text{Σ}}_t\left(r,E\right)$ denotes the macroscopic total cross-section; ${\text{Σ}}_{\text{s}}\left(r,E^{\prime }\to E,\text{Ω"}\to \text{Ω}\right)$ stands for the macroscopic scattering cross-section; $\chi$ indicates the fission spectrum; $\text{ν}$ represents the average secondary neutron number per fission; and ${\text{Σ}}_{\text{f}}\left(r,E^{\prime }\right)$ denotes the macroscopic fission cross-section.

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^{\left(j\right)}$ indicates the initial fission source distribution of the current cycle (j), and $k^{\left(j\right)}$ represents the current estimate of the effective multiplication factor. In conventional MC power iterations, the source neutrons from $f^{\left(j\right)}$ of the current cycle are tracked randomly, according to the transport process, to obtain the flux distribution, ${\psi }^{\left(j\right)}$. Then, the fission source distribution $f^{\left(j+1\right)}$ used for the next cycle, (j+1), can be simulated, based on ${\psi }^{\left(j\right)}$, along with the transport process, using the fission operator, $F$. The process from the initial fission source distribution, $f^{\left(1\right)},$ to the converged fission source distribution, $f^{\left(n\right)}$, as shown in Fig. 1, could be referred to as fission source propagation, and the operator $P=F{\left(T-S\right)}^{-1}$ could be treated as the propagation operator.

Fig. 1.Full-screen View

Schematic of fission source propagation.

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.

Fig. 2.Full-screen View

Standard power iteration and forced propagation processes.

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

where ${\widetilde{fp}}^{\left(j\right)}\left(r,E,\text{Ω}\right)$ represents the proposed FP operator to be determined by the changing trend of the fission source distribution in each iteration cycle; and $f_{\text{FP}}^{\left(j\right)}$ is the FP fission source distribution for cycle (j). As we examine the MC calculation process in more depth, we can determine that the fission source energy ranges are quite narrow, compared with the full energy range of neutrons in the calculation, and that the angle of each fission neutron is sampled randomly. Thus, the FP operator, ${\widetilde{fp}}^{\left(j\right)}\left(r,E,\text{Ω}\right),$ can be simplified as the operator ${\widetilde{fp}}^{\left(j\right)}\left(r\right)$, only related to the position parameter, $r$, so that the FP method description in Eq. (3) can be presented as shown in Eq. (4):

where $f^{\left(j+1\right)}$ indicates the standard fission source distribution calculated using Eq. (1), based on $f_{\text{FP}}^{\left(j\right)}$, from which we could find that the FP fission source distribution, $f_{\text{FP}}^{\left(j\right)}$, can be presented by multiplying the normal fission source distribution, $f^{\left(j\right)}$, and the FP operator, ${\widetilde{fp}}^{\left(j-1\right)}\left(r\right)$.

In the FP process, the operator ${\widetilde{fp}}^{\left(j\right)}\left(r\right)$ 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, ${\widetilde{fp}}^{\left(j\right)}\left(r\right)$ 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 ${\widetilde{fp}}^{\left(j\right)}\left(r\right)$ in the MC method, we need to establish a set of spatial meshes, and then estimate the discrete value of ${\widetilde{fp}}^{\left(j\right)}\left(r\right)$ for each mesh, m, as $f{p}_m^{\left(j\right)}$. 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].

(2) For cycle (j)’s MC transport calculation, obtain the FP fission sources, $f_{\text{FP}}^{\left(j\right)}$, from the previous cycle, (j-1), which includes the position, energy, angle, and weight for each neutron.

(3) Perform the MC transport calculation for each neutron in $f_{\text{FP}}^{\left(j\right)}$, using the standard MC transport method in the RMC.

(4) Store the standard fission source distribution, $f^{\left(j+1\right)}$, and, according to the strategies in the RMC, all the neutrons in $f^{\left(j+1\right)}$ have the same weight, w⁽⁺¹⁾.

(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 $f{p}_m$ using $f{r}_m$. As we obtain the spatial fractions of the standard fission source distribution in mesh m for cycles (j) and (j-1), as $f{r}_m^{\left(j+1\right)}$ and $f{r}_m^{\left(j\right)}$, then $f{p}_m^{\left(j+1\right)}$ can be calculated as shown in (6):

and theoretically, $f{p}_m$ will converge to 1.0 when the fission source distribution has converged, which ensures that the FP method acceleration process is unbiased.

(7) Estimate the FP fission source distribution, $f_{\text{FP}}$, based on $f{p}_m$ and $f$. From Eq. (4), in the FP method, the FP fission sources, $f_{\text{FP}}$, are the product of the FP operator, $\widetilde{fp}$, and the standard fission sources, $f$. We implemented this multiplication process in the RMC (to obtain FP fission sources) by multiplying the weights of standard fission sources in mesh m with $f{p}_m$ and adjusting the total weight, the same as the total weight of $f^{\left(j+1\right)}$. The adjusted fission source distribution with different weights in different meshes and the same total weight as $f^{\left(j + 1\right)}$ constitutes the FP fission source distribution $f_{\text{FP}}^{\left(j+1\right)}$.

(8) Use the FP fission sources, $f_{\text{FP}}^{\left(j+1\right)}$, to perform the MC transport in cycle (j+1), and repeat the processes (2)–(7) to achieve the FP acceleration method.

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.

Fig. 3.Full-screen View

Comparison between the standard iteration and FP iteration processes.

2.2 Techniques to stabilize FP method

The description in section 2.1 shows that the FP method could theoretically accelerate the fission source convergence process. However, when the results of $f{p}_m$ (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 $f{r}_m$ in 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 $f{p}_m$ in 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, $f{p}_m$ may 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 $f{p}_m$, a relaxation factor, r, is implemented after the $f{p}_m$ results are calculated by Eq. (6), as shown in Eq. (7):

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 $f{p}_m$. 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 $f{p}_m$ process could easily derive very large results and the relaxation factor could not perform quite well. Thusly, the upper and lower limits of $f{p}_m$ are 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 $f{p}_m$ was 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.Full-screen View

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), $f{p}_m$ could still converge to 1.0 to ensure that the FP acceleration process is unbiased. Considering the momentum effect properly, the underestimation of $f{p}_m$ caused 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.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 2,000 cycles (1,500 inactive cycles).

Fig. 5.Full-screen View

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

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.

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 ~1,000 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.

Fig. 6.Full-screen View

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

Fig. 7.Full-screen View

Comparison of Shannon entropies achieved using different methods.

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, k_eff, have been listed in Table 1, where it can be observed that the FP method was able to achieve consistent final k_eff results 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

The second test problem involved a typical PWR full-core 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.

Fig. 8Full-screen View

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

Fig. 9Full-screen View

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

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 1,000 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.

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.

Fig. 10.Full-screen View

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

Fig. 11.Full-screen View

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

The resulting calculation times and effective multiplication factors, k_eff, have been listed in Tables 2 and 3; evidently, the FP method achieved unbiased final k_eff results 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.