Research on a Monte Carlo global variance reduction method based on an automatic importance sampling method

NUCLEAR ENERGY SCIENCE AND ENGINEERING

Research on a Monte Carlo global variance reduction method based on an automatic importance sampling method

Yi-Sheng Hao，

Zhen Wu，

Shen-Shen Gao，

Rui Qiu ，

Hui Zhang，

Jun-Li Li

Nuclear Science and Techniques

Vol.35, No.5

Article number 86

Published in print May 2024

Available online 31 May 2024

DOI：10.1007/s41365-024-01404-6

96005

Global variance reduction is a bottleneck in Monte-Carlo shielding calculations. The global variance reduction problem requires that the statistical error of the entire space be uniform. This study proposed a grid-AIS method for the global variance reduction problem based on the AIS method, which was implemented in the Monte Carlo program MCShield. The proposed method was validated using the VENUS-III international benchmark problem and a self-shielding calculation example. The results from the VENUS-III benchmark problem showed that the grid-AIS method achieved a significant reduction in the variance of the statistical errors of the MESH grids, decreasing from 1.08E-02 to 3.84E-03, representing a 64.00% reduction. This demonstrates that the Grid-AIS method is effective in addressing global issues. The results of the self-shielding calculation demonstrate that the grid-AIS method produced accurate computational results. Moreover, the grid-AIS method exhibited a computational efficiency approximately one order of magnitude higher than that of the AIS method and approximately two orders of magnitude higher than that of the conventional Monte Carlo method.

Monte CarloGlobal variance reductionReactor shieldingAutomatic importance sampling

Introduction

Global variance reduction is a classification of shielded computational problems. Global problems require not only the statistics of full-space fluxes or flux-related responses but also low statistical errors for all spaces. The global variance reduction problem is common in the design of new reactors, such as the Chinese Fusion Engineering Testing Reactor (CFETR) [1] and the solid fuel thorium molten salt reactor (TMSR-SF1) [2]. Currently, the geometric structures of the new reactors are relatively complex. When designing the shielding, the statistical error of all positions in the entire space must be small and relatively uniform to ensure that the design requirements for radiation protection are satisfied. If only the radiation parameters of the individual locations are calculated, the location of the radioactive leakage can be ignored. Therefore, there are global issues regarding the shielding design of new reactors. In the calculation space, the magnitude of the flux of the particles at different positions differs. Therefore, a corresponding variance reduction method for the global problem must be adopted to reduce the statistical error of the position with low flux. Finally, a more average statistical error distribution in the whole space is obtained, and simultaneously, it can save computing resources.

In the Monte Carlo calculation of a global problem, the MESH virtual grid is typically set to count the calculation results everywhere in the entire space. Therefore, the goal of Monte Carlo calculations of the global problem can be defined as: for a given relative statistical error $H$ , calculation time $T$ must be low [3], and simulated particle numbers must be few such that the relative statistical error of more grids $R \leq H$ , and the statistical error variance of all grids should be minimal.

In conventional Monte Carlo transport, there are more particles in the region near the source, and the statistical error is small. In other words, the number of particles in the region far from the source is less owing to collisions and other reasons, and the statistical error is large. Although simply increasing the number of simulation particles can reduce the overall variance, the statistical error distribution is not uniform, and it is a waste of resources for particles to continue to be transported near the source, which does not meet the requirements of the global problem. Therefore, a corresponding global variance-reduction technique must be developed.

Principle and Current Situation of Global Variance Reduction

In a study on global variance reduction, Cooper and Larsen et al. [4] noted that a uniform Monte Carlo simulation of particle density can provide a uniform relative statistical error distribution, thus satisfying the goal of global variance reduction. Only particles arriving in the statistical region in the Monte Carlo simulations can contribute to the statistics. The number of Monte Carlo simulation particles in the statistical region determines the magnitude of the relative statistical error. This idea can also be explained by the statistical error calculation Eq. (1) for Monte Carlo simulations. Considering the statistical penetration probability $p$ of the plate shield and assuming that the number of particles arriving in the region behind the plate is $N_{1}$ , the unbiased estimate $\hat{p}$ of $p$ is $\hat{p} = \frac{N_{1}}{N} .$ (1)

Relative statistical error is obtained by substituting the statistical error calculation formula: $\frac{x \hat{σ}}{\hat{p} \sqrt{N}} \approx \frac{x \sqrt{\frac{N_{1}}{N}}}{\frac{N_{1}}{N} \sqrt{N}} = \frac{x}{\sqrt{N_{1}}} .$ (2)

The relative statistical error in this region is entirely determined by the number of particles $N_{1}$ entering the region. If the Monte Carlo simulation particle density is the same in the full space, the relative statistical error is also the same. In an actual Monte Carlo simulation, a numerical analysis of the relative statistical error is more difficult because the statistics are more complex and may use geometric splitting with a roulette, weight window, and other variance reduction techniques.

The various global variance reduction techniques that have emerged are also based on this conclusion, using different methods to achieve uniform Monte Carlo-simulated particle density distributions, including the coupling variance-reduction method based on the forward calculation represented by Cooper and Larsen [4] and the coupling variance reduction technique based on adjoint represented by FW-CADIS. The core idea of this method is to first perform one or two auxiliary calculations to generate a reasonable weight window parameter and realize a uniform Monte Carlo simulation of the particle density distribution.

Cooper and Larsen [4] guided the weight window parameter settings based on deterministic forward transport fluxes. For each mesh $i$ , assuming that $N_{i}$ represents the physical particle density of mesh $i$ , $V_{i}$ represents the mesh volume, $M_{i}$ represents the number of Monte Carlo simulated particles of mesh $i$ , and $w_{i}$ represents the average particle weight of mesh $i$ ; then, for each mesh i: $N_{i} \propto \frac{w_{i} M_{i}}{V_{i}} .$ (3)

