#### Abstract

The neutron diffusion equation is often used to perform core-level neutronic calculations. It consists of a set of second-order partial differential equations over the spatial coordinates that are, both in the academia and in the industry, usually solved by discretizing the neutron leakage term using a structured grid. This work introduces the alternatives that unstructured grids can provide to aid the engineers to solve the neutron diffusion problem and gives a brief overview of the variety of possibilities they offer. It is by understanding the basic mathematics that lie beneath the equations that model real physical systems; better technical decisions can be made. It is in this spirit that this paper is written, giving a first introduction to the basic concepts which can be incorporated into core-level neutron flux computations. A simple two-dimensional homogeneous circular reactor is solved using a coarse unstructured grid in order to illustrate some basic differences between the finite volumes and the finite elements method. Also, the classic 2D IAEA PWR benchmark problem is solved for eighty combinations of symmetries, meshing algorithms, basic geometric entities, discretization schemes, and characteristic grid lengths, giving even more insight into the peculiarities that arise when solving the neutron diffusion equation using unstructured grids.

#### 1. Introduction

The better we engineers are able to solve the equations that model the real physical plants we design and build, the better services we can provide to our customers, and thus, general people can be benefited with better nuclear facilities and installations. The Boltzmann neutron transport equation describes how neutrons move and interact with matter. It involves continuous energy and space-dependent macroscopic cross-sections that should be known beforehand and gives an integrodifferential equation for the vectorial flux as a function of seven independent scalar variables, namely, three spatial coordinates, two angular directions, energy, and time. It represents a balance that holds at every point in space and at every instant in time. Such an equation may be tackled using a variety of approaches; one of them is a simplification that leads to the so-called neutron diffusion approximation that states that the neutron current is proportional to the gradient of the neutron flux by means of a diffusion coefficient, which is a function of the macroscopic transport cross-section. When this approximation—which is analogous to Fick’s law in species diffusion and to the Fourier expression of the heat flux—is replaced into the transport equation, a partial differential equation of second order on the spatial coordinates is obtained. Formally, the neutron diffusion equation may be derived from the transport equation by expanding the angular dependance of the vectorial neutron flux in a spherical harmonics series and retaining both the zero and one-moment terms, neglecting the contributions of higher moments [1, 2]. As crude as it may seem, this diffusion equation gives fairly accurate results when applied under the conditions in which thermal nuclear reactors usually operate. Indeed, it can be shown that in a homogeneous bare critical reactor, the neutron current is proportional to the neutron gradient for every neutron energy [3].

The energy domain is usually divided into a finite number of groups, thus transforming one partial differential equation over space and energy into several coupled equations—one for each group—containing differential operators applied only over the spatial coordinates. This resulting set of second-order PDEs is known as the multigroup neutron diffusion equation and is usually used to model, design, and analyze nuclear reactor cores by the so-called core-level calculation codes. These programs take homogenized macroscopic cross-sections (which may depend on the spatial coordinates through changes of fuel burnup, materials temperature or other properties) computed by lattice-level codes as an input and solve the diffusion equation to obtain the flux (and its related quantities such as power, xenon, etc.) distribution within the core.

Given a certain spatial distribution of materials and its properties inside a reactor core, chances are that the resulting reactor will not be critical. That is to say, in general, the rate of absorptions and leakages will not exactly overcome the neutrons born by fissions sources, and some kind of feedback—either through an external control system or by means of an inherent stability mechanism of the core [4]—is needed to keep the reactor power within a certain interval. Mathematically, this means that the transport—and thus the diffusion—equation does not have a steady-state solution. In practice, given a certain reactor configuration, its steady-state flux distribution is computed by setting all the time derivatives to zero as usual but also by dividing the fission sources by a positive value called the effective multiplication factor, which becomes also one unknown and turns the formulation into a eigenvalue problem. The largest (or smallest, depending on the formulation) eigenvalue is therefore the effective multiplication factor. If , the original configuration was subcritical, and conversely. As successive configurations make , the associated eigenfunctions approach the steady-state critical flux distribution [5].

