Isogeometric analysis for the multigroup neutron diffusion equation with applications in reactor physics Nuclear Energy

Isogeometric Analysis (IGA) has been applied to heterogeneous reactor physics problems using the multi- group neutron diffusion equation. IGA uses a computer-aided design (CAD) description of the geometry commonly built from Non-Uniform Rational B-Splines (NURBS), which can exactly represent complicated curved shapes such as circles and cylinders, common features in reactor design. This work has focused on comparing IGA to ﬁnite element analysis (FEA) for heterogeneous reactor physics problems, including the OECD/NEA C5G7 LWR benchmark. The exact geometry and increased basis function continuity contribute to the accuracy of IGA and an improvement over comparable FEA calculations has been observed. (cid:1) 2016 The Authors. Published by Elsevier Ltd. ThisisanopenaccessarticleundertheCCBYlicense(http://


Introduction
Development of efficient computational methods to solve reactor physics and shielding problems is an ongoing area of research. A linearised version of the Boltzmann transport equation, the neutron transport equation is solved on a seven-dimensional phase space (three spatial dimensions x; y and z, two angular dimensions h and v, energy E and time t. This high dimensionality combined with complicated reactor geometries makes solving threedimensional, heterogeneous, full core problems infeasible, even with high-performance computing (HPC) resources.
Over the past sixty years work in the field has focused on reducing the complexity of the problem to make the equations solvable on contemporary computing resources. Spatial homogenisation of complex geometrical features and utilising the much simpler diffusion approximation are still among the most widely used techniques to solve full-core calculations. Accurately modelling the physical geometry is important if high fidelity solutions to heterogeneous problems are sought. Increasingly over the past decades, the finite element method (FEM) has become one of the most commonly employed methods to represent complicated, arbitrarily shaped geometries in a wide variety of engineering and physical applications, including reactor physics. The flexibility of the FEM has made it possible to solve more complicated, heterogeneous problems and production codes have been developed to solve both first and second-order forms of the transport equation, such as ATILLA (Wareing et al., 2001) and EVENT (de Oliveira, 1986;de Oliveira et al., 1987) respectively. However, the computationally costly mesh generation process can be a significant proportion of the overall design-analysis process (Hughes et al., 2005). In this paper a new method for faithfully modelling the complex full physical geometry is proposed based on Isogeometric analysis (IGA).
The generation of a finite element mesh for finite element analysis (FEA) typically begins with a computer-aided design (CAD) model. The prevailing technology to represent CAD geometries are Non-Uniform Rational B-Splines (NURBS), which can exactly represent complicated curved shapes such as circles and other conic sections and quadric surfaces as well as other more complicated surfaces (Pigel et al., 1997). The principal aim of IGA is to unify the fields of design and analysis using a common mathematical description of the geometry in both processes, eliminating the costly finite element (FE) mesh generation step while also providing an exact geometrical description as opposed to the approximate polygonal or polyhedral meshes used for FEA (Hughes et al., 2005). IGA has been applied to a wide range of physical problems including structural mechanics (Luycker et al., 2011), computational fluid dynamics (Nielsen et al., 2011) and electromagnetics (Buffa et al., 2014). It has been observed that the exact geometrical representation of the domain can contribute significantly to the accuracy of the calculated solution and that the smooth, high continuity NURBS basis functions can better represent solutions to http://dx.doi.org/10.1016/j.anucene.2016.11.015 0306-4549/Ó 2016 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). some problems than the standard Lagrangian basis functions typically employed for FEA. This paper focuses on the application of IGA to the multi-group diffusion approximation of the neutron transport equation and builds upon the work by Hall et al. Hall et al. (2012) who applied IGA to mono-energetic fixed source pincell calculations. IGA has also recently been applied to the first-order form on the transport equation with an S N angular discretisation and discontinuous spatial formulation Owens et al., 2016. The diffusion equation has been chosen due to its relatively cheap computational cost and should provide an indicator to the efficiency of the spatial discretisation when applied to other second-order forms such as the evenparity equation or the self-adjoint angular flux (SAAF) equation. An introduction to the IGA method and its application to the neutron diffusion equation is given in Section 2. In Section 3 the NURBS IGA method is compared to FEA to determine the benefit of an exact geometrical representation. Additionally, the importance of the basis function continuity is studied. The numerical results of the study are presented for two pincell test problems; a single-group pincell with six test cases of increasingly heterogeneous materials and a more representative seven-group UOX pincell from the OECD/NEA C5G7 benchmark. In Section 4 the full quarter-core C5G7 LWR benchmark is studied. Finally, in Section 5 the conclusions of the study are summarised and suggestions for future research in the topic is given.

Isogeometric analysis (IGA)
The principle aim of IGA is to use a common mathematical description of the geometry in both the design and analysis processes. In traditional finite element techniques a CAD model is processed by mesh generation software and a polygonal approximation built of geometric primitives such as triangles, quadrilaterals etc., is generated. These geometric primitives are commonly defined as a linear combination of polynomial functions. In IGA the spatial domain is described by a linear combination of NURBS functions in which the standard finite element machinery can be applied utilising the isoparametric concept. A linear system of equations can then assembled and a set of expansion coefficients, d can be computed. The quantity of interest, which in this study is the scalar neutron flux can then be computed as a linear combination of all (N b ) NURBS shape functions, R multiplied by the expansion coefficients d i .
2.1.1. Basis functions The basis functions for a B-spline or NURBS patch are required for determining the physical shape of the curve, surface or solid. They are also required to perform numerical analysis and provide the shape functions and test functions in IGA by utilising the isoparametric concept.
To define a set of B-spline basis functions a knot vector (N) and the polynomial order (p) of the functions is required. The knot vector is a non-decreasing set of real numbers, for example N ¼ ½n 1 ; n 2 ; n 3 ; n 4 with n i 2 R. In one dimension, the knot vector defines a line in the parameter space.
For a given knot vector N, B-spline basis functions NðnÞ of a given degree p can be calculated as follows; starting with the zeroth-order: where n i is the i th knot point and n is a specified point in the onedimensional parameter space. The p ¼ 0 and p ¼ 1 basis functions are identical to standard finite element piecewise constant and piecewise linear functions respectively which is shown in Fig. 1 If the denominator of the fraction is zero, i.e. n iþp ¼ n i or n iþpþ1 ¼ n iþ1 then the whole fractional term is set to zero. The derivatives of B-spline basis functions are calculated from the p À 1 th B-spline basis function as follows; Sample linear and quadratic B-spline basis functions are displayed in Figs. 1 and 2 respectively.
Multivariate B-spline basis functions are obtained via a tensorproduct construction of the univariate functions up to the required dimension as follows: Knot vectors may be uniform if the difference between each successive knot point (the knot span) is identical or non-uniform where the width of each knot span may be arbitrary. Nonuniform knot vectors allow for knots to be concentrated in certain parts of the domain to provide increased fidelity for the shape. In one-dimension this allows for purely local refinement however, since higher dimensional B-splines are tensor products of onedimensional B-splines local refinements propagates across the full patch. Typically the first and last knot values are repeated p þ 1 times creating what is known as an open knot vector. This ensures that the functions are interpolatory at the beginning and end of the domain which is a useful property for finite element analysis.
Knots may be considered points of reduced continuity. Basis functions within a knot span will possess full C 1 continuity. Basis functions which cross a single knot point will have reduced basis function continuity of order C pÀ1 . Thus for quadratic (p ¼ 2) the first derivative of the basis functions will be continuous however the second derivative will be discontinuous. This is in contrast to the commonly used Lagrangian finite element basis where only the function itself is continuous (C 0 ). The benefit of the increased continuity will be studied in Section 5. Knot points may be repeated, introducing knot-spans of zero width. The effect of this is a further reduction in the continuity of the basis functions. For each additional repeated knot the continuity of basis functions crossing the knots is reduced by 1. If a knot value is repeated p þ 1 times then the basis functions will become discontinuous across that point Hughes et al., 2005. This may be of particular use if the numerical scheme is unstable under a fully continuous Bubnov-Galerkin discretisation, such as the first-order, hyperbolic form of the neutron transport equation.

NURBS basis functions
To construct the geometries required for reactor analysis we are required to use NURBS. NURBS basis functions are derived from the B-spline basis functions which span the same parametric space. The distinction comes that the weighting values for each control point (an additional dimension of the control grid) is required to provide a projective transformation, allowing NURBS to exactly represent conic sections such as circles and ellipses or quadric surfaces such as spheres, ellipsoids and cylinders. In one dimension this can be written in the following form: where R p i ðnÞ is the resulting NURBS basis function evaluated at point n in the parameter space, N p i is the B-spline basis function and w i are the weighting value for each basis function i. It can be seen that NURBS are defined by a division by a B-spline creating a rational function. Two-dimensional bivariate and three-dimensional trivariate NURBS basis function are derived as a tensor-product of the one dimensional B-spline basis functions for all parametric directions n and g in R 2 and n; gand f in R 3 . They are then multiplied by the corresponding control net weight, w and the weighting function W. This is given as R p;q i;j ðn; gÞ ¼ N i;p ðnÞM j;q ðgÞw i;j Wðn; gÞ in two dimensions and as R p;q;r i;j;k ðn; g; fÞ ¼ N i;p ðnÞM j;q ðgÞL k;r ðfÞ Wðn; g; fÞ in three dimensions. A selection of NURBS basis functions mapped to the physical domain for a single pincell geometry are shown in Fig. 3.

The physical, parameter and parent space
In IGA, three separate parametrisations of the domain must be considered. Each NURBS patch is parameterised as a d dimensional Cartesian grid with shape and size determined by a set of knot vectors N n ; N g and N f with parametric directions n; g and f respectively. A NURBS curve in parameter space is a one dimensional line, a NURBS surface is a two-dimensional rectangle and a NURBS solid is a three-dimensional cuboid. Knot points, lines and surfaces in one, two and three dimensions can be considered analogous to element boundaries of a finite element mesh and provide a convenient structure upon which to perform numerical integration to evaluate the integrals when assembling the linear system. The mapping of a function from the parental, integration space through to the physical space is demonstrated in Fig. 4.

Construction of a NURBS patch
To construct a NURBS curve, surface or solid in R d a set of physical points, B must be defined in R dþ1 creating a weighted control net. The d þ 1 th dimension is the corresponding weighting values. Fig. 5 demonstrates the role the control points (analagous to the nodes in finite elements) play in determining the shape of the NURBS surface.
NURBS curves CðnÞ, surfaces Sðn; gÞ and solids Vðn; g; fÞ are constructed as a linear combination of their NURBS basis functions RðnÞ; Rðn; gÞ and Rðn; g; fÞ respectively. . B-splines span across the knots when the order is greater than 1. For quadratics, each basis function is C 1 across knots and C 1 within. When an internal knot is repeated twice, the continuity between elements is reduced to C 0 . Adding an additional knot would make the basis functions discontinuous across the knot, which may be desirable for other applications such as the first-order form of the neutron transport equation to increase the numerical stability of the discretisation. Note that for linear B-splines, the basis functions are identical to those in linear finite elements, with the knot points playing an identical role to element boundaries in finite elements. The knot spans are a convenient way of subdividing the domain, integrating the basis functions using numerical quadrature and assembling the linear system.   4. The mapping of a function in the parameter space to the real space is done in two stages. Gaussian quadrature has been used to integrate the NURBS basis functions and is performed on a ½À1; 1 d parent element. The functions are then mapped through to the element in the real physical space where they can be used as the shape functions for finite element analysis. Fig. 5. The shape of a NURBS patch is determined by a set of coordinates, B in physical space (highlighted in blue) that creates a control polygon. In the above three NURBS surfaces in R 3 have been defined and by altering the position of a control point the surface is perturbed. In this manner and by increasing the number of degrees of freedom complex objects can be modelled exactly.
where i; j and k are the indices of the basis functions in the n; g and f parametric dimension. p; q and r denote the polynomial order of the univariate NURBS basis function used to create the multivariate NURBS shapes.

Discretisation of the Neutron Diffusion Approximation
The standard Galerkin method can be used to discretise the diffusion equation in an identical way to traditional finite element analysis by utilising the NURBS functions described previously instead of the standard polynomial basis functions. The multidimensional, mono-energetic, steady state neutron diffusion equation can be written as where X is the spatial domain and is bounded by the domain boundary C; /ðrÞ is the scalar neutron flux, DðrÞ is the diffusion coefficient, R f ðrÞ is the macroscopic fission cross section, m is the average neutrons produced per fission and k eff is the effective neutron multiplication factor to ensure a steady state solution. The equation is subject to the boundary conditions r/ðrÞ Á n ¼ 0 r 2 C R and where C R and C V are the reflective and vacuum boundaries respectively. n is the outward normal along the boundary. As is well known from the finite element literature the weak form of the equation can be obtained by multiplying Eq. 13 by an arbitrary test function vðrÞ 2 W 1 ðXÞ, integrating over all space X and applying the divergence theorem to the diffusion term of the equation giving ðrv; Dr/Þ þ ðR a /; vÞ þ hDr/ Á n; vi ¼ where ðu; vÞ R X uðrÞvðrÞdX and hu; vi R C uðrÞvðrÞdC are the inner products of two functions over the spatial domain and boundary respectively. In a manner identical to traditional finite element methods, the discretised equations are obtained by substituting Eq. 1 into Eq. 15 and selecting the test functions v to be the same as the shape functions using the Bubnov-Galerkin method. The resulting linear system can be solved using the conjugate gradient method preconditioned by incomplete-LU factorisation.

Numerical results
The main advantages of IGA compared to traditional FEA stem from two unique properties of the NURBS based model; the exact representation of the geometry and the increased in-patch basis function continuity. In the following section the potential advantages of these properties will be systematically studied. In order to facilitate this study two simple test problems have been selected which cover a broad spectrum of physical behaviour. In this section these test cases will be described followed by the computational results and subsequent analysis.
Since the potential future application of this method is expected to be pin-heterogeneous core models for reactor analysis, both test problems are based on the geometry of a Light Water Reactor (LWR) pin cell. The first test problem has been designed to study the relative importance of the approximation in the geometry and the solution function, while the second test case has been selected to demonstrate the method on a problem representative of the fuel pins found in a typical LWR. In both cases the cell disadvantage factor will be used as the quantity of interest, where / mod and / fuel are defined as the average flux in the cell moderator and cell fuel pin respectively. In addition, in the second test case the static eigenvalue (k eff ) is also studied.
3.1. Test case description 3.1.1. One-group pincell The first test problem selected is a one-group pincell benchmark designed to investigate the numerical efficiency of the NURBS basis functions modelling circular pins. Six test cases were studied with increasingly heterogeneous material properties. The prescribed cross-section data for these test cases has been selected so that the challenge for the numerical schemes shifts from the basis functions ability to represent the geometry to their ability to capture sharp features in the solution function. The prescribed cross-section data was generated by varying the relative magnitude of the total macroscropic cross-section in the fuel and moderator regions; whilst maintaining a fixed scattering ratio of c ¼ 0:95. The material properties of the six test cases are detailed in Table 1.
The radius of the pin is 1.0 cm centered inside a square cell of width 4.0 cm as shown in Fig. 6. Reflective boundary conditions are prescribed on all four sides simulating an infinite lattice of pincells. The geometry was discretised with five NURBS patches in an analogous manner to the example shown in Fig. 3 thus each patch has approximately the same area aiming to reduce as much as possible any mesh bias caused by any mesh non-uniformity. A unit source is applied across the full cell. A reference solution was obtained using quartic (p = 4) IGA with 14689 degrees-offreedom, calculated using the IGA compiled in quadruple precision (fp128) and is believed to be accurate up to 14 significant figures. Results for the reference solutions are provided in Table 2.

Seven group UOX pincell
For a more representative problem, the UOX pincell from the OECD/NEA C5G7 LWR benchmark was selected Lewis et al., 2003. The fuel pin has a radius of 0.54 cm centered inside a square cell of width 1.26 cm. Like the single-group pincell, reflective boundary conditions are applied on all four sides. In contrast to the geometry of the one-group pincell the circular fuel pin encompasses a significantly larger fraction of the cell increasing the aspect ratio of the patches representing the moderator. Like in the previous case, a reference solution was obtained using quartic (p = 4) NURBS functions with 32161 spatial degrees-of-freedom. The reference results are provided in Table 3 with the flux normalised to unit fission source. Representative scalar flux contour plots for groups one and seven on a coarse and refined mesh are shown in Fig. 7.

Importance of an exact geometrical representation
The NURBS representation of the geometry is exact at any level of refinement of the IGA model. In contrast, in FEA complex curved shapes are most commonly represented by a polygonal (2-D) and polyhedral (3-D) approximation introducing an underlying error in the numerical discretisation scheme. Higher-order, isoparamet-ric finite elements can also be employed, however these will also generally only approximate complex geometries. Additionally, FE mesh generators will typically not preserve the fissile mass of cir-cular fuel pins, introducing a significant numerical error. This error can be avoided by artificially increasing the radius of the approximate polygon to ensure mass preservation for straight-sided finite element meshes. To mass-preserve, arbitrary geometries with high-order finite elements the meshed volume could be computed by integrating over the resulting mesh and correcting the meshing strategy to compensate. This corrective approach however is not studied in this work and high-order isoparametric meshes are non-mass preserved. To study the benefits of the exact representation of the geometry a series of finite element meshes were generated with as close to identical structure as to those obtained through knot-insertion of the NURBS patches. The similarity in mesh structure ensures that any bias due to preferential element placement is limited. All finite element meshes used in this study were mass preserved. An example of a NURBS mesh and finite element mesh is shown in Fig. 8. Linear, quadratic, cubic and quartic finite element meshes are compared to quadratic, cubic and quartic NURBS meshes.
To ensure that the comparison is focused on the errors introduced by the approximation of the geometry, the continuity of the in-patch NURBS basis functions was restricted to C 0 by repeating internal knot point values. This restriction may limit the potential of the IGA method as this study is focused on the benefits of exact representation of the geometry alone. This problem configuration means that for any given polynomial order the number of elements, degrees-of-freedom and matrix entries are identical for both FEA and IGA schemes. The total number of matrix nonzeroes has been chosen as the metric for comparison as it is considered to be proportional to the cost of a single matrix-vector multiplication. This allows for a fairer performance analysis between schemes of different polynomial order than comparing against the total number of degrees-of-freedom.
Convergence of the cell disadvantage factors for the six singlegroup pincell test cases are shown in Fig. 9. It can be seen that when the ratio of R fuel T to R mod T is low a significant improvement in accuracy is obtained by the IGA schemes compared to the FEA scheme of the same polynomial order. As the cross-section ratio is increased, and the under-resolution of the fuel-moderator interface gradient becomes a dominant source of error then the relative performance advantage of the IGA schemes is reduced as is visible in the case where R fuel T = 100 shown in Fig. 9. When R fuel T is increased to 50 there is no significant advantage by using quadratic NURBS compared to quadratic finite elements. However, if the basis functions are order-elevated then a relative reduction in error is again observed for IGA indicating that the high order finite element schemes are limited by approximation of the geometry rather than the basis function's ability to capture the solution.
In all studied cases, the quadratic, cubic and quartic, straightsided finite element scheme's convergence is restricted to the same rate. As the representation of the geometry does not improve under p-refinement a saturation in the order of convergence is observed. Higher-order IGA schemes produced higher rates of convergence. While the convergence rate of these finite element Table 1 Pincell cross-sections.   schemes is limited by the approximation in the geometry, the exact representation by the NURBS functions yields the full potential of the IGA method. It was observed that the linear finite element scheme yielded the lowest order of convergence illustrating that the ability of the basis functions to represent the solution dominates the error compared to any approximation of the geometry. Convergence of the eigenvalue and the group 1, 4 and 7 cell disadvantage factors for the seven-group UOX pincell are shown in Fig. 10 for the comparison of IGA and mass-preserved, polygonal finite element meshes. Like with the single-group test cases, the NURBS schemes yield a lower error than the FEA scheme of the same polynomial order for all quantities of interest studied. When the variation of the flux profile across the cell is low (indicated by a disadvantage factor close to unity) the improvement in accuracy of the NURBS schemes compared to FEA is lower than when the variation is higher, as observed in group 7. For groups 1 and 4 quadratic IGA and quadratic NURBS produce solutions with similar errors indicating that for these groups the benefits of the exact geometry are low and the ability of the basis to represent the solution is the largest source of error. When the basis functions are orderelevated to cubic and quartic a significant improvement in accu-racy is observed for the IGA schemes confirming the high order finite elements schemes are limited by coarse approximation of the geometry. For high-order finite element schemes, computational work is wasted obtaining an accurate solution to the wrong physical problem. When the quadratic FEA scheme is order elevated a small increase in the rate of convergence for these groups is observed although the increase in the linear system size does not yield improvements in the scheme's efficiency. In contrast for the eigenvalue and the group 7 disadvantage factor a significant improvement in accuracy is observed for the quadratic IGA scheme. For an error of 1 pcm in the eigenvalue a linear system with half the number of matrix non-zeroes (for the same polynomial order both the IGA and FEA schemes have the same number of degrees-of-freedom also) can be used reducing the total memory required and computational time required. Again, it was observed that order-elevation of the NURBS schemes will significantly reduce the error and improve the rate of convergence. However, only when an error significantly below 1 pcm is required are there significant benefits of the higher-order NURBS schemes.
Convergence of the eigenvalue and the group 1, 4 and 7 cell disadvantage factors for the seven-group UOX pincell are shown in Fig. 11 for straight-sided, mass-preserved polygonal finite  elements, non-mass-preserved, high-order, isoparametric finite elements and equivalent NURBS meshes. It can be seen that for all studied quantities the non-mass-preserved, linear finite elements are the least accurate. By preserving the fissile mass the error can be reduced by between one and two orders of magnitude while keeping the same geometric approximation. The group 1 and group 4 cell disadvantage factors show that preserving the fissile mass a more accurate solution can be obtained than the non- Fig. 9. Convergence of the disadvantage factor for the six, single group pincell test cases. Mass-preserved, C 0 -continuous finite elements with a straight-sided polygonal approximation have been compared to C 0 -continuous NURBS isogeometric elements. mass-preserved, isoparametric finite elements for both quadratic and cubic schemes. For both polynomial orders studied the equivalent order isogeometric scheme produces a solution with a lower error compared to both finite element approaches. In contrast, the eigenvalue and group 7 cell disadvantage factor demonstrates that accurate representation of the geometry is more important than mass-preservation for the quadratic and cubic finite element schemes. The improved performance of the quadratic isoparametic finite elements compared to the equivalent isogeometric scheme is believed to be due to a cancellation of errors in this quantity of interest as the same advantage is not observed in the other group fluxes which contribute to the group source term. It was also observed that the isoparametric cubic finite elements do not possess a high-order of convergence for these quantities of interest than the equivalent quadratic scheme and it believed to be due to the rate of the convergence of the geometric approximation being low, dominating the overall solution error.

Importance of increased basis function continuity
NURBS basis functions naturally possess high order (C pÀ1 ) continuity at knot points (analogous to element boundaries in FEA). This is in contrast to the standard finite element techniques where the basis functions are C 0 -continuous at element boundaries. IGA allows for the additional freedom to decrease the basis function continuity by inserting additional knots of the same value, known as k-refinement (Hughes et al., 2005). Each repeated knot value reduces the function continuity by one at the selected point in parametric space. By reducing the basis function continuity additional degrees-of-freedom will be introduced however the average number of basis functions sharing common support will decrease resulting in a sparser linear system to be solved. Fig. 12 shows a representative sample of C 0 -continuous and C 1 -continuous quadratic NURBS basis functions over the UOX pincell geometry. In this study the additional benefit of the increased in-patch basis function continuity of the IGA method has been investigated. Quadratic, cubic and quartic NURBS functions are considered with natural C pÀ1 continuity down to C 0 continuity. Again, the total number of non-zeroes in the linear system has been chosen as the metric for comparison to take into account the increased sparsity of the low-order continuity schemes and to allow for a fairer comparisons of the different polynomial orders. Fig. 10. Convergence of the pincell eigenvalue and groups 1, 4 and 7 disadvantage factor for the UOX pincell. Mass-preserved, C 0 -continuous finite elements with a straightsided polygonal approximation have been compared to C 0 -continuous NURBS isogeometric elements.
Convergence of the cell disadvantage factors for the six singlegroup pincell test cases are shown in Fig. 13. It has been observed that the high in-patch basis function continuity yields a more efficient numerical scheme for the full range of physical problems studied. While all schemes of the same polynomial order reach the same rate of convergence, schemes with high-order continuity reach the asymptotic region more quickly. In addition, it is interesting to note as the basis function continuity is dropped a significant reduction in efficiency is observed after the first set of repeated knots are inserted. The C 1 -continuous and C 0 -continuous cubic schemes are much closer in numerical efficiency compared to the natural C 2 -continuous scheme. This trend is repeated for the quartic schemes where there is little observed difference between the C 0 -continuous, C 1 -continuous and C 2 -continuous schemes.
The convergence of the eigenvalue and group 1, 4 and 7 cell disadvantage factors for the seven-group UOX pincell are shown in Fig. 14. Like for the single-group problem, a significant improvement in accuracy is obtained in all quantity of interest studied when the high-order in-patch basis function continuity is pre-served. If the eigenvalue (k eff ) is required to be known to an error of 1 pcm there is little to choose between any of the studied schemes. As the schemes are further refined and the errors are reduced a significant improvement in efficiency is observed for the schemes with high-order basis functions and high-order inpatch continuity. This trend is also observed in the convergence of the cell disadvantage factors for the three energy groups studied.

OECD/NEA C5G7 benchmark
The OECD/NEA C5G7, seven group quarter-core benchmark is a representative, heterogeneous reactor physics problem routinely used to test core neutronics analysis tools Lewis et al., 2003. The reactor contains two mixed-oxide (MOX) and two uraniumdioxide (UOX) fuel assemblies in a lattice configuration, as shown in Fig. 16. Each assembly contains 289 fuel pins or guide tubes. The heterogeneous fuel pin model represents a significant challenge for finite element codes to obtain accurate solutions efficiently. Group 1, 4 and 7 scalar flux solutions are shown in Fig. 15. Fig. 11. Convergence of the pincell eigenvalue and groups 1, 4 and 7 disadvantage factor for the UOX pincell. Isoparametric C 0 -continuous finite elements (ISO) have been compared to mass-preserved, C 0 -continuous finite elements with a straight-sided polygonal approximation (MP) and also to C 0 -continuous NURBS isogeometric elements.
The NURBS model for the C5G7 benchmark was built using a total of 5849 bi-quadratic patches. For each of the 1156 pincells, five NURBS patches were used requiring a total of 5780 patches for all pincells. For the adjoining reflector, each pincell was connected to a single bi-quadratic patch of size 1.26 cm Â 21.42 cm and a single patch of size 21.42 cm Â 21.42 cm was used to complete the geometry in the lower right corner. Each NURBS patch was assigned two integer values corresponding to the number of knots to be inserted in each parametric direction of the parameter space to perform h-refinement. The coarsest model had zero knots inserted into the fuel assembly patches and 16 knots into the moderator patches in the direction moving away from the assemblies. This resulted in a mesh where the elements in the reflector are square and the same size as an individual pincell. The coarsest mesh was discretised with approximately 25,000 degrees of freedom. It must be noted that this level of mesh refinement in the moderator is heuristic and may not be the most efficient method of solving this problem. Refining the mesh by adding additional elements close to the interface of the reflector and the fuel assemblies may better capture the high flux gradients as the neutrons are thermalised leaving the fuel. However, for this work we wish to study the performance of NURBS based IGA against FEA. An equivalent set of linear and quadratic finite element meshes, with a mass-preserved polygonal approximation were constructed with as close to identical mesh structure as the NURBS model. For the NURBS model we maintain the naturally occurring C 1 continuity of the basis functions as this was previously found to provide increased accuracy.
To study the convergence of the IGA against FEA a high resolution reference solution to the problem was calculated. Convergence plots of the error against the number of matrix non-zeros is given in Fig. 17.
A converged value for k eff of the C5G7 benchmark of 1.183238 (3) for diffusion theory using IGA with a NURBS model was obtained, as shown along with assembly and pinpowers in Table 4. This is consistent with published results from codes including DORT Pautz et al., 2003 andCRONOS F. Moreau et al., 2002. The IGA results were then compared with both linear and quadratic finite element calculations using as close to identically structured meshes as possible to ensure as fair of a comparison as possible. IGA results are detailed in Table 5, with linear and quadratic finite elements results detailed in Tables 6 and 7 respectively. Fig. 17 shows that the eigenvalue of the IGA model converges more rapidly for a given number of matrix non-zeros in the assembled linear system. The difference in the relative accuracy of the two methods becomes greater when the physical mesh is refined.
For an eigenvalue accuracy of 1 pcm (10 À5 ) (compared to the highly refined IGA reference calculation) the required finite element matrix must contain approximately two and half times more matrix entries for C 0 quadratic elements and nearly five and a half times more for C 0 linear elements compared to the IGA calculations. The storage of the preconditioned linear system for each group requires a significant amount of system memory. Reducing the size of the linear system to be stored will allow for more accurate results to be calculated if total system memory is a limiting factor. For high-order even-parity P N transport schemes, for which the diffusion method of this study may form the basis, being able to store the smaller linear systems and associated preconditioners in system memory may result in a significant reduction in overall computing time. It should be noted that the preprocessing, basis function integration and matrix assembly is typically much less than one percent of the total computation time for both IGA and FEA schemes and we do not take this into consideration in our cur- Fig. 12. The continuity of basis functions can be reduced by repeating knot values. In the above five biquadratic patches are used to represent the UOX pincell. In the top row of plots, the central value of each knot vector of all patches is repeated reducing the basis function continuity across knot spans to C 0 -continuous. In the bottom row the natural C 1 continuity of the NURBS basis functions is demonstrated. It can be seen that the basis function is non-zero across multiple elements and posseses a continuous first derivative. rent evaluation of the methods. In the initial preparation for this work the effect of the integration of the NURBS basis functions using Gaussian quadrature was studied. It was found that p þ 1 points per parametric direction provided sufficient accuracy once the mesh was refined. The error from inexact integration of the basis functions was low compared to the mesh's ability to Fig. 13. Convergence of the disadvantage factor for the six, single group pincell test cases. Quadratic, cubic and quartic NURBS basis functions have been studied with all possible in-patch continuity possibilities.
accurately capture the solution and it is believed to be sufficient to integrate the NURBS function with the same number of quadrature points as a polynomial of the same order.
When comparing the calculated pin powers and assembly powers, it has been observed that the IGA simulations performed significantly better than the FEA calculations. For the maximum pin power the isogeometric method used a matrix with five times fewer non-zeros than the quadratic finite elements for an error of 10 À5 . Interestingly, for the maximum pin power the linear finite elements also performed better than the quadratic finite elements. Since straight sided quadratic and linear elements have been used (to allow for preservation of the fissile mass) the boundary of the pin is approximated by a polygon with a different external length compared to the exact circumference of the circular pin. Since each  quadratic element has three nodes per straight side compared to two on a linear element the quadratic mesh has fewer straight sides for the same number of nodes, creating a poorer approximation. Since the NURBS mesh provides an exact geometrical description at the coarsest level (nine degrees-of-freedom for the circle) both the mass and boundary of the pin is modelled exactly allowing for improved accuracy. The minimum pin power and assembly powers show similar trends and a significant improvement in accuracy for the IGA method compared to FEA.
It must be noted that IGA does not appear to show significant improvements over the FEA for coarse meshes in this specific problem. It is believed that this is due to having too few degrees-offreedom to accurately represent the solution so the advantage of the exact geometry is not yet realised. At coarse levels of refinement, the advantage of the exact geometry in the IGA scheme is not yet realised. However, if more accurate solutions are required then IGA can provide significant improvements in accuracy, unlike the linear and quadratic finite elements which vary in relative accuracy for specific physical quantities. The increased complexity of solutions for second-order forms of the transport equation such as the even-parity equations should benefit from the increased continuity of a NURBS model compared to a C 0 finite element mesh potentially allowing for smaller spatial linear systems enabling a higher-order angular approximation for the same amount of computational work and required system memory.

Conclusions
A study of IGA applied to the neutron diffusion equation has been carried out. A comparative study of FEA and IGA applied to a series of single group pincell test cases and a seven group UOX pincell using the diffusion approximation has been performed. Our analysis decomposed the potential advantages of IGA (the increased continuity and exact geometry) to attempt to determine the benefits of each. Our results showed that the C 0 NURBS isogeometric analysis can provide improvements in accuracy compared to both mass preserved, polygonal finite elements and non-mass preserved high-order finite element meshes of the same polynomial order. Compared to these finite element schemes, the NURBS patches exactly represent the geometry allowing for better approximations in the numerical solution for the same number of ele-ments and degrees-of-freedom. In addition, while the rate of convergence of the higher-order finite element schemes saturated at quadratic (p = 2) significant improvements in accuracy can be obtained by order elevating the NURBS meshes. While masspreservation is critical for linear finite elements, the importance was found to vary for higher-order finite element schemes. In the study of the UOX pincell it was found that for some quantities of interest mass-preserved, polygonal, straight-sided finite element meshes produced more accurate solutions than non-masspreserved isoparametric meshes of the same polynomial order. Depending on the quantity of interest the solution accuracy may be dominated either by the error in the fissile mass, the error in the geometric representation or the basis functions ability to represent the solution. As the geometry and fissile mass are both correct for the NURBS meshes, the solution error is determined solely on the functions ability to represent the solution function and reliably out perform both finite element schemes. Our results from the C5G7 benchmark show that an accurate representation of the geometry is beneficial to measuring the maximum pin power and the quadratic, mass preserved, polygonal finite element scheme perform particularly poorly, even compared to the equivalent linear finite elements. The NURBS based description can possess greater than two orders of magnitude improvement in accuracy for a similarly sized linear system.
The increased continuity of the NURBS basis functions (C pÀ1 ) across knot lines was also found to improve the accuracy of the calculations, especially when the mesh was refined. With coarse meshes containing multiple patches the benefit of the increased continuity was reduced since the basis functions are C 0 across patch interfaces. If a geometry is strongly heterogeneous, requiring many smaller patches then this advantage may be reduced.
A key challenge with the NURBS based IGA studied in this work is the propagation of refinement due to the tensor product structure of the patches and the difficulty producing locally refined meshes. A method to circumvent this limitation may be to refine each patch independently based on a given error estimator and then constrain the basis functions on adjacent patches to conform Hughes et al., 2005. The fundamental idea of unifying the geometrical description analysis representation provides significant opportunities for implementing adaptive schemes with significantly reduced complexity compared to FEA.