MATRIX EXPONENTIAL METHODS FOR PARALLEL COMPUTING OF ISOTOPIC DEPLETION AND SPECIES TRANSPORT FOR MOLTEN SALT REACTOR ANALYSIS

Matrix exponential methods have long been utilized for isotopic depletion in nuclear fuel calculations. In this paper we discuss the development of such methods in addition to species transport for liquid fueled molten salt reactors (MSRs). Conventional nuclear reactors work with fixed fuel assemblies in which fission products and fissile material do not transport throughout the core. Liquid fueled molten salt reactors work in a much different way, allowing for material to transport throughout the primary reactor loop. Because of this, fission product transport must be taken into account. The set of partial differential equations that apply are discretized into systems of first order ordinary differential equations (ODEs). The exact solution to the set of ODEs is herein being estimated using the matrix exponential method known as the Chebychev Rational Approximation Method (CRAM).


INTRODUCTION
Matrix exponential methods have long been utilized for isotopic depletion in nuclear fuel calculations. These calculations, in traditional light water reactors, involve solving a set of linear ordinary differential equations (ODEs) for stationary fuel assemblies. In liquid fueled molten salt reactors, the fuel and fission products are allowed to flow around the primary loop, thus requiring that species transport be taken into account. Traditional time discretization methods used for scalar data transport cannot achieve the level of accuracy required for nuclear analysis of fission products. If these methods are utilized then very small time steps must be taken to reach high levels of accuracy, thus leading to increased computational time.

METHODS
Nuclear burnup calculations involve solving a set of first order linear ODEs, written as Equation 1 [2]. n = An + S, n(0) = n 0 Where, n(t) is the nuclide concentration vector, A is the transition matrix, S is a vector of constant source terms and n 0 is the initial condition vector. The transition matrix (A) contains the decay and transmutation coefficients.
In situations when S = 0 Equation 1 has the solution n(t) = e At n 0 where e At is defined by the power series shown in Equation 2[2] [3].
When S is not equal to zero, the transition matrix A can be modified toÂ by adding a dummy species to hold constant source terms [4].
For molten salt reactors, these depletion calculations involve solving scalar linear partial differential equations (PDEs) using matrix exponential methods for high time domain accuracy. This method involves first discretizing the spatial dependence of a PDE to a system of coupled ODEs. The resulting set of ODEs will have the same form of Equation 1 with the addition of new coefficients in A representing the approximation of spatial dependence.

Species Transport
The implementation of species transport is accomplished on a 2D finite volume discretization. Equation 3 is the integral formulation of the conservation of species for a single phase [5].
The left hand side of Equation 3 represents the change in mass concentration of species i with time. On the right hand side the first term is the mass flux across a surface followed by the change in concentration due to volumetric generation rate. In finite volume methods, the integral is taken over a volume element of the mesh and will reduce Equation 3 to volume averaged variables. This results in Equation 4 Species flux is further broken down into two contributions, diffusive and convective. The diffusive term is approximated using a second order central difference scheme. The convective term is calculated via a first order upwind difference scheme, the same method used in [5].

Solutions Based on the Cauchy Integral Formula
Computing the matrix exponential is not easy and throughout time there have been many ways of calculating e At [3]. There are various advantages and disadvantages with the approaches shown in the literature, however, this paper examines methods based on Cauchy's Integral Formula. In general, this problem involves solving contour integrals of the form: where Γ denotes the Hankel contour [6]. When employing this method to transition matrix At it is important to choose contour integrals what enclose the spectrum of the transition matrix. If one does not choose an integral that fully encloses the spectrum, then calculating the resulting contour integral will not be a good approximation to the matrix exponential. It can be shown [2], [7], [1] that the approximation to Equation 5 is equivalent to finding a rational function r(z) that minimizes the error where f (z) = (zI −At) −1 , the singularities of f are the eigenvalues of At. For a rational function with simple poles the approximation to e At n 0 takes the form e At n 0 ≈ r k,k (At)n 0 = α 0 n 0 + 2Re where α 0 is the limit at infinity and α j are the residues at poles θ j [1].
In the Chebychev Rational Approximation Method (CRAM), the rational function r k,k (Z) is chosen to be the best rational approximation of Z on the negative real axis [1]. This is done by finding a unique rational functionr k, where π denotes the set of rational function r k,k (x) = p k (x)/q k (x).
Another method for solving Equation 5 is by using quadratures as rational approximations [6]. Equation 5 can be written for an analytic function φ(θ) that maps the real line onto the contour Γ, shown in Equation 9 .
It can be shown [6] that an Nth order approximation to Equation 9 can be calculated via the trapezoidal rule by where r(z) is of the form, From this approximation, the approximate solution n(t) can be achieved by solving Equation 7 where α 0 equal to zero. The coefficients for α j and θ j are found based on the specific φ contour that is chosen. In this article, one type of contour is examined based on the work of Trefenthen, Weideman, and Schmelzer [6]. The parabolic (Equation 12) contour was derived with optimal parameters by balancing error terms from the approximation of contour integrals by a Nth order quadrature [6].

