2.1. BP Neural Network

An efficient and memory-saving method for evaluating cross-sections is needed for on-the-fly Doppler broadening. Neural networks, which are widely used in machine learning, need to be trained before they are used. Once the parameters of the networks are determined, the networks can be used to calculate the output from the given input. After training, the amount of calculations required is greatly reduced; therefore, this method can be used for on-the-fly Doppler broadening. A combination of the input and expected results is used in the training process, during which the weights and biases of all the neurons are adjusted. The deviation between the outputs of the network and the expected results is gradually reduced during training until the network meets the requirements.

Many studies on neural networks such as convolutional neural networks ^[16] and deep neural networks ^[17] have been conducted. Such networks have many hidden layers and neurons, as well as complex structures, resulting in low computational efficiencies. Because computing speed is important in Monte Carlo codes, a neural network with a simple structure is more suitable.

In this study, the back propagation (BP) network was used. The error back propagation algorithm proposed by Rumelhart ^[18] introduced in the following paragraphs is used for training this network. As shown in Fig. 1, a BP network consists of an input layer, hidden layers, and an output layer. Each layer contains a certain number of neurons. The simplicity of the structure and calculation steps of this type of network ensures its high computational efficiency.

Fig. 1Full-screen View

Structure of BP neural network.

The hidden and output layer nodes of the BP neural network are respectively described by

In the equations, $O_j$ is the output of the hidden layer node; $Y_k$ is the output value of the output layer node. $f_1$ and $f_2$ are the transfer functions of the hidden and output layer nodes, respectively. In general, $f_1$ is a nonlinear function, and $f_2$ can be a linear or nonlinear function. $w_{ij}$ is the weight from the previous node to the hidden layer node. $w_{jk}$ is the weight from the node of the last hidden layer to the node of the output layer. $b_j$ and $b_k$ are the biases of the hidden layer node and output layer node, respectively.

The weights and thresholds of the BP neural network in calculations are generally expressed in matrix form. Equations (1) and (2) are therefore expressed as Eqs. (3) and (4), respectively.

The subscript “pre” means the previous layer, which can be a hidden layer or an input layer. Subscript “h” and “o” means hidden layer and output layer, respectively. In Eq. (3), $O_{\text{h}}$ is the output value of the hidden layer. $W_{\text{pre to h}}$ is the matrix of weights from the previous layer to the hidden layer. $X_{\text{pre}}$ is the matrix of output values of the previous layer and $B_{\text{h}}$ is the matrix of bias of the hidden layer. In Eq. (4), $Y$ is the matrix of output value of the output layer, and it is also the output value of the network. $W$ and $B$ in Eq. (4) have similar meanings to those in Eq. (3). The bias and output data of each layer are column vectors with lengths equal to the number of nodes in the layer. The weight is a matrix where the first and second dimensions are the number of nodes in the previous and current layers, respectively. The matrix elements correspond to those in Eqs. (1) and (2).

The neural networks in this study were trained using MATLAB. The output of the neural networks approached the target value over successive iterations. All the weights and biases in the network were adjusted during the training process, while the structure of the network, including the number of hidden layers and nodes in each layer, and the transfer function, remained unchanged.

The error back propagation algorithm was used for training the network. The error between the result in the training data and the corresponding result calculated by the network for a given input in the training data was propagated back to the parameters of each layer. The process is briefly described as follows.

The set of input vectors and the corresponding target vectors are denoted as

The square of errors between P and T are summed to defined R, which is the error of the network. The factor 1/2 is to simplify the following process.

The goal of the training is to reduce the error. R is therefore expanded using Eq. (4) into the node parameters of the output layer as follows:

Equation 8 can be further expanded using Eq. (1) as follows:

The expansion ends here if the network has one hidden layer. Otherwise, R can be expanded layer by layer. The algorithm is illustrated for a one-hidden-layer network for which the transfer function of the output layer is given by

The partial derivatives of R in Eq. (8) are

The partial derivatives are called the gradient values of the error. The parameters are updated using the gradient values in each training Iteration as follows:

where $\eta$ is the learning rate, which determines the ratio of the correction to the parameters. The amplitude of a single adjustment is directly proportional to $\eta$. The derivation shown in Eqs. (1) to (14) is applied to each individual parameter, and the end result of the derivation can be expressed in the form of a matrix as

The results of Eq. (13) and (14) for all the parameters can be collectively expressed as

The training of the other parameter matrices can be described in a similar manner. The error plays a role in the parameter adjustment of each layer through the result deviation $\Delta$ and propagates to each parameter matrix.

2.2. Training and Calculation

This section describes how the two networks were trained and used in an energy window. The method for dividing the resolved energy is introduced in the next section. Both the energy division and the sequential computation using two networks are designed to reduce the complexity of the training target.

A temperature is chosen as the base temperature for the cross-section of a nuclide. The cross-section at this temperature, denoted as ${\sigma }_{T_0}\left(E\right)$, is a function of the neutron energy $E$. The cross-section broadening coefficient, $K$, is introduced, which is defined as the ratio of the cross-section at a temperature not less than $T_0$ to the cross-section at $T_0$.

The HWN method is similar to the method proposed by Yesilyurt et al. ^[6], which also uses fitting parameters for on-the-fly Doppler broadening. However, the broadening coefficient $K$ in the proposed method applies for all the energies in the corresponding energy window, whereas the method proposed by Yesilyurt et al. requires individual parameters for each energy point.

To reduce the memory required, a network denoted as Network 1 is trained to calculate the cross-section at temperature $T_0$. The upper and lower bounds of the energy window are denoted as $E_{\text{min}}$ and $E_{\text{max}}$, respectively. The corresponding relationship between the input and output of the network is expressed as

