Using Generalized Basis for Functional Expansion

Functional expansion has been rigorously studied as a promising method in stochastic neutron transport and multi-physics coupling. It is a method to represent data specified on a desired domain as an expansion of basis set in a continuous manner. For convenience, the basis set for functional expansion is typically chosen to be orthogonal. In cylindrical PWR pin-cell simulations, the orthogonal Zernike polynomials have been used. The main advantage of using functional expansion in nuclear modeling is that it requires less memory to represent temperature and nuclide variations in fuel then using a fine discretization. Fewer variables are involved in the data storage and transfer process. Each nuclide can have its unique expansion order, which becomes very important for depletion problems. In a recent study, performance analysis was conducted on Zernike-based FETs on a 2D PWR geometry. For reaction rates like the absorption rate of U-238, however, many orders are needed with Zernike-based FETs to achieve a reasonable accuracy. This gap inspires the study in this paper on alternative basis set that can better capture the steep gradient with fewer orders. In this paper, a generalized functional expansion method is established. The basis set can be an arbitrary series of independent functions. To capture the self-shielding effect of U-238 absorption rate, an exponential basis set is chosen. The results show that the expansion order utilizing exponential basis can reduce by half of that from using orthogonal Zernike polynomials while achieving the same accuracy. The integrated reaction rate is also demonstrated to be preserved. This paper also shows that the generalized functional expansion could be a heuristic method for further investigation on continuous depletion problems.


Introduction
High fidelity multiphysics simulations are highly desirable in nuclear engineering applications to improve understanding of next-generation reactors and reduce conservatism in existing reactors. Most of the multiphysics codes are developed individually and coupled together through an external data transfer scheme. These schemes are generally in the form of specific engineered scripts to link one application to another. For example, a lot of Monte Carlo neutron transport codes use the Computational Solid Geometry (CSG) representation of the physical geometry, while thermal hydraulics codes and fuel performance codes usually utilize unstructured mesh to implement finite element type of solvers. To accurately represent data at all points on a domain, physic-specific discretization is needed, which greatly increases computational costs. Creating a unified fine mesh, as many have tried, makes an already very costly problem even more expensive. A coupling method utilizing Functional Expansion Tallies (FETs) to transfer multiphysics data from Monte Carlo codes to other multiphysics tools has been rigorously studied [1][2][3]. This transfer scheme is based on prior work focused on reconstructing fluxes and currents in nuclear systems [4,5]. The FET scheme has been successfully implemented and tested in the MOOSE framework [1,2] and in the SERPENT code [3].
In previous FET implementations, Zernike polynomials are chosen to be the basis set used on 2D cylindrical geometry. In a recent work by the authors, a thorough examination of the Zernike-based FETs against conventional histogram tallies is performed on an LWR pin cell [6]. Comparison was made to identify the optimal expansion order needed for accurate representation of desired reactions rates. The main advantages of using FETs in multiphysics simulation is the memory reduction over finely discretized tallies, the simpler geometrical representation and faster tracking in cell-based codes. Fewer variables are involved in the data transfer process since each nuclide can have its unique expansion order. This feature could potentially reduce the huge cost of on-node storage in nuclide depletion. However, to capture highly varying radial distributions such as that of U-238 absorption, more than 16 expansion terms are needed with Zernike Polynomials whereas in a 10-ring discrete tally, only 10 variables need to be stored [6].
To solve this problem, an alternative basis set that can better capture the steep gradient with fewer orders is desired. In this paper, a generalized functional expansion theory is introduced for arbitrary basis sets. A set of exponential functions is proposed as a better basis for U-238 self-shielding effect. A comparison is made between the new exponentialbased FET and Zernike-based FET for a 2D PWR pin cell problem.

Function Expansion in Nuclear Modeling
The key idea for FET is to represent discrete data as an expansion of continuous function of polynomials [5]. In all previous work, the basis function is chosen to be orthogonal polynomials for their reconstruction simplicity. In the following equations, ψ n (x) is the n th term in an orthogonal basis set ψ(x) and ρ(x) is the weighting function. On the bounded domain Γ, our target distribution can be represented as follows: whereā n is the true n th expansion coefficient defined by: a n = Γ ψ n (x)ρ(x)P(x)dx (2) and k n = 1 ||ψ(x)|| 2 is the normalization factor. In a Monte Carlo simulation,ā n is estimated by an unbiased estimator ofā n as, We can then reconstruct P(x) based onā n and ψ(x).

