A Comparative Analysis of Neutron Transport Calculations Based on Variational Formulation and Finite Element Approaches

: The application of continuous and discontinuous approaches of the finite element method (FEM) to the neutron transport equation (NTE) has been investigated. A comparative algorithm for analyzing the capability of various types of numerical solutions to the NTE based on variational formulation and discontinuous finite element method (DFEM) has been developed. The developed module is coupled to the program discontinuous finite element method for neutron (DISFENT). Each variational principle (VP) is applied to an example with drastic changes in the distribution of neutron flux density, and the obtained results of the continuous and discontinuous finite element (DFE) have been compared. The comparison between the level of accuracy of each approach using new module of DISFENT program has been performed based on the fine mesh solutions of the multi-PN (MPN) approximation. The obtained results of conjoint principles (CPs) have been demonstrated to be very accurate in comparison to other VPs. The reduction in the number of required meshes for solving the problem is considered as the main advantage of this principle. Finally, the spatial additivity to the context of the spherical harmonics has been implemented to the CP, to avoid from computational error accumulation.


Introduction
The Boltzmann equation describes the local density of particles traveling inside a medium with interactions between the particles and the medium, involving the independent variables of position, time, energy, and direction of motion. Applying the well-known approaches to handle the spatial dependence of this equation is one of the central problems of nuclear reactor theory. Among these approaches, the finite element (FE), finite difference, and nodal methods are considered as the most efficient approaches [1][2][3][4][5]. The FE approach as a computer-based method can be applied for numerically solving a wide range of boundary value problems. In recent decades, there has been an increasing tendency in analyzing the neutronic aspects of the reactor, especially numerical solution of neutron transport equation (NTE), to employ techniques based on the families of discontinuous and continuous Galerkin finite element methods (FEMs). The discontinuous finite element (DFE) procedure makes use of the same function space as the continuous procedure, but with relaxed continuity at interface boundaries of each element or each region [6]. It was first introduced by Reed and Hill for the solution of the NTE [7], and recent developments have been reviewed by [8,9]. The essential idea of discontinuous finite element method (DFEM) is derived from the fact that the shape functions can be chosen so that either the field variable or its derivatives, or generally both, are considered discontinuous across the interfaces of elements, while the computational domain continuity is maintained. As a result, the DFEM includes, as its subset, both the FEM and the finite difference method. Therefore, it has the advantages of both the finite difference and the FEMs, in that it can be effectively used in generating various types of meshes with maintaining geometric flexibility and controlling the order of angular approximations through each element. This feature makes it uniquely useful for computational neutron transport calculations. Because of the local nature of a discontinuous formulation, unlike the continuous methods, no global matrix (GM) needs to be assembled; and thus, this reduces the demand on the memory. The effects of the boundary conditions on the interior field distributions then gradually propagate through element-by-element connection. This is another important feature that makes this method useful for neutron transport calculations.
In comparison with the other numerical methods (finite difference and FEs), the discontinuous FE plan has both advantages and disadvantages, which are important for developing specific applications [10,11].
The scalar flux often has a changeable and unstable treatment in different regions with respect to the physical conditions of the system. For example in cases, where the absorption cross-section has a high value, encountering very strong gradients in the scalar flux and discontinuity in angular flux (AF) at interfaces is predictable. This discontinuity mostly can be found in some problems with rapid changes in scalar flux. Therefore, the DFEM can improve these depletions, easily. In addition, the coupling between element variables is achieved through the boundary integrals only. This means that the source terms and stiffness matrixes are fully decoupled between elements, and some special techniques are applied to assembling the local matrixes. With the use of these techniques and considering no dependence between elements, one can find that the order of angular approximation is changeable for each element with respect to the necessity of the problem. Hence, a DFE algorithm will not result in an assembled GM, and thus, the in-core memory demand is not as strong. Furthermore employing the DFEM leads to decreasing the cost of the calculations, while the values of errors are decreasing. In addition, the local formulation makes it very easy and comfortable to parallelize the algorithm, taking advantage of either shared memory parallel computing or distributed parallel computing, which is not the purpose of this paper. Finally, because of the local formulation, both the ℎ-adaptive (mesh size optimization) and -adaptive (optimization the order of angular approximations) refinements are made convenient. The ℎ -adaptive algorithm (mesh and order of angular approximation optimization) based on the discontinuous formulation requires no additional cost associated with node renumbering compared with the continuous finite element (CFE) method. Various approaches have been adapted to interpret the concept of discontinuous and two widely accepted ones, which are based on variational methods in comparison to the CFE method are presented in this paper [12][13][14][15][16].
FEMs (both continuous and discontinuous) have been applied by various techniques to solve NTE. In some techniques, a weighted residual approach is used, while in others, the variational methods are used. In 1997, the classical maximum principle-based on a generalized least-square approach to obtain CFE solutions of the even-parity transport equation (TE) was presented by Ackroyd [17]. The application of variational approaches for solving the NTE has been investigated over recent years.
The numerical solutions to the NTE based on variational formulation and continuous/discontinuous FEM have been successfully applied to our previous investigations: In 2011 an adaptive CFE approach for NTE was applied by Abbasi et al. [18]. As the second study of our group, the spatially adaptive hp refinement approach for the PN NTE using the spectral element method based on the CFE method was carried out [19]. In 2018, an extended half-range spherical harmonics method for first-order NTE based on variational treatment and CFE was presented by Ghazaie et al. [20,21]. In 2019, Yousefi et al. [22] performed an investigation for CFE modeling of the 3D NTE based on the variational methods, and the computational program ENTRANS was developed. Finally, in 2018, Sadeghi et al. [23] applied an adaptive refinement for the PN NTE based on conjoint variational formulation, and DFEM and the new computer code, discontinuous finite element method for neutron transport (DISFENT), in two-dimensional geometry was developed. The main goal of this study is to present a comparative analysis of numerical neutron transport calculations using the variational formulation based on the continuous and discontinuous FEMs. This paper is organized as follows: The theoretical background of the NTE based on even and odd components is presented in the next section. Developing the variational principles (VPs) by considering the direction of motion and spatial dependence to NTE is analyzed in the third section. In Section 3, the maximization of the VP based on changing the penalty parameter is performed. Developing the formulation of the DFE with the element by element neutron conservation (NC) and conjoint maximization is presented in Section 3.
Finally, in Section 4, to compare the performance of both continuous and DFE approaches that are reviewed in the article, each VP is applied to the well-known Reed's edge cell [24] benchmark problem. In order to achieve this goal, the developed computer code DISFENT, in which the module of performing the calculation based on the CFE is also added, has been applied.

