Discontinuous Galerkin spatial discretisation of the neutron transport equation with pyramid finite elements and a discrete ordinate (SN) angular approximation

Abstract In finite element analysis it is well known that hexahedral elements are the preferred type of three dimensional element because of their accuracy and convergence properties. However, in general it is not possible to mesh complex geometry problems using purely hexahedral meshes. Indeed for highly complex geometries a mixture of hexahedra and tetrahedra is often required. However, in order to geometrically link hexahedra and tetrahedra, in a conforming finite element mesh, pyramid elements will be required. Until recently the basis functions of pyramid elements were not fully understood from a mathematical or computational perspective. Indeed only first-order pyramid basis functions were rigorously derived and used within the field of finite elements. This paper makes use of a method developed by Bergot that enables the generation of second and higher-order basis functions, applying them to finite element discretisations of the neutron transport equation in order to solve nuclear reactor physics, radiation shielding and nuclear criticality problems. The results demonstrate that the pyramid elements perform well in almost all cases in terms of both solution accuracy and convergence properties.


Introduction
This paper investigates the use of pyramid elements when computing solutions to the neutron transport equation using a discrete ordinate (S N ) angular discretisation and a discontinuous Galerkin finite element (DGFEM) spatial discretisation (Reed and Hill, 1973;Wareing et al., 2001). An efficient and accurate solution of the neutron transport equation is essential for a variety of nuclear reactor physics, radiation shielding and nuclear criticality safety assessment problems. When creating a finite element mesh most mesh generators will often employ tetrahedral elements, especially for highly complex geometries. This is because the tetrahedron is a simplex granting it maximum flexibility to mesh complex geometries (Frey and George, 2008), meaning that several robust and simple tetrahedral meshing algorithms exist, such as advancing front and Delaunay. When computing solutions to finite element problems hexahedral elements are, in general, more efficient than tetrahedral elements (Cifuentes and Kalbag, 1992;Benzley et al., 1995). However, one of the main issues in using hexahedral elements is that no general robust meshing algorithm exists for hexahedral elements for complex geometries (Schneiders, 2000;Puso and Solberg, 2006).
In order to mesh complex geometries while retaining the favourable properties of hexahedral elements some mesh generators will generate a mixed element mesh where predominantly hexahedral elements are used for the majority of the problem domain with a mixture of tetrahedra, wedges (triangular prisms), and pyramids used to mesh the difficult regions. Pyramid elements are also generated by octree based mesh generators as a part of their mesh refinement process (Dawes et al., 2009) and may be used in order to create a link between regions of high refinement and low refinement hexahedral elements without the need for hanging nodes.
The implementation of nodal finite elements for tetrahedral and hexahedral and wedge elements is well understood. However, pyramid elements will naturally have both polynomial and rational shape functions. Therefore, finding a set of basis functions which provides an effective and stable finite element solution can be challenging. A set of effective basis functions for first-order pyramids has been known for some time (Coulomb et al., 1997), but it was not until a study by (Bergot et al., 2010) that a system of generating correct basis functions for second-order and higher pyramid elements was established.
A previous paper (O'Malley et al., 2017a) examined the performance of pyramid elements generated using the Bergot method when solving the neutron diffusion equation. This paper examines instead their performance when solving the neutron transport equation. The neutron transport equation is hyperbolic whereas the neutron diffusion equation is elliptic, therefore quite different computational behaviour is possible. A three dimensional discontinuous Galerkin finite element (DGFEM) discretisation with discrete ordinates (S N ) angular approximation is used (Wareing et al., 2001). Various nuclear reactor physics, radiation shielding and nuclear criticality problems will be solved in order to demonstrate the computational performance of the pyramid elements. These numerical solutions will then be compared against results obtained with more conventional finite element types such as tetrahedra and hexahedra.

Neutron transport overview
This paper solves the transport equation using a discontinuous (S N ) method established in (Wareing et al., 2001). The monoenergetic steady state transport equation, solving for the angular neutron flux wðr; b XÞ (cm À2 s À1 Sr À1 ) at position r and neutron direc- where R t ðrÞ (cm À1 ) represents the total macroscopic material cross section of the medium at a given position and R s ðr; b X 0 ! b XÞ (cm À1 -Sr À1 ) represents the macroscopic differential scattering cross section from direction b X 0 to direction b X at a given position. Q ðr; b XÞ (cm À3 s À1 Sr À1 ) is a source term representing any neutrons entering the system.
For the case of a neutron transport problem discretised over multiple energy groups the problem must be solved by iterating through the energy groups performing mono-energetic calculations in each one. The removal from each system due to scatter to another energy group is tracked and neutrons entering due to scatter from other energy groups are added to the source term Q. In the case where there is a fission source an iterative eigenvalue method is used (Warsa et al., 2004a). In all cases the basic monoenergetic system defined in Eq. 1 is used.
The scalar neutron flux /ðrÞ may be defined as: XÞ. Let L represent the streaming and removal of neutrons b X Á r þ R t ðrÞ , the inversion of this operator is the equivalent of inverting the local element matrices as described in (Wareing et al., 2001). The operator D represents the conversion of angular flux to scalar flux such that / ¼ Dw and M represents a mapping of scalar flux to angular flux (although w -M/). S is used to represent the action of neutron scatter cross-sections. This notation allows the transport equation shown in Eq.
(1) to be expressed as shown in Eq.
(2) (Warsa et al., 2004b): By multiplying through by the term DL À1 the equation may be expressed in terms of scalar flux: leading to: The expression DL À1 Q is simple to calculate and may be referred to as the swept flux. The operator ðI À DL À1 MSÞ is well defined, meaning that calculating ðI À DL À1 MSÞ/ for a known / is trivial. For an unknown /, letting A ¼ ðI À DL À1 MSÞ and b ¼ DL À1 Q allows the transport equation to be posed as a standard Ax ¼ b matrix problem and it may then be solved using an standard linear algebra techniques (Warsa et al., 2004b). Direct solution methods are impractical for most practical neutron transport problems so indirect or iterative methods are usually employed, this paper uses a GMRES solver for this purpose (Saad, 2003).

Pyramid elements
In order to create a set of basis functions for an element it is first necessary to define a set of expressions which form an orthogonal basis of the finite element space b P r (Bathe, 1996). These are defined here through Jacobi polynomials, where the function P a;b m ðxÞ is a Jacobi polynomial function of order m and a weighting of ð1 À xÞ a ð1 þ xÞ b (Szego, 1975).
For a finite element order r, (Bergot et al., 2010) defines a formula for the orthogonal bases wðx; y; zÞ of a pyramid element with the formula: where: 0 6 i 6 r; 0 6 j 6 r; 0 6 k 6 r À maxði; jÞ ð 6Þ This equation will generate n orthogonal bases where: meaning that n is the number of degrees of freedom of a pyramid element with order r.
Once the orthogonal bases are established it is possible to define a set of basis functions using the nodal positions of a reference element. This method is specified in (Bergot et al., 2010). The typical reference element for a first-order five node pyramid is shown in Fig. 1.

Results
This section presents a variety of neutron transport problems to be solved using pyramid elements of both first and second-order. Some problems will be constructed entirely of pyramid elements in a manner which is useful for testing but not representative of how these elements would be used in real reactor physics applications. Others will use the pyramids in a more realistic way within a mixed element mesh in order to mesh more complex geometries. All problems will compare the computational performance of using pyramid elements against other more conventional element types.
Many of the problems studied here will use the neutron multiplication factor K eff as a benchmark for solution accuracy. K eff is a property of system in which neutrons are being generated by nuclear fission (multiplying systems), representing the ratio of neutrons entering and leaving the system. In multiplying systems the absolute neutron flux is not fixed, only its relative distribution, so in these cases K eff is a useful quantity for benchmarking the accuracy of the solution. The accuracy of problems using pyramid elements will be judged not just by the absolute error of K eff from a reference solution, as this is dependant on a great many factors, but on the error relative to that when using other elements and the convergence towards a value as mesh refinement is increased. It is not the intent of this paper to try to show that pyramids have improved computational performance over other element types. The purpose instead is to demonstrate that they provide a reasonable alternative to other element types in situations where the problem geometry makes their use desirable for meshing purposes. For all timing results in this paper the problems were run in parallel on a computer with two Intel Xeon E5-2680 2.80 GHz processors, each of which possesses ten cores with hyper-threading for a total of forty processes.

L 2 -error analysis
A method of manufactured solutions (MMS) problem is performed for a homogeneous cubic problem of dimensions 1.0 cm Â 1.0 cm Â 1.0 cm. The exact solution of the MMS problem in terms of scalar flux / is chosen to be: The domain is discretised into structured hexahedral and structured pyramid meshes of first and second-order with varying levels of mesh refinement, a transport solution is obtained, and the L 2error is calculated and plotted alongside element characteristic length l c . Here for a problem with N E elements the characteristic length is defined as The L 2 -error results for this problem are plotted in Fig. 2. These results show that when plotted on a logarithmic scale the L 2 -error scales linearly with the characteristic length in all cases. For firstorder elements, both with hexahedra and pyramids, this scaling occurs with a gradient of roughly 2 and for second-order cases the gradient is roughly 3. These results match the expected results for an L 2 -error analysis of a finite element problem. This demonstrates that the pyramid elements are producing the expected computational convergence rates for first and second-order finite elements.

Takeda 2 nuclear reactor physics benchmark problem
The Takeda benchmarks (Takeda and Ikeda, 1991) are a set of three-dimensional nuclear reactor physics benchmarks. This section will study results for the second of the Takeda benchmark problems, Takeda 2, which is a core for a fast breeder reactor. The geometry of this problem is shown in Fig. 3a and b. There are two variations of this problem; in case 1 the control rod (CR) is fully withdrawn leaving only a sodium filled space (CRP), in case 2 the control rod is half inserted and the other half is filled with sodium.
The Takeda benchmark specification (Takeda and Ikeda, 1991) lists criticality results obtained by various individuals for both cases of this problem, displayed in Table 1.  (Bedrosian, 1992) it is impossible to choose polynomial shape functions for a pyramid element while preserving conformity with other elements. The Bergot method (Bergot et al., 2010) outlined above creates an orthogonal base which uses Jacobi polynomials up to order r in the x and y coordinates but is not polynomial in the z coordiante. This enables the generation of a set of non-polynomial basis functions for the pyramid element which conform with standard hexahedral and tetrahedral elements. Due to the simple structure of the geometry it is relatively straightforward to create a mesh of the Takeda problem using structured cubic hexahedral elements. It is also possible, by splitting each hexahedra into six pyramids, to create a fully structured pyramid mesh in the same manner. This splitting was achieved using a python script to create a structured mesh of hexahedra, a node was then added to the centre point of each hexahedra and six pyramids were created each of which has this central node as their apex and a face of the hexahedra as their base. The problem is run with first and second-order hexahedra and pyramids using S 8 angular quadrature.
For the results in Table 2 with hexahedral elements the element length refers to the side length of the elements which are all cubic. For the pyramid problems it represents the side length of a cube which has been split into six identical pyramids. These results demonstrate clearly that the usage of pyramid elements in place of hexahedra has not impacted on the problem solution, for the higher refinement second-order cases the K eff is identical for both hexahedra and pyramids to an accuracy of five significant figures.

Kobayashi dog leg duct radiation shielding problem
In (Kobayashi and Sugimura, 2001) three radiation shielding benchmark solutions are presented. As these problems are monoenergetic fixed source problems they are well suited for studying the spatial convergence of numerical solutions. For the purposes of this spatial convergence study the third Kobayashi problem, which is a dog leg duct radiation shielding problem, is analysed. The dog leg duct featured in this problem makes it challenging for S N neutron transport solvers. Indeed a high level of angular quadrature will be required in order to adequately resolve the anisotropic angular flux through the dog leg duct. The geometry of the Kobayashi problem is shown in Fig. 4a, b and c, with the material properties for two cases shown in Table 3.
As well as defining the problem the Kobayashi paper also defines scalar flux values at various locations for both cases of Fig. 3. Geometry of the Takeda 2 nuclear reactor physics benchmark problem (Takeda and Ikeda, 1991). Table 1 K eff results for Takeda 2 benchmark from various studies using S 8 angular quadrature. If a length is given it refers to the element size (Takeda and Ikeda, 1991). Case 1 is for control rod fully removed and case 2 is for the control rod half inserted.

Study
Case 1 (Nakagawa et al., 1990). For case 2 only Monte Carlo values are provided as an analytical solution is not obtainable. The flux values are defined at 22 discrete points in the system. When the problem is solved in our neutron transport code the scalar flux values at each of these points are taken to determine the accuracy of the solution. By subtracting our value from the reference value and dividing by the reference value a relative error may be obtained for each point, multiplying this value by 100 yields the percentage error for each point. For case 2, the scatter case, these percentage errors are plotted in Fig. 5a, b, and c for various discretisations of the Kobayashi problem using structured hexahedra and pyramids and with a variety of orders of S N quadrature. These figures show how the accuracy of each case varies with S N order and how they compare to the results from GMVP. The presence of a duct and the low scatter within this problem even for case 2 makes it particularly challenging for a S N code, and therefore high angular quadratures are necessary for good accuracy as may be seen from the results for the S 32 case and below. For the S 64 and S 128 cases however there is significantly better accuracy with the error at some points falling within the 1r Monte Carlo error. This is the case for both hexahedra and pyramids, although S 128 for the higher refinement pyramid mesh was not possible due to the large memory requirement exceeding what was available (the computer used had approximately 256 Gb of RAM).
As well as calculating the percentage error at each of the 22 points it is possible to create a more general accuracy metric for the whole problem by taking the square root of the summation of the square of the error at all points. This gives an accuracy metric which will be referred to as the RMS error, explicitly defined in Eq. 9: Using this RMS error formula it is possible to quickly obtain a rough estimate of the accuracy of a large number of solutions of the Kobayashi problem, as seen in Figs. 6 and 7. Here various spatial refinements of the problem are run using structured Fig. 4. Geometry of the Kobayashi dog leg duct radiation shielding problem (Kobayashi and Sugimura, 2001). Table 3 Material properties for both cases of the Kobayashi problem.

Rt
Rs Rs hexahedral, tetrahedral and pyramid elements for S 128 angular quadrature. The high memory requirements of S 128 prevented any higher refinements than those shown from being studied, but some trends may still be observed. Second-order elements perform better than first order and hexahedral elements are generally superior, which is to be expected. The pyramid and tetrahedral elements both perform well however. It is notable that pyramids appear to be trending towards superior performance for higher refinement second-order elements, but more data points for higher spatial refinement cases would be necessary to confirm this trend.

GODIVA nuclear criticality benchmark problem
The GODIVA nuclear criticality benchmark problem, defined in specification 1-A1 (page 12) of (ANL, 1972), is a homogeneous sphere of highly enriched uranium in a vacuum. Since the sphere is fully symmetrical in the x, y and z planes it may be modelled as a single octant of a sphere with reflective boundaries on these planes. The simplest way to mesh a spherical body is with unstructured tetrahedra. It is possible to generate a hexahedral mesh on a sphere problem but doing so requires significant pre-processing effort (Wareing et al., 2001). However, it is possible to create a structured hexahedral mesh for a sphere octant using some tetrahedral, pyramid and wedge elements. Fig. 8 illustrates how the mesh is formed.
The GODIVA benchmark problem was solved using meshes of the kind shown in Fig. 8 with varying refinements and also meshes formed of exclusively unstructured tetrahedra. These will be known as the structured and unstructured cases respectively. A very high spatial refinement variant of the structured case with 137,500 s-order elements was used to calculate a benchmark criticality (K eff ) for the problem, which was calculated as 0.99845. This value is acceptably close to the values for K eff given in (ANL, 1972) for S 8 angular quadrature. Then the GODIVA problem was run with first and second-order elements for the unstructured and structured case, with the error relative to this reference K eff plotted against the solution time. An S 8 angular quadrature was used for all solutions.
Figs. 9 and 10 display the GODIVA benchmark problem results. Fig. 9 shows that as the spatial refinement increases all of the different meshes are converging on the same numerical value. Fig. 10 demonstrates that, while for lower spatial refinement problems the unstructured mesh performs better, when high refinement accuracy is required a structured mesh is much more computationally efficient. This study demonstrates how these pyramid elements may be used as part of a mixed mesh in order to easily generate well performing structured meshes for difficult geometry.

Linking regions with varying spatial refinement
Using a structure of pyramid and tetrahedral elements it is possible to link regions of low and high refinement hexahedral elements without the need for hanging nodes, an example of such a mesh structure is demonstrated in (O'Malley et al., 2017a). If the faces of a low refinement cubic mesh are all to be connected to a high refinement outer cube then an additional structure of pyramids, tetrahedra, and prisms is needed in order to connect the edges. A basic example of where this might be useful for a neutron transport problem is a homogeneous cube with bare boundaries, since the flux profile will change relatively little in the centre. It is therefore feasible to use a low refinement mesh there without a significant loss of accuracy.
An octant of such a problem with reflective boundaries on three faces is generated, illustrated in Fig. 11, which with the reflective boundaries is equivalent to a full cube with the low refinement region in the centre and the high refinement region surrounding it. The octant is 2 cm Â 2 cm Â 2 cm in dimension and the material properties are homogeneous and monoenergetic with R a ¼ 1:0cm À1 ; R s ¼ 1:0cm À1 ; mR f ¼ 1:0cm À1 , and R t ¼ 2:0cm À1 .
The problem is first run for a fixed refinement mesh consisting purely of hexahedra and then for the problem with linked regions of varying refinement. The criticality K eff is calculated and the absolute error between the calculated K eff and a reference K eff determined from a particularly high spatial refinement case is obtained.
The results in Table 4 show that even for this simple example, where there is only a modest reduction in the number of elements present, there is a noticeable saving of computational time. In addition, the accuracy loss is fairly minor, particularly for the secondorder elements.

Highly scattering Kobayashi dog leg problem with reflector
Diffusion synthetic acceleration (DSA) is a vital tool for obtaining computationally efficient convergence of neutron transport problems containing highly diffusive regions (Adams and Martin,  1992). In order to test the computational performance of the pyramid elements with a transport solver using DSA, a variation of the Kobayashi dog leg duct problem is used. The material properties for the source and duct are the same as for Kobayashi Case 2 studied previously but this time the macroscopic scattering crosssection in the body is much higher. In addition to this a 20 cm thick high scatter reflector is added to boundaries of the problem which previously had a reflective boundary condition and all boundaries are treated as bare. These changes make the problem significantly more challenging for standard transport sweep based solvers and therefore necessitate the use of DSA. The material properties for this case are shown in Table 5.
The spatial discretisation scheme used for the DSA acceleration of this problem is a modified interior penalty (MIP) (Wang and Ragusa, 2010). A major reason for choosing this scheme is that it provides a symmetric discontinuous spatial discretisation of the neutron diffusion equation so that computationally efficient preconditioned conjugate gradient (CG) solvers can be used to solve the resulting symmetric positive definite (SPD) system of equations. The conjugate gradient solver is preconditioned using multilevel algorithms based upon a continuous projection of the discontinuous spatial discretisation which have been shown to be effective for neutron diffusion problems (O'Malley et al., 2017b,c). These multilevel preconditioners make use of the algebraic multigrid algorithm AGMG for the low-level correction (Notay, 2010(Notay, , 2012(Notay, , 2014Napov and Notay, 2012).
The GMRES iterations required to solve the transport problem and solution time for the high scatter Kobayashi problem are recorded for various mesh refinements of structured hexahedral, tetrahedral, and pyramid elements using S 32 angular quadrature. The problem is run both with and without DSA and results are tabulated in Table 6. The high scatter Kobayashi is run in parallel in much the same way as for previous results, but the AGMG code used did not support OpenMP parallelism so the low-level correction of the diffusion preconditioner is run in serial. The impact of this on the timings shown would likely be minimal however.
The results provide strong evidence that the pyramid elements are performing well within the DSA framework. The iterations required for the neutron transport solver to complete are being reduced well with a small increase in time required per iteration and the performance of the pyramids matches well that of the hexahedral and tetrahedral elements.
For the lower spatial refinement cases the DSA is not reducing the iteration number by as much as might be expected. This is due to the fact that there are a low number of elements spanning the width of the duct and the penalisation terms at the borders of the duct are therefore leading to a poor diffusion approximation of the transport problem. If a discontinuous spatial discretisation method without any penalisation were to be used for the DSA, such as the Adams-Martin method (Adams and Martin, 1992), then the number of iterations would be significantly lower, particularly for the lower refinements. This would not be a stable method for some cases however so a penalised method such as MIP is necessary. Fig. 8. Visualisation of GODIVA quarter octant meshed with structured hexahedra. Some tetrahedra (green), wedges (yellow) and pyramids (red) must be included. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Fig. 9. Calculated K eff of GODIVA nuclear criticality benchmark problem for all cases as the number of degrees of freedom is increased.

Conclusions
This paper examined the computational performance of pyramid elements, with basis functions generated using the Bergot method, when solving discontinuous Galerkin discrete ordinate discretisations of the neutron transport equation. Through direct comparison with the performance of other element types, in particular tetrahedra and hexahedra, it was demonstrated that these elements perform well in terms of accuracy and convergence for a variety of nuclear reactor physics, radiation shielding, and nuclear criticality benchmarks.
For geometrically simple structured geometries hexahedral elements are generally considered to be the best performing for three dimensional cases. An L 2 error analysis of an MMS problem along with numerical simulations of the Takeda and Kobayashi benchmark problems, on full pyramid structured meshes, demonstrated that pyramid elements have a computational efficiency and accuracy relatively close to hexahedral elements.   . 11. Octant of a cube with a high refinement region exterior linked to a low refinement core using a structure of tetrahedra (green), pyramids (blue) and prisms (yellow). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Table 4 Criticality results for cubic problem with constant refinement and varying refinement linked with pyramids. Ele length refers to the dimensions of the high refinement hexahedral elements in the varied case.

First-order
Second-order Next, numerical simulations of the GODIVA nuclear criticality benchmark case comparing a basic mesh generated using unstructured tetrahedra with an extruded structured mesh which made use of pyramids to mesh a corner demonstrated a more realistic application of pyramid elements. By making use of some pyramid elements this hexahedrally dominated mesh is far simpler to generate than a fully hexahedral mesh would be for this geometry. However, the accuracy and solution times were clearly superior to the unstructured tetrahedral mesh with highly refined spatial meshes.
Finally a high scatter variation of the Kobayashi benchmark problem was used to test the pyramid elements for a more challenging neutron transport problem which made use of DSA. Overall the results presented in this paper demonstrate that the Bergot pyramid elements can be computationally competitive with hexahedral elements for a wide variety of typical neutron transport problems. Future work might focus on investigating the performance of the pyramids in less structured meshes, or when the pyramid elements are highly distorted.