Suppose that the central weight of the window or the survival weight $w_{i} \propto N_{i}$ , then: $\frac{M_{i}}{V_{i}} \approx c .$ (4) where $c$ is a constant. In other words, setting the average weight of the weight window to be proportional to the physical particle density can make the Monte Carlo simulation particle density approximately uniform. The flux is $ϕ_{i} \propto w_{i} v_{i}$ , where $v_{i}$ represents the particle velocity. In practice, Cooper and Larsen [4] chose the forward flux $ϕ_{i}$ to set the weight window parameters as follows: $w_{i} \approx \frac{ϕ_{i}}{\max (ϕ_{i})} .$ (5)

In 2007, Wagner proposed the FW-CADIS method based on CADIS for the global variance reduction of CADIS [5]. The computational goal of the global variance reduction problem is to obtain the calculation results of the global uniform convergence. For the Monte Carlo transport calculations, this can also be understood as obtaining a globally uniform Monte Carlo particle distribution. According to the conjugate transport theory, a conjugate source item $q^{+}$ can be defined such that the conjugate flux calculated by conjugate transport based on the source item can represent the importance of particles to achieve a globally uniform Monte Carlo particle distribution. In the FW-CADIS method, the conjugate source item $q^{+}$ is defined as $q^{+} (r, E) = \frac{σ_{d} (r, E)}{\int ψ (r, E) σ_{d} (r, E) d E} .$ (6)

The equation shows that the FW-CADIS method applies the reciprocal of the total corresponding quantity estimate to weigh the conjugate source item. Therefore, the farther away from the physical source item, the greater the enhancement of the conjugate source item, and the conjugate source item increases accordingly. Additionally, the conjugate calculation process of the FW-CADIS method is consistent with that of the CADIS method. The difference between the two methods is that FW-CADIS requires an additional forward deterministic transport calculation before the conjugate calculation to provide the required weighted source term for the conjugate calculation. Therefore, the calculation effect of the FW-CADIS method depends on the accuracy of the conjugate flux provided by the deterministic program.

Other scholars conducted in-depth research on the global variance-reduction problem. Van Wijk [6] proposed a method for setting the weight window threshold according to the flux and statistical error. This method is entirely based on the I/O function of the MCNP Monte Carlo radiation transfer code family. Yuan [7] improved it by defining a novel "pseudo-source" based on the flux-based variance reduction method and proposed the PS-GVR method. To improve the efficiency of the source analysis, Qingquan [8] proposed a single-step Monte Carlo criticality method that eliminated the calculation of inactive cycles and achieved the highest computational efficiency through a mathematical optimization analysis. This method was used for the numerical analysis of transplutonium isotope production [9]. Meanwhile, the DeGVR method [10] was proposed to improve the efficiency of the shielding calculations. This method adopts the strategy of “density reduction + density extrapolation,” which can obtain global information quickly and is helpful for global variance reduction. Similar methods include PDMC [11], adaptive variance reduction [12] and SP3-coupled methods [13].

For a single probe, the FOM [14] is often used to measure computational efficiency. However, for global problems, because many detectors are used when obtaining MESH statistics, the formula for calculating the FOM factors is no longer applicable. This is a commonly used method to demonstrate the computational efficiency of counting the relative statistical error of each mesh or the probability distribution of the FOM factor and displaying it in the form of two-dimensional images. However, it is still significant to provide a quantitative number for characterizing computational efficiency. Therefore, researchers have attempted to improve the FOM factor formula to satisfy the requirements of global problems. In this study, the following formula is used to measure the computational efficiency of the global variance reduction method: $FOM = \frac{1}{{\bar{R}}^{2} T} .$ (7) where $T$ is the calculation time and $\bar{R}$ is the average relative statistical error of each mesh result. The calculation formula is as follows: $\bar{R} = \frac{Σ R_{i}}{N_{m}},$ (8) where $R_{i}$ represents the relative statistical error of mesh $i$ and $N_{m}$ represents the MESH number.

Limitations of the AIS Method

To solve the deep-penetration problem, Jiajin proposed the AIS method [15]. This method is based on importance sampling and statistical estimation; it introduces a virtual surface, divides the space into multilayer subspaces, generates virtual particles on the virtual surface to be transported to the next layer of the subspace, and performs automatic particle weight adjustment and quantity control. As shown in Fig. 1, the AIS method is as follows.

Fig. 1

AIS method transportation flow chart

(1) $K$ virtual surfaces are introduced, and the transport space is divided into $K + 1$ subspaces;

(2) Virtual particles are generated on the virtual surface at the source event and at each collision event, and the number of virtual particles on the virtual surface is equaled to the number of source particles.

(3) The particles are killed if they reach the current virtual surface during transport.

(4) The virtual particles on the current virtual surface continue to be transported as source particles to the next subspace.

The AIS method uses a statistical estimation to generate virtual particles that continue to be transported as the source term of the subspace in the next layer. In this manner, more particles can be transported to a region far away from the source term, which solves the problem of large statistical errors in the region far away from the source term to a certain extent.

Suppose that the positive particle flow, $J_{+}$ on a surface refers to the number of particles passing through the face from the positive surface. In Monte Carlo simulations, the concept of particle weights is introduced, and $J_{+}$ is calculated by summing the weights of all the particles passing through the surface from a positive surface. Because the number of virtual particles on each virtual surface is constant at $N$ and the weights are all the average weight $\bar{w}$ , the total positive particle flow on the virtual surface is $J_{+} = \frac{\bar{w} N}{N} = \bar{w} .$ (9)

The particle current $Δ J_{+}$ in the infinitesimal element on the virtual surface is: $Δ J_{+} = \frac{\bar{w} Δ N}{N} .$ (10)

