# Analysis of Spatial Discretization Error Estimators Implemented in ARES Transport Code forSN Equations

Authors: Liang Zhang; Bin Zhang; Cong Liu; Yixue Chen
Journal: Science and Technology of Nuclear Installations
 (2018)

---

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


---

## Body

## 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 is(1)Ω→m·∇ψmr→+Σtr→ψmr→=Σsr→ϕr→+Smr→,m=1,⋯,Mwhere ψmr→ is the angular flux at position r→ along the discrete direction Ω→m, Σtr→ and Σsr→ are the macroscopic total and scattering cross sections, respectively, and Smr→ is the fixed-source term. The scalar flux is calculated using(2)ϕr→=∑m=1Mwmψmr→where wm is the quadrature weight corresponding to the discrete direction Ω→m and M is the total number of discrete directions.Weighting equation (1) by a series of test function vli→r→ and integrating over each mesh Qi→=xi-1,xi×yj-1,yj×zk-1,zk give a series of spatial moment balance equations:(3)∮∂Qi→Ω→m·n→Fvli→r→ψmr→dA-∫Qiψmr→Ω→m·∇vli→r→dr→+∫QiΣtr→ψmr→vli→r→dr→=∫QiΣsr→ϕr→+Smr→vli→r→dr→where n→F denotes the outward normal direction of mesh boundary ∂Qi→.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 vli→=1 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 bnr→.(4)ψmr→=∑n=1Nbnr→ψmnUsing 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 1,x,y,z, whereas the trilinear discontinuous (TLD) scheme is defined in the trial function space 1,x,y,z,xy,xz,yz,xyz. 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 theh mesh (coarse) and h/2 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:(5)Eijk,coarseh=ϕi,j,kh-18∑p=01∑q=01∑r=01ϕ2i-p,2j-q,2k-rh/2(6)Eijk,fineh=18∑p=01∑q=01∑r=01ϕ2i-p,2j-q,2k-rh→h/2-ϕ2i-p,2j-q,2k-rh/2where E are the estimated errors in each cell, ϕh and ϕh/2 are scalar fluxes calculated on the h mesh and h/2 mesh, and ϕh→h/2 is the scalar flux projected from the h mesh to the h/2 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:(7)Lψ=Sψ+Qwhere L represents the spatially continuous transport operator and S denotes the scattering operator. ψ and Q are the spatially continuous angular fluxes and fixed source, respectively. Discretizing (7) over the h mesh yields(8)Lhψh=Sψh+Qhwhere Lh represents the discrete transport operator, ψh is the numerical solution obtained on the h mesh, and Qh is the fixed source Q projected onto the discrete domain. The spatial discretization error on the h mesh is defined as εh=ψh-ψh, with ψh being the exact solution to the spatial continuous equation projected onto the discrete domain. Noting the linearity of the transport operator L and the scattering operator S, we obtain(9)Lε=Sε+Sψh+Qh-LψhWith the residual defined as Rh=Sψh+Qh-Lψh, the error can then be obtained by solving (9) with the same iterative method for the transport equation.To calculate residualRh in the zeroth-order spatial discretization schemes, we approximately reconstruct the profiles of angular fluxes on the h mesh with normalized Legendre polynomials:(10)ψmh=ΨmPHere, P=1,P1x,P2x,P1y,P2y,P1z,P2zT, where Plu,u∈x,y,z, are the lth-order normalized Legendre polynomials. And Ψm=ψm0,ψmx,1,ψmx,2,ψmy,1,ψmy,2,ψmz,1,ψmz,2 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 h mesh from:(11)Rmh=0,Rmx,1,Rmx,2,Rmy,1,Rmy,2,Rmz,1,Rmz,2Pwhere(12)Rmu,1=Σs∑m′=1Mwm′ψm′u,1-Σtψmu,1-Ωmu6ψmu,2ΔuRmu,2=Σs∑m′=1Mwm′ψm′u,2-Σtψmu,2,u∈x,y,zFrom (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 h/2 mesh:(13)Lh/2ε=Sε+Rh/2In the zeroth-order scheme, Rh/2 is the average residual value on the h/2 mesh and is obtained by locally integrating Rh on the h/2 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 h/2 mesh to the h mesh:(14)εh=∫εh/2dr(15)εh=∫εh/2drIn scheme 1, the error on the h/2 mesh is integrated over the h mesh first and then the absolute value is taken as the estimated error on the h mesh. In scheme 2, the error on the h mesh is estimated by integrating the absolute value of the error on the h/2 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 byg,ψ:(16)g,ψ=∫D∑m=1Mgmψmdr→To estimate the error in QoI, the following adjoint transport equation must be solved firstly:(17)L∗ψ∗=S∗ψ∗+gwhere L∗ 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 g.As implied in (16), the adjoint source term g 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:(18)g,ψ-g,ψh=g,ψ-ψh=L∗ψ∗-S∗ψ∗,ψ-ψh=ψ∗,Lψ-ψh-Sψ-ψh=ψ∗,Q+Sψh-LψhThis 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 h mesh and reconstructed in the same form as that in (10). (19)ψm∗,h=Ψm∗Pwhere Ψm∗=ψm∗,0,ψm∗,x,1,ψm∗,x,2,ψm∗,y,1,ψm∗,y,2,ψm∗,z,1,ψm∗,z,2. Then, (18) is deduced to(20)g,ψ-g,ψh=∫D∑m=1Mψm∗Rmdr→From 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(21)∫zk-1zkdz∫yj-1yjdy∫xi-1xidxRmψm∗=Vijk3Rmx,1ψm∗,x,1+Rmy,1ψm∗,y,1+Rmz,1ψm∗,z,1+Vijk5Rmx,2ψm∗,x,2+Rmy,2ψm∗,y,2+Rmz,2ψm∗,z,2


