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

Authors: A. Carreño; L. Bergamaschi; A. Martínez; D. Ginestar; A. Vidal-Ferràndiz; G. Verdú
Journal: Computational and Mathematical Methods
 (2022)

---

## Abstract

The time-dependent neutron diffusion equation approximates the neutronic power evolution inside a nuclear reactor core. Applying a Galerkin finite element method for the spatial discretization of these equations leads to a stiff semi-discrete system of ordinary differential 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 different 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 efficiency 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 effectiveness of these methodologies in terms of CPU time and storage for different nuclear benchmark transients, even if the initial preconditioner is not good enough.


---

## Body

## 1. 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].
(1)V−1∂Φ∂t+LΦ=1−βFΦ+∑k=1KλkdχCk,∂Ck∂t=βkF1Φ−λkdCk,k=1,⋯,K,where for the special case of two energy groups and without considering up-scattering, the continuous operators are given by
(2)L=−∇⟶·D1∇⟶+Σa1+Σ120−Σ12−∇⟶·D2∇⟶+Σa2,F=ν1Σf1ν2Σf200,V=1/v1001/v2,F1=ν1Σf1ν2Σf2,χ=10T,Φ=Φ1Φ2T.The vectorsΦ1 and Φ2 denote the fast and thermal flux, respectively. The vectors Ck, 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, Σag, Σfg, Σ12Dg, 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 vg,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 β=∑kKβk. Finally, the neutron precursor delayed constants are denoted as λkd. 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–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 ofSi.The remaining of this paper is organized as follows. In Section2, we describe the spatial and time discretizations that gives rise to the sequences of sparse nonsymmetric linear 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.


## 2. 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](3)VdΦ^dt+LΦ¯=1−βFΦ¯+∑p=1NpλpdXCp,dPCpdt=βkF1Φ^−λpdPCp,p=1,⋯,Np,where Φ^, 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 tn,tn+1 with Δtn=tn+1−tn and the neutron flux at time tn+1 can be approximated by solving the linear system [2].
(4)Tn+1Φ∧n+1=1ΔthVnΦ∧n+∑p=1Npλpde−λpdΔthXCpn,where
(5)Tn+1=1ΔtnVn+1+Ln+1−a^Fn+1,a^=1−β+∑p=1Npβp1−e−λpdΔth.Superindexn denotes the matrices and vectors evaluated at time tn.The neutron precursor equations are approximated by(6)PCpn+1=PCpne−λpdΔtn+βpλpd1−e−λpdΔthF1n+1Φn+1,p=1,⋯,Np.The previous scheme is numerically stable, but due to its implicit character, it requires solving a linear system at each 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–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–13]. Moreover, these preconditioners have a matrix-free implementation if the initial preconditioner can also be applied in a matrix-free 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.


## 3. Updated Preconditioners for a Linear System

Let us consider a linear system(7)Ax=b,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
(8)PAx=Pb.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 isP=A−1 because the matrix of the left preconditioned system is A−1A=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
(9)APy=b,x=Py.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.


### 3.1. 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.
LetP0 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 P0 by a low-rank matrix depending on A and S so as to satisfy the relation
(10)PAS=S.
From the previous definition, the matrixPA 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 toA−1 if we consider the problem Fx=0, for the particular case of Fx=Ax−b [19, 23]. Therefore, they are good candidates to design preconditioners for Krylov methods. In particular, two 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 Pk+1≈A−1 by a low-rank update of the previous approximate matrix as
(11)Pk+1=Pk+sk−PkykykTykTyk,where sk=Δxk=xk+1−xk is the (quasi) Newton step and yk=Fxk+1−Fxk. This is known as the Bad-Broyden variant. To define the preconditioner, yk is taken as yk=Ask to avoid the restriction of choosing the vectors sk as the Newton step. From these assumptions, the sequence would be
(12)Pk+1=PkIN−AskskTATskTATAsk+skskTATskTATAsk.This formula can be generalized to a block rank-k extension if the column vectors of S=s1,⋯sk are in the same subspace that the ones that satisfy the ATA-conjugacy, i.e., siTATAsj=0, i≠j, where i,j=1⋯,k [25]. The authors in [24] called this preconditioner limited memory preconditioner by setting P0=IN. In our case, we do not take this assumption and P0 will be an initial preconditioner.Definition 2.
LetA∈ℝ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
(13)P=P0IN−ASSTATAS−1STAT+SSTATAS−1STAT.
Section4.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 P0=IN).


