NUCL SCI TECH (2017)28:115 DOI 10.1007/s41365-017-0270-3

Global variance reduction method for global Monte Carlo particle transport simulations of CFETR Xing-Chen Nie1 • Jia Li1 • Song-Lin Liu2 • Xiao-Kang Zhang2 • Ping-Hui Zhao1 Min-You Ye1 • German Vogel1 • Xiao Yang1 • Qing-Jun Zhu2

•

Received: 23 June 2016 / Revised: 14 November 2016 / Accepted: 16 November 2016 Ó Shanghai Institute of Applied Physics, Chinese Academy of Sciences, Chinese Nuclear Society, Science Press China and Springer Nature Singapore Pte Ltd. 2017

Abstract It can be difficult to calculate some under-sampled regions in global Monte Carlo radiation transport calculations. The global variance reduction (GVR) method is a useful solution to the problem of variance reduction everywhere in a phase space. In this research, a GVR procedure was developed and applied to the Chinese Fusion Engineering Testing Reactor (CFETR). A cylindrical CFETR model was utilized for comparing various implementations of the GVR method to find the optimum. It was found that the flux-based GVR method could ensure more reliable statistical results, achieving an efficiency being 7.43 times that of the analog case. A mesh tally of the scalar neutron flux was chosen for the GVR method to simulate global neutron transport in the CFETR model. Particles distributed uniformly in the system were sampled adequately through ten iterations of GVR weight window. All voxels were scored, and the average relative error was 2.4% in the ultimate step of the GVR iteration. Keywords Global variance reduction  Weight window  Monte Carlo  MCNP  Neutronics

This work was supported by the National Special Project for Magnetic Confined Nuclear Fusion Energy (Nos. 2013GB108004 and 2015GB108002), and the Chinese National Natural Science Foundation (No. 11175207). & Jia Li [email protected] 1

School of Nuclear Science and Technology, University of Science and Technology of China, Hefei 230027, Anhui, China

2

Institute of Plasma Physics, Chinese Academy of Sciences, Hefei 230031, Anhui, China

1 Introduction There is a growing need for performing high-resolution radiation transport and shutdown dose analyses for fusion reactors [1, 2]. The results of global Monte Carlo (M-C) simulations are statistically acceptable across the entire domain covered. Within the MCNP code [3], the commonly used localized variance reduction (LVR) method [4–9] is a weight window generator, which aims to increase computational efficiency and reduce variance for a localized tally. Nevertheless, the use of a weight window generator in a large and composite target cell implies an obvious risk. The particles that reach one part of the cell much more easily than the other parts will dominate the generated weight windows, and the other parts may not be sampled adequately [2]. Consequently, the LVR can hardly handle the variance reduction problem for a global system. The global variance reduction (GVR) method is an efficient solution to global Monte Carlo particle transport problems. Based on weight windows that uniformly transport non-analog particles throughout the model, which allows all areas of the system to be adequately sampled, GVR improves computational efficiency of hypothetical global detectors and reduces the global variance. Cooper and Larsen [10], who laid theoretical base of the GVR method, utilized a diffusion solution of the forward flux to improve global Monte Carlo particle transport calculations. Davis and Turner [11, 12], and Naish et al. [13], applied this kind of GVR method to specific fission and fusion calculation problems. van Wijk et al. [14] used a GVR procedure to a fission case. In this paper, we estimate the forward flux to establish the weight windows and adopt iteration method in the M-C simulation. A GVR procedure is developed for the global Monte Carlo particle transport

123

115

Page 2 of 7

X.-C. Nie et al.

problems of Chinese Fusion Engineering Test Reactor (CFETR).

2 GVR method For large and complex fusion devices such as CFETR, the dense components of shielding material may well attenuate the neutron population density, causing relatively rare events and poor Monte Carlo estimates in some areas. In order to control the sampling frequency, one can use mesh-based weight windows (WW) to avoid large fluctuations in the Monte Carlo particle weights which usually produce large variances in interesting regions. Setting the parameters of the weight windows in a specified manner allows one to effectively determine which regions should be sampled more or less frequently. Based on the ideas of Ref. [10], we used an M-C estimation of the scalar flux / as a way to set the WW thresholds as Eq. (1) C1 /i WWi ¼ ¼ C2 IMPi Maxð/Þ