Generalized Functional Expansion
For certain reaction rates, like the U-238 absorption rate, the distribution can have very steep gradients due to spatial self-shielding effects. Orthogonal polynomials prove to be a poor choice to capture highly varying distributions [6]. A generalized functional expansion theory is needed for arbitrary alternative basis sets.
Assume the generalized basis set as ψ. Each basis ψ i for i = 0 · · · N does not necessarily need to be orthogonal to each other. However, each basis has to be linearly independent. The quantity of interest is P(x) = ∑ m k=0 a k ψ k (x) for a maximum of m th order expansion. For simplicity, the normalization factor is absorbed in the coefficient a k , and weighting function ρ(x) is assumed to be unity. The Monte Carlo tallying works as an integral on domain Γ: where t i is the tallied value for i th order, j is the tallied event in a total of N events. To correlate tallied values with the expansion coefficients a k , we have the following expression: Note that Γ ψ i (x)ψ k (x)dΓ is the inner product < ψ k , ψ i >. Therefore, the system of equations is given as: The matrix form is Gā =t. G is the matrix for all inner products.ā andt are arrays of true expansion coefficients and arrays of tallied values respectively. Matrix G, by definition, is a Gramian matrix. Since each basis is linearly independent, the Gramian matrix is invertible [7]. The expansion coefficients can then be calculated asā = G −1t .

Conservation of Reaction Rates
For the orthogonal basis set, the matrix G shown in Section 2.2 is a diagonal matrix. Each coefficient is independent of others. Adding up one more term will not affect the lower order terms. However, if the matrix is dense when the basis is not orthogonal, enlarging G will lead to a different solution of coefficient vector a, thus all terms are needed for a faithful reconstruction.
In the original functional expansion theory, the zeroth order term preserves the reaction rates, since zeroth order of the basis set is essentially unity. Since each basis is orthogonal to each other, the total reaction rate is always preserved. For an arbitrary basis set, this property is necessarily not guaranteed. Caution should thus be paid to this issue, since conservation of reaction rate is crucial in nuclear reactor applications.
Suppose the distribution can be represented by an m th order expansion: P(x) = ∑ m k=0 a k ψ k (x). To conserve the reaction rate, the 0 th order basis ψ 0 (x) needs to be 1, which will be imposed in the rest of this study.

Choosing New Basis
The change of basis set can actually be regarded as a change of variable. In conventional functional expansion, we tally the polynomials of true position r to reconstruct a highly varying distribution f (r). Instead of tallying r, the new variable y = h(r) is tallied. With a selected function h, we hope that the reconstructed function f (y) presents a smoother behavior and that y = h(r) is easy to evaluate.
An easy and intuitive guess is to use an exponential function, e.g., y = exp(r C ). By tuning the parameter C, we can change the shape of f (y) such that it can be accurately represented by low order expansion. To determine the optimal constant C, we use the PolynomialFeatures feature in Python Scikit-learn package for polynomial regression on f (y) [8]. The ill-behaved distribution f (r) is obtained by running a 100-ring discrete tally of U-238 on 2D BEAVRS fuel pin model. By varying different C, the smallest order n of regression reaching a R 2 score of 0.995 is recorded in Table 1. Examples of polynomial regression at different C are shown in Figure 1.
From the results shown in Table 1, when C increases, the minimum order to have a good fit with R 2 ≥ 0.995 decreases. When C = 5 and 6, only a 5th order expansion is needed. When C gets even larger, the shape can always be captured by a fifth order polynomial. For a normalized r ∈ [0, 1], the range of y is [1, e]. In Figure 1, the information near y = 1 is very dense because of the exponential transformation. The inner fuel region (circled in the Figure 1b) could not be well represented in this least square polynomial regression. To balance all the effects, C = 6 is chosen to be an optimal guess to capture the U-238 self-shielding effect

Setups
The test problem is to examine the feasibility of the exponential basis function on a 2D PWR pin cell and to compare the performance with conventional Zernike-based FETs. The model is a 2D fuel pin taken from the BEAVRS benchmark [9]. Every fuel rod has a radius of 0.392 cm, with a fresh fuel of 4.5% enrichment.
The Monte Carlo code OpenMC is used as the neutron transport solver in this paper [10]. The new exponential expansion filter has been implemented in OpenMC.
Only the radial dependency is considered in this problem to facilitate comparison with common analysis models. For each pin cell simulation, a total of 10 6 neutron histories are used. For FETs, the tallied results are the coefficients of each basis, which is defined over the whole 2D fuel pin and results are then integrated over the equal volume rings for comparison. To calculate uncertainties on the reconstructed values from FETs, 10 independent tests are conducted to obtain the true variance of the estimates.

