Adaptive discontinuous finite element quadrature sets over an icosahedron for discrete ordinates method

Adaptive discontinuous finite element quadrature sets over an icosahedron for discrete ordinates method 1 first-author https://orcid.org/0000-0002-0347-0791 1 corresp Bin Zhang zhangbin@ncepu.edu.cn Bin Zhang zhangbin@ncepu.edu.cn zhangbin@ncepu.edu.cn https://orcid.org/0000-0001-5351-3860 1 1 School of Nuclear Science and Engineering, North China Electric Power University, Beijing 102206, China School of Nuclear Science and Engineering, North China Electric Power University, Beijing 102206, China The discrete ordinates (S N ) method requires numerous angular unknowns to achieve the desired accuracy for shielding calculations involving strong anisotropy. Our objective is to develop an angular adaptive algorithm in the S N method to automatically optimize the angular distribution and minimize angular discretization errors with lower expenses. The proposed method enables linear discontinuous finite element quadrature sets over an icosahedron to vary their quadrature orders in a one-twentieth sphere so that fine resolutions can be applied to the angular domains that are important. An error estimation that operates in conjunction with the spherical harmonics method is developed to determine the locations where more refinement is required. The adaptive quadrature sets are applied to three duct problems, including the Kobayashi benchmarks and IRI-TUB research reactor, which emphasize the ability of this method to resolve neutron streaming through ducts with voids. The results indicate that the performance of the adaptive method is more efficient than that of uniform quadrature sets for duct transport problems. Our adaptive method offers an appropriate placement of angular unknowns to accurately integrate angular fluxes while reducing the computational costs in terms of unknowns and run times. Shielding calculation Discrete ordinates method Angular adaptivity Discontinuous finite element 1. Introduction 1 High-efficiency shielding calculation centered on the solution of the neutron transport equation is a concern in nuclear computational science. The variables of a steady Boltzmann transport equation include energy, space, and direction of motion, forming a six-dimensional phase space [ 1 1 ] . The discrete ordinates (S N ) method approximates this equation by discretization into a set of algebraic equations coupled only through the scattering operator. The angular fluxes along a finite number of discrete ordinate unit vectors are evaluated in the S N method, and the angular integral is replaced using quadrature sets with weights associated with these ordinates. However, angular discretization errors introduced through quadrature sets in the S N method, in a way, limit its use in many problems, especially ones containing ducts with voids. The duct problem has a highly anisotropic flux throughout the duct regions and features a heavy-streaming region. Angular discretization often requires varying resolutions over the angular domain of the problem under consideration. Uniformly increasing the quadrature directions results in an increase in the size of the calculation, and thus is not the most efficient way to reduce this type of discretization error. This issue has led to the development of local refinement with quadrature sets to optimize the angular distribution and accurately integrate with minimal effort. Many local angular refinement quadrature sets have been studied over the last several decades and applied to many practical problems. The biased quadrature set, first proposed by Jenal and Erickson in 1977, starts with a semi-symmetric Gaussian quadrature set and employs an upward or downward refinement in a hemisphere [ 2 2 ] . In 2001, Longoni developed the longitudinal segmentation (OS) [ 3 3 ] and regional angular refinement (RAR) [ 4 4 ] to add discrete directions only in the focus angular regions for different problems. However, refinements in the quadrature sets mentioned above rely on user intuition or a prior knowledge of certain issues. This technique is unlikely to obtain the optimum distribution and minimize the discretization errors. Therefore, an automated refinement referred to as the angular adaptive method has recently flourished in the S N method. Stone first studied the adaptive method within the S N method in 2007. The focus of his work is based on traditional quadrature sets to automate refinement in different quadrature regions [ 5 5 ] .The adaptivity is restricted to two-dimensional Cartesian geometry and isotropic scattering. In 2010, Jarrell developed regular angular adaptivity and mapping algorithms using linear discontinuous finite-element (LDFE) quadrature sets [ 6 6 ] . However, the non-uniformity in the directions and weights of the LDFE quadrature sets over an octahedron may increase the local integration errors. As an extension of Jarrell’s research, Lau developed linear and quadratic discontinuous finite-element quadrature sets over spherical quadrilaterals and employed them in the corresponding adaptive algorithms [ 7 7 ] . They also studied an optimal mapping algorithm that preserves both the shape and angular moments of interest. In 2018, with the spherical quadrilateral quadrature sets, Zhang evaluated a goal-oriented angular adaptive algorithm based on the adjoint value theory to mitigate ray effects [ 8 8 ] . However, the accuracy of the applied quadrature sets may hinder the development of this angular adaptive method. Linear discontinuous finite element quadrature sets over an icosahedron (ICLDFE quadrature sets) have been developed in previous research [ 9 9 ] . This quadrature set features uniform placement of points and positive weights, local refinement, and high local accuracy for the integration of spherical harmonics with the 4th order convergence. This study aims to evaluate an angular adaptive algorithm built on ICLDFE quadrature sets to produce efficient and accurate solutions for shielding problems with ducts. The remainder of this paper is organized as follows. In Section 2, we briefly introduce the principles of the proposed angular adaptive algorithm. Section 3 presents the numerical results for three benchmark problems. Finally, the conclusions are summarized in Section 4. 2. Angular Adaptivity Algorithm 1 2.1. Adaptive Quadrature Sets 2 The adaptive quadrature sets applied in this study are the ICLDFE quadrature sets. Their construction is shown in Fig. 1 Fig. 1 . The angular domain is first divided into spherical triangles by inscribing an icosahedron into the unit sphere, as shown in Fig. 2 Fig. 2 . Generally, the centroid of the triangles on the icosahedron is chosen as the initial direction and then projected from the regular polyhedron to the unit sphere. The ICLDFE quadrature sets start with a base triangle and add quadrature points by refining it to sub-triangles. The discrete directions in a spherical triangle have an exponential growth by 4 N ( N is the order of the quadrature sets). Linear discontinuous finite element basis functions in the directions are established over spherical triangles. The weights <math id="M1"><mrow><msub><mi>w</mi><mi>i</mi></msub></mrow></math> are then calculated by integrating the underlying basis functions in the corresponding spherical triangles. The basis functions <math id="M2"><mrow><mi>b</mi><mo stretchy="false">(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>Ω</mi></mstyle><mo stretchy="false">)</mo></mrow></math> are represented by a spherical harmonic expansion of the zeroth and first order at each quadrature point denoted as i, where <math id="M3"><mrow><msub><mi>a</mi><mrow><mn>0</mn><mo>,</mo><mi>i</mi></mrow></msub></mrow></math> , <math id="M4"><mrow><msub><mi>a</mi><mrow><mn>1</mn><mo>,</mo><mi>i</mi></mrow></msub></mrow></math> , <math id="M5"><mrow><msub><mi>a</mi><mrow><mn>2</mn><mo>,</mo><mi>i</mi></mrow></msub></mrow></math> and <math id="M6"><mrow><msub><mi>a</mi><mrow><mn>3</mn><mo>,</mo><mi>i</mi></mrow></msub></mrow></math> are unknown coefficients. The angular vector <math id="M7"><mstyle mathvariant="bold" mathsize="normal"><mi>Ω</mi></mstyle></math> is described by the variables <math id="M8"><mi>μ</mi></math> , <math id="M9"><mi>η</mi></math> , and <math id="M10"><mi>ξ</mi></math> , and calculated using the polar angle and azimuthal angle. These coefficients are solved by the property that the basis functions have a value of unity at the corresponding discrete direction but a value of zero for the other three points. <math display="block" id="M11"><mrow><msub><mi>b</mi><mi>i</mi></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>Ω</mi></mstyle><mo stretchy="false">)</mo><mo>=</mo><msub><mi>a</mi><mrow><mn>0</mn><mo>,</mo><mi>i</mi></mrow></msub><mo>+</mo><msub><mi>a</mi><mrow><mn>1</mn><mo>,</mo><mi>i</mi></mrow></msub><mi>μ</mi><mo>+</mo><msub><mi>a</mi><mrow><mn>2</mn><mo>,</mo><mi>i</mi></mrow></msub><mi>η</mi><mo>+</mo><msub><mi>a</mi><mrow><mn>3</mn><mo>,</mo><mi>i</mi></mrow></msub><mi>ξ</mi><mo>,</mo><mi>i</mi><mo>=</mo><mn>1</mn><mo>,</mo><mn>2</mn><mo>,</mo><mn>3</mn><mo>,</mo><mn>4</mn></mrow></math> <math display="block" id="M12"><mrow><msub><mi>w</mi><mi>i</mi></msub><mo>=</mo><mstyle displaystyle="true"><mrow><msub><mo>∫</mo><mi>Δ</mi></msub><mrow><mi>b</mi><mo stretchy="false">(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>Ω</mi></mstyle><mo stretchy="false">)</mo><mtext>d</mtext><mstyle mathvariant="bold" mathsize="normal"><mi>Ω</mi></mstyle></mrow></mrow></mstyle></mrow></math> The basis functions are not completely positive in the corresponding integral regions, which leads to meaningless negative weights as refinement continues. To overcome this, it is necessary to improve the arrangement of the quadrature directions. To avoid duplication, the derivations of the optimization methods have been simplified here, and previous research [ 9 9 ] demonstrated the details. The purpose of this method is to force the weights of the center points equal to their associated sub-region surface areas to guarantee positive weights. The ICLDFE quadrature sets are generated only on one surface of the icosahedron and repeated by rotations on the other 19 surfaces. The number of quadrature points grows exponentially by 20×4 N (i.e., in the case of ICLDFE-S 2 , the number of directions is 20×4 2 =320). They are invariant under 72 degree rotations about the three axes joining the extreme opposite vertices and inversion about the origin. However, it must be noted that reflection in S N calculations requires that the quadrature sets possess 90 degree symmetry and thus ICLDFE quadrature sets can only be available in problems involving vacuum conditions. ICLDFE quadrature sets have many desirable properties that are suitable for adaptive algorithms. They are amenable to local refinement, have strictly positive weights, and have higher local integration accuracy. Our quadrature sets are designed to integrate the zeroth- and first-order spherical harmonic functions exactly in a one-twentieth angular region. Therefore, they integrate non-smooth functions, such as peaky and nearly discontinuous angular flux with higher accuracy, which are commonly encountered in practical transport problems, compared to traditional quadrature sets. 2.2. Adaptive Strategies 2 The common adaptive strategy applies a progressive rule by starting uniformly distributed points and successively increasing the number of quadrature points in a trouble region until the desired accuracy is achieved. Our angular adaptive algorithm with the ICLDFE quadrature sets requires a driving force provided by posterior error estimations to determine the locations that require further refinement. The relative error between the calculated results and interpolations in the testing directions is used as the error estimation in this angular adaptivity. <math display="block" id="M13"><mrow><mi>Δ</mi><msub><mi>ϕ</mi><mrow><mtext>test</mtext></mrow></msub><mo>=</mo><mfrac><mrow><mrow><mo>|</mo><mrow><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>M</mi></munderover><mrow><msub><mi>w</mi><mi>i</mi></msub><msubsup><mi>ψ</mi><mi>i</mi><mrow><mtext>cal</mtext></mrow></msubsup><mo>−</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>M</mi></munderover><mrow><msub><mi>w</mi><mi>i</mi></msub><msubsup><mi>ψ</mi><mi>i</mi><mrow><mi>int</mi></mrow></msubsup></mrow></mstyle></mrow></mstyle></mrow><mo>|</mo></mrow></mrow><mrow><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>M</mi></munderover><mrow><msub><mi>w</mi><mi>i</mi></msub><msubsup><mi>ψ</mi><mi>i</mi><mrow><mtext>cal</mtext></mrow></msubsup></mrow></mstyle></mrow></mfrac></mrow></math> where the testing directions i =1… M are the directions of the higher-order quadrature sets in one spherical triangle. <math id="M14"><mrow><mi>Δ</mi><msub><mi>ϕ</mi><mrow><mtext>test</mtext></mrow></msub></mrow></math> is the maximum relative error between the calculated angular fluxes and interpolated results of all the spatial boundary grids. <math id="M15"><mrow><msubsup><mi>ψ</mi><mi>i</mi><mrow><mtext>cal</mtext></mrow></msubsup></mrow></math> denotes the angular flux through the boundary surfaces calculated by a transport sweep in the testing directions using the latest scattering and fixed source as well as the boundary conditions. <math id="M16"><mrow><msubsup><mi>ψ</mi><mi>i</mi><mrow><mi>int</mi></mrow></msubsup></mrow></math> represents the interpolated angular flux at the testing ordinates determined by the spherical harmonic function. The spherical harmonic function method starts by expanding the angular flux into the spherical harmonic function of the third order and below as follows. <math display="block" id="M17"><mtable><mtr><mtd><mi>ψ</mi><mo stretchy="false">(</mo><mover accent="true"><mi>Ω</mi><mo>→</mo></mover><mo stretchy="false">)</mo><mo>=</mo><msub><mi>c</mi><mn>0</mn></msub><mo>+</mo><msub><mi>c</mi><mn>1</mn></msub><mi>μ</mi><mo>+</mo><msub><mi>c</mi><mn>2</mn></msub><mi>η</mi><mo>+</mo><msub><mi>c</mi><mn>3</mn></msub><mi>ξ</mi><mo>+</mo><msub><mi>c</mi><mn>4</mn></msub><mfrac><mn>1</mn><mn>2</mn></mfrac><mo stretchy="false">(</mo><mn>3</mn><msup><mi>ξ</mi><mn>2</mn></msup><mo>−</mo><mn>1</mn><mo stretchy="false">)</mo><mo>+</mo><msub><mi>c</mi><mn>5</mn></msub><msqrt><mn>3</mn></msqrt><mi>μ</mi><mi>ξ</mi><mo>+</mo><msub><mi>c</mi><mn>6</mn></msub><msqrt><mn>3</mn></msqrt><mi>η</mi><mi>ξ</mi><mo>+</mo><msub><mi>c</mi><mn>7</mn></msub><msqrt><mn>3</mn></msqrt><mi>μ</mi><mi>η</mi></mtd></mtr><mtr><mtd><mo>+</mo><msub><mi>c</mi><mn>8</mn></msub><mfrac><mrow><msqrt><mn>3</mn></msqrt></mrow><mn>2</mn></mfrac><mo stretchy="false">(</mo><msup><mi>μ</mi><mn>2</mn></msup><mo>−</mo><msup><mi>η</mi><mn>2</mn></msup><mo stretchy="false">)</mo><mo>+</mo><msub><mi>c</mi><mn>9</mn></msub><mfrac><mn>1</mn><mn>2</mn></mfrac><mo stretchy="false">(</mo><mn>5</mn><msup><mi>ξ</mi><mn>3</mn></msup><mo>−</mo><mn>3</mn><mi>ξ</mi><mo stretchy="false">)</mo><mo>+</mo><msub><mi>c</mi><mrow><mn>10</mn></mrow></msub><mfrac><mrow><msqrt><mn>6</mn></msqrt></mrow><mn>4</mn></mfrac><mi>μ</mi><mo stretchy="false">(</mo><mn>5</mn><msup><mi>ξ</mi><mn>2</mn></msup><mo>−</mo><mn>1</mn><mo stretchy="false">)</mo><mo>+</mo><msub><mi>c</mi><mrow><mn>11</mn></mrow></msub><mfrac><mrow><msqrt><mn>6</mn></msqrt></mrow><mn>4</mn></mfrac><mi>η</mi><mo stretchy="false">(</mo><mn>5</mn><msup><mi>ξ</mi><mn>2</mn></msup><mo>−</mo><mn>1</mn><mo stretchy="false">)</mo></mtd></mtr><mtr><mtd><mo>+</mo><msub><mi>c</mi><mrow><mn>12</mn></mrow></msub><mfrac><mrow><msqrt><mrow><mn>15</mn></mrow></msqrt></mrow><mn>2</mn></mfrac><mi>ξ</mi><mo stretchy="false">(</mo><msup><mi>μ</mi><mn>2</mn></msup><mo>−</mo><msup><mi>η</mi><mn>2</mn></msup><mo stretchy="false">)</mo><mo>+</mo><msub><mi>c</mi><mrow><mn>13</mn></mrow></msub><msqrt><mrow><mn>15</mn></mrow></msqrt><mi>μ</mi><mi>η</mi><mi>ξ</mi><mo>+</mo><msub><mi>c</mi><mrow><mn>14</mn></mrow></msub><msqrt><mrow><mfrac><mn>5</mn><mn>8</mn></mfrac></mrow></msqrt><mo stretchy="false">(</mo><msup><mi>μ</mi><mn>3</mn></msup><mo>−</mo><mn>3</mn><mi>μ</mi><msup><mi>η</mi><mn>2</mn></msup><mo stretchy="false">)</mo><mo>+</mo><msub><mi>c</mi><mrow><mn>15</mn></mrow></msub><msqrt><mrow><mfrac><mn>5</mn><mn>8</mn></mfrac></mrow></msqrt><mo stretchy="false">(</mo><mn>3</mn><msup><mi>μ</mi><mn>2</mn></msup><mi>η</mi><mo>−</mo><msup><mi>η</mi><mn>3</mn></msup><mo stretchy="false">)</mo></mtd></mtr></mtable></math> where c 0, c 1, …, c 15 are the unknown coefficients to be sought. Equation (4) can be written as a matrix <math display="block" id="M18"><mrow><mrow><mo>[</mo><mrow><mtable><mtr><mtd><mn>1</mn></mtd><mtd><mrow><msub><mi>μ</mi><mn>1</mn></msub></mrow></mtd><mtd><mo>⋯</mo></mtd><mtd><mrow><msqrt><mrow><mfrac><mn>5</mn><mn>8</mn></mfrac></mrow></msqrt><mo stretchy="false">(</mo><mn>3</mn><msub><mi>μ</mi><mn>1</mn></msub><msup><mrow/><mn>2</mn></msup><msub><mi>η</mi><mn>1</mn></msub><mo>−</mo><msub><mi>η</mi><mn>1</mn></msub><msup><mrow/><mn>3</mn></msup><mo stretchy="false">)</mo></mrow></mtd></mtr><mtr><mtd><mn>1</mn></mtd><mtd><mrow><msub><mi>μ</mi><mn>2</mn></msub></mrow></mtd><mtd><mo>…</mo></mtd><mtd><mrow><msqrt><mrow><mfrac><mn>5</mn><mn>8</mn></mfrac></mrow></msqrt><mo stretchy="false">(</mo><mn>3</mn><msub><mi>μ</mi><mn>2</mn></msub><msup><mrow/><mn>2</mn></msup><msub><mi>η</mi><mn>2</mn></msub><mo>−</mo><msub><mi>η</mi><mn>2</mn></msub><msup><mrow/><mn>3</mn></msup><mo stretchy="false">)</mo></mrow></mtd></mtr><mtr><mtd><mo>⋮</mo></mtd><mtd><mo>⋮</mo></mtd><mtd><mo>⋱</mo></mtd><mtd><mo>⋮</mo></mtd></mtr><mtr><mtd><mn>1</mn></mtd><mtd><mrow><msub><mi>μ</mi><mrow><msup><mn>4</mn><mi>n</mi></msup></mrow></msub></mrow></mtd><mtd><mo>…</mo></mtd><mtd><mrow><msqrt><mrow><mfrac><mn>5</mn><mn>8</mn></mfrac></mrow></msqrt><mo stretchy="false">(</mo><mn>3</mn><msub><mi>μ</mi><mrow><msup><mn>4</mn><mi>n</mi></msup></mrow></msub><msup><mrow/><mn>2</mn></msup><msub><mi>η</mi><mrow><msup><mn>4</mn><mi>n</mi></msup></mrow></msub><mo>−</mo><msub><mi>η</mi><mrow><msup><mn>4</mn><mi>n</mi></msup></mrow></msub><msup><mrow/><mn>3</mn></msup><mo stretchy="false">)</mo></mrow></mtd></mtr></mtable></mrow><mo>]</mo></mrow><mo>⋅</mo><mrow><mo>[</mo><mrow><mtable><mtr><mtd><mrow><msub><mi>c</mi><mn>0</mn></msub></mrow></mtd></mtr><mtr><mtd><mrow><msub><mi>c</mi><mn>1</mn></msub></mrow></mtd></mtr><mtr><mtd><mo>⋮</mo></mtd></mtr><mtr><mtd><mrow><msub><mi>c</mi><mrow><mn>15</mn></mrow></msub></mrow></mtd></mtr></mtable></mrow><mo>]</mo></mrow><mo>=</mo><mrow><mo>[</mo><mrow><mtable><mtr><mtd><mrow><msub><mi>ψ</mi><mn>1</mn></msub></mrow></mtd></mtr><mtr><mtd><mrow><msub><mi>ψ</mi><mn>2</mn></msub></mrow></mtd></mtr><mtr><mtd><mo>⋮</mo></mtd></mtr><mtr><mtd><mrow><msub><mi>ψ</mi><mrow><msup><mn>4</mn><mi>n</mi></msup></mrow></msub></mrow></mtd></mtr></mtable></mrow><mo>]</mo></mrow></mrow></math> where <math id="M19"><mrow><mrow><mo>{</mo><mrow><msub><mi>μ</mi><mi>i</mi></msub><mo>,</mo><msub><mi>η</mi><mi>i</mi></msub><mo>,</mo><msub><mi>ξ</mi><mi>i</mi></msub></mrow><mo>}</mo></mrow><mo>,</mo><mi>i</mi><mo>=</mo><mn>1</mn><mo>,</mo><mn>2</mn><mo>,</mo><msup><mrow><mn>...4</mn></mrow><mi>n</mi></msup></mrow></math> are the discrete directions of the current quadrature sets. <math id="M20"><mrow><msub><mi>ψ</mi><mn>0</mn></msub><mo>,</mo><msub><mi>ψ</mi><mn>1</mn></msub><mo>,</mo><mn>...</mn><mo>,</mo><msub><mi>ψ</mi><mrow><msup><mn>4</mn><mi>n</mi></msup></mrow></msub></mrow></math> are the known angular fluxes of the current directions, and n is the quadrature order. A simple inverse is needed to solve these coefficients in the spherical harmonic expansion when the initial quadrature sets default ICLDFE-S 2 quadrature sets (16 directions per angular region). For cases where the quadrature order is more than two, Eq. (5) becomes an overdetermined system solved by the least squares method. The interpolated angular fluxes at the testing directions are then calculated using <math display="block" id="M21"><mrow><mrow><mo>[</mo><mrow><mtable><mtr><mtd><mrow><msubsup><mi>ψ</mi><mn>0</mn><mrow><mi>int</mi></mrow></msubsup></mrow></mtd></mtr><mtr><mtd><mrow><msubsup><mi>ψ</mi><mn>1</mn><mrow><mi>int</mi></mrow></msubsup></mrow></mtd></mtr><mtr><mtd><mo>⋮</mo></mtd></mtr><mtr><mtd><mrow><msubsup><mi>ψ</mi><mrow><msup><mn>4</mn><mrow><mi>n</mi><mo>+</mo><mn>1</mn></mrow></msup></mrow><mrow><mi>int</mi></mrow></msubsup></mrow></mtd></mtr></mtable></mrow><mo>]</mo></mrow><mo>=</mo><mrow><mo>[</mo><mrow><mtable><mtr><mtd><mn>1</mn></mtd><mtd><mrow><msub><mi>μ</mi><mn>1</mn></msub></mrow></mtd><mtd><mo>⋯</mo></mtd><mtd><mrow><msqrt><mrow><mfrac><mn>5</mn><mn>8</mn></mfrac></mrow></msqrt><mo stretchy="false">(</mo><mn>3</mn><msub><mi>μ</mi><mn>1</mn></msub><msup><mrow/><mn>2</mn></msup><msub><mi>η</mi><mn>1</mn></msub><mo>−</mo><msub><mi>η</mi><mn>1</mn></msub><msup><mrow/><mn>3</mn></msup><mo stretchy="false">)</mo></mrow></mtd></mtr><mtr><mtd><mn>1</mn></mtd><mtd><mrow><msub><mi>μ</mi><mn>2</mn></msub></mrow></mtd><mtd><mo>⋯</mo></mtd><mtd><mrow><msqrt><mrow><mfrac><mn>5</mn><mn>8</mn></mfrac></mrow></msqrt><mo stretchy="false">(</mo><mn>3</mn><msub><mi>μ</mi><mn>2</mn></msub><msup><mrow/><mn>2</mn></msup><msub><mi>η</mi><mn>2</mn></msub><mo>−</mo><msub><mi>η</mi><mn>2</mn></msub><msup><mrow/><mn>3</mn></msup><mo stretchy="false">)</mo></mrow></mtd></mtr><mtr><mtd><mo>⋮</mo></mtd><mtd><mo>⋮</mo></mtd><mtd><mo>⋱</mo></mtd><mtd><mo>⋮</mo></mtd></mtr><mtr><mtd><mn>1</mn></mtd><mtd><mrow><msub><mi>μ</mi><mrow><msup><mn>4</mn><mrow><mi>n</mi><mo>+</mo><mn>1</mn></mrow></msup></mrow></msub></mrow></mtd><mtd><mo>⋯</mo></mtd><mtd><mrow><msqrt><mrow><mfrac><mn>5</mn><mn>8</mn></mfrac></mrow></msqrt><mo stretchy="false">(</mo><mn>3</mn><msub><mi>μ</mi><mrow><msup><mn>4</mn><mrow><mi>n</mi><mo>+</mo><mn>1</mn></mrow></msup></mrow></msub><msup><mrow/><mn>2</mn></msup><msub><mi>η</mi><mrow><msup><mn>4</mn><mrow><mi>n</mi><mo>+</mo><mn>1</mn></mrow></msup></mrow></msub><mo>−</mo><msub><mi>η</mi><mrow><msup><mn>4</mn><mrow><mi>n</mi><mo>+</mo><mn>1</mn></mrow></msup></mrow></msub><msup><mrow/><mn>3</mn></msup><mo stretchy="false">)</mo></mrow></mtd></mtr></mtable></mrow><mo>]</mo></mrow><mo>⋅</mo><mrow><mo>[</mo><mrow><mtable><mtr><mtd><mrow><msub><mi>c</mi><mn>0</mn></msub></mrow></mtd></mtr><mtr><mtd><mrow><msub><mi>c</mi><mn>1</mn></msub></mrow></mtd></mtr><mtr><mtd><mo>⋮</mo></mtd></mtr><mtr><mtd><mrow><msub><mi>c</mi><mrow><mn>15</mn></mrow></msub></mrow></mtd></mtr></mtable></mrow><mo>]</mo></mrow></mrow></math> where <math id="M22"><mrow><mrow><mo>{</mo><mrow><msub><mi>μ</mi><mi>i</mi></msub><mo>,</mo><msub><mi>η</mi><mi>i</mi></msub><mo>,</mo><msub><mi>ξ</mi><mi>i</mi></msub></mrow><mo>}</mo></mrow><mo>,</mo><mi>i</mi><mo>=</mo><mn>1</mn><mo>,</mo><mn>2</mn><mo>,</mo><msup><mrow><mn>...4</mn></mrow><mrow><mi>n</mi><mo>+</mo><mn>1</mn></mrow></msup></mrow></math> are the directions of the testing quadrature sets, and <math id="M23"><mrow><msubsup><mi>ψ</mi><mrow><msub><mrow/><mn>1</mn></msub></mrow><mrow><mi>int</mi></mrow></msubsup><mo>,</mo><msubsup><mi>ψ</mi><mrow><msub><mrow/><mn>1</mn></msub></mrow><mrow><mi>int</mi></mrow></msubsup><mo>,</mo><mn>...</mn><mo>,</mo><msubsup><mi>ψ</mi><mrow><msup><mn>4</mn><mrow><mi>n</mi><mo>+</mo><mn>1</mn></mrow></msup></mrow><mrow><mi>int</mi></mrow></msubsup></mrow></math> are the angular fluxes through boundary surfaces at testing directions calculated by interpolation. A user-defined flux error tolerance <math id="M24"><mi>ε</mi></math> is used to determine if refinement is necessary. To refine, the flux difference in the testing directions must be less than the tolerance <math id="M25"><mi>ε</mi></math> on any meshes on the spatial boundary. 2.3. Implementation of Angular Adaptivity 2 The angular adaptivity was implemented in a three-dimensional particle transport code named ARES [ 10 10 ] . The implementation of the angular adaptivity in the transport calculations is illustrated in Fig. 3 Fig. 3 . The source iteration is performed to solve the discrete equations in this study. Each inner iteration is generally composed of normal transport sweeps and angular adaptive sweeps. Some user-defined parameters, such as the initial quadrature set order, number of sweeps before the first adaptive sweep ( N a ), the number of sweeps between two adaptive sweeps, and flux error tolerance, should be defined in advance. At least one normal transport sweep must be completed before any refinement can commence, and the users should determine how often an adaptive sweep is performed. After every user-defined number of adaptive sweeps, the angular fluxes are tested at the spatial boundary to determine whether the angular domain needs to be refined. The angular flux in the first group tends to be the most unsmooth and peaky. As neutrons slow down, the angular distribution becomes isotropic, and the need for integrated accuracy for the quadrature sets is relaxed. Hence, angular adaptivity is only implemented in the first energy group for multi-group shielding problems. An appropriate angular arrangement is then generated and applied to other energy groups for the subsequent calculations. The angular adaptive sweep can be described in detail as follows. All 20 angular regions adopt the same low-order quadrature sets in advance (taking ICLDFE-S 2 as an example). The ICLDFE-S 3 quadrature sets are selected as the testing directions, and the calculated angular fluxes of these directions through the spatial boundary are obtained by transport calculations using the latest scattering and fixed source as well as boundary conditions. The spherical harmonic function method is used to calculate the interpolated angular fluxes in the test directions. The differences between the transport calculations and interpolated results act as error estimations. When refinement is necessary, the quadrature set list adds the testing directions and removes the old directions. It completes an angular adaptive sweep until all the 20 regions are finished, and then the quadrature sets should be updated. This algorithm can acquire different resolutions in every one-twentieth angular region throughout the spatial domain. Since the construction of quadrature sets involves an optimized process and this convergence of optimization tends to be difficult as the quadrature order increases, ICLDFE quadrature sets of S 1 to S 6 have finished calculation and stored by the data structure of “triangle first and order second” in advance to avoid excessive expenses during adaptivity. The adaptive algorithm by necessity uses more computer time per unknown than the uniform quadrature sets because the investment is made to compare the calculated angular fluxes and interpolated values. However, the adaptive process can result in significant computational savings to mitigate this hindrance. Although this adaptivity has the ability to produce accurate results with lower computational cost in terms of unknowns and calculation time, there are some limitations in some areas. A potential downside of this method is that it introduces the angular arrangement throughout the problem domain instead of in different regions. In many practical problems, the angular domain containing the majority of angular fluxes may not be the same for different spatial regions. Additionally, this algorithm is only applicable to problems with full vacuum conditions. In fact, these limitations are caused by the lack of 90 degree symmetry for the ICLDFE quadrature sets. Similar to other adaptive methods in the S N area, the refinement test only performs on the geometry boundaries. Some spatial mesh effects that cause peaky interior angular fluxes but smooth boundary angular fluxes are not captured in this algorithm. 3. Results and Discussion 1 This adaptive algorithm is applied to three benchmarks to demonstrate its effectiveness and shortcomings for shielding problems featured high anisotropy throughout the duct regions with voids. In an effort to understand the behavior of the algorithm in this study as the error tolerance is changed, a convergence analysis is also performed in this section. All numerical calculations were implemented using the three-dimensional particle transport code ARES for assessment. 3.1. Kobayashi II Benchmark 2 The Kobayashi II benchmark is employed to demonstrate the ability of our angular adaptivity in problems with a straight duct. The problem is a rectangular block made of a half-scattering material containing a duct, as shown in Fig. 4 Fig. 4 . The void-duct inlet contains an isotropic source. The material properties and source strengths are listed in Table 1 Table 1 . Vacuum boundary conditions are imposed on all the external boundaries. The mesh used for all the calculations was composed of 60×100×60 Cartesian grids. Our adaptive quadrature sets were produced using an error tolerance of 0.06. The quantity of interest (QOI) is the scalar flux at spatial points along the red line ( X =65 cm, Z =65 cm), as shown in Fig. 4 Fig. 4 . 1 1 Regions 1 1 Source (n·cm −3 ·s −1 ) 1 1 Total cross-section (cm −1 ) 1 1 Scattering cross-section (cm −1 ) 1 1 1 1 1 1 1 1 0.1 1 1 0.05 1 1 2 1 1 0 1 1 10 −4 1 1 0.5×10 -4 1 1 3 1 1 0 1 1 0.1 1 1 0.05 The relative error of the neutron flux calculated by the adaptive method is compared to the results obtained using uniform P N T N and the ICLDFE quadrature sets, as listed in Table 2 Table 2 . The well-vetted published reference solutions [ 11 11 ] for the Kobayashi benchmark are defined as the references. The results demonstrate that the performance of the adaptive method is better than that of uniform quadrature sets with similar directions. The adaptive method with 3392 directions produces a root mean square (RMS) of 3.62×10 -3 , whereas the uniform quadrature sets with approximately the same directions produce an RMS of 1.99×10 -2 . Although the adaptivity method is slightly more expensive than the uniform quadrature sets because of the error estimates and quadrature set updates, its increased accuracy usually compensates for this cost. With a given RMS accuracy of 3.5×10 -3 , the adaptive method shows a significant decrease in the number of directions compared with the ICLDFE-S 5 quadrature sets with over 20000 directions. The difference in directions can result in significant savings by a factor of 8.2 in run times. Figure 5 Figure 5 illustrates the arrangement of the quadrature points in each adaptive process. Only the angular regions oriented to the duct exhibit a strong forward anisotropy, resulting in higher resolution, whereas a lower-order quadrature set is used in other regions with relatively smooth angular fluxes. This indicates that our error estimations can efficiently capture anisotropic angular fluxes and determine the locations that require further refinement. 1 1 Key points (cm) 1 1 Reference   (n·cm -2 ·s −1   ) 1 1 Adaptive (3392)|Relative Error| (%) 1 1 ICLDFE-S 5 (20480) 1 1 P N T N -S 60 (3720) 1 1 (65,105,65) 1 1 8.61×10 0 1 1 0.41 1 1 0.31 1 1 0.06 1 1 (65,115,65) 1 1 2.16×10 0 1 1 0.54 1 1 0.38 1 1 0.23 1 1 (65,125,65) 1 1 8.93×10 -1 1 1 0.18 1 1 0.39 1 1 0.86 1 1 (65,135,65) 1 1 4.77×10 -1 1 1 0.62 1 1 0.37 1 1 0.08 1 1 (65,145,65) 1 1 2.89×10 -1 1 1 0.46 1 1 0.55 1 1 1.33 1 1 (65,155,65) 1 1 1.89×10 -1 1 1 0.20 1 1 0.50 1 1 3.62 1 1 (65,165,65) 1 1 1.31×10 -1 1 1 0.38 1 1 0.16 1 1 4.15 1 1 (65,175,65) 1 1 9.50×10 -2 1 1 0.18 1 1 0.18 1 1 2.45 1 1 (65,185,65) 1 1 7.12×10 -2 1 1 0.01 1 1 0.20 1 1 0.02 1 1 (65,195,65) 1 1 5.45×10 -2 1 1 0.03 1 1 0.06 1 1 0.92 1 1 RMS 1 1 / 1 1 3.62×10 -3 1 1 3.44×10 -3 1 1 1.99×10 -2 1 1 Time/s 1 1 / 1 1 1919 1 1 15765 1 1 1876 One of the key parameters in this adaptive method is the user-defined error tolerance. An excessively small tolerance may result in additional unnecessary quadrature points, in contrast, an excessively large one does not employ sufficient refinement to acquire the desired resolutions. Therefore, a convergence analysis using a varying error tolerance from 0.8 to 0.01 is performed. Figure 6 Figure 6 plots the RMS of the scalar fluxes at key points as a function of the number of directions compared with the adaptive method with uniform ICLDFE quadrature sets. The results illustrate that the adaptive method converges an asymptotic accuracy faster as the angular region is refined compared to the uniform ICLDFE quadrature sets. The adaptive method produces more accurate results and fewer errors with similar directions. It should be noted that the error plateaus around an RMS of approximately 3.5×10 -3 . The RMS does not reduce further even if the number of directions is added because the error is dominated by spatial discretization, while the angular discretization error is sufficiently small. Refined spatial meshes or more accurate spatial discretization schemes should be considered. Overall, the adaptive method reduces the memory requirement and computational time without losing accuracy by refining only the required angular domains. 3.2. Kobayashi III Benchmark 2 The Kobayashi III benchmark is a rectangular block made of a half-scattering material containing complicated bends, as shown in Fig. 7 Fig. 7 . The calculation conditions are the same as those of the Kobayashi II benchmark described in Sect. 3.1. Our adaptive quadrature sets were produced using an error tolerance of 0.03. The quantity of interest (QOI) is the scalar flux at spatial points along the red line ( X =65 cm, Z =65 cm), as shown in Fig. 7 Fig. 7 . A similar analysis was performed in this problem. Table 3 Table 3 lists the relative errors of the neutron fluxes obtained from the adaptive method and uniform quadrature sets. The results also indicate that the adaptive method outperforms uniform quadrature sets with similar directions. The adaptivity with 3392 directions produces an RMS of 6.39×10 -2 , whereas the uniform quadrature sets with 3720 directions produce an RMS of 1.70×10 -2 . The adaptive method shows a significant decrease in the number of directions compared with ICLDFE-S 5 quadrature sets with over 20000 directions for a given RMS accuracy of 6.0×10 -3 . A convergence analysis with a varying error tolerance from 0.1 to 0.005 is performed. Figure 8 Figure 8 plots the RMS of the relative errors of the QOI as the number of directions using different quadrature sets. The results illustrate a similar trend as the Kobayashi II problem. The adaptive method converges faster than the uniform ICLDFE quadrature sets as the number of directions increases, and the error plateaus around an RMS of approximately 5.8×10 -3 . 1 1 Key points (cm) 1 1 Reference (n·cm -2 ·s −1 ) 1 1 Adaptive (3392)|Relative Error| (%) 1 1 ICLDFES 5 (20480) 1 1 P N T N -S 60 (3720) 1 1 (65,105,65) 1 1 8.62×10 0 1 1 0.40 1 1 0.29 1 1 0.35 1 1 (65,115,65) 1 1 2.16×10 0 1 1 0.53 1 1 0.36 1 1 0.56 1 1 (65,125,65) 1 1 8.94×10 -1 1 1 0.21 1 1 0.37 1 1 0.89 1 1 (65,135,65) 1 1 4.78×10 -1 1 1 0.59 1 1 0.34 1 1 0.60 1 1 (65,145,65) 1 1 2.89×10 -1 1 1 0.43 1 1 0.52 1 1 0.82 1 1 (65,155,65) 1 1 1.93×10 -1 1 1 0.23 1 1 0.41 1 1 2.99 1 1 (65,165,65) 1 1 1.05×10 -1 1 1 0.00 1 1 0.23 1 1 3.42 1 1 (65,175,65) 1 1 3.38×10 -2 1 1 0.43 1 1 0.13 1 1 1.36 1 1 (65,185,65) 1 1 1.08×10 -2 1 1 1.07 1 1 0.90 1 1 0.94 1 1 (65,195,65) 1 1 3.40×10 -3 1 1 1.30 1 1 1.26 1 1 1.83 1 1 RMS 1 1 / 1 1 6.39×10 -3 1 1 5.81×10 -3 1 1 1.70×10 -2 1 1 Time (s) 1 1 / 1 1 1927 1 1 16740 1 1 1617 The spatial point at the duct outlet (95 cm, 195 cm, 95 cm) is of particular concern in this problem. Figure 9 Figure 9 plots the RMS of the relative errors at the duct outlet as the number of directions using the adaptive method and uniform quadrature sets. Obviously, our adaptive method does not show benefits at the duct outlet compared with the uniform quadrature sets. The arrangement of the quadrature points in each adaptive process in this problem is the same as the Kobayashi II problem shown in Fig. 5 Fig. 5 . This indicates that our error estimations cannot exactly capture the anisotropic angular fluxes caused by the bend ducts. This adaptivity is introduced throughout the problem domain instead of in different regions, thus, the error estimate is only tested on the geometric boundaries, and peaky interior angular fluxes are not captured in this algorithm. However, the error estimation based on multi-spatial regions will be the focus of future work to attain a more accurate error distribution and more effective angular adaptivity. 3.3. IRI-TUB Research Reactor 2 The IRI-TUB research reactor, located in the Nuclear Technology Department of the Technical University of Budapest, is used to study particle transport in a straight or bent duct [ 12 12 ] . This particular problem is employed to demonstrate the ability of our angular adaptivity to solve complicated multi-group problems containing a straight duct. This reactor consists of a core and a large radiation channel. There are 24 fuel rods with a size of 7.2 cm × 7.2 cm asymmetrically arranged in the active area. The active core is surrounded by graphite and water as the reflective layer, and the shielding layer outside the core is metallic aluminum. A large radiation channel filled with concrete is placed 25 cm from the core, and a cylindrical duct with a radius of 11.8 cm and a length of 187 cm is built in. The outer side of the duct is surrounded by a layer of stainless steel with a thickness of 4.5 mm. Detectors are set up along the central axis at distances of 0, 67, 121, 148, and 175 cm from the entrance of the duct, respectively, as shown by the red squares in Figure 10 Figure 10 . The problem features highly anisotropic fluxes throughout the duct regions and highly anisotropic scattering due to 12 C in graphite. This type of problem poses a difficult challenge for angular discretization in the S N method for obvious reasons and is one of the motivations discussed before. The mesh used for all calculations was composed of 45×153×45 Cartesian grids. The problem uses 3rd-order (P 3 ) scattering along with a 199-group structure for energy discretization [ 13 13 ] . The cross-section messages were obtained from the multi-group cross-section library KASHIL-E70 [ 14 14 ] . Our adaptive quadrature sets were produced using an error tolerance of 0.22. Table 4 Table 4 lists the calculated and experimental values of the <math id="M26"><mrow><mmultiscripts><mtext>I</mtext><mprescripts/><mi/><mrow><mtext>115</mtext></mrow></mmultiscripts><mtext>n(n,n')</mtext><mmultiscripts><mtext>I</mtext><mprescripts/><mi/><mrow><mtext>115m</mtext></mrow></mmultiscripts><mtext>n</mtext></mrow></math> reactive rate using the angular adaptive method and uniform quadrature sets with similar directions, where C/E (AA) is the ratio of the calculation results obtained from our angular adaptive method (5024 directions) to the experimental value, C/E(P N T N ) is the ratio of the calculation results using the P N T N -S 70 quadrature sets (5040 directions) to the experimental values, and C/E (REF) is the ratio of the calculated value in the reference using the two-dimensional discrete ordinates program DOT3.5 [ 15 15 ] coupled with the three-dimensional Monte Carlo program MORSE-SGC/S [ 16 16 ] to the experimental value. The comparison in Table 4 Table 4 indicates that the adaptive method leads to results closer to the measurement among all the three methods. The relative error of the P N T N -S 70 quadrature sets between the calculated result and measurement is within 16% located near the duct outlet. Generally, the angular flux tends to be serious anisotropy as the distance from the duct inlet increases. Although the P N T N -S 70 quadrature sets have 5040 discrete directions, they still have the inability to achieve the necessary resolution to exactly describe the angular flux along the duct. The angular distribution in each adaptive iteration is shown in Figure 11 Figure 11 . As the adaptivity continues, the angular regions along or near the duct adopt higher order quadrature sets to improve the local integral accuracy, while other regions adopt lower quadrature sets because the angular flux is relatively smooth making it easier for the quadrature sets to integrate. The adaptive method can effectively capture the angular flux along the duct and the error near the duct outlet is less than 4%. Overall, the computational results demonstrate that the performance of our adaptive method is better than that of the uniform quadrature sets with similar directions. It should be noted that a larger error of approximately 10% both in our adaptive method and uniform quadrature sets appears near the duct inlet where angular effects are not dominant. Further analysis of the spatial effects should be considered in future work. 1 1 Distances from entrance (cm) 1 1 Measured (s −1 ) 1 1 Calculated (AA) (s −1 ) 1 1 Calculated (P N T N ) (s −1 ) 1 1 C/M(AA) 1 1 C/M (P N T N ) 1 1 C/M (REF) 1 1 0 1 1 5.57×10 -16 1 1 4.95×10 -16 1 1 4.96×10 -16 1 1 0.89 1 1 0.89 1 1 0.77 1 1 67 1 1 3.17×10 -17 1 1 2.79×10 -17 1 1 2.79×10 -17 1 1 0.88 1 1 0.88 1 1 0.63 1 1 121 1 1 7.70×10 -18 1 1 7.17×10 -18 1 1 8.24×10 -18 1 1 0.93 1 1 1.07 1 1 0.54 1 1 148 1 1 4.61×10 -18 1 1 4.46×10 -18 1 1 5.33×10 -18 1 1 0.97 1 1 1.16 1 1 0.54 1 1 175 1 1 3.01×10 -18 1 1 3.12×10 -18 1 1 3.44×10 -18 1 1 1.04 1 1 1.14 1 1 0.54 Normalized angular fluxes profiles at the inner part of the duct in energy Group 1 (1.97×10 7 eV to 1.73×10 7 eV), Group 30 (5.49×10 6 eV to 5.22×10 6 eV), Group 100 (1.43×10 5 eV to 1.36×10 5 eV), and Group 170 (1.3 eV to 1.13 eV), respectively, are exhibited in Figure 12 Figure 12 . For neutrons with high energy, the angular flux in the first group tends to be the most unsmooth or peaky and anisotropy concentrates on similar angular regions in different energy groups. Thermal fluxes tend to be isotropic as the neutrons slow down. Less directions, even uniform lower-order quadrature sets, can integrate angular fluxes exactly. This fact conforms to the previously discussed physical analysis. Figure 13 Figure 13 compares the energy spectrum at the duct inlet (denoted by #1) and 148 cm from the entrance (denoted by #4) of the duct. This figure demonstrates that the computational results are in good agreement with the measurement from 10 -2 eV to 10 7 eV, which indicates the feasibility of our adaptive method in multi-group problems. Because thermal fluxes often exhibit a smoother angular behavior than fast fluxes due to neutron moderation, an angular distribution tailored specifically for fast fluxes is likely too fine for the thermal flux solution, resulting in an unnecessary expense. Hence, considering different refined quadrature sets that would be applicable to different energy groups is an efficient option to further reduce the computational cost. 4. Conclusion 1 A methodology for implementing angular adaptivity in the S N method has been presented. This algorithm is based on discontinuous finite element discrete quadrature sets over an icosahedron to achieve local refinement in the one-twentieth angular region. Numerical results through detailed quantitative analysis indicates that our adaptive method produces a more efficient arrangement of angular unknowns for a given accuracy compared with uniform quadrature sets. In Kobayashi problems with straight ducts, the adaptive method exhibits the same accuracy with a speed-up of approximately eight compared to uniform quadrature sets. For more difficult multi-group problems, our adaptive method also outperforms uniform angular discretization with similar angular unknowns. Although more comparisons with practical problems are necessary before the efficiency of our method is fully demonstrated, the results presented here are illuminating. It should be noted that our adaptivity is introduced throughout the problem domain instead of in different regions, thus, peaky interior angular fluxes caused by bend ducts would not be captured in this algorithm. Future studies will integrate this algorithm with multi-region error estimations for different energy groups and a parallel strategy. REFERENCES E. E. Lewis , W. F. Miller , Computational Methods of Neutron Transport ( John Wiley & Sons Inc , New Jersey , 1984 ) Computational Methods of Neutron Transport J.P. Jenal , P.J. Erickson , Generation of a computer library for discrete ordinates quadrature sets. Report ORNL/TM-6023 , ( 1977 ). doi: 10.2172/5237525 http://doi.org/10.2172/5237525 Generation of a computer library for discrete ordinates quadrature sets. Report ORNL/TM-6023 G. Longoni , A. Haghighat , Development of new quadrature sets with the 'ordinate splitting' technique , ANS International Meeting on Mathematical Methods for Nuclear Applications ( Salt Lake City, USA , 2001 ) Development of new quadrature sets with the 'ordinate splitting' technique G. Longoni , A. Haghighat , Development and application of the regional angular refinement technique and its application to non-conventional problems . Proceedings of PHYSOR 2002 ANS Topical Meeting , ( Seoul, Korea , 2002 ) Development and application of the regional angular refinement technique and its application to non-conventional problems J. C. Stone , Adaptive discrete-ordinates algorithms and strategies , ( Texas A&M University at Texas , 2007 ) Adaptive discrete-ordinates algorithms and strategies J. J. Jarrell , An adaptive angular discretization method for neutral-particle transport in three-dimensional geometries , ( Texas A&M University at Texas , 2010 ) An adaptive angular discretization method for neutral-particle transport in three-dimensional geometries C. Y. Lau , Adaptive discrete-ordinates quadratures based on discontinuous finite elements over spherical quadrilaterals , ( Texas A&M University at Texas , 2016 ) Adaptive discrete-ordinates quadratures based on discontinuous finite elements over spherical quadrilaterals B. Zhang , L. Zhang , C. Liu , et al. , Goal-oriented regional angular adaptive algorithm for the SN equations . Nucl. Sci. Eng. 189 ( 2 ), 120 - 134 ( 2018 ). doi: 10.1080/00295639.2017.1394085 http://doi.org/10.1080/00295639.2017.1394085 Goal-oriented regional angular adaptive algorithm for the SN equations N. Dai , B. Zhang , Y. X. Chen , Discontinuous finite-element quadrature sets based on icosahedron for the discrete ordinates method . Nucl. Eng. Technol. 52 ( 6 ), 1137 - 1147 ( 2020 ). doi: 10.1016/j.net.2019.11.025 http://doi.org/10.1016/j.net.2019.11.025 Discontinuous finite-element quadrature sets based on icosahedron for the discrete ordinates method Y. X. Chen , B. Zhang , L. Zhang , ARES: a parallel discrete ordinates transport code for radiation shielding applications and reactor physics analysis . Sci. Technol. Nucl. Ins. 2017 , pp. 1 - 10 ( 2017 ). doi: 10.1155/2017/2596727 http://doi.org/10.1155/2017/2596727 ARES: a parallel discrete ordinates transport code for radiation shielding applications and reactor physics analysis K. Kobayashi , N. Sugimura , Y. Nagaya , 3D radiation transport benchmark problems and results for simple geometries with void region . Prog. Nuc. Energ. 39 ( 2 ), 119 - 144 ( 2001 ). doi: 10.1016/S0149-1970(01)00007-5 http://doi.org/10.1016/S0149-1970(01)00007-5 3D radiation transport benchmark problems and results for simple geometries with void region J. L. Kloosterman , J. E. Hoogenboom , E. M. Zsolnay , Experiments and calculations on neutron streaming through bent ducts . J. Nucl. Sci. Technol. 30 ( 7 ), 611 - 627 ( 1993 ). doi: 10.3327/jnst.30.611 http://doi.org/10.3327/jnst.30.611 Experiments and calculations on neutron streaming through bent ducts C. Liu , X. L. Hu , B. Zhang , et al , Numerical study of scattering Legendre moments and effect of anisotropic scattering on SN shielding calculation . Nucl. Sci. Tech. 30 ( 11 ), 161 ( 2019 ). doi: 10.1007/s41365-019-0695-y http://doi.org/10.1007/s41365-019-0695-y Numerical study of scattering Legendre moments and effect of anisotropic scattering on SN shielding calculation D. H. Kim , C. Gil , Y. Lee , Validation of an ENDF/B-VII.0-based neutron and photon shielding library in MATXS-format . J. Korean Phys. Soc. 59 , pp. 1199 - 1202 ( 2011 ). doi: 10.3938/jkps.59.1199 http://doi.org/10.3938/jkps.59.1199 Validation of an ENDF/B-VII.0-based neutron and photon shielding library in MATXS-format W. A. Rhoades , F. R. Mynatt , DOT III: two-dimensional discrete ordinates transport code. Report ORNL/TM-4280 , ( 1973 ). doi: 10.2172/4395105 http://doi.org/10.2172/4395105 DOT III: two-dimensional discrete ordinates transport code. Report ORNL/TM-4280 E. A. Straker , P. N. Stevens , D. C. Irving , MORSE code: a multigroup neutron and gamma-ray Monte Carlo transport code. Report ORNL-4585 ( 1970 ). doi: 10.2172/4076246 http://doi.org/10.2172/4076246 MORSE code: a multigroup neutron and gamma-ray Monte Carlo transport code. Report ORNL-4585 All authors contributed to the study conception and design. Material preparation, data collection and analysis were performed by Ni Dai, Bin Zhang, Yi-Xue Chen and Dao-Gang Lu. The first draft of the manuscript was written by Ni Dai and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript. 10.1007/s41365-021-00934-7 National Natural Science Foundation of China (No. 11975097) and the Fundamental Research Funds for the Central Universities (No. 2019MS038). 2021-03-25 2021-08-02 2021-08-05 2021-09-08 2021-09-01 2022-10-11 © China Science Publishing & Media Ltd. (Science Press), Shanghai Institute of Applied Physics, the Chinese Academy of Sciences, ChineseNuclear Society and Springer Nature Singapore Pte Ltd. 2021 This is a post-peer-review, pre-copyedit version of an article published in Nuclear Science and Techniques. The final authenticated version is available online at: http://dx.doi.org/10.1007/s41365-021-00934-7 http://dx.doi.org/10.1007/s41365-021-00934-7 . 2021 Nuclear Science and Techniques 09 32