Weak form of the discrete ordinates method equations

The S_N form was obtained by solving the Boltzmann transport equation along discrete directions and replacing the integrals over the unit sphere with quadrature sets. The steady-state multi-group S_N neutron transport equation can be written as ${\text{Ω}}_m\cdot \nabla {\psi }_m^g\left(r\right)+{\displaystyle {\sum }_{\text{t}}^g\left(r\right)}{\psi }_m^g\left(r\right)=Q_m^g\left(r\right),g=1,\mathrm{...},G;m=1,\mathrm{...},M;r\in D$ (1) Where $D$ is the domain. ${\psi }_m^g\left(r\right)$ is the angular flux in group $g$, direction ${\text{Ω}}_m$ and position $r$. ${\displaystyle \sum {}_t^g}$ is the macroscopic total cross-section in group $g$. $Q_m^g\left(r\right)$ is a source term in group $g$, direction ${\text{Ω}}_m$ and position $r$. For fixed source problems, it can be expressed as: $Q_m^g\left(r\right)={\displaystyle \sum _{g\text{'}=1}^G{\displaystyle \sum _{n=0}^N\frac{2n+1}{4\pi }}}{\displaystyle \sum {}_{s,n}^{\text{g"}\to \text{g}}}\left(r\right){\displaystyle \sum _{k=-n}^n{\varphi }_n^k\left(r\right)}{Y}_n^k\left({\text{Ω}}_m\right)+q_g\left(r\right)$ (2) and for eigenvalue problems, it can be expressed as: $\begin{array}{c}{Q}_m^g\left(r\right)={\displaystyle \sum _{g\text{'}=1}^G{\displaystyle \sum _{n=0}^N\frac{2n+1}{4\pi }}}{\displaystyle \sum {}_{s,n}^{\text{g"}\to \text{g}}}\left(r\right){\displaystyle \sum _{k=-n}^n{\varphi }_n^k\left(r\right)Y_n^k\left({\text{Ω}}_m\right)}\\ +\frac{{\chi }_g}{4\pi k_{\text{eff}}}{\displaystyle \sum _{g\text{'}=1}^G{\upsilon }^{g\text{'}}\left(r\right){\displaystyle \sum {}_f^{\text{g"}}}\left(r\right){\varphi }_{g\text{'}}\left(r\right)},\end{array}$ (3) ${\varphi }_g\left(r\right)={\displaystyle \sum _{m=1}^M{\omega }_m{\psi }_m^g\left(r\right)},$ (4) where $G$ and $M$ is total groups and angles. ${\displaystyle {\sum }_{s,n}^{\text{g"}\to \text{g}}\left(r\right)}$ is the $n$th-order Legendre moment of scattering cross section from group $g\text{'}$ to group $g$ at position $r$. ${\displaystyle {\sum }_{\text{f}}^{\text{g"}}\left(r\right)}$ is the macroscopic fission cross section in group $g\text{'}$ at position $r$. ${\upsilon }^{g\text{'}}\left(r\right)$ is the average number of neutrons produced per fission in group $g\text{'}$. $k_{eff}$ is the eigenvalue. ${\varphi }_n^k\left(r\right)$ represents kth, nth-order flux moment at position $r$. $Y_n^k\left({\text{Ω}}_m\right)$ is $k$th, $n$th-order spherical harmonic in direction ${\text{Ω}}_m$ at position $r$. $q_g\left(r\right)$ is fixed-source in group $g$ at position $r$. ${\varphi }_g\left(r\right)$ is scalar flux in group g at position $r$. ${\omega }_m$ is weight factor of quadrature set.

Eq. (1) is a linear hyperbolic equation: It is necessary to transform this "strong form" differential-integral equation into "weak form" for solving. The spatial domain $D$ is decomposed into $K$ elements, which can be unstructured because the FEM is independent of the space dimension and grid choice. Consider $\partial V_k$ as the surface of the $k$th element and $V_k$ as the volume of the $k$th element, where $k=1,\mathrm{...},K$. To obtain the weak form, Eqs. (1) is multiplied by an arbitrary test function ${\psi }_m^{*}\left(r\right)$ and integrated into any element $k$. Eq. (1) can be rewritten in monoenergetic form as $\begin{array}{l}{\displaystyle {\int }_{V_k}{\text{Ω}}_m\cdot \nabla {\psi }_m\left(r\right){\psi }_m^{*}\left(r\right)}\text{d}V\\ +{\displaystyle {\int }_{V_k}{\text{Σ}}_{\text{t}}\left(r\right){\psi }_m\left(r\right)}{\phi }_m^{*}\text{d}V={\displaystyle {\int }_{V_k}{Q}_m\left(r\right){\psi }_m^{*}}\text{d}V.\end{array}$ (5)

The divergence theorem is then applied to Eq. (5), yielding: $\begin{array}{l}{\displaystyle {\int }_{\partial V_k}{\text{Ω}}_m}\cdot n\left(r\right){\psi }_m^{*}\left(r\right){\psi }_m\left(r\right)\text{d}S-{\displaystyle {\int }_{V_k}{\psi }_m\left(r\right){\text{Ω}}_m\cdot \nabla {\psi }_m^{*}\left(r\right)}\text{d}V+\\ {\displaystyle {\int }_{V_k}{\displaystyle {\Sigma }_{\text{t}}\left(r\right)}{\psi }_m\left(r\right)}{\psi }_m^{*}\left(r\right)\text{d}V={\displaystyle {\int }_{V_k}{Q}_m\left(r\right){\psi }_m^{*}\left(r\right)}\text{d}V\end{array}$ (6) where $n(r)$ is the unit normal vector to the surface of cell $\partial V_k$ at position $r$. Inner products on the volume and surface of any element are introduced to express the bilinear form [23]. ${\langle a,b\rangle }_{V_k}={\displaystyle {\int }_{V_k}a\left(r\right)b\left(r\right)}\text{d}V$ (7) ${\left(a,b\right)}_{\partial V_k}={\displaystyle {\int }_{\partial V_k}\left|{\text{Ω}}_m\cdot n\left(r\right)\right|a\left(r\right)b\left(r\right)}\text{d}S$ (8)

Substituting Eq. (7) and (8) into (6), $\begin{array}{l}{\left({\psi }_m^{*},{\psi }_m\right)}_{\partial V_k}-{\langle {\psi }_m,{\text{Ω}}_m\nabla {\psi }_m^{*}\rangle }_{V_k}\\ +{\langle {\displaystyle {\sum }_t{\psi }_m},{\psi }_m^{*}\rangle }_{V_k}={\langle Q_m,{\psi }_m^{*}\rangle }_{V_k}.\end{array}$ (9)

The general boundary can be divided into three types based on its relative position: interior boundary, outflow domain boundary $\left(\partial V_k^{+}=\left\{r\in \partial V_k,\text{Ω}\cdot n\left(r\right)>0\right\}\right)$, and inflow domain boundary $\left(\partial V_k^{-}=\left\{r\in \partial V_k,\text{Ω}\cdot n\left(r\right)<0\right\}\right)$. The inflow-domain boundary conditions were prescribed as either vacuum or reflective for the incoming directions. So, Eq. (9) can be rewritten as $\begin{array}{l}{\left({\psi }_m^{*},{\psi }_m^{+}\right)}_{\partial V_k^{+}}-{\langle {\psi }_m,{\text{Ω}}_m\nabla {\psi }_m^{*}\rangle }_{V_k}+{\langle {\text{Σ}}_t{\psi }_m,{\psi }_m^{*}\rangle }_{V_k}\\ ={\langle Q_m,{\psi }_m^{*}\rangle }_{V_k}+{\left({\psi }_m^{*},{\psi }_m^{-}\right)}_{\partial V_k^{-}}+{\left({\psi }_m^{*},{\psi }_m^{\text{in}}\right)}_{\partial V_k^{\text{in}}},\end{array}$ (10) where ${\psi }_m^{+}$, ${\psi }_m^{-}$, and ${\psi }_m^{\text{in}}$ are the outflow angular, inflow angular, and interior face fluxes, respectively. Eq. (10) represents the weak form of the S_N equation.

The discontinuous Galerkin finite element discretization scheme

The Galerkin method is theoretically based on the weighted residual method, which multiplies the original partial differential equation by a test function. Seeking a numerical approximation of the angular flux term using Lagrange polynomials of degree P ${\psi }_m={\displaystyle \sum _{j=1}^{J_P}{\phi }_{m,j}{b}_j}\left(r\right),$ (11) where $b_j\left(r\right)$ is the Lagrange shape function, and ${\phi }_{m,j}$ is an unsolved coefficient. The matrix form is convenient for solving Eq. (10), which expresses ${\phi }_{m,j}$ and $b_j\left(r\right)$ in row vector notation, and Eq. (11) can be written as ${\psi }_m=\left[{\phi }_{m,1},\mathrm{...},{\phi }_{m,J_P}\right]{\left[b_1,\mathrm{...},b_{J_P}\right]}^{{\rm T}}={\phi }_m{b}^{{\rm T}}.$ (12)

In Galerkin discretization schemes, test functions share similar function space as basis functions, ${\psi }_m^{*}={\left[b_1,\mathrm{...},b_{J_P}\right]}^{{\rm T}}=b.$ (13)

Substituting Eq. (12) and (13) into (10), $\begin{array}{r}{\left(b,b^{{\rm T}}\right)}_{\partial V_k^{+}}{\phi }_m^{+}-{\langle b^{{\rm T}},{\text{Ω}}_m\nabla b\rangle }_{V_k}{\phi }_m+{\langle {\text{Σ}}_t{b}^{{\rm T}},b\rangle }_{V_k}{\phi }_m\\ ={\langle Q_m,{\psi }_m^{*}\rangle }_{V_k}+{\left({\psi }_m^{*},{\psi }_m^{-}\right)}_{\partial V_k^{-}}+{\left({\psi }_m^{*},{\psi }_m^{in}\right)}_{\partial V_k^{\text{in}}}.\end{array}$ (14)

The upwind scheme is the essence of the DGFEM. Notably, although the test functions on the edges are always taken from within the element, the primal functions on the edges are taken upwind with respect to the streaming direction. For outgoing edges, the primal functions are taken from within the element itself, whereas for incoming edges, the primal functions are taken from the upwind neighbor element.

Multithreaded parallel upwind sweep algorithm

2.3.1

Cell-based sweep algorithm

In the S_N method, a transport sweep was performed for all elements along each discrete direction. The ordering in which the elements and degrees of freedom (DOF) traverse affects the speed of convergence. As mentioned in Section 1, the DGFEM makes it possible to use a cell-based sweeping procedure in S_N methods rather than assembling a global matrix.

The element sweep order for each discrete direction in the ThorSNIPE code was developed using an upwind scheme based on the deal.II library. In our case, an array with all relevant cells is created and sorted in the upwind direction using a "dealii :: DoFRenumbering :: CompareDownstream" object in the deal.II FEM library. A straightforward 2D model with three material regions and various sweep schemes is shown in Fig. 1. Cells in a conventional scheme are often categorized from edge to center to make it easier to facilitate the creation of a global matrix. The downwind system performs renumbering in the direction opposite to the discrete direction, whereas the upwind scheme manages renumbering alongside it. The definitions of "upwind" and "downwind" are determined by multiplying the unit normal vector to the cell's surface by the cosines of the corresponding discrete direction. In reality, the mesh location, discrete direction, and sweep-define the order of each element. The sweep solution technique also complies with neutron transport laws. Not all pathways are mutually exclusive, and each distinct direction results in a unique order in which this solution technique can be applied.

Fig. 1Full-screen View

(Color online) (a) Material region distribution; (b) Conventional scheme; (c) Upwind scheme; (d) Downwind scheme

From the perspective of solving the equations, various sweep schemes also affect the calculation of the incoming angular fluxes ${\psi }_m^{\text{in}}$ in Eq. (14). Primal functions obtained from neighboring elements require numerical values with a certain precision. The calculation error was gradually reduced when the transport equation was solved along the direction of particle motion using the upwind scheme. However, other schemes must process nondirectional data, adapt, and learn based on the data received. This process can be laborious and error-prone, and the accumulation of errors inevitably affects the robustness of the numerical calculation method. Thus, the entire solution procedure required more iterations.

Algorithm 1 displays the cell-based sweep algorithm for each group and angle once the sweep order has been determined. The local equation can be divided into two parts as follows: The transport matrix first makes contributions from the streaming and removal operators determined by the fixed mass matrix and total cross-section. Second, the right-hand side (RHS) vector includes contributions from scattering source operations and fixed source or fission source operators. The in-group and cross-group couplings comprised the contributions of the scattering sources. Lastly, a "dealii :: types :: global dof index" object translates the local solution to the global angular solution. This cell-based sweep technique delivers the DGFEM solution to all cells in a single sweep order without coupling integral terms, fully utilizing the key benefits of DGFEM approaches.

Full-screen View

Algorithm 1 cell-based sweep for each group and angle
Input: g, m, qfixed(r)
for upwind sorted cell do
	for DOF of cell i do
		for DOF of cell $j$ do
			$transport\_matrix\left(i,j\right)=-{\displaystyle {\int }_{V_k}{b}_i^{{\rm T}}\left(r\right){\text{Ω}}_m\cdot \nabla b_j\left(r\right)}dV+{\displaystyle {\int }_{V_k}{\displaystyle {\sum }_t\left(r\right)}{b}_i\left(r\right)}{b}_j^{{\rm T}}\left(r\right)dV$
	for moment n do
		for moment k do
			for DOF of cell i do
				$source\_rhs(i)+={\displaystyle {\int }_{V_k}{b}_i\left(r\right)}\left({\displaystyle \sum {}_r}\left(r\right){\varphi }_n^k\left(r\right)+q_{fixed}\left(r\right)\right)Y_n^k\left({\text{Ω}}_m\right)dV$
	for faces on cell do
		if outflow face then
				for DOF of face i do
					for DOF of face j do
						$boundary\_matrix(i,j)={\displaystyle {\int }_{\partial V_k^{+}}{\text{Ω}}_m}\cdot n\left(r\right)b_i\left(r\right)b_j^{{\rm T}}\left(r\right)dS$
		else if inflow face then
				for DOF of face i do
					$boundary\_rhs\left(i\right)=\{\begin{array}{c}{\displaystyle {\int }_{\partial V_k^{-}}{\text{Ω}}_m}\cdot n\left(r\right)b_i\left(r\right){\psi }_m^{-}\left(r\right)dS\\ {\displaystyle {\int }_{\partial V_k^{\text{in}}}{\text{Ω}}_m}\cdot n\left(r\right)b_i\left(r\right){\psi }_m^{in}\left(r\right)dS\end{array}$
	cell_rhs(i)=source_rhs(i)+boundary_rhs(i)
	cell_matrix(i,j)=transport_matrix(i,j)+boundary_matrix(i,j)
	solution(i)=cell_matrix(i,j)-1cell_rhs(i)
	for DOF of cell i do
		angular_solution[m](local_dof_indices[i])=solution(i)

Show more

The following two steps were employed to compute the local contributions of an individual cell and RHS: First, the integral is transformed from an actual cell $k$ into a unit/reference cell $\widehat{k}$. For example, the streaming operator is transformed into $\begin{array}{l}{\displaystyle {\int }_{V_k}{b}^{{\rm T}}\left(r\right){\text{Ω}}_m\nabla b\left(r\right)}\text{d}V\\ ={\displaystyle {\int }_{V_{\widehat{k}}}\left[J^{-1}\left(\widehat{r}\right)b^{{\rm T}}\left(\widehat{r}\right)\right]\left[{\text{Ω}}_m{J}^{-1}\left(\widehat{r}\right)\widehat{\nabla }b\left(\widehat{r}\right)\right]}\left|\det J\left(\widehat{r}\right)\right|\text{d}\widehat{V},\end{array}$ (15) where the hat indicates the reference coordinates, and $J\left(\widehat{r}\right)$ is the Jacobian of the mapping $r=F_k\left(\widehat{r}\right)$. Second, this integral is approximated using the Gauss-Legendre quadrature. This yields the formula $\begin{array}{l}{\displaystyle {\int }_{V_k}{b}^{{\rm T}}\left(r\right){\text{Ω}}_m\nabla b\left(r\right)}\text{d}V\\ ={\displaystyle {\int }_{V_{\widehat{k}}}\left[J^{-1}\left({\widehat{r}}_q\right)b^{{\rm T}}\left({\widehat{r}}_q\right)\right]\left[{\text{Ω}}_m{J}^{-1}\left({\widehat{r}}_q\right)\widehat{\nabla }b\left({\widehat{r}}_q\right)\right]}\left|\det J\left({\widehat{r}}_q\right)\right|w_q,\end{array}$ (16) where q is the index of the quadrature point, ${\widehat{r}}_q$ is the location of the reference cell, and $w_q$ is the Gaussian quadrature weight. In this study, triangular and tetrahedral grids were used for the 2D and 3D cases, respectively.

2.3.2

Multithread parallel strategy

Parallel transport computation has gained popularity in the nuclear field because of the emergence of high-performance computing clusters [24]. How the processors schedule parallel jobs significantly affects the performance of code using parallel computing. The two basic techniques are spatial-domain decomposition (SDD) [25] and angular-domain decomposition (ADD) [26]. SDD techniques subdivide the spatial domain into subdomains, such that the common mesh-sweep method is performed concurrently in numerous subdomains at diverse angles. However, parallelizing the S_N spatial sweep is more challenging because of the upwind-downwind dependencies between the domain boundaries.

ADD techniques are frequently considered more expandable and require less memory [27]. They entail a mesh sweep along a similar quadrant and are the same for all angles in the chosen quadrature [28]. All operations are statically scheduled to the participating processors without negatively affecting the load balance and are ideally suited for personal computers (PCs). In shared memory machines, the traditional approach to parallelism is to break code down into threads. The deal.II leverages the Threading Building Blocks (TBB) library [29] as basic wrappers to build an object called "tasks.” TBB offers transferable interfaces across many different platforms and abstracts the specifics of run-state signals onto individual threads.

The entire procedure for the multithreaded parallel approach in the ThorSNIPE code is shown in Fig. 2. Both inner and outer iterations contain an angular loop. The number of CPU cores determines the total number of tasks in a PC system. The angular loops were divided into a specific number of subranges that were evenly distributed across the threads. A cell-based sweep sequence is used for each angular loop. A subrange is a task. The tasks are essentially independent during the iterations. They performed the calculations using the same mesh and equation structure as the respective quadrature set. However, the storage of angular fluxes is based on the local DOF indices calculated using sweep schemes. The calculation results for different subranges were collected to update the scattering source. Following the execution of the subranges, threads keep themselves busy by stealing complete portions of the subranges from other threads if they have finished their work. Thus, the code can fully exploit the system thread resource, which requires frequent interruptions by the operating system to allow other threads to execute on the available processor cores.

Fig. 2Full-screen View

Flowchart depicting multithreaded parallel strategy

IAEA benchmark

The IAEA benchmark is a 2D five-region calculation problem consisting of two large source zones and two large absorber zones surrounded by light water [30] and is designed for IAEA advanced reactor neutron transport computation. There is a strong nonuniformity in the fuel region, as shown in Fig. 3. The cross-sections are given in Table 1, and all boundary conditions are vacuum.

Fig. 3Full-screen View

(Color online) Geometry of the IAEA 5-region fixed problem

This benchmark problem is mainly used to analyze the effect of the element sweep order on the cell-based algorithm efficiency and assert the superiority of the upwind sweep scheme. Under the P_NT_N-S₈ angular approximation, comparisons of the numerical results for the fixed source and eigenvalue cases simulated with 2338 elements are presented in Tables 2 and 3, respectively. The relative deviations in the average flux were within 2.50% in the fixed-source case and within 5.00% in the eigenvalue case. All the numerical results show excellent agreement with the references, demonstrating the accuracy and validity of the algorithm and modeling approach for such strongly heterogeneous problems. The change in the element sweep order does not make a large difference to the results. This is because the transport sweep operations were performed individually on each element, and the results were added. Because matrix addition is commutative and these operators in the equation can be described as applying a sum of matrices, the orders of various components will not impact the final result.

As shown in Fig. 4, the number of iterations was compared with the change in the number of elements in the fixed-source calculation. The numerical results indicated that the upwind sweep scheme had a stable convergence speed, regardless of the number of elements. The relationship between the number of iterations and elements grows exponentially in both types of sweep schemes. The primary functions on the edges are determined based on the upwind edges with regard to the streaming direction, and transport problems involve a directional flow of information. The upwind scheme fully uses the initial boundary conditions and provides spatial flux moments for the neighboring elements. The other sweep schemes increased the number of source and power iterations required to meet the convergence criteria in the subsequent matrix calculations to compensate for the lack of constraints. In addition, no overhead cost of assembling and solving the stiffness matrix in the current framework was observed for the three sweep schemes mentioned above. Therefore, the upwind sweep algorithm can save significant computational costs with fewer iterations and is suitable for DGFEM simulations.

Fig. 4Full-screen View

Number of iterations with the change in the number of elements for different sweep schemes

Kobayashi-3i benchmark

The Kobayashi-3i benchmark [31] is a bend duct-type problem with void regions in a highly absorbing medium, as shown in Fig. 5. The cross sections used in this benchmark problem are listed in Table 4. This problem was set up to analyze the features of spatial and angular discretization and demonstrate the parallel efficiency of the multithreaded strategy.

Fig. 5Full-screen View

(Color online) Geometry of the Kobayashi-3i benchmark

Angular flux distribution over full domain with mesh sizes of 5 cm, 4 cm, 3 cm, 2.5 cm, and 2 cm are modeled and simulated. The quality of the quadrature sets was quantified by taking the root mean square (RMS) of the relative scalar flux of all N point detectors. $\text{RMS}=\sqrt{\frac{1}{N}{\displaystyle \sum _{i=1}^N{\left(\frac{{\varphi }_{cal}-{\varphi }_{ref}}{{\varphi }_{ref}}\right)}^2}}$ (17) where the reference solutions ${\varphi }_{\text{ref}}$ were obtained using an analytical method [31]. The RMS of the relative errors in the scalar flux using different quadrature sets versus the number of ordinates is presented in Fig. 6a. It should be noted that the various approximations in the DGFEM with different mesh sizes converged at variable speeds and started to flatten as the number of ordinates increased under S₂₀ angular approximations.

Fig. 6Full-screen View

(a) Root mean square of Scalar flux L2-error versus the order of quadrature sets for the Kobayashi-3i benchmark; (b) Speedup ratio versus the order of quadrature sets for the Kobayashi-3i benchmark

The ThorSNIPE code runs on a PC with 8 cores / 16 threads (one 11th Gen Intel(R) Core(TM) i7-11700F CPU running at 2.5 GHz) to demonstrate the performance and gain of the multithread technology. The speedup ratio with respect to serial mode versus the order of quadrature sets is presented in Fig. 6b, where the parallel implementation efficiency exceeds 100% when the S12 angular approximations are attained, and the mesh size is less than 3 cm. The comparative results show that the proposed method can achieve a better acceleration effect as the number of meshes increases. It can also be observed that the speedup ratio increases more slowly when the S20 angular approximations are attained. In this case, each thread has over 30 parallel tasks, and memory access becomes a bottleneck in the matrix solvers for partial differential equations. The numerical results are shown in Fig. 7, and the neutron fluxes at the key points are listed in Table 5. Compared with MCNP [31], the ThorSNIPE code effectively reduces angular discretization errors and fulfills the transport computational requirement to a certain degree.

Fig. 7Full-screen View

(Color online) Mesh distribution and scalar flux of the Kobayashi-3i benchmark

VENUS-3 benchmark

The VENUS-3 benchmark was proposed to validate the capabilities of the calculation methodologies and cross-section libraries for predicting fluence rates in RPVs because it has a very clean structural geometry representing standard PWR pressure vessel conditions [32]. Fig. 8 illustrates the geometric configuration and material distribution of the VENUS-3 facility. The VENUS-3 model includes two types of fuel rods: the PLSA, core baffle, water reflector, barrel, Pyrex control rods, barrel, neutron pad, three types of grids, bottom support, reactor vessel, and jacket inner walls. Detailed descriptions and qualifications were obtained from the benchmark document. The benchmark arranges 386 dosimeters, with 244 ⁵⁸Ni(n, p)⁵⁸Co dosimeters, 104 ¹¹⁵In(n, n’)^115mIn dosimeters, and 38 ²⁷Al(n, a)²⁴Na dosimeters placed in 268 different spatial locations. The experimental rate sets are provided as a ratio of the calculated reaction rate to the average dosimeter cross-section, named the equivalent fission flux.

Fig. 8Full-screen View

(Color online) Geometry of the Venus-3 benchmark

The numerical results shown in Fig. 9 were modeled and simulated using the P_NT_N-S₈ angular discretization, 199 KASHIL-E70 neutron groups multigroup cross-sections [33], and 407819 tetrahedral elements. Fig. 10 shows a C/E comparison of the equivalent fission fluxes at the three detector positions. For the ¹¹⁵In(n, n’)^115mIn dosimeters, almost all the equivalent fission flux deviations are included within +/-10%, except for a detector position (No. 31 located in the core barrel with 0.7° angle) where the equivalent fission flux numerical value is affected by reflective boundary. For the ²⁷Al(n, a)²⁴Na dosimeters, 92% of the equivalent fission flux deviations were within +/-10%, and the remaining three dosimeters overestimated the corresponding experimental values by approximately 11.9%. For the ⁵⁸Ni(n, p)⁵⁸Co dosimeters, 95% of the equivalent fission flux deviations are limited within +/-10%, except for 11 detector positions where the equivalent fission flux numerical values overestimated the corresponding experimental values by about 10.1%-20.7%. They are concentrated at the junction of the PLSA region, outer baffle, and core barrel. The scalar flux distribution is not smooth because of the strong spatial nonuniformity of the geometry model, and higher requirements are proposed for local mesh refinement, which significantly degrades the simulated accuracy of these regions. The overall agreement between the numerical results from the ThorSNIPE code and the experimental results was very good, which shows that the ThorSNIPE code is highly reliable and stable for complex shielding calculations.

Fig. 9Full-screen View

(Color online) Mesh distribution and fast neutron fluxes of the VENUS-3

Fig. 10Full-screen View

C/E comparison of equivalent fission fluxes at detector positions: (a) ⁵⁸Ni(n,p)⁵⁸Co from No. 1 to 139; (b) ⁵⁸Ni(n,p)⁵⁸Co from No. 140 to 244; (c) ¹¹⁵In(n,n’)^115mIn; (d) ²⁷Al(n,a)²⁴Na

The numerical results [32] from TORT-3.2 with 1004295 hexahedral grids and 26 BUGLE-96 neutron energy groups were selected as comparative references. For the uniformly distributed 3.3% fuel regions, TORT-3.2 with hexahedral structured meshes achieved better results. The ThorSNIPE code with tetrahedral unstructured meshes can decrease the numerical oscillation for complex regions. For the ²⁷Al(n, a)²⁴Na dosimeters, the numerical results of ThorSNIPE tended to be lower than those of TORT-3.2 in the core barrel. This underestimation might have been caused by an error in the average dosimeter cross section.

In Fig. 11, the speedup ratios of the three quadrature sets with the P₃ Legendre expansion order and the four Legendre expansion orders with P_NT_N-S₈ angular discretization are shown. The test conditions were consistent with TORT-3.2 code. Overall, the behavior of the proposed algorithm was similar to that of the Kobayashi-3i benchmark. For P_NT_N-S₁₂, the acceleration ratio was 6.47 with 407819 elements. Comparing the results of different expansion orders, the P₅ expansion order has the steepest upward trend, and the maximum speedup ratio achieved was 6.98. The numerical results indicate that the proposed algorithm achieves a favorable speedup ratio as the scale of the mesh increases, as expected.

Fig. 11Full-screen View

Speedup ratio with the change in the (a) order of quadrature sets; (b) order of Legendre expansion