Arbitrary-Order Bernstein–Bézier Functions for DGFEM Transport on 3D Polygonal Grids

In this paper, we present an arbitrary-order discontinuous Galerkin finite element discretization of the SN transport equation on 3D extruded polygonal prisms. Basis functions are formed by the tensor product of 2D polygonal Bernstein–Bézier functions and 1D Lagrange polynomials. For a polynomial degree p, these functions span {xayb}(a+b)≤p⊗{zc}c∈(0,p) with a dimension of np(p+1)+(p+1)(p−1)(p−2)/2 on an extruded n-gon. Numerical tests confirm that the functions capture exactly monomial solutions, achieve expected convergence rates, and provide full resolution in the thick diffusion limit.


Introduction
The linear Boltzmann transport equation, which we will simply call the transport equation, describes the transport of neutral particles (i.e., neutrons and photons) through background and multiplying media [1]. Some of the most prevalent discretization schemes are the discrete ordinates (S N ) in angle and multigroup in energy. Given an angular quadrature set, {w m , Ω m } M m=1 , the multigroup, S N transport equation for the angular flux of group g traveling in direction Ω m , denoted as ψ m,g ( r) ≡ ψ g ( r, Ω m ), with isotropic scattering within an open, convex spatial domain D, with boundary, ∂D, is written as Ω m · ∇ψ m,g ( r) + σ t,g ( r)ψ m,g ( r) = 1 4π where the source term is given by νσ f ,g ( r)φ g ( r). Eigenvalue (2) σ t,g is the total cross section of group g, σ g →g s is the scattering cross section from group g to group g, νσ f ,g is the fission neutron production term of group g, χ g is the neutron fission spectrum of group g, S g is a distributed source, and the scalar flux of group g, is given by The discontinuous Galerkin finite element method (DGFEM) is a popular spatial discretization for Equation (1) [2][3][4]. Besides 1D, the use of DGFEM for transport has typically been reserved for meshes employing simplicial or tensor mesh cells. Stone and Adams [5] developed a discretization with linear precision for unstructured polygonal cells, which Bailey [6] later extended to unstructured polyhedral cells. Hackemack and Ragusa [7] analyzed quadratic serendipity functions on polygonal grids, which was further extended to extruded polygonal prisms [8]. Recently, utilizing the framework of Floater and Lai [9], Hackemack [10,11] analyzed arbitrary-order DGFEM S N and diffusion solutions on 2D polygonal grids using reduced-space polygonal Bernstein-Bézier functions. This work extends [10] to form arbitrary-order tensor-based functions on extruded polygonal prisms.
The remainder of this paper is as follows. In Section 2, we present the DGFEM discretization of the S N transport equation, followed by the employed 3D extruded basis functions in Section 3. Then, Section 4 presents numerical results demonstrating the expected precision and convergence properties of our DGFEM discretization. Finally, conclusions are presented in Section 5.

The DGFEM Transport Equation
In this section, we review the DGFEM discretization of the S N transport equation. First, we lay down an unstructured mesh, T h ∈ R 3 , over the spatial domain consisting of nonoverlapping spatial elements to form a complete union over the spatial domain: D = K∈T h K. The angular flux unknown of Equation (1) is expanded using discontinuous functions over each mesh cell K: where N K is the number of basis functions on element K, and the functions are nonzero over a single element and discontinuous across the element interfaces. The scalar flux and source term can be expanded analogously. Multiplying Equation (1) by a test function u, integrating over element K, applying the divergence theorem, and using the standard upwind technique yields the following weak form over element K: ( Ω m · n) u, ψ m,g ∂K + − Ω m · ∇u, ψ m,g K + σ t,g u, ψ m,g In Equation (5), ∂K ± represents the outflow/inflow boundary, ∂K ± = r ∈ ∂K | Ω m · n K ( r) ≷ 0 ; with n K ( r) the outward unit normal, ψ ↑ m are the angular flux values on an inflow face taken from the upwind neighbor element, and the inner products correspond to integrations over the cell volume and faces, respectively, where dr ∈ R 3 is within the cell and ds ∈ R 2 is along the cell boundary. The error of the discretized flux solution φ hp ∈ W hp D with polynomial degree p, under the L 2 -norm, can be written as where h is the maximum diameter of all mesh elements, r is the regularity of the transport solution, and C is a constant independent of the mesh. Therefore, the convergence rates are limited by the solution regularity [3].

3D Extruded Polygonal Basis Functions
In Section 2, we gave the DGFEM discretization of the S N transport equation but no definition to the test and trial functions. Now, we provide some details on the employed polygonal prism basis functions. Our extruded polygonal functions are formed from the tensor product of the 2D polygonal Bernstein-Bézier functions (we will refer to them as the (x, y) radial functions) and the 1D Lagrange polynomials (we will refer to them as the z axial functions).

Brief Overview of 2D Polygonal Functions
We now provide a brief overview of the radial polygonal functions. In the interest of brevity, full details are not provided, and the interested reader should see the previous works of Floater and Lai [9] and Hackemack [10]. On a polygonal element, the reducedspace polygonal Bernstein-Bézier functions, { B p }, are formed by a set of functions with support points on vertices, along edges, and within the interior. They reduce to the Lagrange fundamental polynomials on triangles. For a n-gon and a finite element space of degree p, there are n vertex functions, p − 1 functions along each edge, and (p−1)(p−2) 2 interior functions for a total of pn + (p−1)(p−2) 2 functions. They span the {x a y b } (a+b)≤p space of functions. Figure 1 gives the distribution of the degrees of freedom on an example pentagon for polynomial degrees p = 1, 2, 3, 4.

