Accès libre

Multilevel method to compute the lambda modes of the neutron diffusion equation

A. Carreño
A. Carreño
, 
A. Vidal-Ferrándiz
A. Vidal-Ferrándiz
, 
D. Ginestar
D. Ginestar
 et 
G. Verdú
G. Verdú
   | 24 juin 2017
Applied Mathematics and Nonlinear Sciences's Cover Image
À propos de cet article
Article précédent
Article suivant
Citez
Publié en ligne: 24 juin 2017
Pages: 225 - 236
Reçu: 16 avr. 2017
Accepté: 24 juin 2017
DOI: https://doi.org/10.21042/AMNS.2017.1.00019
Mots clés
© 2017 A. Carreño, A. Vidal-Ferrándiz, D. Ginestar, G. Verdú, published by Sciendo
This work is licensed under the Creative Commons Attribution-NonCommercial-NoDerivatives 3.0 License.
Introduction

Most of the simulations performed in nuclear reactors are done using the neutron diffusion equation. This equation is an approximation of the neutron transport equation relying on the assumption that the neutron current is proportional to the gradient of the neutron flux by means of a diffusion coefficient. This approximation is analogous to Fick’s law in species diffusion and to Fourier’s law in heat transfer [1].

For given a configuration of a nuclear reactor core, it is possible to force its criticality dividing the neutron production rate by a positive number, λ, obtaining a neutron balance equation. This equation is known as the λ-modes problem,

ℒϕ=1λℳϕ,$$\begin{array}{} \displaystyle \mathcal{L}\phi=\dfrac{1}{\lambda}\mathcal{M}\phi, \label{eq:lam_pro} \end{array}$$

where ℒ is the neutron loss operator, ℳ is the neutron production operator and φ the neutron flux.

If the two energy groups approximation is used and assuming that fission neutrons are born in the fast group and there is no up-scattering, these operators have the following form,

ℒ=(−∇→(D1∇→)+Σa1+Σ120−Σ12−∇→(D2∇→)+Σa2),ℳ=(νΣf1νΣf200),ϕ=(ϕ1ϕ2),$$\begin{array}{} \displaystyle \begin{array}{*{20}{c}} {{\cal L} = \left( {\begin{array}{*{20}{c}} { - \vec \nabla \left( {{D_1}\vec \nabla } \right) + {\Sigma _{a1}} + {\Sigma _{12}}} & 0 \\ { - {\Sigma _{12}}} & { - \vec \nabla \left( {{D_2}\vec \nabla } \right) + {\Sigma _{a2}}} \\ \end{array}} \right),} & {\$ = \left( {\begin{array}{*{20}{c}} {\nu {\Sigma _{f1}}} & {\nu {\Sigma _{f2}}} \\ 0 & 0 \\ \end{array}} \right),} & {\phi = \left( {\begin{array}{*{20}{c}} {{\phi _1}} \\ {{\phi _2}} \\ \end{array}} \right)} \\ \end{array} \end{array}$$

where φ1 is the fast neutron flux and φ2, the thermal flux. D1 and D2 are the diffusion coefficients, Σa1 and Σa2, the absorption cross sections and νΣf1 and νΣf2, the average number of neutrons produced in each fission multiplied by the fission cross sections. Finally, Σ12 is the scattering cross section from group 1 to group 2. These are coefficients dependent on the position.

The fundamental eigenvalue (the one with the largest magnitude) shows the criticality of the reactor core and its corresponding eigenfunction describes the steady state neutron distribution in the core. Next sub-critical eigenvalues and their corresponding eigenfunctions are useful to develop modal methods to integrate the time dependent neutron diffusion equation [2].