#### 3.1.2. 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 Pk+1 as
(14)Pk+1=Pk+sk−PkykskTPkskTPkyk,where, as before, sk=xk+1−xk and yk=Fxk+1−Fxk. As in the previous case, choosing yk=Ask yields
(15)Pk+1=Pk+sk−PkAskskTPkskTPkAsk,and the previous expression can be extended to a block rank- k update if the column vectors of S=s1,⋯,sk are in the same subspace as the ones that satisfy the PkA-conjugacy, that is, siTPkAsj=0, i≠j, where i,j=1⋯,k.Definition 3.
LetA∈ℝ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
(16)P=P0−P0AS−SSTP0AS−1STP0.
This preconditioner also satisfies the tuning property (10). A selection of S will be studied in Section 4.1.


### 3.2. 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 matrixP0A is diagonalizable, that is,
(17)P0A=VΛV−1,where Λ=diagλi, where λ1≤⋯≤λn are the eigenvalues and V=vi the associated right eigenvectors. Let Vε be the set of right eigenvectors associated with the set of eigenvaluesλi≤ε.Definition 4.
LetS be such that Ac=STAVε has full rank, a left spectral preconditioner is defined as the matrix
(18)P=P0+VεAc−1ST.
As proved in [26], it easy to see that the matrix PA maintains the same right eigenvectors Vε as the matrix P0A 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.


## 3.1. 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.
LetP0 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 P0 by a low-rank matrix depending on A and S so as to satisfy the relation
(10)PAS=S.
From the previous definition, the matrixPA 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 toA−1 if we consider the problem Fx=0, for the particular case of Fx=Ax−b [19, 23]. Therefore, they are good candidates to design preconditioners for Krylov methods. In particular, two 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 Pk+1≈A−1 by a low-rank update of the previous approximate matrix as
(11)Pk+1=Pk+sk−PkykykTykTyk,where sk=Δxk=xk+1−xk is the (quasi) Newton step and yk=Fxk+1−Fxk. This is known as the Bad-Broyden variant. To define the preconditioner, yk is taken as yk=Ask to avoid the restriction of choosing the vectors sk as the Newton step. From these assumptions, the sequence would be
(12)Pk+1=PkIN−AskskTATskTATAsk+skskTATskTATAsk.This formula can be generalized to a block rank-k extension if the column vectors of S=s1,⋯sk are in the same subspace that the ones that satisfy the ATA-conjugacy, i.e., siTATAsj=0, i≠j, where i,j=1⋯,k [25]. The authors in [24] called this preconditioner limited memory preconditioner by setting P0=IN. In our case, we do not take this assumption and P0 will be an initial preconditioner.Definition 2.
LetA∈ℝ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
(13)P=P0IN−ASSTATAS−1STAT+SSTATAS−1STAT.
Section4.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 P0=IN).


### 3.1.2. 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 Pk+1 as
(14)Pk+1=Pk+sk−PkykskTPkskTPkyk,where, as before, sk=xk+1−xk and yk=Fxk+1−Fxk. As in the previous case, choosing yk=Ask yields
(15)Pk+1=Pk+sk−PkAskskTPkskTPkAsk,and the previous expression can be extended to a block rank- k update if the column vectors of S=s1,⋯,sk are in the same subspace as the ones that satisfy the PkA-conjugacy, that is, siTPkAsj=0, i≠j, where i,j=1⋯,k.Definition 3.
LetA∈ℝ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
(16)P=P0−P0AS−SSTP0AS−1STP0.
This preconditioner also satisfies the tuning property (10). A selection of S will be studied in Section 4.1.