Substituted into formula (8), then: $Δ N = N \frac{Δ J_{+}}{\bar{w}} = N \frac{Δ J_{+}}{J_{+}} .$ (11)

In summary, the number of virtual particles on an infinitesimal element is equal to the total number of virtual particles on the virtual surface multiplied by the ratio of the size of the forward particle flow on the infinitesimal element to the size of the total forward particle flow on the virtual surface. Therefore, the number of virtual particles was higher when the level of forward particle flow on the virtual surface was high; thus, the statistical error was smaller. Conversely, at low levels of forward particle flow, the number of virtual particles is lower; thus, the statistical error is larger. The level of the forward particle flow is determined by the problem, independent of the detector location. If the detector is in a region where the horizontal forward particle flow is small, the statistical error of the detector response will be large, and the AIS method will not have a good effect on variance reduction. Therefore, it can be concluded that the AIS method has certain limitations when applied to global problems.

In this study, the grid-AIS method is proposed for the global problem, and the Monte Carlo simulation particle density distribution is made uniform by constructing a corresponding virtual particle number control algorithm.

Grid-AIS Method

4.1

Transport Process

The key to achieving global variance reduction is that the Monte Carlo simulation particle density distribution should be uniform. A virtual surface divides an entire space into several subspaces. If the virtual particle density distribution on each virtual plane is uniform, then the Monte Carlo simulation particle density distribution in the entire space is approximately uniform. The grid-AIS method differs from the AIS method. The AIS method only keeps the number of virtual particles on each virtual surface the same, while the grid-AIS method keeps splitting according to the ratio $α = w / {\bar{w}}_{i}$ of the weight $w$ of virtual particles to the average weight ${\bar{w}}_{i}$ of each grid, and keeps the number of particles in a grid constant, thus ensuring the uniform distribution of the number of particles on the virtual surface.

Based on the analysis, the grid-AIS method achieved a uniform virtual particle density distribution on the virtual surface using a meshing virtual particle number control algorithm. Small meshes of approximately equal areas were evenly divided on the virtual surface, and the number of virtual particles generated by each mesh was maintained. The transport process is as follows.

(1) In the initialization stage, $K$ virtual surfaces were introduced, and the transport space is divided into $K + 1$ subspaces. For each virtual surface $K$ , the virtual surface is divided into $M_{k}$ uniform meshes.

(2) During the source-term sampling and each collision event, the particle produces a virtual particle on the virtual surface. Mesh $i$ to which the virtual particle belongs was determined based on the positions of the virtual particles. For virtual face $K$ , the number of virtual particles in mesh $i$ is maintained as $N_{k, i}$ .

(3) If a particle reaches the current virtual plane during transport, the source particle is eliminated.

(4) The virtual particle on the current virtual surface is returned as the source particle of the next subspace, and (2) continues to be transported.

For virtual surface $K$ , mesh $i$ maintains the corresponding number of particles $N_{k, i}$ . To ensure uniformity in the number of virtual particles, $N_{k, i} = c S_{k, i},$ (12) where $c$ is a constant independent of the mesh and virtual face and $S_{k, i}$ is the mesh area. Thus, the number of particles maintained by each mesh $i$ is proportional to its area. Summing over all the meshes on virtual surface $K$ yields the total number of particles $N_{k}$ on each virtual surface, which is also proportional to its area $S_{k}$ , as follows: $N_{k} = c S_{k} .$ (13)

This differs from the rule in which the number of virtual particles on each virtual surface is the same as that in the traditional AIS method. The traditional AIS method considers the total flux of the space between virtual surfaces and requires only the same number of particles between the virtual surfaces to ensure that the statistical error of each virtual surface is as similar as possible. Global variance reduction must ensure that the statistical error of each mesh flux is as consistent as possible; therefore, the number of virtual particles on the virtual face must be adjusted according to the size of the virtual face.

This distinction is evident in the calculation of item-source global problems. The virtual polygon is generally divided into a series of spheres centered on the point source with radii ranging from small to large. The number of virtual particles in each layer in the traditional AIS method was the same. Although the total number of particles remains the same when the particles are transported to outer space, the number of particles per unit volume is insufficient owing to the increase in the area of the virtual surface. The virtual surface area of the grid-AIS method gradually increased from the inside to the outside, and the number of virtual particles also increased. When particles are transported to outer space, the number of particles per unit volume is consistent with that of the inner layer, which increases the statistical error.

When meshing, this paper requires that the mesh area of each virtual polygon is approximately the same, so if $S_{k, i} \approx S$ is constant, then the number of virtual particles per mesh, $N_{k, i}$ , is also roughly constant. At this point, the number of meshes $M_{k}$ on virtual surface $K$ has the following relation: $M_{k} \approx \frac{S_{k}}{S},$ (14) where $S_{k}$ is the area of each virtual surface and $S$ is the area of each mesh. In other words, the number of meshes on a virtual surface is approximately proportional to its area.

4.2

Meshing and Virtual Particle Adjustment

In the grid-AIS method, virtual particles must be generated on the virtual plane during source item sampling and each collision event, and mesh $i$ is determined according to the position of the virtual particles. First, the mesh division must specify the mesh size parameter $d$ , which determines the granularity of the mesh division; the general value of $d$ is the same as the mesh size of the MESH statistics.

The commonly used planes, cylindrical surfaces, and spherical surfaces must be strictly meshed by area. Figure 2 shows the meshing of the cylindrical virtual surfaces.

Fig. 2

(Color online) Schematic diagram of cylindrical virtual surface meshing

