COMPARISON OF CHEBYSHEV AND ANDERSON ACCELERATIONS FOR THE NEUTRON TRANSPORT EQUATION

This work focuses on the k-eigenvalue problem of the neutron transport equation. The variables of interest are the largest eigenvalue (keff) and the corresponding eigenmode is called the fundamental mode. Mathematically, this problem is usually solved using the power iteration method. However, the convergence of this algorithm can be very slow, especially if the dominance ratio is high as is the case in some reactor physics applications. Thus, the power iteration method has to be accelerated in some ways to improve its convergence. One such acceleration is the Chebyshev acceleration method which has been widely applied to legacy codes. In recent years, nonlinear methods have been applied to solve the k-eigenvalue problem. Nevertheless, they are often compared to the unaccelerated power iteration. Hence, the goal of this paper is to apply the Anderson acceleration to the power iteration, and compare its performance to the Chebyshev acceleration.


INTRODUCTION
This work investigates the solution of the multigroup neutron transport equation with a discrete ordinates (S n ) method. More specifically, it focuses on the k-eigenvalue problem of the equation. The multigroup discrete-ordinate transport equation [1] is written in operator form as: ω n Σ g →g s (r, Ω n · Ω n ) ψ g n (r) ω n νΣ g f ψ g n (r) Therefore, the critical reactor core is obtained by establishing an equilibrium between fission sources and the removal term, without external sources, and is equivalent to solving for the largest eigenvalue k eff such that (k eff , ψ) is the solution to H −1 Fψ = k eff ψ.
In this case, the variables of interest are the largest eigenvalue (k eff ) and the corresponding eigenmode is called the fundamental mode. Mathematically, this problem is usually solved using the power iteration method [2]. However, the convergence of this algorithm can be very slow and, the power iteration method has to be accelerated in some ways to improve its convergence. One such acceleration is the Chebyshev acceleration method [3] which has been applied to legacy codes. In recent years, nonlinear methods have been applied to solve the k-eigenvalue problem. Nevertheless, these are often compared to the unaccelerated power iteration method. Hence, the goal of this paper is to apply the Anderson acceleration method to the power iteration method, and to compare its performance to the Chebyshev acceleration method.

Power iteration
The power iteration method consists in iterating on the fission sources to solve the eigenvalue problem and is as described in Algorithm 1. The approximate eigenvector ψ (n) can be decomposed on where λ i is the i th (in order of decreasing magnitude) eigenvalue, ψ i the associated eigenvector, α i a constant. For thermal reactors, the dominance ratio, κ = λ 2 λ 1 , is usually above 0.95, and the power iteration method converges slowly. Thus, over the years, several acceleration methods have been conceived among which are rebalance techniques, Chebyshev acceleration and more recent nonlinear methods [4].
Apply source iteration for new flux Hψ (n+1) = 1 Compute the new eigenvalue k

Chebyshev acceleration
This consists in modifying the initial problem into an equivalent one with a smaller dominance ratio to improve the convergence of the power iteration algorithm. The Chebyshev acceleration consists in applying power iteration method on a polynomial of H −1 F such that the dominance ratio is as small as possible without changing the eigenvectors of the latter.
Thus, the starting point is to expand ψ (n) on the basis of eigenvectors of A = and by minimising the fractional part of the former expression, the polynomial P l obtained may be applied to accelerate the power iteration [3]. One such polynomial is obtained from the Cheby- where C n (x) = cos(n arccos(x)) if |x| < 1 or C n (x) = cosh(narccosh(x)) otherwise. Following the recurrence property of the Chebyshev polynomials, the acceleration step is hence expressed as: The main drawback of this method is that it requires a priori knowledge of λ 1 and the dominance ratio κ -its efficiency depends on an appropriate estimate for these. Algorithm 2 illustrates an implementation of the Chebyshev acceleration in APOLLO3 R /MINARET [5].

