Eigenvalue problems of neutron diffusion equation

The neutron diffusion equation is expressed as follows: $\frac{1}{v}\frac{\partial \varphi }{\partial t}=-L\varphi +(1-\beta ){\chi }_{\text{p}}F\varphi +{\displaystyle \sum _{i=1}^6}{\chi }_{\text{d}\mathrm{,}i}{\lambda }_i{C}_i\mathrm{.}$ (1) $\frac{\partial C_i}{\partial t}={\beta }_{i}F\varphi -{\lambda }_i{C}_i\mathrm{,}i=\mathrm{1,2,}\cdots \mathrm{,6}$ (2) where ν represents the neutron velocity, ϕ denotes the neutron flux, t stands for time, β represents the total delayed neutron fraction, χ_p signifies the prompt neutron fission spectrum, χ_d,i represents the delayed neutron fission spectrum for the ith group, λi denotes the decay constant of the ith group precursor nucleus, ci represents the concentration of the ith group precursor nucleus, and βi represents the fraction of delayed neutrons for the ith group. The operators L and F are defined as follows: $\begin{array}{c}L\varphi =-\nabla \cdot D(r\mathrm{,}E)\nabla \varphi (r\mathrm{,}E\mathrm{,}t)+{\text{Σ}}_{\text{t}}\left(r\mathrm{,}E\right)\varphi \left(r,E\mathrm{,}t\right)\\ -{\displaystyle {\int }_0^{\infty }}{\text{Σ}}_{\text{s}}\left(r\mathrm{,}{E}^{\text{'}}\to E\right)\varphi \left(r\mathrm{,}{E}^{\text{'}}\mathrm{,}t\right)\text{d}{E}^{\text{'}}\end{array}$ (3) $F\varphi ={\displaystyle {\int }_0^{\infty }}\nu {\text{Σ}}_{\text{f}}\left(r\mathrm{,}{E}^{\text{'}}\right)\varphi \left(r\mathrm{,}{E}^{\text{'}}\mathrm{,}t\right)\text{d}{E}^{\text{'}}$ (4) By introducing a balance factor of 1/λ and adjusting the right-hand side fission neutron source term to make the system pseudocritical, the time derivative term is eliminated, resulting in a steady-state λ-mode neutron diffusion equation, which can be expressed using the following matrices: $L\varphi =\frac{1}{\lambda }\chi F\varphi$ (5) where $\chi =(1-\beta ){\chi }_{\text{p}}+{\displaystyle {\sum }_6^{i=1}{\beta }_i{\chi }_{\text{d}\mathrm{,}i}}$ is the total (or average) fission neutron spectrum. The neutron diffusion equation, which considers only prompt neutrons and disregards delayed neutrons, is as follows: $\frac{1}{v}\frac{\partial \varphi }{\partial t}=-L\varphi +(1-\beta ){\chi }_{\text{p}}F\varphi \mathrm{.}$ (6) Supposing that the solution to Eq. (6) is ${\varphi }_{\alpha \mathrm{,}\text{p}}\left(r\mathrm{,}E\right)e^{\alpha t}\mathrm{.}$ (7) By substituting Eq. (7) into Eq. (6), the neutron diffusion equation for the prompt α-mode can be obtained using Eq. (8): $\frac{\alpha }{v}{\varphi }_{\alpha }=-L{\varphi }_{\alpha }+(1-\beta ){\chi }_{\text{p}}F{\varphi }_{\alpha }\mathrm{.}$ (8) The neutron diffusion equations of λ mode and prompt α mode can be solved mathematically by the eigenvalue problem in the form of Ax=λx or Ax=λBx, yielding a series of discrete eigenvalues ${\lambda }_1\mathrm{,}{\lambda }_2\mathrm{,}\dots \mathrm{,}{\lambda }_n$, ${\alpha }_1\mathrm{,}{\alpha }_2\mathrm{,}\dots \mathrm{,}{\alpha }_n$ along with the corresponding eigenvectors (modes) ${\varphi }_{\lambda \mathrm{,1}}\mathrm{,}{\varphi }_{\lambda \mathrm{,2}}\mathrm{,}\dots \mathrm{,}{\varphi }_{\lambda \mathrm{,}n}$, ${\varphi }_{\alpha \mathrm{,1}}\mathrm{,}{\varphi }_{\alpha \mathrm{,2}}\mathrm{,}\dots \mathrm{,}{\varphi }_{\alpha \mathrm{,}n}$. The largest eigenvalue is referred to as the fundamental eigenvalue. In λ mode, the fundamental eigenvalue carries the physical meaning of the effective multiplication factor k_eff, and the associated fundamental eigenvector represents the neutron flux distribution in precise or pseudo-critical systems. The λ mode is typically used to characterize the reactivity and neutron flux distribution of reactors near criticality. In contrast, the α eigenvalue can be either a positive or negative real number, indicating that the reactor system is in a supercritical or subcritical state. The α-eigenvalue can also represent the evolution behavior of the unsteady neutron flux of the reactor over time [23]. Additionally, for noncritical systems, α and α modes are physically measurable quantities (typically, the fundamental eigenvalue of the prompt α mode is taken as the prompt neutron attenuation constant).

The corresponding adjoint equations for the above λ and α eigenvalue equations are as follows: $L^{+}{\varphi }_{\lambda }^{+}=\frac{1}{\lambda }{(\chi F)}^{+}{\varphi }_{\lambda }^{+}$ (9) $\frac{\alpha }{v}{\varphi }_a^{+}=-L^{+}{\varphi }_a^{+}+(1-\beta ){\left({\chi }_{\text{p}}F\right)}^{+}{\varphi }_a^{+}$ (10) The λ and prompt α adjoint equations have the same eigenvalues as the forward equations, where the λ adjoint flux also has the physical significance of “neutron importance” and can be used to calculate kinetic parameters such as the effective fraction of delayed neutrons [24, 25].

Methods for solving higher-order modes

In practice, when solving the higher-order modes of the neutron diffusion equation, matrix A in the eigenvalue problem (Ax=λx) is often not directly constructed and solved. Instead, iterative methods such as fixed-source iterations combined with subspace iterations as outer iterative algorithms are employed to obtain higher-order modes. Taking IRAM as the outer algorithm, the calculation principle of the higher-order modes is as follows in Algorithm 1[26, 27]. Here, operator A represents the matrix to be solved or a subfunction equivalent to that matrix.

The mth-order eigenvalues and eigenvectors obtained through the IRAM process are the mth-order modes required for the solution. This is described in Algorithm 1, the multiplication of matrices by subspace vectors (Avj) corresponds to a fixed-source inner iteration. The pseudocodes for computing λ and prompt α modes are shown in Algorithms 2 and 3.

In HARMONY2.0 code, the inner iteration is solved using the finite difference method. In Algorithms 2 and 3, ag, bg, … represent the difference coefficients of the leakage terms after discretization. For example, in the x-direction, the inner node formula is given by Eq. (11). $\begin{array}{l}{a}_g\left(i\mathrm{,}j\mathrm{,}k\right)=b_g\left(i-\mathrm{1,}j\mathrm{,}k\right)\\ \text{ =}-\frac{2D_g\left(i\mathrm{,}j\mathrm{,}k\right)D_g\left(i-\mathrm{1,}j\mathrm{,}k\right)\Delta y(j)\Delta z(k)}{D_g\left(i\mathrm{,}j\mathrm{,}k\right)\Delta x(i-1)+D_g\left(i-\mathrm{1,}j\mathrm{,}k\right)\Delta x(i)}\end{array}$ (11) where ag and bg represent the x^- and x⁺ directions, respectively. At the outer boundary with zero-flux boundary conditions, the corresponding formulas are given by Eq. (12). $a_g^{\text{'}}\left(i\mathrm{,}j\mathrm{,}k\right)=-\frac{2D_g\left(i\mathrm{,}j\mathrm{,}k\right)}{\Delta x(i)}\Delta y(j)\Delta z(k)$ (12) For the y and z directions, the subscripts j, k, … in Eq (11) and (12) are changed accordingly.

The discretized central difference equations of the iteration(in line 5 of Algorithm 2 and line 6 of Algorithm 3) are expressed in Eq. 13. Common numerical methods such as the Jacobi iteration and successive over-relaxation iteration can be used to solve these equations, although specific formulations are not provided here. $\begin{array}{l}{a}_g\left(i\mathrm{,}j\mathrm{,}k\right){\varphi }_g\left(i-\mathrm{1,}j\mathrm{,}k\right)+b_g\left(i\mathrm{,}j\mathrm{,}k\right){\varphi }_g\left(i+\mathrm{1,}j\mathrm{,}k\right)+\\ \dots +e_g\left(i\mathrm{,}j\mathrm{,}k\right){\varphi }_g(i\mathrm{,}j\mathrm{,}k)=S_g\left(i\mathrm{,}j\mathrm{,}k\right)\end{array}$ (13) where $e_g\left(i\mathrm{,}j\mathrm{,}k\right)$ denotes the removal rate as shown in Eq. (14): $\begin{array}{c}{e}_g\left(i\mathrm{,}j\mathrm{,}k\right)={\text{Σ}}_{R\mathrm{,}g}\left(i\mathrm{,}j\mathrm{,}k\right)\Delta x(i)\Delta y(j)\Delta z(k)\\ -a_g\left(i\mathrm{,}j\mathrm{,}k\right)-b_g\left(i\mathrm{,}j\mathrm{,}k\right)-\dots \end{array}$ (14) where $S_g\left(i\mathrm{,}j\mathrm{,}k\right)$ denotes the total number of source terms. The λ and prompt α modes have different expressions, as shown in Eq. (15). The corresponding formulas for FS, SS, PFS, and AS are presented in Algorithms 2 and 3. $S_g\left(i\mathrm{,}j\mathrm{,}k\right)=\{\begin{array}{l}F{S}_g\left(i\mathrm{,}j\mathrm{,}k\right)+S{S}_g\left(i\mathrm{,}j\mathrm{,}k\right)\hfill \\ A{S}_g\left(i\mathrm{,}j\mathrm{,}k\right)-PF{S}_g\left(i\mathrm{,}j\mathrm{,}k\right)-S{S}_g\left(i\mathrm{,}j\mathrm{,}k\right)\hfill \end{array}$ (15)

Homogenization by OpenMC

This subsection introduces the MC homogenization method used to obtain the space- and spectrum-homogenized group constants of the reactor.

Traditional two-step methods for reactor core calculations face challenges in accurately describing complex geometric configurations owing to spatial discretization limitations. Additionally, homogenization processes entail addressing intricate issues such as resonance treatment and leakage corrections. To address these concerns, deterministic methods rely on assumptions and approximations. In contrast, MC homogenization methods utilize tallying functions to calculate neutron reaction rates statistically based on continuous-energy cross-sectional data. The group constants are then obtained using principles such as the conservation of reaction rates. This method not only allows for an accurate description of various geometrical configurations but also avoids challenges such as resonance treatment. Several MC codes have functionalities developed to generate group constants for coupling with downstream deterministic codes, including RMC [28], Serpent [29], OpenMC [30, 31], McCard [32], and MVP [33].

In spatial homogenization using traditional approaches, the reactor core is treated as single assemblies or individual fuel rod lattices. Reflective boundary conditions were assumed during homogenization [34]. Because the traditional assembly or lattice homogenization approach neglects the actual boundary conditions of the assembly or lattice, necessary correction methods or special treatments are inevitable when using these homogenized group constants for core calculations [35]. In contrast, due to its geometric flexibility, MC homogenization can not only implement assembly homogenization (infinite lattice homogenization [36]) with reflective boundary conditions but also achieve global homogenization with "real" boundary conditions [37] and "mesh" homogenization for generating group constants based on a mesh with "real" boundary conditions, among other types of homogenization.

Infinite lattice homogenization using the MC method, similar to deterministic assembly homogenization, involves modeling individual assemblies and then calculating the homogenized group constants using the full reflective boundary conditions of each type of assembly, as shown in Fig. 1a. However, global homogenization establishes a model around the assembly or even the entire core to obtain the homogenized group constants for individual assemblies, as illustrated in Fig. 1b. Mesh homogenization, however, directly divides the core space into meshes, calculating the homogenized constants for each mesh under "real" boundary conditions. Given the flexibility in adjusting the mesh parameters, mesh homogenization offers greater applicability in geometric processing. The difference in the steps of the full-core calculation between the infinite lattice or global homogenization and mesh homogenization is demonstrated in Fig. 2.

Fig. 1Full-screen View

(Color online) Different boundary conditions. (a) Reflective boundary conditions; (b) “real" boundary conditions

Fig. 2Full-screen View

(Color online) Core calculation with different assembly homogenized group constants. (a) Infinite lattice homogenization; (b) mesh homogenization

By utilizing mesh homogenization, the critical calculation, tallying, and cross-sectional data treatment required for subsequent computations can be efficiently achieved in a single run of the MC calculation, without the need to perform homogenization calculations for each assembly separately. Furthermore, the group constant data obtained through mesh homogenization can be stored in accordance with the mesh index, which facilitates coupling with downstream deterministic core calculation programs.

Parallel using OpenMP

In this subsection, we discuss the inner iteration acceleration process by using event-driven parallelization and domain decomposition based on OpenMP [38]. The differences between the two types of parallelization were also discussed.

In the context of higher-order mode calculations or subsequent research, improving the computational speed is very important. In the process from MC homogenization to higher-order mode core calculation, several steps can be parallelized for acceleration, such as critical calculation, inner flux iteration, and the last step for eigenpairs solution. Taking the example of solving the prompt α modes of a PHWR model with a mesh division of 18×18×10 and utilizing HARMONY1.0 for serial computation under two-group conditions with a subspace number of 50, the inner iteration process requires approximately 760 s of CPU time. With 200 subspace vectors and an order of 100, it took more than 4000 s. Such rates are unacceptable for studies requiring higher-order mode expansion and reconstruction analyses requiring large number of higher-order modes. (All the calculations in this study were performed on a Linux system equipped with an Intel i5-13600K CPU and 16GB of RAM.)

2.4.1

Event-driven parallelization

One approach to parallelization that is easily achieved using OpenMP is to parallelize individual code blocks (events). For instance, when writing serial code, it is necessary to nest multilayer FOR loop code blocks to calculate the scattering source terms during the inner iteration. However, in a serial execution, each layer loop must be traversed sequentially. In reality, the computation of terms(such as the scatter source) in individual meshes has no inherent dependency, and can be performed simultaneously. This presents a natural opportunity for parallelization.

Similarly, parallel acceleration can be achieved for events with nested loops such as total source calculations and flux iterations. The events and "communication" nodes of the parallel computation achieved using this principle-based approach (hereafter referred to as event-driven parallelization) are depicted as Fig. 3a.

Fig. 3Full-screen View

Two types of parallelization. (a) Event-driven; (b) domain decomposition

The event-driven approach requires a certain amount of time for thread task communication and waiting between multiple parallel-computation steps. Thread waiting primarily occurs because of disparities in workload and performance variations among threads, resulting in a load imbalance [39]. The time consumed by thread-task communication is primarily attributed to the distribution of computational tasks among threads and the aggregation of computational results from different threads. As the number of threads increases, the time required for thread communication inevitably increases and the load imbalance becomes more severe. Both factors contribute to the reduced parallel efficiency. In certain cases of event-driven parallelization, this can result in slower performance compared to serial execution.

2.4.2

Domain decomposition parallelization

Event-driven parallelization divides parallel regions based on each computation event because of the existence of dependencies between different computation events. For instance, the total-source calculation must wait for the completion of the scatter source calculation before proceeding, necessitating a thread waiting between these two computational events. However, although dependencies exist between different computational events, they exist only in the same mesh. For example, the total source calculation at the mesh point (1, 1, 1) must wait for its scatter source calculation to be completed, but the total source at the mesh point (3, 3, 3) does not depend on the completion of the scatter source at (1, 1, 1). This gives rise to the "domain decomposition" parallel method [6, 40].

The principle and schematic diagram of domain decomposition parallel computation is as follows: Consider all spatial meshes (or energy groups) as the "all region" that needs iterative computation. Based on this, divide the "all region" into N subdomains, where N threads are responsible for the iterative computation in each subdomain(as Fig. 4). This approach ensures that thread task communication occurs only at the beginning and after each thread completes its task[41], thereby significantly reducing its frequency. The events and "communication" nodes of the domains decomposition are as Fig. 3b.

Fig. 4Full-screen View

(Color online) Schematic of domain decomposition

The HARMONY2.0 code

Based on the above principles, we developed an MC-deterministic hybrid two-step method for the higher-order mode calculation code HARMONY 2.0, building upon HARMONY1.0. The outer iteration employs the IRAM algorithm, whereas the inner iteration module is recoded and uses Python modules to couple with the mgxs API of OpenMC. The program modules and overall computational process are outlined in Fig. 5.

Fig. 5Full-screen View

(Color online) Schematic diagram of HARMONY 2.0 code system structure and the overall calculation flow

The complete computational process is as follows: Begin by constructing the reactor model using OpenMC and defining the parameters for mesh homogenization, such as mesh structure, energy groups, and cross-section types. Subsequently, execute the critical calculation to generate homogenized cross-sectional data. Afterward, use the Python coupling modules to connect with Fortran modules, which include calculations for difference coefficients, result processing, and the inner iteration module. Finally, compile the Fortran modules for computing higher-order modes. Notably, Python modules require the reading of a simple text file containing information about the geometry mesh/energy bins. Furthermore, a few modules can be invoked to generate HARMONY2.0’s inner iteration modules for parallel computation or to choose the solution mode (λ or prompt α mode; forward, adjoint calculation). Table 1 lists the names and functions of the Python modules.

Comparison of serial and parallel computing

The α mode of the PHWR-3D-2G problem was chosen for computational testing to compare the results and acceleration effects of serial and parallel computations. The number of subspace vectors was set as 200 to compute the first 100 higher-order modes. The PHWR-3D model (19×19×10) shown in Fig. 6 and the first 10 order’s prompt α values are listed in Table 2. A time comparison between the event-driven domain decomposition and serial computations is presented in Table 3. In the case of parallel computation using domain decomposition, the subdomain partitions were as follows: two threads were allocated with 9×18×10, four threads with 9×9×10, eight threads with 9×9×5, 16 threads with 9×9×2 or 3 (z-direction).

Fig. 6Full-screen View

(Color online)PHWR reactor. (a) x-y direction; (b) x-z direction

The results indicate that serial and parallel computations for both methods yielded the same values. The acceleration ratio of event-driven parallelization reached a maximum of 1.95 with five threads but decreased beyond five threads. In contrast, domain decomposition parallelization outperforms event-driven parallelization significantly, with a maximum acceleration ratio of 6.16, which was achieved using 16 threads.

Computational experiments were conducted using a relaxation factor of 1.0. However, adjusting this factor could potentially further accelerate the computation rates. For instance, in the case of domain decomposition with 16 threads, changing the relaxation factor to 1.05 reduced the computational time from 310.07 to 298.62 seconds. In addition, the HARMONY code was rewritten and ported to the Linux platform. This reduced the time for serial computations from over 4000 seconds to 1908.79 seconds, compared with HARMONY1.0. Therefore, with the parallelization implemented in HARMONY2.0, the computational speed increases appreciably.

Validation of C5G7 model

As shown in Table 2, the first four-order results of PHWR agree well with the reference values. However, because the validation of PHWR was obtained directly through deterministic calculations, it is still necessary to validate the hybrid two-step method used in this study before further discussion. Given the limited number of related benchmarks, we selected the C5G7 model for validation.

In the homogenization calculations, the criticality calculation parameters were set as follows: 20,000 generations with the first 500 inactive and 100,000 particles per cycle for C5G7-2D and 500,000 particles per cycle for C5G7-3D. During the homogenization and mode calculations, the mesh divisions for C5G7-2D is 51×51, the mesh division for C5G7-3D is 51×51×10. The models used in OpenMC were 1/4 or 1/8 symmetric, therefore reflective boundary conditions were employed in the mode calculations. The k_eff values of OpenMC were 1.18640(± 5.6 pcm) and 1.18378(± 0.8 pcm), respectively. The C5G7 model and simplified calculation process are shown in Fig. 7.

Fig. 7Full-screen View

C5G7 model and calculation flow

The calculation results are presented in Table 4, where the reference values were obtained from Gupta et al. and calculated using the SSI-based MC method [17]. Here, we observed a difference between the HARMNOY2.0 computed values and reference values. The difference was 267 pcm for the fundamental eigenvalue of C5G7-2D and 253 pcm for C5G7-3D. As the order increases, the differences gradually increase, which could be due to several reasons such as differences in algorithms and discrepancies between diffusion and transport theories. One reason for this is that the settings of the number of meshes lead to differences in the eigenvalue solutions, particularly for higher orders. This will be discussed in the following subsections.

Verification of complex geometry computation

To validate the applicability of HARMONY2.0, to complex geometric reactor cores, we selected the Hoogenboom problem [43], Advanced Test Reactor (ATR)[22, 46], and MUSE-4 experimental facility[47] as test cases. This subsection discusses the computational results for the three models. Models of the Hoogenboom and ATR are shown in Fig. 8, and the MUSE-4 model is illustrated in Figs. 13.

Fig. 8Full-screen View

(Color online) Hoogenboom (a) and ATR (b) model

3.3.1

Hoogenboom

The Hoogenboom problem, proposed by Hoogenboom et al., is an MC performance benchmark comprising 241 fuel assemblies with the same enrichment level. Each assembly consisted of fuel rods arranged in a 17×17 mesh, and included 25 guide tubes. The axial zone of the core is divided into cold and hot regions [43]. In this study, Hoogenboom problem is chosen as a typical square lattice core model for solution.

The computational results for the different spatial mesh configurations and energy group divisions are presented in Table 5. Radially, the model was divided into one or four meshes per assembly. The energy boundaries of the 2-group energy structure were (0, 0.625, 1.96403×10⁷) in eV, and for 6-group, they were (0, 4.53999×10², 9.11822×10³, 6.73795×10⁴, 4.97871×10⁵, 2.23130×10⁶, and 1.96403×10⁷). For the 8-group configuration, the CASMO 8-group [44] energy structure was employed for the eight-group configuration. The MC critical calculation parameters were set to 1,000,000 particles per cycle, totaling 1150 cycles, with the first 150 cycles being inactive. The results showed that the differences between the fundamental eigenvalues and k_eff obtained from OpenMC were within 200 pcm. Specifically, the fundamental eigenvalue calculation for the 34×34×8 8-group configuration is 1.00126, differing by only 1 pcm from OpenMC.s. Furthermore, we obtained the fission matrix using RMC3.5 and solved for higher-order modes using SciPy [45]. Table 5 presents the results along with those from [15].

Table 5Full-screen View

λ mode results of Hoogenboom

Order	OpenMC	HARMONY2.0					RMC&SciPy	Results of [15].
		19×19×8&2 g	19×19×8&6 g	19×19×8&8 g	34×34×8&8 g	36×36×8&8 g	34×34×8	42×42×20
1	1.00127	1.00293	1.00093	1.00031	1.00126	1.00059	1.00193	0.99919
2	±8.5 pcm	0.99402	0.99051	0.98961	0.99078	0.98991	0.98687	0.98483
3		0.99323	0.98861	0.98791	0.98808	0.98840	0.98659	0.98362
4		0.99263	0.98817	0.98766	0.98844	0.98816	0.97444	0.98469
5		0.98437	0.97792	0.97717	0.97763	0.97768	0.96766	0.96956
6		0.98353	0.97773	0.97679	0.97701	0.97729	0.96664	0.96950
7		0.98183	0.97482	0.97377	0.97474	0.97406	0.96011	0.96693
8		0.98032	0.97420	0.97324	0.97250	0.97302	0.95975	0.96591
9		0.97797	0.97138	0.97097	0.97093	0.97241	0.95954	0.96043
10		0.97553	0.96834	0.96706	0.96593	0.96740	0.94412	0.95671

Show more

As mentioned in the previous subsection, the division of meshes affects the calculation results for higher-order modes (similar to the fission matrix method [15]). The discrepancies in the results for the higher-order modes were partly due to differences in the diffusion calculations under different meshes. In addition, from a mathematical perspective, the calculation of higher-order modes is equivalent to an eigenvalue problem. Different mesh divisions affect the order of the matrix A, thereby altering the total number of corresponding eigenvalues. This also resulted in larger discrepancies in the higher-order eigenvalue results as the order increased across different mesh divisions. As shown in Table 5, the comparison indicates that the lack of an energy group dimension leads to progressively larger discrepancies among the higher eigenvalues. By the 10th order, the fission matrix method results from RMC3.5&SciPy are 0.94412, whereas those from HARMONY 2.0, are all above 0.96.

Figure 9 illustrates the mode distribution of Hoogenboom with a 34×34×8-mesh. The eigenvectors of different groups of fundamental modes exhibited similar shapes. Therefore, only the distributions of the eighth group at z=4 for the first ten orders are displayed. As depicted in the Fig. 9, the first, seventh, and tenth orders exhibit fully symmetric distributions, whereas the second, third, fourth, fifth, and sixth orders exhibit a 1/2 symmetric distribution. The eighth and ninth orders displayed a 1/4 symmetric distribution.

Fig. 9Full-screen View

(Color online) λ mode of Hoogenboom (g=8, z=4, order 1-10)

3.3.2

ATR

The ATR is one of the few large research reactors worldwide and is characterized by high flux, numerous irradiation channels, and reactivity control through multiple drums[22, 46]. It utilizes special curved-plate fuel elements, which makes it challenging to calculate higher-order modes using a deterministic method. Hence, we chose ATR as the complex geometry core model for HARMONY2.0 test.

Table 6 presents ATR’s first ten eigenvalues calculation results with a 50×50×1 mesh, equivalent to a mesh size of 2 cm×2 cm×125 cm. The group energy structure follows CASMO 8-group. The OpenMC critical calculation parameters were 1,000,000 particles for 1150 cycles and using the first 150 cycles were inactive. The computed k_eff for ATR was 1.02293. A comparison revealed an error of only 226 pcm for ATR. Additionally, no approximation corrections were applied during MC homogenization, and further adjustments to parameters such as the energy group structure were not made. Figure 10 illustrates the differences in the fundamental flux distributions among the 8th groups in ATR with a mesh division of 50×50×1.

Table 6Full-screen View

λ mode results of ATR

Nev	ATR
	OpenMC	HARMONY2.0
1	1.02293	1.02519
2	±8.3pcm	0.85252
3		0.86857
4		0.77553
5		0.62682
6		0.54278
7		0.53784
8		0.54650
9		0.46770
10		0.45797

Show more

Fig. 10Full-screen View

(Color online) Fundamental λ mode of ATR (g=1-8)

The fundamental neutron flux distribution of the ATR conformed to the characteristic neutron flux distribution of a thermal reactor. This distribution entails fast neutrons generated from fission and predominantly concentrated in the fuel region, whereas the distribution of epithermal neutrons tends to be relatively flat owing to diffusion. By contrast, thermal neutrons exhibit a lower flux in the fuel region because of absorption, which reflects the absorption effect of the fuel material. Notably, this distribution also highlights the neutron absorption by structural materials, as evidenced by the lower neutron flux in regions containing absorbers, such as the surrounding drums.

The 2nd to 9th orders of ATR’s first and 8th group modes are shown in Fig. 11, exhibiting symmetrical distributions similar to 1/2 and 1/4 across the various orders.

Fig. 11Full-screen View

(Color online) λ modes of ATR (group 1&8, order 2-9)

The corresponding adjoint mode distributions for the 1st and 8th groups are shown in Fig. 12. A comparison revealed that the adjoint fundamental mode distribution of the 8th group exhibits a pattern “opposite” to that of the forward fundamental distribution in the fuel region. In other words, the adjoint neutron flux is larger in the fuel region, reflecting the physical significance of the “neutron importance" associated with the adjoint neutron flux.

Fig. 12Full-screen View

(Color online) Adjoint λ modes of ATR (group 1&8, order 1-8)

3.3.3

MUSE-4

Compared with the λ mode, the calculation results of the higher-order modes in the prompt α mode are limited because of fewer experiments compared with critical benchmarks and lower computational efficiency. Therefore, we chose the MUSE-4 experimental facility for validation, which provides experimental data for prompt neutron attenuation constants [47].

The MUSE-4 subcritical experimental apparatus is a modification of the zero-power MASURCA facility that is expanded to accommodate a wider range of experiments. At its core, it utilizes MOX fuels in the active region, with a central aluminum tube extending towards the accelerator. This configuration serves as an external neutron source for simulating the target of an Accelerator-Driven System (ADS) reactor. Surrounding this central feature, the lead components replicate the spallation target of the ADS, coupled with a reflector layer composed of a sodium-stainless steel blend. The outermost layer consists of stainless-steel shielding, which provides radiation protection. This setup allows control of the reactor’s subcriticality by adjusting the loading of the fuel assemblies. When loaded with 976 fuel assemblies, its k_eff corresponds to 0.97 (Fig. 13).

Fig. 13Full-screen View

(Color online) MUSE-4 model in OpenMC (with mesh) (a)x-y direction; (b) x-z direction

The parameters for the OpenMC homogenization calculation were as follows: 500,000 particles, 100 inactive cycles, and 1000 active cycles. The computed k_eff is 0.96817. The parameters for the mesh homogenization are defined with a 6-group energy structure with energy boundaries of (0, 4.53999×10², 9.11822×10³, 6.73795×10⁴, 4.97871×10⁵, 2.23130×10⁶, and 1.96403×10⁷) in eV. The mesh spacing in the xy direction was 5.3 cm, whereas in the z-direction, it was divided into 10.3/12 cm. The numbers of mesh points were 28×30×14 and 28×30×12 for the two configurations, respectively. The computed eigenvalues for both the λ and prompt α modes for the first ten orders are presented in Table 7

Table 7Full-screen View

MUSE-4 eigenvalue results (first 10)

OpenMC	Experiment value/(s-1) [47]	HARMONY2.0
		28×30×14		28×30×12
		λ	α/s-1	λ	α/s-1
keff=0.96817± 3.2pcm	α=-46308 (keff=0.97)	0.96963	-43430	0.96717	-45561
		0.69617	-72951	0.69353	-73006
		0.60302	-73376	0.60063	-73452
		0.60091	-74452	0.59986	-74516
		0.46344	-75652	0.46018	-75803
		0.46288	-75964	0.46095	-76078
		0.44332	-77311	0.44054	-77471
		0.39024	-77433	0.38549	-77650
		0.34971	-77657	0.35699	-77784
		0.34700	-81859	0.34327	-81871

Show more

The results indicate that the computed fundamental eigenvalues for both λ modes under the two different mesh configurations differ from the OpenMC calculation by less than 200 pcm. Similarly, the fundamental prompt α eigenvalue closely aligns with the experimental values. Figure 14 shows the computed first-order eigenvectors of λ and the corresponding α eigenvectors for a mesh configuration of 28×30×12. The shapes of the fundamental distributions for the same energy group are similar for the λ and prompt α modes, differing only in the numerical values.

Fig. 14Full-screen View

(Color online) Fundamental mode of MUSE-4 (group 1–6, z=7). (a)λ mode; (b) prompt α mode

The computed results for the 2–9 orders of the first group’s λ modes and the corresponding prompt α modes are shown in Fig. 15. In the λ-mode, the distributions for each order adhere approximately to symmetric distribution patterns, such as 1/2 and 1/4. Although the prompt α-mode still approximately satisfies the symmetric distribution patterns, its symmetry is less pronounced. In addition to the similarity in the fundamental distribution shapes, there are significant differences in the shapes of the higher-order modes between the λ and prompt α modes.

Fig. 15Full-screen View

(Color online) Higher-order modes of MUSE-4 (group 1, order 2-9, z=7). (a) λ mode; (b) prompt α mode

Computation time evaluation

Based on the cases discussed in the previous section, the execution times for each part of HARMONY2.0 are shown in Fig. 16. In these cases, event-driven parallel examples utilized six threads, whereas domain decomposition employed 16 threads.

Fig. 16Full-screen View

(Color online) Execution time of different cases (first 10 orders)

The most time-consuming part of these cases was the MC homogenization calculation(all cases in Fig. 16 uses 1,000,000 particles with 12 threads). In our tests, some cases using 100,000 particles yielded the same fundamental eigenvalues as those using 1,000,000 particles. For example, the Hoogenboom problem, with a mesh division of 34×34×8 and 8 groups, required only 2710.89 seconds of MC calculation with 100,000 particles (1,000,000 required 25837.21 seconds). However, the use of fewer particles could lead to asymmetrical flux results because of the characteristics of MC calculations. Additionally, because MC homogenization is based on tally functions, if the particle numbers are not sufficiently large, the cross-sections of some meshes might exhibit increased errors owing to insufficient particle counts or even result in a zero cross-section value. These issues may be resolved using variance-reduction techniques.

In addition to the MC calculations, the most time-consuming part is the inner iteration when calculating the first 10 modes with 20 subspace vectors. Taking the C5G7-3D tests as examples, the ratio of the time spent on inner iterations to that spent on other parts(mainly the DNEUPD function [26] and its associated calculations) was approximately 11.85:1. Our tests showed that at the 500th mode with 600 subspace vectors, the inner iteration takes approximately 1446.98 s, while the time for other parts increases to 701.54 s, bringing the ratio to 2.06:1.

Fig. 16 also shows that domain decomposition achieved considerable acceleration in all test cases. Furthermore, the calculations in α mode are significantly slower than those in λ mode. In the case of the MUSE-4 model, despite using domain decomposition, the inner iterations for the first 10 modes required 1518.40 seconds (approximately 60.35 times that of the λ mode), and the time for the other parts increased to 51.87 seconds. Thus, it is evident that acceleration in the α mode is necessary.