COMPARING THE DISCONTINUOUS GALERKIN AND HIGH-ORDER DIAMOND DIFFERENCING METHODS FOR THE TRANSPORT EQUATION ON A LOZENGE-BASED HEXAGONAL GEOMETRY

This paper presents an implementation and a comparison of two spatial discretisation schemes over a hexagonal geometry for the two-dimensional discrete ordinates transport equation. The methods are a high-order Discontinuous Galerkin (DG) finite element scheme and a high-order Diamond Differencing (DD) scheme. The DG method has been, and is being, studied on the hexagonal geometry, also called a honeycomb mesh – but not the DD method. In this research effort, it was chosen to divide the hexagons into (at least) three lozenges. An affine transformation is then applied onto said lozenges to cast them into the reference quadrilaterals usually studied in finite elements. In practice, this effectively means that the equations used in Cartesian geometry have their terms and operators altered using the Jacobian matrix of the transformation. This was implemented in the discrete ordinates solver of the code DRAGON5. Two 2D benchmark problems were then used for the verification and validation, including one based on the Monju 3D reactor benchmark. It was found that the diamond-differencing scheme seemed better. It converged much faster towards the solution at comparable mesh refinements for firstorder expansion of the flux. Even if this difference was not present for second-order, DG was slower, about two to four times slower.


INTRODUCTION
There are several ways of dealing with this geometry. The hexagon can be kept whole and a Discontinuous Galerkin (DG) formulation using fairly exotic basis functions can be built upon it. Or the hexagon can be sub-divided into six triangular or three lozenge elements -which can then be further refined. The two methods have been respectively implemented in the French DG solvers MINARET and ERANOS [1]. However, neither triangular nor lozenge sub-elements of a hexagonal mesh have been investigated for DD. In addition, there does not exist a comparison of the two schemes on hexagons for eigenvalue problems. This paper addresses these two issues.
It is of note that this work has previously been done in Cartesian geometries [2], where it was shown that at low order, DD performs generally better than DG for comparable mesh refinement. Accuracy is comparable at second-order but DD is generally much faster. Beyond posterity's sake, it is quite interesting to compare these two methods since they appear to share quite a few similarities in terms of their formulation and implementation. Because of the upwind nature of the scheme, the problem can be solved locally, i.e., on each unit element -where the matrix to be solved is orders of magnitude smaller, albeit full. However, DG is less demanding on the regularity of the solution because of discontinuities at the boundaries.
For this work, it was chosen to divide the hexagons into lozenges. This brings about a certain elegance and ease of implementation. Indeed, an affine transformation is applied on the lozenges to cast them into reference quadrilaterals (usually squares), leading to a Cartesian reference geometry. In fact, this corresponds to transforming the terms and operators of the transport equation using the Jacobian matrix of the transformation. The fundamental DD and DG equations are not changed substantially, which therefore facilitates implementation. This work is being done in the SNT: solver of the DRAGON5 code -an open-source project from Polytechnique Montreal.

THEORETICAL BACKGROUND
Only the salient points of the theory behind this study will be detailed in this section. The highorder Diamond Difference (DD) and Discontinuous Galerkin (DG) finite element methods have both been covered at length in the literature. The authors invites the reader to peruse [1], [2], [3], and [4], and the references therein.

The Discrete Ordinates Transport Equation
In this angular discretisation method, the angular integral is approximated using a quadrature defined over the unit sphere, with points and weights respectively denoted by Ω n and ω n .
For each point or characteristic of neutron travel, Ω n , the one-speed stationary first-order form of the discrete ordinates transport equation is given by the set of linear hyperbolic equations, ∀r ∈ D, where D is the domain, with boundary Γ, defined in R 2 . Σ is the macroscopic total cross-section, Φ n (r) = Φ(r, Ω) is the neutron flux along that characteristic, and Q n (r) = Q(r, Ω) is the source consisting of fission and external sources as well neutrons scattering from other directions. The boundary condition is given by, ∀r ∈ Γ − , where Γ − is the inflow flow boundary of D, defined by Γ − = {r ∈ Γ : Ω n · n(r) < 0}, with n being the unit outward normal to Γ.
Considering a meshing M h of the domain D into elements K such that K∈M h K = D, the local variational form for the DG method can be expressed by on an element K, where the dependency over r has been made implicit, and the h subscript denotes the fact that these functions are obtained from a broken Sobolev space. The + superscript sign  indicates the restriction of the function taken from within the element K, while the − superscript sign the restriction of the function taken from the neighbouring element to K. Also, ∂K − is defined similarly to Γ − . It should also be noted that the starting point of the DD method is the same as this variational form minus the jump term.

Handling of the Hexagonal Geometry
As mentioned in the introduction, the hexagon is split into (at least) three lozenges. These individual lozenges can then be further sub-meshed, as shown in Fig. 1. These lozenges can then be mapped onto reference quadrilateral elements (denoted byK) by using a transformation F U , as shown in Fig. 2, where U denotes the specific lozenge. The associated Jacobian matrix is given Since the basis functions used in the DG and DD formulations are defined onK, Eq. (3) has to be solved onK. Hence, using the Jacobian matrix previously defined, the following equation is obtained where · is the Euclidean norm in R 2 and t is the unit tangent vector to the inflow boundary, ∂K − . Of particular note here is that while the Jacobian matrix is different for each of the 3 types of lozenges, the determinant and J K,U t are in fact the same for all of them. This greatly helps with the implementation.

