One-dimensional reactor diffusion equation for a single energy group

In this section, a single-group k-eigenvalue problem is introduced for the criticality calculations. The largest value of k (known as the effective neutron multiplication factor or k_eff) should be determined. The single-group neutron diffusion model is given by Eq. (1):(1)where denotes the neutron flux of the energy group E in the r-coordinate at the instant t, v represents the neutron velocity, D denotes the diffusion coefficient, ν represents the neutrons per fission numbers, χ(E) represents the prompt neutron spectra, ∑_t represents total macroscopic cross-section, ∑_s represents the macroscopic scattering cross-section from group E’ to group E, ∑_f represents the fission macroscopic cross-section, and S(r, E, t) represents the neutron source.

In the absence of S(r, E, t), the initial neutron flux density is symmetric along the x-axis. Eq. (1) can be simplified into Eq. (2)[1]:(2)Here, k_∞ represents the infinite multiplication factor, and L denotes the diffusion length. Consider a uniform bare reactor [41] that is shaped as an infinite-plate bare reactor with dimensions of infinite length and width, and a thickness (including the extrapolation distance) of a. This is illustrated in Fig. 1. The analytical solution for the neutron flux can be obtained using the separation-of-variables method, as shown in Eq. (3):(3)

Fig. 1Full-screen View

Infinite Plate Reactor

The critical conditions for a bare reactor in the single-group approximation are given by Eq. (4).(4)where k_eff is the effective neutron multiplication factor. When the system is in a critical state, the neutron flux-density distribution satisfies the wave equation according to the fundamental eigen function corresponding to the minimum eigenvalue . This is given by Eq. (5):(5)If the composition of the system material (i.e., k_∞ and L²) is given, a unique critical size (denoted by a₀) results in k_eff=1. The critical size a₀ corresponds to the critical state of the reactor. For reactor sizes larger than a₀, k_eff>1, which indicates that the reactor is in a supercritical state. Conversely, for reactor sizes smaller than a₀, k_eff<1, which indicates that the reactor is in a subcritical state [42]. However, if the reactor size is given, it is feasible to determine a fuel enrichment (material composition) that satisfies Eq. (4) and ensures that the reactor attains criticality. When the system is in the critical state, the neutron flux density distribution within the reactor can be described as follows:(6)

Two-dimensional reactor diffusion equation for a single energy group

Based on Sect. 2.1, consider a two-dimensional neutron diffusion equation formulated as follows:(7)By discretizing Eq. (7) and approximating the Laplace operator and partial derivatives, the equation can be converted into the following form:(8)where denotes the value of the neutron flux at the spatial grid point (i, j) at time n. Δt is the time step. Δx and Δy are the spatial steps in the x- and y-directions, respectively. D is the diffusion coefficient, and v is the neutron velocity. denotes the infinite multiplication factor, and L is the diffusion length.

According to Eq. (8), the spatial domain meshes use the finite difference method to solve for the entire domain flux magnitude, and the numerically solved data are used as a test set to verify the model accuracy.

Two-dimensional rectangular geometry multi-group multi-material diffusion problem

In a nuclear reactor, the neutron transportation can be described by the multi-group diffusion theory. In this case, the fast- and hot-group neutron fluxes satisfy the following diffusion equations:(9)(10)where D₁ and D₂ are the diffusion coefficients of fast- and hot-group neutrons, respectively. ϕ₁ and ϕ₂ are the fast- and hot-group neutron flux densities, respectively. ν∑_f1 and ν∑_f2 are the neutron production cross sections of the fast- and hot-group neutrons, respectively. ∑_s1→2 is the fast group to hot group fission source term.

For pressurized water reactors, the example is divided into two different material regions, which is shown in Fig. 2. The material parameters for each region are listed in Table 2. Meanwhile, considering that the boundary energy between the fast- and hot-group is low enough, under such circumstances, no hot neutrons are directly produced by nuclear fission. As a result, χ₁=1, χ₂=0, and ∑_s2→1. Eq. (9) and Eq. (10) can be simplified to Eq. (11) and Eq. (12).(11)(12)