(1) The mesh division of the virtual plane surface was relatively simple. Suppose that its length in the $x$ -direction is $L_{x}$ and that in the $y$ -direction is $L_{y}$ ; then, the $x$ -direction is divided into $r o u n d (L_{x} / d)$ meshes of equal length, and the $y$ -direction is divided into $r o u n d (L_{y} / d)$ meshes of equal length. where $r o u n d (x)$ is $x$ rounded to an integer to ensure that the area of each mesh is identical.

(2) The division of the cylindrical surface was similar to that of the plane. Divide evenly in the two directions of the azimuth $φ$ and the height $Z$ of the cylinder.

(3) For the spherical case, if the mesh is divided uniformly in two directions according to the polar angle $θ$ and azimuthal angle $φ$ , it will cause the mesh area to be inconsistent. If the mesh is divided uniformly according to $cos θ$ and azimuthal angle $φ$ , although it can ensure a consistent mesh area, it will cause too much difference in the mesh shape, and a long strip mesh will be generated at the $θ = 0$ position. Thus, within the same mesh, the two ends of the long stripes are far apart, the flux-level difference may be large, and the virtual particles are not uniformly distributed, causing an uneven statistical error distribution that does not meet the requirement of global variance reduction.

In summary, this paper adopts the way that the polar angle $θ$ is uniformly divided while the azimuthal angle $φ$ is nonuniformly divided. Assume that the radius of the sphere is $r$ , the difference between the maximum and minimum of the spherical polar angle is $Δ θ$ , and the difference between the maximum and minimum of the azimuthal angle is $Δ φ$ . Then, the number of meshes in the $θ$ direction is: $n_{θ} = round (\frac{r Δ θ}{d}),$ (15)

For the polar angle mesh $i$ , the $n_{i, φ}$ azimuthal angles $φ$ are further subdivided to ensure that the area of each small mesh is approximately the same. Considering that the spherical area of the polar angle mesh $i$ is $S_{i} = r^{2} (\cos θ_{i} = \cos θ_{i + 1}) Δ φ .$ (16)

Then, $n_{i, φ} = round (\frac{r^{2} (\cos θ_{i} - \cos θ_{i + 1}) Δ φ}{d^{2}}),$ (17) where $θ_{i}$ and $θ_{i + 1}$ are the starting and ending polar angles of the polar angle mesh $i$ , respectively.

Thus, the area of each mesh was approximately the same and the shape of the mesh was closer to a square, which further reduced the flux difference within the mesh.

When mesh $i$ to which the virtual particle belongs is determined according to the position coordinates of the virtual particle on the virtual surface, considering that the virtual surface may have undergone translation, rotation and other operations, the global coordinates of the virtual particle must be converted into local coordinates with the virtual surface as the reference system, and then the virtual surface mesh number must be determined. Consequently, the determination method within local coordinates is relatively simple. If the aforementioned mesh division method is used, only a simple interpolation is required, which is not repeated here.

After determining the mesh $i$ to which the virtual particles belong, roulette and splitting is performed according to the ratio $α = w / \bar{w_{i}}$ of the weight $w$ of the virtual particles to the average weight ${\bar{w}}_{i}$ of each mesh, maintaining the corresponding number of particles $N_{k, i}$ .

For other cases of nonplanar, cylindrical, and spherical surfaces, the modeling of virtual surfaces may include intersection, concatenation, and complementary cases, which are more complex and make it difficult to divide the mesh strictly into equal areas. In this study, an approximate division method was adopted, and statistical MESH was used to divide the virtual surface. The MESH number where the virtual particle is located is used as the virtual face mesh number, and the number of virtual particles per virtual face mesh is maintained in proportion to the volume of the MESH. Here, the volume of the MESH is considered to be approximately proportional to the area of the virtual surface contained in it.

Table 1 shows the comparison between the AIS and grid-AIS methods.

Comparison of AIS and grid-AIS methods

Method	Virtual particle generation	Virtual particle number control	Variance reduction parameter
AIS Method	Generate virtual particles on the virtual surface along the direction of particle motion during source sampling and each collision event	Maintain an equal number of particles on each virtual surface	Only need to set virtual surface position parameters
Grid-AIS Method		Maintain an equal number of virtual particles per mesh on each virtual surface	Need to set virtual surface position parameters and grid size parameters

Currently, the grid-AIS method is implemented in the MCShield program. MCShield is a neutron/photon/electron-coupled transport Monte Carlo program for radiation shielding calculations independently developed by the Radiation Protection and Environmental Protection Laboratory of Tsinghua University. The MCShield program can simulate neutron, photon, and electron coupled transport with massively parallel computation functions and can efficiently and accurately solve the deep penetration problem and the complex shielding problem of the Monte Carlo variance reduction technique. At the same time, the software also contains a variety of powerful pre- processing and post-processing modules, such as CAD geometry conversion, parametric geometry modeling, visualization parameter setting, particle trajectory display and three-dimensional dose display. The program was validated against dozens of international benchmarks [21].

In summary, compared with other methods, the grid-AIS method proposed in this paper for global problems has the following main innovations. (1) For the global variance reduction problem, the weight window method is currently widely used. The weight window method usually divides the grid in the entire space and sets the weight window parameters of each grid. Furthermore, the effect of variance subtraction computation was highly correlated with the quality of the weight window parameters. In contrast to the weight window method, the grid-AIS method must only set a few variance reduction parameters, such as the virtual surface position and grid size, which is convenient for users. (2) In the CADIS method, a deterministic method must be used for the adjoint calculation in a single calculation to calculate the parameters of the weight window. Therefore, the entire calculation process requires two sets of geometric models, Monte Carlo programs, and determinism, which makes the shielding calculation very difficult. The grid-AIS method does not require adjoint calculations and only requires one Monte Carlo forward calculation to obtain the result, which is convenient for users.

4.3

Unbiased Description

