Parallel Jacobian-free Newton Krylov discrete ordinates method for pin-by-pin neutron transport models

A parallel Jacobian-Free Newton Krylov discrete ordinates method (comePSn_JFNK) is proposed to solve the multi-dimensional multi-group pin-by-pin neutron transport models, which makes full use of the good efficiency and parallel performance of the JFNK framework and the high accuracy of the Sn method for the large-scale models. In this paper, the k-eigenvalue and the scalar fluxes (rather than the angular fluxes) are chosen as the global solution variables of the parallel JFNK method, and the corresponding residual functions are evaluated by the Koch–Baker–Alcouffe (KBA) algorithm with the spatial domain decomposition in the parallel Sn framework. Unlike the original Sn iterative strategy, only a “flattened” power iterative process which includes a single outer iteration without nested inner iterations is required for the JFNK strategy. Finally, the comePSn_JFNK code is developed in C++ language and, the numerical solutions of the 2-D/3-D KAIST-3A benchmark problems and the 2-D/3-D full-core MOX/UOX pin-by-pin models with different control rod distribution show that comePSn_JFNK method can obtain significant efficiency advantage compared with the original power iteration method (comePSn) for the parallel simulation of the large-scale complicated pin-by-pin models.

1 Introduction

The efficient parallel algorithms and the acceleration methods of discrete ordinates methods (Sn) (Carlson, 1953) for multi-dimensional large-scale neutron transport models are research hotspots in the reactor simulation field. In the past few decades, many parallel Sn neutron transport codes are developed by using overloading spatial domain decomposition (Bailey and Falgout, 2009) or Koch–Baker–Alcouffe (KBA) algorithms (Baker and Koch, 1998), including DENOVO (Evans et al., 2010), PARTISN (Alcouffe et al., 2005), NECP-hydra (Xu et al., 2018), ARES (Zhang et al., 2017) and so on. To accelerate the calculation of the original Sn iterative strategies based on the power iteration methods (PI), a batch of numerical methods can be employed, such as Wielandt shifts methods, Chebyshev acceleration methods, coarse mesh finite difference methods, diffusion synthetic acceleration method. In this paper, the parallel Jacobian-free Newton-Krylov (JFNK) (Knoll and Keyes, 2004) methods are adopted.

The JFNK methods have good computational efficiency due to the fast and robust convergence, which are widely used to solve the atmospheric convection problem (Hossain and Alam, 2012), the sea ice simulation (Yaremchu and Panteleev, 2022), two-phase flow calculation (Hajizadeh et al., 2018), the neutron diffusion/transport models (Gill, 2009; Knoll et al., 2011), thermal hydraulics simulation (Esmaili et al., 2020) and multi-physics coupled problems (Walker et al., 2019). Meanwhile, some popular multi-physics coupled environments (or platforms), such as MOOSE (Gaston et al., 2009), LIME (Pawlowski et al., 2011) and VERA (Turner et al., 2016), employ the JFNK framework to solve the complicated reactor coupled systems. Recently, we have also developed an efficient JFNK framework based on coarse mesh finite difference and nodal expansion method in the COupling Multiphysics Environment (COME) to solve the steady/transient neutronics and neutronics-thermal hydraulic coupled problems (Zhou, 2022a; Zhou, 2022b; Zhou et al., 2022; Zhou, 2023). To take advantages of the good efficiency and parallel performance of the JFNK framework in COME, in this paper, the parallel JFNK method is combined with the parallel Sn method to solve the multi-dimensional large-scale neutron transport models.

Therefore, a parallel Jacobian-free Newton Krylov discrete ordinates method (comePSn_JFNK) is proposed to solve the multi-dimensional multi-group pin-by-pin neutron transport k-eigenvalue problems, which makes full use of the fast and robust convergence of the parallel JFNK framework and the high accuracy of the Sn method based on the KBA algorithm. The key to implement the JFNK framework into the Sn algorithm is to choose the k-eigenvalue and the scalar fluxes (rather than the angular fluxes) as the global solution variables for parallel JFNK, which are the minimal subset of the non-linear system. The corresponding residual functions are directly constructed from the neutron transport equations and evaluated based on the parallel Sn transport sweeping process by using a “flattened” power iterative process including a single outer iterative step without nested inner iterations.

This paper is organized as follows. The detailed theories of the comePSn_JFNK method are presented in Section 2, including the parallel JFNK method, the parallel Sn method, the residual function evaluation and the solution strategies. In Section 3, the numerical results for the 2-D/3-D full-core pin-by-pin models are analyzed to test the accuracy and efficiency of the comePSn_JFNK method. The conclusions are given in Section 4.

