JFNK method is one of the inexact Newton-Krylov methods, which are the combinations of Newton’s methods for solving non-linear equations, Krylov subspace methods for solving the linear Newton correction equations, and “Jacobian-free” technology for approximately solving Jacobian-vector product (Knoll and Keyes, 2004). The non-linear equations can be generally written by F x = 0 . The Newton iteration step (where n is the index) is to solve the discrete equations as J a x n δ x n = − F x n , x n + 1 = x n + δ x n , (1) where x is the vector of global solution variables, F represents the non-linear discrete system, and J a represents the Jacobian matrix. For Krylov methods, the Jacobian matrix needs not to be formed explicitly because the product of Jacobian matrix and a vector ( v ) rather than the Jacobian matrix is required. Using the first-order finite difference approximation, the Jacobian-vector product can be solved by J a x n v = F x n + ε v − F x n ε (2) where ε is a perturbation parameter which can be specified by x n and v .

For the parallel JFNK method, only the residual functions, the dot products and the norms of vectors are considered to be evaluated synchronously in multiple processors, which are easily developed in the parallel framework. In addition, residual functions are evaluated from the parallel Sn method based on KBA sweeping algorithm as shown in the following sections.

The multi-group neutron Boltzmann transport equations for the k-eigenvalue problem can be written as Ω ^ n ∙ ∇ ψ g , n r → + Σ t , g r → ψ g , n r → = ∑ g ′ = 1 G Σ s , g ′ → g r → ϕ g ′ r → + χ g k e f f ∑ g ′ = 0 G ν Σ f , g ′ r → ϕ g ′ r → , (3) ϕ g r → = ∑ n = 1 N ω n ψ g , n r → , (4) where ψ g , n r → = Angular neutron flux, ϕ g r → = Scalar neutron flux, k e f f = Effective multiplication factor (k-eigenvalue), Σ t , g = Macroscopic total cross section, Σ s , g ′ → g = Macroscopic scattering cross section, χ g = Fission spectrum, ν Σ f , g = Macroscopic production cross section, g = Energy group index and G = Total number of energy group, n = Discrete angle index and N = Total number of discrete angles, r → = Spatial position vector and Ω ^ n = Unit direction vector for discrete angle n , ω n = Quadrature weight for discrete angle n .

Eq. 3 can be written in operator notation form as L ψ = M S ϕ + 1 k e f f M χ f ϕ , (5) where ψ and ϕ are vectors which contain the angular fluxes and scalar fluxes, respectively. L denotes the sum of the streaming operator and removal operator. M is the scalar-to-discrete operator and D is the discrete-to-scalar operator; in the case of isotropic scattering, D represents the action which integrates the angular fluxes to the scalar fluxes as shown is Eq. 4, i.e., ϕ = D ψ . S is the matrix of group-to-group scattering cross sections, f is the matrix of production cross sections and χ is the fission spectrum matrix. The parallel transport sweeping based on the KBA algorithm with the spatial domain decomposition (Baker, 2017) is employed to directly inverse L into L − 1 , which marches across the spatial domain cells in the direction of neutron travel.

Multiplied by D L − 1 , Eq. 5 becomes ϕ = D L − 1 M S ϕ + 1 k e f f D L − 1 M χ f ϕ ⇒ I − D L − 1 M S ϕ = 1 k e f f D L − 1 M χ f ϕ . (6)

It is an eigenvalue problem, whose solution is to get an eigenpair vector ϕ T , k e f f T . Traditionally, it can be solved by power iteration method (PI) which includes two-level nested iterations. The outer iteration can be written as ϕ l + 1 = I − D L − 1 M S − 1 1 k e f f l D L − 1 M χ f ϕ l , (7) k e f f l + 1 = c ^ || f ϕ l + 1 || 1 , (8) where l is the outer iterative index, ∙ 1 denotes the L-1 norm. c ^ is a fixed value during iteration and it can be assigned form the iterative initial value as c ^ ≡ k e f f 0 || f ϕ 0 || 1 . (9)

For each outer iteration, Eq. 7 is equivalent to a fixed source problem, which can generally be calculated using fixed point iteration method (which is the inner iteration) as ϕ m + 1 , l = D L − 1 M S ϕ m , l + 1 k e f f l χ f ϕ l , (10) where m is the inner iterative index. In this paper, this power iterative strategy is called comePSn method. To efficiently get the simultaneous solution of the non-linear system using the parallel JFNK method, the k-eigenvalue and the scalar fluxes (rather than the angular fluxes) are finally chosen as the global solution variables, and x = ϕ T , k e f f T .