Fig. 2Full-screen View

Material Distribution [30]

Based on the multi-group diffusion theory, the distribution of the neutron flux in the core is obtained by iteratively solving the discretized diffusion equation. The obtained dataset was used as a test set in the experiment to evaluate the model prediction accuracy.

2D-IAEA benchmark problem

The 2D-IAEA PWR benchmark problem is a two-dimensional static problem with two neutron groups but without delayed neutron precursors [43]. It is modeled by the following two-dimensional two-group diffusion equations:(13)The reactor has a two-zone core containing 177 fuel assemblies with a width of 20 cm. The core is radially reflected by 20 cm of water. Owing to the symmetry along the x- and y-axes, this one-quarter reactor domain is denoted by Ω. It is composed of four sub-regions of different physical properties Ω_1,2,3,4. The reactor is shown in Fig. 3. Neumann boundary conditions were enforced at the left and bottom boundaries. The group constants for this problem are listed in Table 3.

Fig. 3Full-screen View

Geometric Layout of the 2D-IAEA Benchmark Problem [43]

PINN loss formulation

In PINN, considering the general form of a parameterized and nonlinear PDE,(14)where u represents the latent solution and Ω represents the solution domain. This formula can express PDEs in almost all the fields of mathematical physics.

N_f data points were sampled to measure the physical consistency. This type of loss is collectively referred to as the PDE loss. It has the following form:(15)Subsequently, the initial and boundary losses owing to the known initial and boundary conditions are introduced. N_i and N_b data points are collected to calculate Loss_Initial and Loss_Boundary, respectively. In addition, a few data points should be used for training.(16)(17)(18)Incorporating a small amount of labeled data provides information regarding the correct order of magnitude. This enables the PINN to better calibrate its predictions and prevent unrealistic or divergent solutions. The inclusion of such labeled data significantly contributes to the stability and accuracy of the training process for PINNs. By combining the above losses, networks can be penalized, and the training process is constrained effectively. This ensures that the obtained solutions adhere more closely to the fundamental laws of physics. Therefore, the final network loss formulation is as follows:(19)The weights of the different loss functions and number of data points used to calculate these losses significantly affect the convergence speed and accuracy of the model. The final sampling ratio was determined through multiple experiments. The sampling ratio for the PDE loss, boundary loss, initial loss, and data loss was approximately 30:10:10:1. Multiple losses can be balanced by adjusting the weights w. This prevents the network from prioritizing a single loss, particularly when significant differences in the magnitude between losses exist.

S-CNN architecture

When solving the neutron diffusion equation, the performance of the NN can be affected significantly by the vanishing gradient problem as the network depth increases. This problem can straightforwardly result in training failure and render the results unreliable. Moreover, it significantly limits the increase in the number of layers, reduces the expressive power of the network, and limits the prediction accuracy of the network. Inspired by the effective alleviation of the gradient vanishing problem in the ResNet architecture [44] in the image domain (which introduces the concept of residual learning), networks can learn residual mappings (the difference between the input and desired output). This facilitates the training of very deep neural networks. Thus, a skip-connection mechanism was introduced into the PINN architecture.

The input sampling features are the coordinates, namely, x, y, and t. These are independent. Therefore, a separate filter was applied to each feature in the experiment. A one-dimensional convolution was used as the basis for each network layer. Each hidden layer is expressed as follows:(20)where zl denotes hidden layer l between the input and output layers; Wl and bl are the weight and bias, respectively; and denotes the activation function (e.g., the Tanh function).

On this basis, skip connections were added between different layers. The corresponding hidden layer is expressed as(21)where n represents the skip distance, that is, the number of crossed hidden layers. To determine where a skip connection should be added, the contribution of each layer to the training of the entire network was measured by calculating the gradient norm for each layer. Here, the gradient norm was computed as in [40]:(22)For example, in the 10-layer network, the gradient contribution was calculated for each layer of the network (Fig. 4).