2 Methodologies of the comePSn_JFNK method

2.1 Parallel JFNK

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 $\mathbf{F}\left(\mathbf{x}\right)=0$. The Newton iteration step (where n is the index) is to solve the discrete equations as

$\begin{array}{c}\mathbf{J}\mathbf{a}\left({\mathbf{x}}^{\left(\mathit{n}\right)}\right)\delta {\mathbf{x}}^{\left(\mathit{n}\right)}=-\mathbf{F}\left({\mathbf{x}}^{\left(\mathit{n}\right)}\right),\\ {\mathbf{x}}^{\left(\mathit{n}+1\right)}={\mathbf{x}}^{\left(\mathit{n}\right)}+\delta {\mathbf{x}}^{\left(\mathit{n}\right)},\end{array}(1)$

where $\mathbf{x}$ is the vector of global solution variables, $\mathbf{F}$ represents the non-linear discrete system, and $\mathbf{J}\mathbf{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 ($\mathbf{v}$) rather than the Jacobian matrix is required. Using the first-order finite difference approximation, the Jacobian-vector product can be solved by

$\mathbf{J}\mathbf{a}\left({\mathbf{x}}^{\left(\mathit{n}\right)}\right)\mathbf{v}=\frac{\mathbf{F}\left({\mathbf{x}}^{\left(\mathit{n}\right)}+\epsilon \mathbf{v}\right)-\mathbf{F}\left({\mathbf{x}}^{\left(\mathit{n}\right)}\right)}{\epsilon }(2)$

where $\epsilon$ is a perturbation parameter which can be specified by ${\mathbf{x}}^{\left(\mathit{n}\right)}$ and $\mathbf{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.

2.2 Evaluation of residual functions from parallel Sn

The multi-group neutron Boltzmann transport equations for the k-eigenvalue problem can be written as

${\widehat{\mathrm{\Omega }}}_{\mathit{n}}\bullet \nabla {\psi }_{\mathit{g},\mathit{n}}\left(\overrightarrow{\mathit{r}}\right)+{\Sigma }_{\mathit{t},\mathit{g}}\left(\overrightarrow{\mathit{r}}\right){\psi }_{\mathit{g},\mathit{n}}\left(\overrightarrow{\mathit{r}}\right)={\displaystyle \sum }_{{\mathit{g}}^{\prime }=1}^{\mathit{G}}{\Sigma }_{\mathit{s},{\mathit{g}}^{\prime }\to \mathit{g}}\left(\overrightarrow{\mathit{r}}\right){\mathit{\varphi }}_{{\mathit{g}}^{\prime }}\left(\overrightarrow{\mathit{r}}\right)+\frac{{\chi }_{\mathit{g}}}{{\mathit{k}}_{\mathit{e}\mathit{f}\mathit{f}}}{\displaystyle \sum }_{{\mathit{g}}^{\prime }=0}^{\mathit{G}}\nu {\Sigma }_{\mathit{f},{\mathit{g}}^{\prime }}\left(\overrightarrow{\mathit{r}}\right){\mathit{\varphi }}_{{\mathit{g}}^{\prime }}\left(\overrightarrow{\mathit{r}}\right),(3)$

${\mathit{\varphi }}_{\mathit{g}}\left(\overrightarrow{\mathit{r}}\right)={\displaystyle \sum }_{\mathit{n}=1}^{\mathit{N}}{\omega }_{\mathit{n}}{\psi }_{\mathit{g},\mathit{n}}\left(\overrightarrow{\mathit{r}}\right),(4)$

where ${\psi }_{\mathrm{g},\mathit{n}}\left(\overrightarrow{\mathit{r}}\right)$ = Angular neutron flux, ${\mathit{\varphi }}_{\mathrm{g}}\left(\overrightarrow{\mathit{r}}\right)$ = Scalar neutron flux, ${\mathrm{k}}_{\mathrm{e}\mathrm{f}\mathrm{f}}$ = Effective multiplication factor (k-eigenvalue), ${\Sigma }_{\mathrm{t},\mathrm{g}}$ = Macroscopic total cross section, ${\Sigma }_{\mathrm{s},{\mathrm{g}}^{\prime }\to \mathrm{g}}$ = Macroscopic scattering cross section, ${\chi }_{\mathrm{g}}$ = Fission spectrum, $\nu {\Sigma }_{\mathrm{f},\mathrm{g}}$ = Macroscopic production cross section, $\mathrm{g}$ = Energy group index and $\mathrm{G}$ = Total number of energy group, $\mathit{n}$ = Discrete angle index and $\mathrm{N}$ = Total number of discrete angles, $\overrightarrow{\mathrm{r}}$ = Spatial position vector and ${\widehat{\Omega }}_{\mathit{n}}$ = Unit direction vector for discrete angle $\mathit{n}$, ${\omega }_{\mathit{n}}$ = Quadrature weight for discrete angle $\mathit{n}$.

Eq. 3 can be written in operator notation form as

$\mathbf{L}\mathbf{\psi }=\mathbf{M}\mathbf{S}\mathbf{\varphi }+\frac{1}{{\mathit{k}}_{\mathit{e}\mathit{f}\mathit{f}}}\mathbf{M}\mathbf{\chi }\mathbf{f}\mathbf{\varphi },(5)$

where $\mathbf{\psi }$ and $\mathbf{\varphi }$ are vectors which contain the angular fluxes and scalar fluxes, respectively. $\mathbf{L}$ denotes the sum of the streaming operator and removal operator. $\mathbf{M}$ is the scalar-to-discrete operator and $\mathbf{D}$ is the discrete-to-scalar operator; in the case of isotropic scattering, $\mathbf{D}$ represents the action which integrates the angular fluxes to the scalar fluxes as shown is Eq. 4, i.e., $\mathbf{\varphi }=\mathbf{D}\mathbf{\psi }$. $\mathbf{S}$ is the matrix of group-to-group scattering cross sections, $\mathbf{f}$ is the matrix of production cross sections and $\mathbf{\chi }$ 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 $\mathbf{L}$ into ${\mathbf{L}}^{-1}$, which marches across the spatial domain cells in the direction of neutron travel.

Multiplied by $\mathbf{D}{\mathbf{L}}^{-1}$, Eq. 5 becomes

$\mathbf{\varphi }=\mathbf{D}{\mathbf{L}}^{-1}\mathbf{M}\mathbf{S}\mathrm{\varphi }+\frac{1}{{\mathit{k}}_{\mathit{e}\mathit{f}\mathit{f}}}\mathbf{D}{\mathbf{L}}^{-1}\mathbf{M}\mathbf{\chi }\mathbf{f}\mathbf{\varphi }\Rightarrow \left(\mathbf{I}-\mathbf{D}{\mathbf{L}}^{-1}\mathbf{M}\mathbf{S}\right)\mathbf{\varphi }=\frac{1}{{\mathit{k}}_{\mathit{e}\mathit{f}\mathit{f}}}\mathbf{D}{\mathbf{L}}^{-1}\mathbf{M}\mathbf{\chi }\mathbf{f}\mathbf{\varphi }.(6)$

It is an eigenvalue problem, whose solution is to get an eigenpair vector ${\left({\mathrm{\varphi }}^{\mathrm{T}},{\mathrm{k}}_{\mathrm{e}\mathrm{f}\mathrm{f}}\right)}^{\mathrm{T}}$. Traditionally, it can be solved by power iteration method (PI) which includes two-level nested iterations. The outer iteration can be written as

${\mathbf{\varphi }}^{\left(\mathit{l}+1\right)}={\left(\mathbf{I}-\mathbf{D}{\mathbf{L}}^{-1}\mathbf{M}\mathbf{S}\right)}^{-1}\frac{1}{{{\mathit{k}}_{\mathit{e}\mathit{f}\mathit{f}}}^{\left(\mathit{l}\right)}}\mathbf{D}{\mathbf{L}}^{-1}\mathbf{M}\mathbf{\chi }\mathbf{f}{\mathbf{\varphi }}^{\left(\mathit{l}\right)},(7)$

${{\mathit{k}}_{\mathit{e}\mathit{f}\mathit{f}}}^{\left(\mathit{l}+1\right)}=\widehat{\mathit{c}}||\mathbf{f}{\mathbf{\varphi }}^{\left(\mathit{l}+1\right)}{||}_1,(8)$

where $\mathit{l}$ is the outer iterative index, ${‖\bullet ‖}_1$ denotes the L-1 norm. $\widehat{\mathrm{c}}$ is a fixed value during iteration and it can be assigned form the iterative initial value as

$\widehat{\mathit{c}}\equiv \frac{{{\mathit{k}}_{\mathit{e}\mathit{f}\mathit{f}}}^{\left(0\right)}}{||\mathbf{f}{\mathrm{\varphi }}^{\left(0\right)}{||}_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

${\mathbf{\varphi }}^{\left(\mathit{m}+1,\mathit{l}\right)}=\mathbf{D}{\mathbf{L}}^{-1}\mathbf{M}\left(\mathbf{S}{\mathbf{\varphi }}^{\left(\mathit{m},\mathit{l}\right)}+\frac{1}{{{\mathit{k}}_{\mathit{e}\mathit{f}\mathit{f}}}^{\left(\mathit{l}\right)}}\mathbf{\chi }\mathbf{f}{\mathbf{\varphi }}^{\left(\mathit{l}\right)}\right),(10)$

where $\mathrm{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 $\mathbf{x}={\left[{\mathbf{\varphi }}^{\mathrm{T}},{\mathrm{k}}_{\mathrm{e}\mathrm{f}\mathrm{f}}\right]}^{\mathrm{T}}$.

Meanwhile, $\mathbf{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

$\mathbf{F}\left(\mathbf{x}\right)\equiv \left[\begin{array}{l}{\mathbf{F}}_{\mathrm{\varphi }}\left(\mathbf{x}\right)\\ {\mathbf{F}}_{\mathit{k}}\left(\mathbf{x}\right)\end{array}\right]=\left[\begin{array}{l}\mathbf{\varphi }-{\left(\mathbf{I}-\mathbf{D}{\mathbf{L}}^{-1}\mathbf{M}\mathbf{S}\right)}^{-1}\frac{1}{{\mathit{k}}_{\mathit{e}\mathit{f}\mathit{f}}}\mathbf{D}{\mathbf{L}}^{-1}\mathbf{M}\mathbf{\chi }\mathbf{f}\mathbf{\varphi }\\ {\mathit{k}}_{\mathit{e}\mathit{f}\mathit{f}}-\widehat{\mathrm{c}}||\mathbf{f}{\left(\mathbf{I}-\mathbf{D}{\mathbf{L}}^{-1}\mathbf{M}\mathbf{S}\right)}^{-1}\frac{1}{{\mathit{k}}_{\mathit{e}\mathit{f}\mathit{f}}}\mathbf{D}{\mathbf{L}}^{-1}\mathbf{M}\mathbf{\chi }\mathbf{f}\mathbf{\varphi }{||}_1\end{array}\right],(11)$

where ${\mathbf{F}}_{\mathrm{\varphi }}\left(\mathbf{x}\right)$ and ${\mathbf{F}}_{\mathrm{k}}\left(\mathbf{x}\right)$ denotes the residual functions for $\mathrm{\varphi }$ and ${\mathrm{k}}_{\mathrm{e}\mathrm{f}\mathrm{f}}$, respectively. However, as mentioned above, the solution of the ${\left(\mathbf{I}-\mathbf{D}{\mathbf{L}}^{-1}\mathbf{M}\mathbf{S}\right)}^{-1}$ operation needs nested inner iterations with a certain amount of transport sweeping steps (multiple calculation of ${\mathbf{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

${\mathbf{F}}_{\mathrm{\varphi }}\left(\mathbf{x}\right)=\mathbf{\varphi }-\mathbf{D}{\mathbf{L}}^{-1}\mathbf{M}\mathbf{S}\mathbf{\varphi }-\frac{1}{{\mathit{k}}_{\mathit{e}\mathit{f}\mathit{f}}}\mathbf{D}{\mathbf{L}}^{-1}\mathbf{M}\mathbf{\chi }\mathbf{f}\mathbf{\varphi }=\mathbf{\varphi }-\mathbf{D}{\mathbf{L}}^{-1}\mathbf{M}\left(\mathbf{S}+\frac{1}{{\mathit{k}}_{\mathit{e}\mathit{f}\mathit{f}}}\mathbf{\chi }\mathbf{f}\right)\mathbf{\varphi },(12)$

${\mathbf{F}}_{\mathrm{k}}\left(\mathbf{x}\right)={\mathit{k}}_{\mathit{e}\mathit{f}\mathit{f}}-\widehat{\mathrm{c}}||\mathbf{f}\mathbf{D}{\mathbf{L}}^{-1}\mathbf{M}\left(\mathbf{S}+\frac{1}{{\mathit{k}}_{\mathit{e}\mathit{f}\mathit{f}}}\mathbf{\chi }\mathbf{f}\right)\mathbf{\varphi }{||}_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.

2.3 Parallel code design for evaluation of residual functions

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 $\mathbf{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 $\mathbf{x}$ to $\mathrm{\varphi }$ and ${\mathrm{k}}_{\mathrm{e}\mathrm{f}\mathrm{f}}$ in each MPI rank. 2) Calculate the new scalar fluxes $\stackrel{\sim }{\mathrm{\varphi }}$ by using the KBA algorithm where $\stackrel{\sim }{\mathrm{\varphi }}=\mathbf{D}{\mathbf{L}}^{-1}\mathbf{M}\left(\mathbf{S}+1/{\mathrm{k}}_{\mathrm{e}\mathrm{f}\mathrm{f}}\mathbf{\chi }\mathbf{f}\right)\mathrm{\varphi }$, and then get the new k-eigenvalue as ${\stackrel{\sim }{\mathrm{k}}}_{\mathrm{e}\mathrm{f}\mathrm{f}}=\widehat{\mathrm{c}}||\mathbf{f}{\stackrel{\sim }{\mathbf{\varphi }}||}_1$. 3) Calculate the residual functions $\mathbf{F}\left(\mathbf{x}\right)$ according to Eqs 12, 13, where $\mathbf{F}\left(\mathbf{x}\right)=\left[\begin{array}{l}\mathbf{\varphi }-\stackrel{\sim }{\mathbf{\varphi }}\\ {\mathrm{k}}_{\mathit{e}\mathit{f}\mathit{f}}-{\stackrel{\sim }{\mathrm{k}}}_{\mathrm{e}\mathrm{f}\mathrm{f}}\end{array}\right]$.

2.4 Solution strategy

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 ${\mathbf{x}}^{\left(0\right)}$ 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

$\frac{{‖\mathbf{F}\left({\mathbf{x}}^{\left(\mathit{n}\right)}\right)‖}_2}{{‖\mathbf{F}\left({\mathbf{x}}^{\left(0\right)}\right)‖}_2}<{10}^{-6},(14)$

where ${‖\bullet ‖}_2$ denotes the L-2 norm and $\mathbf{F}\left({\mathbf{x}}^{\left(0\right)}\right)$ denotes the residuals of the iterative initial value. In each Newton step ($\mathit{n}$), the adaptive convergence criterion with the Eisenstat-Walker forcing term ${\eta }^{\left(\mathit{n}\right)}$ (Eisenstat and Walker, 1996) is

$‖\mathbf{J}\mathbf{a}\left({\mathbf{x}}^{\left(\mathit{n}\right)}\right)\delta {\mathbf{x}}^{\left(\mathit{n}\right)}+\mathbf{F}\left({\mathbf{x}}^{\left(\mathit{n}\right)}\right){\Vert }_2<{\eta }^{\left(\mathit{n}\right)}‖\mathbf{F}\left({\mathbf{x}}^{\left(\mathit{n}\right)}\right){\Vert }_2,(15)$

${\eta }^{\left(\mathit{n}\right)}=0.1{\left(\frac{\mathbf{F}\left({\mathbf{x}}^{\left(\mathit{n}\right)}\right)}{\mathbf{F}\left({\mathbf{x}}^{\left(\mathit{n}-1\right)}\right)}\right)}^{2.0}.(16)$

FIGURE 1

FIGURE 1. Iterative strategies and flowchart for the comePSn_JFNK and comePSn methods.

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.

3 Numerical results and analysis

In this section, we apply the comePSn_JFNK and comePSn codes to two series of complicated multi-dimensional pin-by-pin neutron transport models. First, the popular 2-D/3-D KAIST-3A benchmark problems with Gd burnable absorbers are simulated to analyze the numerical accuracy, the parallel performance and the computational efficiency of the comePSn_JFNK and comePSn codes. Then, the more complicated 2-D/3-D full-core MOX/UOX pin-by-pin models with different insertion forms of the control rods are further employed to test the high numerical efficiency of the comePSn_JFNK code compared with the comePSn code in an order of 1000 CPU cores on the computer clusters with Hygon C86 7185 32-core processors. In all of our numerical simulations, the S₈ Carlson quadrature set (Longoni, 2004) is adopted to discretize the angular direction space.

3.1 KAIST-3A benchmark problem

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.

FIGURE 2

FIGURE 2. Radial geometry of the KAIST-3A benchmark problem.

FIGURE 3

FIGURE 3. Fuel assembly configuration of the KAIST-3A benchmark problem: (A) UOX fuel assembly and (B) MOX fuel assembly.

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.

TABLE 1

TABLE 1. Results of k-eigenvalues on different mesh sizes per pin for 2-D KAIST-3A benchmark problems.

FIGURE 4

FIGURE 4. Relative error (%) of power density between comePSn_JFNK and comePSn codes on different mesh sizes per pin for 2-D KAIST-3A benchmark problems: (A) 1 × 1 for ARO, (B) 2 × 2 for ARO, (C) 4 × 4 for ARO, (D) 32 × 32 for ARO, (E) 1 × 1 for ARI, (F) 2 × 2 for ARI, (G) 4 × 4 for ARI and (H) 32 × 32 for ARI.

FIGURE 5

FIGURE 5. Relative error (%) of pin-wise power density on different mesh sizes per pin compared with those on 32 × 32 meshes per pin for 2-D KAIST-3A benchmark problems using the comePSn_JFNK code: (A) 1 × 1 for ARO, (B) 2 × 2 for ARO, (C) 4 × 4 for ARO, (D) 1 × 1 for ARI, (E) 2 × 2 for ARI and (F) 4 × 4 for ARI.

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.

TABLE 2

TABLE 2. Comparison of the computational cost for 2-D KAIST-3A benchmark problems.

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.

FIGURE 6

FIGURE 6. Relative error (%) of radial average power density between comePSn_JFNK and comePSn codes for 3-D KAIST-3A benchmark problems: (A) ARO and (B) ARI.

FIGURE 7

FIGURE 7. Radial average power density and 3-D power density distribution for 3-D KAIST-3A benchmark problems using the comePSn_JFNK code: (A) radial power density for ARO, (B) 3-D power density for ARO, (C) radial power density for ARI and (D) 3-D power density 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.

TABLE 3

TABLE 3. Comparison of the computational cost for 3-D KAIST-3A benchmark problems (100 CPU cores).

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.

FIGURE 8

FIGURE 8. Parallel efficiency for 3-D KAIST-3A ARO case of the comePSn_JFNK code.

3.2 Full-core MOX/UOX pin-by-pin models

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.

FIGURE 9

FIGURE 9. Fuel assembly configuration of the full-core MOX/UOX pin-by-pin models: (A) UOX fuel assembly and (B) MOX fuel assembly.

FIGURE 10

FIGURE 10. Radial core configuration of the full-core MOX/UOX pin-by-pin models: (A) model A, (B) model B, (C) model C and (D) model D.

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%.

TABLE 4

TABLE 4. Results of k-eigenvalues for 2-D full-core MOX/UOX pin-by-pin models A–D.

FIGURE 11

FIGURE 11. Power density distribution using the comePSn_JFNK code for 2-D full-core MOX/UOX pin-by-pin models: (A) model A, (B) model B, (C) model C and (D) model D.

FIGURE 12

FIGURE 12. Relative errors (%) of radial average power density between comePSn_JFNK and comePSn codes for 2-D full-core MOX/UOX pin-by-pin models: (A) model A, (B) model B, (C) model C and (D) model D.

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.

TABLE 5

TABLE 5. Comparison of the computational cost for 2-D Full-core MOX/UOX pin-by-pin models A–D.

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.

FIGURE 13

FIGURE 13. Radial average power and 3-D power density distribution for 3-D Full-core MOX/UOX pin-by-pin models: (A) radial power density for model A, (B) 3-D power density for model A, (C) radial power density for model B, (D) 3-D power density for model B, (E) radial power density for model C, (F) 3-D power density for model C, (G) radial power density for model D and (H) 3-D power density for model D.

TABLE 6

TABLE 6. Comparison of the k-eigenvalues and the computational cost for 3-D full-core MOX/UOX pin-by-pin models (1024 CPU cores).

4 Conclusion

In this paper, a parallel Jacobian-free Newton Krylov discrete ordinates method (comePSn_JFNK) is developed to solve the multi-dimensional multi-group pin-by-pin neutron transport problems, which combines the parallel JFNK framework and the parallel Sn method based on KBA algorithm. The comePSn_JFNK code exhibits the good efficiency compared with the traditional parallel Sn code with power iterative strategy (comePSn). Furthermore, by making full use of transport sweeping iterative framework, the corresponding residual functions of the parallel JFNK framework can be easily evaluated by only a “flattened” power iterative process which includes a single outer iteration without nested inner iterations. The comePSn_JFNK and the comePSn codes are developed in the unified COupled Multiphysics Environment (COME).

By simulating the 2-D/3-D KAIST-3A pin-by-pin benchmark problems and the 2-D/3-D full-core pin-by-pin MOX/UOX models, the speedups of the comePSn_JFNK code are 6–15 for the KAIST-3A benchmark problems and 15–30 for the full-core MOX/UOX pin-by-pin models, which indicates the comePSn_JFNK code has significant efficiency advantage for complicated 3-D full-core pin-by-pin models. Further study is to improve the parallel efficiency on the massively parallel computers, to apply the parallel coarse mesh finite difference methods and the parallel diffusion synthetic acceleration method as preconditioners, and to extend the comePSn_JFNK method to solve the pin-by-pin multi-physics coupled problems.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author contributions

YZ: Software, modeling, validation, data analysis and original draft preparation. XZ: Conception, methodology, software, review, editing and supervision.

Funding

This research was funded by National Natural Science Foundation of China: 12005073; National Key Research and Development Program of China: 2018YFE0180900; National Key Research and Development Program of China: 2020YFB1901600; Project of Nuclear Power Technology Innovation Center of Science Technology and Industry for National Defense: HDLCXZX-2021-HD-033.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Alcouffe, R. E., Baker, R. S., Dahl, J. A., Davis, E. J., Saller, T. G., and Turner, S. A. (2005). Partisn: A time-dependent parallel neutral Particle transport code system. LAUR-05-3925 (Los Alamos, NM, USA: Los Alamos National Laboratory).

Google Scholar

Bailey, T. S., and Falgout, R. D. “Analysis of massively parallel discrete-ordinates transport sweep algorithms with collisions,” in Proceedings of the International Conference on Mathematics, Computational Methods & Reactor Physics, Saratoga Springs, NY, USA, May 2009.

Google Scholar

Baker, R. S. (2017). An SN algorithm for modern architectures. Nucl. Sci. Eng. 185, 107–116. doi:10.13182/nse15-124

CrossRef Full Text | Google Scholar

Baker, R. S., and Koch, K. R. (1998). An Sn algorithm for the massively parallel CM-200 computer. Nucl. Sci. Eng. 128, 312–320. doi:10.13182/nse98-1

CrossRef Full Text | Google Scholar

Carlson, B. G. (1953). Solution of the transport equation by Sn approximations. LA-1599 (Los Alamos, NM, USA: Los Alamos Scientific Laboratory Report).

Google Scholar

Cho, N. Z. (2000). Daejeon, South Korea: Korea Advanced Institute of Science and Technology benchmark report. Available at: https://github.com/nzcho/Nurapt-Archives/tree/master/KAIST-Benchmark-Problems. Benchmark problem 3A: MOX fuel-loaded small PWR coreMOX Fuel Zoning, 7 Group Homog. Cells.

Google Scholar

Eisenstat, S. C., and Walker, H. F. (1996). Choosing the forcing terms in an inexact Newton method. SIAM J. Sci. Comput. 17, 16–32. doi:10.1137/0917003

CrossRef Full Text | Google Scholar

Eisenstat, S. C., and Walker, H. F. (1994). Globally convergent inexact Newton methods. SIAM J. Optim. 4, 393–422. doi:10.1137/0804022

CrossRef Full Text | Google Scholar

Esmaili, H., Kazeminejad, H., Khalafi, H., and Mirvakili, S. (2020). Subchannel analysis of annular fuel assembly using the preconditioned Jacobian-free Newton Krylov methods. Ann. Nucl. Energy 146, 107616. doi:10.1016/j.anucene.2020.107616

CrossRef Full Text | Google Scholar

Evans, T. M., Stafford, A. S., Slaybaugh, R. N., and Clarno, K. T. (2010). Denovo: A new three-dimensional parallel discrete ordinates code in scale. Nucl. Technol. 171, 171–200. doi:10.13182/nt171-171

CrossRef Full Text | Google Scholar

Gaston, D., Newman, C., Hansen, G., and Lebrun-Grandie, D. (2009). Moose: A parallel computational framework for coupled systems of nonlinear equations. Nucl. Eng. Des. 239, 1768–1778. doi:10.1016/j.nucengdes.2009.05.021

CrossRef Full Text | Google Scholar

Gill, D. F. (2009). “Newton-Krylov methods for the solution of the k-eigenvalue problem in multigroup neutronics calculations,”. PhD thesis (Pennsylvania, PA, USA: Pennsylvania State University).

Google Scholar

Hajizadeh, A., Kazeminejad, H., and Talebi, S. (2018). Formulation of a fully implicit numerical scheme for simulation of two-phase flow in a vertical channel using the Drift-Flux Model. Prog. Nucl. Energy 103, 91–105. doi:10.1016/j.pnucene.2017.11.009

CrossRef Full Text | Google Scholar

Hossain, M. A., and Alam, J. M. (2012). Assessment of a symmetry-preserving JFNK method for atmospheric convection. Comput. Phys. Commun. 269, 108113. doi:10.1016/j.cpc.2021.108113

CrossRef Full Text | Google Scholar

Knoll, D. A., and Keyes, D. E. (2004). Jacobian-free Newton–Krylov methods: A survey of approaches and applications. J. Comput. Phys. 193, 357–397. doi:10.1016/j.jcp.2003.08.010

CrossRef Full Text | Google Scholar

Knoll, D. A., Park, H., and Newman, C. (2011). Acceleration of k-eigenvalue/criticality calculations using the jacobian-free Newton-Krylov method. Nucl. Sci. Eng. 167 (2), 133–140. doi:10.13182/nse09-89

CrossRef Full Text | Google Scholar

Longoni, G. (2004). Advanced quadrature sets and acceleration and preconditioning techniques for the discrete ordinates method in parallel computing environments. PhD Thesis (Gainesville, FL, USA: University of Florida).

Google Scholar

Pawlowski, R., Bartlett, R., Belcourt, N., Cannon, S. R., and Warren, H. R. (2011). A theory manual for multi-physics code coupling in LIME, version 1.0. SAND2011-2195 (Albuquerque, NM, USA: SAND Report).

Google Scholar

Saad, Y., and Schultz, M. H. (1986). Gmres: A generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM J. Sci. Stat. Comput. 7, 856–869. doi:10.1137/0907058

CrossRef Full Text | Google Scholar

Smith, M. A., Lewis, E. E., and Na, B. C. (2005). Benchmark on deterministic transport calculations without spatial homogenization: MOX fuel assembly 3-D extension case. Paris, France: Nuclear Energy Agency Organisation for Economic Co-Operation and Development, NEA/NSC/DOC.

Google Scholar

Turner, J. A., Clarno, K., Sieger, M., Bartlett, R., Collins, B., Pawlowski, R., et al. (2016). The virtual environment for reactor applications (VERA): Design and architecture. J. Comput. Phys. 326, 544–568. doi:10.1016/j.jcp.2016.09.003

CrossRef Full Text | Google Scholar

Walker, E. D., Collins, B., and Gehin, J. C. (2019). Low-order multiphysics coupling techniques for nuclear reactor applications. Ann. Nucl. Energy 132, 327–338. doi:10.1016/j.anucene.2019.04.022

CrossRef Full Text | Google Scholar

Xu, L. F., Cao, L. Z., Zheng, Y. Q., and Wu, H. (2018). Hybrid MPI-communication for the multi-angular SN parallel sweep on 3-D regular grids. Ann. Nucl. Energy 116, 407–416. doi:10.1016/j.anucene.2018.03.003

CrossRef Full Text | Google Scholar

Yaremchu, M., and Panteleev, G. (2022). On the Jacobian approximation in sea ice models with viscous-plastic rheology. Ocean. Model. (Oxf). 177, 102078. doi:10.1016/j.ocemod.2022.102078

CrossRef Full Text | Google Scholar

Zhang, L., Zhang, B., Liu, C., Zheng, J., Zheng, Y., and Chen, Y. (2017). Calculation of the C5G7 3-D extension benchmark by ARES transport code. Nucl. Eng. Des. 318, 231–238. doi:10.1016/j.nucengdes.2017.04.011

CrossRef Full Text | Google Scholar

Zhou, X. F. (2023). Jacobian-free Newton Krylov coarse mesh finite difference algorithm based on high-order nodal expansion method for three-dimensional nuclear reactor pin-by-pin multiphysics coupled models. Comput. Phys. Commun. 282, 108509. doi:10.1016/j.cpc.2022.108509

CrossRef Full Text | Google Scholar

Zhou, X. F. (2022). Jacobian-free Newton Krylov two-node coarse mesh finite difference based on nodal expansion method. Nucl. Eng. Technol. 54, 3059–3072. doi:10.1016/j.net.2022.02.005

CrossRef Full Text | Google Scholar

Zhou, X. F. (2022). Operator split, Picard iteration and JFNK methods based on nonlinear CMFD for transient full core models in the coupling multiphysics environment. Ann. Nucl. Energy. revised.

Google Scholar

Zhou, X. F., Zhong, C. M., and Zhang, Y. Y. (2022). Jacobian-free Newton Krylov two-node coarse mesh finite difference based on nodal expansion method for multiphysics coupled models. Ann. Nucl. Energy 168, 108915. doi:10.1016/j.anucene.2021.108915

CrossRef Full Text | Google Scholar