Research basics

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-7]. Therefore, there are various definitions of entropy. In this study, the entropy concept used is the one proposed by Shannon in [8], which has some novel applications in various fields [9, 10]. For example, Ref. [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-14]. Suppose there are N cells in total for some nuclear systems and M fission neutrons at some moment. Let mi be the number of fission neutrons in cell i. Hence, ∑imi=M. Let $p_i=\frac{m_i}{M}$ and $H_{\text{ori}}=-{\displaystyle \sum _{i=1}^M{p}_i}\ln (p_i)\mathrm{,}$ (7) we obtain the Shannon entropy H_ori. For the power iteration method of eigenvalue calculation, if there are K iteration steps and $H_{\text{ori}}^k\mathrm{,}k=\mathrm{1,}\cdots \mathrm{,}K$ 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-18] for the convergence of the fission source distribution, which belongs to the time-independent particle transport problem. It has already been noted [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

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

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 $V_i\left(i=\mathrm{1,}\cdots \mathrm{,}N\right)$ be the volume of cell i and mij be the number of surviving particles in cell i and energy group j. If defining $p_{ij}=\frac{m_{ij}/V_i}{{\displaystyle {\sum }_{ij}{m}_{ij}/V_i}}$ (8) and $H_{\text{t}}=-{\displaystyle \sum _{ij}{p}_{ij}}\ln \left(p_{ij}\right)\mathrm{,}$ (9) then there are L Ht values in total for L finished batches. Here, H_t is the entropy corresponding to survival particle distribution. Note that pij ≥ 0 and ${\displaystyle {\sum }_{ij}{p}_{ij}}=1$; thus, pij 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 Vi is not trivial because it represents the concept of number density. Note that mij in batch k should contain mij in all batches k₁ satisfying k₁ lt;k.

Indicator for judging whether global tallies have converged sufficiently

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 Fi is the flux in cell i, let $p_i=\frac{F_i}{{\displaystyle {\sum }_{i=1}^N{F}_i}}$ (10) and $H_{\text{f}}=-{\displaystyle \sum _{i=1}^N{p}_i}\ln (p_i)\mathrm{.}$ (11) 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, Fi in batch k should contain the contributions of the particles in all batches k₁ that satisfy k₁ 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

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]. The stochastic oscillator indicator of the n^th iteration step of criticality calculation is defined as follows: $K^n=\frac{H^n-H_{\text{min}}^{n\mathrm{,}p}}{H_{\text{max}}^{n\mathrm{,}p}-H_{\text{min}}^{n\mathrm{,}p}}\mathrm{,}$ (12) where $H_{\text{min}}^{n\mathrm{,}p}$ and $H_{\text{max}}^{n\mathrm{,}p}$ are the minimum and maximum Shannon entropies, respectively, over the last p steps, and Hn 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]. The posterior rule for accessing the convergence is given by the following inequality: $\left|K^n-0.5\right|<ϵ\text{ and }\left|\frac{1}{m}{\displaystyle \sum _{i=0}^{m-1}{K}^{n+i}}-0.5\right|<ϵ\mathrm{,}$ (13) In other words, this criterion requires that both Kn 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] 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 Hi $(i=\mathrm{1,}\cdots \mathrm{,}S)$ and Ki $(i=p\mathrm{,}\cdots \mathrm{,}S)$. When completing the n^th $(n\ge m+p-1)$ batch and the following inequality is satisfied: $\left|K^{n-m+1}-0.5\right|<ϵ\text{ and }\left|\frac{1}{m}{\displaystyle \sum _{i=0}^{m-1}{K}^{n-i}}-0.5\right|<ϵ\mathrm{,}$ (14) 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

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]; 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], 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] 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] 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, 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.

Model one

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²³⁵ and U²³⁸ 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 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 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: $\text{Global\_FOM}=\frac{N}{T{\displaystyle {\sum }_{i=1}^N{R}_i^2}}$ (15) to evaluate the efficiency of global flux tallies [27]. Here, N=600 is the total cell number, T is the calculation time, and Ri 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, 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 and 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 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.

Fig. 1Full-screen View

Comparison of neutron number

Fig. 2Full-screen View

Global efficiency of all flux tallies in some step

Fig. 3Full-screen View

Shannon entropy H_t in some step

Fig. 4Full-screen View

Shannon entropy H_f in some step

Fig. 5Full-screen View

The real sample size in different cases

Model two

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 (²³⁴U, ²³⁵U, ²³⁶U, ²³⁸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 $\lambda =\frac{1}{N(t)}\frac{\text{d}N(t)}{\text{d}t}\mathrm{.}$ (16) 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, 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.