DEVELOPMENT AND TESTING OF THE MODIFIED SIMPSON’S RULE FOR DISCRETE ORDINATES TRANSPORT APPLICATIONS

A new angular quadrature type termed Modified Simpson Trapezoidal (MST) is developed based on the conventional Simpson’s 1/3 rule where the angular pattern over polar levels has a trapezoid shape. An adaptive coefficient correction scheme is developed to enable our new quadrature to integrate the angular flux over subintervals separated by the interior jump irregularities. A two-dimensional test problem is employed to verify the angular discretization error in the uncollided SN scalar flux computed with our new quadrature sets, as well as conventional angular quadrature types. Numerical results show that the MST quadrature error in the point-wise scalar flux converges with second order against increasing number of discrete angles, while the error obtained with other conventional quadrature types converges slower than first order depending on the regularity of the exact point-wise uncollided angular flux. In order to reduce the number of discrete points needed, a variant of the MST quadrature, namely MSTP30, is developed by using the Quadruple Range [1] polar quadrature with fixed 30 polar angles and applying the MST quadrature to the azimuthal dependence in each polar level. The angular discretization error in the point-wise SN scalar flux obtained with MSTP30 sets converges with fourth order because the polar discretization error is sufficiently reduced that MSTP30 behaves like a one-dimensional quadrature. Furthermore, because MSTP30 computes the integral over subintervals that keep the true solution’s irregularity at the boundaries, this fourth order convergence rate is unaffected by such inevitable irregularities.


INTRODUCTION
Accurate numerical solution of the steady-state monoenergetic particle transport equation in the deterministic scope is one of the key research topics in the nuclear computational science community. Both spatial and angular variables need to be discretized, and the applied discretization methods introduce discretization errors into the numerical solution. In this work, we focus on the angular discretization error caused by the widely-applied Discrete Ordinates (S N ) method in the numerical solution of the transport equation. Our previous work [2,3] considered the angular discretization error in the region-averaged scalar flux, and therefore here we switch our research objective to the angular discretization error in the point-wise scalar flux, assuming, as before, vanishing spatial discretization error. Under certain assumptions, the uncollided point-wise S N angular flux is exact, and thus, the angular discretization error in its scalar counterpart equals the angular quadrature error. Motivated by the conclusion from that work, we realized that it is the regularity of the exact angular flux that governs the angular quadrature error convergence rate. We have established that the exact uncollided point-wise angular flux's regularity with respect to the azimuthal angle is either C (1) Z (Ω ϕ ) if the boundary condition is continuous at the incoming corner, or C (0) Z (Ω ϕ ) if the boundary condition is discontinuous, where Ω ϕ = [0, 2π] is the range of azimuthal angle, and the jump singularities locate on the inside of each angular quadrant. Hence, the error of traditional angular quadrature types, e.g. Level Symmetric (LS) [4], Legendre-Chebyshev Quadrangular (LCQ) [5], Legendre-Chebyshev Triangular (LCT) [5], Quadruple Range (QR) [1] and Quadruple Range Spence type [6] (QRS), cannot converge rapidly. The reason is that the true angular flux usually is not sufficiently smooth within each angular quadrant, while these quadrature types' integration intervals are either the unit sphere or one quadrant. Based on the regularity of the exact angular flux and the conventional Simpson's rule, we developed a new quadrature type that can treat the interior singularities inside the integration interval to achieve faster error convergence rate than the traditional angular quadrature types.

DEVELOPMENT OF MODIFIED SIMPSON'S RULE
The conventional M -point one-dimensional Simpson's 1/3 rule is naturally considered since its error converges like O(M −4 ), such that a two-dimensional product Simpson's rule error is expected to converge with second order [7]. However, the conventional one-dimensional Simpson's rule requires the usage of the two endpoints of the integration interval [8]. This feature may cause an essential singularity in the transport solution because sin θ sin ϕ and sin θ cos ϕ approach zero for angles approaching quadrant boundaries, while they appear in the denominator of the characteristic path length that governs the attenuation term [2]. Also to achieve fourth order convergence, Simpson's rule requires the integrand to belong to C (r) Z [a, b], r ≥ 4, over the integration interval [a, b]. Simpson's rule error converges with less-than-fourth order if the integrand has a low regularity order caused by irregularities inside the integration interval. Therefore, we modified Simpson's rule to avoid these singular points while retaining the same error convergence order.

