Abstract
For the verification test of some scientific calculation programs, different comparison methods are commonly applied to ensure the correctness of the computations. However, it is difficult to verify whether the testing output is correct, because the oracles which include the expected output are not always available or too hard to get. For this reason, the authors focus on using the Richardson Extrapolation to estimate the convergences of the numerical solution on different levels of mesh refinement. These numerical convergence properties can be applied to verification test, without the need for giving the oracles. In the present study, the authors take the program test of the multigroup neutron diffusion equations as a study case and propose the Richardson Extrapolation-based verification method. Three verification criterions are obtained based on our approach. In addition, a test experiment is conducted demonstrating the validity of our theoretical results.
1. Introduction
The importance of scientific software in mission-critical and safety-critical applications has increased dramatically during the last three decades. The engineers, developers, and users of scientific software face a critical problem: How should confidence in computational science and engineering be critically assessed [1]? Since any subtle defects in critical system software will cause catastrophic consequences, the credibility of the computational results must be raised to a higher level.
Verification test is the major process for improving this confidence of the software quality. Great progress has been made in the aspects of verification in scientific software [2]. There remains a great challenge in this verification community; that is, there is no reliable test oracle to indicate what the correct output should be for arbitrary input, which was called oracle problem. This problem commonly arises when conducting verification testing of scientific software; many scientific applications fall into the category of nontestable programs [3] where an oracle is unavailable or too difficult to implement. However, many testers of scientific software generally adopt the following mechanisms as test oracles: (a) comparing with analytical solutions, simulation results, tabulated values, or hand-calculations [4]; (b) verifying with standard mathematical libraries or reference software packages [5]. Additionally, other methods for solving the oracle problem have been developed, including formal proof [6], analytic solutions for simplified physics, method of manufactured solutions [7], PDE benchmark solutions, and metamorphic testing [8, 9]. For example, around 30 verification benchmarks have been developed by the National Agency for Finite Element Methods and Standards (NAFEMS), which have been used widely in verification. In addition, ANSYS [10] and ABAQUS have roughly 270 formal verification tests. Readers are referred to [11–13] for more details on verification. The commonality of these methods is to use the analytical solution, simulation results, and reference values as test oracle. However these methods are unavailable for every newly investigated program.
In order to solve the oracle problem in scientific applications, the mathematical models and numerical methods are factors that should be investigated again. The a posteriori estimation methods based on Richardson Extrapolation are examined [15]. Emphasis is placed on discretization error estimation methods based on Richardson Extrapolation [16]. In the process of using Richardson Extrapolation to analyze the order of error accuracy of the numerical solution, the exact solution of the partial differential equation is actually an auxiliary role and does not need to be given. In other words, Richard Extrapolation methods can be applied to verify the no-test oracle problem. Moreover, the relevant verification criterions could be constructed by the analysis process based on Richardson Extrapolation method [16].
Since the quality of the nuclear reactor design software is closely related to the safety and economy of reactors, the verification test is very important for improving the credibility of programs. Actually, the program for solving multigroup diffusion equations is the primary procedure of the nuclear reactor design software [14]. In the critical state of the reactor, the multigroup diffusion equations are usually described as the solution of Lambda Modes problem or λ-eigenvalue problem [17], which provides the fundamental eigenvalue (the maximal eigenvalue) that is called the effective multiplication factor (k-effective). In this case, the neutron flux is a corresponding eigenfunction. Computations of the Lambda Modes problem provide the results of the k-effective and neutron flux distribution for reactor design calculations. Hence, if the program is calculated incorrectly, it will lead to error prediction and affect the safety and economy of the reactor. Though traditional comparison methods are commonly applied in verification test of partial differential equations program, it is difficult to give expected values to verify whether the calculation results are correct.
In the present study, the program of the multigroup neutron diffusion equations is performed as a study case. We investigate the mathematical models and the finite difference numerical algorithm of the multigroup diffusion equations. Based on Richardson Extrapolation, the convergence properties of numerical solution in and are obtained, and then the rigorous verification criterions are constructed for verification. However, we would not need to give exact solutions to verification test.
The rest of the paper is organized as follows. In Section 2 we introduce the mathematical model. Section 3 presents the main numerical method and convergence analysis process based Richardson Extrapolation. Test experiments are conducted, and the results are reported in Section 4. Finally, we give a summary and conclude the paper in Section 5.
2. Background
2.1. Mathematical Model
The steady state multigroup neutron diffusion equations are the basic mathematical model in reactor design, where the governing equations are given by (1) and (2).where is the energy group, is the diffusion coefficient of the energy group , is the removal cross section of the energy group , is the scattering cross section from energy group to , is the average number of neutrons emitted by fission of the energy group , is the fission cross section of the energy group , is the neutron scalar flux of the energy group , is the integrated fission spectrum of the energy group , and is the effective multiplication factor of reactor. is the fission term. In the multigroup formulation, the neutron diffusion equations are represented by a coupled system of partial differential equations for the flux.
Assuming a bounded 2D or 3D domain with a convex boundary , the domain is composed of the finite number of subassemblies , . For with an outer boundary conditions, we havewhere is the outer normal to the boundary ; at an interface between dissimilar regions with continuous conditions, we haveandwhere and represent the value of both sides of the interface, respectively. The quantities , , , , are constant in each .
To simplify the manipulations, (1) and (2) are expressed in a matrix system by the following form:where the vector of the neutron flux is , the coefficient matrix is defined withand the coefficient matrix is defined with .
2.2. Power Iteration
The power iteration (PI) method is the most used method to obtain the fundamental eigenvalue and its eigenvector [18]. This iterative source method is described in detail by [19]. The power iteration algorithm process is as follows.(I)Starting with a positive initial estimate for and , they are substituted into the right hand of (6) to obtain the neutron flux . The upper labels represent the number of iterations.(II)For the i-th iteration, (6) is written as where .(III)Updating the value for new and the new flux vector , and are calculated by the following equations, respectively, with the stop criterion as follows: where and are small constants.(IV)If the result does not correspond to the stop criterion, return to step (II).
The iteration process of step (II) is called inner iteration, which usually applies the finite differences, finite elements, and other methods to solve (1). Readers are referred to [17, 20, 21] for more details on the mathematical theory. In particular, the finite difference method is a widely used method.
At each step of the power iteration procedure, if we omit the superscript and subscript, (1) is expressed as G uncoupled self-consistency elliptic boundary value problems. Equation (1) is rewritten by
Here, contains the fission and scattering terms. Actually, the fission source is not known so that the inner iteration must be embedded in step (III), which is called outer iteration. Therefore, we fill the values of the inner iteration to the outer iteration.
3. Scientific Software and Numerical Method
Many large and complex reactor physics calculation software packages include programs for solving the multigroup neutron diffusion equations, which widely adopted the differential method to carry on spatial discretization [22]. In the current study, we design and compile a calculation program called HARMONY-2, which also adopts the finite difference method to discrete the multigroup neutron diffusion equation, and invokes ARPACK which is an open source package developed by the Department of Computing and Applied Mathematics of RICE University. The ARPACK implements the IRAM (Implicitly Restarted Arnoldi Method) algorithm, which is presented in [23]. HARMONY-2 is used to calculate high order harmonics of the neutron diffusion equations.
3.1. The Finite Difference Method
We consider an example of one-dimensional and multigroup neutron diffusion equations:
where , , in is zero, and is continuous inside of .
Assuming is divided into a limited number of subregions. Space discretization in one-dimensional X geometry is shown in Figure 1.
Let be the split size and be the split number in the X geometry. The medium in each subdivision area is uniform, so the various cross sections and diffusion coefficients are constant. Using a backward difference format, we obtain the equations as formulation (14). Equation (14) is written as an overrelaxation iterative format: where superscript is the inner iteration count and is the outer iteration count of calculation algorithm, where
Boundary conditions are as follows:Calculation flow of the program algorithm for (13) is shown in Figure 2.
3.2. Richardson Extrapolation
Based on the Richardson Extrapolation method, the convergence of the numerical solution in (1) and (2) was investigated. Let be the numerical solution of (13) in a suitably fine mesh and be the exact solution of this problem. The Richardson Extrapolation procedure assumes that the numerical scheme is first-order accurate or second-order accurate. Considering the Sobolev space , we assume , and we have
Specifically, the numerical solution is calculated by the same discretization scheme with three different grid sizes: , , and , respectively.
If , the error can be rewritten byEquations (24) to (26) yieldThen, we getIf , we have and that the mesh is refined or coarsened by a factor of two. Similarly, we have the following.
Based on the above derivation, the numerical solution convergence of and was obtained by Richardson Extrapolation method. During the above analysis process, it is very clear that the exact solutions of neutron diffusion equations are eliminated in the derivation, so it is not needed to provide the exact solution as the test oracle. Moreover, it is very convenient to verify the program by checking whether the numerical solution satisfies the convergence properties. Moreover, as the mesh is encrypted, the neutron flux and Max-eigenvalues (effective multipliers) converge, respectively, which is also a property that can be used for verification. Based on the convergence properties, three verification criterions are constructed. The first and second one are the numerical solution convergence of and , respectively. The last one is the Max-eigenvalues convergence with the mesh refinement.
4. Test Experiments
4.1. SLAB_1D_2G Benchmark Problem
The purpose of the test experiment is to apply the three verification criterions presented above to verify the program HARMNOY-2. In this experiment, we select the SLAB_1D_2G benchmark problem as the computational case, which is the one-dimensional slab reactor and double-group diffusion problem [14]. The thickness of the slab is 450cm with uniform mesh. The reactor core geometry is shown in Figure 3, and the cross section parameters are shown in Table 1.
(1)The HARMONY-2 program performs the SLAB_1D_2G benchmark calculation. The grid is divided into 18 equal parts in one-dimensional Cartesian coordinate system. For comparison with [14], the output results of the eigenvalues are set as four modes eigenvalues.(2)Verification testing by three verification criterions: Let neutron flux be , where is the neutron scalar flux for 1st group energy, is the neutron flux for 2nd group energy, and denotes the size of grid. We assume that , , , , , and are presented as different refined grids, respectively. The point in X axis is selected as test point. HARMONY-2 is executed multiple times depending on the different meshes. The calculated results are postprocessed, and the value of the flux on observed point is taken out and normalized, where the and denote the numerical solution convergence of and , respectively, as follows:
4.2. Test Results and Discussions
The calculation results of the four modes eigenvalues are shown in Table 2 and compared with the results of [14]. Table 2 indicates that the maximum absolute error between the results of our program and lecture [14] is very small. In addition, the neutron flux of , in 18 equal elements mesh are shown in Figures 5 and 6, respectively, which indicate that the distribution of the neutron scalar flux is geometrically symmetric. The three verification criterions are applied to verification test for the numerical solution of the multigroup neutron diffusion equations. Table 3 shows that the convergence of on the neutron flux value of point in different meshes. In Table 3, the convergence order approximates to 4. The derivative values of neutron flux on point in different refined grids are listed in Table 4. As the mesh is encrypted, the results are about 2.0 order accuracy in Table 4. The calculation results of Max-eigenvalue in different refined grids are shown in Table 5. Figure 4 shows that Max-eigenvalue converges with the mesh encryption. When the mesh is divided into 72 equal parts, the converges to around 9.86943E-01. Those results indicate that calculation of HARMONY-2 satisfies the three verification criterions constructed by the proposed method in this paper.
Table 2 actually reflects a comparison method, which is commonly used in verification tests. Although this method is easy to perform by comparing with the reference values, it is difficult to get reference values for some newly investigated programs. Our method not only is more rigorous than traditional comparison methods, but also does not need reference solutions, exact solutions, or reference values etc. Though our method is based on convergence properties of numerical solution, it is not difficult to get these convergence properties from the other second-order elliptic partial differential equations by our method. Clearly, our method can be applied to verify the programs of these equations. Therefore, our method is very helpful in alleviating oracle problem and improving the quality of scientific software.
5. Conclusion
In this paper, we proposed a Richardson Extrapolation-based verification method, which used the Richardson Extrapolation on different levels of mesh refinement to estimate the convergences of the numerical solution, and we constructed the verification criterions based on these properties. Compared with traditional methods, our method did not need the exact solution, analysis solutions, and reference solutions as expected values. The calculation program of multigroup neutron diffusion equation was used as an example of verification test. The results show that the numerical convergences accord with the verification criterion. Though the test experiment indicated the validity of our method for alleviating the oracle problem, further research might explore the application of our method to test other programs for solving second-order elliptic partial differential equations.
Data Availability
The data of HARMONY-2 code used to support the findings of this study are available from the corresponding author upon request. The data of Table 1 supporting the calculation of HARMONY-2 in this paper are from previously reported studies and datasets, which have been cited in [14]. The processed data are available from the corresponding author upon request. The data of Tables 2, 3, 4, and 5 are the calculation results of HARMNOY-2, which are used to support the findings of this study and are available from the corresponding author upon request. The data of Figures 1–6 are the calculation results of HARMNOY-2, which are used to support the findings of this study and are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported by the National Natural Science Foundation of China under Grant No. 11805093. The authors would like to thank Dr. Xie Jinsen of University of South China for providing advice on program of HARMONY-2.