# Richardson Extrapolation-Based Verification Method of Scientific Calculation Program without the Oracles: A Case Study

Authors: Shiyu Yan; Xiaohua Yang; Guodong Cheng; Hua Liu
Journal: Mathematical Problems in Engineering
 (2019)

---

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


---

## Body

## 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 inH0 and H1 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 Section2 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).(1)-∇·Dgr∇ϕgr+ΣR,grϕgr=∑g′=1,g′≠gGΣg′→grϕg′r+χgkeffQr,g=1,2,…,G,(2)Qr=∑g′=1GυΣfrg′ϕg′r,where g is the energy group, Dg(r) is the diffusion coefficient of the energy group g, ΣR,g(r) is the removal cross section of the energy group g, Σg′→g(r) is the scattering cross section from energy group g′ to g, υg is the average number of neutrons emitted by fission of the energy group g, (∑f(r))g′ is the fission cross section of the energy group g′, ϕg(r) is the neutron scalar flux of the energy group g, χg is the integrated fission spectrum of the energy group g, and keff is the effective multiplication factor of reactor. Q(r) 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 Ωi,  i=1,2,…,M. For r∈Ω with an outer boundary conditions, we have(3)Dg∂ϕgr∂n+ϕgr=0,r∈∂Ω,where n is the outer normal to the boundary ∂Ω; at an interface between dissimilar regions Ωi with continuous conditions, we have(4)ϕgr+=ϕgr-,and(5)Dg∂ϕgr+∂n=Dg∂ϕgr-∂n,r∈∂Ωi,where r+ and r- represent the value of both sides of the interface, respectively. The quantities ΣR,g(r),Dg(r),Σg′→g(r),χg,υΣf(r) are constant in each Ωi.To simplify the manipulations, (1) and (2) are expressed in a matrix system by the following form:(6)Mφ=1keffFφ,where the vector of the neutron flux is φ=[ϕ1,ϕ2,…,ϕG]T, the coefficient matrix M=(ai,j)G,G is defined with(7)ai,j=-Di∇2+ΣR,i,i=j-Σj′,i,i≠j′,and the coefficient matrix F=(fi,j)G,G is defined with fi,j=χiυjΣfj.


### 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φ(0) and keff(0), they are substituted into the right hand of (6) to obtain the neutron flux φ(1). The upper labels represent the number of iterations.(II)
For the i-th iteration, (6) is written as(8)Mφi+1=1keffiSi,whereS=Fφ.(III)
Updating the value for newkeff and the new flux vector φ, keff(i+1) and S(i+1) are calculated by the following equations, respectively,(9)Si+1=Fφi+1,keffi+1=keffi∫ΩSi+1dΩ∫ΩSidΩ,with the stop criterion as follows:(10)keffi+1-keffikeffi+1<ε1,∫ΩSi+1dΩ-∫ΩSidΩ∫ΩSi+1dΩ<ε2,whereε1 and ε2 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(11)-∇Dr∇ϕr+ΣRrϕr=Sr,r∈Ω(12)Dr∂ϕr∂n+ϕr=0,r∈∂Ω.Here,S(r) 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.


