Monte Carlo fission matrix acceleration method with adaptive mesh

NUCLEAR ENERGY SCIENCE AND ENGINEERING

Monte Carlo fission matrix acceleration method with adaptive mesh

PAN Liu-Jun，

WANG Rui-Hong，

JIANG Song

Nuclear Science and Techniques

Vol.26, No.5

Article number 050603

Published in print 20 Oct 2015

Available online 20 Oct 2015

DOI：10.13538/j.1001-8042/nst.26.050603

46601

In large-scale, loosely coupled systems, Monte Carlo criticality calculation always suffers from slow fission source convergence resulting from the high dominance ratio. The fission matrix acceleration method has shown its potential to accelerate the convergence of the fission source in many numerical simulations. In practice, however, instability of this method may be caused by imbalanced precisions of elements of the fission matrix. Hence, an improved method, in which the space mesh used to compute the fission matrix is defined adaptively based on the fission bank in each cycle, is introduced. The proposed method ensures balanced precisions of elements of the fission matrix, so is more stable than the existing fission matrix method.

Monte CarloCriticalFission matrixAdaptive

I. INTRODUCTION

Monte Carlo (MC) method is widely used to compute the multiplication factor and fundamental mode eigenfunction of critical systems in reactor physics. The power iteration method is the most important technique for criticality calculation in general Monte Carlo codes. However, this method always suffers from slow fission source convergence in large, loosely coupled systems. To achieve a sufficiently converged fission source, a large number of inactive cycles must be performed before active cycles start. This slow convergence may result in unacceptable computational cost, especially for whole core analysis of commercial reactors. To overcome this issue, various acceleration methods have been suggested in MC criticality calculations. The super-history method [1] and the Wielandt’s method [2] decrease the dominance ratio by modifying the power iteration scheme, and then reduce the necessary number of inactive cycles to reach the convergence of the fission source. However, they increase the time expense of computation during each inactive cycle. The functional method [3] and the coarse-mesh finite difference accelerating Monte Carlo method [4] employ standard Monte Carlo techniques to estimate nonlinear functionals, which are used in low-order functional Monte Carlo (FMC) equations and coarse-mesh finite dfference (CMFD) equations to obtain the eigenvalues and discrete representations of the eigenfunction. The modified power iteration method [5] can be used for convergence acceleration in addition to the capability of calculating the second eigenpair. There are still no examples or illustrations to demonstrate that the modified power iteration method will work for continuous energy problems.

Compared with the methods mentioned above, the fission matrix acceleration (FM) method [6] is more intuitive. When the system is spatially discretized, the fission matrix can be tallied during the Monte Carlo simulation at each cycle. The FM method adjusts the weights of fission neutrons according to the fundamental mode eigenvector of the fission matrix at each cycle. This method has shown its potential to accelerate the fission source in many numerical simulations. However, some instability issues restrict the applications of this method. Actually, the FM method constructs a space mesh before the iteration process and uses this fixed mesh to tally the fission matrix in each cycle. The elements of the fission matrix may be computed with various precisions, as the elements correspond to zones with different source intensities are tallied with a different number of samples. The fission source intensities may become very small in some zones after several cycles. This results in large statistical errors of the corresponding columns of the fission matrix. The imbalanced precisions of elements of the fission matrix may cause the FM method to be instable. This instability issue of the FM method was recently addressed in the semi-fixed-source FM method [7]. This method samples all elements of the fission matrix using an equal number of particles to improve efficiency. But the semi-fixed-source FM method increases the computational cost, as it samples many more particles than the existing method, especially when the real fission source of the system is severely imbalanced. In this paper, the space mesh used for computing the fission matrix is defined adaptively based on the fission bank in each cycle. This method ensures balanced precisions of elements of the fission matrix without increasing the number of sampled particles.

It is worth noting that "accumulation" [8] can be adapted to address similar issues. This method improves the precision of the fission matrix, but it is not able to balance the precision of elements of the fission matrix.

II. FISSION MATRIX METHOD WITH ADAPTIVE MESH

A. Fission matrix method

The main purpose of the criticality calculation is to find the fundamental mode solution of the eigenvalue equation [7]