General Modified Simpson's Rule
An M -point standard Simpson's rule requires evenly-spaced discrete points over its integration interval [a, b] where the two endpoints are included in these M discrete points. If singularities occur at the endpoints, then the Standard Simpson's rule is not able to achieve the theoretical convergence rate even if the integrand over the integration interval is sufficiently smooth, and therefore we have to choose two alternative points for x 1 and x M .
Let us start with the simplest case of one subinterval. For a quadratic function where a, b, c 0 , c 1 , c 2 and c 2 are some constants. First we select three points within (a, b) and avoid the two endpoints: where h = (b − a)/2, α and β are constants. Then we build up a 'three-point system' in terms of h, and set the middle point to be the local origin for convenience: where γ and δ are functions of a, b, α and β, and it is easy to transform the integration interval from [a, b] to [−αh,βh]. The exact integration of f (x) over [−αh,βh] is Then we construct a linear system using points from Eq.(2) to solve for the corresponding coefficients: Thus, we have the modified Simpson's formula: Note that the traditional Simpson's rule is a special case with α = β = 0 andα =β = γ = δ = 1.
For the composite Simpson's rule, we apply Eq.(5) to every subinterval. Let f (x) be the integrand, x ∈ [a, b], and N be the number of subintervals in the integration interval, then the total number of discrete points is M = 2N + 1, and h = (b − a)/(M − 1). Assume a, b are points of essential singularities, and the interior singularity is not considered in this first step where the discrete points are evenly spaced, except the two endpoints: Next, for each subinterval i, i = 1, 2, · · · , N , we split the resulting discrete points into the threepoint systems, and also transform the corresponding integration interval to compute Simpson's coefficients w 1,i , w 2,i , w 3,i via Eq.(4). The composite modified Simpson's rule is given by Note that in this work, we assume the interior singularities are known in advance. These interior singularities can be determined from the configuration of problem, and they depend on both the spatial and angular variables. The method for identifying the solution's singularity will be reported in the future. Thus, to properly deal with the interior singularities, the Simpson coefficients in each subinterval should be determined on the fly. The coefficients of each fixed discrete ordinate changes according to the adaptive coefficient correction scheme that we develop in the following section. It is important to emphasize here that the added computational burden here is minimal because the set of discrete ordinates is fixed for all spatial points in the problem. To compute the scalar flux, and more generally angular moments of the flux, at a given point we first determine the locations of the internal azimuthal-angle irregularities, then perform partial integrals separately in each subinterval using the Modified Simpson Rule that avoids the end-points essential irregularity, then add the contributions for the full-quadrant integral as described in the next section.

Adaptive Coefficient Correction Scheme
For simplicity we assume there is only one singular point within the integration interval; the same method can be applied to cases with multiple singularities. Let ϕ s ∈ (a, b) be the only singularity of f (x) inside the integration interval, and x j < ϕ s < x j+1 , j = 1, 2, · · · , M , where M is the total number of discrete points. Then the integral can be split into two parts: For integral I L [f ], there are j discrete points located over (a, ϕ s ). If j is odd, and j ≥ 3, then the number of subintervals in (a, ϕ s ) is N L = (j−1)/2, we need to adjust Simpson's coefficients of the as depicted in Figure 1(a). Here again, Eq.(4) is applied to determine Simpson's coefficients once we represent ϕ s =β h. Note thatβ also depends on the number of points: As a special case, a flat approximation is employed when j = 1. If j is even, and j ≥ 4, then there is one isolated point that cannot define Simpson's rule for the corresponding subinterval. Let x j be the isolated point, and the (j − 1) point forms Simpson's rule with N L subintervals, where N L = (j −2)/2. To implement the modified Simpson's rule, x j−1 and x j−2 are reused to integrate the function over (x j−1 , ϕ s ). These three points form the (N L + 1)-th subinterval, as illustrated in Figure 1(b). We represent ϕ s =β h, and apply the same procedure again to obtain Simpson's coefficients for the three points. Finally adding w 1,N L +1 and w 2,N L +1 to w 2,N L and w 3,N L respectively we obtain the corrected Simpson's coefficients for x j−2 and x j−1 . As a special case, the trapezoidal rule is employed when j = 2.
The integral I R [f ] can be treated using analogous modifications to Simpson's rule, but x j+1 is chosen to be the isolated point if (M − j) is even. The one-dimensional modified Simpson's rule with the adaptive coefficient correction scheme This was achieved by creating subinterval interfaces where points of irregularity exist.