## 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).(1)-∇·Dgr∇ϕgr+ΣR,grϕgr=∑g′=1,g′≠gGΣg′→grϕg′r+χgkeffQr,g=1,2,…,G,(2)Qr=∑g′=1GυΣfrg′ϕg′r,where g is the energy group, Dg(r) is the diffusion coefficient of the energy group g, ΣR,g(r) is the removal cross section of the energy group g, Σg′→g(r) is the scattering cross section from energy group g′ to g, υg is the average number of neutrons emitted by fission of the energy group g, (∑f(r))g′ is the fission cross section of the energy group g′, ϕg(r) is the neutron scalar flux of the energy group g, χg is the integrated fission spectrum of the energy group g, and keff is the effective multiplication factor of reactor. Q(r) 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 Ωi,  i=1,2,…,M. For r∈Ω with an outer boundary conditions, we have(3)Dg∂ϕgr∂n+ϕgr=0,r∈∂Ω,where n is the outer normal to the boundary ∂Ω; at an interface between dissimilar regions Ωi with continuous conditions, we have(4)ϕgr+=ϕgr-,and(5)Dg∂ϕgr+∂n=Dg∂ϕgr-∂n,r∈∂Ωi,where r+ and r- represent the value of both sides of the interface, respectively. The quantities ΣR,g(r),Dg(r),Σg′→g(r),χg,υΣf(r) are constant in each Ωi.To simplify the manipulations, (1) and (2) are expressed in a matrix system by the following form:(6)Mφ=1keffFφ,where the vector of the neutron flux is φ=[ϕ1,ϕ2,…,ϕG]T, the coefficient matrix M=(ai,j)G,G is defined with(7)ai,j=-Di∇2+ΣR,i,i=j-Σj′,i,i≠j′,and the coefficient matrix F=(fi,j)G,G is defined with fi,j=χiυjΣfj.


## 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φ(0) and keff(0), they are substituted into the right hand of (6) to obtain the neutron flux φ(1). The upper labels represent the number of iterations.(II)
For the i-th iteration, (6) is written as(8)Mφi+1=1keffiSi,whereS=Fφ.(III)
Updating the value for newkeff and the new flux vector φ, keff(i+1) and S(i+1) are calculated by the following equations, respectively,(9)Si+1=Fφi+1,keffi+1=keffi∫ΩSi+1dΩ∫ΩSidΩ,with the stop criterion as follows:(10)keffi+1-keffikeffi+1<ε1,∫ΩSi+1dΩ-∫ΩSidΩ∫ΩSi+1dΩ<ε2,whereε1 and ε2 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(11)-∇Dr∇ϕr+ΣRrϕr=Sr,r∈Ω(12)Dr∂ϕr∂n+ϕr=0,r∈∂Ω.Here,S(r) 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:(13)-ddrDgrdϕgdr+ΣR,grϕgr=Sgr,r∈Ωg=1,2,…,G,whereΩ={r:0<r<R}, ∂Ω={0,R}, ϕg in ∂Ω is zero, and ϕg 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.Figure 1
Space discrete in one-dimensional X geometry.LetΔx be the split size and NX 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). (14)agiϕg,mi-1+bgiϕg,mi+cgiϕg,mi+1-SSg,mi=FSg,mig=1,2,…,G,i=1,2,…,NXEquation (14) is written as an overrelaxation iterative format: (15)ϕg,kK=ωbgiFSg,m-1i+SSg,mK-1i-agiϕg,mKi-1-cgiϕg,mK-1i+1+1-ωϕg,mK-1i,where superscript K is the inner iteration count and m is the outer iteration count of calculation algorithm, where(16)agi=-Dgi-1/2Δxi-1,(17)bi=12Σt,gi+12Δxi+12Σt,gi-12Δxi-1-agi-cgi,(18)cgi=-Dgi+1/2Δxi,(19)FSgi=χgi∑g′=1GυΣfg′ϕg′,m-1iΔxi,(20)SSg,mi=∑g′=1GΣg′→giϕg′,miΔx.Boundary conditions are as follows:(21)agiB-=-2DgiB-ΔxiB-,ϕgiB--1=0,(22)agiB+=-2DgiB+ΔxiB+,ϕgiB++1=0.Calculation flow of the program algorithm for (13) is shown in Figure 2.Figure 2
Calculation flow chart of the algorithm of HARMNOY-2.


### 3.2. Richardson Extrapolation