ð1Þ

The lower WW threshold (WWi ) is set inversely proportional to the importance (IMPi ) of a certain mesh element i. C1 and C2 are constants defined by different applications. Equation (1) ensures that neutrons are more likely to be rouletted in mesh elements of low importance (high flux), hence the increase in computational efficiency. It also allows particles to be split in the mesh elements of high importance (low flux), so that a considerable number of daughter neutrons lead to better statistic estimates. As a result, the global region can be sampled uniformly. The weights of particles may vary greatly from one region to another, which prevents the distant regions from being underpopulated, and the close-to-source regions from being over-populated. Figure 1 shows the flowchart of GVR method based on neutron flux. Global tallies (GVR tallies) are set over the whole model before calculation. An initial analog is run to obtain the neutron distribution data and the WW file format. The iteration goes on when the sampling is insufficient in the global system. A WW value will be generated according to Eq. (1), if the neutron flux result is qualified (the relative error is less than a certain threshold). Otherwise, the WW bound will be set to a certain constant (0.001 is chosen as the constant). Then, the value of WW lower bound is processed by a code that gives play to the functions of data reordering and replacing. A new GVR WW is valuable for the next iteration. The neutron distribution data may be poor in highly attenuated regions, but it can provide a basis for further iterations. Each iteration improves the calculation accuracy. The iteration ends when

123

Fig. 1 Flowchart of GVR method based on neutron flux

enough neutrons spread through the whole model and the simulation obtains credible global results.

3 Tests and discussion 3.1 Comparison of GVR methods for cylindrical model A cylindrical neutronics model of CFETR, with the top and bottom planes being refectory boundaries, was employed to assess the efficiency of various forms of GVR techniques (Fig. 2, R = 964.5 cm, H = 2000 cm, B1 = 214.4 cm, B2 = 234.5 cm). An initial analog run

Fig. 2 Cylindrical CFETR model (in cm)

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

was performed to obtain neutron distribution data, which was utilized to generate different cell-based GVR WW. Most kinds of distribution data were inversely related to the cell importance. Given this relation, different types of data were selected from a wider range of possible candidates. For example, Davis and Turner [11, 12] compared the flux, weight and population. van Wijk et al. [14] used relative error as one of the particle data. In the present work, neutron energy and track were also selected for GVR comparisons, as thoroughly decrease when neutrons penetrate deeper. They are important among the reference data for Monte Carlo transport. For estimating GVR calculation efficiency, we introduced a global figure of merit FOMG from Ref. [11] to represent the efficiency degree of the GVR WW. It is P defined as FOMG = N/(Tcpu e2i ), (i = 1 ? N), where ei is the error in the ith cell, Tcpu is the CPU time in minutes and N is the number of total mesh elements. The spatial distribution of relative errors throughout the mesh should be as flat as possible. The standard relative error P P rerr = [( e2i )/N - ( ei)2/N2]1/2 (i = 1 ? N) should remain as small as possible to indicate the uniformity of the samples. Based on various neutron distribution data, only one step of iteration was taken under conditions of the same calculation time (One iteration is enough for the cylindrical model). Table 1 shows that the best method to generate a global importance map is the neutron flux of the highest FOMG (54.32). Regarding global variance reduction efforts, the flux-based GVR method ensures more statistically reliable results with the lowest average error (0.108), the highest number of scoring voxels (100%) and a lower rerr (0.16). The flux-based GVR method achieved an efficiency of 7.43 times higher than that of the analog method. Figure 3 shows the analog results yielded with an obvious uncertainty in the area far away from the first wall, especially in the outboard blanket. In comparison, the flux-based GVR method

Table 1 Comparison of different distribution data Distribution data

FOMG

Average error

rerr

Scoring (%)

Flux

54.32

0.108

0.16

100

Weight

51.04

0.109

0.17

91.03

Energy Population

44.69 41.16

0.154 0.110

0.14 0.19

100 100

Track

39.51

0.108

0.2

97.44

Error

19.02

0.162