Problem (1) is discretized using a high order Galerkin Finite Element Method (FEM), transforming this differential problem into an algebraic eigenvalue problem. The matrices of the problem obtained from the discretization with a FEM are large and iterative eigenvalue solvers have to be used. Different methods have been used to solve the generalized eigenvalue problem such as Krylov methods [3, 4]. However, if it is necessary to compute an accurate solution, the degree of the polynomial expansions have to be increased. As consequence, the matrix sizes in the eigenvalue problem increase considerably and efficient methods should be used. Multilevel strategies are based on the discretization of the problem using different meshes and have been successfully used for the integration of the time dependent neutron diffusion equation [5]. In this work, a multilevel method is proposed to solve the λ-modes problem combining two levels of meshes with a Modified Block Newton method [8]. The Krylov-Schur method is used to initiate the multilevel method.

The structure of the rest of the paper is as follows. In section 2 the spatial discretization of the lambda modes problem using a high order finite element method is briefly presented. In section 3, the proposed multilevel method and the eigenvalue solvers are described. Numerical results for the analysis of two different benchmark problems are presented in section 4. Finally, the main conclusions of the paper are summarized in section 5.

Spatial discretization

To solve the problem (1), a spatial discretization of the equations has to be selected. In this work, a high order continuous Galerkin finite element method is used, [6], leading to an algebraic eigenvalue problem associated with the discretization of (1) with the following structure,

Lϕ˜=1λMϕ˜,$$\begin{array}{} \displaystyle \mathbf{L}\tilde{\phi}=\dfrac{1}{\lambda}\mathbf{M}\tilde{\phi}, \label{disc_mat} \end{array}$$

where ϕ˜=(ϕ˜1,ϕ˜2)T$\begin{array}{} \displaystyle \tilde \phi = {({\tilde \phi _1},{\tilde \phi _2})^T} \end{array}$ is the algebraic vector of weights associated with the fast and thermal neutron fluxes. The matrices elements are given by

Lij=Σe=1Nt(D1∫Ωe∇→N1i⋅∇→N1jdV−D1∫ΓeN1i∇→N1jdS→)+D2∫Ωe∇→N2i⋅∇→N2jdV−D2∫ΓeN2i∇→N2jdV+(Σa1+Σ12)∫ΩeN1iN1jdV+Σa2∫ΩeN2iN2jdV−Σ12∫ΩeN2iN1jdV),Mij=Σe=1Nt(νΣf1∫ΩeN1iN1jdV+νΣf2∫ΩeN1iN2jdV),$$\begin{array}{} \displaystyle \begin{array}{*{20}{c}} {{L_{ij}}} \hfill & = \hfill & {\mathop \Sigma \limits_{e = 1}^{{N_t}} \left( {{D_1}{\smallint _\Omega }_{_e}\vec \nabla {N_{1i}} \cdot \vec \nabla {N_{1j}}{\rm{d}}V - {D_1}{\smallint _{{\Gamma _e}}}{N_{1i}}\vec \nabla {N_{1j}}{\rm{d}}\vec S} \right)} \hfill \\ {} \hfill & + \hfill & {{D_2}{\smallint _\Omega }_{_e}\vec \nabla {N_{2i}} \cdot \vec \nabla {N_{2j}}{\rm{d}}V - {D_2}{\smallint _{{\Gamma _e}}}{N_{2i}}\vec \nabla {N_{2j}}{\rm{d}}V} \hfill \\ {} \hfill & + \hfill & {\left( {{\Sigma _{a1}} + {\Sigma _{12}}} \right){\smallint _\Omega }_{_e}{N_{1i}}{N_{1j}}{\rm{d}}V + {\Sigma _{a2}}{\smallint _\Omega }_{_e}{N_{2i}}{N_{2j}}{\rm{d}}V} \hfill \\ {} \hfill & - \hfill & {{\Sigma _{12}}{\smallint _\Omega }_{_e}{N_{2i}}{N_{1j}}{\rm{d}}V),} \hfill \\ {{M_{ij}}} \hfill & = \hfill & {\mathop \Sigma \limits_{e = 1}^{{N_t}} \left( {\nu {\Sigma _{f1}}{\smallint _\Omega }_{_e}{N_{1i}}{N_{1j}}{\rm{d}}V + \nu {\Sigma _{f2}}{\smallint _\Omega }_{_e}{N_{1i}}{N_{2j}}{\rm{d}}V} \right),} \hfill \\ \end{array} \end{array}$$