Based on the Richardson Extrapolation method, the convergence of the numerical solution in (1) and (2) was investigated. Let φh 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 Hα(Ω), we assume φh∈Hα(Ω), and we have (23)φ-φhHαΩ=oh2-α.Specifically, the numerical solutionφh is calculated by the same discretization scheme with three different grid sizes: h1=h, h2=h/2, and h3=h/4, respectively.Ifα=1, the error can be rewritten by(24)φh=φexact+c0h+oh2,(25)φh/2=φexact+c0h2+oh22,(26)φh/4=φexact+c0h4+oh42.Equations (24) to (26) yield(27)φh-φh/2H1=c0h2,(28)φh/2-φh/4H1=c0h4.Then, we get(29)φh-φh/2H1φh/2-φh/4H1≅2.If α=0, we have (30)φh=φexact+c1h2+oh3,(31)φh/2=φexact+c1h22+oh23,(32)φh/4=φexact+c1h42+oh43,and that the mesh is refined or coarsened by a factor of two. Similarly, we have the following.(33)φh-φh/2H0φh/2-φh/4H0≅4Based on the above derivation, the numerical solution convergence ofH0 and H1 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 H0 and H1, respectively. The last one is the Max-eigenvalues convergence with the mesh refinement.


## 3.1. The Finite Difference Method

We consider an example of one-dimensional and multigroup neutron diffusion equations:(13)-ddrDgrdϕgdr+ΣR,grϕgr=Sgr,r∈Ωg=1,2,…,G,whereΩ={r:0<r<R}, ∂Ω={0,R}, ϕg in ∂Ω is zero, and ϕg 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.Figure 1
Space discrete in one-dimensional X geometry.LetΔx be the split size and NX 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). (14)agiϕg,mi-1+bgiϕg,mi+cgiϕg,mi+1-SSg,mi=FSg,mig=1,2,…,G,i=1,2,…,NXEquation (14) is written as an overrelaxation iterative format: (15)ϕg,kK=ωbgiFSg,m-1i+SSg,mK-1i-agiϕg,mKi-1-cgiϕg,mK-1i+1+1-ωϕg,mK-1i,where superscript K is the inner iteration count and m is the outer iteration count of calculation algorithm, where(16)agi=-Dgi-1/2Δxi-1,(17)bi=12Σt,gi+12Δxi+12Σt,gi-12Δxi-1-agi-cgi,(18)cgi=-Dgi+1/2Δxi,(19)FSgi=χgi∑g′=1GυΣfg′ϕg′,m-1iΔxi,(20)SSg,mi=∑g′=1GΣg′→giϕg′,miΔx.Boundary conditions are as follows:(21)agiB-=-2DgiB-ΔxiB-,ϕgiB--1=0,(22)agiB+=-2DgiB+ΔxiB+,ϕgiB++1=0.Calculation flow of the program algorithm for (13) is shown in Figure 2.Figure 2
Calculation flow chart of the algorithm of HARMNOY-2.


## 3.2. Richardson Extrapolation