## 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 theh mesh (coarse) and h/2 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:(5)Eijk,coarseh=ϕi,j,kh-18∑p=01∑q=01∑r=01ϕ2i-p,2j-q,2k-rh/2(6)Eijk,fineh=18∑p=01∑q=01∑r=01ϕ2i-p,2j-q,2k-rh→h/2-ϕ2i-p,2j-q,2k-rh/2where E are the estimated errors in each cell, ϕh and ϕh/2 are scalar fluxes calculated on the h mesh and h/2 mesh, and ϕh→h/2 is the scalar flux projected from the h mesh to the h/2 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:(7)Lψ=Sψ+Qwhere L represents the spatially continuous transport operator and S denotes the scattering operator. ψ and Q are the spatially continuous angular fluxes and fixed source, respectively. Discretizing (7) over the h mesh yields(8)Lhψh=Sψh+Qhwhere Lh represents the discrete transport operator, ψh is the numerical solution obtained on the h mesh, and Qh is the fixed source Q projected onto the discrete domain. The spatial discretization error on the h mesh is defined as εh=ψh-ψh, with ψh being the exact solution to the spatial continuous equation projected onto the discrete domain. Noting the linearity of the transport operator L and the scattering operator S, we obtain(9)Lε=Sε+Sψh+Qh-LψhWith the residual defined as Rh=Sψh+Qh-Lψh, the error can then be obtained by solving (9) with the same iterative method for the transport equation.To calculate residualRh in the zeroth-order spatial discretization schemes, we approximately reconstruct the profiles of angular fluxes on the h mesh with normalized Legendre polynomials:(10)ψmh=ΨmPHere, P=1,P1x,P2x,P1y,P2y,P1z,P2zT, where Plu,u∈x,y,z, are the lth-order normalized Legendre polynomials. And Ψm=ψm0,ψmx,1,ψmx,2,ψmy,1,ψmy,2,ψmz,1,ψmz,2 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 h mesh from:(11)Rmh=0,Rmx,1,Rmx,2,Rmy,1,Rmy,2,Rmz,1,Rmz,2Pwhere(12)Rmu,1=Σs∑m′=1Mwm′ψm′u,1-Σtψmu,1-Ωmu6ψmu,2ΔuRmu,2=Σs∑m′=1Mwm′ψm′u,2-Σtψmu,2,u∈x,y,zFrom (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 h/2 mesh:(13)Lh/2ε=Sε+Rh/2In the zeroth-order scheme, Rh/2 is the average residual value on the h/2 mesh and is obtained by locally integrating Rh on the h/2 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 h/2 mesh to the h mesh:(14)εh=∫εh/2dr(15)εh=∫εh/2drIn scheme 1, the error on the h/2 mesh is integrated over the h mesh first and then the absolute value is taken as the estimated error on the h mesh. In scheme 2, the error on the h mesh is estimated by integrating the absolute value of the error on the h/2 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 byg,ψ:(16)g,ψ=∫D∑m=1Mgmψmdr→To estimate the error in QoI, the following adjoint transport equation must be solved firstly:(17)L∗ψ∗=S∗ψ∗+gwhere L∗ 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 g.As implied in (16), the adjoint source term g 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:(18)g,ψ-g,ψh=g,ψ-ψh=L∗ψ∗-S∗ψ∗,ψ-ψh=ψ∗,Lψ-ψh-Sψ-ψh=ψ∗,Q+Sψh-LψhThis 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 h mesh and reconstructed in the same form as that in (10). (19)ψm∗,h=Ψm∗Pwhere Ψm∗=ψm∗,0,ψm∗,x,1,ψm∗,x,2,ψm∗,y,1,ψm∗,y,2,ψm∗,z,1,ψm∗,z,2. Then, (18) is deduced to(20)g,ψ-g,ψh=∫D∑m=1Mψm∗Rmdr→From 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(21)∫zk-1zkdz∫yj-1yjdy∫xi-1xidxRmψm∗=Vijk3Rmx,1ψm∗,x,1+Rmy,1ψm∗,y,1+Rmz,1ψm∗,z,1+Vijk5Rmx,2ψm∗,x,2+Rmy,2ψm∗,y,2+Rmz,2ψm∗,z,2