where Ni is the prescribed shape function for the ith node. For simplicity, the shape functions used are part of Lagrange finite elements. Ωe, with e = 1,...,Nt, are the reactor subdomains (cells) in which the reactor domain is divided. In the same way, Γe are the corresponding subdomain surfaces which are part of the reactor frontier. The boundary conditions considered, for the reactors analyzed in this work, are zero-flux conditions, and then, the shape functions of the corresponding nodes are fixed to zero. More details on the spatial discretization used and more general boundary conditions can be found in [4].

The finite element method has been implemented using the open source finite elements library Deal.II [7].

Multilevel method

Let us consider the eigenvalue problem (2) rewritten as

Mϕ=λLϕ$$\begin{array}{} \displaystyle M\phi=\lambda L \phi \label{eq:eig_problem} \end{array}$$

where L and M are n × n matrices arising from the discretization on a given domain Ω. For simplicity, we have dropped the tilde in the notation of the algebraic vectors and we assume that Ω is discretized on a uniform grid Ωf.

To solve (3) one typically calculates the successive approximations to the exact solution, φ, using an iterative technique, of the form

ϕnew=G(L,M,ϕold),$$\begin{array}{} \displaystyle \phi^{\text{new}}=G(L,M,\phi^{\text{old}}), \end{array}$$

where G is some expression involving the old assignment of φ. One natural way to improve this process is to obtain better initial guesses (closer to the solution) at a computational cost as low as possible. To accomplish this aim, a coarse grid Ωc is used to generate improved initial guess to the solution in the finest grid Ωf.

In the coarse level, an eigenvalue problem of the form

Mcϕc=λcLcϕc,$$\begin{array}{} \displaystyle M^{{c}}\phi^{c}=\lambda^{c} L^{c} \phi^{c}, \label{eq:eig_problem_{c}} \end{array}$$

is solved, where Mc and Lc are the matrices associated with the spatial discretization Ωc.

To use a vector φc as an initial guess for the problem associated with the finer grid Ωf, we define the linear operator called interpolation operator Icf:Ωc→Ωf$\begin{array}{} \displaystyle I_{c}^{f}:\Omega^{c}\rightarrow \Omega^{f} \end{array}$. Thus, ϕf=Icfϕc$\begin{array}{} \displaystyle \phi^{{f}}=I_{c}^{f}\phi^{{c}} \end{array}$. Similarly, to represent a vector φf on a coarser grid Ωc, we use the restriction operator Ifc:Ωf→Ωc$\begin{array}{} \displaystyle I_{f}^{c}:\Omega^{f}\rightarrow \Omega^{c} \end{array}$. In the restriction process, besides the coarsening of the initial spatial discretization, the materials and their corresponding cross section must be redefined in the coarse mesh with a homogenization method. In each coarse cell, d that is equal to the join of cells d1,...,dn, i.e. d=∪i=1mdi$\begin{array}{} \displaystyle d=\cup_{i=1}^md_i \end{array}$, the value of each cross section XSd is computed with a volume averaged, as,

XSd=1V∑i=1mViXSdi,$$\begin{array}{} \displaystyle XS_d=\dfrac{1}{V}\sum_{i=1}^m V_i XS_{d_i}, \end{array}$$

where XSdi is the value of the cross section in cell di. V is the volume of the coarse cell d and Vi the volume of the cell di.

The multilevel method can be summarized in the following steps:

Restrict the problem in a fine mesh to a coarse mesh ϕc=Ifcϕf$\begin{array}{} \displaystyle \phi^{{c}}=I_{f}^{c}\phi^{{f}} \end{array}$.