Based on the Richardson Extrapolation method, the convergence of the numerical solution in (1) and (2) was investigated. Let φh 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 Hα(Ω), we assume φh∈Hα(Ω), and we have (23)φ-φhHαΩ=oh2-α.Specifically, the numerical solutionφh is calculated by the same discretization scheme with three different grid sizes: h1=h, h2=h/2, and h3=h/4, respectively.Ifα=1, the error can be rewritten by(24)φh=φexact+c0h+oh2,(25)φh/2=φexact+c0h2+oh22,(26)φh/4=φexact+c0h4+oh42.Equations (24) to (26) yield(27)φh-φh/2H1=c0h2,(28)φh/2-φh/4H1=c0h4.Then, we get(29)φh-φh/2H1φh/2-φh/4H1≅2.If α=0, we have (30)φh=φexact+c1h2+oh3,(31)φh/2=φexact+c1h22+oh23,(32)φh/4=φexact+c1h42+oh43,and that the mesh is refined or coarsened by a factor of two. Similarly, we have the following.(33)φh-φh/2H0φh/2-φh/4H0≅4Based on the above derivation, the numerical solution convergence ofH0 and H1 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 H0 and H1, 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.Table 1
Parameters of SLAB_1D_2G one-dimensional slab benchmark problem.
Energy group D / c m Σt/cm-1 v Σ f / c m - 1 χ p / χ d Σ 1 → 2 / c m - 1 1 1.264E-00 8.154E-03 0.000E-00 1.00E-00 7.368E-03 2 9.328E-01 4.100E-03 4.562E-03 0.00E-00 7.368E-03Figure 3
Reactor core configuration of SLAB_1D_2G.(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φh=(ϕ1h,ϕ2h), where ϕ1h is the neutron scalar flux for 1st group energy, ϕ2h is the neutron flux for 2nd group energy, and h denotes the size of grid. We assume that h1=450/18=25cm, h2=h1/2, h3=h1/4, h4=h1/8, h5=h1/16, and h6=h1/32 are presented as different refined grids, respectively. The point xmid=225cm 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 Ratio1 and Ratio2 denote the numerical solution convergence of H0 and H1, respectively, as follows:(34)Ratio1=φhi-φhi+1H0φhi+1-φhi+2H0,Ratio2=φhi-φhi+1H1φhi+1-φhi+2H1.


### 4.2. Test Results and Discussions

The calculation results of the four modes eigenvalues are shown in Table2 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 ϕ1, ϕ2 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 H0 on the neutron flux value of xmid point in different meshes. In Table 3, the convergence order approximates to 4. The derivative values of neutron flux on xmid 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 keff 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
Comparison of results by the four modes eigenvalues.
n λ-eigenvalue Absolute error HARMONY-2 Lecture [14] 1 9.86988E-01 9.8698E-01 0.0000 2 9.35159E-01 9.3516E-01 1.000E-05 3 8.58979E-01 8.5898E-01 1.000E-05 4 7.69751E-01 7.6976E-01 1.000E-05Table 3
The neutron flux values ofxmid in different refined grids.
Mesh Element The neutron flux value ofxmid for 1th group energy. Ratio1 The neutron flux value ofxmid for 2th group energy. Ratio1 h 1 18 6.641297145594463E-002  6.641294098874310E-002 h 2 36 6.660342485457200E-002  6.660331582277153E-002 h 3 72 6.665076278267033E-002 4.02327280213385 6.665082494407264E-002 4.00712176556248 h 4 144 6.666264561430789E-002 3.98372454834075 6.666253389102599E-002 4.05750589616633 h 5 288 6.666565458515177E-002 3.94913485510786 6.666562905904372E-002 3.78297620234630 h 6 576 6.666633289577738E-002 4.43597775157830 6.666641383170335E-002 3.94403140799427Table 4
The derivative values of neutron flux onxmid in different refined grids.
Mesh Element The derivative of neutron flux value ofxmid for 1th group energy. Ratio2 The derivative of neutron flux value ofxmid for 2th group energy. Ratio2 h1 18 0.80710858E-04  0.8071363955E-04 h2 36 0.42252071E-04 0.19102225E+01 0.4224298523E-04 0.19106992E+01 h3 72 0.21143797E-04 0.19983199E+01 0.2118749068E-04 0.19937701E+01 h4 144 0.10581871E-04 0.19981151E+01 0.1061129701E-04 0.19966919E+01 h5 288 0.79545852E-05 0.13302858E+01 0.7929059044E-05 0.13382795E+01Table 5
The calculation results of Max-eigenvalue in different refined grids.
Mesh h 1 h 2 h 3 h 4 h 5 h 6 Element 18 36 72 144 288 576  Max-eigenvalue 9.86988E-01 9.86954E-01 9.86945E-01 9.86943E-01 9.86943E-01 9.86942E-01Figure 4
The convergence of the Max-eigenvalue.Figure 5
ϕ 1 in grid with 18 equal elements.Figure 6
ϕ 2 in grid with 18 equal elements.Table2 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.


## 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.Table 1
Parameters of SLAB_1D_2G one-dimensional slab benchmark problem.
Energy group D / c m Σt/cm-1 v Σ f / c m - 1 χ p / χ d Σ 1 → 2 / c m - 1 1 1.264E-00 8.154E-03 0.000E-00 1.00E-00 7.368E-03 2 9.328E-01 4.100E-03 4.562E-03 0.00E-00 7.368E-03Figure 3
Reactor core configuration of SLAB_1D_2G.(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φh=(ϕ1h,ϕ2h), where ϕ1h is the neutron scalar flux for 1st group energy, ϕ2h is the neutron flux for 2nd group energy, and h denotes the size of grid. We assume that h1=450/18=25cm, h2=h1/2, h3=h1/4, h4=h1/8, h5=h1/16, and h6=h1/32 are presented as different refined grids, respectively. The point xmid=225cm 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 Ratio1 and Ratio2 denote the numerical solution convergence of H0 and H1, respectively, as follows:(34)Ratio1=φhi-φhi+1H0φhi+1-φhi+2H0,Ratio2=φhi-φhi+1H1φhi+1-φhi+2H1.


## 4.2. Test Results and Discussions

The calculation results of the four modes eigenvalues are shown in Table2 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 ϕ1, ϕ2 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 H0 on the neutron flux value of xmid point in different meshes. In Table 3, the convergence order approximates to 4. The derivative values of neutron flux on xmid 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 keff 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
Comparison of results by the four modes eigenvalues.
n λ-eigenvalue Absolute error HARMONY-2 Lecture [14] 1 9.86988E-01 9.8698E-01 0.0000 2 9.35159E-01 9.3516E-01 1.000E-05 3 8.58979E-01 8.5898E-01 1.000E-05 4 7.69751E-01 7.6976E-01 1.000E-05Table 3
The neutron flux values ofxmid in different refined grids.
Mesh Element The neutron flux value ofxmid for 1th group energy. Ratio1 The neutron flux value ofxmid for 2th group energy. Ratio1 h 1 18 6.641297145594463E-002  6.641294098874310E-002 h 2 36 6.660342485457200E-002  6.660331582277153E-002 h 3 72 6.665076278267033E-002 4.02327280213385 6.665082494407264E-002 4.00712176556248 h 4 144 6.666264561430789E-002 3.98372454834075 6.666253389102599E-002 4.05750589616633 h 5 288 6.666565458515177E-002 3.94913485510786 6.666562905904372E-002 3.78297620234630 h 6 576 6.666633289577738E-002 4.43597775157830 6.666641383170335E-002 3.94403140799427Table 4
The derivative values of neutron flux onxmid in different refined grids.
Mesh Element The derivative of neutron flux value ofxmid for 1th group energy. Ratio2 The derivative of neutron flux value ofxmid for 2th group energy. Ratio2 h1 18 0.80710858E-04  0.8071363955E-04 h2 36 0.42252071E-04 0.19102225E+01 0.4224298523E-04 0.19106992E+01 h3 72 0.21143797E-04 0.19983199E+01 0.2118749068E-04 0.19937701E+01 h4 144 0.10581871E-04 0.19981151E+01 0.1061129701E-04 0.19966919E+01 h5 288 0.79545852E-05 0.13302858E+01 0.7929059044E-05 0.13382795E+01Table 5
The calculation results of Max-eigenvalue in different refined grids.
Mesh h 1 h 2 h 3 h 4 h 5 h 6 Element 18 36 72 144 288 576  Max-eigenvalue 9.86988E-01 9.86954E-01 9.86945E-01 9.86943E-01 9.86943E-01 9.86942E-01Figure 4
The convergence of the Max-eigenvalue.Figure 5
ϕ 1 in grid with 18 equal elements.Figure 6
ϕ 2 in grid with 18 equal elements.Table2 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.



---

*Source: 5782146-2019-01-30.xml*