Parallel Computing
Both CRAM and quadrature contour methods involve solving k/2 sets of independent linear systems. Because of the way in which these algorithms work, the ability exist to easily compute these linear systems independently on multiple processors. The best way to divide up the work is to evenly divide the linear systems and send the set to each processor. In this work, MPI is utilized to demonstrate both easy and efficient parallelization of these algorithms.

Parallel Performance
A test matrix is constructed of various sizes to access both the parallel performance and computational efficiency of CRAM. The number of processors is increased by powers of two, up to eight, so that the number of linear systems in evenly divided for each processor. Each matrix size is evaluated 20 times and the average run time is recorded. Figures 1 and 2 show the parallel speedup and efficiency for CRAM. Figure 1 shows good linear speedup for CRAM, especially when the matrix size is large. In Figure 2, there is a point at four processors where a slightly superlinear efficiency is observed. This may point to sub-optimal serial performance for systems of large equations.

Neutron Precursor Analysis
The results presented in this section revolve around a single practice problem which involves the modeling of six delayed neutron precursor groups represented by Equation 13 where ρ k is the mass concentration of each group k. Starting with the first term on the right hand side of Equation 13, the first term is the transport term, the second is generation from fission and the last is decay. These set of six equations are discretized to form a coupled set of first order linear ODEs.
The sample problem in is article models the molten salt reactor experiment (MSRE), however, the results shown are in no way intended to reproduce or mimic the operation of the MSRE. For now, the results shown are only meant to demonstrate the methodological capability, not to model the actual performance of the MSRE. The source terms for fission and decay were obtained using the MPACT neutronics code. The fission source terms were obtained for a single MSRE channel near the middle of the core and follows a cosine shape in the axial direction. Because this problem is a 2-D axial/radial representation, the source terms are scaled with a cosine shape in the radial direction to model a more realistic reactor flux shape. The resulting 2-D axial/radial slice is split up into 14 radial levels and 32 axial levels. Sixteen of the axial levels are located inside the core and 16 are located outside of the core to model the flow loop. The 14 channels start from the center of the core and proceed outward to the core boundary. MSRE contained flow channels but these channels are not modeled, the core is modeled as an open vat. Periodic boundary conditions are placed at the top and bottom of the model however, Dirichlet boundary conditions are also examined to show their impact on the imaginary eigenvalue parts of the system. Flow velocity is assumed constant in the axial direction. In the radial direction, a velocity shape is assumed and taken from work by Andreas van Wijk [8]. Fluid velocity starts low in the center of the core, increases to each a maximum then decays to the reactor boundary.
The main assumption when using either CRAM or the contour denoted by Equation 12 is that the eigenvalues of the transition matrix must be clustered around the negative real axis, or at least in a region where the rational approximation is still valid. For the test previously described, eigenvalues are plotted for nominal flow (black x), 50% nominal flow (red x) and 150% nominal flow (blue x) in Figure 3 for periodic boundary conditions and in Figure 4 for Dirichlet boundary conditions.  Examining Figures 3 and 4 shows the impact that a change in boundary conditions has on the spectrum of the transition matrix. This, in fact, makes sense knowing that the imaginary parts come from damped oscillations raising from the periodic boundary conditions. The fact that the imaginary components of the eigenvalues are quite large seemed quite concerning at first. Both CRAM and the parabolic contour methods require the eigenvalues to be cluster around the negative real axis. To understand the impact this has on the errors, the error for both approximations on the complex plane is shown in Figures 5 and 6, both figures plot three different things. The first is a contour of the log of the relative error on the complex plane. The seconds is each of the quadrature points denoted by black x's. The last is the worst behaving eigenvalues of Figure 3 where the eigenvalue is the slope of the line multiplied by the time step of the problem. Black dots correspond nominal flow, red to 50% flow and blue to 150% flow. Both of these figures contain a lot of information about the relative accuracy of your solution based on how long of a time step you take. For time steps corresponding to one second the eigenvalues of the transition matrix will fall into the plateau where your error is a minimum. If the time step is between 0.1 and 0.6 seconds the 16 order parabolic contour will contain error due to the spectrum falling in a region of high relative error (1e-5 to 1e-2). For the same time steps with CRAM, the error is lower. Results for the simulation are presented in Figure 7, where CRAM is the solution method. These results were generated by taking a single time step of 140 seconds. Cases were also evaluated by taking smaller time steps and examining how the solution behaves in transient scenarios. It was found that by taking a single step to a predetermined time vs taking many small steps to this same time resulted in slightly different answers. Although most of these differences were at or near machine double precision (10 −16 ) they are concerning and more research needs to be done to explain this phenomenon.

CONCLUSIONS
In this work matrix exponential methods were applied on the basis of contour integrals for the solution of partial differential equations relevant to molten salt reactor analysis. It is shown that the spectrum of the transition matrix has a big influence on the accuracy of these solution methods and it is shown that it is possible to utilize CRAM for accurate time domain solutions.