Solve with Krylov-Schur method the eigenvalue problem Mcϕc = λcLcϕc.

Interpolate the solution of coarse mesh ϕ0f=Icfϕc$\begin{array}{} \displaystyle \phi_0^{{f}}=I_{c}^{f}\phi^{{c}} \end{array}$.

Solve with the modified block Newton method the problem in a fine mesh Mfφf = λfLfφf using as initial guess the interpolated solution ϕ0f$\begin{array}{} \displaystyle \phi_0^{f} \end{array}$.

Krylov-Schur method

Krylov-Schur method was introduced in 2002 by Stewart [9] and can be seen as an improvement on traditional Krylov subspace methods such as Arnoldi and Lanczos for computing a subset of eigenvalues and their corresponding eigenvectors of a large and sparse matrix.

Given an ordinary eigenvalue problem

Ax=λx,$$\begin{array}{} \displaystyle Ax=\lambda x, \end{array}$$

the Krylov-Schur method is defined by generalizing the Arnoldi decomposition of order m,

AVm=VmHm+hm+1,mνm+1e*m,$$\begin{array}{} \displaystyle A{V_m} = {V_m}{H_m} + {h_{m + 1,m}}{\nu _{m + 1}}e{*_m}, \end{array}$$

where νm+1 is the result of orthonormalizing Aνm with respect to previous columns and Hm is an upper Hessenberg matrix, to obtain a so-called Krylov decomposition of order m,

AVm=VmBm+νm+1b*m+1,$$\begin{array}{} \displaystyle A{V_m} = {V_m}{B_m} + {\nu _{m + 1}}b{*_{m + 1}}, \end{array}$$

in which matrix Bm is not restricted to be upper Hessenberg and bm+1 is an arbitrary vector.

Krylov decompositions are invariant under (orthogonal) similarity transformations, so that

AVmQ=VmQ(QTBmQ)+vm+1bm+1TQ,$$\begin{array}{} \displaystyle AV_mQ=V_mQ(Q^{T}B_mQ)+v_{m+1}b_{m+1}^{T}Q, \end{array}$$

with QTQ = I, is also a Krylov decomposition. In particular, one can choose Q in such a way that Sm = QTBmQ is in a (real) Schur form, that is, upper (quasi-)triangular with the eigenvalues in the 1 × 1 or 2 × 2 diagonal blocks. This particular class of relation, called Krylov-Schur decomposition, can be written in block form as

A[V˜1V˜2]=[V˜1V˜2][S11S120S22]+vm+1[b˜1Tb˜2T],$$\begin{array}{} \displaystyle A\begin{bmatrix} \tilde{V}_1& \tilde{V}_2 \end{bmatrix}=\begin{bmatrix} \tilde{V}_1& \tilde{V}_2 \end{bmatrix}\begin{bmatrix} S_{11}&S_{12}\\0& S_{22} \end{bmatrix}+v_{m+1}\begin{bmatrix} \tilde{b}^{T}_1 & \tilde{b}^{T}_2 \end{bmatrix}, \end{array}$$

and has the nice feature that it can be truncated, resulting in a smaller Krylov-Schur decomposition,

AV˜1=V˜1S11+vm+1b˜1T,$$\begin{array}{} \displaystyle A\tilde{V}_1=\tilde{V}_1S_{11}+v_{m+1}\tilde{b}^{T}_1, \end{array}$$

that can be extended again to order m.

When addressing a generalized eigenvalue problem, Mx = λLx, a reformulation of the problem is used to obtain an ordinary eigenvalue problem. In particular, L−1Mx = λx, is solved, if L is non-singular. The iterative methods to solve the eigenvalue problem rely on the matrix-vector multiplication. The computation of L−1Mx is done solving linear systems. These linear systems are solved with BiCGStab method, together a Cuthill-McKee reordering to optimize the fill-in and the incomplete LU factorization for preconditioning the matrices.

