Sample size adaptive strategy for time-dependent Monte Carlo particle transport simulation

Sample size adaptive strategy for time-dependent Monte Carlo particle transport simulation 1 first-author https://orcid.org/0000-0002-5968-3665 1 1 1 1 1 corresp † sgji0906@qq.com † sgji0906@qq.com sgji0906@qq.com Institute of Applied Physics and Computational Mathematics, Beijing 100094, China Institute of Applied Physics and Computational Mathematics, Beijing 100094, China When multiphysics coupling calculations contain time-dependent Monte Carlo particle transport simulations, these simulations often account for the largest part of the calculation time, which is insufferable in certain important cases. This study proposes an adaptive strategy for automatically adjusting the sample size to fulfil more reasonable simulations. This is realized based on an extension of the Shannon entropy concept and is essentially different from the popular methods in time-independent Monte Carlo particle transport simulations, such as controlling the sample size according to the relative error of a target tally or by experience. The results of the two models show that this strategy can yield almost similar results while significantly reducing the calculation time. Considering the efficiency, the sample size should not be increased blindly if the efficiency cannot be enhanced further. The strategy proposed herein satisfies this requirement. Time-dependent Monte Carlo particle transport simulation Shannon entropy Adaptive strategy Introduction 1 Multiphysics processes have been extensively studied in several important fields. Owing to its complexity, an analytical solution is hardly possible, and numerical calculations are treated as efficient tools instead of expensive experiments. When time-dependent processes governed by different physical rules are coupled, one apparent solution is to divide the entire calculation into many steps, in which different control equations should be solved individually. Some data must be transmitted from one process to another to mimic the coupling effect. In certain important cases, the time-dependent particle transport is an indispensable part of the multiphysics process. The coupling between fluid mechanics and neutron transport was used as an example, which was used to demonstrate the sample-size adaptive strategy. In many cases, other physical processes must also be considered. However, considering only these two processes is not trivial and is meaningful in certain important situations. The governing equations for fluid mechanics are as follows: <math display="block" id="M1"><mrow><mfrac><mrow><mo>∂</mo><mi>ρ</mi></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mo>∇</mo><mo>⋅</mo><mi>ρ</mi><mover accent="true"><mi>u</mi><mo>→</mo></mover><mo>=</mo><mn>0</mn><mo>,</mo></mrow></math> <math display="block" id="M2"><mrow><mfrac><mrow><mo>∂</mo><mi>ρ</mi><mover accent="true"><mi>u</mi><mo>→</mo></mover></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mo>∇</mo><mo>⋅</mo><mi>ρ</mi><mover accent="true"><mi>u</mi><mo>→</mo></mover><mover accent="true"><mi>u</mi><mo>→</mo></mover><mo>=</mo><mo>−</mo><mo>∇</mo><mi>p</mi><mo>,</mo></mrow></math> <math display="block" id="M3"><mrow><mfrac><mrow><mo>∂</mo><mi>ρ</mi><mi>E</mi></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mo>∇</mo><mo>⋅</mo><mi>ρ</mi><mi>E</mi><mover accent="true"><mi>u</mi><mo>→</mo></mover><mo>=</mo><mo>−</mo><mo>∇</mo><mo>⋅</mo><mi>p</mi><mover accent="true"><mi>u</mi><mo>→</mo></mover><mo>+</mo><mi>Q</mi><mn>.</mn></mrow></math> In the above equations, ρ is the density, p is the pressure, <math id="M4"><mover accent="true"><mi>u</mi><mo>→</mo></mover></math> is the fluid velocity, <math id="M5"><mrow><mi>E</mi><mo>=</mo><mfrac><mn>1</mn><mn>2</mn></mfrac><mover accent="true"><mi>u</mi><mo>→</mo></mover><mover accent="true"><mi>u</mi><mo>→</mo></mover><mo>+</mo><mi>I</mi></mrow></math> is the energy, I is the material-specific internal energy, and Q is the feedback quantity from neutron transport. The fluid pressure p , density ρ , and internal energy I are related by the equation of state, whose general form is <math display="block" id="M6"><mrow><mi>p</mi><mo>=</mo><mi>f</mi><mo stretchy="false">(</mo><mi>ρ</mi><mn>,</mn><mi>I</mi><mo stretchy="false">)</mo><mn>.</mn></mrow></math> When neutrons fly in a medium, nuclear energy is released by reactions between the neutrons and the nucleus. These reactions change the temperature, density, and velocity of the medium. However, because these reactions rely on their relative velocities, the thermal and hydrodynamic motions of the nucleus influence the neutron transport process. Therefore, the correct transport equation should include the distribution function of the nucleus, except for the neutron flux, which should be determined by the corresponding transport equation. This created a complex problem involving a group of simultaneous transport equations. Fortunately, because the local thermal equilibrium hypothesis is suitable, the use of Maxwell distribution to describe the thermal motion of a nucleus is sufficiently precise. After averaging, the coefficients of the neutron transport equation rely only on the relative velocity of the neutron and the hydrodynamic motion of the medium. Therefore, the neutron transport equation that considers the effect of fluid motion is given by [ 1 1 ] <math display="block" id="M7"><mtable columnalign="left"><mtr><mtd><mfrac><mrow><mo>∂</mo><mi>N</mi></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mover accent="true"><mi>v</mi><mo>→</mo></mover><mo>⋅</mo><mo>∇</mo><mi>N</mi><mo>+</mo><mfrac><mrow><mtext>d</mtext><mover accent="true"><mi>N</mi><mo>→</mo></mover></mrow><mrow><mtext>d</mtext><mi>t</mi></mrow></mfrac><mo>⋅</mo><msub><mo>∇</mo><mi>V</mi></msub><mi>N</mi><mo>+</mo><mi>V</mi><mstyle displaystyle="true"><mo>∑</mo></mstyle><mrow><mo>(</mo><mrow><mover accent="true"><mi>r</mi><mo>→</mo></mover><mn>,</mn><mi>V</mi></mrow><mo>)</mo></mrow><mi>N</mi><mrow><mo>(</mo><mrow><mover accent="true"><mi>r</mi><mo>→</mo></mover><mn>,</mn><mover accent="true"><mi>V</mi><mo>→</mo></mover><mn>,</mn><mi>t</mi></mrow><mo>)</mo></mrow></mtd></mtr><mtr><mtd><mtext>   </mtext><mo>=</mo><mstyle displaystyle="true"><mrow><mo>∫</mo></mrow></mstyle><msup><mstyle displaystyle="true"><mo>∑</mo></mstyle><mo>′</mo></msup><msup><mi>f</mi><mo>′</mo></msup><msup><mi>V</mi><mo>′</mo></msup><msup><mi>N</mi><mo>′</mo></msup><mi>d</mi><msup><mover accent="true"><mi>V</mi><mo>→</mo></mover><mo>′</mo></msup><mo>+</mo><mi>q</mi><mrow><mo>(</mo><mrow><mover accent="true"><mi>r</mi><mo>→</mo></mover><mn>,</mn><mover accent="true"><mi>V</mi><mo>→</mo></mover><mn>,</mn><mi>t</mi></mrow><mo>)</mo></mrow><mn>,</mn></mtd></mtr></mtable></math> where <math id="M8"><mover accent="true"><mi>v</mi><mo>→</mo></mover></math> is the neutron velocity, <math id="M9"><mrow><mover accent="true"><mi>V</mi><mo>→</mo></mover><mo>=</mo><mover accent="true"><mi>v</mi><mo>→</mo></mover><mo>−</mo><mover accent="true"><mi>u</mi><mo>→</mo></mover></mrow></math> is the relative velocity, <math id="M10"><mrow><mi>N</mi><mo stretchy="false">(</mo><mover accent="true"><mi>r</mi><mo>→</mo></mover><mn>,</mn><mover accent="true"><mi>V</mi><mo>→</mo></mover><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo><mtext>d</mtext><mover accent="true"><mi>r</mi><mo>→</mo></mover><mtext>d</mtext><mover accent="true"><mi>V</mi><mo>→</mo></mover></mrow></math> is the mean neutron number in the volume element <math id="M11"><mrow><mtext>d</mtext><mover accent="true"><mi>r</mi><mo>→</mo></mover><mtext>d</mtext><mover accent="true"><mi>V</mi><mo>→</mo></mover></mrow></math> around <math id="M12"><mrow><mo stretchy="false">(</mo><mover accent="true"><mi>r</mi><mo>→</mo></mover><mn>,</mn><mover accent="true"><mi>V</mi><mo>→</mo></mover><mo stretchy="false">)</mo></mrow></math> at time t , and <math id="M13"><mrow><msup><mrow><mstyle displaystyle="true"><mo>∑</mo></mstyle></mrow><mo>′</mo></msup><msup><mi>f</mi><mo>′</mo></msup><msup><mi>V</mi><mo>′</mo></msup><msup><mi>N</mi><mo>′</mo></msup><mo>=</mo><mstyle displaystyle="true"><mo>∑</mo></mstyle><mrow><mo>(</mo><mrow><mover accent="true"><mi>r</mi><mo>→</mo></mover><mn>,</mn><msup><mi>V</mi><mo>′</mo></msup></mrow><mo>)</mo></mrow><mi>f</mi><mrow><mo>(</mo><mrow><mover accent="true"><mi>r</mi><mo>→</mo></mover><mn>,</mn><msup><mover accent="true"><mi>V</mi><mo>→</mo></mover><mo>′</mo></msup><mo>→</mo><mover accent="true"><mi>V</mi><mo>→</mo></mover></mrow><mo>)</mo></mrow><msup><mi>V</mi><mo>′</mo></msup><mi>N</mi><mrow><mo>(</mo><mrow><mover accent="true"><mi>r</mi><mo>→</mo></mover><mn>,</mn><msup><mover accent="true"><mi>V</mi><mo>→</mo></mover><mo>′</mo></msup><mn>,</mn><mi>t</mi></mrow><mo>)</mo></mrow></mrow></math> , <math id="M14"><mrow><mstyle displaystyle="true"><mo>∑</mo></mstyle><mo stretchy="false">(</mo><mover accent="true"><mi>r</mi><mo>→</mo></mover><mn>,</mn><mi>V</mi><mo stretchy="false">)</mo></mrow></math> is the total macro cross section. In the total macro cross section, <math display="block" id="M15"><mrow><mstyle displaystyle="true"><mo>∑</mo></mstyle><mrow><mo>(</mo><mrow><mover accent="true"><mi>r</mi><mo>→</mo></mover><mn>,</mn><msup><mi>V</mi><mo>′</mo></msup></mrow><mo>)</mo></mrow><mi>f</mi><mrow><mo>(</mo><mrow><mover accent="true"><mi>r</mi><mo>→</mo></mover><mn>,</mn><msup><mover accent="true"><mi>V</mi><mo>→</mo></mover><mo>′</mo></msup><mo>→</mo><mover accent="true"><mi>V</mi><mo>→</mo></mover></mrow><mo>)</mo></mrow><mo>=</mo><mstyle displaystyle="true"><munder><mo>∑</mo><mi>i</mi></munder><mrow><mstyle displaystyle="true"><munder><mo>∑</mo><mi>x</mi></munder><mrow><mstyle displaystyle="true"><msubsup><mo>∑</mo><mi>x</mi><mi>i</mi></msubsup><mrow><mrow><mo>(</mo><mrow><mover accent="true"><mi>r</mi><mo>→</mo></mover><mn>,</mn><msup><mi>V</mi><mo>′</mo></msup></mrow><mo>)</mo></mrow></mrow></mstyle></mrow></mstyle></mrow></mstyle><msubsup><mi>f</mi><mi>x</mi><mi>i</mi></msubsup><mrow><mo>(</mo><mrow><mover accent="true"><mi>r</mi><mo>→</mo></mover><mn>,</mn><msup><mover accent="true"><mi>V</mi><mo>→</mo></mover><mo>′</mo></msup><mo>→</mo><mover accent="true"><mi>V</mi><mo>→</mo></mover></mrow><mo>)</mo></mrow><mn>,</mn></mrow></math> where <math id="M16"><mrow><mstyle displaystyle="true"><msubsup><mo>∑</mo><mi>x</mi><mi>i</mi></msubsup><mrow><mrow><mo>(</mo><mrow><mover accent="true"><mi>r</mi><mo>→</mo></mover><mn>,</mn><msup><mi>V</mi><mo>′</mo></msup></mrow><mo>)</mo></mrow></mrow></mstyle></mrow></math> is the probability of x reaction taking place with i nucleus when a neutron travels a unit distance, <math id="M17"><mrow><msubsup><mi>f</mi><mi>x</mi><mi>i</mi></msubsup><mrow><mo>(</mo><mrow><mover accent="true"><mi>r</mi><mo>→</mo></mover><mn>,</mn><msup><mover accent="true"><mi>V</mi><mo>→</mo></mover><mo>′</mo></msup><mo>→</mo><mover accent="true"><mi>V</mi><mo>→</mo></mover></mrow><mo>)</mo></mrow></mrow></math> is the probability of generating neutrons with velocity <math id="M18"><mover accent="true"><mi>V</mi><mo>→</mo></mover></math> when a neutron takes place x reaction with i nucleus at point <math id="M19"><mover accent="true"><mi>r</mi><mo>→</mo></mover></math> , and <math id="M20"><mrow><mi>q</mi><mrow><mo>(</mo><mrow><mover accent="true"><mi>r</mi><mo>→</mo></mover><mn>,</mn><mover accent="true"><mi>V</mi><mo>→</mo></mover><mn>,</mn><mi>t</mi></mrow><mo>)</mo></mrow></mrow></math> is an external source. At the beginning of each step, the fluid calculation provides a new transport environment for neutron transport, and the time-dependent particle transport problem is calculated to obtain feedback quantities for the fluid calculation. This cycle is repeated several times to obtain the required quantities. Based on the above control equations, the Jointed numerical simulation software for mUlti-PhysIcs and multi-maTerials problems under ExtReme conditions (JUPITER) platform was developed by a team at the Institute of Applied Physics and Computational Mathematics. It utilizes mature framework technology to support software development and has been confirmed and verified by many models. Owing to the accuracy of the Monte Carlo method for geometry [ 2 2 ], physics modeling, and other merits, this method is often used to simulate particle transport [ 3 3 ]; however, it is extremely time-consuming. When there is a problem with multi-materials with large deformations and the scale of the model's grid number is large, several steps are required for a reliable calculation, and each step contains a Monte Carlo simulation of a large model. Therefore, in this case, the calculation time is insufferable even with large-scale parallelization using a powerful supercomputer. Similar to the traditional Monte Carlo simulation of the time-independent particle transport problem, the core objective is to achieve the highest possible efficiency. This means finding a method with high calculation accuracy but not time-consuming. In other words, when efficiency flattens with an increase in sample size, there is no need to perform more calculations if efficiency, not calculation accuracy, is the goal. For the Monte Carlo method, it is well known that when the sample size tends to infinity, the calculation accuracy can reach any level because of the law of large numbers, which is impossible. Apparently, it is difficult to set different reasonable sample sizes for each step based on experience because each step is a different particle transport problem. If a fixed, sufficiently large sample size is set for all steps, which is the original sample size setting method of the JUPITER platform, many calculations may be wasted because the efficiency factor will be flat as the sample size is increased to a certain level. Therefore, experience-based methods are not considered to be efficient. Moreover, because the Monte Carlo simulation provides several feedback quantities with very different fluctuation levels, the sample size controlled by the relative error of one tally may not be sufficient or may be excessively large for another tally. Therefore, error-based methods [ 4 4 ] are unsuitable. Because the calculation of time-dependent particle transport is expensive, Monte Carlo simulations still need to develop an adaptive strategy to scientifically set the sample size, even with the rapid progress of supercomputers. The remainder of this paper is organized as follows: The second section describes all aspects of a sample size adaptive strategy based on an extension of the Shannon entropy concept. The third section contains some active numerical results from the calculations of the two models. The final section provides the concluding remarks. Sample size adaptive strategy 1 Research basics 2 Entropy is a notion existing in both physics and mathematics. Clausius is considered the father of the concept of entropy. The important contributions to the development of this concept were made by Boltzmann, Gibbs, and Planck. Moreover, mathematicians, including Shannon, Kolmogorov, Sinai, Rényi, Holevo, and Tsallis, were also preoccupied with the concept. On the history of different paths to entropy concept, see, e.g., [ 5 5 - 7 7 5-7 ]. Therefore, there are various definitions of entropy. In this study, the entropy concept used is the one proposed by Shannon in [ 8 8 ], which has some novel applications in various fields [ 9 9 , 10 10 ]. For example, Ref. [ 9 9 ] showed that Shannon entropy is an efficient dynamic indicator that provides a direct measure of the diffusion rate and, thus, a timescale for the instabilities arising when dealing with chaos. It is worth pointing out that the Shannon entropy was first introduced into Monte Carlo particle transport simulation as a posterior indicator for convergence of fission source distribution in eigenvalue calculation of nuclear systems [ 11 11 - 14 14 11-14 ]. Suppose there are N cells in total for some nuclear systems and M fission neutrons at some moment. Let m i be the number of fission neutrons in cell i . Hence, ∑ i m i = M . Let <math id="M21"><mrow><msub><mi>p</mi><mi>i</mi></msub><mo>=</mo><mfrac><mrow><msub><mi>m</mi><mi>i</mi></msub></mrow><mi>M</mi></mfrac></mrow></math> and <math display="block" id="M22"><mrow><msub><mi>H</mi><mrow><mtext>ori</mtext></mrow></msub><mo>=</mo><mo>−</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>M</mi></munderover><mrow><msub><mi>p</mi><mi>i</mi></msub></mrow></mstyle><mi>ln</mi><mo stretchy="false">(</mo><msub><mi>p</mi><mi>i</mi></msub><mo stretchy="false">)</mo><mn>,</mn></mrow></math> we obtain the Shannon entropy H ori . For the power iteration method of eigenvalue calculation, if there are K iteration steps and <math id="M23"><mrow><msubsup><mi>H</mi><mrow><mtext>ori</mtext></mrow><mi>k</mi></msubsup><mn>,</mn><mi>k</mi><mo>=</mo><mn>1,</mn><mo>⋯</mo><mn>,</mn><mi>K</mi></mrow></math> are entropies of each step, then the convergence of this entropy sequence is a necessary condition for the convergence of the fission source distribution. It is primarily used as follows: When the entire calculation has been completed and the Shannon entropy of each step has been obtained, a step number n is selected by judging the entropy sequence after this step is perceived to have converged. If the number of the first active steps is greater than n , this calculation is considered reasonable. Otherwise, it is judged that the inactive iteration steps are insufficient, and the results are systematically biased. In recent years, there have been some attempts to explore on-the-fly judgment methods [ 15 15 - 18 18 15-18 ] for the convergence of the fission source distribution, which belongs to the time-independent particle transport problem. It has already been noted [ 19 19 ] that the Shannon entropy concept can be problematic when it indicates the convergence of the fission source distribution in a slow-convergence system. Fortunately, our concerns did not include this type of problem. Analysis for an reasonable adaptive strategy for sample size 2 Two important facts should be considered. (1) The surviving particles of the current step should be transferred to the next step as an external source. The attributes of these particles are samples of unknown distribution that change with the number of steps. This distribution can be called survival particle distribution for convenience; (2) Monte Carlo simulation of particle transport should provide global tallies (one tally per cell) for mechanical calculation as a source item. These facts have a general sense of many multiphysics coupling processes. It would be interesting to develop a sample size adaptive strategy to balance small errors and calculation costs. The basic considerations of this study are as follows: Suppose a sufficiently large sample size of one step is divided into batches, and each batch contains a fixed number of samples, the Monte Carlo simulation should finish all batches one at a time. However, if the number of surviving particles is sufficient, and global tallies have converged according to some judgment rules after some of these batches have been finished, the remaining batches can be omitted because the safely finished batches are sufficient to obtain reasonable results. This point can be reasonably verified by the fact that when stops in advance in one step, the global efficiency of global tallies, according to some index, is flat with batches that have already been finished. To realize these ideas, three questions must be answered. The first question is whether the surviving particles are sufficient. The second one is by which indicator the global tallies can be judged to converge sufficiently. The last one is which on-the-fly diagnostic can be used to stop the Monte Carlo simulation in advance. The solutions to these problems are presented in the following three subsections. These methods are based on the Shannon entropy concept. Indicator for judging whether survival particles are enough 2 Based on the definition of Shannon entropy, an indicator can be designed to show if the number of surviving particles are sufficient. Suppose there are N cells in total for the system and G energy groups so that the total number of phase–space grids is NG . Let <math id="M24"><mrow><msub><mi>V</mi><mi>i</mi></msub><mrow><mo>(</mo><mrow><mi>i</mi><mo>=</mo><mn>1,</mn><mo>⋯</mo><mn>,</mn><mi>N</mi></mrow><mo>)</mo></mrow></mrow></math> be the volume of cell i and m ij be the number of surviving particles in cell i and energy group j . If defining <math display="block" id="M25"><mrow><msub><mi>p</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub><mo>=</mo><mfrac><mrow><msub><mi>m</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub><mo>/</mo><msub><mi>V</mi><mi>i</mi></msub></mrow><mrow><mstyle displaystyle="true"><msub><mo>∑</mo><mrow><mi>i</mi><mi>j</mi></mrow></msub><mrow><msub><mi>m</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub><mo>/</mo><msub><mi>V</mi><mi>i</mi></msub></mrow></mstyle></mrow></mfrac></mrow></math> and <math display="block" id="M26"><mrow><msub><mi>H</mi><mtext>t</mtext></msub><mo>=</mo><mo>−</mo><mstyle displaystyle="true"><munder><mo>∑</mo><mrow><mi>i</mi><mi>j</mi></mrow></munder><mrow><msub><mi>p</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub></mrow></mstyle><mi>ln</mi><mrow><mo>(</mo><mrow><msub><mi>p</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub></mrow><mo>)</mo></mrow><mn>,</mn></mrow></math> then there are L H t values in total for L finished batches. Here, H t is the entropy corresponding to survival particle distribution. Note that p ij ≥ 0 and <math id="M27"><mrow><mstyle displaystyle="true"><msub><mo>∑</mo><mrow><mi>i</mi><mi>j</mi></mrow></msub><mrow><msub><mi>p</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub></mrow></mstyle><mo>=</mo><mn>1</mn></mrow></math> ; thus, p ij can be considered as a discrete probability distribution. When batches are individually simulated, the H t values form a sequence. If the sequence converges, it can be said that these L batches have made the surviving particles sufficient. Apparently, finer phase-space grids represent stricter standards. In Eq. (8), dividing by V i is not trivial because it represents the concept of number density. Note that m ij in batch k should contain m ij in all batches k 1 satisfying k 1 lt; k . Indicator for judging whether global tallies have converged sufficiently 2 When particle transport is coupled with other processes, the Monte Carlo simulation must provide feedback quantities to those processes as an input. The precision of these inputs, also known as tallies in the Monte Carlo method, significantly influence the confidence of the coupling calculation. In most cases, these quantities can be expressed as integrals of the flux and the corresponding functions. Therefore, if global flux tallies converge sufficiently, these feedback quantities are credible to a certain degree. Focusing on global flux tallies can also avoid the problem of dealing with different types of tallies in different multiphysics processes. If F i is the flux in cell i , let <math display="block" id="M28"><mrow><msub><mi>p</mi><mi>i</mi></msub><mo>=</mo><mfrac><mrow><msub><mi>F</mi><mi>i</mi></msub></mrow><mrow><mstyle displaystyle="true"><msubsup><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></msubsup><mrow><msub><mi>F</mi><mi>i</mi></msub></mrow></mstyle></mrow></mfrac></mrow></math> and <math display="block" id="M29"><mrow><msub><mi>H</mi><mtext>f</mtext></msub><mo>=</mo><mo>−</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></munderover><mrow><msub><mi>p</mi><mi>i</mi></msub></mrow></mstyle><mi>ln</mi><mo stretchy="false">(</mo><msub><mi>p</mi><mi>i</mi></msub><mo stretchy="false">)</mo><mn>.</mn></mrow></math> Here, H f is the entropy corresponding to global flux. The convergence of H f sequence formed by finished batches indicates the global convergence of the flux tallies. In addition, F i in batch k should contain the contributions of the particles in all batches k 1 that satisfy k 1 lt; k . To the best of our knowledge, this is the first time an indicator like Eq. (11) has been introduced into the Monte Carlo simulation of particle transport to indicate the convergence of global tallies. This has the advantage that there is no need to calculate the relative error. This advantage saves significant memory and calculation costs, particularly for large-scale models. On-the-fly diagnostic for judging whether entropy sequence has converged 2 Based on the stochastic oscillator indicator, which is commonly used in the technical analysis of financial markets, a posterior diagnostic method has been proposed to assess fission source convergence in Monte Carlo criticality calculations [ 18 18 ]. The stochastic oscillator indicator of the n th iteration step of criticality calculation is defined as follows: <math display="block" id="M30"><mrow><msup><mi>K</mi><mi>n</mi></msup><mo>=</mo><mfrac><mrow><msup><mi>H</mi><mi>n</mi></msup><mo>−</mo><msubsup><mi>H</mi><mrow><mtext>min</mtext></mrow><mrow><mi>n</mi><mn>,</mn><mi>p</mi></mrow></msubsup></mrow><mrow><msubsup><mi>H</mi><mrow><mtext>max</mtext></mrow><mrow><mi>n</mi><mn>,</mn><mi>p</mi></mrow></msubsup><mo>−</mo><msubsup><mi>H</mi><mrow><mtext>min</mtext></mrow><mrow><mi>n</mi><mn>,</mn><mi>p</mi></mrow></msubsup></mrow></mfrac><mn>,</mn></mrow></math> where <math id="M31"><mrow><msubsup><mi>H</mi><mrow><mtext>min</mtext></mrow><mrow><mi>n</mi><mn>,</mn><mi>p</mi></mrow></msubsup></mrow></math> and <math id="M32"><mrow><msubsup><mi>H</mi><mrow><mtext>max</mtext></mrow><mrow><mi>n</mi><mn>,</mn><mi>p</mi></mrow></msubsup></mrow></math> are the minimum and maximum Shannon entropies, respectively, over the last p steps, and H n is the entropy of the current n th step. It has been proven that when the entropy sequence converges, this stochastic oscillator index fluctuates around 0.5 stochastically [ 18 18 ]. The posterior rule for accessing the convergence is given by the following inequality: <math display="block" id="M33"><mrow><mrow><mo>|</mo><mrow><msup><mi>K</mi><mi>n</mi></msup><mo>−</mo><mn>0.5</mn></mrow><mo>|</mo></mrow><mo><</mo><mi>ϵ</mi><mtext> and </mtext><mrow><mo>|</mo><mrow><mfrac><mn>1</mn><mi>m</mi></mfrac><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>0</mn></mrow><mrow><mi>m</mi><mo>−</mo><mn>1</mn></mrow></munderover><mrow><msup><mi>K</mi><mrow><mi>n</mi><mo>+</mo><mi>i</mi></mrow></msup></mrow></mstyle><mo>−</mo><mn>0.5</mn></mrow><mo>|</mo></mrow><mo><</mo><mi>ϵ</mi><mn>,</mn></mrow></math> In other words, this criterion requires that both K n and its average over the next m steps should be within ϵ of 0.5. If the active cycle starts after the iteration step, satisfying the aforementioned inequality, the entire calculation should be accepted. However, this rule can be modified to be an on-the-fly diagnostic [ 20 20 ] for the convergence of the entropy sequence corresponding to the survival particle distribution and global flux tally, as explained in the last two subsections. This point is as follows: If we decompose all samples into S batches and simulate all batches individually, we can obtain a sequence of entropy values and stochastic oscillator indicators called H i <math id="M34"><mrow><mo stretchy="false">(</mo><mi>i</mi><mo>=</mo><mn>1,</mn><mo>⋯</mo><mn>,</mn><mi>S</mi><mo stretchy="false">)</mo></mrow></math> and K i <math id="M35"><mrow><mo stretchy="false">(</mo><mi>i</mi><mo>=</mo><mi>p</mi><mn>,</mn><mo>⋯</mo><mn>,</mn><mi>S</mi><mo stretchy="false">)</mo></mrow></math> . When completing the n th <math id="M36"><mrow><mo stretchy="false">(</mo><mi>n</mi><mo>≥</mo><mi>m</mi><mo>+</mo><mi>p</mi><mo>−</mo><mn>1</mn><mo stretchy="false">)</mo></mrow></math> batch and the following inequality is satisfied: <math display="block" id="M37"><mrow><mrow><mo>|</mo><mrow><msup><mi>K</mi><mrow><mi>n</mi><mo>−</mo><mi>m</mi><mo>+</mo><mn>1</mn></mrow></msup><mo>−</mo><mn>0.5</mn></mrow><mo>|</mo></mrow><mo><</mo><mi>ϵ</mi><mtext> and </mtext><mrow><mo>|</mo><mrow><mfrac><mn>1</mn><mi>m</mi></mfrac><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>0</mn></mrow><mrow><mi>m</mi><mo>−</mo><mn>1</mn></mrow></munderover><mrow><msup><mi>K</mi><mrow><mi>n</mi><mo>−</mo><mi>i</mi></mrow></msup></mrow></mstyle><mo>−</mo><mn>0.5</mn></mrow><mo>|</mo></mrow><mo><</mo><mi>ϵ</mi><mn>,</mn></mrow></math> the surviving particle distribution and global flux tallies can then be perceived to have converged. The parameters p,m were set to 20 and 50, respectively, but the ϵ value was set by the user. The smaller the ϵ is, the stricter the rule is. Therefore, a conservative or radical strategy can be developed by using different ϵ values. Description of the adaptive strategy 2 The entire mechanism for automatically adjusting the sample size is summarized as follows: In each step of the time-dependent Monte Carlo particle transport simulation, a sufficiently large sample size is initially set and divided into M batches. After finishing one batch, new H t and H f values are obtained. When a finished batch number exceeds m + p -1, the on-the-fly diagnostic is invoked to decide if it needs to stop the Monte Carlo simulation in advance according to both the entropy sequences. Only when both judgments return a true value, the simulation of this step is stopped in advance. If stopped in advance in step L , there will not be enough survival particles to be used in step L +1 although the current survival particles can be seen as sufficiently sampled from the survival particle distribution. In this case, we randomly select from the current surviving particle bank to supply more surviving particles as source particles in step L +1. This operation is considered reasonable based on the probability theory. Note that in parallel computing, frequently calculating these entropies should use the method proposed in [ 21 21 ]; otherwise, the time cost of obtaining these entropies will offset the benefit of the sample size adjustment. In other words, this method uses local data in each process to calculate the local entropy and the average of the entropies of all processes as the final entropy value. Thus, massive data broadcasts can be avoided. It can be verified that when the sample size approaches infinity, the entropy is the same as the original entropy. Even with a finite sample size, this new entropy can be used to indicate the convergence of the corresponding survival particle distribution or global tallies. Although there is a statistical estimation method for the Shannon entropy [ 22 22 ], we have not adopted it because of its complexity; however, this will be the focus of our future studies. It is worthwhile to provide four annotations. (1) When an under-sampling problem [ 23 23 ] exists for a specific model, an insufficient history number causes the Shannon entropy corresponding to the survival particle distribution to fluctuate more drastically. Therefore, the inequality Eq. (14) cannot be satisfied easily, and the entire calculation will not be stopped in advance. (2) Theoretically, the method in this study may fail for some symmetric and oscillating models. However, for time-dependent particle transport simulation, strict symmetric and oscillating models are rare; hence, it is not a serious issue. (3) The Fourier fundamental mode coefficient works similar to the traditional Shannon entropy; therefore, the recently proposed Fourier diagnosis [ 24 24 ] can be a substitute for the Shannon entropy diagnosis. Because we adopted an efficient method for calculating Shannon entropy in MPI parallel environments and did not perform a numerical comparison, we do not know which method is better. (4) The recently proposed batch-size growth method [ 25 25 , 26 26 ] is consistent with our method. However, our work aims to improve the efficiency of time-dependent particle transport simulations using multiphysics calculations. The two aforementioned studies aim to improve the efficiency of time-independent criticality calculations. The existing difficulties are different for these two types of problems. Theoretically, the concepts proposed in these studies are suitable for transplantation. These performances can only be compared by future numerical experiments. Numerical results 1 It is worth mentioning again that our goal in designing this sample size adaptive strategy is to avoid unnecessary calculations when the results are sufficiently good when evaluated by a standard. In each step, when the original sample size is divided into several batches, and each batch is simulated individually according to the original plan, we expect that this strategy can stop in advance when the efficiency evaluated by some standard becomes flat. Apparently, because of the inevitable stochastic nature of this strategy, we could not precisely capture the stop point. However, capturing it approximately can still save considerable calculation time while keeping the results approximately unchanged when compared to the original case. The two models discussed in the next sections were used to demonstrate this effect. The Arbitrary Lagrange-Euler (ALE) algorithm was used to perform mechanical calculations. Model one 2 Model one was a sector piece with an open angle of 2 and a radius of 6.45 cm. There are two concentric layers. Each layer is composed of a mixture of U 235 and U 238 but with different densities. In total, 600 cells were used. The initial sample size was 12,800,000 samples per step, which was sufficiently large for this model. These samples were divided into 640 batches with 20,000 samples per batch. The calculations used a total of 32 cores. We perform two independent calculations with ϵ = 0.1 and ϵ = 0.05. All the other conditions were identical to those described above. Figure 1 Figure 1 shows how the neutron number( N p ) in the model changes with time for three cases, with and without sample size adaptation but with two different ϵ values). Table 1 Table 1 shows the comparison of neutron numbers ( N p ) at the end of the calculations for these three cases and the time costs. It can be observed that the proposed sample size adaptive strategy significantly decreases the calculation time while maintaining an almost unchanged result. When the judging rule is stricter, similar to the ϵ =0.05 case, the result is closer to the original case; however, the reduction in the calculation time is smaller. For a more detailed comparison, we randomly chose a step to determine whether the adaptive strategy could stop the calculation in advance when the efficiency flattened. We used the following equation: <math display="block" id="M38"><mrow><mtext>Global_FOM</mtext><mo>=</mo><mfrac><mi>N</mi><mrow><mi>T</mi><mstyle displaystyle="true"><msubsup><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></msubsup><mrow><msubsup><mi>R</mi><mi>i</mi><mn>2</mn></msubsup></mrow></mstyle></mrow></mfrac></mrow></math> to evaluate the efficiency of global flux tallies [ 27 27 ]. Here, N =600 is the total cell number, T is the calculation time, and R i is the relative error of the flux in cell i . In this step, the ϵ =0.1 case stops at the 416 th batch, whereas the ϵ =0.05 case cannot stop in advance; that is, all 640 batches are simulated in this case. From Fig. 2 Fig. 2 , we see that the sample size adaptive strategy stops in advance when the global efficiency is almost flat with batch numbers. Considering the efficiency, when global efficiency becomes flat, it is not necessary to do more simulations. Figure 3 Figure 3 and Fig. 4 Fig. 4 show the two Shannon entropies ( H t and H f ) for the same step. As expected, the sample size adaptive strategy stopped the entire calculation when these entropies flattened. Figure 5 Figure 5 shows the real sample size that was simulated in different cases. The reduction in calculation time was evidently due to a decrease in the sample size. However, the results were almost the same and the efficiency remained almost unchanged because the simulated sample size was already sufficient. We cannot say which is the best ϵ because of inevitable randomness; however, it can be said that a small ϵ means a conservative strategy, whereas a large ϵ means a radical strategy. This is a characteristic of most of the algorithms. However, the results indicate that ϵ can maintain the concerned physical quantities. Therefore, if used appropriately, the method proposed in this study can definitely improve the calculation efficiency and provide a tool for different users. 1 1 Different cases 1 1 N p 1 1 Time cost (hours) 1 1 without sample size adaptive 1 1 2.058× 10 24 1 1 1.75 1 1 sample size adaptive with ϵ =0.1 1 1 2.049× 10 24 1 1 0.93 1 1 sample size adaptive with ϵ =0.05 1 1 2.059× 10 24 1 1 1.73 Model two 2 Model two is a cuboid which contains 200,000 hexahedron cells. It is divided into three layers that contain four types of isotopes of U ( 234 U, 235 U, 236 U, 238 U) but with different densities. The mechanical calculation is done in one step, and the Monte Carlo particle transport simulation is iterated for 100 steps to calculate <math display="block" id="M39"><mrow><mi>λ</mi><mo>=</mo><mfrac><mn>1</mn><mrow><mi>N</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow></mfrac><mfrac><mrow><mtext>d</mtext><mi>N</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow><mrow><mtext>d</mtext><mi>t</mi></mrow></mfrac><mn>.</mn></mrow></math> We performed 20 calculations without sample size adaptation but with different random number seeds to obtain the reference answer. Each calculation used 1,600 batches with 30,000 samples per batch, and 32 cores were utilized. From Table 2 Table 2 , we can see that the sample size adaptive strategy obtains a reasonable λ value while significantly decreasing the calculation time at approximately only 22% of the original case without adaptive sample size adjustment. Note that the number in parentheses is the standard deviation of 20 λ values, and the λ value of the sample size adaptive strategy is located in the 95% confidence interval obtained from these 20 λ values. 1 1 Different cases 1 1 λ 1 1 Time cost (hours) 1 1 without sample size adaptive 1 1 3.7727(0.0238) 1 1 3.04 1 1 sample size adaptive with ϵ =0.1 1 1 3.7709 1 1 0.67 Conclusion 1 Based on the extension of the Shannon entropy concept, two indicators are designed to show whether samples of survival particle distribution are sufficient and whether global flux tallies have converged sufficiently. An on-the-fly diagnostic method based on the stochastic oscillator indicator for the convergence of these two indicators is proposed to avoid manual intervention, thus forming a complete adaptive strategy. This is the first time that it has become possible to adjust the sample size according to certain definite and reasonable rules. Based on the results of the two models displayed in the previous section, this self-adjustment mechanism maintains reasonable precision while reducing the calculation time. This indicator and mechanism can be easily generalized to other multiphysics calculations, including time-dependent Monte Carlo particle transport simulations. References S.H. Du , S.F. Zhang , T.G. Feng et al. , Computer simulation of transport problems , ( Hunan Science and Technology Press , China , 1989 ) Computer simulation of transport problems X. Wang , J.L. Li , Z. Wu et al. , CMGC - A CAD to Monte Carlo geometry conversion code, Nuc . Sci. Tech. 31 , 82 ( 2020 ). doi: 10.1007/s41365-020-00793-8 http://doi.org/10.1007/s41365-020-00793-8 CMGC - A CAD to Monte Carlo geometry conversion code, Nuc L. Deng , G. Li , B.Y. Zhang et al. , A high fidelity general purpose 3-D Monte Carlo particle transport program JMCT3.0 . Nuc. Sci. Tech. 33 , 108 ( 2022 ). doi: 10.1007/s41365-022-01092-0 http://doi.org/10.1007/s41365-022-01092-0 A high fidelity general purpose 3-D Monte Carlo particle transport program JMCT3.0 X-5 Monte Carlo Team , MCNP-A General Monte Carlo N-Particle Transport Code , Version 5 , LA-UR-03-1987 , 2003 MCNP-A General Monte Carlo N-Particle Transport Code F. Balibrea , On Clausius, Boltzmann and Shannon notions of entropy . J. Mod. Phys. 7 , 219 - 227 ( 2016 ). doi: 10.4236/jmp.2016.72022 http://doi.org/10.4236/jmp.2016.72022 On Clausius, Boltzmann and Shannon notions of entropy L. Benguigui , The different paths to entropy . Euro. J. Phys. 34 , 303 - 321 ( 2013 ). doi: 10.1088/0143-0807/34/2/303 http://doi.org/10.1088/0143-0807/34/2/303 The different paths to entropy D. Hammer , A. Romashchenko , A. Shen et al. , Inequalities for Shannon entropy and Kolmogorov complexity . J. Compu. Sys. Sci. 60 , 442 - 464 ( 2000 ). doi: 10.1006/jcss.1999.1677 http://doi.org/10.1006/jcss.1999.1677 Inequalities for Shannon entropy and Kolmogorov complexity C.E. Shannon , A mathematical theory of communication . Mob. Compu. Comm. Rev. 5 3–55,( 2001 ). doi: 10.1145/584091.584093 http://doi.org/10.1145/584091.584093 A mathematical theory of communication P.M. Cincotta , C.M. Giordano , R.A. Silva et al. , The Shannon entropy: An efficient indicator of dynamical stability . Phys. D 417 , 132816 ( 2021 ). doi: 10.1016/j.physd.2020.132816 http://doi.org/10.1016/j.physd.2020.132816 The Shannon entropy: An efficient indicator of dynamical stability R. M. Yulmetyev , N. A. Emelyanova , F. M. Gafarov , Dynamical Shannon entropy and information Tsallis entropy in complex systems . Phys. A 341 , 649 - 676 ( 2004 ). doi: 10.1016/j.physa.2004.03.094 http://doi.org/10.1016/j.physa.2004.03.094 Dynamical Shannon entropy and information Tsallis entropy in complex systems T. Ueki , Stationarity diagnostics with relative entropy and Wilcoxon signed rank initerated-source Monte Carlo methods . Nucl. Sci. Eng. 160 , 242 ( 2008 ). doi: 10.13182/NSE160-242 http://doi.org/10.13182/NSE160-242 Stationarity diagnostics with relative entropy and Wilcoxon signed rank initerated-source Monte Carlo methods T. Ueki , F.B. Brown , Stationarity modeling and informatics-based diagnostics in Monte Carlo criticality calculations . Nucl. Sci. Eng. 148 , 38 ( 2005 ). doi: 10.13182/NSE04-15 http://doi.org/10.13182/NSE04-15 Stationarity modeling and informatics-based diagnostics in Monte Carlo criticality calculations F.B. Brown , On the use of Shannon entropy of the fission distribution for assessing convergence of Monte Carlo criticality calculations . LA-UR-06-3737 ( 2006 ) On the use of Shannon entropy of the fission distribution for assessing convergence of Monte Carlo criticality calculations Z.G. Li , K. Wang , Y.C. Guo et al. , Forced propagation method for Monte Carlo fission source convergence acceleration in the RMC . Nuc. Sci. Tech. 32 , 27 ( 2021 ). doi: 10.1007/s41365-021-00868-0 http://doi.org/10.1007/s41365-021-00868-0 Forced propagation method for Monte Carlo fission source convergence acceleration in the RMC T. Ueki , On-the-fly diagnostics of particle population in iterated-source Monte Carlo Methods . Nucl. Sci. Eng. 158 , 15 - 27 ( 2008 ). doi: 10.13182/NSE08-A2735 http://doi.org/10.13182/NSE08-A2735 On-the-fly diagnostics of particle population in iterated-source Monte Carlo Methods Y. Naito , J. Yang , The sandwich method for determining source convergence in Monte Carlo calculatio . J. Nucl. Sci. Tech. 41 , 559 ( 2004 ). doi: 10.1080/18811248.2004.9715519 http://doi.org/10.1080/18811248.2004.9715519 The sandwich method for determining source convergence in Monte Carlo calculatio I. Kim , H. Kim , Y. Kim , Deterministic truncation of the Monte Carlo transport solution for reactor eigenvalue and pinwise power distribution . Nucl. Sci. Eng. 194 , 14 - 31 ( 2020 ). doi: 10.1080/00295639.2019.1654815 http://doi.org/10.1080/00295639.2019.1654815 Deterministic truncation of the Monte Carlo transport solution for reactor eigenvalue and pinwise power distribution G. Aldrich , S. Dutta , J. Woodring , OpenMC in situ source convergence detection . LA-UR-16-23217 ( 2016 ). doi: 10.2172/1253484 http://doi.org/10.2172/1253484 OpenMC in situ source convergence detection T. Ueki , F.B. Brown , D.K. Parsons et al. , Autocorrelation and dominance ratio in Monte Carlo criticality calculations . Nucl. Sci. Eng. 145 , 279 - 290 ( 2003 ). doi: 10.13182/NSE03-04 http://doi.org/10.13182/NSE03-04 Autocorrelation and dominance ratio in Monte Carlo criticality calculations D.H. ShangGuan , Z.C. Ji , L. Deng et al. , New strategy for global tallying in Monte Carlo criticality calculation . Acta. Phys. Sin. 68 , 122801 ( 2019 ). doi: 10.7498/APS.68.20182276 http://doi.org/10.7498/APS.68.20182276 New strategy for global tallying in Monte Carlo criticality calculation D.H. ShangGuan , L. Deng , B.Y. Zhang et al. , Efficient method of calculating Shannon entropy of non-static transport problem in message passing parallel programming environment . Acta. Phys. Sin. 65 , 142801 ( 2016 ). doi: 10.7498/aps.65.142801 http://doi.org/10.7498/aps.65.142801 Efficient method of calculating Shannon entropy of non-static transport problem in message passing parallel programming environment A. Bulinski , D. Dimitrov , Statistical estimation of the Shannon entropy . Acta. Math. Sin.(English Series) 35 , 17 - 46 ( 2019 ). doi: 10.1007/s10114-018-7440-z http://doi.org/10.1007/s10114-018-7440-z Statistical estimation of the Shannon entropy R.J. Brissenden , A.R. Garlick , Biases in the estimation of K eff and its error by Monte Carlo methods . Ann. Nucl. Energy 13 , 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 M.R. Omar , J.A. Karim , Fission source stationarity diagnostics using the Fourier fundmental mode coefficient . Prog. Nucl. Energy (English Series). 146 , 104164 ( 2022 ). doi: 10.1016/j.pnucene.2022.104164 http://doi.org/10.1016/j.pnucene.2022.104164 Fission source stationarity diagnostics using the Fourier fundmental mode coefficient Q.Q. Pan , N. An , T.F. Zhang et al. , Single-step Monte carlo criticality algorithm . Comp.Phys. Comm. 279 , 108439 ( 2022 ). doi: 10.1016/j.cpc.2022.108439 http://doi.org/10.1016/j.cpc.2022.108439 Single-step Monte carlo criticality algorithm Q.Q. Pan , T.F. Zhang , X.J. Liu et al. , Optimal batch size growth for wielandt method and superhistory method, Nucl . Sci. Eng. 196 , 183 - 192 ( 2022 ). doi: 10.1080/00295639.2021.1968223 http://doi.org/10.1080/00295639.2021.1968223 Optimal batch size growth for wielandt method and superhistory method, Nucl Q.Q. Pan , Redevelopment of shielding module and research on advanced variance reduction methods based on RMC code. Ph.D thesis, Tsinghua University ( 2020 ) Redevelopment of shielding module and research on advanced variance reduction methods based on RMC code. All authors contributed to the study conception and design. Material preparation, data collection and analysis were performed by Dan-Hua ShangGuan,Wei-Hua Yan, Jun-Xia Wei, Zhi-Ming Gao, Yi-Bing Chen, and Zhi-Cheng Ji. The first draft of the manuscript was written by Dan-Hua ShangGuan and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript. 10.1007/s41365-023-01202-6 This work was supported by the CAEP Found (No. CX20200028) and Youth Program of National Natural Science Foundation of China (No. 11705011). 2022-12-03 2023-01-29 2023-02-14 2023-04-25 2023-06-07 © 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. 2023 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-023-01202-6 http://dx.doi.org/10.1007/s41365-023-01202-6 . Springer Nature or its licensor (e.g. a society or other partner) holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law. 2023 38408216 The data that support the findings of this study are openly available in Science Data Bank at https://www.doi.org/10.57760/sciencedb.j00186.00061 and http://resolve.pid21.cn/31253.11.sciencedb.j00186.00061. Data Availability Statement Nuclear Science and Techniques 4 34