Theoretical Background
The full NTE is an integro-differential equation in seven independent variables (three in space, two in angle, one in energy, and one in time), and as a consequence, the number of degrees of freedom that would be required, for the accurate discretization of all the independent variables can grow extremely quickly. However, the monoenergetic approximation, called one-speed approximation, which is most widely considered neutron transport model both analytically and numerically, is described as [25]: where ( ⃗, Ω ⃗ ) is the distributed extension source of the neutron, ( ⃗) is the total cross-section, and ( ⃗, Ω ⃗ . Ω ⃗ ) is the differential scattering cross-section at ⃗ for scattering from the direction Ω ⃗ into the direction Ω ⃗ . Equation (1) is the first order form of NTE in its steady-state formulation. However, for generating the VPs, the second-order of TE is required [17]. Hence, the even and odd parity (EAOP) TE is derived.
Now, these parity fluxes and sources are substituted into Equation (1) to obtain the mixed parity equations: where the leakage and removal operators for the parities of an arbitrary function ⃗, Ω ⃗ are defined as: It has been demonstrated that both of the operators and and their inverse and are both self-adjoint and positive definite [16].
From Equation (10), the odd parity of AF can be obtained as: and finally, by substituting Equation (13) in to Equation (9), the second-order equation for the evenparity flux ⃗, Ω ⃗ can be written as: − This equation is called the even-parity TE. In a similar way, the odd-parity TE can be obtained as:

VP Based on the FE Approach
Traditionally, a partial differential equation is said to have a VP if there is a function for which the vanishing of the first variation yields the partial differential equation. It does not directly provide VPs for the first-order Boltzmann equations for steady-state and time-dependent transport. Thus, the second-order even-parity equation and the odd-parity equation can be used to construct VPs. The extremum principles for both the steady-state and time-dependent can be obtained from a generalized least-square method. A feature of classical VPs is that admissible TFs have to satisfy continuity at the interfaces of regions. However, this can be wasteful in the use of elements for those places where the AF is expected to be changing rapidly with , because a fine spatial mesh of elements has to be used in such a situation with a classical principle. However, if utilizing the discontinuities in the TF was permissible at the interfaces of elements, an acceptable approximate solution with fewer elements would be possible.
For generating a VP based on DFE a volume of system with various types of the boundary such as bare, perfect reflector, or albedo is considered. If in the identity (8), the approximation value is used instead of the exact value , then one can write: Now by substituting this definition into the pair Equations (9) and (10) the following consequence yields: where , are the residuals that provide a measure of the errors made by a trial function in approximating the exact solution . For obtaining a VP, which is based on continuous TF, two measures of error are made by a trial function : Failure to satisfy the TE in volume of the system and is denoted by ( ).

II.
Failure to satisfy the boundary condition on boundary , which is specified as ( ).
However, as regards that the purpose is to achieve a principle that admits discontinuous TF, hence another error should be added too. The nature of this error is based on failing to satisfy the interface continuity condition between two elements having volumes and . According to Figure  1, when two elements are considered separated and it is assumed that there is no connection between them, then the error is appeared. Hence, an overall error measure is obtained by adding the three error measures and is given by: with definitions: and where, with employing the even and odd parts of AF, the error reduced to: Now, by collecting the introduced errors, the functional ( ), which is a generalized leastsquares measure of the amounts by which the trial function fails to satisfy the first-order Boltzmann equation and the interface and boundary condition, can be obtained.
Using the generalized least-square method and with respect to the special properties of the operator's leakage and removal (self-adjoint and definite positive) and the divergence theorem, the important maximum principle can be derived as In addition, the formulation of classical principle has been derived in [18], which is based on only the volume component of error and the boundary error functional and again using the generalized least-square method and algebraic operators, the maximum principles that admit the continuous TFs in space and angular variables, can be written as: where and and the complementary of this principle is: where and Now, through comprehensive observation of Equation (25) and making a comparison between the first three terms of Equations (27) and (28), the maximum principle [ ] can be reduced to Here, ⃗, Ω ⃗ − ⃗, Ω ⃗ is the jump in the TF at ⃗ on the interface ( ∩ ), and is a positive number called penalty parameter and it has the role of controlling the jump in the approximate solution at interfaces.

The Directional Dependence of the AF
The spherical harmonics, or , is one of the most popular methods for describing the directional dependence of AF, which is used in this study. They can produce a high-quality approximation for the angular dependence even when a relatively small number of basis functions are employed, particularly when the solution is closed to isotropic. Thus, the method can be used to accurately model a wide variety of neutron transport problems.
In one dimensional geometry, the angular dependence of the EAOP flux is expressed in terms of even and odd order Legendre polynomials, respectively. Thus, the EAOP flux can be approximated by the TFs with definitions: and also, the even parity TF for two dimensional can be written as where The series expansions of the AF into angular basis functions can be written in matrix notation as: and where the transposed vectors and represent the EAOP angular basis functions, respectively, for corresponding one-and two-dimensional geometries.

Finite Element TFs for Spatial Domain
There are some approaches for verifying spatial discretization of the TE. However, in this paper, the spatial dependence of the even-parity AF is presented in terms of the nodal values of the flux and the FE shape functions. For this purpose, the even parity flux for each element can thus be approximated by the trial functions ( , Ω) as where is the number of nodes in each element, which can vary from region to region. The Equation (39) can be written in compact notations as: and with a similar procedure, the odd parity flux appears as Again, the matrix form of this expression is given by

