Strategies of Preconditioner Updates for Sequences of Linear Systems Associated with the Neutron Diffusion

The time-dependent neutron di ﬀ usion equation approximates the neutronic power evolution inside a nuclear reactor core. Applying a Galerkin ﬁ nite element method for the spatial discretization of these equations leads to a sti ﬀ semi-discrete system of ordinary di ﬀ erential equations. For time discretization, an implicit scheme is used, which implies solving a large and sparse linear system of equations for each time step. The GMRES method is used to solve these systems because of its fast convergence when a suitable preconditioner is provided. This work explores several matrix-free strategies based on di ﬀ erent updated preconditioners, which are constructed by low-rank updates of a given initial preconditioner. They are two tuned preconditioners based on the bad and good Broyden ’ s methods, initially developed for nonlinear equations and optimization problems, and spectral preconditioners. The e ﬃ ciency of the resulting preconditioners under study is closely related to the selection of the subspace used to construct the update. Our numerical results show the e ﬀ ectiveness of these methodologies in terms of CPU time and storage for di ﬀ erent nuclear benchmark transients, even if the initial preconditioner is not good enough.


Introduction
The safety and the design of a nuclear power plant require, among others, fast and accurate codes that simulate the behaviour of the neutrons inside of the core of a nuclear reactor. The time-dependent multigroup neutron diffusion equation approximates the time evolution of the neutronic power. This equation is an approximation of the neutron transport equation that assumes that the neutron current is proportional to the gradient of the scalar neutron flux by means of a diffusion coefficient. This equation can be expressed as [1].
where for the special case of two energy groups and without considering up-scattering, the continuous operators are given by The vectors Φ 1 and Φ 2 denote the fast and thermal flux, respectively. The vectors C k , k = 1, ⋯, K are the concentration of the neutron delayed precursors, where K is the total number of precursor groups. The cross-sections and diffusion coefficients, Σ a g , Σ f g , Σ 12 D g , g = 1, 2, are functions that depend on the materials of the reactor. ν g , g = 1, 2, is the average number of neutrons released per fission while v g , g = 1, 2, is the neutron velocity associated with the energy group g. The spectrum of the prompt and the delayed neutrons are denoted by χ. The fraction of the delayed neutrons is β k such that the total delayed neutron fraction is β = ∑ K k β k . Finally, the neutron precursor delayed constants are denoted as λ d k . The rest of the coefficients in this work are considered constants related to the neutron delayed precursors. The problem (1) depends on both the position and time; thus, spatial and time discretization must be selected in order to compute the numerical solution. Applying a Galerkin finite element method for the spatial discretization of (1) leads to a stiff semi-discrete system of ordinary differential equations. For time discretization, an implicit scheme is used, which implies solving a linear system of equations for each time step. In a realistic nuclear reactor, the size of the resulting matrices is huge, and Krylov subspace methods are chosen to solve the systems, but these methods may exhibit slow convergence unless a good preconditioner is provided.
In many works, preconditioners for these unstructured linear systems are based on incomplete factorizations of the coefficient matrix [2]. Incomplete LU factorization works well for this type of problem, but it requires the storage of both the full matrix and the preconditioner, and this can be very expensive from the point of view of computational cost and memory requirement. Since the computation of a good preconditioner can be expensive, while solving a sequence of many linear systems, it is advantageous to recycle preconditioners, that is, update a previous preconditioner and reuse the updated version. In, e.g., [3][4][5][6][7] the authors describe updating strategies for incomplete factorizations or approximate inverse preconditioners in the solution of sequences of linear systems. Unfortunately, these techniques do not avoid the need to store the preconditioner. In this work, we focus on updating strategies that can be applied to any initial preconditioner and that can be constructed and applied in a matrix-free regime.
This work explores several matrix-free strategies based on different updated preconditioners, which are constructed by low-rank corrections of a given initial preconditioner. We analyze the performance of two tuned preconditioners based on the bad and good Broyden's methods, initially developed for nonlinear equations and optimization problems, and spectral preconditioners. We study the influence of the selection of the subspace used to construct the update on the efficiency of the resulting preconditioner. The efficiency is measured by both the decrease in the number of iterations needed to reach convergence and the reduction of the overall solution time. We present numerical results showing the effectiveness of these methodologies in terms of CPU time and storage for different nuclear benchmark transients, even when starting from a poor-quality initial preconditioner. Summarizing, the aims of this contribution are (1) to analyze the different low-rank updates, (2) to apply this tool to the sequence of linear systems arising from the discretization of the neutron diffusion equation, and (3) to perform an experimental study on the choice of S i .
The remaining of this paper is organized as follows. In Section 2, we describe the spatial and time discretizations that gives rise to the sequences of sparse nonsymmetric lin-ear systems. Section 3 describes the different low-rank preconditioner updates we consider while in Section 4, we describe the sequences of low-rank preconditioners and analyze different choices for selecting the basis of the subspace used to construct the low-rank update. In Section 5, we present the numerical experiments and discuss the results. The paper is ended by stating the main conclusions of this study in Section 6.