Core-level codes traditionally use regular grids to discretize the differential operators over the space. Depending on the characteristics and symmetry of the reactor core, either squares or hexagons are used as the basic shape of the mesh. Usual discretization schemes involve cell-centered finite differences or two-step coarse-mesh/coupling coefficient methods [6–8], which are accurate and efficient enough for most applications. There are, nevertheless, certain limitations that are not inherently related to structured grids but to the way their coarseness is utilized and how they are computed and applied to the reactor geometry. These issues can be overcome by discretizing the spatial operators using a scheme which could be applied to nonstructured grids, namely, finite volumes or finite elements.

This way, unstructured grids may be used to study, analyze, and understand the numerical errors introduced by the discretization of the leakage operator with a difference-based scheme over a coarse structured grid by successively refining the mesh whilst comparing the solutions with the structured one, as depicted in Figure 1. Another example of application may be the improvement of the response matrix of boundary conditions over cylindrical surfaces, which is the case for most of the nuclear power plants cores. When a coarse structured mesh is applied to a curved geometry, there appears a geometric condition known as the staircase effect. Even though arbitrary shapes cannot be perfectly tessellated with unstructured grids, for the same number of unknowns, they better represent the geometry than structured meshes, as Figure 2 illustrates. Again, by refining the mesh and comparing solutions, the matrix responses used to set the boundary conditions on the coarse meshes can be optimized. Finally, other complex geometries such as slanted control rods (Figure 3) or irradiation chambers can be directly taken into account by unstructured grids, possibly by refining the mesh in the locations around said complexities.

(a) Coarse structured grid commonly used in diffusion codes such as in [9]. Each lattice cell is divided into a 2 2 grid for solving the neutron diffusion equation |

**(b) Discretization of the lattice cell using a fine unstructured grid as proposed in this work**

**(c) Further refinement of the unstructured grid**

**(a) Continuous two-dimensional domain**

**(b) Structured grid**

**(c) Unstructured grid**

#### 2. The Steady-State Multigroup Neutron Diffusion Equation

In the present work, we take the steady-state multigroup neutron diffusion equation for granted. That is to say, we focus on the mathematical aspects of the eigenvalue problem and make no further reference to its derivation from the transport equation nor to the validity of its application to reactor problems, as these subjects that are extensively discussed in the classic literature [1, 5, 11].

The differential formulation of the steady-state multigroup neutron diffusion equation over an -dimensional domain with groups of energy is the set of coupled differential equations: where are the unknown flux distributions and the s and the s are the macroscopic cross-sections and diffusion coefficients, which ought to be computed by a lattice-level code and taken as an input to core-level codes. However, for the purposes of fulfilling the objectives of this work, we take the macroscopic cross-sections as functions of the spatial coordinate which is known beforehand. The coefficient represents the fission spectrum and is the fraction of the fission neutrons that are born into group . As stated above, the -fissions are divided by a positive effective multiplication factor and all the time derivatives, that is, the right hand of (1) is set to zero. We note that, in order to define a consistent nomenclature, we use the absorption cross section of the group , which is equal to the total cross section minus the self-scattering cross section . Therefore, the sum over the scattering sources excludes the term corresponding to . If we had used the total cross section in the second term of (1), we would have had to include the self-scattering term as a source.

If the cross sections depend only in an explicit way on the spatial coordinate , then (1) is linear. If, as is the general case, the cross sections depend on through the flux itself—such as by means of the xenon concentration or by local temperature distributions—then (1) is nonlinear. Nevertheless, this general case can be solved by successive linear iterations so the basic problem can be regarded as being purely linear.