Product Modified Simpson's Rule
For solving the two-dimensional transport equation, the product Modified Simpson's rule is easy to obtain. First, we generate the N p -point polar quadrature over polar range [0, π/2]. Realizing that in 2D configurations there exists no interior singularities over the polar integration interval [0, π/2], we only need to treat the essential singularity occurring at θ = 0. Note that we restrict the number of polar angles N p per octant to be even for generating the azimuthal quadratures over each polar level as detailed in the next paragraph. The modified Simpson's rule is employed, and the discrete polar angles θ i , i = 1, 2, · · · , N p + 1, are determined by Eq.(6) with fixed parameters: a = 0, b = π/2, α = 1/2, β = 0, M = N p + 1. Note that the modified Simpson's rule gives N p + 1 discrete points. Since N p is even, we have to abandon the endpoint θ Np+1 = π/2. The adaptive coefficient correction scheme is applied to determine the corresponding polar quadrature weights w p i , i = 1, 2, · · · , N p .
For each polar level i, we set the number of discrete azimuthal angles N i a = N p , and thus the total number of the discrete angles per octant is M o = N 2 p . The procedures of generating the azimuthal quadrature is different from the generation of the polar quadrature. First, we equally divide the azimuthal range [0, π/2] into N 2 p subintervals, then all the azimuthal angles ϕ j , j = 1, 2, · · · , N 2 p , are determined using the midpoint of each resulting subinterval. Next, let all azimuthal angles in every polar level be symmetric with respect to π/4, then for polar level i, the azimuthal angles from (0, π/4) are selected as ϕ i k = ϕ kNp−i+1 , k = 1, 2, · · · , N p /2. The azimuthal angles from π/4 to π/2 are simply the mirror image of the azimuthal angles in (0, π/4). In this way, the azimuthal angle distribution has a trapezoidal shape, as illustrated by the example of N p = 4 in Figure 2. Once all discrete azimuthal angles are determined and fixed, the adaptive coefficient correction scheme is implemented to generate the modified Simpson's coefficients for each discrete azimuthal angles. The conventional product quadrature formula [7] is used to produce the two-dimensional product angular quadrature sets. Therefore we call this quadrature type the Modified Simpson's rule with the adaptive coefficient correction scheme and a Trapezoid shape of azimuthal angles (MST-N ). Note that N in this notation indicates the number of the polar angles per octant N p .
Finally, we develop another version of the MST quadrature sets to improve the error's convergence order by replacing the MST polar quadrature with the QR polar quadrature fixed at 30 discrete polar angles per octant. The QR polar quadrature is a Gaussian Christoffel quadrature, therefore the quadrature error converges rapidly as long as the integrand is smooth enough. Note that the point-wise angular flux has infinite continuous weak derivatives with respect to the polar angle for two-dimensional problems. We have verified that if the azimuthal angle is fixed, then the polar quadrature error is of order 10 −16 when using QR polar quadrature with 30 polar angles per octant. That's why we employ the 30-point QR polar quadrature to replace the MST polar quadrature sets.
Since for this new quadrature type, the number of the polar angles is fixed, the number of azimuthal angles per polar level and per octant is set to N and therefore the total number of angles per octant is M o = 30N , and we follow the same procedure of generation of MST azimuthal quadrature sets to obtain the azimuthal angles per polar level. Similarly, the adaptive coefficient correction scheme is implemented to generate the modified Simpson's coefficient for each discrete azimuthal angle, and again, the conventional product quadrature formula is used to produce the two-dimensional product angular quadrature sets. We name this quadrature type as MSTP30-N , where N represents the number of azimuthal angles per polar level and per octant.