H s (\vec{r}, E) = k s (\vec{r}, E),

(1)

where k is the eigenvalue, $s (\vec{r}, E)$ is the corresponding fission source distribution function, and

H s (\vec{r}, E) = \int_{0}^{\infty} d E^{'} \int_{V} d \vec{r}' f (\vec{r}', E^{'} \to \vec{r}, E) s (\vec{r}', E^{'}),

(2)

in which $f (\vec{r}', E^{'} \to \vec{r}, E)$ is the fission kernel and represents the expected number of first generation fission neutrons produced per unit volume at $\vec{r}$ per unit energy at E from a fission neutron born at $\vec{r}'$ with energy E’. The fission neutrons are assumed to emit isotropically (this assumption is not necessary for Monte Carlo simulation, only for simplifying the formulation).

As we know, the standard power method is used to solve problem (1) in most Monte Carlo codes. A guessed fission source distribution or a distribution obtained by other calculations is used as an initial fission source distribution, and then the distribution is iteratively recalculated until it approaches to convergence. The cycles before obtaining a convergent fission source are called inactive cycles, while active cycles are the cycles after achieving a convergent fission source. Inactive cycles are used to converge to a fission source and active cycles are used for tallying the final results.

When the system is spatially discretized, the real fission source distribution satisfies the following equation:

H S = k S,

(3)

where S indicates a fission source distribution vector, its elements represent the fraction of the fission source in each zone. k is the multiplication factor, and the fission matrix H is defined as follows:

H_{i, j} = \frac{\int_{0}^{\infty} d E \int_{V_{i}} d \vec{r} \int_{0}^{\infty} d E^{'} \int_{V_{j}} d \vec{r}' f (\vec{r}', E^{'} \to \vec{r}, E) s (\vec{r}', E^{'})}{\int_{0}^{\infty} d E^{'} \int_{V_{j}} d \vec{r}' s (\vec{r}', E^{'})} .

(4)

The "cycle fission matrix" [9] can be constructed in each cycle as

H_{i, j}^{(n)} = \frac{\int_{0}^{\infty} d E \int_{V_{i}} d \vec{r} \int_{0}^{\infty} d E^{'} \int_{V_{j}} d \vec{r}' f (\vec{r}', E^{'} \to \vec{r}, E) s_{m c}^{(n - 1)} (\vec{r}', E^{'})}{\int_{0}^{\infty} d E^{'} \int_{V_{j}} d \vec{r}' s_{m c}^{(n - 1)} (\vec{r}', E^{'})} .

(5)

As the fission source distribution evaluated from the fission matrix converges more rapidly than ordinary Monte Carlo simulations [7], the fission source convergence can be accelerated by adjusting the weights of fission neutrons according to the fundamental mode eigenvector of H⁽ⁿ⁾ in each cycle. If the m-th fission neutron is produced in region i, its weight, $w_{m}^{(n)}$ , is corrected to [10]

\frac{{(s_{f m}^{(n)})}_{i}}{{(s_{m c}^{(n)})}_{i}} ω_{m}^{(n)},

(6)

where $S_{f m}^{(n)}$ is the fundamental eigenvector of H⁽ⁿ⁾, and ${(S_{f m}^{(n)})}_{i}$ is the relative source intensity at zone i obtained from an ordinary MC simulation (normalized to the same number with $S_{f m}^{(n)}$ ). In practice, the fission matrix acceleration is only used in inactive cycles. It is recommended that a few inactive cycles without the weight correction should be set before moving into active cycles to eliminate the bias of region-averaged correction factors.

B. Adaptive mesh

In Monte Carlo simulation, the number of histories sampled in each zone is decided by the source intensity. When a fixed mesh is used for computing the fission matrix during all cycles in the FM method, the fission source intensities may become very small in some zones after several cycles. Therefore, the number of histories sampled in those zones are less. Imbalanced precisions of elements of the fission matrix may occur and result in instability of the FM method. In this paper, to improve efficiency and overcome instability, adaptive mesh is proposed as follows:

(1) A relatively fine mesh is defined, consisting of a relatively large number of zones.

