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

NUCLEAR ENERGY SCIENCE AND ENGINEERING

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

Xing-Chen Nie，

Jia Li，

Song-Lin Liu，

Xiao-Kang Zhang，

Ping-Hui Zhao，

Min-You Ye，

German Vogel，

Xiao Yang，

Qing-Jun Zhu

Nuclear Science and Techniques

Vol.28, No.8

Article number 115

Published in print 01 Aug 2017

Available online 07 Jul 2017

DOI：10.1007/s41365-017-0270-3

48300

It can be difficult to calculate someunder-sampled regions in global Monte Carlo radiation transport calculations. The global variance reduction (GVR) method is a useful solution to the problem of variancereduction everywhere in aphase space. In this research, a GVR procedure was developed and applied to the Chinese Fusion Engineering Testing Reactor (CFETR). A cylindricalCFETR modelwas utilized for comparingvarious implementations of the GVR methodto find the optimum.It was found that the flux-based GVR methodcouldensure more reliable statistical results, achievingan efficiency being7.43 times that of theanalog case. A mesh tally of the scalar neutron ﬂux was chosen for the GVR method to simulate global neutron transport in the CFETR model. Particles distributed uniformly in the system were sampled adequatelythrough teniterations ofGVR weight window. All voxels were scored, and the average relative error was2.4% in the ultimate step of the GVR iteration.

Global variance reductionWeightwindowMonteCarloMCNPNeutronics

1 Introduction

There is a growing need for performinghigh-resolution radiation transport and shutdown dose analyses for fusion reactors [1,2]. The results of global Monte Carlo (M-C)simulationsarestatistically 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, whichaims to increase computational efﬁciency and reduce variance for alocalized tally. Nevertheless, the use of a weight window generator in a large and composite target cell implies an obvious risk. Theparticles thatreach 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) methodis an efficient solution to global Monte Carlo particle transport problems. Based on weight windowsthat uniformly transport non-analog particles throughoutthe model, which allows all areas of the system to be adequately sampled, GVR improves computational efﬁciency of hypothetical global detectors, and reducesthe global variance. M. A. Cooper and E. W. 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. A. Davis and A. Turner [11-12], and J. Naish et al. [13], appliedthis kind of GVR method to specific fission and fusion calculation problems. A.J. van Wijket al.[14] used a GVR procedure to a fission case.In this paper,we estimate the forward flux to establish the weight windows, and adoptiterationmethod in the M-C simulation.A GVR procedure isdeveloped for the global Monte Carlo particle transport problems of CFETR(Chinese Fusion Engineering Test Reactor).

2 GVR method

For large and complex fusion devices such as CFETR, the dense componentsof shielding material may wellattenuate the neutronpopulation density, causingrelatively 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)toavoidlarge 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 ﬂux $ϕ_{}$ as a way to set the WW thresholdsas Eq. (1)

W W_{i} = \frac{C_{1}}{I M P_{i}} = C_{2} \frac{ϕ_{i}}{Max (ϕ_{})} .

(1)

The lower WW threshold ( $W W_{i}$ ) is set inversely proportional to the importance ( $I M P_{i}$ ) of a certain mesh element i. $C_{1}$ and $C_{2}$ are constants defined by different applications. Eq.(1) ensuresthat neutrons are more likely to be rouletted in mesh elements of low importance (high ﬂux),hence the increase ofcomputational efficiency. It also allows particles to be split in the mesh elements of high importance (low ﬂux),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 particlesmay vary greatly from one region to another, whichprevents the distant regions from being under-populated, and the close-to-sourceregions from being over-populated.