0.28

93.59

7.31

0.302

0.43

74.36

Analog

Page 3 of 7 115

provides more valuable flux data distributed across the whole model. 3.2 Applications to a three-dimensional CFETR model By comparing various forms of GVR methods, a mesh tally of the scalar neutron flux was chosen for the GVR method to simulate global neutron flux of a 3-DCFETR model (water-cooled ceramics breeding blanket) [16–18] of 22.5° (Fig. 4). It contained all the representative inboard and outboard components: the plasma chamber, breeding blanket, shield blanket, divertor, vacuum vessel, thermal shield, magnets, and ports. Radial sizes of the blanket (breeding plus shielding) were 0.7 and 1.2 m in the inboard and outboard zones, respectively. Three ports were filled with shielding materials of water (in volume ratio of 70%) and steel (in volume ratio of 30%). Thus, global deep penetration is a typical problem in this model. For calculating the neutron flux and GVR WW, a mesh tally of 50 divisions along the X and Y directions and 80 divisions along the Z direction was utilized to cover the entire reactor model. The voxel dimension was about 10 cm 9 25 cm 9 25 cm. For this global radiation simulation, the WW lower bound was normalized to 0.5 in the source region. The mesh-based GVR WW was divided into the thermal and fast neutrons ranges. Instead of using all flux information, flux results with relative errors greater than a certain criterion were filtered out, because previous simulations indicated significant statistic uncertainties fora large change in flux gradient. Therefore, a region with a poor flux estimate region could lead neutrons to overly split due to a sharp decrease of the WW bounds, and result in a less effective simulation. With a too low threshold of error, less flux values would be used for generating GVR WW, hence an insufficient variance reduction for the global problem. In this simulation, the tolerance of flux relative error was 0.5. For generating a GVR WW, a stable procedure was compiled by C?? language. This procedure had the two functions: (1) reordering the mesh tally data according to the right order of the WW file; and (2) replacing previous lower bound WW with a new bound GVR WW. Figure 5 shows neutron flux distributions and relative error maps at different GVR WW iteration stages. Without GVR WW, the neutron flux produced by the analog method occupied only a part of the blanket, and most of the areas presented poor neutron flux results or no recorded results in blank voxels. This under-sampling downgraded the response to the global detectors placed far away from the plasma source. Although a vast number of starting particles

123

115

Page 4 of 7

X.-C. Nie et al.

Fig. 3 (Color online) Neutron flux distribution in the outboard blanket (a) and inboard blanket (b) of a cylindrical CFETR model

Fig. 4 (Color online) Neutronics model of CFETR

were tried, the dense component of the reactor material greatly attenuated the neutrons. In the subsequent iterations with the GVR WW set by the scalar flux, the distribution of particles in the system became more extensive and uniform. The flux spread across a larger region of the model, resulting in a larger range of spatial WW, which facilitated the GVR function in the next iteration. We note that neutrons have little possibility of being transported into the three ports filled with the shielding materials of water and steel. That deep penetration problem was solved gradually, when the steps of the GVR WW iteration were increased. Suitable iteration time and steps, which should be decided by the type of model and user experience, could improve the efficiency of the integrated iteration calculation. In Table 2, the rerr decreases with the iteration times and reaches 0.09 in the final iteration, which means that the GVR method is capable of flattening the relative error distribution of a global mesh tally. As expected, the average error decreases with increasing iteration times, while the scoring voxels increases. The average relative error

123