It should be noted that, if at least one of the diffusion coefficients is discontinuous over space, then the divergence operator is not defined at the discontinuity points. Therefore, the differential formulation—also known as the strong formulation—is not complete when there exist material discontinuities that involve the diffusion coefficients. At these material interfaces, the differential equation has to be replaced by a neutron current conservation condition: where the plus and minus sign denote both sides of the interface. As the diffusion coefficients are different, the resulting flux distribution ought to have a discontinuous gradient at the interface in order to conserve the current.

When transforming the strong formulation into a weak formulation—not just into an integral formulation—both (1) and (2) can be taken into account by a single expression. In effect, let be arbitrary functions of for . Multiplying each of the equations (1) by , integrating over the domain , and applying Green’s formula [12] to the leakage term, we obtain

These coupled equations should hold for any arbitrary set of functions . Making an analogy between (3) and the principle of virtual work for structural problems, we call these functions virtual fluxes. This formulation does not involve any differential operator over the diffusion coefficient and yet, at the same time, can be shown to be equivalent to the strong formulation given by (1).

In any case, both formulations involve the computation of unknown functions of and one unknown real value . Except for homogeneous cross sections in canonical domains , the multigroup diffusion equation needs to be solved numerically. As discussed below, any nodal-based discretization scheme replaces a continuous unknown function by discrete values for . If we arrange these unknowns into a vector such as then the continuous eigenvalue problem—in either formulation—can be transformed into a generalized matrix eigenvalue/eigenvector problem casted in either of the following forms: where and are square matrices. We call the fission matrix, as it contains all the -fission terms which are the ones that we artificially divided by . The rest of the terms are grouped into the removal or transport matrix , which includes the rest of the neutron-matter interaction mechanisms. It can be shown [11] that, for any real set of cross sections, exists and that the pairs of eigenvalue/eigenvector solutions of (5) satisfy that(1)there is a unique real positive eigenvalue greater in magnitude than any other eigenvalue, (2)all the elements of the eigenvector corresponding to that eigenvalue are real and positive, (3)all other eigenvectors either have some elements that are zero or have elements that differ in sign from each other.

##### 2.1. Boundary Conditions

Being a differential equation over space, the neutron diffusion equation needs proper boundary conditions to conform a properly defined mathematical problem. These can be imposed flux (Dirichlet), imposed current (Neumann), or a linear combination (Robin). However, due to the fact that in the linear problem in absence of external sources—such as (1) or (3)—the problem is homogeneous; if is a solution, then any multiple is also a solution. That is to say, the eigenvectors are defined up to a multiplicative constant whose value is usually chosen as to obtain a certain total thermal power. Thus, the prescribed boundary conditions should also be homogeneous and be defined up to a multiplicative constant. Therefore, the allowed Dirichlet conditions are zero flux at the boundary; the Neumann conditions should prescribe that the derivative in the normal direction should be zero (i.e., symmetry condition), and the Robin conditions are restricted to the following form: being the outward unit normal to the boundary of the domain .

##### 2.2. Grids and Schemes

One way of solving the neutron diffusion equation—and in general any partial differential equation over space—is by discretizing the differential operators with some kind of scheme that is applied over a certain spatial grid. Given an -dimensional domain, a grid defines a partition composed of a finite number of simple geometric entities. In structured grids, these elementary entities are arranged following a well-defined structure in such a way that each entity can be identified without needing further information than the one provided by the intrinsic structure. On the other hand, the geometric entities that compose an unstructured grid are arranged in an irregular pattern, and the identification of the elementary entities needs to be separately specified, for example, by means of a sorted list. For instance, Figure 2(b) shows a structured grid that approximates the continuous domain of Figure 2(a). Each square may be uniquely identified by means of two integer indexes that indicate its relative position in each of the horizontal and vertical directions. In the case of Figure 2(c) that shows an unstructured grid, there is no way to systematically make a reference to a particular quadrangle without any further information.