Figure 1 shows the flow chart 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 (therelative 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 enough neutrons spread through the whole model and the simulation obtains credible global results.

Fig. 1

Flow chart of GVR method based on neutron flux.

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 beingreﬂectory boundaries,was employed to assess the efﬁciency 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 was performed to obtain neutron distributiondata,which wasutilized to generate different cell-based GVR WW.Most kinds of distribution data were inversely related to the cell importance. Given this relation, differenttypes of data were selected from a wider range of possible candidates.For example, A. Turner and A. Davis [11-12] compared the flux,weight and population.A.J. 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 theyroughly decrease when neutrons penetrate deeper. They are important among the reference data for Monte Carlo transport.

Fig. 2

Cylindrical CFETR model (in cm).

For estimating GVR calculation efﬁciency, we introduced a globalﬁgure of meritFOM_Gfrom Ref.[11] to represent the efficiency degree of the GVR WW. It is defined as FOM_G =N/(T_cpu· $\sum ε_{i}^{2}$ ), (i =1→N), where ε_i is the error in the i^th cell, T_cpu is the CPU time in minutesandNis the number of total mesh elements.The spatial distribution of relative errors throughout the mesh should be as ﬂat as possible. The standard relative error σ_err= [( $\sum ε_{i}^{2}$ )/N− (∑ε_i)²/N²]^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 ﬂux of the highestFOM_G (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 σ_err (0.16). The flux-based GVR method achieved an efficiency of 7.43 times higher than that of the analog method. Fig.3 shows the analog results yielded with an obvious uncertainty in the area far away from the first wall, especially in the outboard blanket. Incomparison, the flux-based GVR method providesmore valuable flux data distributed across the whole model.

Comparison of different distribution data.

Distribution data	FOM_G	Average error	σ_err	Scoring (%)
Flux	54.32	0.108	0.16	100
Weight	51.04	0.109	0.17	91.03
Energy	44.69	0.154	0.14	100
Population	41.16	0.110	0.19	100
Track	39.51	0.108	0.2	97.44
Error	19.02	0.162	0.28	93.59
Analog	7.31	0.302	0.43	74.36

Fig. 3

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

3.2 Applicationsto a three-dimensional CFETR model

Bycomparingvarious forms of GVR methods, a mesh tally of the scalar neutron ﬂux was chosen for the GVR methodto 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) were0.7 m 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 of50 divisions along the Xand Ydirections and 80 divisions along the Zdirectionwas utilized to cover the entire reactor model. The voxel dimension was about10cm× 25cm× 25 cm.

Fig. 4

(Color online) Neutronics model of CFETR.

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 lesseffective 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 WWwitha 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 ﬂux producedbythe analogmethod 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 particleswere tried, the dense componentof the reactor material greatly attenuated the neutrons.In the subsequent iterations with the GVR WW set by the scalar ﬂux, thedistribution of particles in the system becamemore extensive and uniform.The flux spreadacross a larger region of the model, resulting in a larger range of spatial WW,which facilitatedthe GVR function in the next iteration.Wenote that neutrons have little possibilityof being transported into the three ports filled with the shielding materials of water andsteel. 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 efﬁciency of the integrated iteration calculation.

Fig. 5

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

In Table 2, the σ_err decreases with the iteration times, and reaches 0.09 in the final iteration, which means that the GVR method is capable of ﬂattening the relative error distribution of a global mesh tally. As expected, the average error decreases with increasing iteration times, while the scoring voxelsincreases.The average relative error changes from 78.8% under the analog method (Step 0) to 5.7% in the 8^th step. The percentage of scoring voxels increases from 27.79% to 100%.

Summary of GVR weight window iteration (time in min)*

Iterations	CTME	CTM	Computer time	Av. time perhistory	Parallel efficiency (%)	σ_err	Average error	Scoring(%)
0	1000	1002	1480	3.24×10⁻⁵	67.66	0.38	0.788	27.79
1	1000	1030	1232	3.32×10⁻⁴	83.57	0.42	0.732	42.77
2	1000	1040	1259	4.95×10⁻⁴	82.54	0.42	0.657	56.87
3	1000	1026	3633	1.14×10⁻³	28.24	0.34	0.390	88.56
4	1000	2908	19159	2.90×10⁻²	15.18	0.24	0.324	95.57
5	5000	5069	27503	1.27×10⁻²	18.43	0.22	0.215	97.62
6	10000	12178	74467	2.03×10⁻²	16.35	0.19	0.189	99.02
7	20000	20666	138240	2.95×10⁻²	14.95	0.17	0.132	99.77
8	40000	40517	207360	9.42×10⁻³	19.54	0.09	0.057	100.00

* CTME, computer time spent on the transport portion of the problem, which is set before the simulation.

In theory, the average time per history should increase as the WW covers more areas in the latter iterations. Nonetheless, Table 2 showsthat this value fluctuates through this series of iterations, whichindicates that some special histories have an abnormal lifespan.Also,the CTMs are greaterthan 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 timesthat of the CTME. Regardless of how we adjusted the CTME (whichwas below 2907.95 min), the CTM was constantly 2907.95 min. This suggests that certain neutron history may bevery 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 whenrunning MCNP in parallel mode. Inthis simulation, the MCNP was adapted to parallel processing with 24 processors. As a result, the computer time was much longer than thatfor 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, thoughit was not detailed enough. Fewerdevice gaps may make the fewer high statistical weight particles stream to the regions where a large amount ofWW splitting would take place.Bulk shielding and void gap are main reasons for the long history problem. When the weights of the particleshave 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 in the later iterations, the entire iteration processwas completedwithin an acceptable time range.

Figure 6 shows the fraction of mesh voxels under a certain error threshold as a function ofthe relative error. The more iterations are simulated, the larger the fraction of meshes that present an acceptableerror in the results. The lines move rapidly to the top left, where greater fractions of meshvoxels are plotted atsmall statistical errors.As expected, the analogsimulation (AN) performed the poorest, and only16.3% of the total voxelshad an error below 10%. Neutrons can scarcely be transported beyond the shielding material and filled ports; while 88.7% of the total voxels had an error under 10% in the 8^th step of adequately sampled iterations. So, the GVR methodcan improve the variance reduction for the global model.

Fig. 6

(Color online) Cumulative distribution of relative errors for Iteration 0–8.

In Fig. 7, the GVRand analogresultsare compared at the final step of iteration (the 10^th iteration)in the same computer time (10000 min).The flux and distribution of analog MonteCarlo particles vary by many orders of magnitude far from the source,withlarge statistical errors in regions far from the source.Using theGVR method,MonteCarlo particles are uniformly distributed throughout thesystem. All meshes are tallied (scoring=100%). The relative error in the scalar fluxremains considerably flat (σ_err=0.05). In an average relative error of 0.024, few errors in each mesh are greater than 0.05.

Fig. 7

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

4 Conclusion

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

References

A. Davis,

Radiation shielding of fusion systems. Ph.D. Thesis

, University of Birmingham, 2010.