## 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 Pk+1≈A−1 by a low-rank update of the previous approximate matrix as
(11)Pk+1=Pk+sk−PkykykTykTyk,where sk=Δxk=xk+1−xk is the (quasi) Newton step and yk=Fxk+1−Fxk. This is known as the Bad-Broyden variant. To define the preconditioner, yk is taken as yk=Ask to avoid the restriction of choosing the vectors sk as the Newton step. From these assumptions, the sequence would be
(12)Pk+1=PkIN−AskskTATskTATAsk+skskTATskTATAsk.This formula can be generalized to a block rank-k extension if the column vectors of S=s1,⋯sk are in the same subspace that the ones that satisfy the ATA-conjugacy, i.e., siTATAsj=0, i≠j, where i,j=1⋯,k [25]. The authors in [24] called this preconditioner limited memory preconditioner by setting P0=IN. In our case, we do not take this assumption and P0 will be an initial preconditioner.Definition 2.
LetA∈ℝ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
(13)P=P0IN−ASSTATAS−1STAT+SSTATAS−1STAT.
Section4.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 P0=IN).


## 3.1.2. 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 Pk+1 as
(14)Pk+1=Pk+sk−PkykskTPkskTPkyk,where, as before, sk=xk+1−xk and yk=Fxk+1−Fxk. As in the previous case, choosing yk=Ask yields
(15)Pk+1=Pk+sk−PkAskskTPkskTPkAsk,and the previous expression can be extended to a block rank- k update if the column vectors of S=s1,⋯,sk are in the same subspace as the ones that satisfy the PkA-conjugacy, that is, siTPkAsj=0, i≠j, where i,j=1⋯,k.Definition 3.
LetA∈ℝ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
(16)P=P0−P0AS−SSTP0AS−1STP0.
This preconditioner also satisfies the tuning property (10). A selection of S will be studied in Section 4.1.


## 3.2. 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 matrixP0A is diagonalizable, that is,
(17)P0A=VΛV−1,where Λ=diagλi, where λ1≤⋯≤λn are the eigenvalues and V=vi the associated right eigenvectors. Let Vε be the set of right eigenvectors associated with the set of eigenvaluesλi≤ε.Definition 4.
LetS be such that Ac=STAVε has full rank, a left spectral preconditioner is defined as the matrix
(18)P=P0+VεAc−1ST.
As proved in [26], it easy to see that the matrix PA maintains the same right eigenvectors Vε as the matrix P0A 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.


## 4. Sequence of Preconditioners for Sequences of Linear Systems

Now, let us consider a sequence of linear systems(19)Aixi=bi,1≤i≤imax,where Ai∈ℝN×N are nonsymmetric sparse matrices; xi∈ℝN are the solution vectors, which in this application are the neutron flux in the different time steps; and bi∈ℝN are the independent terms.We are interested in computing a sequence of preconditionersPi for the sequence of linear systems (19). To construct these preconditioners, we use the updated preconditioners of Section 3, where the initial preconditioner is kept constant and matrix S is updated along the sequence.Therefore, the Bad-Broyden preconditioner for thei-th linear system is
(20)Pi=P0IN−AiSiSiTAiTAiSi−1SiTAiT+SiSiTAiTAiSi−1SiTAiT.The Good-Broyden preconditioner is(21)Pi=P0−P0AiSi−SiSiTP0AiSi−1SiTP0.Finally, the spectral preconditioner is(22)Pi=P0+VεSiTAiVε−1SiT.In selectingP0, 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 A0 is used as P0 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, P0 corresponds to the inverse of the diagonal matrix Ai.


### 4.1. Choice ofSi

