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,
where ℒ is the neutron loss operator, ℳ is the neutron production operator and
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,
where
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
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.
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,
where
where
The finite element method has been implemented using the open source finite elements library Deal.II [7].
Let us consider the eigenvalue problem (2) rewritten as
where
To solve (3) one typically calculates the successive approximations to the exact solution,
where
In the coarse level, an eigenvalue problem of the form
is solved, where
To use a vector
where
The multilevel method can be summarized in the following steps:
Restrict the problem in a fine mesh to a coarse mesh
Solve with Krylov-Schur method the eigenvalue problem
Interpolate the solution of coarse mesh
Solve with the modified block Newton method the problem in a fine mesh
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
the Krylov-Schur method is defined by generalizing the Arnoldi decomposition of order
where
in which matrix
Krylov decompositions are invariant under (orthogonal) similarity transformations, so that
with
and has the nice feature that it can be truncated, resulting in a smaller Krylov-Schur decomposition,
that can be extended again to order
When addressing a generalized eigenvalue problem,
To solve the ordinary and generalized eigenvalue problems by Krylov-Schur method, the library SLEPc [10] has been used.
Given a partial generalized eigenvalue problem of the form
where
where
This problem is undetermined since the eigenvectors are defined up to a constant. To determine the problem, the biorthogonality condition
Thus, a new iterated solution arises as,
where Δ
that is obtained substituting (12) into (11) and removing second order terms.
The system (13) is coupled, since the matrix
Defining
At each iteration, the matrix
with
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
MBNM 1: Orthonormalize ▹ M 2: Compute the Ritz approximations ▹ R 3: conv=0 (conv indicates the number of converged eigenvalues) 4: 5: Compute ▹ C 6: 7: Orthonormalize( ▹ M 8: Compute the Ritz approximatisons ▹ R 9: 10: ▹ C 11: Λ = [Λ, ▹ L 12: 14: 15: 16: 17: conv=i 18: 19: 20:
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.
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.
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 reactorDoFs CPU Time FED Fine Coarse Krylov-Schur Multilevel (CM) Improv. 2 137842 18522 106s 87s (3.8s) 17% 3 453962 59582 794s 555s (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.
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.
The computational times computed with multilevel method and Krylov-Schur method are shown in Table 2.
Computational times of Ringhals reactorDoFs CPU Time FED Fine Coarse Krylov-Schur Multilevel (CM) its. MBNM Improv. 2 334510 125730 237s 180s (35s) 5 24% 3 1106180 411992 1864s 1596s (463s) 4 14%
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%.
To make a better analysis of the multilevel method, the computational times with the residual errors, computed as
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.
The
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