Almost all of the grid-based schemes—which are known as nodal schemes, which are to be differentiated from modal schemes based on series expansions—are based in either the finite differences, volumes, or elements method. Finite differences schemes provide the most simple and basic approach to replace differential (i.e., continuous) operators by difference (i.e., discrete) approximations. However, they are not suitable for unstructured meshes and may introduce convergence problems with parameters that are discontinuous in space, which is the case for any reactor core composed of at least two different materials. Moreover, boundary conditions are hard to incorporate and usually give rise to incorrect results.

Methods of the finite volumes family involve the integration of the differential equation over each elementary entity, applying the divergence theorem to transform volume integrals into surface integrals and providing a mean to estimate the fluxes through the entity’s surface using information contained in its neighbors. In this context, each elementary entity is called a cell, and finite volumes methods give the mean value of each of the group fluxes at the th cell, which may be associated to the value of where is the location of the cell barycenter. Being an integral-based method, spatial-discontinuous parameters are handled more efficiently than in finite differences, and boundary conditions can be easily incorporated as forced cell fluxes. Nonetheless, even though these methods may be applied to unstructured grids, the estimation of the fluxes on the cells’ surfaces is usually performed by using some geometric approximations that may lose validity as the quality of the grid is worsened. Moreover, the simple integral approach cannot take into account discontinuous diffusion coefficients, so when two neighbors pertain to different materials the flux has to be computed in a certain particular way to conserve the neutron current according to (2).

Finally, finite elements methods rely on a weak formulation of the differential problem similar to (3) that maintains all the mathematical characteristics of the original strong formulation plus its boundary conditions. The method is based on shape functions that are local to each elementary entity—now called element, defined by nodes as corners—and on finding a set of nodal values such that a certain condition is met, which is usually that the residual of the solution has to be orthogonal to each of the shape functions. This condition is known as the Galerkin method and implies that the error committed in the approximate solution of the continuous problem is confined into a small subset of the original vector space of the continuous functions. These mathematical properties make finite elements schemes very attractive. However, these features depend on a large number of integrations that ought to be performed numerically, so a computational effort/desired accuracy tradeoff has to be considered. Not only does the finite elements method give the values of the flux at the th node but also the shape functions provide explicit expressions to interpolate and to evaluate the unknown functions at any location of the domain . Boundary conditions are divided into essential and natural. The first group comprises the Dirichlet boundary conditions which are satisfied exactly—within the precision of the eigenvalue problem solver—by the obtained solution. The latter include the Neumann and the Robin conditions that are satisfied only approximately by the derivatives of the interpolated solution through the shape functions with finer meshes giving better agreement with the prescribed values.

The same unstructured grid may be used either for the finite volumes or for the finite elements method. In the first case, the unknowns are the mean value of the fluxes over each cell, whilst in the latter the unknowns are the fluxes evaluated at each node. Therefore, the number of unknowns is different for each method, even when using the same grid. Figure 4 shows this situation, and in Section 3.1, we further illustrate these differences. A fully detailed mathematical description of the actual algorithms for both volumes and elements-based proposed discretizations can be found in an academic monograph written by the author of this paper [13].

#### 3. Results

We now proceed to show two illustrative results that are to be taken as an overview of the possibilities that unstructured grids can provide in order to tackle the multigroup neutron diffusion problem. The examples are two-dimensional problems, as they contain some of the complexities a real three-dimensional reactor posse, yet the reported results are not so complicated as to be easily understood and analyzed. In particular, we state some basic differences between the finite volumes and the finite elements methods by solving a two-group homogeneous bare reactor with the Robin boundary conditions over the very same grid, although a rather coarse one, so the differences can be observed directly into the resulting figures. We then solve the classical two-dimensional LWR problem, also known as the 2D IAEA PWR benchmark. Not only do we show again the differences between the finite volumes and elements formulation but also we solve the problem using different combinations of meshing algorithms, basic shapes, and characteristic lengths of the mesh.