The selection ofSi 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 Si as the eigenvectors of the matrix AiP0 or P0Ai (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.
LetA be a nonsingular matrix and Q=q1,⋯,qk be a matrix whose columns form a set of orthonormalized vectors. A Ritz vector of A associated with the subspace span q1,⋯,qk is defined as the vector w=Qy, where y is an eigenvector of QTAQ. 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 matrixAiP0 or P0Ai. However, numerical results will show that the Ritz vectors associated with matrix Ai can be used instead, without any significant detriment of the convergence rate of the iterative method.The Ritz vectors can be complex because the matricesAi, P0Ai, or AiP0 are nonsymmetric. In the case that a Ritz vector (w) with a complex Ritz value (θ) appears, the pair (θ¯,w¯) is also a Ritz pair (Ai, P0Ai, and AiP0 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, the columns of Si 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 vectorsQ.If a sequence of linear systems is solved with the GMRES(m) in k iterations, it is assumed that the relation QlTAiQl=Hl holds with Ql∈ℝN×k and Hl∈ℝk×k (k≤m). In this case, the Ritz vectors are computed with the basis Ql constructed during the solution of the previous linear system [28, 9, 10].Note that if the Ritz vectors associated with the matrixAi are computed, these are obtained as W=QlZ, where the columns of the matrix Z are the eigenvectors associated with the smallest eigenvalues of Hl.


#### 4.1.1. Previous Solutions of the Systems

The second option we consider is to use the orthonormal basisQ using solutions of the previous linear systems [29, 20]. In that sense, to solve the i-th linear system, the matrix Q~ is defined as
(23)Q~=xi−p,⋯,xi−1.The column vectors ofQ are orthonormal vectors that span the same subspace as Q~ computed at each iteration before constructing the Ritz vectors.This selection is based on the reasonable assumption that the components of a solutionxi are in the directions of the leftmost eigenvectors of Ai, because if the independent terms are expanded in terms of the eigenvectors of Ai, bi=∑j=1Nαjuji, then xk=∑j=1Nαj/λjiuji with the largest weights provided by the smallest eigenvalues.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 ofAi closely approximate those eigenvectors of Ai+1, for a given i.


## 4.1. Choice ofSi

The selection ofSi 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 Si as the eigenvectors of the matrix AiP0 or P0Ai (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.
LetA be a nonsingular matrix and Q=q1,⋯,qk be a matrix whose columns form a set of orthonormalized vectors. A Ritz vector of A associated with the subspace span q1,⋯,qk is defined as the vector w=Qy, where y is an eigenvector of QTAQ. 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 matrixAiP0 or P0Ai. However, numerical results will show that the Ritz vectors associated with matrix Ai can be used instead, without any significant detriment of the convergence rate of the iterative method.The Ritz vectors can be complex because the matricesAi, P0Ai, or AiP0 are nonsymmetric. In the case that a Ritz vector (w) with a complex Ritz value (θ) appears, the pair (θ¯,w¯) is also a Ritz pair (Ai, P0Ai, and AiP0 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, the columns of Si 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 vectorsQ.If a sequence of linear systems is solved with the GMRES(m) in k iterations, it is assumed that the relation QlTAiQl=Hl holds with Ql∈ℝN×k and Hl∈ℝk×k (k≤m). In this case, the Ritz vectors are computed with the basis Ql constructed during the solution of the previous linear system [28, 9, 10].Note that if the Ritz vectors associated with the matrixAi are computed, these are obtained as W=QlZ, where the columns of the matrix Z are the eigenvectors associated with the smallest eigenvalues of Hl.


### 4.1.1. Previous Solutions of the Systems

The second option we consider is to use the orthonormal basisQ using solutions of the previous linear systems [29, 20]. In that sense, to solve the i-th linear system, the matrix Q~ is defined as
(23)Q~=xi−p,⋯,xi−1.The column vectors ofQ are orthonormal vectors that span the same subspace as Q~ computed at each iteration before constructing the Ritz vectors.This selection is based on the reasonable assumption that the components of a solutionxi are in the directions of the leftmost eigenvectors of Ai, because if the independent terms are expanded in terms of the eigenvectors of Ai, bi=∑j=1Nαjuji, then xk=∑j=1Nαj/λjiuji with the largest weights provided by the smallest eigenvalues.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 ofAi closely approximate those eigenvectors of Ai+1, for a given i.


## 4.1.1. Previous Solutions of the Systems

The second option we consider is to use the orthonormal basisQ using solutions of the previous linear systems [29, 20]. In that sense, to solve the i-th linear system, the matrix Q~ is defined as
(23)Q~=xi−p,⋯,xi−1.The column vectors ofQ are orthonormal vectors that span the same subspace as Q~ computed at each iteration before constructing the Ritz vectors.This selection is based on the reasonable assumption that the components of a solutionxi are in the directions of the leftmost eigenvectors of Ai, because if the independent terms are expanded in terms of the eigenvectors of Ai, bi=∑j=1Nαjuji, then xk=∑j=1Nαj/λjiuji with the largest weights provided by the smallest eigenvalues.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 ofAi closely approximate those eigenvectors of Ai+1, for a given i.


## 5. Numerical Results

### 5.1. 3D Langenbuch Transient

This transient was defined in [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.Figure 1
Geometry and materials distribution of the Langenbuch reactor.
(a)
Radial profile(b)
Axial profileMaterials 4 and 6 represent control rods. It has been initiated by the withdrawal of a bank of four partially inserted control rods (C1 in Figure1) 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.Figure 2
Global power of the Langenbuch transient.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 bi−Aixi2<tol. In this problem, MATLAB is used to compute the numerical results.First, we present the number of iterations to solve the first linear systemA1x1=b1 by using the LGMRES and updated preconditioners from 5 eigenvectors of A and P0A associated with the smallest eigenvalues (Table 1). The eigenvalues are computed by using the command eigs of MATLAB. The initial preconditioner P0 is set equal to the ILU (0) of the matrix A1 obtained with the MATLAB command ilu. This table shows that all the updated preconditioners with the eigenvectors of P0A1 and A improve the results obtained when the initial P0 (with no update) is used. Interestingly, one can observe that using the eigenvectors of A1 provides similar results as those obtained when the eigenvectors of P0A1 are used.Table 1
Number of iterations of LGMRES to solve the linear problemAx=b with the different updated preconditioners and choices of S. P0=ILUA.
SNo updateBad-BroydenGood-BroydenSpectraleigs(P0A)45313432eigs(A)45323435Now, we test the performance of the updated preconditioners for the sequence of linear systems. In the following, we select the initial preconditionerP0 as the ILU (0) preconditioner of matrix A1 for all computations.First, we test the type of vectors used to construct matrixSi. 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 Ai (eigs), the Ritz vectors associated with the matrix Ai from the basis obtained from GMRES (gmres), or the Ritz vectors associated with matrix Ai 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 P0 kept fixed without updating for all the linear systems in the sequence.Figure 3
Number of iterations of the LGMRES method to solve the sequence of linear systems with the updated preconditioners with rank 5 corrections andP0 equal to ILU (0).
(a)
Bad-Broyden(b)
Good-Broyden(c)
SpectralBroyden’s preconditioners (Figures3(a) and 3(b)) with the “exact” eigenvectors (eigs) provide a constant reduction of the number of iterations, with respect to using P0 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.We investigate in Figure4 the effect of using different maximum ranks (2,5,10) to construct the matrix Si to correct the initial preconditioner. In this experiment, the columns of matrix Si 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).Figure 4
Number of iterations of the LGMRES method to solve the sequence of linear systems with the updated preconditioners and different rank corrections whereP0= is ILU (0).
(a)
Bad-Broyden-sols(b)
Good-Broyden-sols(c)
Spectral-solsTable2 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(A1) 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.Table 2
CPU times (s) of LGMRES to solve the sequence of 10 linear systems with the updated preconditioners andP0 equal ILU (0). “Pi construction” collects the CPU time to construct the ILU (0) preconditioner of the first matrix and the CPU time to compute the subspace Si. “Pi application” contains the total CPU time in the computation of the 10 linear systems devoted to apply the Pi. “Total” is the total CPU time spent to obtain the solution of the 10 linear systems.
PreconditionerRankTotal its.CPU time (s)Pi constructionPi applicationPi application per it.TotalILU(A1)—4561.38.10.01828.1Bad-Broyden-sols23651.47.50.02125.653281.56.80.02123.3103131.56.60.02122.5Good-Broyden-sols23561.47.50.02123.853211.57.10.02123.1103081.57.00.02222.6Spectral-sols24011.48.10.02026.954101.58.10.02027.0104121.58.20.02027.2Figure5 compares the number of iterations of LGMRES to solve the sequence of linear systems with the updated preconditioners constructed using rank-5 corrections and P0 equal ILU (0). The Ritz vectors associated with the previous solutions are chosen as set Si. Broyden’s preconditioners (with similar results for both strategies) provide a smaller number of iterations than the spectral preconditioner.Figure 5
Number of LGMRES iterations to solve the sequence of linear systems with the updated preconditioners, 5-rank corrections andP0 equal ILU (0).We note that the studied updated preconditioners require only matrix-vector products if the application of the initial preconditionerP0 relies on this operation. If P0 is taken as the ILU (0) factorization, then the complete assembly of the matrix A1 is required. P0 can be set as the inverse of the diagonal of the matrices Ai to avoid assembling this matrix. The performance of the updated preconditioners with LGMRES when P0 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.Figure 6
Number of iterations of LGMRES to solve the sequence of linear systems with the updated preconditioners, 5-rank corrections, andP0 equal to the inverse of the diagonal of Ai.Moreover, note that, for this application, Broyden’s preconditioners are better updating strategies than the spectral preconditioner, independently of the type ofP0 used. 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 P0 as a diagonal preconditioner can be suitable for huge problems with high memory demands.


#### 5.1.1. 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. Figure7 reports the number of iterations obtained with the updated preconditioners by using different maximum ranks (2,5,10) to construct the matrix Si. The columns of matrix Si are set using the previous solutions. The initial preconditioner used is the ILU (0) of the matrix A0. 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.Figure 7
Number of iterations of the LGMRES method to solve the sequence of linear systems with the updated preconditioners and different rank corrections whereP0= is ILU (0).
(a)
Bad-Broyden-sols(b)
Good-Broyden-sols(c)
Spectral-sols


### 5.2. 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.Figure 8
Geometry and materials distribution of the VVER-440 reactor.
(a)
Radial profile(b)
Axial profileThe transient is defined as follows. The control rod denoted by number 26 is ejected in the first0.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.Figure 9
Global power of the AER-DYN-001 transient.In this problem, we use the Good-Broyden preconditioner, where the columns ofSi are the five Ritz vectors associated with the smallest Ritz values of Ai 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 P0 the Block Gauss-Seidel preconditioner of A1 (BGSA1), unless the number of iterations of a linear system j is greater than 50. In such case, the new P0 is computed as the Block Gauss-Seidel preconditioner of the matrix Aj. 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 P0 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.Figure 10
Number of iterations and CPU time using the Good-Broyden preconditioner versus the Block Gauss-Seidel preconditioner.
(a)
N. iterations(b)
CPU time (s)However, as in the previous case, we propose an alternative to setP0, which can be implemented using only matrix-vector 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 Ai. For that purpose, in this work, we use ten iterations of the GMRES method preconditioned with the inverse diagonal of Ai. 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 P0 either without updating or updated with the Good-Broyden method. The total CPU time is 21.2 hours, when the P0 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 P0 (=chebyshev) is corrected.Figure 11
Number of iterations and CPU time using the Good-Broyden preconditioner versus the Chebyshev preconditioner.
(a)
N. iterations(b)
CPU time (s)Now, we compare these results (whereP0=chebyshevAi) with the previous (where P0=BGSA1). 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.Table 3
CPU time and mean number of iterations of GMRES to reach convergence with different preconditioners for the AER-DYN-001 transient.
P0=BGS(A1)P0=Chebyshev(Ai)No updateGood-BroydenNo updateGood-BroydenCPU time (h)19.918.621.29.5Avg. GMRES iterations50.343.8391.5162.6


## 5.1. 3D Langenbuch Transient

This transient was defined in [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.Figure 1
Geometry and materials distribution of the Langenbuch reactor.
(a)
Radial profile(b)
Axial profileMaterials 4 and 6 represent control rods. It has been initiated by the withdrawal of a bank of four partially inserted control rods (C1 in Figure1) 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.Figure 2
Global power of the Langenbuch transient.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 bi−Aixi2<tol. In this problem, MATLAB is used to compute the numerical results.First, we present the number of iterations to solve the first linear systemA1x1=b1 by using the LGMRES and updated preconditioners from 5 eigenvectors of A and P0A associated with the smallest eigenvalues (Table 1). The eigenvalues are computed by using the command eigs of MATLAB. The initial preconditioner P0 is set equal to the ILU (0) of the matrix A1 obtained with the MATLAB command ilu. This table shows that all the updated preconditioners with the eigenvectors of P0A1 and A improve the results obtained when the initial P0 (with no update) is used. Interestingly, one can observe that using the eigenvectors of A1 provides similar results as those obtained when the eigenvectors of P0A1 are used.Table 1
Number of iterations of LGMRES to solve the linear problemAx=b with the different updated preconditioners and choices of S. P0=ILUA.
SNo updateBad-BroydenGood-BroydenSpectraleigs(P0A)45313432eigs(A)45323435Now, we test the performance of the updated preconditioners for the sequence of linear systems. In the following, we select the initial preconditionerP0 as the ILU (0) preconditioner of matrix A1 for all computations.First, we test the type of vectors used to construct matrixSi. 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 Ai (eigs), the Ritz vectors associated with the matrix Ai from the basis obtained from GMRES (gmres), or the Ritz vectors associated with matrix Ai 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 P0 kept fixed without updating for all the linear systems in the sequence.Figure 3
Number of iterations of the LGMRES method to solve the sequence of linear systems with the updated preconditioners with rank 5 corrections andP0 equal to ILU (0).
(a)
Bad-Broyden(b)
Good-Broyden(c)
SpectralBroyden’s preconditioners (Figures3(a) and 3(b)) with the “exact” eigenvectors (eigs) provide a constant reduction of the number of iterations, with respect to using P0 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.We investigate in Figure4 the effect of using different maximum ranks (2,5,10) to construct the matrix Si to correct the initial preconditioner. In this experiment, the columns of matrix Si 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).Figure 4
Number of iterations of the LGMRES method to solve the sequence of linear systems with the updated preconditioners and different rank corrections whereP0= is ILU (0).
(a)
Bad-Broyden-sols(b)
Good-Broyden-sols(c)
Spectral-solsTable2 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(A1) 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.Table 2
CPU times (s) of LGMRES to solve the sequence of 10 linear systems with the updated preconditioners andP0 equal ILU (0). “Pi construction” collects the CPU time to construct the ILU (0) preconditioner of the first matrix and the CPU time to compute the subspace Si. “Pi application” contains the total CPU time in the computation of the 10 linear systems devoted to apply the Pi. “Total” is the total CPU time spent to obtain the solution of the 10 linear systems.
PreconditionerRankTotal its.CPU time (s)Pi constructionPi applicationPi application per it.TotalILU(A1)—4561.38.10.01828.1Bad-Broyden-sols23651.47.50.02125.653281.56.80.02123.3103131.56.60.02122.5Good-Broyden-sols23561.47.50.02123.853211.57.10.02123.1103081.57.00.02222.6Spectral-sols24011.48.10.02026.954101.58.10.02027.0104121.58.20.02027.2Figure5 compares the number of iterations of LGMRES to solve the sequence of linear systems with the updated preconditioners constructed using rank-5 corrections and P0 equal ILU (0). The Ritz vectors associated with the previous solutions are chosen as set Si. Broyden’s preconditioners (with similar results for both strategies) provide a smaller number of iterations than the spectral preconditioner.Figure 5
Number of LGMRES iterations to solve the sequence of linear systems with the updated preconditioners, 5-rank corrections andP0 equal ILU (0).We note that the studied updated preconditioners require only matrix-vector products if the application of the initial preconditionerP0 relies on this operation. If P0 is taken as the ILU (0) factorization, then the complete assembly of the matrix A1 is required. P0 can be set as the inverse of the diagonal of the matrices Ai to avoid assembling this matrix. The performance of the updated preconditioners with LGMRES when P0 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.Figure 6
Number of iterations of LGMRES to solve the sequence of linear systems with the updated preconditioners, 5-rank corrections, andP0 equal to the inverse of the diagonal of Ai.Moreover, note that, for this application, Broyden’s preconditioners are better updating strategies than the spectral preconditioner, independently of the type ofP0 used. 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 P0 as a diagonal preconditioner can be suitable for huge problems with high memory demands.


### 5.1.1. 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. Figure7 reports the number of iterations obtained with the updated preconditioners by using different maximum ranks (2,5,10) to construct the matrix Si. The columns of matrix Si are set using the previous solutions. The initial preconditioner used is the ILU (0) of the matrix A0. 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.Figure 7
Number of iterations of the LGMRES method to solve the sequence of linear systems with the updated preconditioners and different rank corrections whereP0= is ILU (0).
(a)
Bad-Broyden-sols(b)
Good-Broyden-sols(c)
Spectral-sols


## 5.1.1. 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. Figure7 reports the number of iterations obtained with the updated preconditioners by using different maximum ranks (2,5,10) to construct the matrix Si. The columns of matrix Si are set using the previous solutions. The initial preconditioner used is the ILU (0) of the matrix A0. 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.Figure 7
Number of iterations of the LGMRES method to solve the sequence of linear systems with the updated preconditioners and different rank corrections whereP0= is ILU (0).
(a)
Bad-Broyden-sols(b)
Good-Broyden-sols(c)
Spectral-sols


## 5.2. 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.Figure 8
Geometry and materials distribution of the VVER-440 reactor.
(a)
Radial profile(b)
Axial profileThe transient is defined as follows. The control rod denoted by number 26 is ejected in the first0.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.Figure 9
Global power of the AER-DYN-001 transient.In this problem, we use the Good-Broyden preconditioner, where the columns ofSi are the five Ritz vectors associated with the smallest Ritz values of Ai 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 P0 the Block Gauss-Seidel preconditioner of A1 (BGSA1), unless the number of iterations of a linear system j is greater than 50. In such case, the new P0 is computed as the Block Gauss-Seidel preconditioner of the matrix Aj. 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 P0 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.Figure 10
Number of iterations and CPU time using the Good-Broyden preconditioner versus the Block Gauss-Seidel preconditioner.
(a)
N. iterations(b)
CPU time (s)However, as in the previous case, we propose an alternative to setP0, which can be implemented using only matrix-vector 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 Ai. For that purpose, in this work, we use ten iterations of the GMRES method preconditioned with the inverse diagonal of Ai. 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 P0 either without updating or updated with the Good-Broyden method. The total CPU time is 21.2 hours, when the P0 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 P0 (=chebyshev) is corrected.Figure 11
Number of iterations and CPU time using the Good-Broyden preconditioner versus the Chebyshev preconditioner.
(a)
N. iterations(b)
CPU time (s)Now, we compare these results (whereP0=chebyshevAi) with the previous (where P0=BGSA1). 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.Table 3
CPU time and mean number of iterations of GMRES to reach convergence with different preconditioners for the AER-DYN-001 transient.
P0=BGS(A1)P0=Chebyshev(Ai)No updateGood-BroydenNo updateGood-BroydenCPU time (h)19.918.621.29.5Avg. GMRES iterations50.343.8391.5162.6


## 6. 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×105 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.



---

*Source: 3884836-2022-06-26.xml*