The Maximization of the VP
The maximum principles, given by [ ] and [ ] functionals, which admit continuous TFs, have been derived in [26]. The TFs used in these principles are forced to obey the continuity conditions at the interfaces between elements. However, in the DFE approach, there is another situation. The maximum principle, given by [ ] functional, admits both continuous TFs and discontinuous TFs. Unlike the continuous case, the discontinuous TFs do not have to obey the continuity conditions at interfaces between elements. Hence, generally, there are two different approaches for implementation of the FEM on the TF. According to the first approach, all the elements are completely connected with the adjacent elements. It means that the continuity condition between the adjacent elements is established. Hence, this approach is used for the classical principle, which admits continuous TF. However, the main objective of this paper is based on analyzing the various discontinuous aspects of the FEM. In the discontinuous approach, two types of discontinuities can be considered between adjacent regions or adjacent elements of each region, which are called region-wise technique and elementwise technique, respectively. In the first case, a continuous flux profile is assumed to exist within the interior of the regions of the system, while discontinuities are allowed at the interfaces between these regions. In the second case, a totally discontinuous mesh is considered where the TFs are allowed to be discontinuous at all interfaces between adjacent elements. These two approaches for twodimensional geometries are shown schematically in Figure 2. In global terms, the reduced functional of [ ] can be written as where the interfacial functional is expressed as In addition, is the number of interfaces around each element, which is illustrated for various types of elements in Figure 3. By inserting the even parity TF into Equation (43), the local stiffness matrix for each element and an interface matrix for each interface can be obtained. After performing the assembling procedure, the following discretized functional holds: where the local stiffness matrix , obtained in this case, is the same as was obtained in principle.
Requiring [ ] to be stationary with respect to variations of each component of [ ] yields: where the GM + is symmetric, positive-definite, and because of the local compact support, also sparse and banded.
The essential difference from the calculations viewpoint between principle and principle is to form the GM. In the continuous case, all the local matrixes are assembled as ordinary, which was introduced in Cassiano et al. [27]. However, in the principle, which is based on the DFE, the stiffness matrixes of interfacial elements are located next to each other without any assembling, and then, the interfacial matrix couple local matrixes. According to this, different orders of angular approximations can be considered in various regions.

Influence of the Penalty Parameter
A simple treatment is given here for the influence of the penalty parameter on a discontinuous variational approximation for the AF, which is given by the discontinuous variational approximations for the even and odd parity AFs. In fact, the penalty parameter is used to control the discontinuity range of elements. As a result of this, when the penalty parameter is increased, the jump at the interfaces is reduced. This remark suggests that there may be some specific values of for which a DFE solution is superior to a classical continuous variational solution [17].
For specifying the penalty parameter in applying [ ] , some appropriate cases are recommended. Generally, by using a high value (physical infinite) for , the solutions tend to continuous solutions. However, in multi-region problems, the value of is defined with respect to the properties of each region, and hence = , where h is element width, can be considered as a choice for the penalty parameter. Another possible choice for the penalty parameter is = 0.

DFE with the Element by Element NC
The maximum principle based on the DFE can be imposed with only a single discontinuous trial function (TF). However, in the sequel a general maximum principle ( , , , ) can be considered with the penalty parameter non-negative, so that the degree of control on discontinuities is varied from none for = 0 to their extinction as → ∞. In this approach of DFE, two TFs are used. The TFs and are continuous across interfaces, and they couple the elements into each other. In addition, and are the helping TFs to ensure element by element NC and, within an element, are functions of neutron direction only. In this section, the new maximum principle based on DFE and element by element NC is derived.

Specification of TFs
Over the whole system, the expression used for classical components ± ⃗, Ω ⃗ of the discontinuous TFs ± ⃗, Ω ⃗ are of the kind ± ⃗, Ω ⃗ = ± ± ± ⃗, Ω ⃗ = 1,2, … for an FE spatial mesh with M nodes. In the following, a general representation for the discontinuous components ± of the TF ± ( , Ω) can be written as where outside of element the basis functions ± ⃗, Ω ⃗ vanish.