(2) At the beginning of each cycle, construct a coarse mesh by combining the zones of the fine mesh according to the source distribution, making the source intensities of different zones of the coarse mesh roughly equal. Each zone of coarse mesh may include several zones of the fine mesh. A different number of zones in the coarse meshes for different cycles is allowed.

(3) During each cycle, the fission matrix is calculated based on the coarse mesh. Correspondingly, the order of the fission matrixes may be different in different cycles.

(4) During each cycle, count the number of the first generation fission neutrons in each zone of the fine mesh and construct the fission bank for the next cycle.

(5) After each cycle, adjust the weights of the neutrons in the fission bank according to the fundamental eigenvector of the fission matrix.

Step (2) is the key step of this method. This step is easy to implement for a 1 dimensional model. For this purpose, the fine mesh contains N zones can be defined by

a = x_{0} < x_{1} < \dots < x_{N - 1} < x_{N} = b .

(7)

The source intensity (according to the fission source bank given by the previous cycle) of zone Ii=[x_i-1,xi] is Δsi. If we want to construct a coarse mesh containing about K zones (K should be much smaller than N), the source intensity of each zone of the coarse mesh should be about

t = \frac{\sum_{i = 1}^{N} Δ s_{i}}{K} .

(8)

We first select k₁ satisfying the following equation:

| \sum_{i = 1}^{k_{1}} Δ S_{i} - t | = \min_{k} {| \sum_{i = 1}^{k} Δ S_{i} - t |},

(9)

let $\cup_{i = 1}^{k_{1}} I_{i}$ be the first zone of the coarse mesh, then find k₂ satisfying

| \sum_{i = k_{1}}^{k_{2}} Δ S_{i} - t | = \min_{k} {| \sum_{i = k_{1}}^{k} Δ S_{i} - t |},

(10)

let $\cup_{i = k_{1} + 1}^{k_{2}} I_{i}$ be the second zone of the coarse mesh, and so on. After this, if the intensity of the last zone of the coarse mesh is less than a settled value (for example, t/2), the last two zones of the coarse mesh are combined. This method is not able to make the intensity of each zone of coarse mesh equal to t exactly, so the number of zones in the coarse meshes may be a little different than K.

III. NUMERICAL RESULTS

We consider a one dimensional problem of large heterogeneous fissile regions surrounded by a thin reflector. The neutrons are assumed to be scattered isotropically, and the physical data is given in Table 1.

Data for test problem

Region	Location (cm)	∑_t (cm^-1)	∑_s (cm^-1)	v∑_f (cm^-1)
1	0 < x < 5	1.0	0.86	0.00
2	5 < x < 10	1.0	0.86	0.15
3	10 < x < 30	1.0	0.86	0.10
4	30 < x < 35	1.0	0.86	0.149
5	35 < x < 40	1.0	0.86	0.00

Monte Carlo simulations were performed with 500000 histories at each cycle. The initial fission source distribution is spatially flat. The FM method is performed with the fissile regions discretized to 10 zones. The 30 cm fissile region is equally discretized to 300 zones for the fine mesh. The amounts of zones of coarse mesh constructed for the fission matrix acceleration method with the adaptive mesh (ADFM) method are also about 10, but may be a little different at some cycles.

Figure 1 shows the number of fission neutrons in each zone corresponding to the initial source, cycle 5 and cycle 20, respectively, when the fixed mesh is used. From this figure, it can be seen that the number of fission neutrons becomes very imbalanced after several cycles. This would result in very imbalanced precisions of elements of the fission matrix, because only the histories started from zone j contribute to the tallies of the elements Hij of the fission matrix.

Fig. 1.

The numbers of fission neutrons of each zone corresponding to different cycles when fixed mesh is used.

Figure 2 shows the coarse meshes of the ADFM method corresponding to different cycles. The mesh of the FM method is the same as the mesh of the ADFM method at the first cycle. Notice that the meshes are obviously different at different cycles.

Fig. 2.

The coarse meshes for ADFM method corresponding to different cycles.

100 cycles are done with the FM method and the ADFM separately, the Shannon entropy [11] evolution is shown in Fig. 3 (the fine mesh is used for calculating the Shannon entropy). The average of the Shannon entropy over cycles 3001 to 5000 of the standard Monte Carlo power iteration (PI) is taken as the reference.