Meanwhile, x is also the minimal subset of all the variables for the discrete non-linear system based on the parallel Sn method, which indicates that the JFNK solution has the least number of global solution variables and can obtain good numerical efficiency.

Then, according to the PI strategy for the parallel Sn method as shown in Eq. 7 and Eq. 8, the corresponding residual functions for the parallel JFNK method can be easily expressed as F x ≡ F ϕ x F k x = ϕ − I − D L − 1 M S − 1 1 k e f f D L − 1 M χ f ϕ k e f f − c ^ || f I − D L − 1 M S − 1 1 k e f f D L − 1 M χ f ϕ || 1 , (11) where F ϕ x and F k x denotes the residual functions for ϕ and k e f f , respectively. However, as mentioned above, the solution of the I − D L − 1 M S − 1 operation needs nested inner iterations with a certain amount of transport sweeping steps (multiple calculation of L − 1 ), which can lead to significantly expensive cost when evaluating residual functions. In fact, only one transport sweeping step is enough for the evaluation of the residual functions for JFNK. As a result, in this paper, the residual functions are constructed according to Eq. 6 and 8 as F ϕ x = ϕ − D L − 1 M S ϕ − 1 k e f f D L − 1 M χ f ϕ = ϕ − D L − 1 M S + 1 k e f f χ f ϕ , (12) F k x = k e f f − c ^ || f D L − 1 M S + 1 k e f f χ f ϕ || 1 , (13) where the evaluation is equivalent to a “flattened” power iterative process (Gill, 2009) including a single outer iteration without nested inner iterations as shown in Eq. 10. In this paper, the combined strategy of the parallel JFNK method and the parallel Sn method is called comePSn_JFNK method.

In this paper, the comePSn and comePSn_JFNK codes are parallelly designed on spatial domain decomposition with the Message Passing Interface (MPI) standard. For the comePSn_JFNK code, the message communication between MPI ranks is mainly required for the evaluation of the residual functions, the dot products and the norms. The parallel dot products and norms of vectors can be easily calculated by using the MPI_Allreduce function. And the evaluation of residual functions employs the KBA algorithm, which can directly calculate the inversion of L in the parallel Sn framework.

KBA is constructed on the (N-1)-dimensional spatial domain decomposition for N-dimensional models. For 3-D rectangular coordinate system, the spatial grid is divided into domains on the XY layout of CPU cores (or MPI ranks). Each domain is further decomposed into computational chunks, which are solved in the “wavefront” ordering. Each computational chunk, in a single octant with N/8 discrete angles for a single energy group, is solved after receiving the incoming angular fluxes from the upstream chunks (or from the boundary conditions), and then the outer angular fluxes are sent to the downstream chunks if required. The “simultaneous in angle, successive in quadrants” pipelining method described in the reference (Baker and Koch, 1998) are chosen in this paper. In addition, to reduce the idle time, the sweeps start at all the four corners of the spatial grid at the same time (Bailey and Falgout, 2009).

As mentioned above, the KBA transport sweeping process can be easily applied to evaluate the residual functions for the parallel JFNK framework. In this work, the evaluation of the residual functions is a three-step process. 1) Map x to ϕ and k e f f in each MPI rank. 2) Calculate the new scalar fluxes ϕ ∼ by using the KBA algorithm where ϕ ∼ = D L − 1 M S + 1 / k e f f χ f ϕ , and then get the new k-eigenvalue as k ∼ e f f = c ^ || f ϕ ∼ || 1 . 3) Calculate the residual functions F x according to Eqs 12, 13, where F x = ϕ − ϕ ∼ k e f f − k ∼ e f f .

Based on the above methodologies, the parallel comePSn_JFNK and comePSn codes are respectively developed in the unified COME using C++ language to study the numerical accuracy, the efficiency and the parallel performance for the 3-D/2-D pin-by-pin neutron transport models. The iterative strategies and the flowchart are shown in Figure 1. The good initial guess x 0 of the comePSn_JFNK code comes from the solution after 5 power iterative steps. The CPU time and iterative properties of the comePSn code after 5 iterative steps are also compared with those of the comePSn_JFNK code. In addition, the number of inner steps is fixed as one per outer step for the comePSn code, which ensures that a single power iterative step for the comePSn code is equivalent to one calling calculation of residual functions for the comePSn_JFNK code. Then, the relative convergence criterion for both comePSn_JFNK and comePSn codes is F x n 2 F x 0 2 < 10 − 6 , (14) where ∙ 2 denotes the L-2 norm and F x 0 denotes the residuals of the iterative initial value. In each Newton step ( n ), the adaptive convergence criterion with the Eisenstat-Walker forcing term η n (Eisenstat and Walker, 1996) is J a x n δ x n + F x n ‖ 2 < η n F x n ‖ 2 , (15) η n = 0.1 F x n F x n − 1 2.0 . (16)

