DAPHNE-3D: A NEW TRANSPORT SOLVER FOR UNSTRUCTURED TETRAHEDRAL MESHES

A new Discrete Ordinates transport solver for unstructured tetrahedral meshes is presented. The solver uses the Discontinuous Galërkin Finite Element Method with linear or quadratic expansion of the flux within each cell. The solution of the one-group problem is obtained with non-preconditioned fixed-point or GMRES iterations. Precision and performances of the solver are evaluated on the 3D Radiation Transport Benchmark Problems proposed by Kobayashi, showing very good agreement with the reference and good computing times in serial execution.


INTRODUCTION
The Service d'Études des Réacteurs et de Mathématiques Appliquées (SERMA) is currently developing a computing platform for radiation shielding applications, called OPERA. The aim is to provide an integrated environment containing tools for data pre-and post-processing, particle transport and material activation. OPERA already contains simplified transport solvers based on point-kernel, a 3D discrete ordinates solver for Cartesian geometries, and a reference Monte Carlo transport code. To extend OPERA capabilities, we are currently developing a new solver, designed to treat arbitrary 3D geometries. The new solver is called DAPHNE (for Deep Attenuation PHotons and NEutrons) and it solves the first-order form of the transport equation using the multigroup approximation for the energy and the discrete ordinates approximation for the angle. In its first implementation, the geometry representation is based on a tetrahedral mesh partitioning the space.
Solvers with similar characteristics are not new in the transport community. Among the codes that solve the first-order form of the transport equation we can quote Split-Cell [1], ATTILA [2], AHOT-C-UG [3] and AETIUS [4]. The first one uses a linear or exponential expansion of the flux within tetrahedral cells, and a short-characteristics closure. In a similar fashion, AHOT-C-UG uses the short-characteristics formulation but with a robust arbitrary polynomial expansion of the flux within each cell. The ATTILA code uses the Discontinuous Galërkin Finite Element Method (DGFEM) with linear expansion of the flux within tetrahedral and a trilinear expansion in hexahedral elements. The AETIUS code mimics ATTILA in that it uses linear DGFEM for tetrahedral meshes. The main difference between the characteristic formulation and DGFEM is in the choice of the closure relating boundary fluxes: in short-characteristics, this relation is retrieved from the integral form of the transport equation, while DGFEM uses the "upwind" approximation. As a consequence, the characteristics formulation better approximates, in principle, problems with pure streaming with the disadvantage that the coupling coefficients involve the evaluation of computationally expensive exponential functions of the total cross-section. On the contrary, the DGFEM coefficients have a simple linear dependence on the total cross-section but pure transmission is not explicitly accounted for.
In the first implementation of DAPHNE-3D, we chose to explore the DGFEM formulation. In addition to the linear DGFEM that can be found in ATTILA and AETIUS, we also implemented a quadratic expansion of the flux within each tetrahedral cell. The objective is to explore whether a higher order formulation helps to achieve similar precision with reduced computational resources.
The paper is organized as follows: in next section we will briefly recall the DGFEM formulation as proposed by Wareing [2] and provide some algorithmic details on the construction of DGFEM matrices and on the solution of the one-group equation. The following section is dedicated to the results obtained on the benchmarks proposed by Kobayashi [5].