Fig. 4Full-screen View

(Color online) Gradient norm of each layer

When the gradient is significant, the parameters of the layer can be updated conveniently. A large gradient results in an unstable model, and a small gradient indicates that the layer may need to undergo more iterations to attain the optimal state and may encounter a gradient vanishing problem. Thus, the network increases the shallow gradient paths by establishing shallow-to-deep skip connections. This provides more gradient paths in the shallow layers and prevents vanishing gradients in the deep layers. The detailed 10-layer S-CNN architecture used in Sect. 4 is presented in Table 4. Here, B, C, and L represent the batch size, channel, and length, respectively.

In the experiment, the skip distance n was set to 3 for larger gradient propagation spans. It is the best parameter for shallow networks. When n < 3, the jump distance is insufficient. This may limit the network’s capability to learn complex physical fields and does not prevent the gradient vanishing problem. When n > 3, the number of jumping layers is excessively large. This may pose a challenge to gradient propagation and, thus, increase the difficulty of optimizing the network. Experiments were conducted on S-CNN models with varying depths, and high-precision results were obtained consistently. This validated the effectiveness of the residual module in addressing the gradient problem.

RAR Mechanism

Owing to specific regions with large gradients in the overall solution domain, if sampled uniformly, the network displays inadequate fitting of the local regions. This issue can be addressed by increasing the weights of specific sampling points (as shown in Eq. (23)) or by resampling using Eq. (24):(23)(24)where X_raw represents the original dataset, X_new represents the new dataset, X_f represents the set of sampling points with more significant PDE residuals, X_resample represents the set of resampling points, and w represents the penalty weight.

Meanwhile, referring to the grid division of conventional numerical computation methods, fine-grained samples should be collected from regions with large gradients to improve the prediction accuracy of the network in these regions.

Based on this, we consider using Eq. (24) to update our dataset, that is, introducing RAR [38] to improve the distribution of residual points during the training process of PINNs to address the bottleneck phenomenon that originates from the difficulty of reducing the PDE residuals in certain regions. This ultimately enhances the predictive accuracy of the model.

By selectively sampling more points in regions where the PDE residuals are significant, this approach enables the network to focus on challenging areas and adjust the sampling density accordingly. This, in turn, results in improved learning and prediction capabilities. The method is particularly effective for capturing the complex behavior of PDE solutions and identifying sharp gradient regions.

The pseudocode for the RAR method is as follows:

This algorithm divides the solution domain into α² subintervals with α subdivisions in both x- and t-directions. The Latin hypercube sampling (LHS) method [39] was used for random sampling. This ensured that sample points covered the entire space without clustering or bias.

As shown in Fig. 5, by adaptively adding points to regions with more significant residuals, more intensive sampling was conducted in one of the subdomains of the entire domain. This enabled the network to better capture the PDE solution’s behavior. This adaptive sampling approach improved the distribution of residual points and mitigated the bottleneck phenomenon. Ultimately, it enhanced the accuracy and performance of the model and accelerated its convergence.

Fig. 5Full-screen View

(Color online) PDE points distribution. (a) Before resampling; (b) After resampling

Proposed Framework: R²-PINN

Our study proposed an R²-PINN architecture. It combines the PINN structure, S-CNN, and RAR mechanisms to solve PDEs with improved accuracy and computational efficiency.

The structure of R²-PINN is shown in Fig. 6. In R²-PINN, S-CNN is used as the backbone of the network. It can better capture the features of the PDE solution and improve the training efficiency. The RAR mechanism is employed to balance Loss_PDE across regions of the domain. This ensures that the network focuses on regions with large errors and adaptively refines the mesh.

Fig. 6Full-screen View

(Color online) Overall architecture

