An Extended Linear Discontinuous Method for One-group Fixed Source Discrete Ordinates Problems with Isotropic Scattering in Slab Geometry

Nowadays, the obtainment of an accurate numerical solution of fixed source discrete ordinates problems is relevant in many areas of engineering and science. In this work, we extend the hybrid Finite Element Spectral Green’s Function method (FEM-SGF), originally developed to solve eigenvalue diffusion problems, for fixed source problems using as a mathematical model, the discrete ordinates formulation in one energy group with isotropic scattering in slab geometry. This new method, Extended Linear Discontinuous Discrete Ordinates (ELD-SN), is based on the use of neutron balance equations and the construction of a hybrid auxiliary equation. This auxiliary equation combines a linear discontinuous approximation and spectral parameters to approximate the neutron angular flux inside the cell. Numerical results for benchmark problems are presented to illustrate the accuracy and computational performance of our methodology. ELD-SN method is free from spatial truncation errors in S2 quadrature, and generate good results in the other quadrature sets. This method is more accurate than the conventional Diamond Difference (DD) and Linear Discontinuous (LD) methods, but surpassed by the Spectral Green ́s Function (SGF) method, for quadrature order greater than two.


INTRODUCTION
The neutron population calculations in non-multiplying medium is relevant in different areas such as nuclear physics and technology, mathematics, radiation protection and material sciences.
The neutron transport equation [20] is the mathematical model used to estimate the neutron density related quantities. As neutrons are not charged particles, they cause ionization through a complicated mechanism involving energetic secondary particles emission. For that reason, neutrons are often called "indirectly ionizing particles". They may be quite penetrating and the shielding required may be massive and expensive [17]. Therefore, accurate solutions for neutron transport problems are crucial in nuclear science.
Analytic solutions for neutron transport equation are limited to a few particular simple cases, involving simplified geometries and rough approximations. This is why numerical approaches are frequently used. There are two main approaches in computational methods to solve neutron transport equation: the probabilistic and the deterministic. In this work, we use the deterministic approach based on the discrete ordinates formulation (S N ) [8] as mathematical model. The computational requirements used by S N calculations are expensive in memory and CPU time. For that reason, researchers have dedicated much effort to obtain accurate numerical solutions in coarse spatial meshes.
An alternative to performs coarse mesh neutron transport calculations are the Spectral Green's Functions (SGF) nodal methods. The first spectral nodal method was proposed in [5] to solve fixed source one-dimensional S N problems. Since then, several authors extended the SGF methods to different geometries and formulations for the transport equation, both for fixed source [6,12,13] and eigenvalue problems [11,3,4]. Among the advantages of SGF methods are the high numerical precision and weak dependence on the spatial mesh. However, it has a high algebraic outlay to obtain the sweep equations and also a high computational cost. Recently, some authors developed numerical approaches that exploit the advantages of SGF methods, minimizing its deficiencies. Among these, stand out the composite spatial grid methods [14] and the finite element spectral nodal hybrid methods [24].
In neutron transport problems with azimuthal symmetry, the one-dimensional Cartesian geometry offers good results. For that reason, many authors used this geometry in recent works [21,10,19,25]. In this paper, we develop the hybrid Extended Linear-Discontinuous method, for discrete ordinates formulation (ELD-SN), to solve fixed source problems in slab geometry, considering one energy group approximation and an isotropic scattering. This method extends the prime work in hybrid methods, proposed by [24] to solve eigenvalue diffusion problems in one-dimensional Cartesian geometry. The ELD-SN combines a linear discontinuous approximation inside the spatial cell, common in finite element method, and a quasi-analytic approach, characteristic in SGF methods for preserving the analytical spectral solutions inside the cell. This new method offers accurate numerical results and minimizes the drawbacks of SGF methods. We validate this numerical formulation by solving benchmark problems.
As we know, there are several efficient methodologies to solve fixed source problems in S N formulation for slab geometry [28,26,5]. Nevertheless, as there are no previous researches from the use of hybrid spectral nodal methods in neutron transport formulations, the main goals of this work are to propose a new numerical approach and to evaluate its performance in the solution of neutron transport problems.
In sections 2, 3 and 4 we present the numerical formulation of ELD-SN method. First, we obtain the spatial balance equations, then we describe the spectral analysis, and finally we present the auxiliary and sweep equations. Section 4, offers numerical results for model problems in homogeneous and heterogeneous domains. At the end, we conclude our work and give some suggestions for future work.