Extrusion to Polygonal Prisms
With the radial functions given, we can now provide details on converting them into an arbitrary-order functional space compatible with extruded polygonal prisms. For polynomial degree p, the extruded 3D functions, { E p }, are formed from the tensor product of the radial functions { B p } and the (p + 1) 1D Lagrange polynomials (denoted as { L p }): For a polygonal prism with lower and upper bounds of z B and z T , respectively, the axial functions have evenly-spaced support points given by and the i-th axial function is given by Since these extruded functions are taken as the tensor product of our radial and axial functions, their functional space is given by span( E p ) = {x a y b } ⊗ {z c }, where dim( E p ) = np(p + 1) + (p + 1)(p − 1)(p − 2)/2 (n still represents the number of vertices of the 2D polygon). The support points of the degrees of freedom are given in Figure 2 for an example pentagonal prism for polynomial degrees p = 1, 2, 3. An analogous p = 4 case compared to Figure 1 was not provided due to the excessive cluttering of the nodes.

Basis Function Verification
We first provide numerical experiments that verify properties of our extruded basis functions and the corresponding DGFEM transport solutions. We test three properties: (1) verify the DGFEM captures the proper interpolation space, (2) verify the DGFEM achieves proper convergence rates in h-type and p-type refinement, and (3) verify that the solutions maintain full resolution in the thick diffusion limit (TDL) [12].
We first test the interpolation properties and convergence rates of the S N solutions using the method of manufactured solutions (MMS). Our domain is [0, 1] 3 , with σ t = 1, σ s = 0, S 8 quadrature, and all vacuum boundaries. We define two functional forms to test interpolation, ψ p , and convergence, ψ s , properties given by x a y b z c , ψ s (x, y, z) = sin(k x πx) sin(k y πy) sin(k z πz), (11) with k x = 1, k y = 2, and k z = 3. Three extruded mesh types are utilized: Cartesian (hexahedra), triangular, and polygonal. The 2D polygonal meshes were generated with the PolyMesher software [13] before extrusion. Table 1 gives the L 2 -norm of the error for DGFEM orders of p = 1, ..., 6. From the machine precision results, we can clearly see that the transport solutions fully-live within the functional space of the polygonal Bernstein-Bézier functions. The loss of exact precision is due to the magnitude of the functions to grow past unity with increasing polynomial degree. Next, Figures 3 and 4 provide the convergence of the sinusoidal solution, ψ s , under uniform h-refinement and p-refinement, respectively. Since N do f ∝ h −3 in 3D, the slopes in Figure 3 approach −(p + 1)/3 as expected. Additionally, since ψ s is analytic within r ∈ D, the exponential convergence of p-refinement in Figure 4 is expected.       Having verified the interpolation and convergence properties of our DGFEM scheme, we now confirm that our DGFEM transport solutions retain full resolution in the thick diffusion limit. From the work of Adams [12], it is known that a poor discretization of the S N transport equation (with full upwinding) leads to 'locking' into absurd transport solutions. Adams states that two properties of the DGFEM functions are needed for full resolution: (1) locality of the functions and (2) surface-matching (see Adams work for further details). The extruded polygonal prism functions of this work do possess these properties, and as such we expect full resolution.
In the TDL, the domain mean free path approaches infinity. If we fix the physical dimensions of the problem, we can then scale the cross sections and source to reflect the properties of the TDL. Introducing a scaling parameter, , we scale σ t as 1/ and σ a and Q as . This leads to a scaled transport equation and corresponding diffusion equation given by respectively. The diffusion equation is discretized using a conforming Galerkin finite element method (CGFEM) with strongly-enforced homogeneous Dirichlet boundary con-ditions. As → 0, the discretized transport solution, φ T , will approach the discretized diffusion solution, φ D , at a rate of under the L 2 -norm (||φ D − φ T || L 2 ). This is confirmed in Figure 5.

Multigroup Pincell Problem
Our final numerical example is a 7 group pincell problem (with reflecting radial boundary conditions) that utilizes the geometry and cross sections of the UO 2 pincells from the C5G7 benchmark problem [14]. The fuel region is formed from a 16-sided regular polygon where the fine mesh cells are formed by square tiling until the pin boundary is intersected. The convergence history of the error in k (given in units of pcm) are given in Figure 6. Superior convergence is realized for increasing polynomial degree but with diminishing returns. A significant reduction in error is realized going from linear to quadratic FEM.

Conclusions
In this paper, we have presented an arbitrary-order DGFEM discretization of the S N transport equation compatible with extruded polygonal grids. The basis functions are formed by a tensor product of 2D radial (xy) polygonal Bernstein-Bézier functions and 1D axial (z) Lagrange polynomials. On an extruded n-gon with polynomial degree p, there are np + (p − 1)(p − 2)/2 radial functions and p + 1 axial functions, with their tensor product resulting in np(p + 1) + (p + 1)(p − 1)(p − 2)/2 total functions. These functions can exactly interpolate the {x a y b } ⊗ {z c } space of functions, where a + b ≤ p and c = 0, . . . , p. Due to the locality and surface matching properties of the functions across edges, these functions allow for transport solutions to maintain full resolution in the thick diffusion limit. Numerical testing also demonstrated the functions' interpolation properties and appropriate convergence rates were observed. The methodology presented in this work can thus be used to compute S N solutions on 3D extruded polygonal prisms (e.g., reactor analysis calculations).
Funding: This research received no external funding.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.