DGFEM Formulation
For the derivation of the DGFEM equations, we start from the first-order form of the one-group, discrete ordinates transport equation assuming the standard Legendre expansion of the scattering cross-sections: with ψ is the angular flux, q is the total angular source, A l is the real Spherical Harmonic of order l, and φ l is the angular moment of the flux with respect to the Spherical Harmonic of order l, computed with an angular quadrature formula. We have also assumed that the external source, Q l , shares the same Spherical Harmonics expansion of the flux. Boundary conditions (BCs) are associated to Eq. (1). In our implementation we allow inhomogeneous BCs and reflective BCs. The latter are expressed as the mapping between the outgoing flux for a given direction and the incoming flux for the reflected one. The derivation of DGFEM equations is done following the formulation proposed by Wareing et al. [2]. For simplicity of the notation, we will omit the angular dependence and write the equations for one single angle. We start assuming that the spatial domain of the problem is partitioned into an unstructured mesh composed of K tetrahedral cells of volume V k and surfaces Γ kf for k = 1, . . . , K and f = 1, . . . , 4. For each cell, we define a basis θ k (r) and write the expansion of the angular flux: A similar expansion is assumed for the angular moments of the flux, the external source and the angular source. In our implementation we use basis of the Lagrangian type for linear (4-nodes) and quadratic elements (10-nodes). By projecting Eq. (1) onto the basis function, and by making use of the divergence theorem we obtain: Here, r i is either x, y or z and Ω i is the projection of Ω onto these axes. The symbol n kf represents the outgoing normal of face f of element k. Notice that for linear tetrahedral elements, this is constant for each face. The boundary fluxes, ψ kf , are provided by the upwind closure: where with the symbol k − we indicate the neighbor region through the face f . For the external boundaries, the incoming flux is given by BCs.
Eq. (3) is written in the global spatial coordinate system. As shown in [2], the calculation of DGFEM coefficients can be efficiently done in the local coordinate system of a unit tetrahedron . By applying this transformation, the volume integrals in Eq. (3) become: with J k the Jacobian matrix of the transformation. Matrix M and the three K j can be computed once and used for all elements. As for volume integrals, the integral on surfaces in Eq. (3) is done using a local coordinate system mapping any triangular surface to the unit triangle of surfacẽ The integration is done by noticing that the trace of the Lagrangian basis function onto the surfaces is identically zero for the components that do not belong to the surface. We can thus write: with P kf the rectangular permutation matrix N v × N s , with N v the number of volume nodes (4 for linear, 10 for quadratic elements) and N s the number of surface nodes (3 for linear and 6 for quadratic). The permutation matrix is written for the face f and element k to take into account the local numbering of nodes of each region. The vectorθ s contains the N s non-zero entries of the trace of the basis. By inserting Eq. (4) and Eq. (7) into Eq. (3) we can finally obtain the surface integral as: with |J f | being the magnitude of the surface Jacobian of the face f . Also note that we have distinguished the right-and left-permutation matrices to account for the cases of "outgoing" and "incoming" fluxes due to the upwind closure. In practice, since local numbering is known within each cell, the U f = P kf U P T kf is pre-computed, while for the incoming flux, the permutation is done on-the-fly using pre-computed U . The final form of Eq. (3) is: For each region and angle, we build the matrix on the l.h.s of the equation on-the-fly and solve the problem using the Eigen library [6]. This strategy favors the reduction of the memory since the matrices do not have to be stored, but it can be computationally intensive especially for higherorder elements.

Void Limit
Some discretization schemes fail to provide the solution of transport in void regions. This is for example the case of schemes based on the second-order form of the transport equation [7], or when using cell-balance in the characteristics formulation [8]. The solution of Eq. (9) in void regions requires the matrix on the l.h.s. to be invertible for σ t = 0. The analysis of U f and K j alone reveals that these matrices are singular for both linear and quadratic elements. However, it is easy to show that their kernels are orthogonal, thus insuring that the matrix resulting from their sum is always invertible. It must be noticed that for cells with high aspect ratio, the resulting matrix can be ill-conditioned.

One-Group Solution Algorithm
The solution of Eq. (1) is done by "sweeping" regions downstream for all the angles independently. This procedure does not allow an exact inversion of the one-group equation because of the presence of the scattering source and the reflective BCs. To solve Eq. (1) we implemented non-preconditioned fixed-point iterations and GMRES as proposed by Warsa et al. [9]. No modifications are done to their algorithm so details will not be given. The only difference worth noticing concerns the surface angular fluxes included in the Krylov space: while in [9] they include incoming and outgoing fluxes for reflective surfaces only, we chose to include outgoing fluxes for all the surfaces. This allows obtaining the outgoing flux for every surface without any additional sweep. If needed, incoming fluxes at reflective boundaries may be obtained by applying BCs.