By combining the PINN structure, S-CNN, and RAR mechanisms, the R²-PINN can effectively solve PDEs with improved accuracy and computational efficiency. It leverages the capability of deep learning and adaptive refinement to capture the complex features of a PDE solution. Thereby, it enables a more precise representation of the physical phenomena. Section 4 describes the verification of the superiority of the proposed R²-PINN network and each mechanism through a series of experiments.

Dataset preparation

S-CNN depth and kernel size ablation Experiment

All the parameters except for the base network were maintained consistent to compare the FCN and S-CNN architectures. The total number of sampled data points was 5000, with 3000 data points used to calculate Loss_PDE, 1000 data points sampled for Loss_Initial and 1000 data points sampled for Loss_Boundary. Network training does not involve feeding the data to calculate Loss_Data. The Tanh activation function and LBFGS optimizer were used. The Gaussian distribution random sampling method was used to initialize the network weights and biases. The number of hidden neurons per layer in the network was set to 26.

To investigate the impact of the layer number and kernel size on the model performance, ablation experiments were conducted on S-CNN networks with different depths and kernel sizes. The padding was set to zero for kernel size one and to one for kernel size three. The detailed experimental results are listed in Table 5 and Fig. 7.

Fig. 7Full-screen View

Comparative accuracy with different model parameters. (a) ϕ₀=ϕ₁; (b) ϕ₀=ϕ₂

In Table 5 and Fig. 7, Ω MSE refers to the mean squared error (MSE) score over the entire domain, and Ω₁ MSE refers to the MSE score at t=0.015 s, which is calculated separately to assess the capability of the model to extrapolate in time. The baseline refers to the results of FCN.

From Table 5 and Fig. 7, it can be observed that when the number of layers is less than 10, the accuracy shows a positive correlation with the number of layers. However, as the number of layers increases further, the accuracy does not increase. Hence, it is preferable to select fewer layers to make the NN train faster and achieve a high accuracy. This is because although the expressive capability of the neural network increases with the number of layers, beyond a certain threshold, an increase in the number of layers does not significantly improve the expressive capability of the network. An increase in the number of layers only causes an increase in the training time, and the model experiences overfitting problems. Moreover, because the input features are coordinates and not related to each other, the kernel size set to one can consider each channel as an independent feature to be addressed. This operation is more in line with reality. Therefore, a kernel size set to one can achieve a higher accuracy than a kernel size set to three. Based on the above analysis, the 10-layer S-CNN with a kernel size of one exhibited the best predictive accuracy with the minimum number of parameters. Consequently, this configuration was used in the subsequent experiments.

Ablation experiment on resampling parameters in S-CNN

According to the RAR algorithm, the resampling granularity was set to α to define the granularity of the subdomains in the solution domain. The solution domain was divided into α subintervals along the x and t dimensions. This resulted in α² subintervals. Following the RAR algorithm, in every 1000 epochs, the network calculates and compares the MSEs of each subdomain and performs resampling in the subinterval with the largest MSE. The number of resampling points is denoted by m.

The initial PDE sampling point was set to 3000 to prevent excessive sampling and training. Moreover, 2000 data points were sampled, with 1000 points to calculate the initial loss and 1000 points to calculate the boundary loss. The maximum number of PDE samples was set to 5000. Resampling was stopped when the total number of PDE samples attained 5000 after multiple iterations. The LHS method was used to resample.

Ablation experiments were conducted in two dimensions (resampling granularity and number of resampling points) while maintaining an identical S-CNN architecture (10 layers) and the other hyperparameters. The S-CNN network was compared with an FCN (baseline) using identical initial conditions, that is, ϕ₀=ϕ₁. The detailed experimental results are presented in Tables 6 and 7.

From Table 6 and Table 7, it is evident that under different sampling hyperparameters, our network consistently achieved a better Ω₁ MSE result. It was two orders of magnitude higher than that of the FCN baseline. When using the optimal hyperparameters, the R²-PINN achieved an accuracy of E-08or even E-09. This indicates that our method has a significant advantage in determining whether the flux values attain a steady-state at a specific instant. This advantage was demonstrated during the k_∞ search described in Sect. 4.4.