SPATIAL BALANCE EQUATIONS
In order to develop the proposed ELD-SN method, we consider as a mathematical model, the one-speed S N formulation in one-dimensional Cartesian geometry with isotropic scattering [20] in a generic domain D of length L is the angular flux of particles traveling in the discrete ordinates direction µ m , σ t (x) is the total cross section, σ s0 (x) is the zeroth order component of the differential scattering cross section, ω n is the Gauss-Legendre angular quadrature weight for each direction, Q(x) is the fixed external source, f m and g m are prescribed functions, and N is the number of discrete ordinates.
In order to obtain a solution for the one-speed S N equations, we divide the domain into I spatial cells of length h i = x i+1/2 − x i−1/2 as represented in figure 1, and solve eq. (2.1) in each of these cells. The cross sections and the external sources are piecewise constant functions inside the cells. The first step is to obtain the S N spatial balance equations, then we can use the operator where P l is the Legendre polynomial of l order, and x i is the midpoint of spatial cell. Considering l = 0 in operator (2.3) and applying in eq. (2.1), we obtain the zeroth-order balance equation The first-order balance equation is obtained by considering l = 1 In spatial balance equations (2.4 and 2.5) we define the zero'th and the first moment of the angular flux and the scattering source as average quantities in the spatial cell 9) and the auxiliary constants (2.10) Several numerical methods have been reported in the literature to solve one-energy group, one-dimensional, fixed source discrete ordinates problems with isotropic scattering. Lewis and Miller (1984) describes the traditional Diamond Difference (DD) method conforming a marching scheme that follow the direction of neutron travel using the zeroth-order balance equation (2.4) together with the Diamond Difference relation as auxiliary equation [20]. Also, Hill (1974) and Barros (1990) studied the Linear Discontinuous (LD) method that forms an iterative linear system, using the zeroth and first-order balance equations (2.4), (2.5) and two auxiliary sweep equations [15,2]. The main difference between the DD and LD methods is the solution representation inside each cell. The DD method aproximates the solution by continuous linear functions across the spatial grid. The LD method approximates the solution by discontinuous linear function across the spatial grid. Therefore, the DD method has one degree of freedon, while the LD method has two.
These two methods implement a marching scheme, yielding an iterative algorithm based on the source term, called Source-Iteration (SI) scheme [20]. This scheme can be used because the auxiliary equations are decoupled in the angular directions. This approximation is weak, because the contribution of other directions to the average flux is neglected. Nevertheless, the sweeping equations are simpler. This SI scheme consists of three stages: first, a sweep from left to right for µ m > 0 is made, then, a similar sweep is made, now, from right to left for µ m < 0 and finally, the stopping criterion is verified.
More recently Barros and Larsen (1990) proposed the Spectral Green´s Function (SGF) method describing a convergent numerical scheme wich generates numerical solutions completely free of spatial truncation errors in slab geometry. This is because, it preserves the analytic general solution [5]. This method uses zeroth-order balance equation (2.4) and an auxiliary equation, that computes the average angular flux inside the cell, considering the contribution of all incoming discrete directions. The authors solve the SGF equations using the one-Node Block Inversion (NBI) iterative scheme. This scheme uses the most recent available estimates for the node-edge angular incident fluxes on a given node, to calculate the exiting angular fluxes in the upwind directions [6]. So, the NBI sweep involves more difficult algebraic manipulation and needs higher computational cost per iteration than standard SI sweeps.
The discretization process and the spatial balance equations are common to all numerical approaches to solve S N formulation, the differences between them, appear in the auxiliary equations. In the next section we developed auxiliary equations for the ELD-SN method. The proposed method is compared with DD, LN and SGF methods in section 5.