The general transport architecture of the grid-AIS method was the same as that of the AIS method, whereas the unbiased nature of the AIS method was demonstrated and is not repeated here. The number of virtual particles in the virtual surface mesh was controlled using roulette and splitting, and the process was unbiased. The amount of particle transport between virtual surfaces varies with the virtual surface area, which is different from the AIS method. However, this process is achieved by gambling splitting, the weights are adjusted accordingly, and the overall process is unbiased. Overall, the Grid-AIS method was unbiased.

VENUS-III International Benchmark Problem

5.1

Description of Benchmark

The VENUS-III [22] experimental reactor is a low-flux thermal neutron reactor that was built in 1964 to study reactor core designs and the irradiation damage of nuclear materials. This provides benchmark experiments for numerical simulation programs. The 1/4 core model of the VENUS-III reactor is shown in Fig. 3, and it consists of 639 fuel rods and 11 control rods, with dimensions of 73.65 cm length, 37.80 cm width, and 37.80 cm height. In the figure, the light-grey area represents the control rods, whereas the other parts are fuel rods. The grey parts represent fuel rods enriched to 4% burnable poison, the white parts represent fuel rods enriched to 3% burnable poison, and the dark grey parts represent fuel rods enriched to 3% burnable poison.

Fig. 3

(Color online) (a) Top view and (b) side view of the VENUS-III experimental reactor benchmark

The VENUS-III reactor plane diagram is shown in Fig. 3. It is divided into nine horizontal regions from the center to the left as follows: (1) the central water channel outside the core; (2) the inner layer baffle outside the core (stainless steel thickness of 2.858 cm); (3) the core fuel area (4% fuel rods and 3.3% fuel rods); (4) the outer layer baffle outside the core (stainless steel thickness of 2.858 cm); (5) the water reflector outside the core (minimum water thickness of 2.169 cm); (6) the outer core basket (stainless steel thickness of 4.99 cm); (7) the cooling water layer outside the basket (water thickness of 5.80 cm); (8) the thermal shield (stainless steel thickness of 6.72 cm); and (9) the reactor pressure vessel (stainless steel). From the center, it was divided into eight regions, from bottom to top, as follows: (1) the lower water region inside the pressure vessel, (2) the lower core support structure, (3) the outer bottom support frame of the core, (4) the lower internal reflector of the core, (5) the core fuel area, (6) the upper internal reflector of the core, (7) the outer upper support frame of the core, and (8) the upper water region inside the pressure vessel.

5.2

Calculation Results

The source term in the core region consists of fission neutrons with a Watt spectrum and isotropic distribution. The source term is distributed within the cylindrical fuel rods, and the sampling probabilities for different positions within each fuel rod vary. In this case, the neutron flux in the computational region was calculated using MESH tallies. The MESH grid size is 20.00 mm × 20.00 mm × 20.00 mm, with 46, 46, and 85 MESH grids in the X, Y, and Z directions, respectively. In this study, three methods were compared: AIS, grid-AIS, and traditional Monte Carlo methods. Nine cylindrical virtual surfaces were set up for the grid-AIS and AIS methods. The center of the virtual surface of the cylinder moved down from the center of the core (0, 0, 0) to (0, 0, -400 mm), and the radius of the virtual surface of the cylinder increased from 500 to 900 mm. Virtual particles were generated on the top, bottom, and lateral surfaces of the cylinder. In the Grid-AIS method, grid size parameter d=40, 80, and 160 mm. Table 2 shows the calculation results of the VENUS-III international benchmark.

Calculation Results of the VENUS-III International Benchmark

Method	NPS	CPU Time (min)	Mean absolute error	MESH average statistical error (%)
Conventional Monte Carlo	1.00 ×10⁸	460	6.29×10^-8	24.00
Conventional Monte Carlo	1.00×10⁷	42	2.50×10^-7	36.89
AIS	1.00×10⁷	330	1.41×10^-7	7.58
Grid-AIS	1.00×10⁷	352	1.42×10^-7	6.80

Table 3 and Fig. 4 show that both the grid-AIS and AIS methods showed significant improvements in computational results, compared to the traditional Monte Carlo algorithm. The average statistical error of the MESH grid decreased from 0.24 to 0.07, demonstrating the effectiveness of the AIS method. Compared to the AIS method, the Grid-AIS method showed the following improvements in the computational results: the percentage of MESH grids with a statistical error greater than 0.80 decreased from 0.31% to 0.02%; the percentage of MESH grids with a statistical error greater than 0.50 decreased from 1.56% to 0.33%; the percentage of MESH grids with a statistical error less than 0.30 increased from 95.98% to 98.40%; and the percentage of MESH grids with a statistical error less than 0.10 increased from 83.14% to 87.58%. Furthermore, the variance in the statistical errors for the MESH grids in the grid-AIS method decreased from 1.08% to 0.38%, representing a 64.00% reduction. This indicates that the grid-AIS method reduced the global variance problem. In addition, by comparing the calculation results for different grid size parameters, we found that the calculation results for d=80 and 160 mm were better than those for d=40 mm. Generally, the grid size parameter should match the size of the MESH grid and the number of simulated particles. A grid-size parameter that is either too large or too small affects the effectiveness of the grid-AIS method.

Statistical Results of the VENUS-III International Benchmark

Method	Variance of MESH statistical error (%)	Proportion of grids with statistical error > 30%	Proportion of grids with statistical error < 10%
Conventional Monte Carlo	4.36	30.67%	27.43%
Conventional Monte Carlo	5.98	62.55%	20.65%
AIS	1.08	4.02%	83.14%
Grid-AIS	0.38	1.60%	87.58%

Fig. 4