Results and Discussion
Comparison with Zernike-Based FET From Section 2.4, the new functional expansion basis at order n refers to ψ n = (exp (r 6 )) n = exp (nr 6 ). A direct approach to compare Exponential-based and Zernikebased FET is to compare the relative difference of FETs against histogram estimator tallies in a given ring. As with the discrete benchmark, the tracklength-based estimator is chosen which typically provides better statistics in a radially discretized fuel pin [11]. The relative difference between tracklength-based estimator and collision-based estimator is the target accuracy for this study.
Again, the relative difference at a certain ring is given below, and the standard deviation of this relative difference is calculated using the following error propagation formula.
where f stands for the reaction rate at the given ring. Because of the self-shielding effect, the reaction rates at the outer most ring is theoretically larger, and the slope of change is larger. For a larger first order derivative, it generally requires higher order of expansions to capture the change in shape. Therefore, the outer most ring is chosen for comparison since it represents a tougher test for the FET tallies. Figure 2 is the comparison between Exponential and Zernike-based FETs. Figure 2a shows the relative difference between the FET results and the tracklength results in the outer most ring of the 10-ring model. The error bars are too small to be visible, which means the fluctuation on the plots cannot be explained by stochastic behavior. The dips are from the change in sign inside the absolute values given in Equations (8) and (9). This figure clearly shows that the relative difference improves continually as expansion order increases. For Zernike-based FET, even at 50th order (26 terms for radial-only Zernike Polynomials in even order mode), the relative difference is still larger than the target standard (the green line in Figure 2a). The relative difference converges to the magnitude of 10 −3 after 16 terms. However, for Exponential-based FET, the relative difference achieves the tolerable range at about 5 terms and then fluctuates around the target accuracy. Figure 2b shows the reconstructed radial distribution from FETs using different basis at different orders against tracklength histogram tally. A 6-term (5th order) of exponential basis in the form of ψ n = exp (nr 6 ) can well represent the discrete results from tracklength estimator. Each point is within the standard deviation of discrete results. As with Zernikebased FETs, expansion with 6 terms is far from sufficient. There is an obvious oscillation as shown in the plot. Comparing to more than 16 terms to be stored in conventional Zernikebased FETs, only a third of that is needed for the new exponential basis. For conservative expansion of using 12 terms, there are still fewer variables to be stored.
Another quick check when generating Figure 2b shows the conservation of total reaction rates on the disk. Applying numerical integration on the functional expansion functions gives us the same total reaction rates no matter what basis to choose and what order to use.

Conclusions
Functional expansion tallies provide the possibility of simplifying the geometrical representation and reducing tally storage needs. However, for reaction rates like the absorption rate of U-238, excessive expansion orders are needed for Zernike-based FETs to achieve a reasonable accuracy. In this work, a generalized functional expansion is introduced to mitigate this problem by using an alternative basis set. There is no requirements on orthogonality of the basis set, and by keeping the zeroth order as unity, the integral reaction rate is preserved. By converting the problem from choosing a better basis to a change of variables, a set of exponential functions in the form of ψ n = exp (nr 6 ) is proposed as a better basis for tallying U-238 absorption reaction rates. The results show that this basis set can well represent the same distribution as traditional discrete tallies. A comparison is made between the new exponential-based FET and Zernike-based FET in a 2D PWR pin cell problem. The variables needed to be stored in this new exponential basis set is approximately half of the orthogonal Zernike-based FET, which can help in data storage and transfer. In addition, fewer expansion could potentially reduce the computational cost in continuous depletion problems [1].
This work provides a heuristic discussion on generalized functional expansion theory. By using fabricated basis sets, steep gradients that occur in neutronics simulation can be well represented by a limited order of expansion. While orthogonal Zernike-based FETs work well for most reactions, the fabricated basis sets are essential in depletion analysis where steep gradients often occur due to spatial self-shielding effects. Further investigation will focus on applying the generalized functional expansion on continuous depletion problems. Funding: This research was supported by the Exascale Computing Project (17-SC-20-SC), a joint project of the U.S. Department of Energy's Office of Science and National Nuclear Security Administration, responsible for delivering a capable exascale ecosystem, including software, applications, and hardware technology, to support the nation's exascale computing imperative.