To solve the ordinary and generalized eigenvalue problems by Krylov-Schur method, the library SLEPc [10] has been used.

Modified Block Newton method for generalized eigenvalue problems

Given a partial generalized eigenvalue problem of the form

MV=LVΛ,$$\begin{array}{} \displaystyle MV = LV\Lambda , \end{array}$$

where V∈ℝn×q$\begin{array}{} \displaystyle V \in {^{n \times q}} \end{array}$ is a matrix of eigenvectors and Λ∈ℝq×q$\begin{array}{} \displaystyle \Lambda \in {^{q \times q}} \end{array}$ is a diagonal matrix whose diagonal elements are the dominant eigenvalues, n denotes the degree of freedoms in finite element method and q the number of desired eigenvalues. It is assumed that the eigenvectors can be factorized as

V=ZS,$$\begin{array}{} \displaystyle V = Z S, \end{array}$$

where ZTZ = Iq. Problem (8) can be rewritten as

MV=LVΛ⇒MZS=LZSΛ⇒MZ=LZSΛS−1⇒MZ=LZK.$$\begin{array}{} \displaystyle MV=LV\Lambda \Rightarrow MZS=LZS\Lambda \Rightarrow MZ=LZS\Lambda S^{-1}\Rightarrow MZ=LZK. \label{sis:rel} \end{array}$$

This problem is undetermined since the eigenvectors are defined up to a constant. To determine the problem, the biorthogonality condition WTZ = Iq is introduced, where W is a fixed matrix of rank q. Newton’s method is used to solve the following problem

F(Z,Λ):=(MZ−LZKWTZ−Iq)=(00) .$$\begin{array}{} \displaystyle F(Z,\Lambda) := \left( \begin{array}{c} M Z-LZ K \\W^{T}Z-I_q \end{array} \right) = \left( \begin{array}{c} 0 \\ 0 \end{array} \right) \ . \label{eq:MatFunctW} \end{array}$$

Thus, a new iterated solution arises as,

Z(k+1)=Z(k)−ΔZ(k),    K(k+1)=K(k)−ΔK(k),$$\begin{array}{} \displaystyle Z^{(k+1)}=Z^{(k)}-\Delta Z^{(k)},\qquad K^{(k+1)}=K^{(k)}-\Delta K^{(k)}, \label{inc_new} \end{array}$$

where ΔZ(k) and ΔK(k) are solutions of the system