(Color online) Statistical error result of (a) the AIS method (Y= -900 mm), (b) grid-AIS method (Y= -900 mm, d = 80), (c) AIS Method (Z = 0 mm), (d) grid-AIS method (Z = 0 mm, d = 80), (e) grid-AIS method (Z = 0 mm, d = 40), and (f) grid-AIS method (Z = 0 mm, d = 160) results of the VENUS-III benchmark

In this section, the effectiveness of the grid-AIS method is validated using the VENUS-III international benchmark problem. However, in this benchmark problem, the radius of the computational region was relatively small compared to the core radius, resulting in minimal differences in the surface areas of the different cylindrical virtual surfaces. In addition, because of the relatively uniform distribution of geometry and materials in the VENUS-III international benchmark questions, the calculation results of the AIS method were already good, with little room for improvement. Therefore, based on an actual reactor structure, this study proposes a self-designed reactor shielding example in which pipes and other structures are added to increase the geometric anisotropy, and the grid-AIS method is further verified.

Self-Design Reactor Shielding Example

6.1

Example Description

In this section, the grid-AIS method is validated using a self-designed reactor-shielding algorithm. The geometric model is shown in Fig. 5. The core is a cylinder with a radius of 0.50 m and height of 1 m, and its lower surface is flush with the bottom of the water layer. Outside the core is a water layer with a radius of 1.50 m and height of 4.00 m. Outside of the water layer is a concrete layer with a thickness of 0.50 m, and outside the vacuum boundary. Radial and axial air pipes with a radius of 0.10 m pass through the water layer and concrete. The source was uniformly distributed in the core area, the energy distribution was a Watt fission spectrum, and the direction was isotropic. The calculation example has a thick shielding body, pipes, and other structures with strong anisotropy, which is a typical shielding calculation problem with difficult convergence of global counting, and can be used to test the computational capability of the Monte Carlo method under a more anisotropic geometry.

Fig. 5

(Color online) Three views of the geometry of the self-design reactor shielding example

The subbody fluxes in full space were counted using MESH with a mesh size of 4.40 cm × 4.40 cm × 4.40 cm and a total of 568,800 meshes. Calculations were performed separately using the original MCShield, AIS, and Grid-AIS methods. The AIS virtual surfaces were 18 spherical surfaces with the center of the core as the center of the sphere, from the inside to the outside, and the radii of the spherical surfaces were spaced 20.00 cm apart. Owing to the limitation of the overall geometric boundary, the spheres are not closed spheres, and the starting and ending polar angles are determined by the spheres and vacuum geometric boundary.

The computing platform used in this study is a multinode computing cluster with two 18-core processors per node: Intel Xeon E5-2699 V3 @2.3 GHz CPU with 64 GB of memory.

The calculation results are listed in Tables 4~5. To facilitate comparison of the statistical errors, the number of simulated particles was controlled for each method. The simulation times of the methods are relatively close; therefore, the magnitude of the statistical error represents the computational efficiency. Using 100 processes for parallel computation, the CPU time listed in Table 1 is the sum of the running times of each process. The computational efficiency was measured using the FOM factor applied to a global problem.

Simulation parameters for the self-design reactor shielding example

Method	NPS	CPU time (min)	Mean absolute error	Mean relative statistical error (%)	FOM (min^-1)
Conventional Monte Carlo	1.00×10¹⁰	6.54×10⁴	8.17×10^-10	62.20	3.95×10^-5
AIS	6.00×10⁸	7.37×10⁴	1.15×10^-9	11.90	9.56×10^-4
Grid-AIS	1.00×10⁸	7.81×10⁴	2.78×10^-9	3.30	1.15×10^-2

Calculation Results of the Self-Design Reactor Shielding Example

Method	Variance of MESH statistical error (%)	Proportion of grids with statistical error > 30%	Proportion of grids with statistical error < 10%
Conventional Monte Carlo	7.13	33.99%	45.09%
AIS	3.80	14.81%	69.28%
Grid-AIS	1.67	5.12%	86.95%

It can be seen from the parameter results that the average relative statistical error of the conventional Monte Carlo method was as high as 62.20%, and the results were unreliable. The grid-AIS method has an average statistical error of 3.30%, which is less than 5.00% and is acceptable for shielding calculations. In addition, the FOM factor of the grid-AIS method was approximately 12 times higher than that of the AIS method, and 290 times higher than that of the conventional Monte Carlo method.

6.2

Example Results

Sectional views of the statistical MESH flux are shown in Figs. 6, 7, and 8. The top of each figure represents the X-Z plane, which is the plane through the vertical air duct. The figure below shows the YZ-plane perpendicular to the horizontal air duct. The left and right images of the two-dimensional images in Figs. 6, 7, and 8 have the same meaning, and the flux is displayed in logarithmic coordinates. The colors represent the fluxes. Red indicates the highest, blue indicates the lowest, and the span between each contour is one order of magnitude.

Fig. 6

(Color online) Conventional Monte Carlo flux cross-section: (a) front view and (b) side view

Fig. 7

(Color online) AIS flux cross-section: (a) front view and (b) side view

Fig. 8

(Color online) Grid-AIS flux cross-section: (a) front view and (b) side view

The maximum flux attenuation was approximately 10^-19. The flux decreased exponentially with the distance from the core in the range near the core. In the upper area, far away from the core, the flux was higher in the air duct and near the air duct because of the influence of the two air ducts, and the flux was lower in other areas. The conventional Monte Carlo method has many empty meshes at positions far away from the source and only has values near the air pipe, which indicates that the conventional Monte Carlo method cannot transport particles to positions far away from the core. By contrast, the AIS and grid-AIS methods essentially have no empty meshes.