Anderson acceleration
In this research effort, the Anderson acceleration [6] is applied to the power iteration algorithm. The power iteration works with two iterates ψ (n+1) and ψ (n) -it corresponds to a Krylov subspace of dimension 1. The method can be improved by using information from previous iterates. Chebyshev acceleration or Successive Over-Relaxation update the current iteration using a linear combination of previous iterates. In this sense, Anderson acceleration (Anderson mixing in quantum chemistry [7]) attempts to benefit from past iterations. It has been applied in [8] to accelerate neutron-thermal hydraulics problems. Anderson acceleration is an extrapolation method which has been successfully applied recently to neutron transport by [4,9].
Anderson belongs to the class of Newton-Krylov methods except that the subspace projection is orthogonalised using only the number of vectors prescribed by the history length of the method. Unlike Newton methods, it does not require the forming of the Jacobian matrix (similar to Jacobian-Free Newton-Krylov). No formal proof of convergence exists in literature but [7] showed that under some assumptions, convergence is observed in practice. The power iterations described in Algorithm 1 (lines 2 and 3) can be recast as fixed point iterations under the following form: associated with the non-linear problem u = f (u) where u is the ratio of the fission source to the dominant eigenvalue F ψ k ef f . Anderson acceleration can be applied to accelerate such sequences as described in Algorithm 3. Each Anderson iteration requires the solution to a constrained least- Determine α (n) = (α squares problem that is solved in our implementation in a naive way using the normal equations. On the difference of the Chebyshev algorithm, Anderson acceleration offers only two free parameters: the damping coefficient β and the history length m n . Besides, on the computational side, it should be noted that increasing the order of the method leads to higher memory considerations to store the "longer" history of vectors.

NUMERICAL RESULTS
The Chebyshev and Anderson acceleration methods were implemented in a mock-up in-house S n code and several benchmarks have been considered. First, 1D cases with reflector/active core/reflector were studied to obtain a general trend of the acceleration methods. Then, 2D Cartesian benchmarks have been set up with varying core sizes up to sizes with 10 × 10 assemblies of sidelength 20 cm (with larger cores, the power iteration method converges more and more slowly) and one reflector assembly of similar dimensions. All calculations were carried out with S 2 levelsymmetric quadratures [1]. Two-group cross sections with upscattering have been employed as given in Table 1 In forthcoming sections, the results illustrated are for the 2D Cartesian geometry Table 1: Benchmark cross sections. Fuel and reflector share the same cross sections, except the reflector has a nil νΣ f . with 10 × 10 assemblies as it is the most challenging case considered. The reference solution is a well-converged calculation, with a tolerance of 10 −13 with free iterations to avoid any effect of the acceleration method. The plots given next correspond to the error on the fission source compared to the reference solution given against the total number of inner iterations (transport sweeps). These figures are for a maximum of 5 spatial sweep iteration (Richardson), and 5 multigroup iterations (Gauss-Seidel type). The maximum number of power iterations was arbitrarily set to 500. In each plot, the total number of transport sweeps for unaccelerated (or free) iterations are given as N f ree total (if N f ree total = 5 × 5 × 500 = 12500, the unaccelerated calculation is deemed unconverged.)

Anderson acceleration results
Here, some general trends for the Anderson acceleration are highlighted.

Impact of the convergence of the inner iterations
The transport equation is solved using a series of nested iterations (Gauss-Seidel, Richardson). As shown in Fig. 1 (left figure), if unsufficient transport sweeps are carried out, the calculation does not converge at all, as the initial solution vectors are inaccurate. Allowing for more transport sweeps ( Fig. 1 right figure) mitigates this problem.
Besides, other works like [10] concluded that for Newton-like methods, the accuracy of the inner iterations must be increased as the outer iteration approaches the solution. For nested iterations with outer Krylov such as Anderson acceleration, the initial iterates must be computed with good accuracy, which can be relaxed later on in the process [11]. Even when affected, the Chebyshev acceleration is less sensitive to the inner loop convergence compared to the Anderson acceleration.

Impact of the history length
Starting the acceleration method as soon as possible (when sufficient iterations have been carried out to extrapolate the source vector using past iterations for a given order) is not always a right strategy. From Fig. 2, it can be observed that imposing sufficient free iterations beforehand leads to better eigenvectors to start building the solution space, and thus increases the rate of convergence. Besides, it has been observed that for less strict stopping criteria, the impact was more significant as the accuracy is greatly improved. In addition, as the order increases, the effect of the history length is less appreciable.   3 illustrates the effect of the damping coefficient β n on the rate of convergence. Two sets of results are provided: β n of 0.5 and 1, and for minimal history lengths and with a fixed number of initial free iterations. In any case, a damping coefficient of 1 is optimal with respect to the rate of convergence. Furthermore, it appears that as the order increases the impact of the damping coefficient is less significant (∼1000 sweeps difference for order 5, whereas > 5000 for order 1).

Anderson vs. Chebyshev
In any case as shown in Fig. 4, even for less stringent stopping criteria, Anderson 3 -which employs the same number of previous iterations for extrapolation as the Chebyshev acceleration -is more efficient than Chebyshev acceleration. It converges with almost three times less iterations than required by Chebyshev acceleration. With strict tolerances, the Chebyshev algorithm fails to accelerate the source to the required criteria.

CONCLUSIONS
The Anderson acceleration method has been applied successfully to the power iteration as a good candidate to substitute the Chebyshev acceleration. The Anderson acceleration relies on a strong mathematical background with two parameters left to user: the history length and the damping factor. Triggering the acceleration too early is not desirable, sufficient free iterations are required to obtain better initial histories and thus more precise initial vectors for the subspace. As for the damping parameter, it is observed that although it can be varied, a value of unity (1) is optimal for the application of Anderson acceleration to the power iteration method for the neutron transport equation. Anderson acceleration is also sensitive to the convergence of inner iterations and further work must be carried out to investigate the optimisation of the computational cost by introducing relaxation techniques for the convergence of nested loops.