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

Ni Dai ? Bin Zhang ? Yi-Xue Chen ? Dao-Gang Lu

Abstract The discrete ordinates (SN) 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 SN 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 the 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.

Keywords Shielding calculation · Discrete ordinates method·Angular adaptivity·Discontinuous finite element

1 Introduction

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].The discrete ordinates (SN) 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 SNmethod,and the angular integral is replaced using quadrature sets with weights associated with these ordinates.However,angular discretization errors introducedthroughquadraturesetsintheSNmethod,inaway,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.Uniformlyincreasingthequadraturedirectionsresults inanincreaseinthesizeofthe calculation,andthusitisnot 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 of 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 refinementina hemisphere[2].In2001,Longonidevelopedthe longitudinal segmentation (OS) [3] and regional angular refinement (RAR) [4] to add discrete directions only in the focus angular regions for different problems.However,refinements inthe 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 SNmethod.Stone first studied the adaptive method withinthe SNmethodin 2007.The focusof his work isbasedon traditional quadrature sets to automate refinement in different quadrature regions[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].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].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 algorithmbasedonthe adjointvaluetheorytomitigate rayeffects[8].However,the accuracy of the appliedquadrature 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].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 Sect.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 Sect.4.

2 Angular Adaptivity Algorithm

2.1 Adaptive Quadrature Sets

The adaptive quadrature sets applied in this study are the ICLDFE quadrature sets.Their construction is shown in Fig.1.

Fig.1 (Color online) Construction of ICLDFE quadrature sets

The angular domain is first divided into spherical triangles by inscribing an icosahedron into the unit sphere,as shown in 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 4N(Nis the order of the quadrature sets).Linear discontinuous finite element basis functions in the directions are established over spherical triangles.The weightswiare then calculated by integrating the underlying basis functions in the corresponding spherical triangles.The basis functionsb（Ω）are represented by a spherical harmonic expansion of the zeroth and first order at each quadrature point denoted asi,wherea0，i,a1，i,a2，ianda3，iare unknown coefficients.The angular vector Ω is described by the variables μ,η,and ξ,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.

Fig.2 (Color online) Projection from the icosahedron to the unit sphere

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 the previous research[9]has 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 × 4N(i.e.,in the case of ICLDFE-S2,the number of directions is 20 × 42=320).They are invariant under 72-degree rotations about the three axes joining the extreme opposite vertices and inversions about the origin.However,it must be noted that reflection in SNcalculations requires that the quadrature sets possess 90-degree symmetry and thus the 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

The common adaptive strategy applies a progressive rule by starting with 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.

where the testing directionsi=1…Mare the directions of the higher-order quadrature sets in one spherical triangle.Δφtestis the maximum relative error between the calculated angular fluxes and interpolated results of all the spatial boundary grids.denotes the angular flux through the boundary surfaces calculated by a transport sweep in the testing directions using the latest scattering and the fixed source as well as the boundary conditions.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.

wherec0，c1，...，c15are the unknown coefficients to be sought.Equation (4) can be written as a matrix

where ｛μi，ηi，ξi｝，i=1，2，...4nare the discrete directions of the current quadrature sets.ψ0，ψ1，...，ψ4nare the known angular fluxes of the current directions,andnis the quadrature order.A simple inverse is needed to solve these coefficients in the spherical harmonic expansion when the initial quadrature sets default ICLDFE-S2quadrature 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 the equation:

where ｛μi，ηi，ξi｝，i=1，2，...4n+1are the directions of the testing quadrature sets,andare the angular fluxes through boundary surfaces at testing directions calculated by interpolation.A user-defined flux error tolerance ε is used to determine if refinement is necessary.To refine,the flux difference in the testing directions must be less than the tolerance ε on any meshes on the spatial boundary.

2.3 Implementation of Angular Adaptivity

The angular adaptivity was implemented in a three-dimensional particle transport code named ARES [10].The implementation of the angular adaptivity in the transport calculations is illustrated in 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 userdefined parameters,such as the initial quadrature set order,number of sweeps before the first adaptive sweep (Na),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 the 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.