Furthermore, the results were obtained by adopting another initial condition, where ϕ₀=ϕ₂. A boxplot analysis of the resampling hyperparameters for all the results is shown in Fig. 8 and Fig. 9. When the number of resamples is 500 and the resample granularity is set to two, the model achieves the best accuracy.

Fig. 8Full-screen View

MSE Comparison Between Parameters. (a) m Comparison; (b) α Comparison

Fig. 9Full-screen View

Resample parameter 3D visualization

Given the introduction of the RAR mechanism, it is necessary to consider whether there is a significant time overhead. Hence, time comparisons were performed for different numbers of samples. The results are listed in Table 8. It can be observed that when the number of resamplings was set to 500, the training time required by the model reduced significantly compared with that when it was set to zero. By calculating the time consumed per cycle, we determined that this trend did not increase significantly. Thus, setting an appropriate number of resamplings effectively reduced the overall network training time. This shows that the RAR method involves an increase in the number of sampling points. This results in a different density of samples in each region, as well as a relatively large number of sampling points in complex regions. This contributes substantially to the convergence speed of the network.

The model uses the optimal resampling parameter configuration to predict the flux distribution under different k_∞ values where ϕ₀=ϕ₁. The error plot is shown in Fig. 10. Except for the relatively large errors at t=0 s, the errors in the other regions were relatively flat. Significantly, as the time increased, the errors increased marginally. In addition, Fig. 11 shows the distribution of the values predicted by the model. When t approaches zero, there is a specific deviation between the predicted result in the boundary region and zero. This guided us to appropriately increase the weights of the Loss_Boundary and Loss_Initial.

Fig. 10Full-screen View

(Color online) Error field. (a) k_∞=1.0001; (b) k_∞=1.0041

Fig. 11Full-screen View

(Color online) Predict result. (a) k_∞=1.0001; (b) k_∞=1.0041

To evaluate the performance of the model and measure the difference between the losses, Fig. 12 shows the training and testing losses, respectively. R²-PINN converged in 2000 epochs. In Fig. 12, Loss_PDE is larger than Loss_Boundary and Loss_Initial. This may result in Loss_PDE dominating the optimization process, whereas the other losses cannot be optimized effectively. To address this issue, w in Eq. (19) was set to 100 to impose higher penalties on the boundary and initial region error.

Fig. 12Full-screen View

Loss over iterations. (a) Training loss (k_∞=1.0001); (b) Training loss (k_∞=1.0041); (c) Test loss

Based on the experimental results in Sect. 4.2 and Sect. 4.3, an optimized S-CNN architecture with an optimized RAR mechanism was used to compose the R²-PINN for further use in searching critical parameters. This is discussed in Sect. 4.4 and Sect. 4.5.

k_∞ Search with R²-PINN

For a given geometric shape and volume of the reactor core, k_∞ and L² can be adjusted by modifying the reactor core size or modifying the composition of the materials within the reactor such that k_eff is one. When the system attains a steady state after a sufficient period, the neutron flux density follows the distribution described by Eq. (6), and the reactor is in a critical state.

In this experiment, the L² parameter size was not altered. Moreover, different networks were trained by adjusting k_∞ to predict the evolution of the neutron flux at this parameter. By continuously adjusting the value of k_∞, we attempted to adjust k_eff to one, thereby attaining the critical state. Initially, the parameter search range was set as [1.0001, 1.0041]. It has been verified that when k_∞ is 1.0001, k_eff<1. This indicates a subcritical state in which ϕ(x,t) decays exponentially with time t. When k_∞ is 1.0041, k_eff>1 indicates a supercritical state in which ϕ(x,t) increases continuously. The specific distributions are shown in Fig. 13.

Fig. 13Full-screen View