## 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 Figure1. The geometry configuration has dimension of 20 cm × 20 cm × 20 cm. In the homogeneous problem, the total cross section Σt=0.5cm-1 is adopted for the two materials, whereas in the heterogeneous problem, the two different total cross sections Σt1=0.5×10-4cm-1 and Σt2=0.5cm-1 are employed for the two materials. A vacuum boundary condition is applied to the +x, +y, and +z faces of the domain and a reflecting boundary condition is applied to the -x, -y, and -z faces. Uniform 1 cm × 1 cm × 1 cm meshes and PNTN S16 quadrature sets [21] are used in the calculation performed in this paper.Figure 1
Geometric configuration for the test problem suites.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.(22)ψmMMSr→0=ψmMMSr→Be-τr→0,r→B+∫r→0r→Bdr→Qr→e-τr→0,r→where r→B is on the boundary of the domain and can be located by ray-tracing along the -Ω→m direction from the point r→0. The combined source Qr→, including the fixed external source and the scattering source, is fixed and listed in Table 1. Also, the optical distance τr→1,r→2 is represented by integrating the total cross section from r→1 to r→2.(23)τr→1,r→2=∫r→1r→2dr→Σtr→With the purpose of generating results that correspond to the spatial discretization, we integrated the manufactured angular fluxes over each mesh.(24)ψm,Qi→MMS=∫Qi→dr→ψmMMSr→To construct the manufactured solution with scattering, the scattering source portion is subtracted from Qr→.(25)qQi→MMS=QQi→-Σs,Qi→∑m=1Mwmψm,Qi→MMSThe manufactured solution to the SN equation with scattering can now be stated as follows: the angular flux, ψm,Qi→MMS, is the solution to the SN equations when the fixed external source equals qQi→MMS, total cross section is Σt,Qi→, and the scattering cross section is Σs,Qi→. In our implementation, the integral over the spatial mesh, (24), is calculated numerically.Table 1
Total cross section and combined source of the test problem suites.
Case Total cross sectionΣt (cm−1) Combined sourceQ (cm−3s−1) Material 1 Material 2 Material 1 Material 2 Homogeneous 0.5 0.5 1 1 Heterogeneous 0.5 0.5×10−4 1 1×10−4Having 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.(26)θijk=log10εijkEijkwhere εijk is the estimated error and Eijk is the true error, both associated with cells  i,j,k. 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.(27)θ-1=1N∑θijk(28)θ-2=∑θijkEijk∑EijkThe θ-1 value gives a simple average of the absolute values of the effectivity index, whereas the θ-2 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 ratioc=0.0,0.1,0.3,0.5,0.7,0.9.Figure2 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 h/2 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.Figure 2
Histogram of effectivity for two-mesh estimator.
(a) (b)Figure3 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.Figure 3
Histogram of effectivity for residual-based estimator.
(a) (b)To globally compare the effectiveness of these estimators, Figure4 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 θ-1 and θ-2 are calculated using (27) and (28), respectively.Figure 4
Comparison of the averaged effectivity for the various zeroth-order spatial discretizations.
(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 h mesh and the h/2 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 Figure4, the true error weighted averaged effectivity index θ-2 is obviously smaller than θ-1 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 inx,y,z∈0,10cm. 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).Figure5 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.Figure 5
Ratio of the estimated error to the true error for the homogeneous case.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 Figure6. 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.Figure 6
Effects of the various adjoint spatial discretizations on the accuracy of the DWR estimator.
(a) (b)