To solve the two examples shown below, we employed the milonga code, which was written from scratch by the author of this paper and is currently still under development within his ongoing PhD thesis. There exists a first public release [14] under the terms of the GNU General Public License—that is, it is a free software—that can only handle structured grids. There is a second release being prepared, whose main relevance is that it can work with nonstructured grids as well, which is the main feature of the code. It uses a general mathematical framework—coded from scratch as well—which provides input file parsing, algebraic expressions evaluation, one- and multidimensional interpolation of scattered data, shared-memory access, numerical integration facilities, and so forth. The milonga code was designed with four design basis vectors in mind—which are thoroughly discussed in the documentation [15]—which includes the type of problems it can handle, the code scalability, which features are expected, and what to do with the obtained results.

It works by first reading an input file that, using plain-text English keywords and arguments, defines the number of spatial dimensions, the number of group energies, and optionally a mesh file. Currently, only grids generated with the code gmsh [16] are only supported, mainly because it is also a free software and it suits perfectly well milonga’s design basis in the sense that the continuous geometry can be defined as a function of a number of parameters. Afterward, the physical entities defined in the grid are mapped into materials with macroscopic cross sections, which may depend on the spatial coordinates by means of intermediate functions such as burn-up or temperatures distributions, which in turn may be defined by algebraic expressions, by interpolating data located in files, by reading shared-memory objects or by a combination of them. In the same sense, boundary conditions are applied to appropriate physical entities of dimension . The problem matrices and are then built and stored in an appropriate sparse format using the free PETSc [17, 18] library, and the eigenvalue problem is solved using the free SLEPc [19] library. The results are stored into milonga’s variables and functions, which may be evaluated, integrated—usually using the free GNU Scientific Library [20] routines—and of course written into appropriate outputs. Milonga can also solve problems parametrically and be used to solve optimization problems. As stated above, the code is a free software so corrections and contributions are more than welcome.

##### 3.1. A Coarse Bare Homogeneous Circle

Figures 5 and 6 show the results of solving a bare two-dimensional circular homogeneous reactor of radius over an unstructured grid which is deliberately coarse, using the finite volumes method and the finite elements method, respectively. Two group energies were used and a null-flux boundary condition was fixed at the external surface. Figures 5(a) and 6(a) compare the obtained fast flux distributions. In the first case, the numerical solution provides mean values for each neutron flux group in each cell, while in the latter, the solution is computed at the nodes, and continuous functions are evaluated by means of the shape functions used in the formulation. Figures 5(b) and 6(b) illustrate the fast fluxes unknowns and its relative position in space. It can be noted that the mesh coarseness gives results that differ substantially in both cases. Finally, the structure of the sparse eigenvalue problem matrices is shown for each case with blue, red, and cyan representing positive, negative, and explicitly inserted zero values. In the finite volume case, there are 184 unknowns (92 cells 2 groups), while in the second case there are 218 unknowns (109 nodes 2 groups). The volumes’ fission matrix is almost diagonal: it has a bandwidth equal to the number of energy groups, which in this case is two. The off-diagonal values appear right next to the diagonal elements because of the chosen ordering of the unknowns in the flux vector . Other orderings may be used, but the rate of convergence of the eigenvalue problem can be deteriorated. On the other hand, the elements’ fission matrix has a nontrivial structure, because the grid is unstructured, and the net fission rate at each element depends on the fluxes at the nodes whose location inside the matrix depends in turn on how the grid was generated. This effect is similar to the one that appears in structural analysis where mass matrices need to be lumped in order to simplify transient computations [21]. The diagonal block of element’s matrix and the null block in correspond to the discrete equations that set the null-flux boundary conditions on the nodes located at the external surface. Other types of boundary conditions lead to different kinds of structures within the problem matrices.

**(a) Fast flux distribution**

**(b) Fast flux unknowns**

(c) Matrix structure |

(d) Matrix structure |

**(a) Fast flux distribution**

**(b) Fast flux unknowns**

(c) Matrix structure |

(d) Matrix structure |