(Color online) ϕ Distribution. (a) ϕ₀=ϕ₁, k_∞=1.0001; (b) ϕ₀=ϕ₁, k_∞=1.0041; (c) ϕ₀=ϕ₂, k_∞=1.0001; (d) ϕ₀=ϕ₂, k_∞=1.0041

In the absence of delayed neutrons, when the system approaches criticality it converges to a critical state within a significantly short time (5-8 ms) regardless of the boundary conditions. Here, ϕ does not vary. At this point, the ϕ-distribution is a steady-state analytical solution. The parameter search interval is partitioned into n equal parts, yielding multiple k_∞ values. The network is trained for each value, calculates ϕ_t after a specific time tτ, and performs a search until ϕ_t approaches zero within a reasonable accuracy range. The process is illustrated in Fig. 14.

Fig. 14Full-screen View

Search flowchart [29]

In each network, Eq. (2) is used as the equation for Loss_PDE. Eq. (25) is used as the equation for Loss_Initial. The MSE between the boundary predicted values and zero is used to compute Loss_Boundary.

Using the automatic differentiation mechanism of the NN, the partial derivative of ϕ for t was calculated after obtaining the predicted values. Given the network’s capability to predict flux variations within a time interval of 0.015 s, the experiment used the last five time points to calculate and to determine whether the system was in a critical state [29].

Using Fig. 14, searches for k_∞ when the initial distribution follows ϕ₁ and ϕ₂ were conducted. Each search result was recorded. ϕ_t as a function of k_∞ is shown in Fig. 15. When n=2 and the grid method degenerates into a binary search, only approximately 20 iterations are required to obtain the search results with an accuracy level of E-05. The total runtime of the program was approximately 30 min.

Fig. 15Full-screen View

(Color online) Search result of ϕ_t. (a) ϕ₀=ϕ₁; (b) ϕ₀=ϕ₂

From Fig. 15, ϕ_t varies exponentially with k_∞. This further indicates that gradient-based methods such as gradient descent or Newton’s method can be evaluated to search for parameter values more rapidly and efficiently. Alternatively, using curve-fitting techniques to fit the scattered data and the intersection of the resulting curve with the x-axis yields the value of k_∞ at the critical state. This process can be completed in a small number of search iterations. The k_∞ value corresponding to the critical state is identified by progressively refining the interval. The search results are presented in Table 9. Among these, Δϕ at the last five time points is recorded to determine whether ϕ attains a steady state. Compared with the results of the FCN [29], the R²-PINN has a smaller Δϕ. This implies that the R²-PINN search has a higher accuracy than FCN. Furthermore, the search results were validated by generating the ϕ distribution for the obtained k_∞ values, as shown in Fig. 16.

Fig. 16Full-screen View

(Color online) Steady-state verification. (a) ϕ₀=ϕ₁; (b) ϕ₀=ϕ₂

From Fig. 16, the flux tends to stabilize as t increases. This implies that the system has attained a critical state. This indicates that our network can effectively search for optimal parameters corresponding to the critical state.

k_∞ Search efficiency improvement

When searching for k_∞, the initial search interval is large. The k_∞ value to be searched differs considerably from the critical state k value. Thus, the accuracy of prediction is not highly required. Only the prediction is used to compute ϕ_t to serve as a priori information for the next interval refinement. Therefore, during the training section, we examined whether the rate of variation in the fluxes converged at regular intervals. When ϕ_t converges, the training section is stopped, and the next k-value network training begins. This results in a faster parameter search without loss of accuracy. The equation for determining when to stop network iteration is as follows:(27)The parameter λ was set to 0.01. At every 200 iterations, Eq. (27) was computed to determine if the network had stopped iterating. The comparative results for ϕ₀=ϕ₁ are listed in Table 10. The early termination mechanism can significantly reduce the search time (0.56 times the original one).

Furthermore, as shown in Fig. 15, ϕ_t is exponentially related to k_∞. Therefore, it can be used to determine the critical value of k_∞ by fitting a quadratic function. To optimize the search algorithms, an experiment was conducted to compare three search methods: binary, grid, and quadratic fitting. Using R²-PINN with an early termination mechanism, the results are shown in Fig. 17.