changes from 78.8% under the analog method (Step 0) to 5.7% in the eighth step. The percentage of scoring voxels increases from 27.79 to 100%. In theory, the average time per history should increase as the WW covers more areas in the latter iterations. Nonetheless, Table 2 shows that this value fluctuates through this series of iterations, which indicates that some special histories have an abnormal lifespan. Also, the CTMs are greater than corresponding CTMEs, since the MCNP time cutoff would be delayed to the next rendezvous to assure consistent data. At Iteration 4, the CTM is about three times that of the CTME. Regardless of how we adjusted the CTME (which was below 2907.95 min), the CTM was constantly 2907.95 min. This suggests that certain neutron history may be very long (even longer than 1907.95 min) in this iteration. This phenomenon is called as a long history problem [15], which can drastically shorten the computer time when running MCNP in parallel mode. In this simulation, the MCNP was adapted to parallel processing with 24 processors. As a result, the computer time was much longer than that for CTME, and the parallel efficiency [15] (defined as the transport CPU time in minutes divided by the wall clock time and the number of processors) dropped in the later steps of the iteration containing more GVR WWs. In this case, the 3-D CFETR model was utilized for preliminary design, though it was not detailed enough. Fewer device gaps may make the fewer high statistical weight particles stream to the regions where a large amount of WW splitting would take place. Bulk shielding and void gap are main reasons for the long history problem. When the weights of the particles have a magnitude higher than the lower GVR WW threshold, excessive splitting would appear near the area [15]. Hence, the long history problem is not so detrimental in the whole iteration process. Although some computer time was significantly extended

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

Page 5 of 7 115

Fig. 5 (Color online) Neutron flux (top) and relative error (bottom) maps obtained using a normalized neutron source

Table 2 Summary of GVR weight window iteration (time in min) Iterations

CTME

CTM

Computer time

Av. time per history

Parallel efficiency (%)

rerr

Average error

Scoring (%)

0

1000

1002

1480

3.24 9 10-5

67.66

0.38

0.788

27.79

1

1000

1030

1232

3.32 9 10-4

83.57

0.42

0.732

42.77

-4

2

1000

1040

1259

4.95 9 10

82.54

0.42

0.657

56.87

3

1000

1026

3633

1.14 9 10-3

28.24

0.34

0.390

88.56

4

1000

2908

19,159

2.90 9 10-2

15.18

0.24

0.324

95.57

5

5000

5069

27,503

1.27 9 10-2

18.43

0.22

0.215

97.62

6

10,000

12,178

74,467

2.03 9 10-2

16.35

0.19

0.189

99.02

7

20,000

20,666

138,240

2.95 9 10-2

14.95

0.17

0.132

99.77

207,360

-3

19.54

0.09

0.057

100.00

8

40,000

40,517

9.42 9 10

CTME, computer time spent on the transport portion of the problem, which is set before the simulation CTM, the maximum amount of computer time spent on the Monte Carlo calculation, which is shown in the output file

in the later iterations, the entire iteration process was completed within an acceptable time range. Figure 6 shows the fraction of mesh voxels under a certain error threshold as a function of the relative error. The more iterations are simulated, the larger the fraction of meshes that present an acceptable error in the results.

The lines move rapidly to the top left, where greater fractions of mesh voxels are plotted at small statistical errors. As expected, the analog simulation (AN) performed the poorest, and only 16.3% of the total voxels had an error below 10%. Neutrons can scarcely be transported beyond the shielding material and filled

123

115

Page 6 of 7

Fig. 6 (Color online) Cumulative distribution of relative errors for iteration 0–8

X.-C. Nie et al.