NUMERICAL RESULTS
Using the equations obtained in the preceding section, the high-order Diamond Differencing (DD) method and the Discontinuous Galerkin (DG) finite element method were implemented in the DRAGON5 code [5]. The implementation was validated using a simple 2D one-group problem taken from [6]. A more complicated 2D reactor mock-up four-group (without up-scattering) test case was then computed.
The results presented in the following subsections are k ef f values and errors on the absorption rate  Figure 3: Domain of the 2D one-group benchmark, featuring a 30 • symmetry, the material mixtures, and vacuum boundary conditions. Hexagon side is 19 cm.  Table 1: Cross-section data for the 2D one-group benchmark shown in Fig. 3. Note the anisotropy in mixture 2.
distribution. The latter is calculated using where i denotes the cell element of calculation, and V i its volume. The maximum and average errors are then respectively worked out on these reaction rate values using the equations below: where the asterisk denotes the reference value, and V core is the volume of the whole domain.
Various parameters are investigated for the results, namely, the • angular discretisation order, N , which corresponds to N (N +2) 2 allowed directions in 2D.
• spatial discretisation order, with values of 0, 1, 2, and 3, respectively corresponding to flat, linear, parabolic and cubic expansions of the polynomials used to approximate the flux in the discretisation method. DD does not have cubic expansions, and DG does not have flat expansion.
Also, no DSA or GMRES(m) acceleration are available for the DG method as of the time of writing this paper. Hence, only a Livolant acceleration was used for the calculations for both methods. Finally, Level-Symmetric (LS) quadratures were used in all computations, and the power iteration is stopped when the maximum relative discrepancy between flux unknowns is less than 10 −6 .

One-Group Benchmark
The domain is described in Fig. 3, with the cross sections for the different material mixtures given in Table 1. It should be noted that the geometry is first unfolded in the solver before computation.
This benchmark is not at all characteristic of a real world reactor. However, transport and anisotropy effects are greatly emphasised over the whole domain and at the void boundary. This allows for a numerical verification and validation of the method, and its implementation. It was obtained from the paper [6], which studied the implementation of an SP n solver, and contrasted it with another discrete-ordinates solver, ERANOS/SNATCH -both used a lozenge submeshing of the geometry.
The asymptotic "reference" value for the ERANOS/SNATCH solver was 1.002313. This was obtained using a Gauss-Chebyshev quadrature allowing for 180 directions, a DG method with a flux Figure 4: Domain of the Monju-2D benchmark. Hexagon side is 6.67417 cm.  Table 3: Cross section data for the Monju-2D benchmark shown in Fig. 4. M. denotes the material Mixture, and G. the energy Group.
expansion of parabolic order using a Legendre hierarchical basis, and 27 lozenges per hexagons.
The results from the current implementation in the SNT: solver are printed in Table 2, and the obtained asymptotic "reference" value is 1.002337, a margin of 2.4 pcm. An S 18 also allows for 180 directions; therefore, the differences are the choice of basis and the order of expansion. Taking that into consideration, it can be said that convergence towards the same solution has been achievedand the solver validated.

The Monju-2D Reactor Benchmark
The Monju-2D reactor benchmark was formulated by selecting the second axial plane from the Monju 3D reactor benchmark. The latter is a three-group Fast Breeder Reactor (FBR) mockup proposed by Komano et al. [7]. The description of the domain considered in this paper is given in Fig. 4 and the material cross sections in Table 3. Mixtures 1, 2, 3 and 4 respectively correspond to the active inner core, active outer core, blanket and control rod. Also, again, while the domain features a 120 • symmetry, it is first unfolded before computation. The results are printed in Table 4.

Discussion of Results
From both Table 2 and Table 4, looking at each method independently, it can be seen that they perform mostly as expected with each varying parameter. There are a few inconsistencies. For example, for DD, when looking at low angular orders, especially for the linear order, the error on the k ef f can be seen to increase as the mesh is refined. This can be attributed to spatial and polynomial convergence being reached at that point, but angular order, N , still being the limiting factor. Indeed, as N is increased, this disappears. Comparing the two methods is more interesting. Even if DD was better in Cartesian, there was still some expectation among the authors that DG might be more precise and faster than DD on hexagons. However, the results indicate quite the contrary still. At comparable orders, for lower lozenge submeshes, DD does generally better, while DG is persistently slower -about 1.3× to 1.5× slower for the one-group benchmark and about 3× to 4× slower for Monju. This is consistent with what was observed for the Cartesian geometry [2], albeit to a stronger effect. While it could be ascribed to the source being slightly more complex to build for DG (because of the upwind sources), it is not only that. DG just requires more inner iterations to converge. This has not been illustrated here for lack of space, but it is consistently the case. More investigation is required into their behaviour with preconditioning and a GMRES(m) acceleration.

CONCLUSIONS
High-order Discontinuous Galerkin (DG) and high-order Diamond Differencing (DD) schemes for the discrete ordinates transport equation were implemented for hexagonal geometries. This was achieved by dividing the hexagons into lozenges, and transforming them into squares. While it had been done before for DG -but not for DD, the ease of implementation has been more formally shown. This had been lacking in previous literature. The methods were then compared on two eigenvalue problems. It was shown that at low order, DD has smaller margins of error at lower mesh-refinements. It is also consistently significantly faster. In the future, more work still needs to be done with acceleration techniques such as the DSA for better evaluation.