Fig. 3.

(Color online) Shannon entropies of first 100 cycles of different methods for one dimensional slab problem.

In Fig. 3, it can be seen that both the Shannon entropies of the FM method and the ADFM method cross the reference more rapidly than the standard PI and the Shannon entropies of the FM oscillate much more wildly than the ADFM method.

In some situations, it may need many cycles to eliminate this oscillation after terminating the weight correction when the FM method is used. For example, the Shannon evolution when the weight correction is terminated at cycle 48 is shown in Fig. 4.

Fig. 4.

(Color online) Shannon entropy evolution when the weight correction is terminated at cycle 48.

Table 2 shows the cycle numbers when the Shannon entropy crosses the reference for the first time after terminating the weight correction by FM and ADFM at different cycles (by standard PI, the Shannon entropy crosses the reference for the first time at cycle 851). It can be seen from this table that the ADFM method is more efficient and more stable than the existing FM method.

Cycle numbers when the Shannon entropy crosses the reference at first time after terminating the weight correction

The cycle terminating the weight correction		41	42	43	44	45	46	47	48	49	50
The first cycle when the Shannon entropy crosses the reference after terminating the weight correction	FM	690	473	51	283	69	257	237	786	515	249
	ADFM	81	188	170	64	114	93	57	62	130	227

The eigenfunction estimated during cycles 50 to 100 by the FM method, the ADFM method, and standard power iteration are shown in Fig. 5, in which the eigenfunction estimated during cycles 3001 to 3050 by the standard PI is taken as a reference.

Fig. 5.

(Color online) Eigenfunction estimated during cycles 50 to 100 of one dimensional slab problem.

IV. DISCUSSION

The convergence of source distribution can be diagnosed by checking whether the Shannon entropy crosses the reference [11]. It can be seen that both the Shannon entropies of the FM method and the ADFM method cross the reference rapidly from the numerical results (Fig. 3 and Fig. 4). It seems that the convergence rate of the ADFM method is the same as the FM method. But it can also be seen that the FM method can be too oscillatory to use in this model. To eliminate this oscillation, extra cycles without weight correction should be set before moving into active cycles. Actually, when the FM method is used, terminating the weight correction at different cycles may lead to different effects. In some cases, it requires many cycles to make the Shannon entropy cross the reference again after terminating the weight correction. The amount of extra cycles may even be comparable to the amount of cycles of standard PI to make the Shannon entropy cross the reference. But when ADFM is used, it always needs much fewer cycles to make the Shannon entropy cross the reference again after terminating the weight correction than standard PI (Table 2). Therefore, the ADFM method is more efficient and stable than the FM method.

Compared to the FM method, the extra time required for the ADFM method can be neglected in the tested 1D model. The same situation can be expected for 2D or 3D problems. Actually, the main extra cost of the ADFM method is to calculate the intensity of the fission source in each zone of the fine mesh. When ADFM is used, regular fine mesh can be defined (the union of the zones of fine mesh must contain the fissile regions but not need to be the same as the fissile regions). Then it is easy to find where a fission neutron in the fission source bank is located, so it is easy to count the number of fission neutrons (which is the intensity of the fission source) in each zone of the fine mesh.

The FM method has been implemented into the MCNP code to simulate 2D and 3D continuous-energy problems [12]. The only difference between our ADFM method and the FM method lies in how to construct the space mesh. Although the ADFM method is only tested on a multi-group problem for simplicity, there is no difficulty in applying the ADFM method to continuous-energy problems.

There is also no difficulty in applying the ADFM method to 2D and 3D problems. But for complex 2D or 3D systems, the coarse meshes must be constructed elegantly to make the ADFM method efficient. The reason is as follows. The weight corrections of the FM method and the ADFM method can only accelerate the balance of the fission source between different zones; they are not able to accelerate the balance of the local fission source in each zone. In the ADFM method, as the source intensities in each zone of the coarse mesh are roughly equal, local balance of the fission source in large zones are not important for their low source densities. For 1D problems, the balances of local fission source in small zones are easy to reach. But for 2D or 3D problems, the small area or volume of a zone does not certainly lead to a quick balance of the local fission source, because even if the area or volume of a zone is small, it is possible that the scale of one direction of that zone is large (for example, a long and narrow rectangle). To make the ADFM method efficient, the zones of coarse mesh must meet one extra requirement: the scale of one direction can not be much larger than others.