Spatial and Time Discretizations
A high-order finite element method (FEM) is considered for the spatial discretization of the neutron diffusion equation to get a semi-discrete system of ordinary differential equations. This methodology can be applied to any type of reactor geometry. In this work, a continuous Galerkin method with Lagrange polynomials is used. It yields to the system of ordinary differential equations [2] where b Φ, C, L, F, V, and X are the discrete version of the previous differential operators. The implementation of the FEM has been made by using the structures of the open source library Deal.ii [8]. For a realistic nuclear reactor problem, the system (3) is, in general, stiff. Hence, it is necessary to use implicit methods to avoid artificial stepsize restriction due to stability and not accuracy reasons. In this work, a one-step implicit scheme is used, such that the time domain, ½0, T is divided into several time intervals ½t n , t n+1 with Δ t n = t n+1 − t n and the neutron flux at time t n+1 can be approximated by solving the linear system [2]. where Superindex n denotes the matrices and vectors evaluated at time t n .
The neutron precursor equations are approximated by The previous scheme is numerically stable, but due to its implicit character, it requires solving a linear system at each 2 Computational and Mathematical Methods time step. In a realistic nuclear reactor, the size of the matrices is huge, and Krylov subspace methods are chosen to solve the systems, but, to accelerate convergence, a good preconditioner is mandatory. In many works, preconditioners for these unstructured linear systems are based on incomplete factorizations of the coefficient matrix [2]. Incomplete LU factorization works well for this type of problem, but it requires the storage of the full matrix together with the preconditioner, and this can be very expensive from the point of view of memory requirement.
To alleviate the cost of computing the preconditioners, other alternatives, based on low-rank updates of a given incomplete factorization, have been studied, see, e.g., [3][4][5][6], which can be implemented using only vector-matrix products and without requiring additional storage.
In our application, although the matrix operators change at each time step, they are not very different, and the coefficient matrices related to successive time steps share similar spectral properties. Therefore, it is reasonable to explore methodologies to improve the convergence of the Krylov methods by using some information generated when the solutions of the previous linear systems are computed. A technique which takes into account information of the Krylov process to accelerate the iterative solution of subsequent linear systems has been also investigated, e.g., in [9]. In [10], spectral information on the system matrix (which is kept constant throughout the sequence) is extracted and refined from the conjugate gradient solution of previous linear systems and used to deflate the subsequent systems in the sequence.
In this work, an initial preconditioner is updated from this spectral information to improve the convergence of the linear solver along a given transient. This technique has been studied in previous works [3,[11][12][13]. Moreover, these preconditioners have a matrix-free implementation if the initial preconditioner can also be applied in a matrixfree regime, as for instance the polynomial preconditioner recently resumed in [14], reducing the CPU resources of the matrix allocation considerably. Another strategy to update preconditioners for sequences of linear systems, based on the sparse minimization of the Frobenius norm of a suitable error function, has been recently studied in [15].
We describe in the next section different preconditioner updates for a linear system. Next, the methodology is extended to a sequence of linear systems.

Updated Preconditioners for a Linear System
Let us consider a linear system where A ∈ ℝ N×N is a nonsymmetric sparse matrix; x, b ∈ ℝ N is the solution vector and the independent term, respectively. Under suitable conditions, the rate of convergence of Krylov methods such as the GMRES [16] can depend on the distribution of the eigenvalues of A [17]. Even if, in the nonsymmetric instances, this is not directly related to extremal eigenvalues, however, eigenvalues with small modulus are usually responsible of slow convergence of these methods. In this context, we define a left preconditioned system as The solution of (8) is the same as that of system (7) if P is nonsingular; thus, P is designed such that the condition number of matrix PA is smaller than the condition number of A.
Theoretically, the optimal preconditioner is P = A −1 because the matrix of the left preconditioned system is A −1 A = I, which has optimal condition number 1, requiring a single iteration for convergence. However, applying the optimal preconditioner is as difficult as solving the original system, and preconditioners that approximate A −1 are used. Moreover, we are also interested in obtaining a preconditioner that is cheap to construct and apply. Alternatively, the system (19) can be also preconditioned on the right as If the systems are associated with a symmetric and positive definite (SPD) matrix, we can apply the conjugate gradient method, and the convergence does not depend on using right or left preconditioning. However, this behaviour does not generally happen when the GMRES is applied for nonsymmetric matrices. The results can be different if the preconditioner is applied on the left or on the right. Typically, the right preconditioner is preferred since, in this case, the exit test is based on the true residual.

Tuned Preconditioner.
The concept of tuned preconditioner was introduced in [18], and other papers of the same authors, in the framework of iterative eigensolvers. The link between these preconditioner and the quasi-Newton updates, presented in [19], is described in the survey [20]. In this work, we use a block definition of these updates as discussed in [21,22]. Definition 1. Let P 0 be an initial preconditioner of the matrix A and S ∈ ℝ N×k a matrix with full column rank, a tuned preconditioner is a matrix P obtained by updating P 0 by a low-rank matrix depending on A and S so as to satisfy the relation From the previous definition, the matrix PA has the eigenvalue 1 with (at least) multiplicity k associated with the eigenvectors corresponding to each column of S.
Different types of tuned preconditioners can be developed. The following tuned preconditioners are derived from the quasi-Newton methods developed in numerical optimization. These methods provide an approximation of the inverse of the Jacobians, which is equal to A −1 if we consider the problem FðxÞ = 0, for the particular case of FðxÞ = Ax − b [19,23]. Therefore, they are good candidates to design preconditioners for Krylov methods. In particular, two 3 Computational and Mathematical Methods Broyden algorithms are used for this purpose, the Bad-Broyden method and the Good-Broyden method.
3.1.1. The Bad-Broyden Preconditioner. This strategy has been developed in [24], where the preconditioner is designed from rank-1 updates from either the "bad" Broyden's method or Eirola and Nevalinna's method. The Sherman-Morrison formula constructs the sequence of approximations of the inverse of the Jacobian P k+1 ≈ A −1 by a low-rank update of the previous approximate matrix as where s k = Δx k = x k+1 − x k is the (quasi) Newton step and y k = Fðx k+1 Þ − Fðx k Þ. This is known as the Bad-Broyden variant. To define the preconditioner, y k is taken as y k = As k to avoid the restriction of choosing the vectors s k as the Newton step. From these assumptions, the sequence would be This formula can be generalized to a block rankk extension if the column vectors of S = ½s 1 , ⋯s k are in the same subspace that the ones that satisfy the A T A-conjugacy, i.e., [25]. The authors in [24] called this preconditioner limited memory preconditioner by setting P 0 = I N . In our case, we do not take this assumption and P 0 will be an initial preconditioner.
Definition 2. Let A ∈ ℝ N×N be a general nonsingular matrix and S ∈ ℝ N×k of full rank k, with k ≤ N. The Bad-Broyden preconditioner P ∈ ℝ N×N for nonsymmetric systems is defined as Section 4.1 will study the selection of S. This preconditioner satisfies the tuning property (10). The authors [25] show spectral properties of the operator AP (where P 0 = I N ).

The Good-Broyden Preconditioner.
Additionally, one can update the preconditioner by using the Good-Broyden method [20]. This algorithm updates the approximate inverse of the Jacobian P k+1 as where, as before, s k = x k+1 − x k and y k = Fðx k+1 Þ − Fðx k Þ. As in the previous case, choosing y k = As k yields and the previous expression can be extended to a block rankk update if the column vectors of S = ½s 1 , ⋯, s k are in the same subspace as the ones that satisfy the P k A-conjugacy, that is, s T i P k As j = 0, i ≠ j, where i, j = 1 ⋯ , k: Definition 3. Let A ∈ ℝ N×N be a general nonsingular matrix and S ∈ ℝ N×k of full rank k, with k ≤ N. The Good-Broyden preconditioner P ∈ ℝ N×N for nonsymmetric systems is defined as This preconditioner also satisfies the tuning property (10). A selection of S will be studied in Section 4.1.

Spectral Preconditioners.
The spectral preconditioners are based on low-rank corrections obtained from the spectral data associated with the smallest eigenvalues of the preconditioned matrix [26,27].
We start with a preconditioner such that the preconditioned matrix P 0 A is diagonalizable, that is, where Λ = diag ðλ i Þ, where jλ 1 j ≤ ⋯≤jλ n j are the eigenvalues and V = ðv i Þ the associated right eigenvectors. Let V ε be the set of right eigenvectors associated with the set of eigenvaluesjλ i j ≤ ε.
Definition 4. Let S be such that A c = S T AV ε has full rank, a left spectral preconditioner is defined as the matrix As proved in [26], it easy to see that the matrix PA maintains the same right eigenvectors V ε as the matrix P 0 A with eigenvalues η i = λ i + 1, showing that the spectral preconditioner shifts the smallest eigenvalues of the preconditioned matrix away from zero. This is expected to improve the convergence of the Krylov methods. To select the matrix S, different options are proposed in Section 4.1. In [26], a suitable choice of S is proved to leave unchanged all the other eigenvalues and eigenvectors.

Sequence of Preconditioners for Sequences of Linear Systems
Now, let us consider a sequence of linear systems where A i ∈ ℝ N×N are nonsymmetric sparse matrices; x i ∈ ℝ N are the solution vectors, which in this application are the neutron flux in the different time steps; and b i ∈ ℝ N are the independent terms. We are interested in computing a sequence of preconditioners P i for the sequence of linear systems (19). To construct these preconditioners, we use the updated preconditioners of 4 Computational and Mathematical Methods Section 3, where the initial preconditioner is kept constant and matrix S is updated along the sequence. Therefore, the Bad-Broyden preconditioner for the i-th linear system is The Good-Broyden preconditioner is Finally, the spectral preconditioner is In selecting P 0 , two different possibilities are studied in function of the initial preconditioner chosen. If we use the ILU factorization for the initial preconditioner, the ILU of the matrix A 0 is used as P 0 for all systems of the transient. However, if we choose the inverse of the diagonal of the matrix as the initial preconditioner, whose computation is negligible, P 0 corresponds to the inverse of the diagonal matrix A i .

4.1.
Choice of S i . The selection of S i will influence the efficiency of the updated preconditioner in terms of both the number of iterations and the cost of applying the preconditioner. The best option is to set the columns of matrix S i as the eigenvectors of the matrix A i P 0 or P 0 A i (depending on whether right or left preconditioning is used) associated with the eigenvalues of the smallest modulus. However, in practical computations, this is not feasible, and other approximations are commonly used. A practical solution is to approximate these eigenvectors from the Ritz vectors that are defined as follows: Definition 5. Let A be a nonsingular matrix and Q = ½q 1 , ⋯ , q k be a matrix whose columns form a set of orthonormalized vectors. A Ritz vector of A associated with the subspace span fq 1 , ⋯, q k g is defined as the vector w = Qy, where y is an eigenvector of Q T AQ. The eigenvalue θ, corresponding to the eigenvector y, is called Ritz value, and the pair ðθ, wÞ is called the Ritz pair.
It seems reasonable to update the initial preconditioners using the Ritz vectors associated with matrix A i P 0 or P 0 A i . However, numerical results will show that the Ritz vectors associated with matrix A i can be used instead, without any significant detriment of the convergence rate of the iterative method.
The Ritz vectors can be complex because the matrices A i , P 0 A i , or A i P 0 are nonsymmetric. In the case that a Ritz vector (w) with a complex Ritz value (θ) appears, the pair (θ, w) is also a Ritz pair (A i , P 0 A i , and A i P 0 are real). Note that the real and imaginary parts of the Ritz vector w generate the same subspace as the two conjugate Ritz vectors. Therefore, (a) Radial profile    5 Computational and Mathematical Methods the columns of S i are set equal to the real and imaginary parts of the Ritz vector w. Finally, in case any linearly dependent (or nearly) column appears, the last vector is removed.
In the following, we describe two ways to get a set of vectors Q.
If a sequence of linear systems is solved with the GMRES(m) in k iterations, it is assumed that the relation Q T l A i Q l = H l holds with Q l ∈ ℝ N×k and H l ∈ ℝ k×k (k ≤ m). In this case, the Ritz vectors are computed with the basis Q l constructed during the solution of the previous linear system [28,9,10].
Note that if the Ritz vectors associated with the matrix A i are computed, these are obtained as W = Q l Z, where the columns of the matrix Z are the eigenvectors associated with the smallest eigenvalues of H l .

Previous
Solutions of the Systems. The second option we consider is to use the orthonormal basis Q using solutions of the previous linear systems [29,20]. In that sense, to solve the i-th linear system, the matrixQ is defined asQ The column vectors of Q are orthonormal vectors that span the same subspace asQ computed at each iteration before constructing the Ritz vectors.
This selection is based on the reasonable assumption that the components of a solution x i are in the directions of the leftmost eigenvectors of A i , because if the independent terms are expanded in terms of the eigenvectors of On the other hand, in our application, two consecutive matrices in the sequence of linear systems display small changes since the time-step is usually considered small. In principle, this property suggests, but does not guarantee, that also eigenvectors of the two consecutive matrices are close. Experimentally, we found that the eigenvectors of A i closely approximate those eigenvectors of A i+1 , for a given i.   Computational and Mathematical Methods [30]. It corresponds to an operational transient of the Langenbuch reactor. This reactor is a small LWR composed of 77 fuel assemblies and two types of fuel. Figure 1 shows the geometry of the reactor model, whose spatial discretization is composed of 1170 cells. The spatial discretization of the problem is made with a degree of the polynomials in the finite element method equal to 3 to obtain a problem of size N = 69440. Materials 4 and 6 represent control rods. It has been initiated by the withdrawal of a bank of four partially inserted control rods (C1 in Figure 1) at a rate of 3 cm/s over 0 < t < 26:7 s. A second bank of control rods (C2 in Figure 1) is inserted at the same rate over 7:5 < t < 47:5 s. The transient is followed during 60 s. Figure 1(b) represents the axial profile at t = 0:0 s. The evolution of global power can be observed in Figure 2.

3D Langenbuch Transient. This transient was defined in
As a first case, a sequence of 10 linear systems is used to test the performance of the updated preconditioners. They are obtained from the backward difference method (Equation (4)) with Δt = 1:0 s. The positions corresponding to the different systems are marked in Figure 2 with red dots. GMRES with left preconditioning (LGMRES) is applied to solve the linear systems. The initial vectors in the iterative process are set equal to a vector of zeros for all systems. The tolerances for the stopping criterion have been set to tol = 10 −8 and the exit test is based on the unpreconditioned residual kb i − A i x i k 2 < tol. In this problem, MATLAB is used to compute the numerical results.
First, we present the number of iterations to solve the first linear system A 1 x 1 = b 1 by using the LGMRES and updated preconditioners from 5 eigenvectors of A and P 0 A associated with the smallest eigenvalues ( Table 1). The eigenvalues are computed by using the command eigs of MATLAB. The initial preconditioner P 0 is set equal to the ILU (0) of the matrix A 1 obtained with the MATLAB command ilu. This table shows that all the updated preconditioners with the eigenvectors of P 0 A 1 and A improve the results obtained when the initial P 0 (with no update) is used. Interestingly, one can observe that using the eigenvectors of A 1 provides similar results as those obtained when the eigenvectors of P 0 A 1 are used.

Computational and Mathematical Methods
Now, we test the performance of the updated preconditioners for the sequence of linear systems. In the following, we select the initial preconditioner P 0 as the ILU (0) preconditioner of matrix A 1 for all computations.
First, we test the type of vectors used to construct matrix S i . Rank-5 corrections are taken in all cases in the updated preconditioners. Figure 3 shows the number of iterations of the LGMRES method needed to reach convergence with different choices of the columns of matrix S, namely, the eigenvectors of the matrix A i (eigs), the Ritz vectors associated with the matrix A i from the basis obtained from GMRES (gmres), or the Ritz vectors associated with matrix A i and the basis obtained from the solutions of the previous linear systems (sols). Each subfigure reports the results obtained for each type of updated preconditioner (good or bad Broyden or spectral). The figures also show the number of iterations obtained if we apply only the preconditioner P 0 kept fixed without updating for all the linear systems in the sequence.
Broyden's preconditioners (Figures 3(a) and 3(b)) with the "exact" eigenvectors (eigs) provide a constant reduction of the number of iterations, with respect to using P 0 fixed, from the very first system. In contrast, the selection of the previous solutions (sols) and the vectors of the GMRES method (gmres) improves the convergence after a few linear systems are solved, even though from the 8th linear system, the improvement deteriorates. The 8th linear system corresponds to the transient at t = 8 s. At t = 7:5 s, a bank of control rods is inserted to shut down the reactor. Therefore, the subspace obtained with the 1-7th linear systems is not as useful to construct the preconditioner for the 8th linear system as we have observed for the systems 1-7. In the spectral preconditioner (Figure 3(c)), the reduction for all vectors of the number of iterations is more constant. For all updated preconditioners, the selection of the previous solutions (sols) speeds up the convergence of the LGMRES more than when the Ritz vectors constructed from GMRES (gmres) are used.    Table 2: CPU times (s) of LGMRES to solve the sequence of 10 linear systems with the updated preconditioners and P 0 equal ILU (0). "P i construction" collects the CPU time to construct the ILU (0) preconditioner of the first matrix and the CPU time to compute the subspace S i . "P i application" contains the total CPU time in the computation of the 10 linear systems devoted to apply the P i . "Total" is the total CPU time spent to obtain the solution of the 10 linear systems. We investigate in Figure 4 the effect of using different maximum ranks (2,5,10) to construct the matrix S i to correct the initial preconditioner. In this experiment, the columns of matrix S i have been set using the previous solutions. For Broyden's preconditioners, the results show that corrections of maximum rank-5 improve the convergence. Corrections of maximum rank-10 only slighlty reduce the number of iterations because in two cases, the new solutions are linearly dependent from the previous ones, and the rank of the subspaces for these linear systems is quite similar (see Subsection 5.1.1 for an analysis of the effect of the maximum rank on the number of iterations to compute 20 linear systems). Table 2 shows the total number of iterations and CPU times for solving the sequence of systems with different preconditioners and maximum ranks. We also report the CPU times devoted to construct and apply the preconditioners. We observe that the updated preconditioners are more convenient than using the fixed ILU(A 1 ) preconditioner, both in terms of number of iterations and CPU time, due also to the fact that the cost per iteration is only slightly increased by the preconditioner update. The best performance is obtained with the Broyden updates. In particular, the Good-Broyden updates are slightly more efficient when considering the number of iterations. This improvement will be more evident for more difficult cases, as documented in Section 5.2. Figure 5 compares the number of iterations of LGMRES to solve the sequence of linear systems with the updated preconditioners constructed using rank-5 corrections and P 0 equal ILU (0). The Ritz vectors associated with the previous solutions are chosen as set S i . Broyden's preconditioners (with similar results for both strategies) provide a smaller number of iterations than the spectral preconditioner.
We note that the studied updated preconditioners require only matrix-vector products if the application of the initial preconditioner P 0 relies on this operation. If P 0 is taken as the ILU (0) factorization, then the complete assembly of the matrix A 1 is required. P 0 can be set as the inverse of the diagonal of the matrices A i to avoid assembling this matrix. The performance of the updated preconditioners with LGMRES when P 0 is equal to the inverse of the diagonal is displayed in Figure 6. This figure shows that the updated preconditioners decrease the number of iterations.
Moreover, note that, for this application, Broyden's preconditioners are better updating strategies than the spectral preconditioner, independently of the type of P 0 used.

Computational and Mathematical Methods
However, one can observe that if the inverse of the diagonal is used, the number of iterations to reach the convergence is much higher than when the ILU (0) factorization is used. Selecting P 0 as a diagonal preconditioner can be suitable for huge problems with high memory demands.

Case 2:
A Sequence of 20 Linear Systems. To conclude this section, we study the performance of the updated preconditioners for solving a sequence of 20 linear systems associated with the previous transient. Figure 7 reports the number of iterations obtained with the updated preconditioners by using different maximum ranks (2, 5, 10) to construct the matrix S i . The columns of matrix S i are set using the previous solutions. The initial preconditioner used is the ILU (0) of the matrix A 0 . For this sequence, we can observe the improvement, in terms of reduction of the number of iterations, provided by the rank-10 Broyden preconditioners. Regarding the spectral preconditioner, increasing the rank does not translate in further reduction of the number of iterations.

AER-DYN-001 Problem.
To test the efficiency of the updated preconditioner in a more realistic nuclear benchmark, the AER-DYN-001 transient is used. This problem was introduced in [31]. It corresponds to an asymmetric control rod ejection accident without feedback in a VVER440 reactor. The core has 250 cm height, with two reflector layers of 25 cm each added, one to the top and the other one to the bottom of the core. The assembly pitch is 14.7 cm. The core is a VVER-440 core type, with 25 fuel elements across the diameter. The disposition of materials together with the initial position of the control rods is shown in Figure 8. For the spatial discretization, the deal. II library cannot handle hexagonal cells, so each hexagon is subdivided into three quadrilaterals to obtain a total number of 15156 cells. Cubic polynomials are used for the finite element method to obtain algebraic problems of size 857882 degrees of freedom.
The transient is defined as follows. The control rod denoted by number 26 is ejected in the first 0:08 s with a velocity 25 m/s. Then, scram is initiated, inserting the safety rods 23 and 25 at t = 1:0 s with a velocity 0:25 m/s, so that the bottom position is reached at t = 11:0 s. The drop of control rod group 21 is also started at 1:0 s with the same velocity. The transient is followed during 4 s. Figure 9 represents the global power of the transient. For the time discretization, Δt = 0:01 s is used to obtain a squence of 400 linear systems to be solved. The methodology for this reactor has been implemented in C++ by using Deal.II and PETSc structures. The computer for the calculations is an Intel Core i7-4790 @ 3:60 GHz × 8 processor with 32 Gb of RAM running on Ubuntu GNU/Linux 18:04 LTS.
In this problem, we use the Good-Broyden preconditioner, where the columns of S i are the five Ritz vectors associated with the smallest Ritz values of A i and the previous solutions. The left GMRES of the PETSc library is applied to solve the linear systems. For the initial preconditioner, first, we use as P 0 the Block Gauss-Seidel preconditioner of A 1 (BGSðA 1 Þ), unless the number of iterations of a linear system j is greater than 50. In such case, the new P 0 is computed as the Block Gauss-Seidel preconditioner of the matrix A j . This preconditioner permits the use of a semi matrix-free implementation of the matrices (see [32] for more details). Figure 10 shows the number of iterations (Figure 10(a)) and the CPU time (Figure 10(b)) needed to reach convergence to the required accuracy for each linear system. The   total CPU time of the computation is equal to 19.9 hours (no update) whereas it is equal to 18.7 hours if P 0 is updated by the Good-Broyden formula. In this case, the updated preconditioner does not provide a significant reduction neither of the number of iterations nor the CPU time.
However, as in the previous case, we propose an alternative to set P 0 , which can be implemented using only matrixvector products. Chebyshev polynomial preconditioners are recommended for parallel computations and matrix-free implementations [33]. In this work, five iterations of the Chebyshev preconditioner, whose implementation is provided by the Deal.II library [8], are applied. The application of the Chebyshev method requires estimating the largest eigenvalue of A i . For that purpose, in this work, we use ten iterations of the GMRES method preconditioned with the inverse diagonal of A i . Figure 11 shows the number of iterations (on the left) and the CPU time (on the right) that each linear system needs to reach convergence using P 0 either without updating or updated with the Good-Broyden method. The total CPU time is 21.2 hours, when the P 0 is not updated, while it reduces to 9.5 hours, with the Good-Broyden preconditioner. We observe a large improvement in the results if the P 0 (=chebyshev) is corrected. Now, we compare these results (where P 0 = chebyshevð A i Þ) with the previous (where P 0 = BGSðA 1 Þ). Table 3 collects the total CPU times and the number of iterations that GMRES needed with each preconditioner to reach convergence. If the preconditioners are not updated, the Chebyshev preconditioner is not as efficient as Block Gauss-Seidel in terms of CPU time. However, if we use the Good-Broyden to correct the initial setting, the Chebyshev preconditioner, in addition to its negligible memory requirements, reduces the CPU time by more than a half.

Conclusions
This work describes and compares some strategies of updated preconditioners to speed up the computation of the sequence of linear systems obtained from the time discretization of the neutron diffusion equation. In particular, three types of techniques are considered: two tuned-type preconditioners, the Bad-Broyden and the Good-Broyden, and a spectral preconditioner. All of them are based on low-rank corrections of an initial preconditioner. Theoretically, they are constructed by matrix-vector multiplications that involve a set of eigenvectors of the preconditioned matrix. Moreover, based on different subspaces, cheaper options to define these preconditioners are tested.
Numerical results show that the low-rank updates presented in this work accelerate the convergence of the GMRES, with tuned preconditioners a better option than the spectral preconditioner. We found Broyden preconditioners to be less sensitive to the selection of the basis vectors used to construct the update in comparison with the spectral preconditioner. In fact, the vectors forming the basis for the updates for Broyden's preconditioner can be obtained from the solutions to the previous linear systems in the sequence. The cost-free choice proposed in this work, namely, the Ritz vectors extracted from the GMRES iteration, does not always guarantee a sufficient acceleration of the spectral preconditioner. In any case, low-rank corrections are generally enough to obtain efficient preconditioners for our problems.
The CPU times needed to apply the updated preconditioners are higher than those for applying the initial preconditioner only. However, these preconditioners enjoy a matrix-free implementation if the initial preconditioner can be applied in a matrix-free regime such as the Chebyshev preconditioner. In this case, low-rank updated preconditioners allow for great savings in storage and also on the CPU time needed to assemble the matrices of the sequence.
Results on a realistic nuclear benchmark requiring the solution of a sequence of 400 linear systems, each one of size of more than 8.5 × 10 5 degrees of freedom, show that the Good-Broyden update applied to the matrix-free Chebyshev preconditioner provides an impressive reduction of the CPU time taking 9.5 hours to complete the simulation vs. the 21.2 hours needed without updating.

Data Availability
The data used to support the findings of this study are included or referenced within the article.

Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.