In case a part of the domain contained a nonmultiplicative material such as a reflector, then there would appear sections of the fission matrix with null values in both methods, rendering singular in both methods. Care should be taken when dealing with the numerical schemes for the eigenvalue problem solution. It can be seen in the elements’ matrices a particular structure that implements the boundary conditions at the external surfaces. This structure does not appear in the volumes’ matrices because the boundary conditions appear as flux terms which are summed up over all the surfaces of each cell, so they are masked inside the volumetric discretization of the divergence and gradient operators.

As the two-group neutron diffusion equation with uniform cross sections over a circle subject to null-flux boundary conditions has an analytical solution, it is adequate to compare how the two proposed numerical schemes relate to it. In the studied problem, we ignored upscattering and fast fissions. Then, the analytical effective multiplication factor is being, , the smallest root of Bessel’s first-kind function of order zero .

Figure 7 shows how the two numerical effective multiplication factors compare to the analytical solution given by (7) as a function of the mesh refinement, indicated by the quotient between the radius of the circle and the characteristic length of the cell/elements. We can see that the computed by the finite volumes (elements) method is always greater (less) than the analytical solution. Indeed, it can be proven that for a bare one-dimensional slab this is always the case [13]. However, this result does not hold for problems with nonuniform cross sections, even in simple one-dimensional reflected reactors.

We may draw two other conclusions from Figure 7. First, that the finite element method seems to provide a better solution than the volumes-based scheme and, second, that even though the error committed tends to decrease with finer grids, its behavior is not monotonic for finite volumes. In fact, Figure 8 shows the difference between the numerical and analytical solutions using a logarithmic vertical scale, where both conclusions are even more evident. We defer the discussion of such differences between the discretization scheme until the next section. It is worth to note, however, that the fact that a finite-element-based scheme throws better results for the particular bare homogeneous circular reactor under study than the finite volumes does not imply that finite volumes ought to be discarded as a valid tool for solving the neutron diffusion equation in general cases. The combination of lattice and core-level computations is usually performed using cell-based results which when fed into node-based methods to solve the few-group neutron diffusion equation may introduce errors which potentially can lead to unacceptable solutions. Nevertheless, this analysis is far beyond the scope of this paper which focuses on solving a mathematical equation over unstructured grids.

##### 3.2. The 2D IAEA PWR Benchmark Problem

This is a classical two-group neutron diffusion problem, first designed in the early 1970s and taken as a reference benchmark for computational codes. A number of codes were used to solve either this problem or its three-dimensional formulation [22, 23], including milonga using a structured grid [24]. The original formulation can be found in the reference [10], and there is a reproduction that may be easily found online in reference [24]. The geometry consists of a one-quarter of a PWR core depicted in Figure 9, and the homogeneous macroscopic cross sections are listed in Table 1. An axial buckling term should be taken into account. The external surface should be subject to a zero incoming current condition, which may be written as

The expected results are as follows.(1)Maximum eigenvalue. (2)Fundamental flux distributions: (a)Radial flux traverses and . Note: the fluxes shall be normalized such that (b)Value and location of maximum power density. This corresponds to maximum of in the core. It is recommended that the maximum values of both in the inner core and at the core/reflector interface be given. (3) Average subassembly powers where is the volume of the th subassembly and designates the fuel subassemblies as shown in Figure 9.(4) Number of unknowns in the problem, number of iterations, and total and outer.(5) Total computing time, iteration time, IO-time, and computer used.(6) Type and numerical values of convergence criteria.(7) Table of average group fluxes for a square mesh grid of 20 20 cm.(8) Dependence of results on mesh spacing.