The relative statistical errors of the MESH fluxes are shown in the cross-sectional plots in Figs. 9, 10, and 11, and the relative statistical errors are displayed in linear coordinates. For display purposes, the width of each color was 2.50%, from 0.00 to 20.00%, and colors above 20.00% were red, for a total of nine colors. As shown in the illustration, the conventional Monte Carlo method yields poor results in the upper part of the space away from the core, where the statistical error is lower than 20.00% only near the vertical air duct and higher than 20.00% everywhere else. The AIS method has a smaller statistical error near the air ducts in the upper part of the space away from the core, and a larger statistical error at a location away from the ducts. Combined with the flux cut diagram, the statistical error changes in size with the same pattern of change in the flux size. The statistical error was smaller when the flux level was high on the virtual surface, and vice versa. This is consistent with the results of a previous analysis conducted using the AIS method. As shown in Fig. 11, the statistical errors of the Grid-AIS method were smaller and more uniformly distributed than those of the AIS method, and the statistical errors were also smaller at locations far from the air ducts. This indicates the strategy of the grid-AIS method to ensure that the virtual particles are uniformly distributed on the virtual surface.

Fig. 9

(Color online) Conventional Monte Carlo statistical error cross-section: (a) front view and (b) side view

Fig. 10

(Color online) AIS statistical error cross-section: (a) front view and (b) side view

Fig. 11

(Color online) Grid-AIS statistical error cross-section: (a) front view and (b) side view

To verify the correctness of the grid-AIS method, the relative deviation results of the grid-AIS method and results of the conventional Monte Carlo and AIS methods are plotted separately, as shown in Figs. 12-13; the colors used in the plots have the same meaning as in the statistical error plots. The figure shows that the relative deviation of the Grid-AIS method results from the conventional Monte Carlo and AIS method results near the core and inside the vertical air ducts is low, below 2.50%, which further verifies the correctness of the Grid-AIS method. By contrast, regions with higher relative deviations also have higher statistical errors and are therefore not meaningful for comparison.

Fig. 12

(Color online) Sectional diagram of flux relative deviation between conventional Monte Carlo and Grid-AIS: (a) front view and (b) side view

Fig. 13

(Color online) Sectional diagram of flux relative deviation between AIS and Grid-AIS (a) front view and (b) side view

The cumulative probability distribution functions with frequency histograms of the statistical errors for each MESH for the three methods are shown in Fig. 14. The horizontal coordinates represent the relative statistical error, which ranged from 0.00% to 100.00%. The vertical coordinate is the cumulative probability distribution function of the relative error, representing the proportion of MESH below a certain relative statistical error relative to the total mesh. As shown in Fig. 14a, the statistical error of the grid-AIS method increases most steeply, and the meshes with statistical errors below 10.00% account for approximately 94.00%. The AIS method, on the other hand, has approximately 77.00% of the meshes with statistical errors below 10.00%. The conventional Monte method was even lower, accounting for approximately 20.00%. In addition, the statistical errors of the conventional Monte Carlo and AIS methods show a step increase of approximately 100.00% because these two methods have partially empty meshes.

Fig. 14

(a) Cumulative probability distribution function and (b) frequency histogram of MESH statistics results

Figure 14b presents the frequency histograms of the relative statistical errors for each MESH. Relative statistical errors were divided into 1.00% intervals, and the number of MESHs within each relative statistical error segment was counted as a percentage of the total number of MESH. To make the graph more concise, the vertical boundaries of each statistical error segment in the histogram are hidden. As shown in Fig. 14b, the statistical errors of the Grid-AIS method were relatively concentrated, and the proportion of higher statistical errors was much lower than those of the conventional Monte Carlo and AIS methods. The spatial distribution of the Monte Carlo simulation particle density in the grid-AIS method is more uniform; therefore, the statistical error is more uniform in the whole space and is reflected in the frequency histogram, and the statistical error distribution shows a trend of being more concentrated.

In addition, it can be found from the data that the statistical error distribution of the conventional Monte Carlo and AIS methods has a more consistent trend in the areas with larger relative statistical errors. Although the AIS method reduces the overall statistical error by introducing a virtual surface, it does not change the Monte Carlo simulation particle density distribution on the virtual surface; therefore, the error change trend is close to that of the conventional Monte Carlo method. A relative statistical error of 100.00% represents no statistical value, that is, an empty mesh. In the figure, the vertical coordinate at a 100.00% relative statistical error represent the ratio of the number of empty meshes to the total number of meshes. This ratio was approximately 50.00% for the conventional Monte Carlo method, approximately 2.00% for the AIS, and 0.00% for the grid-AIS, representing no empty mesh.

Conclusion

Based on the basic AIS method, a Grid-AIS method applicable to global variance reduction is proposed in this paper and implemented in the Monte Carlo program MCShield.

In this study, the grid-AIS method was validated using the VENUS-III international benchmark problem and a self-shielding reactor calculation. The results demonstrated that the grid-AIS method produced accurate computational results for addressing the global variance reduction problem, thereby confirming the correctness of the grid-based AIS method. In both validation cases, there was a reduction of approximately 60.00% in the variance of the statistical errors of the computed results. The grid-AIS method also exhibited a computational efficiency approximately one order of magnitude higher than that of the AIS method and approximately two orders of magnitude higher than that of the conventional Monte Carlo method. These results validate the effectiveness of the grid-AIS method.

In summary, the Grid-AIS method can better solve the global problem in shielding calculations and is an extension and supplement to the AIS method.

References

X. C. Nie, J. Li, S. L. Liu et al.,

Global variance reduction method for global Monte Carlo particle transport simulations of CFETR. Nucl

. Sci. Tech. 28, 115 (2017). https://doi.org/10.1007/s41365-017-0270-3

Baidu Scholar

Google Scholar

P. Yang, Y. Dai, Y. Zou et al.,

Application of global variance reduction method to calculate a high-resolution fast neutron flux distribution for TMSR-SF1