Fig. 17Full-screen View

Comparison of search method

The quadratic fitting search method determined k_∞ rapidly in 250 s and attained an accuracy of E-04. Meanwhile, the grid search method determined k_∞ in 522 s and attained an accuracy of E-05. Different methods can be selected for parameter searches based on the actual requirement for time and accuracy.

Network prediction for different k_∞: R²-PINN’s robustness

The network was trained at different k_∞ values. The results were compared with the analytical solution generated to obtain the accuracy at each k_∞ value, as shown in Fig. 18, to validate the robustness of the R²-PINN.

Fig. 18Full-screen View

MSE Accuracy Comparison for Different k_∞

The robustness of R²-PINN is illustrated in Fig. 18. Across the tested values, the network consistently maintained a high accuracy of at least E-07. This indicated its capability to handle various scenarios and maintain reliable predictions. The observed differences in MSE across parameter values were relatively small. This further highlights the stability and effectiveness of the network. The robustness of the R²-PINN network is evident from its consistently high accuracy for different parameter values.

R²-PINN for solving a two-dimensional reactor diffusion equation for a single energy group

In this experiment, an 11-layer S-CNN network was used for training. The remaining hyperparameters were selected as described in Sect. 4.4.

To validate the generalizability of the models, an experiment was conducted between different k_∞. The results are shown in Table 11. Here, the network achieved an accuracy of E-05 under different parameters, and the MSE for the extrapolation region (that is, the region with t = 1 s) still achieved an accuracy of E-05. This demonstrates that the model performs well in terms of temporal extrapolation capability and can be used to search for parameters corresponding to the critical state.

Searching for the k_∞ parameter in the manner described in Sect. 4.4 yields k_∞ = 1.1378. Using the source iteration method, we determined that k_∞=1.1202. The calculation error is approximately 1.5%. The results of the prediction and reference truths are presented in Fig. 19. The error plots are shown in Fig. 20.

Fig. 19Full-screen View

(Color online) Comparison between result and truth. (a) Prediction result (t=0 s); (b) Prediction result (t=1 s); (c) Truth (t=0 s); (d) Truth (t=1 s)

Fig. 20Full-screen View

(Color online) Loss field. (a) t=0 s; (b) t=1 s

R²-PINN for solving two-dimensional rectangular geometry multi-group multi-material diffusion problem

In this two-group diffusion problem, each neutron group’s flux is predicted independently to enhance the accuracy. Specifically, we separately model the fast group neutron flux (ϕ₁) and hot group neutron flux (ϕ₂) using dual PINNs with a criticality factor of k_eff = 0.9693. Given the interconversion relationship between the fluxes of these two energy groups, dual PINNs are designed to share loss functions and are optimized sequentially. This shared loss structure ensures that the interaction between the energy groups is captured effectively while allowing for customized optimization paths for each flux. To train the models, a four-layer S-CNN was employed. The other PINN-related parameters, such as the sampling ratios, activation functions, and optimization strategies, remained consistent with those detailed in Sect. 4.4.

Using R²-PINN, the MSE of ϕ₁ attains 1.36×10^-6 and that of ϕ₂ attains 2.49×10^-7. To better illustrate the flux distribution, particularly the abrupt variations in the neutron flux at the material interfaces for different neutron groups, cross-sectional views of the neutron flux were adopted with slices extracted from the x-z and y-z planes. The results and reference solutions are presented in Figs. 21 and 22.

Fig. 21Full-screen View

(Color online) Predict result. (a) ϕ₁ Distribution in x-z plane; (b) ϕ₂ distribution in x-z plane; (c) ϕ₁ Distribution in y-z plane; (d) ϕ₂ distribution in y-z plane

Fig. 22Full-screen View

(Color online) Truth result. (a) ϕ₁ distribution in x-z plane; (b) ϕ₂ distribution in x-z plane; (c) ϕ₁ distribution in y-z plane; (d) ϕ₂ distribution in y-z plane