The restarted GMRES (Generalized Minimal REsidual) method (Saad and Schultz, 1986) is adopted to solve Eq. 15 and the maximum Krylov subspace dimension is set to 20. At last, the global convergent backtracking algorithm (Eisenstat and Walker, 1994) is employed to ensure the global convergence of come_PSn_JFNK.

KAIST-3A benchmark problems are published by KAIST Nuclear Reactor Analysis and Particle Transport laboratory (Cho, 2000), whose model is a simplified PWR core. The quarter symmetric radial geometry is illustrated in Figure 2. The core consists of two types of fuel assemblies, Uranium (UOX, including UOX-1 with 2.0% enrichment and UOX-2 with 3.3% enrichment) and Plutonium (MOX), which are surrounded by the baffle and reflector. Each fuel assembly follows a 17 × 17 lattice design of 21.42 cm width with a 1.26 cm pin pitch, which includes 264 fuel rods (or burnable absorber rods), 24 guide tubes (or control rods inserted) and one instrument guide tube as shown in Figure 3. There are two cases: all rods out (ARO) and all rods in (ARI), based on the configurations with or without the control rods (CR) inserted into the four UOX-2 fuel assemblies. Seven-group macroscopic cross sections are provided in the benchmark report, including UO₂ fuel, MOX fuel, guide tube, BA (Burnable absorber) rod, control rod, baffle and reflector.

At first, the 2-D quarter symmetric models are simulated by comePSn_JFNK and comePSn codes. The pin cell is divided into 1 × 1, 2 × 2, 4 × 4 and 32 × 32 rectangular meshes, respectively.

Table 1 and Figure 4 respectively show the comparison of the k-eigenvalue and the relative errors (%) of power density for ARO and ARI with different mesh sizes per pin using comePSn_JFNK and comePSn codes. It can be seen that there is no obvious difference between the results from the two codes. To further analyze the accuracy of the two codes with different mesh sizes per pin, the relative errors (%) of pin-wise power density on 1 × 1, 2 × 2 and 4 × 4 meshes per pin are presented in Figure 5 compared with the results on 32 × 32 meshes per pin. For 1 × 1 mesh per pin, the maximum value of the relative errors reaches 12% for ARO and 18% for ARI, and the errors of the k-eigenvalues are about 200 pcm (for ARO) and 750 pcm (for ARI). While for 2 × 2 and 4 × 4 meshes per pin, the errors of the power density and the k-eigenvalues are respectively less than 2% and 50 pcm, and it indicates that the solution on 2 × 2 meshes per pin can obtain acceptable numerical accuracy.

To further analyze the numerical efficiency of comePSn_JFNK and comePSn codes, the CPU times and iteration numbers are presented in Table 2. It should be noted that the CPU times do not include the cost of the initialization, reading the input files and writing the output files. The results indicate that the different mesh sizes per pin make no obvious difference on the convergence rate for both comePSn_JFNK and comePSn codes. The computational efficiency of the comePSn_JFNK code is about 6–8 times as much as the comePSn code for ARO and about 8–10 for ARI.

Then the KAIST-3A 2-D quarter symmetric models are extend to 3-D full-core models with the same radial configuration, and axial reflectors are set at the top/bottom of the core. The active fuel length is 365.76 cm and the width of axial reactors is 21.42 cm. The control rods can be inserted from the top of the upper reflector to the bottom of the active fuel region in the UOX-2 (CR) assemblies. ARO and ARI cases are still simulated by comePSn_JFNK and comePSn codes.

The 3-D full-core model is divided into 340 × 340 × 340 meshes. Specifically, one pin cell is divided into 2 × 2 meshes in the radial direction; in the axial direction, the active region is divided into 300 meshes and the top/bottom reflectors are respectively divided into 20 meshes. Both comePSn_JFNK and comePSn codes solve the problems using 100 CPU cores. Figure 6 presents the relative errors (%) of radial average power density between the two codes and Figure 7 shows the radial average power density and 3-D power density distribution in the active regions calculated the comePSn_JFNK code. The numerical results of radial power density and the k-eigenvalues from the two codes also agree well with each other. The k-eigenvalues are respectively 1.12612 for ARO and 0.97063 for ARI.