. Nucl. Sci. Tech. 30, 125 (2019). https://doi.org/10.1007/s41365-019-0650-y

Baidu Scholar

Google Scholar

C. B. Kiedrowski, Ibrahim et al., Evaluating the efficiency of estimating numerous Monte Carlo tallies. 2011 ANS Annual Meeting (2011).

M. A. Cooper, E. W. Larsen,

Automated weight windows for global Monte Carlo particle transport calculations

. Nucl. Sci. Eng. 137, 1-13 (2001). https://doi.org/10.13182/NSE00-34

Baidu Scholar

Google Scholar

J. C. Wagner, D. E. Peplow, S. W. Mosher,

FW-CADIS method for global and regional variance reduction of Monte Carlo radiation transport calculations

. Nucl. Sci. Eng. 176, 37-57 (2014). https://doi.org/10.13182/NSE12-33

Baidu Scholar

Google Scholar

A.J. van Wijk, G. Van den Eynde, J.E. Hoogenboom,

An easy to implement global variance reduction procedure for MCNP

. Ann. Nucl. Energy 38, 2496-2503 (2011) https://doi.org/10.1016/j.anucene.2011.07.037

Baidu Scholar

Google Scholar

Y. Hu, S. Yan, Y. F. Qiu et al.,

Implementation and benchmarking of an automatic global variance reduction method on OpenMC

. Fusion Eng. Des. 173, 112829 (2021). https://doi.org/10.1016/j.fusengdes.2021.112829

Baidu Scholar

Google Scholar

Q.Q. Pan, N. An, T.F. Zhang et al.,

Single-step Monte Carlo criticality algorithm

. Comput. Phys. Commun. 279, 108439 (2022). https://doi.org/10.1016/j.cpc.2022.108439

Baidu Scholar

Google Scholar

Q.Q. Pan, Q.F. Zhao, L.J. Wang et al.,

Rapid diagnosis method for transplutonium isotopes production in high flux reactor

. Nucl. Sci. Tech. 34, 44 (2023). https://doi.org/10.57760/sciencedb.j00186.00039

Baidu Scholar

Google Scholar

10.

Q.Q. Pan, L.J. Wang, Y. Cai et al.,

Density-extrapolation global Variance Reduction (DeGVR) method for large-scale radiation field calculation

. Comput. Math. Appl. 143, 10-22 (2023). https://doi.org/10.1016/j.camwa.2023.04.024

Baidu Scholar

Google Scholar

11.

Q.Q. Pan, H.W. Lv, S.Q. Tang et al.,

Pointing probability driven semi-analytic Monte Carlo Method (PDMC) - Part I: Global Variance Reduction for large-scale radiation transport analysis

. Comput. Phys. Commun. 291, 108850 (2023). https://doi.org/10.1016/j.cpc.2023.108850

Baidu Scholar

Google Scholar

12.

Q.Q. Pan, J.J. Rao, S.F. Huang et al.,

Improved adaptive variance reduction algorithm based on RMC code for deep penetration problems

. Ann. Nucl. Energy 137, 107113 (2020). https://doi.org/10.1016/j.anucene.2019.107113

Baidu Scholar

Google Scholar

13.

Q.Q. Pan, T.F. Zhang, X.J. Liu et al.,

SP3-coupled global variance reduction method based on RMC code

. Nucl. Sci. Tech. 32, 122 (2021). https://doi.org/10.1007/s41365-021-00973-0

Baidu Scholar

Google Scholar

14.

A.C. Olivieri, G.M. Escandar,

Analytical figures of merit

. Introduction to Multivariate Calibration: A Practical Approach, 159-177 (2018). https://doi.org/10.1016/b978-0-12-410408-2.00006-5

Baidu Scholar

Google Scholar

15.

J. J. Fan,

Research on Monte Carlo method for solving problems with small probability and big contribution

. Ph.D. thesis, Tsinghua University (2004) (in Chinese)

Baidu Scholar

Google Scholar

16.

C.Y. Li,

Research on Monte Carlo method for solving point detector flux and deep penetration problem

. Ph.D. thesis, Tsinghua University (2008) (in Chinese)

Baidu Scholar

Google Scholar

17.

X. Wang,

Research on key methods and program development of Monte Carlo calculation for high-efficiency radiation shielding

. Ph.D. thesis, Tsinghua University (2016) (in Chinese)

Baidu Scholar

Google Scholar

18.

S. S. Gao,

Research on Variance Reduction method for Monte Carlo simulation of reactor Shielding

. Ph.D. thesis, Tsinghua University (2019) (in Chinese)

Baidu Scholar

Google Scholar

19.

X. Wang, J.L. Li, Z. Wu et al.,

Improved algorithms and coupled neutron-photon transport for auto-importance sampling method

. Chin. Phys. C 41, 014103 (2017). https://doi.org/10.1088/1674-1137/41/1/014103

Baidu Scholar

Google Scholar

20.

Y.S. Hao, R. Qiu, Z. Wu et al.,

Research on the source-detector variance reduction method based on the AIS adjoint Monte Carlo method

. Ann. Nucl. Energy 191, 109916 (2023). https://doi.org/10.1016/j.anucene.2023.109916

Baidu Scholar

Google Scholar

21.

S.S. Gao, Z. Wu, X. Wang et al.,

Development of a radiation shielding Monte Carlo Code: RShieldMC

. Presented at International Conference on Mathematics and Computational Methods Applied to Nuclear Science and Engineering 2017, Jeju, Korea, April 16-20 (2017).

Baidu Scholar

Google Scholar

22.

R. Rulko, Computing radiation dose to reactor pressure vessel and internals. Italy: the ENEA-Bologna Nuclear Data Centre, 1997.

Footnote

The authors declare that they have no competing interests.