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

Forced propagation method for Monte Carlo fission source convergence acceleration in the RMC 1 corresp first-author Corresponding author, lizeguang@mail.tsinghua.edu.cn Corresponding author, lizeguang@mail.tsinghua.edu.cn lizeguang@mail.tsinghua.edu.cn https://orcid.org/0000-0002-3457-8471 2 https://orcid.org/0000-0003-0704-8069 2 2 https://orcid.org/0000-0001-7633-9391 Institute of Nuclear and New Energy Technology, Tsinghua University, Beijing 100084, China Institute of Nuclear and New Energy Technology, Tsinghua University, Beijing 100084, China Department of Engineering Physics, Tsinghua University, Beijing 100084, China Department of Engineering Physics, Tsinghua University, Beijing 100084, China 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. Fission source convergence acceleration Monte Carlo method Forced propagation method RMC code 1. Introduction 1 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, k eff , and flux; however, it consumes a ample computing time, especially for loosely coupled or large-scale problems with high dominance ratios [ 1 1 , 2 2 ], such as full-core problems. Therefore, accelerating fission s ource 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 3 , 4 4 ], stratified source sampling method [ 5 5 ], Wielandt method [ 6 6 ], functional MC (FMC) method [ 7 7 ], coarse-mesh finite difference (CMFD) acceleration method [ 8 8 , 9 9 ], asymtotic Wielandt method (AWM), asymptotic super-history method (ASM) [ 10 10 ], and others [ 11 11 - 13 13 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 3 , 4 4 , 6 6 , 14 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 multi-group 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 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 1 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 1 , 2 2 ], such as a loosely coupled case with few neutron interactions between fissile media or large-scale 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 2 The neutron transport equation of criticality is calculated using the equations listed below: <math display="block" id="M1"><mtable columnalign="left"><mtr><mtd><mstyle mathvariant="bold" mathsize="normal"><mi>T</mi></mstyle><mi>ψ</mi><mo>−</mo><mstyle mathvariant="bold" mathsize="normal"><mi>S</mi></mstyle><mi>ψ</mi><mo>=</mo><mfrac><mn>1</mn><mrow><msub><mi>k</mi><mrow><mtext>eff</mtext></mrow></msub></mrow></mfrac><mstyle mathvariant="bold" mathsize="normal"><mi>F</mi></mstyle><mi>ψ</mi><mo>,</mo></mtd></mtr><mtr><mtd><mtext>Transport: </mtext><mstyle mathvariant="bold" mathsize="normal"><mi>T</mi></mstyle><mi>ψ</mi><mo>=</mo><mtext>Ω</mtext><mo>⋅</mo><mo>∇</mo><mi>ψ</mi><mrow><mo>(</mo><mrow><mi>r</mi><mo>,</mo><mi>E</mi><mo>,</mo><mtext>Ω</mtext></mrow><mo>)</mo></mrow><mo>+</mo><msub><mtext>Σ</mtext><mi>t</mi></msub><mrow><mo>(</mo><mrow><mi>r</mi><mo>,</mo><mi>E</mi></mrow><mo>)</mo></mrow><mi>ψ</mi><mrow><mo>(</mo><mrow><mi>r</mi><mo>,</mo><mi>E</mi><mo>,</mo><mtext>Ω</mtext></mrow><mo>)</mo></mrow><mo>,</mo></mtd></mtr><mtr><mtd><mtext>Scattering: </mtext><mstyle mathvariant="bold" mathsize="normal"><mi>S</mi></mstyle><mi>ψ</mi><mo>=</mo><mo>∬</mo><mtext>d</mtext><mi>E</mi><mo>'</mo><msub><mtext>dΩ"Σ</mtext><mi>s</mi></msub><mrow><mo>(</mo><mrow><mi>r</mi><mo>,</mo><msup><mi>E</mi><mo>′</mo></msup><mo>→</mo><mi>E</mi><mo>,</mo><mtext>Ω"</mtext><mo>→</mo><mtext>Ω</mtext></mrow><mo>)</mo></mrow><mi>ψ</mi><mrow><mo>(</mo><mrow><mi>r</mi><mo>,</mo><mi>E</mi><mo>'</mo><mo>,</mo><mtext>Ω"</mtext></mrow><mo>)</mo></mrow><mo>,</mo></mtd></mtr><mtr><mtd><mtext>Fission: </mtext><mstyle mathvariant="bold" mathsize="normal"><mi>F</mi></mstyle><mi>ψ</mi><mo>=</mo><mfrac><mrow><mi>χ</mi><mrow><mo>(</mo><mi>E</mi><mo>)</mo></mrow></mrow><mrow><mn>4</mn><mi>π</mi></mrow></mfrac><mo>∬</mo><mtext>d</mtext><mi>E</mi><mo>'</mo><msub><mtext>dΩ"νΣ</mtext><mi>f</mi></msub><mrow><mo>(</mo><mrow><mi>r</mi><mo>,</mo><msup><mi>E</mi><mo>′</mo></msup></mrow><mo>)</mo></mrow><mi>ψ</mi><mrow><mo>(</mo><mrow><mi>r</mi><mo>,</mo><mi>E</mi><mo>'</mo><mo>,</mo><mtext>Ω"</mtext></mrow><mo>)</mo></mrow><mo>,</mo></mtd></mtr></mtable></math> where <math id="M2"><mi>ψ</mi></math> indicates the neutron flux distribution; k eff stands for the effective multiplication factor of the system; <math id="M3"><mrow><mrow><mo>(</mo><mrow><mi>r</mi><mo>,</mo><mi>E</mi><mo>,</mo><mtext>Ω</mtext></mrow><mo>)</mo></mrow></mrow></math> represent the 3-D position, energy, and angle parameters, respectively; <math id="M4"><mrow><msub><mtext>Σ</mtext><mi>t</mi></msub><mrow><mo>(</mo><mrow><mi>r</mi><mo>,</mo><mi>E</mi></mrow><mo>)</mo></mrow></mrow></math> denotes the macroscopic total cross-section; <math id="M5"><mrow><msub><mtext>Σ</mtext><mtext>s</mtext></msub><mrow><mo>(</mo><mrow><mi>r</mi><mo>,</mo><msup><mi>E</mi><mo>′</mo></msup><mo>→</mo><mi>E</mi><mo>,</mo><mtext>Ω"</mtext><mo>→</mo><mtext>Ω</mtext></mrow><mo>)</mo></mrow></mrow></math> stands for the macroscopic scattering cross-section; <math id="M6"><mi>χ</mi></math> indicates the fission spectrum; <math id="M7"><mtext>ν</mtext></math> represents the average secondary neutron number per fission; and <math id="M8"><mrow><msub><mtext>Σ</mtext><mtext>f</mtext></msub><mrow><mo>(</mo><mrow><mi>r</mi><mo>,</mo><msup><mi>E</mi><mo>′</mo></msup></mrow><mo>)</mo></mrow></mrow></math> 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): <math display="block" id="M9"><mrow><msup><mi>f</mi><mrow><mrow><mo>(</mo><mrow><mi>j</mi><mo>+</mo><mn>1</mn></mrow><mo>)</mo></mrow></mrow></msup><mo>=</mo><mstyle mathvariant="bold" mathsize="normal"><mi>F</mi></mstyle><msup><mi>ψ</mi><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msup><mo>=</mo><mfrac><mn>1</mn><mrow><msup><mi>k</mi><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msup></mrow></mfrac><mstyle mathvariant="bold" mathsize="normal"><mi>F</mi></mstyle><msup><mrow><mrow><mo>(</mo><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>T</mi></mstyle><mo>−</mo><mstyle mathvariant="bold" mathsize="normal"><mi>S</mi></mstyle></mrow><mo>)</mo></mrow></mrow><mrow><mo>−</mo><mn>1</mn></mrow></msup><msup><mi>f</mi><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msup><mo>,</mo></mrow></math> where <math id="M10"><mrow><msup><mi>f</mi><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msup></mrow></math> indicates the initial fission source distribution of the current cycle ( j ), and <math id="M11"><mrow><msup><mi>k</mi><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msup></mrow></math> represents the current estimate of the effective multiplication factor. In conventional MC power iterations, the source neutrons from <math id="M12"><mrow><msup><mi>f</mi><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msup></mrow></math> of the current cycle are tracked randomly, according to the transport process, to obtain the flux distribution, <math id="M13"><mrow><msup><mi>ψ</mi><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msup></mrow></math> . Then, the fission source distribution <math id="M14"><mrow><msup><mi>f</mi><mrow><mrow><mo>(</mo><mrow><mi>j</mi><mo>+</mo><mn>1</mn></mrow><mo>)</mo></mrow></mrow></msup></mrow></math> used for the next cycle, ( j +1), can be simulated, based on <math id="M15"><mrow><msup><mi>ψ</mi><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msup></mrow></math> , along with the transport process, using the fission operator, <math id="M16"><mstyle mathvariant="bold" mathsize="normal"><mi>F</mi></mstyle></math> . The process from the initial fission source distribution, <math id="M17"><mrow><msup><mi>f</mi><mrow><mrow><mo>(</mo><mn>1</mn><mo>)</mo></mrow></mrow></msup><mo>,</mo></mrow></math> to the converged fission source distribution, <math id="M18"><mrow><msup><mi>f</mi><mrow><mrow><mo>(</mo><mi>n</mi><mo>)</mo></mrow></mrow></msup></mrow></math> , as shown in Fig. 1 Fig. 1 , could be referred to as fission source propagation, and the operator <math id="M19"><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>P</mi></mstyle><mo>=</mo><mstyle mathvariant="bold" mathsize="normal"><mi>F</mi></mstyle><msup><mrow><mrow><mo>(</mo><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>T</mi></mstyle><mo>−</mo><mstyle mathvariant="bold" mathsize="normal"><mi>S</mi></mstyle></mrow><mo>)</mo></mrow></mrow><mrow><mo>−</mo><mn>1</mn></mrow></msup></mrow></math> could 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, <math id="M20"><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>P</mi></mstyle><mo>,</mo></mrow></math> as shown in Fig. 2 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: <math display="block" id="M21"><mrow><msubsup><mi>f</mi><mrow><mtext>FP</mtext></mrow><mrow><mrow><mo>(</mo><mrow><mi>j</mi><mo>+</mo><mn>1</mn></mrow><mo>)</mo></mrow></mrow></msubsup><mo>=</mo><mfrac><mn>1</mn><mrow><msup><mi>k</mi><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msup></mrow></mfrac><msup><mrow><mover accent="true"><mrow><mi>f</mi><mi>p</mi></mrow><mo stretchy="true">˜</mo></mover></mrow><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msup><mrow><mo>(</mo><mrow><mi>r</mi><mo>,</mo><mi>E</mi><mo>,</mo><mtext>Ω</mtext></mrow><mo>)</mo></mrow><mo>⋅</mo><mstyle mathvariant="bold" mathsize="normal"><mi>P</mi></mstyle><msubsup><mi>f</mi><mrow><mtext>FP</mtext></mrow><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msubsup><mo>,</mo></mrow></math> where <math id="M22"><mrow><msup><mrow><mover accent="true"><mrow><mi>f</mi><mi>p</mi></mrow><mo stretchy="true">˜</mo></mover></mrow><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msup><mrow><mo>(</mo><mrow><mi>r</mi><mo>,</mo><mi>E</mi><mo>,</mo><mtext>Ω</mtext></mrow><mo>)</mo></mrow></mrow></math> represents the proposed FP operator to be determined by the changing trend of the fission source distribution in each iteration cycle; and <math id="M23"><mrow><msubsup><mi>f</mi><mrow><mtext>FP</mtext></mrow><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msubsup></mrow></math> 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, <math id="M24"><mrow><msup><mrow><mover accent="true"><mrow><mi>f</mi><mi>p</mi></mrow><mo stretchy="true">˜</mo></mover></mrow><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msup><mrow><mo>(</mo><mrow><mi>r</mi><mo>,</mo><mi>E</mi><mo>,</mo><mtext>Ω</mtext></mrow><mo>)</mo></mrow><mo>,</mo></mrow></math> can be simplified as the operator <math id="M25"><mrow><msup><mrow><mover accent="true"><mrow><mi>f</mi><mi>p</mi></mrow><mo stretchy="true">˜</mo></mover></mrow><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msup><mrow><mo>(</mo><mi>r</mi><mo>)</mo></mrow></mrow></math> , only related to the position parameter, <math id="M26"><mi>r</mi></math> , so that the FP method description in Eq. (3) can be presented as shown in Eq. (4): <math display="block" id="M27"><mrow><msubsup><mi>f</mi><mrow><mtext>FP</mtext></mrow><mrow><mrow><mo>(</mo><mrow><mi>j</mi><mo>+</mo><mn>1</mn></mrow><mo>)</mo></mrow></mrow></msubsup><mo>=</mo><msup><mrow><mover accent="true"><mrow><mi>f</mi><mi>p</mi></mrow><mo stretchy="true">˜</mo></mover></mrow><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msup><mrow><mo>(</mo><mi>r</mi><mo>)</mo></mrow><mo>⋅</mo><msup><mi>f</mi><mrow><mrow><mo>(</mo><mrow><mi>j</mi><mo>+</mo><mn>1</mn></mrow><mo>)</mo></mrow></mrow></msup><mo>=</mo><mfrac><mn>1</mn><mrow><msup><mi>k</mi><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msup></mrow></mfrac><msup><mrow><mover accent="true"><mrow><mi>f</mi><mi>p</mi></mrow><mo stretchy="true">˜</mo></mover></mrow><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msup><mrow><mo>(</mo><mi>r</mi><mo>)</mo></mrow><mo>⋅</mo><mstyle mathvariant="bold" mathsize="normal"><mi>P</mi></mstyle><msubsup><mi>f</mi><mrow><mtext>FP</mtext></mrow><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msubsup><mo>,</mo></mrow></math> where <math id="M28"><mrow><msup><mi>f</mi><mrow><mrow><mo>(</mo><mrow><mi>j</mi><mo>+</mo><mn>1</mn></mrow><mo>)</mo></mrow></mrow></msup></mrow></math> indicates the standard fission source distribution calculated using Eq. (1), based on <math id="M29"><mrow><msubsup><mi>f</mi><mrow><mtext>FP</mtext></mrow><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msubsup></mrow></math> , from which we could find that the FP fission source distribution, <math id="M30"><mrow><msubsup><mi>f</mi><mrow><mtext>FP</mtext></mrow><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msubsup></mrow></math> , can be presented by multiplying the normal fission source distribution, <math id="M31"><mrow><msup><mi>f</mi><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msup></mrow></math> , and the FP operator, <math id="M32"><mrow><msup><mrow><mover accent="true"><mrow><mi>f</mi><mi>p</mi></mrow><mo stretchy="true">˜</mo></mover></mrow><mrow><mrow><mo>(</mo><mrow><mi>j</mi><mo>−</mo><mn>1</mn></mrow><mo>)</mo></mrow></mrow></msup><mrow><mo>(</mo><mi>r</mi><mo>)</mo></mrow></mrow></math> . In the FP process, the operator <math id="M33"><mrow><msup><mrow><mover accent="true"><mrow><mi>f</mi><mi>p</mi></mrow><mo stretchy="true">˜</mo></mover></mrow><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msup><mrow><mo>(</mo><mi>r</mi><mo>)</mo></mrow></mrow></math> 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, <math id="M34"><mrow><msup><mrow><mover accent="true"><mrow><mi>f</mi><mi>p</mi></mrow><mo stretchy="true">˜</mo></mover></mrow><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msup><mrow><mo>(</mo><mi>r</mi><mo>)</mo></mrow></mrow></math> 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 <math id="M35"><mrow><msup><mrow><mover accent="true"><mrow><mi>f</mi><mi>p</mi></mrow><mo stretchy="true">˜</mo></mover></mrow><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msup><mrow><mo>(</mo><mi>r</mi><mo>)</mo></mrow></mrow></math> in the MC method, we need to establish a set of spatial meshes, and then estimate the discrete value of <math id="M36"><mrow><msup><mrow><mover accent="true"><mrow><mi>f</mi><mi>p</mi></mrow><mo stretchy="true">˜</mo></mover></mrow><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msup><mrow><mo>(</mo><mi>r</mi><mo>)</mo></mrow></mrow></math> for each mesh, m, as <math id="M37"><mrow><mi>f</mi><msubsup><mi>p</mi><mi>m</mi><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msubsup></mrow></math> . 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 16 ]. (2) For cycle ( j )’s MC transport calculation, obtain the FP fission sources, <math id="M38"><mrow><msubsup><mi>f</mi><mrow><mtext>FP</mtext></mrow><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msubsup></mrow></math> , 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 <math id="M39"><mrow><msubsup><mi>f</mi><mrow><mtext>FP</mtext></mrow><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msubsup></mrow></math> , using the standard MC transport method in the RMC. (4) Store the standard fission source distribution, <math id="M40"><mrow><msup><mi>f</mi><mrow><mrow><mo>(</mo><mrow><mi>j</mi><mo>+</mo><mn>1</mn></mrow><mo>)</mo></mrow></mrow></msup></mrow></math> , and, according to the strategies in the RMC, all the neutrons in <math id="M41"><mrow><msup><mi>f</mi><mrow><mrow><mo>(</mo><mrow><mi>j</mi><mo>+</mo><mn>1</mn></mrow><mo>)</mo></mrow></mrow></msup></mrow></math> 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): <math display="block" id="M42"><mrow><mi>f</mi><msubsup><mi>r</mi><mi>m</mi><mrow><mrow><mo>(</mo><mrow><mi>j</mi><mo>+</mo><mn>1</mn></mrow><mo>)</mo></mrow></mrow></msubsup><mo>=</mo><mfrac><mrow><msub><mstyle mathsize="140%" displaystyle="true"><mo>∑</mo></mstyle><mrow><msup><mi>f</mi><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msup><mo>∈</mo><mi>m</mi></mrow></msub><msup><mi>w</mi><mrow><mrow><mo>(</mo><mrow><mi>j</mi><mo>+</mo><mn>1</mn></mrow><mo>)</mo></mrow></mrow></msup></mrow><mrow><msub><mstyle mathsize="140%" displaystyle="true"><mo>∑</mo></mstyle><mrow><msup><mi>f</mi><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msup><mo>∈</mo><mi>M</mi></mrow></msub><msup><mi>w</mi><mrow><mrow><mo>(</mo><mrow><mi>j</mi><mo>+</mo><mn>1</mn></mrow><mo>)</mo></mrow></mrow></msup></mrow></mfrac><mo>.</mo></mrow></math> (6) Calculate the discrete <math id="M43"><mrow><mi>f</mi><msub><mi>p</mi><mi>m</mi></msub></mrow></math> using <math id="M44"><mrow><mi>f</mi><msub><mi>r</mi><mi>m</mi></msub></mrow></math> . As we obtain the spatial fractions of the standard fission source distribution in mesh m for cycles ( j ) and ( j -1), as <math id="M45"><mrow><mi>f</mi><msubsup><mi>r</mi><mi>m</mi><mrow><mrow><mo>(</mo><mrow><mi>j</mi><mo>+</mo><mn>1</mn></mrow><mo>)</mo></mrow></mrow></msubsup></mrow></math> and <math id="M46"><mrow><mi>f</mi><msubsup><mi>r</mi><mi>m</mi><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msubsup></mrow></math> , then <math id="M47"><mrow><mi>f</mi><msubsup><mi>p</mi><mi>m</mi><mrow><mrow><mo>(</mo><mrow><mi>j</mi><mo>+</mo><mn>1</mn></mrow><mo>)</mo></mrow></mrow></msubsup></mrow></math> can be calculated as shown in (6): <math display="block" id="M48"><mrow><mi>f</mi><msubsup><mi>p</mi><mi>m</mi><mrow><mrow><mo>(</mo><mrow><mi>j</mi><mo>+</mo><mn>1</mn></mrow><mo>)</mo></mrow></mrow></msubsup><mo>=</mo><mfrac><mrow><mi>f</mi><msubsup><mi>r</mi><mi>m</mi><mrow><mrow><mo>(</mo><mrow><mi>j</mi><mo>+</mo><mn>1</mn></mrow><mo>)</mo></mrow></mrow></msubsup></mrow><mrow><mi>f</mi><msubsup><mi>r</mi><mi>m</mi><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msubsup></mrow></mfrac><mo>,</mo></mrow></math> and theoretically, <math id="M49"><mrow><mi>f</mi><msub><mi>p</mi><mi>m</mi></msub></mrow></math> 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, <math id="M50"><mrow><msub><mi>f</mi><mrow><mtext>FP</mtext></mrow></msub></mrow></math> , based on <math id="M51"><mrow><mi>f</mi><msub><mi>p</mi><mi>m</mi></msub></mrow></math> and <math id="M52"><mi>f</mi></math> . From Eq. (4), in the FP method, the FP fission sources, <math id="M53"><mrow><msub><mi>f</mi><mrow><mtext>FP</mtext></mrow></msub></mrow></math> , are the product of the FP operator, <math id="M54"><mrow><mover accent="true"><mrow><mi>f</mi><mi>p</mi></mrow><mo stretchy="true">˜</mo></mover></mrow></math> , and the standard fission sources, <math id="M55"><mi>f</mi></math> . 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 <math id="M56"><mrow><mi>f</mi><msub><mi>p</mi><mi>m</mi></msub></mrow></math> and adjusting the total weight, the same as the total weight of <math id="M57"><mrow><msup><mi>f</mi><mrow><mrow><mo>(</mo><mrow><mi>j</mi><mo>+</mo><mn>1</mn></mrow><mo>)</mo></mrow></mrow></msup></mrow></math> . The adjusted fission source distribution with different weights in different meshes and the same total weight as <math id="M58"><mrow><msup><mi>f</mi><mrow><mrow><mo>(</mo><mrow><mi>j</mi><mo> </mo><mo>+</mo><mo> </mo><mn>1</mn></mrow><mo>)</mo></mrow></mrow></msup></mrow></math> constitutes the FP fission source distribution <math id="M59"><mrow><msubsup><mi>f</mi><mrow><mtext>FP</mtext></mrow><mrow><mrow><mo>(</mo><mrow><mi>j</mi><mo>+</mo><mn>1</mn></mrow><mo>)</mo></mrow></mrow></msubsup></mrow></math> . (8) Use the FP fission sources, <math id="M60"><mrow><msubsup><mi>f</mi><mrow><mtext>FP</mtext></mrow><mrow><mrow><mo>(</mo><mrow><mi>j</mi><mo>+</mo><mn>1</mn></mrow><mo>)</mo></mrow></mrow></msubsup></mrow></math> , 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 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 2 The description in section 2.1 shows that the FP method could theoretically accelerate the fission source convergence process. However, when the results of <math id="M61"><mrow><mi>f</mi><msub><mi>p</mi><mi>m</mi></msub></mrow></math> (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 <math id="M62"><mrow><mi>f</mi><msub><mi>r</mi><mi>m</mi></msub></mrow></math> 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 <math id="M63"><mrow><mi>f</mi><msub><mi>p</mi><mi>m</mi></msub></mrow></math> 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, <math id="M64"><mrow><mi>f</mi><msub><mi>p</mi><mi>m</mi></msub></mrow></math> 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 <math id="M65"><mrow><mi>f</mi><msub><mi>p</mi><mi>m</mi></msub></mrow></math> , a relaxation factor, r , is implemented after the <math id="M66"><mrow><mi>f</mi><msub><mi>p</mi><mi>m</mi></msub></mrow></math> results are calculated by Eq. (6), as shown in Eq. (7): <math display="block" id="M67"><mrow><mi>f</mi><msubsup><mi>p</mi><mi>m</mi><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msubsup><mo>=</mo><mi>r</mi><mo>⋅</mo><mi>f</mi><msubsup><mi>p</mi><mi>m</mi><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msubsup><mo>+</mo><mrow><mo>(</mo><mrow><mn>1</mn><mo>−</mo><mi>r</mi></mrow><mo>)</mo></mrow><mo>⋅</mo><mn>1.0</mn><mo>,</mo></mrow></math> 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 <math id="M68"><mrow><mi>f</mi><msub><mi>p</mi><mi>m</mi></msub></mrow></math> . 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 <math id="M69"><mrow><mi>f</mi><msub><mi>p</mi><mi>m</mi></msub></mrow></math> process could easily derive very large results and the relaxation factor could not perform quite well. Thusly, the upper and lower limits of <math id="M70"><mrow><mi>f</mi><msub><mi>p</mi><mi>m</mi></msub></mrow></math> 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 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 Fig. 4 . In the RMC, we also used momentum consideration in the FP method, whereby <math id="M71"><mrow><mi>f</mi><msub><mi>p</mi><mi>m</mi></msub></mrow></math> 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): <math display="block" id="M72"><mtable columnalign="left"><mtr><mtd><mi>m</mi><msubsup><mi>t</mi><mi>m</mi><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msubsup><mo>=</mo><mrow><mo>(</mo><mrow><mn>1</mn><mo>−</mo><mi>t</mi></mrow><mo>)</mo></mrow><mo>⋅</mo><mi>f</mi><msubsup><mi>p</mi><mi>m</mi><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msubsup><mo>+</mo><mi>t</mi><mo>⋅</mo><mi>m</mi><msubsup><mi>t</mi><mi>m</mi><mrow><mrow><mo>(</mo><mrow><mi>j</mi><mo>−</mo><mn>1</mn></mrow><mo>)</mo></mrow></mrow></msubsup><mo>,</mo></mtd></mtr><mtr><mtd><mi>f</mi><msubsup><mi>p</mi><mi>m</mi><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msubsup><mo>=</mo><mi>m</mi><msubsup><mi>t</mi><mi>m</mi><mrow><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow></mrow></msubsup><mo>,</mo></mtd></mtr></mtable></math> where t indicates the momentum rate set between 0.0 and 1.0, and, according to Eq. (8), <math id="M73"><mrow><mi>f</mi><msub><mi>p</mi><mi>m</mi></msub></mrow></math> could still converge to 1.0 to ensure that the FP acceleration process is unbiased. Considering the momentum effect properly, the underestimation of <math id="M74"><mrow><mi>f</mi><msub><mi>p</mi><mi>m</mi></msub></mrow></math> 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 2 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 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 1 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 2 The first test problem was a modified OECD / NEA source convergence benchmark [ 18 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 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). 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 Figs. 6 and 7 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 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, k eff , have been listed in Table 1 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. 1 1 Method 1 1 Final k eff 1 1 Converged cycles (Estimation) 1 1 Convergence time (min) 1 1 Power iteration 1 1 0.921268 ± 0.000040 1 1 1,000 1 1 21.49 1 1 FP method, r = 0.8 1 1 0.921296 ± 0.000040 1 1 80 1 1 1.94 1 1 FP method, r = 0.6 1 1 0.921255 ± 0.000039 1 1 200 1 1 4.12 3.2 Hoogenboom full-core problem 2 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 19 ]. The Hoogenboom problem con sists 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 Figures 8 and 9 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 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 Figs. 10 and 11 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, k eff , have been listed in Tables 2 Tables 2 and 3 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. 1 1 Method 1 1 Final k eff 1 1 Converged cycles (estimate) 1 1 Convergence time (min) 1 1 Power iteration 1 1 1.000508 ± 0.000029 1 1 300 1 1 35.01 1 1 FP method, r = 0.8 1 1 1.000461 ± 0.000028 1 1 50 1 1 5.86 1 1 FP method, r = 0.6 1 1 1.000487 ± 0.000028 1 1 50 1 1 5.83 1 1 Method 1 1 Final k eff 1 1 Converged cycles (estimate) 1 1 Convergence time (min) 1 1 Power iteration 1 1 1.000535 ± 0.000028 1 1 300 1 1 34.83 1 1 FP method, r = 0.8 1 1 1.000496 ± 0.000029 1 1 50 1 1 5.79 1 1 FP method, r = 0.6 1 1 1.000507 ± 0.000029 1 1 50 1 1 5.91 4. Conclusion 1 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. References S. Gorodkov , Using core symmetry in Monte Carlo dominance ratio calculations , Nucl. Sci. Eng . 168 ( 3 ), 242 - 247 ( 2011 ). doi: 10.13182/NSE10-37 http://doi.org/10.13182/NSE10-37 Using core symmetry in Monte Carlo dominance ratio calculations T. Ueki , F.B. Brown , D.K. Parsons et al ., Autocorrelation and Dominance Ratio in Monte Carlo Criticality Calculations , Nucl. Sci. Eng . 145 ( 3 ), 279 - 290 ( 2003 ). doi: 10.13182/NSE03-04 http://doi.org/10.13182/NSE03-04 Autocorrelation and Dominance Ratio in Monte Carlo Criticality Calculations R.J. Brissenden , A.R. Garlick , Biases in the estimation of keff and its error by Monte Carlo methods , Ann. Nucl. Energy . 13 ( 2 ), 63 - 83 ( 1986 ). doi: 10.1016/0306-4549(86)90095-2 http://doi.org/10.1016/0306-4549(86)90095-2 Biases in the estimation of keff and its error by Monte Carlo methods R.N. Blomquist , E.M. Gelbard , Alternative implementations of the Monte Carlo power method , Nucl. Sci. Eng . 141 , 85 - 100 ( 2002 ). doi: 10.13182/NSE01-30 http://doi.org/10.13182/NSE01-30 Alternative implementations of the Monte Carlo power method A. Mohamed , E.M. Gelbard , Stratified Source-Sampling Techniques for Monte Carlo Eigenvalue Analysis , Proc. Int. Conf. Physics of Nuclear Science and Technology , Hauppage, New York , October 5–8, 1998 . Stratified Source-Sampling Techniques for Monte Carlo Eigenvalue Analysis T. Yamamoto , Y. Miyoshi , Reliable method for fission source convergence of Monte Carlo criticality calculation with Wielandt’s method , J. Nucl. Sci. Technol . 41 , 99 - 107 ( 2004 ). DOI: 10.1080/18811248.2004.9715465 http://doi.org/10.1080/18811248.2004.9715465 Reliable method for fission source convergence of Monte Carlo criticality calculation with Wielandt’s method E.W. Larsen , J. Yang , A functional Monte Carlo method for k-eigenvalue problems , Nucl. Sci. Eng . 159 ( 2 ), 107 - 126 ( 2008 ). doi: 10.13182/NSE07-92 http://doi.org/10.13182/NSE07-92 A functional Monte Carlo method for k-eigenvalue problems M.J. Lee , H. Joo , D. Lee et al ., Investigation of CMFD ACCELERATED Monte Carlo Eigenvalue Calculation with Simplified Low Dimensional Multigroup Formulation , International Conference on the Physics of Reactors 2010, PHYSOR 2010 , 2010 . Investigation of CMFD ACCELERATED Monte Carlo Eigenvalue Calculation with Simplified Low Dimensional Multigroup Formulation Y.G. Jo , N.Z. Cho , Acceleration and real variance reduction in Continuous-Energy Monte Carlo whole-core calculation via p-CMFD feedback , Nucl. Sci. Eng . 189 , 26 - 40 ( 2018 ). doi: 10.1080/00295639.2017.1373517 http://doi.org/10.1080/00295639.2017.1373517 Acceleration and real variance reduction in Continuous-Energy Monte Carlo whole-core calculation via p-CMFD feedback D. She , K. Wang , G.L. Yu , Asymptotic Wielandt method and superhistory method for source convergence in Monte Carlo criticality calculation , Nucl. Sci. Eng . 172 ( 2 ), 127 - 137 ( 2012 ). doi: 10.13182/NSE11-44 http://doi.org/10.13182/NSE11-44 Asymptotic Wielandt method and superhistory method for source convergence in Monte Carlo criticality calculation J. Dufek , W. Gudowski , Stability and convergence problems of the Monte Carlo fission matrix acceleration methods , Ann. Nucl. Energy . 36 , 1648 - 1651 ( 2009 ). doi: 10.1016/j.anucene.2009.07.020 http://doi.org/10.1016/j.anucene.2009.07.020 Stability and convergence problems of the Monte Carlo fission matrix acceleration methods M.R. Omar , J.A. Karim , T.L. Yoon , An efficient Monte Carlo fission source convergence acceleration strategy adapted from the survival biasing technique , Ann. Nucl. Energy . 138 , 107 - 164 ( 2020 ). doi: 10.1016/j.anucene.2019.107164 http://doi.org/10.1016/j.anucene.2019.107164 An efficient Monte Carlo fission source convergence acceleration strategy adapted from the survival biasing technique B. Toth , D. Griesheimer , A novel source convergence acceleration scheme for Monte Carlo criticality calculations, part II: implementation and results , Proceedings of Joint International Topical Meeting on Mathematics & Computation and Supercomputing in Nuclear Applications (M&C + SNA 2007) , Monterey, California , Apr 15-19, 2007 . A novel source convergence acceleration scheme for Monte Carlo criticality calculations, part II: implementation and results J. Dufek , W. Gudoski , Analysis and Improvement of the Fission Matrix Method for Source Convergence Acceleration in Monte Carlo Eigenvalue Calculations, Technical Report, NURESIM Work Package No. T1.1.2., Milestone No. 4c , 2006 . Analysis and Improvement of the Fission Matrix Method for Source Convergence Acceleration in Monte Carlo Eigenvalue Calculations, Technical Report, NURESIM Work Package No. T1.1.2., Milestone No. 4c K. Wang , Z. Li , D. She et al ., RMC - A Monte Carlo code for reactor core analysis , Ann. Nucl. Energy . 82 , 121 - 129 ( 2015 ). doi: 10.1016/j.anucene.2014.08.048 http://doi.org/10.1016/j.anucene.2014.08.048 RMC - A Monte Carlo code for reactor core analysis T. Ueki , F.B. Brown , Stationarity Diagnostics Using Shannon Entropy in Monte Carlo Criticality Calculation I: F Test , Trans. Am. Nucl. Soc . 87 , 156 - 158 ( 2002 ). Stationarity Diagnostics Using Shannon Entropy in Monte Carlo Criticality Calculation I: F Test N. Qian . On the momentum term in gradient descent learning algorithms . Neural networks . 12 ( 1 ), 145 - 151 ( 1999 ). doi: 10.1016/S0893-6080(98)00116-6 http://doi.org/10.1016/S0893-6080(98)00116-6 On the momentum term in gradient descent learning algorithms R. Blomquist , A. Nouri , M. Armishaw et al ., OECD/NEA Source Convergence Benchmark Program: Overview and Summary of Results , Proc. Int. Conf. Nuclear Criticality Safety (ICNC2003) , Tokai-mura, Ibaraki, Japan , October 20–24, 2003 . OECD/NEA Source Convergence Benchmark Program: Overview and Summary of Results J.E. Hoogenboom , W.R. Martin , B. Petrovic , The Monte Carlo performance benchmark test - aims, specifications and first results, Int . Conf. Mathematics and Computational Methods Applied to Nuclear Science and Engineering , Rio de Janeiro, Brazil , 2011 . The Monte Carlo performance benchmark test - aims, specifications and first results, Int All authors contributed to the study conception and design. Material preparation, data collection and analysis were performed by Ze-Guang Li, Yu-Chuan Guo and Xiao-Yu Guo. The first draft of the manuscript was written by Ze-Guang Li and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript. 10.1007/s41365-021-00868-0 This research was supported by the National Natural Science Foundation of China (Nos. 11775126, 11545013, 11605101), the Young Elite Scientists Sponsorship Program by CAST (No. 2016QNRC001), Science Challenge Project by MIIT of China (No. TZ2018001), and Tsinghua University, Initiative Scientific Research Program. 2020-10-07 2020-12-24 2021-01-16 2021-03-16 2021-03-01 2021-03-31 © China Science Publishing & Media Ltd. (Science Press), Shanghai Institute of Applied Physics, the Chinese Academy of Sciences, ChineseNuclear Society and Springer Nature Singapore Pte Ltd. 2021 This is a post-peer-review, pre-copyedit version of an article published in Nuclear Science and Techniques. The final authenticated version is available online at: http://dx.doi.org/10.1007/s41365-021-00868-0 http://dx.doi.org/10.1007/s41365-021-00868-0 . 2021 Nuclear Science and Techniques 03 32