## 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 Figure1. The geometry configuration has dimension of 20 cm × 20 cm × 20 cm. In the homogeneous problem, the total cross section Σt=0.5cm-1 is adopted for the two materials, whereas in the heterogeneous problem, the two different total cross sections Σt1=0.5×10-4cm-1 and Σt2=0.5cm-1 are employed for the two materials. A vacuum boundary condition is applied to the +x, +y, and +z faces of the domain and a reflecting boundary condition is applied to the -x, -y, and -z faces. Uniform 1 cm × 1 cm × 1 cm meshes and PNTN S16 quadrature sets [21] are used in the calculation performed in this paper.Figure 1
Geometric configuration for the test problem suites.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.(22)ψmMMSr→0=ψmMMSr→Be-τr→0,r→B+∫r→0r→Bdr→Qr→e-τr→0,r→where r→B is on the boundary of the domain and can be located by ray-tracing along the -Ω→m direction from the point r→0. The combined source Qr→, including the fixed external source and the scattering source, is fixed and listed in Table 1. Also, the optical distance τr→1,r→2 is represented by integrating the total cross section from r→1 to r→2.(23)τr→1,r→2=∫r→1r→2dr→Σtr→With the purpose of generating results that correspond to the spatial discretization, we integrated the manufactured angular fluxes over each mesh.(24)ψm,Qi→MMS=∫Qi→dr→ψmMMSr→To construct the manufactured solution with scattering, the scattering source portion is subtracted from Qr→.(25)qQi→MMS=QQi→-Σs,Qi→∑m=1Mwmψm,Qi→MMSThe manufactured solution to the SN equation with scattering can now be stated as follows: the angular flux, ψm,Qi→MMS, is the solution to the SN equations when the fixed external source equals qQi→MMS, total cross section is Σt,Qi→, and the scattering cross section is Σs,Qi→. In our implementation, the integral over the spatial mesh, (24), is calculated numerically.Table 1
Total cross section and combined source of the test problem suites.
Case Total cross sectionΣt (cm−1) Combined sourceQ (cm−3s−1) Material 1 Material 2 Material 1 Material 2 Homogeneous 0.5 0.5 1 1 Heterogeneous 0.5 0.5×10−4 1 1×10−4Having 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.(26)θijk=log10εijkEijkwhere εijk is the estimated error and Eijk is the true error, both associated with cells  i,j,k. 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.(27)θ-1=1N∑θijk(28)θ-2=∑θijkEijk∑EijkThe θ-1 value gives a simple average of the absolute values of the effectivity index, whereas the θ-2 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 ratioc=0.0,0.1,0.3,0.5,0.7,0.9.Figure2 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 h/2 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.Figure 2
Histogram of effectivity for two-mesh estimator.
(a) (b)Figure3 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.Figure 3
Histogram of effectivity for residual-based estimator.
(a) (b)To globally compare the effectiveness of these estimators, Figure4 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 θ-1 and θ-2 are calculated using (27) and (28), respectively.Figure 4
Comparison of the averaged effectivity for the various zeroth-order spatial discretizations.
(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 h mesh and the h/2 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 Figure4, the true error weighted averaged effectivity index θ-2 is obviously smaller than θ-1 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 inx,y,z∈0,10cm. 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).Figure5 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.Figure 5
Ratio of the estimated error to the true error for the homogeneous case.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 Figure6. 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.Figure 6
Effects of the various adjoint spatial discretizations on the accuracy of the DWR estimator.
(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.



---

*Source: 3965309-2018-07-09.xml*