{MΔZ(k)−LΔZ(k)K(k)−LZ(k)ΔK(k)=MZ(k)−LZ(k)K(k),WTΔZ(k)=WTZ(k)−Iq,$$\left\{ \matrix{ {M\Delta {Z^{(k)}} - L\Delta {Z^{(k)}}{K^{(k)}} - L{Z^{(k)}}\Delta {K^{(k)}} = M{Z^{(k)}} - L{Z^{(k)}}{K^{(k)}},} \cr {{W^T}\Delta {Z^{(k)}} = {W^T}{Z^{(k)}} - {I_q},}} \right. $$

that is obtained substituting (12) into (11) and removing second order terms.

The system (13) is coupled, since the matrix K(k) is not necessarily a diagonal matrix. To decouple the system, the Modified Block Newton method applies two previous steps. The first step consists of an orthogonalization to the matrix Z(k) using the modified Gram-Schmidt Orthogonalization. Once Z(k) is a matrix whose columns are orthogonal, i.e., Z(k)TZ(k) = Iq, as a second step, a Rayleigh-Ritz procedure for generalized eigenvalue problems is applied [11, 12], which consists of obtaining the eigenvectors S(k) and their corresponding eigenvalues Λ(k) that satisfy

Z(k)TMZ(k)S(k)=Z(k)TLZ(k)S(k)Λ(k).$$\begin{array}{} \displaystyle Z^{(k)^{T}} MZ^{(k)}S^{(k)}=Z^{(k)^{T}}LZ^{(k)}S^{(k)}\Lambda^{(k)}. \label{sis:rrgen} \end{array}$$

Defining Z¯(k):=Z(k)S(k)$\begin{array}{} \displaystyle \bar{Z}^{(k)}:=Z^{(k)}S^{(k)} \end{array}$, we have, from (14), that Λ(k) is a diagonal matrix whose elements, λi are the Ritz values and Z¯(k)$\begin{array}{} \displaystyle \bar{Z}^{(k)} \end{array}$ are the approximated Ritz eigenvectors, satisfying the equation

Z(k)T(MZ¯(k)−LZ¯(k)Λ(k))=0.$$\begin{array}{} \displaystyle Z^{(k)^{T}}( M\bar{Z}^{(k)}-L\bar{Z}^{(k)}\Lambda^{(k)})=0. \end{array}$$

At each iteration, the matrix W in equation (13) is chosen as the previous approximation for the invariant subspace, that is, W=Z¯(k)$\begin{array}{} \displaystyle W=\bar{Z}^{(k)} \end{array}$. Using the definition of K(k) on (10), system (13) is decoupled into the q linear systems

(M−Lλi(k)LZ¯(k)Z¯(k)†0)(Δz¯i(k)−Δλi(k))=(Mz¯i(k)−Lz¯i(k)λi(k)0)$$\begin{array}{} \displaystyle \begin{pmatrix} M-L\lambda_i^{(k)} & L\bar{Z}^{(k)}\\ \bar{Z}^{(k)^{\dagger}} &0 \end{pmatrix} \begin{pmatrix}\Delta \bar{z}_i^{(k)} \\ -\Delta \lambda_i^{(k)} \end{pmatrix}= \begin{pmatrix} M\bar{z}_i^{(k)}-L\bar{z}_i^{(k)}\lambda_i^{(k)} \\ 0 \end{pmatrix} \label{sis:corrnewton} \end{array}$$

with i = 1,...,q, and Δz¯i(k)$\begin{array}{} \displaystyle \Delta \bar{z}_i^{(k)} \end{array}$ the i-th column of ΔZ¯(k)$\begin{array}{} \displaystyle \Delta \bar{Z}^{(k)} \end{array}$. Vectors Z(k + 1) are updated according to equation (12) and the eigenvalues λ(k)i are obtained from the small problem (14).

The Modified Block Newton method can be summarized as in Algorithm 1. In this Algorithm, a deflation technique is implemented which consists of working only with the non-converged eigenvectors, thus “locking” those eigenvectors that have already converged. Furthermore, to solve the system (16) a dynamic procedure is used to set the tolerance. Also, an expansion of the dimension of this method can be considered. In this case, the initial approximation has a size equal to the dimension, that is denoted in the Algorithm by d. When the dimension is not indicated this value is set to d = q (number of desired eigenvalues).

MBNM

Input: Initial approximation U = [u1,...,ud] for the eigenvectors.
Output: Diagonal matrix of eigenvalues Λ and matrix V with the eigenvectors as its columns.
 1:Orthonormalize U▹ MODIFIED GRAM-SCHMIDT
 2:Compute the Ritz approximations U, W▹ RAYLEIGH-RITZ GEN.
 3:conv=0 (conv indicates the number of converged eigenvalues)
 4:while conv < q do
 5:   Compute ΔU=[Δu1λ,…,Δuendλ]$\begin{array}{} \displaystyle \Delta U =[\Delta u^{\lambda}_{1},\dots,\Delta u^{\lambda}_{\text{end}}] \end{array}$▹ CORRECTION NEWTON EQ. (16)
 6:   U = U − ΔU
 7:   Orthonormalize(U)▹ MODIFIED GRAM-SCHMIDT
 8:   Compute the Ritz approximatisons U, W▹ RAYLEIGH-RITZ GEN.
 9:   fori =conv → q do
10:       if‖Mui−wiLui‖‖ui‖<toli${{M{u_i} - {w_i}L{u_i}} \over {{u_i}}} < to{l_i}$then▹ CONVERGENCE CHECKING
11:         Λ = [Λ,wi]▹ LOCKING CONVERGED EIGENPAIRS
12:         V = [V,ui]
13:       else
14:         break
15:       end if
16:   end for
17:   conv=i
18:   W = [wconv,...,wend]
19:   U = [uconv,...,uend]
20:end while

Numerical results

In this section, two different benchmark problems are presented to illustrate the efficiency of the multilevel scheme proposed in this paper.

The options for multilevel method have been: Krylov-Schur method to solve the coarse mesh with tolerance equal to 10−3. The number of eigenvalues required in each method has been 4. The dimension of the Krylov subspace set has been equal to 11, which is a reasonable value for the number of eigenvalues required. In the modified Block Newton algorithm, the tolerance has been set to 10−6. The linear systems of the Newton correction (Eq. (16)) have been solved with GMRES method and incomplete LU factorization for preconditioning. The tolerances in the dynamic procedure have been, in the first iteration 10−2, 10−3 in the second one, 10−5 in the third one and in the following iterations, 10−8.

For each problem, the computational time of the multilevel method is compared with Krylov-Schur method for the finest mesh. The Krylov-Schur tolerance to solve directly the fine mesh has been 10−6. The dimension has been set equal to the one used for the coarse mesh, 11. The time to solve the coarse mesh (CM) has been added to observe the amount of time in multilevel method spent in this part. The improvement in time of the multilevel method (improv.) has been computed dividing the difference between the CPU time with Krylov-Schur and multilevel method by the time with Krylov-Schur.

Cuboid reactor

This theoretical reactor consists of 3D cuboid with different materials placed heterogeneously in the different cells. Even though this problem is completely theoretical, it is relevant to analyze the computational time of multilevel scheme in a mesh with a large difference between the coarse and the fine meshes (8 cells in fine mesh are joint to form 1 cell in the coarse one). It has been designed with 8000 different assemblies (20 cells in each axis). The dimensions are 20 cm × 20 cm × 20 cm. The fine and coarse meshes considered in this problem are represented in Figure 1. Also, the number of cells in each mesh (n. cells) are indicated.

Fig. 1

Meshes for 3D cuboid reactor

The computational times observed in this problem and the improvement of multilevel method with respect to Krylov-Schur method are shown in Table 1. The results are computed for finite element degrees (FED) equal to 2 and 3. In DoFs, big differences are observed between fine and coarse meshes. The matrices in the fine mesh are approximately seven times bigger than the matrices in the coarse one. So, the time to solve the problem in the coarse mesh for the multilevel method is relatively small. The number of iterations in MBNM to obtain the tolerances desired have been 4 and 5, for FE= 2 and FE= 3, respectively. This indicates that the solution obtained in the coarse mesh is close to the solution in the fine mesh.

Computational times of cuboid reactor

DoFsCPU Time
FEDFineCoarseKrylov-SchurMultilevel (CM)Improv.
213784218522106s87s (3.8s)17%
345396259582794s555s (41s)30%

If we compare the computational times between Krylov-Schur and multilevel method, we observe a difference of 19s with FED= 2 and a difference of 239s with FED= 3, that are improvements of 17% and 30% with respect to Krylov-Schur method. So, the improvement, in this case, increases when the degree of polynomials used in the finite element method is increased. Note that if the dimension in Krylov subspace is increased these differences are reduced. However, there are not improvements in time with expansions of dimension in the MBNM.

Ringhals C9 Case

For a practical application of the multilevel method, we have chosen a configuration of Ringhals rector. Particularly, we have chosen the C9 point of the Ringhals I stability benchmark, which corresponds to a point of operation that degenerated in an out-of-phase oscillation [13]. Its geometry (fine mesh) and the coarse mesh considered are represented in Figure 2.

Fig. 2

Meshes for Ringhals reactor

The computational times computed with multilevel method and Krylov-Schur method are shown in Table 2.

Computational times of Ringhals reactor

DoFsCPU Time
FEDFineCoarseKrylov-SchurMultilevel (CM)its. MBNMImprov.
2334510125730237s180s (35s)524%
311061804119921864s1596s (463s)414%

As in the previous case, the finite element degrees (FED) considered have been 2 and 3. In DoFs, we observe that the relation in the number of cells between the coarse and the fine mesh is a little higher than 1/3. The number of iterations in MBNM to obtain the tolerances desired have been relatively low for both cases. So, the solution obtained with the coarse mesh is also close to the solution in the fine mesh. This fact can be also observed in Figure 3, where the fast flux corresponding to the first eigenvalue has been represented in the coarse and the fine mesh with FED= 2. In this figure, we observe that both eigenfunctions have a similar shape. Similar results are obtained to FED= 3. If the computational times between Krylov-Schur and multilevel method are compared, similar conclusions to the ones obtained in the cuboid case are obtained. For FED= 3, a difference of 268 seconds and an improvement of 14% are obtained and for FED= 2, a difference of 57 seconds is observed that is an improvement of 24%.

Fig. 3

Fast fluxes associated with the first eigenvalue

To make a better analysis of the multilevel method, the computational times with the residual errors, computed as

maxi‖Mϕi−λiLϕi‖‖ϕi‖,$$\mathop {{\rm{max}}}\limits_i {{||M{\phi _i} - {\lambda _i}L{\phi _i}||} \over {||{\phi _i}||}},$$

are represented in Figure 4 and Figure 5 for FED= 2 and FED= 3, respectively. The iterations are indicated with a dot. Furthermore, the computational times to solve the coarse problem and to build the matrices of fine problem are plotted. In these figures, it is observed that the initial approximations are near to 1 and the tolerance is reached with a small number of iterations since the block Newton method has a high rate of convergence. The time used in each iteration is large, being smaller in the first iterations, due to the dynamic tolerance implemented. The computational time to solve the coarse mesh represents an important part of the time in multilevel method, so a greater coarsening of initial mesh means less time in this part. The time to interpolate the solution and to build the matrices to the fine problem is very small, in both cases.

Fig. 4

CPU time of multilevel method for Ringhals reactor with FED= 2

Fig. 5

CPU time of multilevel method for Ringhals reactor with FED= 3

Conclusions

The λ-modes problem has been discretized using a high order finite element method. A multilevel procedure based on a block Newton method has been proposed to solve the algebraic eigenvalue problem obtained from the discretization. This method has been based on using two grids and a modified block Newton method.

To investigate the performance of the methodology two benchmark problems have been analyzed: a theoretical cuboid and the realistic Ringhals reactor. In both problems we have used different finite element degrees, and the computational times necessary to solve the problem with the multilevel method have been smaller than the ones used by Krylov-Schur method. For the reactor with prismatic geometry, it is possible to make a further coarsening of the mesh and to obtain smaller number of cells, so the relative time necessary to solve the problem in the coarse mesh is smaller than in the realistic reactor. In this last type of problems, it is necessary to study different strategies to coarse the mesh according to the geometry reactor. Regarding to the modified block Newton method, it has a high rate of convergence when the initial approximation is close enough to the solution in the fine mesh.

In summary, the multilevel method is a good strategy to compute a set of λ-modes of a nuclear power reactor, specially when the geometry of the reactor allows to define a coarse mesh with a large reduction of the degrees of freedom of the problem.

eISSN:
2444-8656
Langue:
Anglais
Périodicité:
2 fois par an
Sujets de la revue:
Life Sciences, other, Mathematics, Applied Mathematics, General Mathematics, Physics
RSS Feed de la revue