ports; while 88.7% of the total voxels had an error under 10% in the 8th step of adequately sampled iterations. So, the GVR method can improve the variance reduction for the global model. In Fig. 7, the GVR and analog results are compared at the final step of iteration (the 10th iteration) in the same computer time (10,000 min). The flux and distribution of analog Monte Carlo particles vary by many orders of magnitude far from the source, with large statistical errors in regions far from the source. Using the GVR method, Monte Carlo particles are uniformly distributed throughout the system. All meshes are tallied (scoring = 100%). The relative error in the scalar flux remains considerably flat (rerr = 0.05). In an average relative error of 0.024, few errors in each mesh are [0.05.

Fig. 7 (Color online) Simulation results by the analog and GVR methods in 10,000 min of computer time at the final iteration (the 10th), using a normalized neutron source

123

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

4 Conclusion In typical global neutron transport simulation of CFETR, one could rarely obtain precise results. As a solution, a global variance reduction procedure was demonstrated and applied. For a 1-D cylindrical model, various GVR methods were compared. The results show that scalar flux-based GVR is the most efficient one among them. It can achieve an efficiency of 7.43 times higher than the analog simulations. Applying the GVR method to a 3D model, a good global importance map was obtained through several steps of GVR WW iteration, and an acceptable smoothness of errors across a large space was achieved. In the final WW iteration, all meshes scored an average error of 2.4%. The mitigation of the long history problem and the application of GVR to more complex models are still under ongoing research.

8.

9.

10.

11.

12.

13.

References 14. 1. A. Davis, Radiation shielding of fusion systems. Ph.D. thesis, University of Birmingham (2010) 2. F. Wasastjerna, Using MCNP for fusion neutronics. Ph.D. thesis, Helsinki University of Technology (2008) 3. X-5 Monte Carlo Team, A general Monte Carlo N-particle transport code version 5. Los Alamos National Laboratory, LACP-03-0245 (2003) 4. T.E. Booth, J.S. Hendricks, Importance estimation in forward Monte Carlo calculations. Fus. Sci. Technol. 59, 90–100 (1984) 5. T.E. Booth, Genesis of the weight window and the weight window generator in MCNP—a personal history. Los Alamos National Laboratory, LA-UR-06-5807 (2006) 6. S. Sato, H. Iida, T. Nishitanie, Evaluation of shutdown gammaray dose rates around the duct penetration by three-dimensional Monte Carlo decay gamma-ray transport calculation with variance reduction method. J. Nucl. Sci. Technol. 39, 1237–1246 (2002). doi:10.1080/18811248.2002.9715316 7. F. Wasastjerna, Application of importances or weight windows in MCNP4C to a geometry similar to an ITER equatorial port. Ann.

15.

16.

17.

18.

Page 7 of 7 115

Nucl. Energy 35, 425–431 (2008). doi:10.1016/j.anucene.2007. 07.011 B.E. Thomas, MCNP variance reduction examples. Los Alamos National Laboratory, LA-UR-12-25907 (2004). doi:10.2172/ 1054246 J.E. Eggleston, M.A. Abdou, M.Z. Youssef, The use of MNCP for neutronics calculations within large buildings of fusion facilities. Fus. Eng. Des. 42, 281–288 (1998). doi:10.1016/S09203796(98)00215-4 M.A. Cooper, E.W. Larsen, Automated weight windows for global Monte Carlo particle transport calculations. Nucl. Sci. Eng. 137, 1–13 (2001). doi:10.13182/nse00-34 A. Davis, A. Turner, Comparison of global variance reduction techniques for Monte Carlo radiation transport simulations of ITER. Fus. Eng. Des. 86, 2698–2700 (2011). doi:10.1016/j. fusengdes.2011.01.059 A. Davis, A. Turner, Application of novel global variance reduction methods to fusion radiation transport problems, in International Conference on Mathematics and Computational Methods Applied to Nuclear Science and Engineering, Rio de Janeiro, RJ, Brazil, May 8–12 (2011). http://www.ccfe.ac.uk/ assets/Documents/INTERN1.pdf J. Naish, F. Fox, Z. Ghani et al., Radiation mapping at JET and ITER using advanced computational acceleration techniques and tools. Nucl. Technol. 192, 299–307 (2015). doi:10.13182/nt14132 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). doi:10.1016/j.anucene.2011. 07.037 A. Davis, A. Turner, Improving computation efficiency of MonteCarlo simulation with variance reduction, in International Conference on Mathematics and Computational Methods Applied to Nuclear Science and Engineering, Sun Valley, Idaho, USA, May 5–9, (2013). https://arxiv.org/abs/1309.6166 S.L. Liu, Y. Pu, X.M. Cheng et al., Conceptual design of a water cooled breeder blanket for CFETR. Fus. Eng. Des. 89, 1380–1385 (2014). doi:10.1016/j.fusengdes.2014.01.065 J. Li, X.K. Zhang, F.F. Gao et al., Neutronics comparison analysis of the water cooled ceramics breeding blanket for CFETR. Plasma Sci. Technol 18, 179–183 (2016). doi:10.1088/10090630/18/2/14 Q.J. Zhu, J. Li, S.L. Liu, Neutronics analysis of water-cooled ceramic breeder blanket for CFETR. Plasma Sci. Technol 18, 775–780 (2016). doi:10.1088/1009-0630/18/7/13

123