NUMERICAL RESULTS
A two-dimensional problem is employed to verify the angular discretization error obtained by our new modified Simpson's rule. The configuration of this problem is depicted in Figure 3. This is a rectangular region (x ∈   First, we derive the analytical form of the exact point-wise uncollided angular flux at an arbitrary point by using the integral form of the transport equation, and we integrate (using MATLAB) the resulting exact angular flux with respect to angle over the unit sphere to obtain the exact scalar flux, then we apply the quadrature formula to the exact angular flux to obtain the S N scalar flux with no spatial discretization error. The angular discretization error is determined by subtracting the latter from the former for each quadrature type and number of angles per octant, M o . Figure 4(a) shows the absolute value of the angular discretization error in the point-wise uncollided S N scalar flux against increasing number of discrete angles per octant for Case 1. We observe that the scalar flux errors obtained with the LC and QR class quadrature sets converge with first order, and the error computed with MST sets converges with second order. These error trends are consistent with the comprehensive quadrature error theory reported in Ref. [2]. In this case, the regularity of the exact solution with respect to the azimuthal angle is C where Ω q ϕ = [(j − 1)π/2, jπ/2], j = 1, 2, 3, 4, is the azimuthal angle range per quadrant, thus explaining the LC and QR class' first order convergence. In contrast, for the MST quadrature sets, although their designed integration interval is one angular quadrant, the adaptive coefficient correction scheme allows them to treat the interior singularities exactly. We verified that the exact point-wise uncollided angular flux is analytic with respect to the azimuthal angle in ((j − 1)π/2, ϕ s ) ∪ (ϕ s , jπ/2). Thus, the product MST quadrature error converges with second order. Compared with the MST sets, the MSTP30 quadrature error converges with fourth order. This phenomenon can be explained as follows. Since we fixed the number of discrete polar angles, the increasing number of discrete angles in MSTP30 sets is only due to increasing the number of azimuthal angles, and thus it can be considered as a special one-dimensional quadrature set. Also, the polar quadrature error in this case is of order 10 −16 , essentially numerical zero for double-precision computations, and far smaller than the error of the corresponding azimuthal quadrature error. Hence, if the error is dominated by the azimuthal quadrature, then the MSTP30 sets yield the standard fourth order convergence rate expected for the one-dimensional quadrature sets.  , while the MST quadrature error still converges with second order. In this case, the boundary condition is discontinuous, so the exact angular flux is discontinuous inside the angular quadrant, and therefore its regularity with respect to the azimuthal angle is C (0) Z (Ω q ϕ ), which leads to the order one-half convergence of the quadrature error for the LC and QR class. Here again, the interior singularities cannot retard the error convergence rate of the MST sets due to the adaptive coefficient correction scheme. Note that irregular oscillation occurs in the MST quadrature error when the number of discrete angles is greater than 10 4 . The reason behind this phenomenon remains an open question. The MSTP30 quadrature error, here again, converges with fourth order.  ). The azimuthal quadrature employed in our LC quadrature sets is the Gauss-Chebyshev quadrature (first kind) that is equally spaced and weighted, and the quadrature points do not include points of angular irregularity due to the discontinuous dependence of the total cross section on the azimuthal angle at this point, namely at ϕ = nπ/2, n = 0, · · · , 4. As a result, the true solution can be considered continuous in C It is worthing mentioning that the MSTP30 error will converge to some non-zero limit as the number of azimuthal angles keeps increasing. Because for any finite N p the fixed polar quadrature error is not perfectly zero, once the azimuthal quadrature error in one polar level is smaller than the polar quadrature error, the total MSTP30 error observed in Figure 4 will become equal to N a times polar quadrature errors, where N a is the number of azimuthal angles. This indicates that the total MSTP30 error will increase beyond a critical value of M o where the azimuthal quadrature error is smaller than the polar quadrature error.

CONCLUSIONS
We conclude that our angular quadrature rule based on the traditional Simpson's rule is able to avoid the essential singularity at the endpoints of the angular quadrant, and the adaptive coefficient correction scheme is able to directly treat the interior singularities. Thus the quadrature error convergence rate is only dominated by the dimensionality and the intrinsic fourth order convergence of the one-dimensional quadrature error. Numerical results of the point-wise uncollided scalar flux show that our new quadrature error converges faster than conventional angular quadrature sets, unaffected by various solution irregularities.