Even though the original problem is based on a quarter-core situation, the problem has an eighth-core symmetry which cannot be taken into account by structured grids which are the main target of the benchmark. However, nonstructured grids can take into consideration any kind of symmetry almost without loss of accuracy and at the same time reducing roughly the number of unknowns by a half and the associated computational effort needed to solve the problem by a factor of four. Answers to items can be given completely by milonga using a single input file (see Supplementary data in Supplementary Material available online at http://dx.doi.org/10.1155/2013/641863). As asked in item , how results depend not only on the mesh spacing but also on the meshing algorithm, on the grid’s elementary geometric shape, and on the discretization scheme may shed lights on the subject which may be even more interesting than the numerical results themselves.

Taking advantage of milonga’s capability of reading and parsing command-line arguments, the selection of the core geometry (quarter or eighth), the meshing algorithm (delaunay [16] or delquad [25]), the shape of the elementary entities (triangles or quadrangles), the discretization scheme (volumes or elements), and the characteristic length of the mesh can be provided at run time. Fixing five values for gives possible combinations, which we solve with successive invocations to milonga with the same input file but with different arguments from a simple script. Figures 10, 11, 12, 13, 14, 15, 16, 17, and 18 show the results corresponding to items (1)–(7) for some illustrative cases. The complete set of figures and the code used to generate them may be provided upon request. Table 2 compiles the answers to the problem for every case studied.

**(a)**

**(b)**

**(c)**

**(d)**

**(a)**

**(b)**

**(c)**

**(d)**

**(a)**

**(b)**

**(c)**

**(d)**

**(a)**

**(b)**

**(c)**

**(d)**

**(a)**

**(b)**

**(c)**

**(d)**

**(a)**

**(b)**

**(c)**

**(d)**

**(a)**

**(b)**

**(c)**

**(d)**

**(a)**

**(b)**

**(c)**

**(d)**

**(a)**

**(b)**

**(c)**

**(d)**

As the milonga code is still under development, its numerical routines are not yet fully optimized nor designed for parallel computation. Therefore, the reported times are only rough estimates and should be taken with care. The solution comprises five steps:(1)generate the grid with the requested geometry, meshing algorithm, basic shape, and characteristic length by calling to gmsh; (2)read the generated mesh; (3)build the matrices;(4)solve the eigenvalue problem; (5)compute the requested results.

The CPU time reported in Figures 10–18 is thus the sum of all of these five steps, but not the time needed to generate the figures—which in some cases with coarse meshes was considerably larger than the solution time itself. The eigenvalue problem was solved using a multilayer iterative Krylov-Schur method [26] with a tolerance relative to the matrices norm

The reported residual norm and relative error are respectively. The fields marked as outer, linear, and inner iterations refer to the number of steps needed to attain the requested tolerance in each layer of the Krylov-Schur algorithm. The computer used to solve the problem has an Intel i7 920 @ 2.67 GHz processor with 4 Gb of RAM running Debian GNU/Linux Wheezy.

When using a finite volumes-based scheme over an unstructured mesh, the solver has to gather information about which cells are neighbors and which are not. Currently gmsh does not write this kind of lists in its output files, so milonga has to explicitly solve the neighbors problem. Performing a linear search is an task, which is unacceptable for problem sizes of interest. The code uses a search based on a -dimensional tree [27], which can in principle be performed in steps. Still, for large values of , the time needed to read and parse the mesh (number two previously mentioned) is the bottleneck of the solution. This step is not needed in finite elements, although the construction of the elementary matrices involves the computation of the multidimensional Jacobians and integrals, which then have to be assembled. Again, for large , this step (number three) takes up most of the time needed to solve the problem.

The Delaunay algorithm is a standard method for generating two-dimensional grids [16] by tessellating a plane with triangles. If the mesh needs to be based on quadrangles instead, a recombination algorithm can be used to transform two adjacent triangles in one quadrangle, whenever is possible. However, for geometries which are based on rectangular shapes there exist other algorithms [25] both for meshing and for recombining the triangles that give rise to elementary entities with right angles almost everywhere, which may be a desired property of the resulting grid. For finite elements, the steps needed to build the matrices depend on the selection of triangles or quadrangles as the basic elementary geometry because the shape functions change. However, the number of unknowns is the same as the number of nodes that does not change after a recombination procedure. On the other hand, the number of unknowns in the finite volumes schemes with triangles is roughly twice as the number of unknowns with quadrangles for the same grid.

Figure 19 shows how the computed static reactivity (i.e., ) depends on the number of unknowns for each of the sixteen combinations of geometry algorithm shape scheme. The four accepted results published in the original reference [10] are included for reference, although it should be taken into account that said reactivities were computed almost forty years ago. Figure 20 shows the total wall time needed to solve the problem as a function of , while Figures 21, 22, 23, and 24 show the times needed for each of the first four steps involved in the solution, maintaining the same logarithmic scale for both the abscissas and the ordinates. Green data represent finite volumes, whilst blue bullets correspond to finite elements. Solid lines indicate quarter-core and dashed lines eighth-symmetry. Fillet bullets are results obtained by the Delaunay triangulation, and empty bullets were computed with the delquad algorithm. Finally, triangle-shaped data correspond to triangles and squares, and diamonds denote quadrangles as the basic geometry of the grid.

It can be seen that finite elements produce a much smaller dispersion of eigenvalues than finite volumes with the refinement of the mesh. This can be explained because finite volumes methods rely on a geometric condition of the mesh which may change abruptly if the meshing algorithm decides to allocate the cells in a rather different form for small changes in . Finite elements methods are less influenced by these discontinuous lay-out changes of the elements. As expected, eighth-core symmetries give almost the same results as the quarter-core geometries with roughly half the unknowns. For small problems, finite volumes run faster that finite elements because the time needed to solve the neighbor problem is negligible. When the problem size grows, this time increases and exceeds the overhead implied in the construction and assembly of the finite elements matrices. Also, at least for this configuration, it is seen that the eigenvalue problem is solved faster for finite volumes than for finite elements.

#### 4. Conclusions

Unstructured grids provide the cognizant engineer with a wide variety of possibilities to deal with the design or analysis of nuclear reactor cores. These kinds of grids can successfully reproduce continuous geometries commonly found in reactor cores such as cylinders, and therefore, not only can the diffusion equation be better approximated inside the domain of definition but also the fulfillment of boundary conditions is improved. A free computer code was written from scratch that is able to completely solve the 2D IAEA PWR Benchmark using unstructured grids for sixteen combinations of geometry, meshing algorithm, basic shape, and discretization scheme plus any value of the grid’s characteristic length by using a single input file. The complete set of input files and code—executable and source—is available either online or upon request, with comments, experiences, suggestions, and corrections being more than welcome. Further development should include tackling full three-dimensional geometries with complete thermal hydraulic feedback in order to analyze how the solutions of the coupled neutronic-thermal problem depend on the spatial discretization scheme of the neutron leakage term. Parallelization of the computation and assembly of the matrices and of the solution of the eigenvalue problem and its implementation using GPUs are also desired features to implement. A problem with direct applications that the future versions of milonga ought to solve is the analysis of how the geometry of the absorbing materials should be taken into account in structured coarse grids in order to mitigate effects such as the rod-cusp problem.

Suitable schemes for approximating the continuous differential operators by discrete matrix expressions include finite volumes and finite elements families. Finite volumes methods compute cell mean values, whilst finite elements give functional values at the grid’s nodes. In general, finite elements are less sensitive to changes in the mesh so the results they provide do not change significantly for different meshing algorithms or elementary shapes. Small problems are best solved by finite volumes as the neighbor-finding problem is faster than the process of building and assembling the eigenvalue-problem matrices. For a large number of unknowns, the process of finding which cell is neighbor of which—even based on a -dimensional tree—overwhelms the computation and assembly of elementary matrices, and finite-element methods perform better.

#### Supplementary Materials

The supplimentary material consists of the input file for the free neutronic code milonga used to solve eighty cases of the 2D PWR IAEA benchmark problem disucussed in the article.