Fig.3 (Color online)Implementation of the angular adaptive algorithm

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-S2as an example).The ICLDFE-S3quadrature 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 the 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 S1to S6have 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 usesmore computer time per unknown than the uniform quadrature sets because the investment is made to compare the calculated angular fluxes and the 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 costs 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 SNarea,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

This adaptive algorithm is applied to three benchmarks to demonstrate its effectiveness and shortcomings for shielding problems that 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

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.The void-duct inlet contains an isotropic source.The material properties and source strengths are listed in 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 (Color online) Configuration of Kobayashi II problem

The relative error of the neutron flux calculated by the adaptive method is compared to the results obtained using uniform PNTNand the ICLDFE quadrature sets,as listed in Table 2.The well-vetted published reference solutions[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-S5quadrature sets with over 20,000 directions.The difference in directions can result in significant savings by a factor of 8.2 in run times.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.

Table 1 Cross-section and source strength

Table 2 Relative error of QOI for the Kobayashi II problem

Table 3 Relative error of QOI for the Kobayashi III problem

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 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 requirements and computational time without losing accuracy by refining only the required angular domains.

3.2 Kobayashi III Benchmark

The Kobayashi III benchmark is a rectangular block made of a half-scattering material containing complicated bends,as shown in 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.6 (Color online) RMS of relative errors of QOI along 2A as number of directions

Fig.7 (Color online) Configuration of Kobayashi III problem

Fig.8 (Color online) RMS of relative errors of QOI along 3A as number of directions

Fig.9 (Color online) RMS of relative error at the duct outlet as number of directions

A similar analysis was performed in this problem.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-S5quadrature sets with over 20,000 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 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.

Table 4 Calculated and experimental values of 115In（n，n′）115mIn reactive rate

The spatial point at the duct outlet(95 cm,195 cm,and 95 cm) is of particular concern in this problem.Figure 9 plots the RMS of the relative errors at the duct outlet as the number of directions using the adaptive method and the 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.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.

Fig.5 (Color online) Arrangement of directions in each adaptive process:a initial,b first,c second,and d third

3.3 IRI-TUB Research Reactor

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].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 Fig.10.The problem features highly anisotropic fluxes throughout the duct regions and highly anisotropic scattering due to12C in graphite.This type of problem poses a difficult challenge for angular discretization in the SNmethod 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 (P3) scattering along with a 199-group structure for energy discretization [13].The cross-section messages were obtained from the multigroup cross-section library KASHIL-E70 [14].Our adaptive quadrature sets were produced using an error tolerance of 0.22.

Fig.10 (Color online) Configuration of IRI-TUB problem at Z=0 cm

Table 4 lists the calculated and experimental values of the115In（n，n′）115mIn 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(PNTN) is the ratio of the calculation results using the PNTN-S70quadrature 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] coupled with the three-dimensional Monte Carlo program MORSE-SGC/S [16] to the experimental value.The comparison in Table 4 indicates that the adaptive method leads to results closer to the measurement among all the three methods.The relative error of the PNTN-S70quadrature 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 PNTN-S70quadrature 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 Fig.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.

Fig.11 (Color online) Distribution of directions in each adaptive iteration:a initial,b first,c second,d third,and e fourth

Fig.12 (Color online)Normalized angular fluxes profiles at inner part of duct in energy a Group 1,b Group 30,c Group 100,and d Group 170

Fig.13 (Color online) Comparison with measured and calculated energy spectrum at #1 and #4

Normalized angular fluxes profiles at the inner part of the duct in energy Group 1 (1.97 × 107eV to 1.73 × 107eV),Group 30 (5.49 × 106eV to 5.22 × 106eV),Group 100 (1.43 × 105eV to 1.36 × 105eV),and Group 170 (1.3 eV to 1.13 eV),respectively,are exhibited in Fig.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 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–2to 107eV,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 costs.

4 Conclusion

A methodology for implementing angular adaptivity in the SNmethod has been presented.This algorithm is based on discontinuous finite element discrete quadrature sets over an icosahedron to achieve local refinement in the onetwentieth angular region.Numerical results through detailed quantitative analysis indicate 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.

Author contributionsAll 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.