where ${\sigma }_{T_0}^{\text{"}}$ is the cross-section calculated by Network 1 at $T_0$ and $F_1$ is the mapping relationship between the input energy and output cross-section of Network 1. It should be noted that some errors exist between the results and the actual cross-section.

Network 2 is trained to calculate the broadening coefficient, $K$. If $K$ is used directly as the input data for training Network 2 and Eq. (18) is used to calculate the broadened cross-section, the error of the calculation result will be the superposition of the errors of the two networks.

To reduce the error, the $K$ value calculated by the following formula is used as the training target.

The difference between Eqs. (20) and (18) is that the cross-section at $T_0$ in Eq. (18) is replaced by the output of Network 1 in Eq. (20). Therefore, the influence of the accuracy of Network 1 on the final result can be eliminated. However, the complexity of $K$ will increase if the error of Network 1 is extremely large, which may affect the accuracy of Network 2. It is therefore necessary to control the error of Network 1.

$K$ in Eq. (20) is used as the input data to train Network 2. The temperature range of the network fitting is limited by the input data. A temperature range denoted as $T_{\text{min}}\le T\le T_{\text{max}}$ that could meet the actual requirements is selected. $T_{\text{min}}$ should not be lower than $T_0$. The corresponding relationship between the input and output can therefore be expressed as

where $K^{\prime }$ is the broadening coefficient calculated by Network 2, and $F_2$ is the mapping relationship between the input energy and the output coefficient of Network 2. $K^{\prime }$ is a function of the independent variables $E$ and $T$ and is denoted as $K^{\prime }\left(E,T\right)$.

By rewriting Eq. (20), the equation for on-the-fly Doppler broadening is obtained as

Equation (22) describes the Doppler broadening process. Network 1 is first used to calculate the cross-section at the given energy $E$ and basic temperature $T_0$. Then, the cross-section is broadened to the temperature $T$ using the broadening coefficient $K$ calculated by Network 2.

The values of $T_0$, $T_{\text{min}}$, and $T_{\text{max}}$ were determined. The temperature ranges and their corresponding fields of study were summarized by Yesilyurt et al. ^[6] as shown in Table 1. Monte Carlo codes for reactors should be able to handle benchmarking calculations and reactor operation problems. Thus, $T_{\text{min}}$ and $T_{\text{max}}$ were set as 250 K and 1650 K, respectively. It should be noted $K$ approaches 1 as $T$ approaches $T_0$, which affects the accuracy of the neural network in the vicinity of $T_0$. As a result, training using data at $T=T_0$ should be avoided, and $T_0$ should not be close to $T_{\text{min}}$. Considering the lower complexity of the cross-section curve at higher temperatures, the base temperature $T_0$ was chosen as 200 K to reduce the difficulty of the network training.

2.3. Window Division and Parameter Determination

The cross-section curve in the resolved resonance energy range of heavy nuclides at a single temperature is very complex because of the large number of resonance peaks. The addition of the temperature dimension further increases the complexity. Therefore, it is impractical to train the neural network directly using the ACE data over the entire resolved resonance energy range as the input.

Dividing the resolved resonance energy range into energy windows can significantly reduce the difficulty of training in each window. Because similar processes are carried out for each window, the data should be divided evenly according to the training difficulty. Dividing the resonant peaks equally between the different windows is an easy and effective method. Each window has the same number of maximum points that correspond mainly to the positions of the resonance peaks within the resolved resonance energy range. The edges of the windows are set as the minimum points of the cross-section curve.

Because of the Doppler broadening effect, the resonant peaks are broadened at temperatures above 0 K. At $T_0=200\text{ K}$, some adjacent resonant peaks are combined into single peaks, leading to a reduction in the number of maximum points. The maximum points, rather than the resonant peaks at 0 K, were used for the deviation of the cross-section curve at $T_0$. Because the distributions of the resonant peaks differ for different heavy nuclides, the number of windows and the positions of the edges determined by this process will also be different.

The number of maximum points in each window was carefully determined. A smaller number of points for a given network would result in a higher number of divided windows and a corresponding increase in the memory required. In addition, if large number of peaks are included, the amount of data in each window will be extremely large, and the accuracy of the networks will be decreased. The number of points was set to 15 based on experience and testing. The aforementioned process was applied to the total cross-section curve at the base temperature of $T_0=200\text{ K}$ for window division.

Relatively large deviations near the boundary of the windows were often observed during network training. To avoid this, the windows were extended along both edges, as is shown in Fig. 2. The data in the extended window were used for neural network training, whereas the acquired neural network was used only within the range of the original window. The windows next to the unresolved resonance energy range and thermal energy range were only extended in the direction of the resolved resonance energy range.

Fig. 2Full-screen View

Extension of window.

The training results were greatly affected by the structure of the neural networks and the training parameters. The network parameters were determined according to the test results.

The training function training was used in network training. This approach has the advantages of a faster convergence speed and better accuracies. The transfer function of the hidden layer nodes is the tansig function in MATLAB, which is defined as

A BP network with one hidden layer is sufficient for solving many problems. Networks with multiple hidden layers may show better performance in some cases ^[19]. A one-hidden-layer BP network was used for Network 1 (described in Sect. 2.2) to reduce the memory requirements. For Network 2, a two-hidden-layer network was used to avoid overfitting in the training of K(E, T).

Tests were performed to determine the number of neurons in each hidden layer. The number of neurons in the hidden layer of Network 1 was determined first. The total cross-section data of U-235 at 0 K were divided into 181 windows, and 19 equally spaced windows were selected. Networks with 80, 90, 100, 110, 120, and 130 nodes in the hidden layer were trained using the data in the windows. Each network was trained with the data from each window 10 times, and the minimum value of the maximum absolute relative error was calculated. The geometric means of the error values for the 19 windows are shown in Fig. 3.

Fig. 3Full-screen View

Relationship between the number of nodes in the hidden layer and the mean minimum absolute relative error.

In general, the accuracy of the network increased with the number of nodes, and it was at an acceptable level when the number of nodes was 130. As shown in Fig. 3, a further increase in the number of nodes had a limited effect on the improvement of accuracy. Networks with 130 nodes in the hidden layers were used for most windows. A larger number of nodes was used in a few windows in which the relative error was extremely large.

Training and comparisons were performed to determine the number of nodes in the two hidden layers of Network 2. A reasonable maximum epoch number, training time limit, and accuracy target were set for training. Data from the second window of the U-235 total cross-section were used to evaluate networks with different combinations of node numbers. The performance was determined using the percentage of data points where the absolute relative error with respect to the input data was less than 0.1%. The results are compared in Table 2. The results indicate that the best combination has 40 nodes in the first hidden layer and 20 nodes in the second hidden layer.

3.1. Microscopic Accuracy and Memory Requirement Comparison

The HWN method was applied to U-235 and U-238, which are representative heavy nuclides with important resonances. The total cross-section, elastic scattering cross-section, absorption cross-section of U-235, and total cross-section of U-238 were compared with the ACE data at the same temperature. The cross-sectional results and absolute relative errors with respect to the ACE data are plotted in Fig. 4, where the comparisons were made in the energy regions with strong resonances at 300 K and 700 K. It can be observed that the errors are low for most data points and that the cross-sections evaluated are accurate at both temperatures. The relative errors fluctuate with the energy. The fluctuations are caused by the characteristics of neural networks with many nodes.

Fig. 4Full-screen View

(Color online) Comparison of microscopic cross-sections at 300 K and 700 K.

The HWN method and the windowed multipole method are both methods with new ways, which do not use continuous-energy data, of storing the cross-section data. Therefore, their memory requirements are lower than that of the ACE data. The relative errors of the HWN and windowed multiple ^[1] methods are compared in Table 3. The results show that the maximum relative errors are similar and that the average relative errors of both methods are within 0.1%.

The theoretical memory consumption of the network parameters and ACE cross-section data are compared. In the continuous-energy Monte Carlo code using the ACE data, the cross-sections are stored in the form of double-precision floating-point numbers. Each double-precision floating-point number occupies eight bytes of memory. Both the energy grid and cross-section values are needed for cross-section evaluation. For the total cross-section, elastic scattering cross-section, and absorption cross-section, four double-precision floating-point numbers, which occupy 32 bytes of memory, are needed for each energy point. In comparison, most of the parameters in the HWN method are stored in the form of double-precision floating-point numbers, and a few parameters are integers. The memory consumptions of the two methods are calculated using the number of parameters and the memory space needed for each parameter.

The ACE data at 0 K processed by NJOY were used for comparison. If on-the-fly Doppler broadening is not introduced to the Monte Carlo code, point-wise cross-sections at more than a dozen temperatures will be needed for thermal-hydraulics coupled analysis. The method proposed by Yesilyurt et al. uses parameters that require several times the memory needed for the ACE data at 0 K. Point-wise data are also needed in the TMS method. As a result, the HWN method shows a significant reduction of the memory requirement compared to the 0 K ACE data.

The comparison results are listed in Table 4. The results show that for the three cross-sections of U-235, the HWN method could reduce the memory requirements of cross-section data in the resolved resonance energy range by 66.1% as compared with the case of the 0 K ACE data. The memory requirement reduction was 65.9% over the entire energy range.

The networks for each window can be used independently because the training process of each window is independent. It is easy to set the scope of the method when the method is implemented in a Monte Carlo code. It is not necessary to store the 0 K ACE data within the selected windows if the HWN method is used. However, the speed of the cross-section evaluation in these windows will decrease. The efficiency drop is described in Sect. 4. Table 4 clearly shows that the resolved resonance energy range accounts for most the nuclear data. The memory optimization ratio is significant when the HWN method is used in some of the windows. Users can decide the efficiency to be compromised for saving memory.

3.2. Comparison of Macroscopic Test Results

The HWN method was used for Doppler broadening in all windows of the resolved resonance energy in the HWN cases, while ACE data from the processed ENDF/B-VII.0 database at the same temperature were used in the ACE case. All calculations were performed using an Intel i7-9750H CPU without parallelism.

3.2.1. Concentric spheres

The first comparison example consisted of two concentric spheres, as shown in Fig. 5. The radii of the inner and outer spheres were 6.782 cm and 11.862 cm, respectively. The inner sphere was filled with Material 1. The space between the two spheres was filled with Material 2. The outside of the outer sphere was set to vacuum. The nuclide compositions of the materials are listed in Table 5. The temperature of each material was set to 300 K.

Fig. 5Full-screen View

(Color online) Concentric spheres example.

Table 5Full-screen View

Nuclide composition of materials in concentric spheres example

Materials	Nuclides	Atomic density (${10}^{24}\text{atoms}\cdot {\text{cm}}^{-3}$)
Material 1	U-235	4.4917×10-2
	U-238	2.5993×10-3
	U-234	4.9210×10-4
Material 2	U-235	3.4428×10-4
	U-238	4.7470×10-2
	U-234	2.6299×10-6

Show more

The calculation parameters are listed in Table 6, and the results are presented in Table 7. The deviation in $k_{\text{eff}}$ is very small. There was a slight increase in the calculation time. The accuracy of the HWN method was therefore confirmed.

Table 6Full-screen View

Calculation parameters for the example of concentric spheres

Parameters	Value
Neutrons per cycle	100000
Inactive cycles	400
Active cycles	1600

Show more

Table 7Full-screen View

Calculation results of the concentric spheres example

Case	$k_{\text{eff}}$	Standard deviation	Calculation time (min)	Time ratio
ACE	0.995747	0.000048	597.3738	1.00
HWN	0.995748	0.000048	601.7461	1.01

Show more

3.2.2. PWR assembly

The second comparison example is the PWR assembly model shown in Fig. 6. The model comprised an infinite cylinder with a cross-section of 21.42 cm × 21.42 cm. There were 264 fuel rods and 25 pipes arranged in a 17 × 17 square. The fuel rods were cylinders with diameters of 0.8192 cm. Each fuel rod was surrounded by a 0.082-mm air layer and a 0.572-mm zirconium wall. The inner and outer diameters of the pipes were 1.138 and 1.2294 cm, respectively. The pipes were filled with water, and the material of their walls was zirconium. The remainder of the assembly was filled with water. The nuclide compositions of the fuels are listed in Table 8. The temperature of all the materials was set to 700 K.

Fig. 6Full-screen View

(Color online) Schematic diagram of the PWR assembly.

Table 8Full-screen View

Nuclide composition of fuel in the PWR assembly

Nuclide	Atomic density (1024 atoms/cm-3)
U-235	6.9100 × 10-3
U-238	2.2062×10-1
O-16	4.5510×10-1

Show more

The calculation parameters are listed in Table 9, and the results are presented in Table 10. The difference in $k_{\text{eff}}$ between the two cases is within twice the standard deviation, indicating that there is no significant difference. The calculation time of the HWN case is 1.39 times that of the ACE case, which indicates that using the HWN method will prolong the calculation time. The calculation time is shorter if the method is not used in all the windows.

Table 9Full-screen View

Calculation parameters of PWR assembly.

Parameters	Value
Neutrons per cycle	40000
Inactive cycles	400
Active cycles	1600

Show more

Table 10Full-screen View

Calculation results of PWR assembly.

Case	keff	Standard deviation	Calculation time (min)	Time ratio
ACE	1.349355	0.000068	2174.0543	1.00
HWN	1.349255	0.000062	2961.4330	1.39

Show more

The neutron flux spectrum of the fuel in a fuel rod adjacent to the central pipe and the neutron flux inside each fuel rod and pipe were calculated for both the ACE and HWN cases. The results are presented in Fig. 7. The blocks in Fig. 7c do not represent the actual geometry, rather they the corresponding positions.

Fig. 7Full-screen View

(Color online) Comparison results of ACE and HWN cases.

Figure 7a shows that the neutron flux spectra of the fuel rods are the same. Figure 7b shows that the relative deviations of the fluxes in most statistical intervals are within three times the standard deviation of the ACE case. Figure 7c compares the neutron fluxes of all the fuel rods in the assembly. There is no significant difference between the ACE and HWN cases. The blank grid squares represent the pipes whose fluxes are not shown in this figure. A numerical comparison shows that the deviation of most fuel rods and pipes is within twice the standard deviation, and the deviation of the remaining fuel rods and pipes is within three times the standard deviation except for a few fuel rods and pipes. The comparison results therefore demonstrate the accuracy of the proposed HWN method.