New Formulated of [ ]
The [ ] maximum principle for the first order Boltzmann equation, which is given by [16] can be expressed as where the bilinear functional ( , ) for arbitrary ⃗, Ω ⃗ and ⃗, Ω ⃗ is defined as: Equation (48) is derived from the identity where the definition of the error functional is and similarly

Conjoint Maximization of ( + ) and ( + )
With > 0, the interfacial term in the ( ) functional ensures that the vanishing of ( ) leads to = , the exact solution. However, with choosing the zero value for the penalty parameter, a different situation arises. It can be demonstrated that the vanishing of functional ( + ) does not necessarily imply + = [15]. Hence, in constructing the VP, the following property is used: if, and only if, = and = 0 = .
This property and the identity which is the consequence of Equations (58) and (59), provides the motivation for the conjoint maximum principle: The new maximum principle ( , , ) can be simplified by letting the discontinuous component of the odd parity flux vanish so that = 0, and hence: For a one-dimensional slab geometry case, the discontinuous components and can be expanded as (66) Hence, the optimized solution for in an interior element, is obtained by maximizing principle with respect to the coefficients of , which yields < , > = − ( , ) + ( , ) where the angular direction (μ = cosθ) comprises the polar angle between the z-axis and the unit vector Ω ⃗ in the direction of the neutron motion.
With the particular choice of unity for , Equation (67) reduces to the conservation condition for an element, which states the supply rate of neutrons entering , from the interior source and the inward flow across the surface of the element, which is balanced by the absorption rate for and the leakage rate of neutrons outwards across . Therefore, the conservation conditions for an element, which is imposed as desirable physical constraints on a TF can arise also as the consequences of maximizing principle. By taking Legendre polynomials as the basis functions of the directional dependence of the continuous and discontinuous components for the EAOP flux and using the orthogonality property of the Legendre polynomials in the interval [−1,+1], the various terms of Equation (67) after tedious algebraic calculations can be determined and eventually the reason is expressed as where the continuous components ± are optimized by means of the [ ] and [ ] principles respectively. Since Equation (68) can be solved for each element independently of the other elements, the unknown coefficients can be obtained very efficiently. Finally, the average scalar flux is approximated by summation of the various components of the TF, which is written as In addition, Equation (67) is developed for two-dimensional problems and simple square elements, which the result in an element at the bare surface boundary giving:

Computational Results
In the previous sections, it was showed that the various approaches of the FEM with arbitrary boundary conditions can be applied to VPs. In addition, a new variational scheme based on the element-by-element NC in which neutron AF is approximated by its decomposition into continuous and discontinuous EAOP components was described. The comprehensive developed computer code based on continuous and DFEM in one dimensional and − geometry using multiple VPs is used to generate an approximate solution for AF for any order approximation. To check the capabilities of various schemes, the solutions of the scalar flux resulting from each VP versus the exact transport solution are compared. An analytic solution for the Boltzmann equation is available only for some limiting cases. In order to evaluate the obtained results of the explained methods, solutions of the multi-PN (MPN) approximation are considered as the reference solutions [20,21]. The multi-PN (MPN) approximation to the first-order NTE is based on applying a variational approach. The considered reference solution in this study uses the 16P15 with a uniform spatial mesh of 2000 linear elements. The angular domain in the interval of [−1,+1] is divided into 16 folds and, in each fold, the order of Legendre polynomial is 15. The most important purpose is to demonstrate the advantages and disadvantages of different types of VPs.
To illustrate the effects of continuous and DFE approaches on the accuracy of numerical solutions, each VP is applied to the well-known Reed's edge cell benchmark problem. This example is formed of five regions and diverse properties are considered for each one of them. With respect to the material properties of the edge cell problem, which is given in Table 1, it is predictable to be encountered with very strong gradients at interfaces in the scalar flux and consequently, sharp changes in the distribution of neutron flux are extremely expected. Hence, this problem can be considered as an excellent example to compare the capability of each method. The geometry of the problem is given in Figure 4.    Tables 2-6, the  relative errors for , approximations with a set of elements are indicated. The number of elements in each region initially started by coarse meshes and the ascending procedure for the number of elements is considered. The applied boundary conditions in calculations are based on the formulas of (67) and (21). The graphical representation of obtained errors for approximation is presented in Figure 5A-E. The most striking results to emerge from all these figures are that: Table 2. Relative errors for , approximations with a set of elements for region 1.  a. These figures and tables analyze and compare the obtained results from various methods of the FE approach. According to this comparison, it is evident that the conjoint principle (CP) appears to behave better than the methods using CFEs or principle. In addition, as follows from these plots, the obtained errors from and principles almost overlap. As mentioned earlier, by increasing the penalty parameter, the results tend to continuous solutions. b. The graph indicates that one of the great advantages of the conjoint approach is that by using the coarse meshes in every region yields the solutions with more agreement in comparison to other approaches. However, in some regions, the obtained results from the CP at first tend to zero, but by increasing the number of elements, the calculated errors are raised. In fact, it is surprising to find that the errors are decreased and again are increased. There is a general justification for describing this situation . As mentioned previously, the CP is based on the NC in each element. Hence, in some cases, by implying some coarse elements, NC is satisfied and increasing the number of elements not only cannot improve the accuracy but also decreases, because of the accumulation of computational errors. This difficulty is partly compensated by employing the spatial additivity to the context of the spherical harmonics which mainly is based on two common error estimation algorithms. The first algorithm is applied by using some information about the exact solution, which is called prior estimator [28]. The second approach is based on a posterior error estimator [29][30][31] and is developed by a particle balance relationship over an element, which is derived by integrating the Boltzmann equation over an element and can be written as:

Order of Angular
Imposing the parity parts of the trial function φ (consisting of the continuous and discontinuous components), the balance equation reduced to and finally, the required conservation condition is expressed as by introducing = Ω the source in the element e = [ + ] Ω removal rate in the element = Ω. Ω , net current leakage from the element e also, is an assessment criterion which is evaluated the difference between approximate and reference solutions. According to expression (72), when the residual is greater than the specified tolerance, the spatial element is refined by dividing an element into two sub-elements. By implementing this technique, only the elements, which cannot satisfy the NC will be divided into two smaller parts. The main advantage of this algorithm is to prevent increasing the number of elements for transport calculation during refining process (in comparison to the case that all the elements should be divided into smaller parts). Figure 6 indicates how the residual in each region is affected by increasing the number of elements in both uniform and adaptive case. According to this curve yielding, the residual approaches a value of zero at low number of elements in the uniform case, but as decreasing the size of the elements uniformly, the NC will be far away from balance situation, while with respect to the adaptive refinement, this difficulty can be solved as is shown in Figure 6b. c. The obtained results showed the efficiency of various principles in spatial analysis. However, for the visual representation of the ability of VPs about angular dependence, the reader is referred to Table 7. In this table the total errors of the entire domain for various angular approximations, which are evaluated by a fine mesh and the solution, are summarized. The most conspicuous observation emerges from the data comparison is that the CP again has lower errors in comparison to other approaches.

Conclusions
This paper presents the comparative analysis of numerical neutron transport calculations based on variational formulation and DFE method. The formulations of three main approaches based on the continuous and DFE methods are developed. In order to perform a comparison between the various methods, some extra modules have been attached to the first version of the program DISFENT [21]. The main focus of this paper has been the demonstration of the ability of each described method for modeling a complex example and comparing the accuracy of the obtained results. It has been shown that the developed calculator program DISFENT can accurately model various examples both using continuous and DFE approaches.
The obtained results show that the CP appears to behave better than the methods using continuous FEs or principle. In addition, the CP leads to more accurate results in comparison to other principles, when the coarse meshes are applied to all the suggested principles. The main point in using the CP is that by implying coarse elements, the NC is satisfied in some elements and increasing the number of elements just leads to error accumulation and oscillation in results. Hence, the spatial adaptivity to the context of the spherical harmonics has been employed in order to avoid the stated difficulty.