SPECTRAL ANALYSIS
In order to develop the Extended Linear Discontinuous method for discrete ordinate formulation, it is essential to discuss spectral theory for the analytic solution of the S N equations. The spectral theory for neutron transport problems began in the late 60's with Zweifel and Case's work to solve analytically the one-speed transport equation for homogeneous infinite medium [9]. Later, several authors extended these mathematical fundamentals to other formulations and geometries [5,4,23].
First, let's consider that the general form of local analytic solution of eq. (2.1) inside one given spatial cell can be written as where ψ p m,i (x) and ψ h m,i (x) are respectively, the homogeneous and particular components of the solution. For constant external source, the non-homogeneous solution is x-independent, and appears as In order to obtain an expression for the homogeneous component, we consider the ansatz The eigenvectors are normalized using the quadrature parameters in the form By substituting eq. (3.3) into the homogeneous problem corresponding to eq. (2.1) and considering normalization restriction, eq. (3.4), we obtain the solution for the m-th eigenvector component .
Now, in order to calculate the eigenvalue associated with each eigenvector, we substitute eq. (3.5) into the normalization condition eq. (3.4). We obtain the dispersion relation When the ratio between scattering cross section and total cross section eq. (2.10), satisfies 0 < C 0,i < 1, the roots of the eq. (3.6) are all simple and real numbers and also, are the eigenvalues of our system. In symmetric angular quadrature sets, these N-roots appear in pairs, with opposite signs. Finally, considering the homogeneous [eq. (3.
The main idea of the present ELD-SN method is to modify the linear discontinuous auxiliary equations for preserving the analytic solution of angular flux [eq. (3.7)] inside the spatial cell. For that reason, two spectral parameters are added to auxiliary equations. These parameters are used to keep the characteristics of the spectral problem. The auxiliary equations of ELD-SN method are presented in the next section.

AUXILIARY EQUATIONS AND ITERATIVE SCHEME
The first development in the hybrid nodal methods, combining linear discontinuous finite element approximations and spectral nodal methods was introduced by Rocha et. al. (2016) [24]. Specifically, that work described the Finite Element Spectral Green's Function method (FEM-SGF) to solve the neutron diffusion equation in non-multiplying media for slab geometry. The present work extends the FEM-SGF method, for one-dimensional S N fixed source problem, in one energy-group, following a similar methodology used by Larsen (1986) [18] to improve the standard DD method.
For ELD-SN method we proposed auxiliary equations as a linear approximation involving, zeroth and first order moments of the angular flux. Then the auxiliary equations appear as where A and B are constant, called spectral parameters, wich are to adjust the angular flux inside the spatial cell to preserve the general analytic solution (3.7). At this point, we proceed to calculate the spectral parameters. The zeroth-order moment of angular flux is defined in eq. (2.6).
Substituting it into eq. (3.7), and solving the integral, we obtain Using a similar process for the first-order moment given in eq. (2.8), the result is In eq. (4.1) we consider the analytical solution (3.7) for the left term, and substitute eqs. (4.2 and 4.3) into the right side. After some algebraic manipulations, we obtain For each k-term in the summation, eq. (4.4) is (4.5) Equation (4.5) must be satisfied at the spatial cell edges. So, evaluating it in the left edge (x i−1/2 ) and in the right edge (x i+1/2 ), and after some manipulations we obtain a linear algebraic system of two equations in the form (4.7) Solving this system, we obtain expressions for the spectral parameters A and B: (4.9) As it can be seen in eqs. (4.8 and 4.9), A and B depends on the eigenvalues ν k . So, the basic idea in the ELD-SN method is that the higher eigenvalue is independent of the spatial cell width (h i ). To attend this constraint, we fix k = 1 in eqs. (4.8 and 4.9), and we substitute the spectral parameters in eq (4.1) thus obtaining (4.10) The decision to preserve the dominant eigenvalue guarantees that all eigenvalues does not tend to infinity or become complex for any h i value. In this sense, the proposed method presents a weak dependence to spatial mesh dimensions if compared with another numerical approaches. Also, it is important to note that numerical solution is completely free of spatial truncation errors for S 2 angular quadrature order. This is because, in this case, we have two eigenvalues, ν 1 = −ν 2 , then in S 2 problems all eigenvalues are preserved.
At this point, we obtain the iterative sweep equations of the source iteration scheme, for the ELD-SN method. Note that, in our spatial mesh (see Fig. 1), the set of zeroth and first-order balance eqs. (2.4 and 2.5) with boundary conditions eqs. (2.2), form an undetermined linear algebraic system with 3NI + N unknowns (ψ m,i±1/2 , ψ m,i , ψ m,i ) and 2NI + N equations. In order to obtain a full determined linear system, we need to add NI auxiliary equations (4.1) for each sweep. In the µ m > 0 sweep, when we compute the emergent angular fluxes in the right spatial cell edge, eq. (4.10) is evaluated at x i+1/2 . For the µ m < 0 sweep, the computed emergent angular fluxes appear in the left spatial cell edge, and eq. (4.10) is evaluated at x i−1/2 . The result for each sweep is considering the upper-sign for the right-to-left sweep, and the lower-sign for the left-to-right sweep. To simplify the notation in eq. (4.11) we use the spectral parameters A and B from eqs.
. Finally, we get the iterative sweep equation for first-order moment of angular flux using eq. (4.11) in the first-order balance equation (2.5). That is, . In these equations the upper-sign is used for the right-to-left sweep (µ m > 0), and the lower-sign for the left-to-right sweep (µ m < 0). Using an approximation for scattering source we compute the exiting fluxes in the right/left sweeps for each arbitrary spatial cell i, and also, for each angular direction m, by following these three steps: 1. In eqs. 3. Finally, after finish the right and left sweeps, the stopping criterion is checked to finish the iterative process or the scattering source moments eqs. (2.7 and 2.9) are updated for the next iteration.
In order to quantify the computational cost of each method, we use the definition made by De Barros (1990) [1] for the flop (floating point operation). In the DD method, the amount of work required by the SI-scheme for one spatial cell is approximately 4(N +2) flops and for LD method, is about 8 (N + 3). This flop count only considers the floating-point operations necessary to evaluate the exiting cell-edge angular flux in one direction, including the scattering source calculation as the transport sweeping is performed. We notice, that as the quadrature order N increases, the amount of work required by the LD method per direction and per iteration is approximately two times greater than the amount of work required by the DD method [1].
De Barros (1990) also reported that the amount of work required by the SGF method with the NBI scheme is 2[N(N + 2) + 4] flops. It is important to notice that this flop count does not include the preliminary calculations necessary in this method, for example, calculation of the eigenvalues and the node-average Green´s function matrices.
On the other hand, the amount of computational work demanded by the proposed ELD-SN method with the SI scheme is the same that the LD method, which is 8 (N + 3) flops. So, the computational cost per direction and per iteration is significantly less than the amount of work required by the SGF method.

NUMERICAL RESULTS AND DISCUSSION
In this section, we solved a numerical benchmark to illustrate the accuracy and performance of the developed ELD-SN hybrid method. Then, we compared this method with the DD, LN and SGF methods on a heterogeneous domain.
The benchmark problem was adapted from Barros (1990) [7]. It is a heterogeneous slab 0 ≤ x ≤ 100 divided into three regions without external source (Q = 0). The boundary conditions, geometry and material regions parameters are shown in figure 2. We use four quadrature sets (S 2 , S 8 , S 16 and S 32 ) and different spatial meshes to obtain the numerical solutions. In tables 1 and 2 we show the scalar flux at the region interfaces, the iteration number, the maximum deviation (δ ), execution time (CPU) and the relative cost efficiency (κ). To calculate the relative deviation we used as the reference solution, the results of SGF method for a fine mesh (2000 × 5000 × 3000 cells). The relative efficiency was proposed by [5] as a ratio between CPU time consumed by each method and CPU time used to generate SGF comparable results (differences less than 5% at any mesh point) on the same computer. The convergence number to stop the iterative process was 1E − 08. In Table 1, we can observe that SFG and ELD-SN methods are completely free of spatial truncation errors for S 2 quadrature, this means that the numerical results do not depend on the spatial mesh. In this quadrature order the spectrum has only two elements, and these are the eigenvalues and their eigenfunctions, which are preserved for all h i values. In addition, analyzing tables 1 and 2 for higher quadrature orders, we see that the DD, LD and ELD-SN results converge to the reference solution when the spatial mesh becomes finer. However, the SGF results are free of spatial truncation errors. The proposed method is more accurate than the DD and LD methods but is not free from spatial truncation error. This is because the two higher eigenvalues are preserved for all h i but the N − 2 remaining eigenvalues are not. Similar behavior has been observed in Larsen (1986), who implemented the Extended Diamond Difference method [18]. We remark that for higher order quadrature calculations the ELD-SN generates more accurate results than the DD and the LD methods, with respect to the SGF results. We observe that regarding the CPU time reported in tables 1 and 2, the proposed ELD-SN generates accurate results using less CPU time than the other methods. This behavior is confirmed by its relative efficiency parameter that reports values lower than 1 for coarser meshes and higher quadratures. If we consider computational performance for one-sweep iteration (See section 4), ELD-SN has a similar computational cost as DD and LD methods, also, it is cheaper than the SGF method. Besides, it can be noticed in these tables, that the SI scheme requires more transport sweeps than the NBI scheme, the same behavior reported by De Barros and Larsen (1990) [6] when they analyzed these two iteration schemes. If we use convergence acceleration methods as it is described in [19,16], the iteration sweep number in SI scheme would be significantly reduced, consequently, the computational performance of the ELD-SN method could be improved.
Finally, the ELD-SN method shows better computational performance that the SGF method for higher order quadratures on coarse mesh calculations, considering the same spatial mesh.

CONCLUSIONS AND FUTURE WORKS
In this work, we present the ELD-LN method to solve fixed source problems of neutron transport using, as the mathematical model, the one-group discrete ordinates formulation in slab geometry with isotropic scattering. The fundamental idea behind this method is to partially preserve the analytic general solution inside the spatial cell by introducing two spectral parameters in the linear auxiliary equations. Thus, this method can combine the spectral nodal methods' accuracy with the linear approximation's simplicity.
Analysing the results, we observe that for S 2 quadrature the ELD-SN method is free of spatial truncation errors, however, for other quadratures is not. Even so, the method is more accurate than the DD and LD methods and also, generates good results in problems with strong flux gradients. If compared with the SGF method, the ELD-SN method generates the same numerical results for S 2 quadrature, and for the other quadratures, the results are very close.
The proposed method, based on the SI sweep scheme, is simpler than the SGF method. Also, it presents a satisfactory computational performance, surpassing the SGF method for the same spatial grid. Nevertheless, the ELD-SN method is not free of spatial truncation error, as is in the SGF method. The mesh independence of the SGF allows it to generate accurate solutions with only a spatial cell per region. For that reason, in slab geometries, SGF offers a better computational performance than the ELD-SN method.
Knowing that there is several analytical, semi-analytical and numerical methodologies to solve discrete ordinates problems in slab geometry [28,26,5,27,22], the main contribution of this work is to explore the hybrid spectral nodal methods in the solution of S N formulation using an elementary problem. So, we establish the bases to extend this simpler approach to more complex scenarios as multi-group problems, anisotropic mediums, and multidimensional geometries, where its simplicity and computational performance improvement would be relevant.
In the multi-group approach, we will consider slab geometries and in the multidimensional case, the X,Y -Cartesian geometry. In both cases, the spectrum cardinality of the S N equations increases. As the ELD-SN method just preserves two spectral elements, the method approximation is rougher. We can improve the approximation by increasing the finite element polynomial order. That is, extending the method to cubic discontinuous (CD); fifth-degree discontinuous (5D), etc. The CD method will be able to preserve four elements of the kernel, the 5D should preserve six; and so forth. To extend ELD-SN method to anisotropic problems we should introduce the corresponding terms in the scattering source and reproduce the mathematical procedure described in sections 2, 3 and 4.

ACKNOWLEDGMENTS
Authors acknowledge to Fundação de Amparoà Pesquisa do Estado da Bahia (FAPESB, BRAZIL) for the partial financial support provided.