V. CONCLUSION

The fission matrix acceleration method has shown its potential to accelerate the convergence of the fission source in many numerical simulations. In practice, however, instability of this method may be caused by the imbalanced precision of elements of the fission matrix. In this paper, the space mesh used to compute the fission matrix is defined adaptively based on the fission bank in each cycle. This method ensures a balanced precision of elements of the fission matrix, so is more stable. Although the ADFM method is only tested on one dimensional multi-group problem for simplicity, there is no difficulty in using the ADFM method in two dimensional or three dimensional continuous-energy problems, on the condition of properly constructed adaptive mesh.

References

[1]

R N Blomquist and E M Gelbard.

Alternative implementations of the Monte Carlo power method

. Nucl Sci Eng, 2002, 141: 85-100. DOI: 10.13182/NSE01-30

Baidu Scholar

Google Scholar

[2]

T Yamamoto and Y Miyoshi.

Reliable method for fission source convergence of Monte Carlo criticality calculation with Wielandt’s method

. J Nucl Sci Technol, 2004, 41: 99-107. DOI: 10.1080/18811248.2004.9715465

Baidu Scholar

Google Scholar

[3]

E W Larsen and J Yang.

A functional Monte Carlo method for k-eigenvalue problems

. Nucl Sci Eng, 2008, 159: 107-126. DOI: 10.13182/NSE07-92

Baidu Scholar

Google Scholar

[4]

E R Wolters, E W Larsen and W R Martin.

Hybrid Monte Carlo-CMFD method for accelerating fission source convergence

. Nucl Sci Eng, 2013, 174: 286-299. DOI: 10.13182/NSE12-72

Baidu Scholar

Google Scholar

[5]

T Booth.

Power iteration method for the several largest eigenvalues and eigenfunctions

. Nucl Sci Eng, 2006, 154: 48-62. DOI: 10.13182/NSE05-05

Baidu Scholar

Google Scholar

[6]

T Kitada and T Takeda.

Effective convergence of fission source distribution in Monte Carlo simulation

. J Nucl Sci Technol, 2001, 38: 324-329. DOI: 10.1080/18811248.2001.9715036

Baidu Scholar

Google Scholar

[7]

J Dufek and W Gudowski.

Fission matrix based Monte Carlo criticality calculations

. Ann Nucl Energy, 2009, 36: 1270-1275. DOI: 10.1016/j.anucene.2009.05.003

Baidu Scholar

Google Scholar

[8]

S Yun and N Z Cho.

Acceleration of source convergence in Monte Carlo k-eigenvalue problems via anchoring with a p-CMFD deterministic method

. Ann Nucl Energy, 2010, 37: 1649-1658. DOI: 10.1016/j.anucene.2010.07.018

Baidu Scholar

Google Scholar

[9]

L L Carter and N J McCormick.

Source convergence in Monte Carlo calculations

. Nucl Sci Eng, 1969, 36: 438-441. DOI: 10.13182/NSE69-1

Baidu Scholar

Google Scholar

[10]

H Kadotani, Y Hariyama, M Shiota, et al. Acceleration of fission distribution convergence using eigenvectors from matrix k calculations in the KENO code. ICNC 91, Oxford, UK, paper II-1, Sep. 9-13, 1991.

[11]

T Ueki and F B Brown.

Stationarity modeling and informatics-based diagnostics in Monte Carlo criticality calculations

. Nucl Sci Eng, 2005, 149: 38-50. DOI: 10.13182/NSE04-15

Baidu Scholar

Google Scholar

[12]

F Brown, S Carney, B Kiedrowski, et al.

Fission matrix capability for MCNP Monte Carlo

. Joint International Conference on Supercomputing in Nuclear Applications and Monte Carlo. LA-UR-13-26962, La Cité des Sciences et de l’Industrie, Paris, France, Oct. 27-31, 2013. 10.1051/snamc/201403502

Baidu Scholar

Google Scholar