Benchmark Results
The aim of the paper is to test the precision and the effectiveness of our DGFEM implementation for shielding applications regardless of the errors arising from the multi-group cross-section generation. To that end, we choose to test the solver on the Kobayashi Benchmark Problems [5]. They consist of three regions (source, void, and shield) with total cross section σ t = 0.1 cm −1 for source and shield, and σ t = 1E−4 cm −1 for the void region. Three geometry configurations are proposed, each one having two cases: Case i with no scattering, and Case ii with scattering ratio σ s /σ t = 0.5 for all regions. The benchmark problems are modelled and meshed using the GEOM and SMESH modules of SALOME-8 [10] respecting the initial geometry and by applying reflective boundary conditions.  fixed-point iterations, and 1E−7 residual in GMRES. Several tests have been performed for different mesh sizes and orders of angular quadrature aiming to the identification of the numericallyconverged solution. Table 1 summarizes the chosen mesh and the angular quadratures, as well as the Root Mean Square (RMS), memory usage and computing time for a serial execution. The RMS here is computed on the calculation points required by the benchmark. Detailed results are given in Tables 2 and 3. Punctual values are obtained reconstructing the scalar flux in the desired point using the spatial basis function of the DGFEM approximation.
As a general comment, we see that both schemes have small relative error even for fluxes with attenuation of several orders of magnitude. Notice that the solution for Case ii (50% scattering) converges with lower order of quadrature as compared to Case i. This is due to reduced ray-effect in problems with scattering. Fig. 1 shows the RMS for the three tests for varying number of degrees of freedom. We notice that in tests 1 and 3, for comparable RMS, less degrees of freedom are required by quadratic elements with consequent reduced memory usage. This behavior is not observed in Test 2. However, in terms of computing time, the reduced number of degrees of freedom does not imply faster solution, as it can be seen from Table 1. As an example, Figure 2 shows a 3D plot of the mesh and of the scalar flux obtained for Case 3ii with linear DGFEM.
A comparison of our results with those published by AETIUS in [4] shows comparable precision. However, AETIUS results are computed using much coarser angular quadrature formula. This is probably due to the lack of the First Collision Source treatment in our solver. On the other hand, AETIUS solutions use much more refined spatial meshes. A comparison of our results with those provided by Ferrer with AHOT-C-UG [8] cannot be done since they model void regions using σ t = 0, leading to higher fluxes. In addition, they provide results only for Case ii and for a subset of the points required by the benchmark. We can only notice that the convergence of the spatial mesh seems to be achieved with equal or coarser meshes. No information is provided about computing times in both references. Comparisons with the ATTILA code are not possible since results about Kobayashi benchmark have not been published.

Conclusions and Perspectives
DAPHNE-3D is a new discrete ordinates transport solver for three-dimensional unstructured tetrahedral meshes. In this work we have implemented a Discontinuous Galërkin Finite Element Method for space discretization using linear and quadratic expansions of the angular flux. Inhomogeneous and reflective boundary conditions are currently available. The one-group solution algorithm uses non-preconditioned fixed-point iterations or GMRES. The precision and the performance of the DGFEM solver have been tested on the 3D Radiation Transport Benchmark Problems proposed by Kobayashi for shielding applications. Results show very good agreement between the  Comparisons between linear and quadratic elements show that the higher-order expansion helps reducing the total number of degrees of freedom, with consequent reduction of the memory. However, the use of higher-order elements does not always provide a reasonable reduction of computing time. Results also confirmed that a very large quadrature order is required to eliminate the ray-effect arising in this kind of configuration.
As a perspective here, we are interested in the investigation and comparisons of additional spatial discretization methods such as those proposed by [1] and [3], in particular for shielding applications. Future work will be also dedicated to the construction of suitable low-order preconditioner of the one-group equation for accelerating the fixed-point and GMRES iterative schemes. The implementation of ray-effect mitigation techniques such as the First Collision Source method is also a necessary step for the effective solution of shielding problems. Finally, since the solver is currently implemented using only serial execution, significant reduction of the computing time can be obtained by applying shared-memory parallelism.