The error fields for each energy group are presented in Fig. 23. Owing to the small value of the flux, to clarify the error representation, we calculated the loss rate. This is shown in Fig. 23. The specific equations are given in Eq. (28).(28)

Fig. 23Full-screen View

(Color online) Loss field and loss rate field

The results show that we can effectively determine the flux distribution under the specified k_eff according to the steady-state equations. Additionally, according to the experiments presented in Sects. 4.4 and 4.7, we can determine the critical parameters according to the transient equations. Thus, the model can be adapted well to solve the two major problems of nuclear reactors, that is, solving for the fluxes and searching for the steady-state parameters.

R²-PINN for solving 2D-IAEA problem

Using a 6-layer R²-PINN structure, we incorporated k_eff as a parameter for iterative optimization in neural network training. The model utilized 18000 points for computing Loss_PDE, 500 points for Loss_Boundary, and 76 points for Loss_Data.

Finally, the relative error and relative L_∞ error were used to evaluate the accuracy of the R²-PINN. Herein, the relative L_∞ error is particularly important in the nuclear engineering domain. The specific equations are expressed in Eqs. (29) and (30). The prediction results and reference truth are presented in Fig. 24. The absolute error plots are shown in Fig. 25.(29)(30)

Fig. 24Full-screen View

(Color online) Comparison between results and reference. (a) R²-PINN result of ϕ₁; (b) R²-PINN result of ϕ₂; (c) FreeFem++ result of ϕ₁; (d) FreeFem++ result of ϕ₂

Fig. 25Full-screen View

(Color online) Absolute loss field. (a) Absolute error of ϕ₁; (b) Absolute error of ϕ₂

Considering the engineering acceptance criteria for the 2D IBP, the flux calculation error in fuel assemblies with a relative flux higher than 0.9 should be less than 5%. In fuel assemblies with a relative flux less than 0.9, the flux calculation error should be less than 8%. In addition, the relative error of k_eff should be less than 0.005 [47]. The predicted results shown in Table 12 satisfy these acceptance criteria. This indicates that R2-PINN also holds practical engineering application value.

The strategy for hyperparameter selection

The hyperparameter selection strategy was formulated meticulously to balance the model expressiveness, training stability, and prediction accuracy across the computational domain. This subsection explains the important hyperparameter selections for neural networks.

First, a 10-layer network was selected to achieve an optimal tradeoff between complexity and expressiveness. Compared with shallower networks, deeper networks more effectively approximate complex nonlinear functions. This is typical in physical modeling. As shown in Table 5, our evaluation across configurations indicated that a 10-layer S-CNN structure consistently outperforms the others. This makes it the optimal option. This selection method was similarly applied to different benchmark cases to ensure robustness.

Second, the Tanh activation function was selected over ReLU and Sigmoid because of its smoother nonlinearity, which is critical for approximating continuous functions in physical problems. Although ReLU is effective in mitigating vanishing gradients, it has insufficient stability in capturing smooth physical fields and generally encounters difficulty with highly nonlinear PDEs. In contrast, Tanh enhances numerical stability. This makes it more suitable for PINNs.

Third, the LBFGS optimization algorithm was selected for its faster convergence in high-dimensional problems and capability to prevent gradient explosions. Compared with first-order optimizers such as Adam and SGD, LBFGS provides a second-order approximation. This results in a more stable training and better performance, particularly in data-limited scenarios. During testing, Adam and SGD were vulnerable to early convergence and local minima, thereby yielding suboptimal predictions. Meanwhile, LBFGS provided better stability and training completeness.

Finally, the selection of 26 hidden neurons per layer achieved a balance between the model capacity and generalization. Excessively few neurons result in underfitting, whereas excessively many neurons risk overfitting and instability. Through extensive experimentation, 26 neurons per layer were observed to provide an optimal trade-off. This ensured accurate and stable predictions.