Abstract
The discrete ordinates method (SN) is one of the mainstream methods for neutral particle transport calculations. Assessing the quality of the numerical solution and controlling the discrete error are essential parts of large-scale high-fidelity simulations of nuclear systems. Three error estimators, a two-mesh estimator, a residual-based estimator, and a dual-weighted residual estimator, are derived and implemented in the ARES transport code to evaluate the error of zeroth-order spatial discretization for SN equations. The difference in scalar fluxes on coarse and fine meshes is adopted to indicate the error in the two-mesh method. To avoid zero residual in zeroth-order discretization, angular fluxes within one cell are reconstructed by Legendre polynomials. The error is estimated by inverting the discrete transport operator using the estimated directional residual as an anisotropic source. The inner product of the forward directional residual and the adjoint angular flux is employed to quantify the error in quantities of interest which can be denoted by a linear functional of forward angular flux. Method of Manufactured Solutions (MMS) is adopted to generate analytical solutions for SN equation with scattering and the determined true error is used to evaluate the effectivity of these estimators. Promising results are obtained in the numerical results for both homogeneous and heterogeneous cases. The larger error region is well captured and the average effectivity index for the local error estimation is less than unity. For the series test problems, the estimated goal quantity error can be contained within an order of magnitude around the exact error.
1. Introduction
The Boltzmann transport equation [1] is central to the neutral particle transport problems as encountered in the field of nuclear engineering community. However, the numerical solution of this transport equation for realistic three-dimensional radiation shielding problems still remains a challenge because of the high-dimensional nature of the solution space and the nonsmoothness of the exact solution. In the deterministic transport calculation, a large part of errors stem from the numerical discretization of all its variables; therefore it is important to develop an approximate measure of these discretization errors. With the estimated errors, one can get the confidence level of the numerical solution and avoid unreasonable discretization. Also, the estimated error can be used to drive an adaptive mesh refinement (AMR) algorithm.
Over the past 20 years, the concept of adaptivity has been widely employed to obtain accurate solutions to the transport equation with fewer unknowns. As an essential part of adaptivity, error estimation remains active research area. In 1973, Madsen developed a method for calculating rigorous a posteriori error bounds for steady-state neutron transport problems [2]. In 1998, Jessee adopted the normalized gradient of the radiant intensity to measure the local truncation error in the radiative transport equation and developed an AMR algorithm with this measurement [3]. In 2003, the value of the mean free path in a certain cell was adopted to drive mesh refinement in the code Styx [4]. For the discrete ordinates nodal transport methods, Azmy derived a posteriori error estimator as a function of volume residual and surface jumps, which were approximated by solving local problems on a cell by cell basis. Also, the global norm error is split into a local error indicator to account for the error of each element [5, 6]. Based on the arbitrary high-order discontinuous Galerkin method, Fournier presented an error estimator for the SN transport equation and investigated its adoption to nonconforming meshes and a variable polynomial order [7]. To drive the AMR process, Wang used an error estimator based on a two-mesh approach and a jump-based error indicator. Also, he evaluated the error in the quantity of interest (QoI) with the direct and adjoint errors [8]. For the DGFEM-0 spatial discretization, Azmy estimated the residual using Taylor expansion and inverted the discrete transport operator using the residual as a source to estimate the error [9]. Recently, a dual-weighted residual (DWR) method was employed to estimate the error of the QoI. In 2011, Lathouwers combined the forward residuals and adjoint error to give accurate error measures for the effective multiplication factor in two-dimensional reactor problems [10]. For the diamond difference (DD) scheme, Jeffers derived DWR method to estimate the error in the QoI arising from the spatial discretization for fixed source and eigenvalue neutron transport problems [11].
In this paper, we assess the accuracy of the error estimator implemented in the ARES [12] transport code. The reminder of this paper is organized as follows. In Section 2, we briefly present the theory of the error estimator implemented in ARES transport code, including the two-mesh estimator, the residual-based estimator, and DWR method. Test problem specifications are described and numerical results are discussed in Section 3. Also, we provide concluding remarks in Section 4.
2. Theory of Error Estimator
The discrete ordinate method (SN) is employed to deal with the discretization of angular variables. And the monoenergetic steady-state SN transport equation with isotropic scattering iswhere is the angular flux at position along the discrete direction , and are the macroscopic total and scattering cross sections, respectively, and is the fixed-source term. The scalar flux is calculated usingwhere is the quadrature weight corresponding to the discrete direction and is the total number of discrete directions.
Weighting equation (1) by a series of test function and integrating over each mesh give a series of spatial moment balance equations:where denotes the outward normal direction of mesh boundary .
Multiple spatial discretization schemes are available in ARES. The traditional low-order schemes, including DD, weighted diamond difference (WD) [13], and step characteristic (SC) [14, 15], are based on the cell-balance equation, which can be obtained by setting in (3). Additional auxiliary equations, which assume the behavior of the angular fluxes across the spatial cell, are required to close the cell-balance equation. For example, the DD scheme assumes that the cell-averaged flux is a linear average of opposing surface-averaged fluxes, whereas the SC scheme analytically inverts the transport operator along the discrete directions within one cell. The finite-element method [16, 17] approximates the angular fluxes within one cell by a linear combination of trial functions .Using identical test spaces and trial spaces yields the Galerkin approximation. For example, the linear discontinuous (LD) scheme expands the angular flux in the trial function space , whereas the trilinear discontinuous (TLD) scheme is defined in the trial function space . In this paper, we mainly focus on the error estimators implemented to the zeroth-order spatial discretization.
2.1. Two-Mesh Estimator
Considering the justification that spatial discretization errors decrease as the spatial cell is refined in the asymptotic regime, one can use the difference in the solution obtained on the mesh (coarse) and mesh (refined) to estimate spatial discretization errors. Generally, this estimator underestimates the true error but captures the error distribution accurately. Another drawback is that the accuracy of this estimator is also affected by the occurrence of irregularities [18] in angular fluxes, which excludes some spatial cells from the asymptotic regime.
In the transport sweep process, angular fluxes are not stored; therefore the angular integrated quantity, namely, the scalar flux, is chosen as the basis of this estimator. In this paper, two different calculation schemes are investigated:where are the estimated errors in each cell, and are scalar fluxes calculated on the mesh and mesh, and is the scalar flux projected from the mesh to the mesh. The projection process is based on zeroth spatial moment conservation for low-order discretization. In scheme 1, the absolute value is taken based on the coarse mesh, whereas the absolute value is taken on fine mesh in scheme 2.
2.2. Residual-Based Estimator
We rewrite (1) in operator form:where represents the spatially continuous transport operator and denotes the scattering operator. and are the spatially continuous angular fluxes and fixed source, respectively. Discretizing (7) over the mesh yieldswhere represents the discrete transport operator, is the numerical solution obtained on the mesh, and is the fixed source projected onto the discrete domain. The spatial discretization error on the mesh is defined as , with being the exact solution to the spatial continuous equation projected onto the discrete domain. Noting the linearity of the transport operator and the scattering operator , we obtainWith the residual defined as , the error can then be obtained by solving (9) with the same iterative method for the transport equation.
To calculate residual in the zeroth-order spatial discretization schemes, we approximately reconstruct the profiles of angular fluxes on the mesh with normalized Legendre polynomials:Here, , where , are the th-order normalized Legendre polynomials. And is the vector of coefficients and is obtained through the cell-averaged and six face-averaged angular fluxes. With these assumed profiles in hand, we can calculate the residual on the mesh from:whereFrom (9), we can see that inverting the transport operator on the directional residual returns the estimated error. Note that the integral of the directional residual over each mesh is zero. To avoid the zero-residual term, we solve (9) on the mesh:In the zeroth-order scheme, is the average residual value on the mesh and is obtained by locally integrating on the mesh. Considering that the residual may be negative or positive, the DD scheme is employed to solve (13). The following two schemes are investigated to project the error from the mesh to the mesh:In scheme 1, the error on the mesh is integrated over the mesh first and then the absolute value is taken as the estimated error on the mesh. In scheme 2, the error on the mesh is estimated by integrating the absolute value of the error on the mesh.
2.3. Dual-Weighted Residual Estimator
In realistic shielding calculations, an accurate solution may not be required throughout the entire computational domain but only some linear functionals of the angular fluxes, such as the detector reaction rate and the average scalar flux in the subdomain. The linear functionals of the solution can be denoted by :To estimate the error in QoI, the following adjoint transport equation must be solved firstly:where is the adjoint transport operator and denotes adjoint angular fluxes. The adjoint angular flux is commonly referred to as the importance function of some goals, which is dependent on the adjoint source term .
As implied in (16), the adjoint source term is dependent on the goal of our calculation. With the relationship of forward and adjoint transport equation, the error in QoI can be derived as follows:This expression indicates that the error in the quantity of interest can be estimated by the inner product of the adjoint fluxes and the residual of forward transport equation. In our zeroth-order schemes implementation, the adjoint angular fluxes are also solved on the mesh and reconstructed in the same form as that in (10). where . Then, (18) is deduced toFrom the orthogonality of the Legendre polynomials, the inner product of the residual and adjoint fluxes over one cell along one direction are found to be
3. Numerical Results and Discussion
For an assessment of the accuracy of these error estimators, two series test problems are set up and the reference solutions are obtained by the MMS [19, 20] method. The two-mesh estimator and residual-based estimator are used to obtain local errors, while the DWR estimator is employed to estimate the error of region-integrated fluxes. Also, some measurements are introduced to quantify the effectiveness of the estimated errors.
3.1. Test Problem Description and Estimator Rubrics
In the test problem suites, a box-shaped domain including two alternate variant materials is considered and depicted in Figure 1. The geometry configuration has dimension of 20 cm × 20 cm × 20 cm. In the homogeneous problem, the total cross section is adopted for the two materials, whereas in the heterogeneous problem, the two different total cross sections and are employed for the two materials. A vacuum boundary condition is applied to the , , and faces of the domain and a reflecting boundary condition is applied to the , , and faces. Uniform 1 cm × 1 cm × 1 cm meshes and PNTN S16 quadrature sets [21] are used in the calculation performed in this paper.
To evaluate the quality of the estimate methods, the MMS is adopted to generate the spatial continuous solution. The exact solution of the SN equations without scattering can be determined along one angular direction.where is on the boundary of the domain and can be located by ray-tracing along the direction from the point . The combined source , including the fixed external source and the scattering source, is fixed and listed in Table 1. Also, the optical distance is represented by integrating the total cross section from to .With the purpose of generating results that correspond to the spatial discretization, we integrated the manufactured angular fluxes over each mesh.To construct the manufactured solution with scattering, the scattering source portion is subtracted from .The manufactured solution to the SN equation with scattering can now be stated as follows: the angular flux, , is the solution to the SN equations when the fixed external source equals , total cross section is , and the scattering cross section is . In our implementation, the integral over the spatial mesh, (24), is calculated numerically.
Having obtained the reference solution, we quantify the effectiveness of the error estimators by defining the following measurements. The effectivity index is defined locally in every cell.where is the estimated error and is the true error, both associated with cells . When the estimated error agrees with the true error, the value of equals zero. A positive value means an overestimation, whereas a negative value signifies an underestimation. To evaluate the accuracy of the error estimator in the global sense, the average effectivity is constructed and considered.The value gives a simple average of the absolute values of the effectivity index, whereas the value weights the effectivity index with the true error, thereby capturing mainly the accuracy for regions with large true errors.
3.2. Local Error Estimator Results
We estimated the spatial discretization errors using the two-mesh estimator and the residual-based estimator for all combinations of the following cases: homogeneous and heterogeneous cases, various zeroth-order spatial discretization schemes including DD, DZ, DTW, EDW, and SC, and scattering ratio .
Figure 2 depicts the histogram plots of the scalar fluxes effectivity for two-mesh estimator. Two calculation schemes, (5) and (6), were adopted and investigated for the homogeneous case without scattering when DTW spatial discretization is employed. In scheme 1, the effectivity index peaks slightly less than zero, indicating that the estimator underestimates the error but captures the error distribution well. The underestimation is expected because of the nature of two-mesh methods. By taking the absolute difference on the mesh, scheme 2 becomes more conservative. From Figure 2(b), the effectivity index spans more orders of magnitude, meaning that it performs poorly in producing the exact error distribution. Color-coding is employed to quantify the magnitude of the true error. For regions where the true error is relatively large, scheme 1 is able to estimate the error accurately, whereas scheme 2 evidently overestimates the error.
(a)
(b)
Figure 3 shows the scalar fluxes effectivity results for the residual-based estimator. The heterogeneous case without scattering is considered and DTW spatial discretization scheme is employed. A similar effectivity index distribution is obtained in the two schemes. By comparing (14) and (15), we conjecture that the error in each cell maintains the same positivity for most cells, although the residual in each cell does not have the same positivity. For large true error regions, scheme 2 is slightly more conservative than scheme 1.
(a)
(b)
To globally compare the effectiveness of these estimators, Figure 4 depicts the averaged effectivity index as a function of the scattering ratio for the zeroth-order spatial discretization implemented in ARES. In Figure 4, Case1 (C1) denotes the heterogeneous case, whereas Case2 (C2) is the homogeneous case. Estimator 1 (E1) refers to the two-mesh estimator, and Estimator 2 (E2) is the residual-based estimator. The averaged effectivities and are calculated using (27) and (28), respectively.
(a)
(b)
(c)
(d)
Except for the high scattering region, both of these estimators generally perform better in the homogeneous case than in the heterogeneous case. This is because the strong discontinuity in the cross section increases the singularity of the angular fluxes, posing a challenge for both spatial discretization and estimators. For the two-mesh estimator, it has been pointed out that the local accuracy does not increase with mesh refining in the regions where angular flux is discontinuous [18, 22]. Thus, the effectivity of the two-mesh estimator is influenced by the unphysical oscillations or numerical diffusion occurring in the solution on the mesh and the mesh. For the residual-based estimator, the unsmooth angular flux increases the difficulty in approximating the residual with polynomials and the strong discontinuity in the cross section may also decrease the accuracy of inverting the transport operator on the residual.
From Figure 4, the true error weighted averaged effectivity index is obviously smaller than except for the SC discretization scheme, implying that these estimators capture larger errors more accurately. This property is preferable in AMR algorithm. When comparing these two estimators, we find that the two-mesh estimator outperforms residual-based estimator in the DD scheme, whereas the residual-based estimator performs better in the SC scheme. In the EDW scheme, these two estimators are competitive in terms of the effectivity. In the DTW scheme, two-mesh estimator is better for the homogeneous case and the residual-based estimator is slightly better in the heterogeneous case.
With increasing scattering ratio, the average effectivity increases significantly in DD, DTW, and EDW scheme. In the SC, there is no evident relation between the effectivity and scattering ratio. Overall, estimators employed in the EDW achieve relatively better accuracy, with the maximum average effectivity less than 0.4. When these estimators are used in the DD or DTW schemes, the average effectivity reaches nearly 0.8 for high scattering case.
3.3. DWR Estimator Results
For the DWR estimator, we estimate the error of the region-integrated scalar fluxes and compare the estimated error with the true error. The goal region is located in . By setting the adjoint source equal to 1 in this region, we can obtain the adjoint angular fluxes by solving adjoint transport equation. Then, the goal error is obtained by calculating the inner product of the forward directional residual and the adjoint angular fluxes with (20) and (21).
Figure 5 depicts the ratio of the estimated goal error and the true goal error for the homogeneous case. The same spatial discretization scheme is used in both the forward and adjoint calculations. The DD is not included in Figure 5 because its performance in predicting the adjoint angular fluxes is poor, especially in the zero adjoint source regions. For low scattering problems, the ideal ratio values of unity, meaning the estimated error is close to the true error, were obtained in the DTW, EDW, and SC schemes. Nonetheless, with increasing scattering ratio, the error is overestimated in the SC scheme and underestimated in the DTW and EDW schemes. Conservative results, the estimated error is apparently larger than the exact error, are produced in the DZ scheme. Overall, the estimated error is contained within an order of magnitude around the true error.
To investigate the effect of spatial discretization when solving the adjoint transport problems, we solved the forward transport problem with SC spatial discretization and the adjoint transport problem with the five spatial discretization schemes. The results for the homogeneous and heterogeneous cases are presented in Figure 6. Overall, the estimator captures the error well, although the exact relative error of the region-integrated scalar fluxes is extremely small, ranging from 0.006% to 0.349% in the homogeneous case and from 0.118% to 0.672% in the heterogeneous case. Except for the DD adjoint spatial discretization in the homogeneous case, similar trends are obtained. With increasing scattering ratio, the true relative error increases monotonically and the estimator becomes more conservative. From comparing the performances of the DD and DZ adjoint spatial discretizations in the homogeneous case, we conjecture that negative adjoint angular fluxes occur and induce oscillation of adjoint angular fluxes. Therefore, the accuracy of the inner product of the adjoint angular fluxes and the directional residual can be influenced. However, in the heterogeneous cases, the DD adjoint spatial discretization has a negligible effect on the accuracy of the estimator, which is attributed to the negligible contribution that the inaccurate adjoint angular fluxes in the source-free region makes to the estimated goal error. Also, the exact relative error of the heterogeneous case is significantly larger than the homogeneous case, especially for low scattering problems. Therefore, the effect of spatial discretization is almost imperceptible.
(a)
(b)
4. Conclusion
This paper briefly introduces the spatial discretization schemes implemented in ARES transport code to treat the SN equations and presents error estimators employed in ARES to estimate the spatial discretization error for zeroth-order schemes. The profiles of angular fluxes within each cell were approximately reconstructed by the Legendre polynomials, so that the residual could be obtained. MMS was adopted to generate analytical solutions to the SN equations with scattering, so that the effectivity of these estimators could be investigated by comparing the estimated error and exact true error.
These estimators achieved promising results for the series test problems, which contain homogeneous and heterogeneous cases with varying scattering ratio. For the local error estimators, the larger error region was captured effectively and the average effectivity index was contained to less than unity. Also, the estimated goal error was contained within an order of magnitude around the true error. Generally, best results occurred in the homogeneous case without scattering, and with increasing scattering ratio, the accuracy of the estimators diminished. Moreover, the estimators were influenced by the exact angular flux profiles and the employed spatial discretization techniques; for example, the two-mesh estimator performs better when DD spatial discretization was used, and residual-based estimator is better in SC.
The proposed estimators have the potential to drive the AMR process and control error when the zeroth-order spatial discretization is adopted. For further improvements of these estimators, the representation of the residual should be investigated for nonsmooth angular fluxes and the numerical oscillation needs to be ameliorated in the residual-to-error transport calculations and adjoint transport calculations. And another major challenge is to investigate the couple spatial and angular discretization error, develop estimator for the coupled error, and test with international shielding benchmark problems.
Data Availability
Because our deterministic transport code is not public, the data is not suitable to share online. However, all the data related to this paper can be obtained by contacting the authors reasonably for academic purposes.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (11505059 and 11575061) and the Fundamental Research Funds for the Central Universities (2017XS091).