To further analyze the efficiency of comePSn_JFNK and comePSn codes for 3-D models, the iterative properties and CPU times (on 100 CPU cores) are listed in Table 3. Compared with the comePSn code, the speedups of the comePSn_JFNK code are 12.41 for ARO and 14.64 for ARI. It is worth noting that the acceleration rates of the comePSn_JFNK code for 3-D models are higher than those for 2-D models, which shows the advantage of numerical efficiency of the comePSn_JFNK code for the complicated 3-D neutron transport problems.

In addition, the parallel efficiencies of the comePSn_JFNK code for the 3-D KAIST-3A ARO case using 4, 16, 100 and 400 CPU cores are shown in Figure 8. The “Solving” curve means the parallel efficiency of the total CPU time and the “Sweeping” curve indicates the parallel efficiency of KBA sweeping. The efficiencies of “Solving” and “Sweeping” using 400 CPU cores are respectively about 63% and 68%, which can be further improved by optimizing the parallel performance of the Krylov subspace methods and the strategies of the KBA algorithm.

A series of 3-D/2-D PWR full-core MOX/UOX pin-by-pin models with different insertion positions of control rods are simulated to further study the numerical properties of comePSn_JFNK and comePSn codes. The reactor core also consists of UOX and MOX fuel assemblies surrounded by reflector. As shown in Figure 9, each fuel assembly is 21.42 cm wide and is made up of 17 × 17 lattice with a 1.26 cm pin pitch, which includes 264 fuel pins, 24 guide tubes (or control rods inserted) and one instrument tube for fission chamber. The UOX and MOX assemblies are distributed in a checker-board configuration and four core models (A, B, C and D) are designed with different configurations of the UOX assemblies with control rods inserted (UOA) as illustrated in Figure 10. The active fuel length is 365.76 cm and there are 21.42 cm high axial reflectors at the top/bottom of the core just as the 3-D KAIST-3A models. It should be noted that the control rods are inserted form the top of the upper reflector to the bottom of the active fuel regions in the models B, C and D; and the control rods are still in the upper reflector if withdrawn in the models A, C and D. In this paper both 2-D and 3-D models are simulated by comePSn_JFNK and comePSn codes. The seven-group cross sections for pin cells are chosen form the C5G7 benchmark report (Smith et al., 2005). The pin cells are divided into 2 × 2 rectangular meshes in the radial direction; in the axial direction, the active fuel region is divided into 580 meshes and the top/bottom reflectors are respectively divided into 30 meshes.

Table 4 shows the results of the k-eigenvalues from comePSn_JFNK and comePSn codes and the reference solution is from the paper (Zhou, 2022a). It can be observed that the errors of k-eigenvalues from the two codes are about 150 pcm compared with the reference. Figure 11 shows the power density distribution from the comePSn_JFNK code and Figure 12 shows the corresponding relative errors (%) compared with the comePSn code. It can be observed that the power density distribution from the two codes agree well with each other, and even in the model B where the control rods are all inserted in, the maximum of the errors is less than 0.002%.

To analyze the numerical efficiency, Table 5 presents the comparation of the computational cost for the 2-D models from comePSn_JFNK and comePSn codes (using 1 core). Compared with the comePSn code, the speedups of the comePSn_JFNK code are respectively 15.40, 16.84, 29.90, and 29.67 for models A–D, which indicates that the acceleration rate of the comePSn_JFNK code for the asymmetric models C and D is higher than that for the symmetric models A and B. In addition, there is obvious difference of the iterative steps between the four models from the comePSn code, which range from 7093 to 27041. While the comePSn_JFNK code takes similar non-linear iterative steps (range from 5 to 10) due to the robust convergence of the JFNK methods.

To further test the computational properties of the comePSn_JFNK code for 3-D models, the radial average power density and the 3-D power density distribution calculated by the comePSn_JFNK code using 1024 CPU cores are presented in Figure 13, and the k-eigenvalues and the computational cost compared with the comePSn code is shown in Table 6. It can be seen that the acceleration rate of the comePSn_JFNK code are respectively 20.48, 22.57, 24.93 and 26.45 for models A– D, which indicates that the comePSn_JFNK code still has significant efficiency advantage for the complicated 3-D full-core pin-by-pin models with different control rod distributions.