An adaptive, hanging-node, discontinuous isogeometric analysis method for the first-order form of the neutron transport equation with discrete ordinate (SN) angular discretisation

In this paper a discontinuous, hanging-node, isogeometric analysis (IGA) method is developed and applied to the first-order form of the neutron transport equation with a discrete ordinate (SN) angular discretisation in two-dimensional space. The complexities involved in upwinding across curved element boundaries that contain hanging-nodes have been addressed to ensure that the scheme remains conservative. A robust algorithm for cycle-breaking has also been introduced in order to develop a unique sweep ordering of the elements for each discrete ordinates direction. The convergence rate of the scheme has been verified using the method of manufactured solutions (MMS) with a smooth solution. Heuristic error indicators have been used to drive an adaptive mesh refinement (AMR) algorithm to take advantage of the hanging-node discretisation. The effectiveness of this method is demonstrated for three test cases. The first is a homogeneous square in a vacuum with varying mean free path and a prescribed extraneous unit source. The second test case is a radiation shielding problem and the third is a 3 × 3 “supercell” featuring a burnable absorber. In ∗ Corresponding author. E-mail address: a.owens12@imperial.ac.uk (A.R. Owens). http://dx.doi.org/10.1016/j.cma.2017.01.036 0045-7825/ c ⃝ 2017 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/ licenses/by/4.0/). 216 A.R. Owens et al. / Comput. Methods Appl. Mech. Engrg. 318 (2017) 215–241 the final test case, comparisons are made to the discontinuous Galerkin finite element method (DGFEM) using both straight-sided and curved quadratic finite elements. c ⃝ 2017 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons. org/licenses/by/4.0/).

equations of electromagnetics. The work presented here is an extension of [19], in which a discontinuous Galerkin IGA method was applied to a hyperbolic system of equations, specifically the first-order form of the neutron transport equation with a discrete ordinates (S N ) angular discretisation.
A major drawback of the spatial discretisation presented in [19], as well as in NURBS-based IGA in general is that the tensor product structure of the basis functions leads to the propagation of mesh refinement throughout the entire physical domain. This can lead to computational effort being wasted by over-refining parts of the physical domain that are of little importance to the accuracy of the solution. For continuous discretisations, existing techniques for dealing with this refinement propagation include T-splines [20] and LR B-splines [21].
In a discontinuous setting, the simpler hanging-node approach [22] can be taken to stop this refinement propagation. This has the additional advantage of enabling adaptive mesh refinement (AMR) techniques to refine the mesh, as each element in the mesh can be refined independently of its neighbours. This allows mesh refinement to be prioritised in regions of the physical domain that have a large error associated with them, as calculated by a selection of heuristic error indicators.
The remainder of this paper is organised as follows. Section 2 explains the underlying concepts of isogeometric analysis methods and the NURBS basis functions the method uses. Section 3 derives the spatially discretised system of equations from the mono-energetic neutron transport equation and details the upwinding methodology in the presence of hanging-nodes. Section 4 describes the method used to determine a unique sweep ordering through the mesh for each discrete ordinates direction. Section 5 describes the AMR strategy employed as well as the error indicators used. Section 6 presents numerical results of applying the adaptive hanging-node method to a variety of verification test cases.

Isogeometric analysis
Isogeometric analysis was introduced in an attempt to unify the geometric description of physical problems, utilised in CAD programs, with the computational analysis and numerical discretisation of the governing equations [12]. This unification was achieved by using the NURBS representation for both the underlying geometry and the numerical discretisation method. A large proportion of the design and analysis cycle time is typically taken up with generating a gap free NURBS representation of the geometry that can then be meshed using a mesh generator [12]. This is particularly the case if adaptivity of the spatial mesh is required, or if the geometry changes in a time dependent problem. In these cases the CAD geometry may need to be revisited every time the mesh changes if there are curved boundaries. This is due to the fact that conventional finite elements cannot exactly represent conic sections such as circles and ellipses or quadric surfaces such as cylinders, spheres and ellipsoids. These are all commonly found in reactor physics, nuclear criticality safety assessment and radiation shielding applications. In particular in reactor physics problems, the fissile mass of the system must be conserved in order to obtain an accurate solution [10]. Even if the fissile mass is conserved, the geometry is often still poorly approximated, and so the leakage of neutrons through surfaces will contain an error. The impact of this imperfect geometrical representation will be covered in Section 6.4.
As IGA uses the same NURBS representation as the CAD geometry to discretise the weak form of the S N equations, it can exactly represent the aforementioned geometric primitives, as well many others [23]. As a direct result of this, the geometry used by the analysis program is exact from the coarsest level of refinement. This allows an ancillary mesh generator to be completely bypassed as any refinement can be performed directly in the analysis program. In contrast, the finite element method requires a mesh generator to be employed every time a more accurate representation of the geometry is required. Finally, the parameterisation of the mapping from the parametric to the physical space does not change under h or p refinement when using IGA [12]. This simplifies the implementation of a hanging-node formulation with an arbitrary distribution of nodes on a given element face.

B-splines
Three properties completely define a B-spline curve, also known as a patch [23] 3. A knot vector Ξ , which is a sequence of non-decreasing real numbers. The length of a knot vector is defined to be the number of values, including repeats, in Ξ . Open knot vectors will be used in the remainder of this paper, in which the first and last knot value are repeated p + 1 times. In this case, the knot vector will have length n + p + 1, Ξ = (ξ 1 , ξ 2 , . . . , ξ n+ p+1 ).
B-Spline basis functions in parametric space are completely defined by their polynomial order and knot vector. A uniform knot vector has equally spaced knot points, otherwise it is non-uniform [23]. Knot points separate the parametric space into knot spans. As these knot points are points of reduced continuity, the knot spans can be considered to be analogous to elements in traditional FEM. If a knot point is repeated m times, it has multiplicity m and the basis functions are C p−m continuous at this point. In particular, if a knot point has multiplicity m = p + 1 then the basis functions will be discontinuous at this point, and the patch is split into two new patches. Denoting basis function i of polynomial order p by N i, p (ξ ), then the basis functions are given by the following Cox-de Boor recursion formula [24,25] for p = 0: and for p ≥ 1 Due to the recursive nature of their construction, the derivatives of B-splines can be represented in terms of lower order B-splines as follows [23]: Terms of the form 0/0 can arise in Eqs. (2) and (3), in which case the term is defined to be zero. These basis functions form a partition of unity at all values of ξ within the knot vector [23], that is: The basis functions are pointwise non-negative, which means that every term in a mass matrix is also non-negative, which can be useful for mass lumping schemes [26]. The basis functions are also interpolatory at the end points of the knot vector, such that N i, p (ξ 1 ) = δ i,1 and N i, p (ξ n+ p+1 ) = δ i,n .
B-Spline curves are then defined by associating each basis function with a control point in physical space and summing the product of the two: Higher dimensional B-Spline objects are constructed from one dimensional B-splines in a tensor product fashion. In two dimensions, polynomial orders p and q, knot vectors Ξ = (ξ 1 , ξ 2 , . . . , ξ n+ p+1 ) and H = (η 1 , η 2 , . . . , η m+q+1 ) and a control net {B i, j : i = 1, . . . , n and j = 1, . . . , m} are used to construct B-Spline surfaces as follows: where the N i, p are as above and the M j,q are defined in the same way using polynomial order q and knot vector H . In three dimensions, we add another parametric direction with polynomial order r , knot vector Z = (ζ 1 , ζ 2 , . . . , ζ l+r +1 ) and changing to a control lattice {B i, j,k : i = 1, . . . , n, j = 1, . . . , m and k = 1, . . . , l} to get the volume: The partition of unity, pointwise non-negativity and interpolatory properties extend into 2 and 3 dimensions due to the tensor product structure. In particular, the interpolatory nature of the basis functions limits the number of nonzero basis functions on each face of an element. This can improve the memory efficiency of a sweep based solver by reducing the size of the mass matrix that must be stored on each face of each element for upwinding.

NURBS
While B-splines can be used to represent a wide variety of geometries, NURBS are needed to exactly represent conic sections. To construct NURBS from B-splines, a weight is associated with each control point, denoted by ω i , ω i, j and ω i, j,k in one, two or three dimensions respectively. One dimensional NURBS basis functions of order p are then defined by [23]: Nˆi , p (ξ )ωˆi (8) with NURBS curves given by: The rational part of the name NURBS comes from the division by B-splines. Two and three dimensional basis functions are then defined as follows: Nˆi , p (ξ )Mˆj ,q (η)Lk ,r (ζ )ωˆi ,ĵ,k . (10b) NURBS surfaces and solids are then constructed by associating each basis function with a control point, analogously to Eq. (9).
All of the problems considered in this work are two-dimensional, and in practice can be constructed from the union of biquadratic ( p = 2) NURBS surface patches using the open knot vectors Ξ = H = (0, 0, 0, 1, 1, 1). As the length of an open knot vector is n + p + 1, there will be 3 basis functions in each parametric direction, for a total of 9 basis functions per patch. An example of a pincell geometry and the associated control points for two of the patches are given in Fig. 1 and Table 1 respectively. Fig. 1. The 5 biquadratic patches used to create a pincell geometry with a pin of radius 1 cm and a cell pitch of 4 cm. The control points for the circle (red) and one of the moderator patches (blue) have been marked on, with control points used by both patches marked in purple. The origin is taken to be at the centre of the circle, and the physical coordinates and weights of these control points are given in Table 1. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Table 1 The physical coordinates and weights for the control points of two of the patches shown in Fig. 1 The interpolatory nature of the basis functions at the end points of the knot vectors forces the corner control points ((1, 1), (3, 1), (1, 3) and (3,3)) to coincide with the corners of the physical patch. The mid-side control points ((2, 1), (1, 2), (3, 2) and (2, 3)) of the circle are then fixed by the equation for a circle arc using quadratic NURBS [26], as is control point (2, 3) of the moderator patch. The remaining control points all have some freedom available in their (x, y) coordinates and weight. The values given in Table 1 have been chosen based on symmetry considerations and, in the case of control point (2, 2) of the moderator patch, an attempt to give a smooth variation in the Jacobian over the patch.
When these patches are refined, knot points are inserted with multiplicity p + 1 [23], so that the basis functions are discontinuous between each element as in the fully discontinuous Galerkin scheme described in [19].

Discrete ordinates equations & weak form
In the following derivation vector quantities are in bold while matrices are double underlined. The steady state, mono-energetic neutron transport equation with isotropic scattering and fixed (extraneous) source is given by [27]: where r is a point in the physical domain V , Ω is the direction of neutron flight given by a vector on the unit sphere, ψ is the angular neutron flux (cm −2 s −1 Sr −1 MeV −1 ), φ is the scalar neutron flux (cm −2 s −1 MeV −1 ), Σ t is the macroscopic total cross-section (cm −1 ), Σ s is the macroscopic scattering cross-section (cm −1 ) and Q is the isotropic extraneous source strength (cm −3 s −1 MeV −1 ). This work will only consider problems in two-dimensional space V ⊂ R 2 , although the derivation is identical in one or three-dimensional problems. The discrete ordinates angular discretisation then splits the angular dependence into N discrete directions Ω k , along which Eq. (11a) is solved for the ψ k (r) = ψ(r, Ω k ). Each ordinate also has a corresponding quadrature weight w k for performing angular integrals, which gives an approximation to the scalar flux: The steady state, mono-energetic, fixed source discrete ordinates equations for neutron transport with isotropic scattering are then given by: As with traditional DGFEM, Eq. (13) is multiplied by a test function v(r) ∈ H 1 (V e ), as both the gradient and the trace of v on the element boundary will be needed [28], and then integrated over an element V e : The divergence theorem is then applied to the first term of Eq. (14), yielding: where Γ e is the boundary of element V e , and n(r) is the outward pointing unit normal to this boundary. Defining the two inner products (u, v) :=  V e dru(r)v(r) and ⟨u, v⟩ :=  Γ e d Su(s)v(s), combining this notation with Eq. (14), substituting in (15) and utilising the fact that the cross-sections and external source term are constant within an element yields: where Q e is the source strength in element e. To deal with the boundary integral term, it is necessary to split the element boundary Γ e into inflow and outflow sections; the inflow boundary, where Ω k · n < 0 is denoted by Γ − e , while the outflow boundary where Ω k · n > 0 is denoted by Γ + e . By defining two additional inner products ⟨u, v⟩ , the boundary integral term can be rewritten as: Ω k · ⟨ψ k , vn⟩ = Ω k · ⟨ψ k , vn⟩ + + Ω k · ⟨ψ k , vn⟩ − .
As in traditional DGFEM S N [29], the outflow integral over Γ + e is performed using the flux from within the element. The inflow integral over Γ − e is fully upwinded using the flux ψ ad j from the adjacent element(s). This process is more complicated with the hanging-node formulation than in the conforming mesh case outlined in [19], and will be covered in detail below. If Γ − e is on the physical domain boundary, ψ ad j is instead determined by the appropriate boundary condition. Eq. (16) can then be rewritten as: A solution ψ k ∈ H 1 (V e ) is then sought such that Eq. (18) holds ∀v ∈ H 1 (V e ). The isoparametric concept is followed, in which the solution representation uses the basis functions that represent the geometry. A Galerkin projection is then performed within the element, such that: where there are nnp basis functions within the element. The weight functions R A and the trial functions R B are taken from the same space, the NURBS used to describe the geometry. The source vector, along with mass, streaming and boundary matrices are then defined to be: The volumetric integrals are performed using p + 1 Gaussian quadrature points in each parametric direction over a parent element defined to be with Jacobian matrix J corresponding to this mapping, entry {A, B} of the mass matrix is given by: where the W i are the Gaussian quadrature weights and the fractional term represents the determinant of the Jacobian of the mapping from Defining solution vectors ψ k and φ in the obvious way from Eqs. (19b) and (19c) respectively, Eq. (18) can be written in the compact form: where Ω k = (µ k , η k ) in two spatial dimensions and element A of the vector ψ ad j is defined to be Ω k · ⟨ψ ad j , R A n⟩ − . Defining: This system of equations can then be solved for each element in an order that depends on the discrete ordinates direction Ω k . The difference between the conforming mesh case outlined in [19] and the hanging-node formulation presented here lies in the calculation of the incoming upwinded flux ψ ad j . In the conforming mesh case, on the element boundary the adjacent flux ψ ad j is represented in terms of the same basis functions as the flux ψ k within the element to be solved for. This is no longer the case if an element face contains a hanging-node, such as in Fig. 2. In this case, projected mass matrices are needed of the form  where the R A are the basis functions defined over element e that are non-zero along Γ − e , and the R B are the basis functions defined over neighbouring elements that are non-zero along Γ − e . In DGFEM these element faces are typically straight, or approximated as such [29], and so a single projected mass matrix is needed on each of the faces Γ 2 , Γ 3 and Γ 4 in Fig. 2. The NURBS case is complicated by the fact that, in general, the use of an average normal vector on an element face limits the convergence of the spatial scheme to first order [19]. The product R B R A along with a normal vector must be stored at each quadrature point being used for integration on each element face. This represents the only fundamental difference between the method presented here and existing hanging-node DGFEM methods for the discrete ordinates equations [22,30], and incurs a performance penalty of ≈5%-10% in the transport solve, depending on the polynomial order of the basis functions employed.
In addition, the sign of Ω k · n can change over a curved boundary, and so integration must be performed over partial element boundaries. This is performed as in [19], where the integral over a partial element boundary, e.g. Γ + e is treated as an integral over the whole boundary in which the integrand vanishes where Ω k · n < 0. I.e.
For example in the quadratic case depicted in Fig. 2, p = 2 and therefore 3 quadrature points are required on an element face as demonstrated in [19]. This means integrals over Γ ± 2 , Γ ± 3 and Γ ± 4 must be performed using 3 quadrature points each. Similarly, integrals over Γ − 1 must be split up into integrals over Γ 2 , Γ 3 and Γ 4 and summed, leading to a total of 9 quadrature points being used for this integration. It needs to be investigated when performing integrals over Γ + 1 whether we need to use these same 9 quadrature points, or whether 3 are enough. The method of integration over partial element faces given by Eq. (24) dictates that the method with 9 quadrature points must be used. If not, in cases where the sign of Ω k · n changes over an element face, the outgoing integral over Γ + 1 will make one approximation to this partial boundary integration, while the incoming integrals over Γ − 2 , Γ − 3 and Γ − 4 will make a different approximation. The number of neutrons leaving e 1 across Γ 1 will then not in general be equal to the sum of the neutrons entering e 2 , e 3 and e 4 over Γ 2 , Γ 3 and Γ 4 respectively. The scheme would then not be conservative, and so integration must always be performed using the same quadrature points on any given element face. This simplifies the logic of the scheme, in that integration over any inter-element boundary must always be performed over the finest subdivision of that boundary by the elements on either side of it.

Sweep ordering
As with DGFEM, an order must be chosen in which to solve the elements for each discrete ordinates direction to update the solution at each iteration. This is done using the directed graph (digraph) associated with each direction. For each direction, each element in the mesh has a vertex in the digraph. There is an edge from vertex i to vertex j if neutrons can move from element i to element j in the mesh for this direction. If we define Γ i, j to be the boundary between element i and element j and n to be the outward pointing unit normal on Γ i, j with respect to element i, then there is an edge from vertex i to vertex j if Ω k · n > 0 at any point on Γ i, j .
In traditional DGFEM, element faces are usually either straight/planar, or approximated as such [29]. In this case, meshes of two-dimensional geometries produce acyclic digraphs [29], and a perfect order can be found in which to sweep the elements for each discrete ordinate direction. This is not the case when curved element faces are explicitly taken into account [19]. A digraph G is said to contain a cycle if there exists a vertex v such that by traversing the edges of G, v can be returned to, whilst visiting at least one other vertex inbetween. Cycles can cause problems in sweep ordering, as it becomes unclear which element should be solved for first for a given ordinate if neutrons can move from element i to element j and vice versa. This can also reduce the efficiency of the scheme, as the solution for a direction cannot be obtained by solving for each element once at each iteration.
Strongly connected components (SCCs) are defined to be a set of vertices {v 1 , v 2 , . . . , v n } such that there is a path from v i to v j and vice versa for all i, j ∈ {1, 2, . . . , n}, and are directly correlated to the existence of cycles in the mesh. The partitioning of vertices into SCCs is unique, as these sets of vertices cannot overlap. The condensation of a digraph G is formed by collapsing all the vertices in an SCC into a single new vertex, with edges in the condensation between these new vertices if there were any edges between any of the vertices in the original graph G. The condensation digraph is then acyclic, and so a perfect ordering of these is possible [31].
In the conforming mesh case, it was found in [19] that solving each element within an SCC once is more efficient than iteratively converging the solution within an SCC before moving on to the next one. This approach will be followed here. Therefore a sweep ordering of the elements within the SCC must be found. The method employed for this in [19] involves assigning a weight to each edge in the SCC such that the weight of the edge from vertex i to vertex j α i, j is given by: where Γ + i, j is defined to be the portion of Γ i, j where Ω k · n is positive. This is then a measure of how "outgoing" the direction Ω is on this boundary. Where there was an edge both from vertex i to vertex j and vice versa the edge with the lower weight of the two was deleted. In the absence of hanging-nodes this method worked for all examples considered.
However, when there are hanging-nodes present, this method is no longer sufficient, as demonstrated in Fig. 3. In the SCC consisting of vertices 2, 3 and 4, following this method would delete the edge from vertex 2 to vertex 4. The resulting graph would still contain a cycle 2 → 3 → 4 → 2.
An alternative approach used here is to delete the weakest edge in the entire SCC. This is in contrast with the previous method where edges would only be deleted if there was an edge from vertex i to vertex j and vice versa. If the resulting graph is acyclic, an ordering can be found. If not, the next weakest edge is deleted, and so on until the graph is acyclic and an ordering can be found. Note that unlike the previous method, this is guaranteed to provide an ordering, as even in a pathological example eventually there will be no edges left in the SCC and so finding an ordering is trivial.

Hanging-node adaptive mesh refinement algorithm
While it is possible to specify hanging-node meshes manually by specifying the knots to insert in each parametric direction in each patch, the number of possible permutations grows rapidly with the number of patches. In addition, this requires a priori knowledge of which parts of the geometry will need a finer mesh and which parts can be left with a coarse mesh, which is not always available. In order to take advantage of the discretisation, an adaptive scheme is used here to choose how to refine the patches specified in the geometry.
Here, the coarsest mesh was chosen to be the geometry in an AMR scheme, so that each patch in the geometry becomes an element in the mesh. This geometry is created to be conforming, in the sense that each face of each element that is not on a domain boundary exactly matches the face of another element, i.e. there are no hanging-nodes in the coarsest mesh. This condition could be relaxed, at the expense of complicating the procedure of identifying which elements are adjacent to each other in the coarsest mesh, as the control points on the faces of adjacent elements would no longer coincide.
The solution to the discrete ordinates equations is then calculated on this mesh via the procedure outlined in Section 3. Each element e in the mesh is then assigned an error indicator β e , designed to estimate the error in the numerical solution in this element. An element is then flagged for refinement if it meets the following criteria [22]: with the parameter 0 ≤ α ≤ 1. In this study α was chosen to be 0.3 as in [22]. Elements flagged for refinement are bisected in each parametric direction, creating four new elements. This procedure is then repeated until a given exit criteria is reached, such as the number of AMR iterations or the total number of elements in the mesh.
As the initial mesh is conforming and elements are only ever refined by halving them in each parametric direction, this ensures that element faces are always either fully contained by an adjacent face or contain an integer number of element faces. This is deliberate to simplify the integration over element faces for upwinding.
Three different error indicators will be used here to calculate β e .

Scalar Jump Error Indicator
For this error indicator a variation on the jump-based indicator presented in [30] is used, which measures the jump in the scalar flux around an element boundary. In [30] the square of the scalar flux was used, whereas here the absolute value is used, such that: where Γ e is the boundary of the element, φ e is the scalar flux within the element and φ ad j is the scalar flux in adjacent elements.

Angular Jump Error Indicator
This is similar to the jump-based error indicator, only the jump in the angular flux is used, such that: Residual Error Indicator This more complex error indicator is based on that described in Duo et al. [32]. Defining the residual in an element e and direction k to be: The volumetric residual term for element e is then given by: The jump in the incoming current is also taken into account, through the term: The error indicator β e is then given by: where h e is a measure of the size of the element e, taken here to be the square root of the area of the element.

Method of manufactured solutions
In order to check that the hanging-node formulation attains the correct rate of convergence, the Method of Manufactured Solutions was used. For a full description of the method, see [33]. Here, a modified version of the MMS solution presented in [34] is used to ascertain the order of convergence of the hanging-node formulation. The manufactured solution chosen is for the discrete ordinates equations rather than the full transport equation, and can therefore be used to test the convergence properties of only the spatial discretisation.
The spatial domain is given by the unit centimetre square [0, 1] × [0, 1]. The manufactured solution is then given by: where N is the number of directions in the angular quadrature set used. All MMS calculations were performed with an S 4 level symmetric angular quadrature set, and the scalar flux solution is presented in Fig. 4. The solution is infinitely differentiable in space, allowing a spatial discretisation to attain its full rate of convergence. The manufactured solution is also equal to zero on the domain boundaries, allowing a vacuum boundary condition to be prescribed. The anisotropic source term in each direction can then be calculated:  Σ t = 1 cm −1 and Σ s = 0.5 cm −1 as in [34]. Once the numerical scalar flux solution φ num has been calculated, the L 2 error ϵ 2 is then given by: Expression (35) is then evaluated using 16 2 Gaussian quadrature points in each element. ϵ 2 will then be exact (up to machine precision) for numerical solutions with polynomial order p ≤ 15 on B-Spline meshes. While this is not the case for NURBS meshes, p ≤ 4 in the following examples and so the error introduced by this integration is expected to be negligible [35].
In order to generate meshes with hanging-nodes to test the convergence rate of the scheme, an adaptive procedure was used. In this rare case when the exact discrete ordinates solution is known, the L 2 error in the scalar flux in each element can be calculated. The error indicator β e in element e is then given by: The first mesh types generated were orthogonal grids, in which case the scheme is equivalent to DGFEM. This is achieved by beginning with the square geometry as a single patch. When the AMR iteration procedure is applied to this geometry, the resulting element boundaries are always aligned with the coordinate axes. This allows both the volumetric and boundary integrals to be performed exactly using Gaussian quadrature. Example meshes generated using this procedure are given in Fig. 5.
The results of applying this procedure with first, second and third order basis functions are presented in Fig. 6(a). Fourth order basis functions cannot be used to test convergence rates on orthogonal grids as the manufactured solution itself is fourth order. Fourth order basis functions are therefore capable of exactly representing the solution, and so give an error only due to machine precision, which has been verified to be the case. Fig. 6(a) shows that the hangingnode formulation attains approximately p + 1 order convergence with order p basis functions. The convergence rates are not perfectly straight lines as they are for equivalent formulations without hanging-nodes [19]. However this is to be expected, as the element size is no longer a constant within a given mesh, so h is not a perfect metric of element size.
A more stringent test of the convergence rate attained by the hanging-node formulation was performed by applying the above methodology but using an initial geometry represented by a pincell style mesh shown in Fig. 7. The radius of the circle was given by 1 √ 5π so that the area of each of the 5 patches in the original geometry is identical. This is a much tougher test of the spatial discretisation. The basis functions are now rational, so that both the volumetric and boundary integrals performed using Gaussian quadrature are not exact. The element boundaries can now also be curved, so that integration must be performed along partial element boundaries when discrete ordinate directions are tangent to the boundary.
The convergence rates for second, third and fourth order NURBS are shown in Fig. 6(b). Note that fourth order NURBS on these meshes cannot exactly represent the manufactured solution and so convergence tests can be (a) Orthogonal meshes.
(b) Pincell meshes.  performed using them. Fig. 6(b) shows that with rational basis functions and curved element boundaries, the hangingnode formulation still attains order p + 1 convergence with order p basis functions. The L 2 error is limited to ≈10 −13 due to a combination of machine precision and the condition number of the system.

Homogeneous square in a vacuum
While the solution in Section 6.1 is useful for testing theoretical convergence rates, typical transport solutions are not infinitely differentiable. With piecewise constant cross sections, the solution belongs in H 3 2 at best [36]. This restricts the attainable rate of convergence of even high order spatial discretisations to 3 2 . A problem that demonstrates this reduced convergence order is that of a homogeneous square with vacuum boundary conditions prescribed [22,30].
The problem consists of a 10 cm × 10 cm square with a uniform isotropic unit source, a scattering ratio c = Σ s Σ t = 0.9 and vacuum boundary conditions. This was solved with an S 4 level symmetric angular quadrature set for Σ t = 1 cm −1 and Σ t = 100 cm −1 . The scalar flux solutions to this problem are shown in Fig. 8(a) for Σ t = 1 cm −1 and in Fig. 8(b) for Σ t = 100 cm −1 . In both cases, uniform refinement of the spatial domain was compared to AMR driven by each of the 3 error indicators outlined in Section 5.
For both values of Σ t , L 2 errors were calculated using a reference solution. Two different reference solutions were generated for each Σ t value. In both cases one of these used a uniform refinement of the domain with 4 11 = 4,194,304 quartic ( p = 4) elements. The other was generated using the Residual error indicator driven AMR technique with p = 4, with 1,040,173 elements for Σ t = 1 cm −1 and with 524,656 elements in the Σ t = 100 cm −1 case.
(b) Σ t = 100 cm −1 . Fig. 9. The convergence of the L 2 error calculated using both uniform and adaptive quartic ( p = 4) reference solutions for the Residual error indicator with p = 3 for the Σ t = 1 cm −1 case (a) and the Σ t = 100 cm −1 case (b) of the homogeneous square problem. The errors calculated using the uniform references plateau around 10 −5 and 10 −6.7 for the Σ t = 1 cm −1 and Σ t = 100 cm −1 cases respectively. In both cases this corresponds to the L 2 difference between the uniform and adaptive references. Fig. 9(a) shows that for errors larger than ≈10 −4.5 the L 2 error computed using either of the references give indistinguishable results for the Σ t = 1 cm −1 case. The errors computed using the uniform reference then diverge from those computed using the adaptive reference and plateau at 10 −5 . This is due to the error in the uniform reference. While it has more elements than the adaptive reference, they are evenly distributed, and so many are used to overrefine areas of the geometry of little significance. Fig. 10 shows that the adaptive scheme has concentrated elements along the ordinate directions emanating from the corners of the geometry. This is to be expected as these are the lines across which the true S N solution can be discontinuous [22,30].
Similarly, Fig. 9(b) shows indistinguishable errors computed using both references above ≈10 −6 for the Σ t = 100 cm −1 case, with the error compared to the uniform reference plateauing around 10 −6.7 . Fig. 11 demonstrates the adaptive scheme concentrating elements along the boundary of the domain as there is an extremely sharp gradient in the solution here as demonstrated in Fig. 8(b). The magnified portion of Fig. 11 shows that near the boundary the elements are also being concentrated along the ordinate directions emanating from the domain corners. For these reasons, the remaining results use the L 2 error computed using the adaptive reference solutions.  Figs. 12(a), 12(c) and 12(e) show the convergence in the Σ t = 1 cm −1 case of uniform refinement and all 3 AMR schemes for p = 1, 2 and 3 respectively. For all polynomial orders all 3 AMR schemes behave very similarly, and in all cases have a higher rate of convergence on a per element basis than the uniform refinement method. In the case of the most refined cubic cases, this results in efficiency savings of two orders of magnitude, either obtaining the same level of accuracy using 100x fewer elements, or using the same number of elements and reducing the L 2 error by a factor of 100.
Figs. 12(b), 12(d) and 12(f) show the convergence in the Σ t = 100 cm −1 case of uniform refinement and all 3 AMR schemes for p = 1, 2 and 3 respectively. Again a higher rate of convergence on a per element basis is observed using (e) Σ t = 1 cm −1 , p = 3.
(f) Σ t = 100 cm −1 , p = 3. Fig. 12. The L 2 error convergence of uniform refinement along with all 3 AMR schemes for the Σ t = 1 cm −1 case (left) and Σ t = 100 cm −1 case (right) of the homogeneous square problem with p = 1 (top), 2 (middle) and 3 (bottom). In all cases the AMR schemes have a higher rate of convergence on a per element basis than uniform refinement. This difference in rate markedly increases in the Σ t = 1 cm −1 case as the polynomial order of the basis is increased, which is not observed in the Σ t = 100 cm −1 case. all 3 of the AMR schemes. All 3 AMR schemes behave qualitatively similarly, although are not in as close agreement as in the Σ t = 1 cm −1 case. In particular, the Scalar Jump error indicator in the p = 1 case shows a plateau in the error in Fig. 12(b) around 10 −3 before resuming its convergence. Fig. 13 shows the meshes generated by this scheme at the AMR steps involved. The elements on the boundary of the domain are not being refined in the AMR steps when the L 2 error is plateaued. The error then reduces dramatically once these elements are refined at AMR step 12. One possible cause for this behaviour could be that there is no jump in the scalar flux along the face of the elements that form part of the domain boundary, which would artificially suppress the error indicators in boundary elements using this method. Figs. 14(a) and 14(b) show comparisons between spatial discretisations using the same error indicator but different polynomial orders for Σ t = 1 cm −1 and Σ t = 100 cm −1 respectively. As elements of different polynomial orders have different numbers of basis functions, total degrees of freedom are used as the independent variable. For Σ t = 1 cm −1 there is a significant increase in accuracy on a per degree of freedom basis when moving from linear to quadratic basis functions, and a slight improvement again in using cubics. In the Σ t = 100 cm −1 case all 3 orders behave very similarly until the error reaches ≈10 −2.5 . At that point, the rates of convergence differ significantly, with higher convergence rates observed the greater the polynomial order of the spatial discretisation.

Shielding problem
This problem is an adaptation of Example 2 used in [22], and is representative of a shielding problem [30]. A central source region is surrounded by a highly scattering material, which in turn is surrounded by an absorber. In order to test the NURBS and curved boundaries in a more taxing regime than that in Section 6.1, the Cartesian geometry specified in [22] has been modified so that the central source region is a quarter circle, the surrounding scattering media is a quarter annulus, and the absorber fills the rest of the geometry as shown in Fig. 15. The material properties in each region are given in Table 2.
A level symmetric S 4 angular quadrature set was used to solve the problem. A reference solution was generated using the Residual error indicator with 202,532 quartic elements, the scalar flux of which is presented in Fig. 15.
All three AMR schemes and uniform refinement were compared using quadratic elements. Fig. 16 shows the convergence in the L 2 error of all three AMR schemes and uniform refinement with quadratic elements compared to the reference solution. All of the AMR schemes behave very similarly, and converge significantly faster than the (a) Σ t = 1 cm −1 .
(b) Σ t = 100 cm −1 . For the Σ t = 1 cm −1 case the rate of convergence is very similar for all polynomial orders, with an increase in accuracy per degree of freedom observed as the polynomial order is increased. In the Σ t = 100 cm −1 case the rates of convergence on a degree of freedom basis differ significantly once the error is below ≈10 −2.5 , improving as the polynomial order is increased.    16. The convergence of the L 2 error for all three AMR schemes and uniform refinement with quadratic elements for the shielding problem. All 3 AMR schemes behave very similarly, and converge at a higher rate than uniform refinement. uniform refinement method on a per element basis. Fig. 17 shows both uniform and Residual error indicator meshes with similar numbers of elements (all 3 AMR schemes give qualitatively similar meshes). The AMR schemes refine the mesh along the characteristics across which the S 4 solution can be discontinuous, as well as around the sharp gradient in the flux near the source region. The uniform refinement under-resolves these areas, while over-refining areas of the domain where the solution is flat, such as near the vacuum boundaries.

Supercell with burnable absorber
This test problem consists of a 3 × 3 array of cells, adapted from [37] by only considering isotropic scattering. The outer 8 cells are identical fuel pins, while the central cell contains a burnable absorber, all surrounded by a moderator. The geometry is given in Fig. 18 and the material properties are given in Table 3. The problem was solved with an S 32 triangular-Chebyshev angular quadrature set to attempt to mitigate ray effects as much as possible. The strongly absorbing central pin causes a severe depression in the scalar flux, as demonstrated in Fig. 18. A reference solution was generated using the Residual error indicator with 684,064 quartic elements. Uniform refinement here refers to selecting the number of knot insertions so as to keep the average element size in each patch as similar as possible. The absorption rate in the fuel pins in the middle of one side of the geometry was used as the detector response to compare this uniform refinement to each of the three AMR schemes with quadratic elements in all cases. Fig. 19 shows the convergence in the detector response for quadratic elements using uniform refinement and all three AMR schemes. The AMR schemes are more accurate than uniform refinement for this problem by approximately one order of magnitude when using the same number of elements. The Residual and Angular Jump error indicators give reasonably smooth convergence plots, while the Scalar Jump error indicator has more erratic behaviour. Fig. 20 shows the Residual scheme concentrating the mesh refinement around the sharp gradient induced by the strongly absorbing central pin (all 3 AMR schemes give qualitatively similar meshes). Uniform meshes with a similar number of elements under-resolve this region while over-refining areas of the geometry which are less important, such close to the domain boundary.    19. The convergence of the detector response for all 3 AMR schemes and a uniformly refined method used quadratic elements for the supercell problem. The AMR schemes are approximately an order of magnitude more accurate than uniform refinement for a given number of elements. In order to demonstrate the advantage of IGA over FEM in the context of AMR, the same problem was solved using finite elements. This can be performed using the same program, as there is a one-to-one mapping between tensor-product (bilinear, biquadratic etc.) finite elements and two-dimensional B-splines. As the finite element and B-spline basis functions span the same function space in the parent element, the transformation between the two is simply a change of basis. The B-spline control points B i can therefore be written as a linear combination of the finite element nodes N j : with the A i, j calculated by Galerkin projection [38]. The issue with DGFEM AMR is that unless the CAD geometry is revisited at every iteration, the representation of the geometry does not change as the AMR proceeds. Therefore if this representation is not exact, as with circular pincells, this does not improve throughout the AMR process. In order to demonstrate this, the supercell geometry was set up using biquadratic finite elements using two different approaches.
In the first, all elements were straight-sided, and so the circles and annuli in the geometry are approximated by regular polygons. This approach is frequently taken in practice, as fuel regions can then easily be mass preserved using a priori knowledge of the circle radius and the number of sides in the regular polygon, and fissile mass preservation is crucial to obtaining accurate criticality solutions [10].
In the second, curved elements were used. These are formed from the (non mass preserved) straight-sided meshes described above by moving the mid-side nodes of element faces on a circle boundary onto the centre point of the circle arc between the two corner nodes on the same boundary. This produces a quadratic, rather than linear, approximation to the geometry, but it is still not exact. In particular, the fissile mass will always be underestimated using this approach. Conserving region volumes is no longer simple in this case, as regions are bounded by piecewise quadratic curves rather than regular polygons. To do so requires the numerical integration of the area of the elements approximating the circular geometry, followed by moving the nodes on the circumference of the circle away from the centre of the circle by an appropriate scaling factor, which is based on the ratio of this area to the exact area of the circle. Due to this complexity, this approach would not usually be used in practice.
However, in order to facilitate comparison between these two finite element approaches and IGA, all regions will be mass preserved in the following calculations. The supercell geometry was set up with the circle boundaries represented by regular polygons and piecewise quadratic curves with between 4 and 80 sides, some examples of which are shown in Fig. 21. The Residual error indicator driven AMR scheme was then used to refine this initial mesh, without changing the original geometry. Fig. 22 shows the results of these methods, along with the quadratic IGA scheme using the Residual error indicator for comparison. With the straight-sided elements presented in Fig. 22(a), increasing the accuracy in the geometry representation decreases the error in the detector response. However, for a fixed representation of the geometry, the error in the detector response eventually plateaus as the mesh is adaptively refined. This suggests that the converged detector response for a given polygonal representation of the geometry is not the same as the converged detector response for the true circular geometry, and that this is not a viable method of reducing the error when using straightsided finite elements. In contrast, the error in the detector response in the IGA scheme continues to decrease on the finest meshes considered, and is approximately an order of magnitude more accurate than any of the straight-sided finite element results.
The results when using the curved elements, presented in Fig. 22(b), are more complex. The coarsest representation of the geometry, denoted by DGFEM 4 in the legend, plateaus after a few AMR steps, similarly to the straight-sided element cases presented in Fig. 22(a). The AMR procedure was initially set to exit once a solution had been calculated on a mesh with at least 10 5 elements. These results implied that the detector response for the less accurate geometries (DGFEM 8-20) plateaued at a lower level than those of the more accurate geometries (DGFEM 40-80).  However, the convergence behaviour of the curved elements, and to a lesser extent the straight-sided and IGA elements, are "stepped", in that many AMR steps produce little or no improvement in the detector response error, while others produce a step change. This could be caused by a phenomenon similar to that described in Fig. 13, in which the detector response accuracy only improves when certain regions of the mesh are refined, although this is more challenging to verify for this problem due to the relative complexity of the geometry. This hypothesis is, however, supported by the fact that the values of the error in the detector response at which the various schemes plateau coincide at a few discrete values, e.g. 1.2 × 10 −5 .
In order to verify that the DGFEM 80 curved elements could achieve a detector response error as low as the DGFEM 8-20 schemes, the calculation was repeated with more AMR steps, so that the final mesh had ≈4 × 10 5 elements. It can be seen in Fig. 22(b) that the error does indeed decrease again.
A likely cause of this seemingly counterintuitive behaviour could be the initial finite element meshes used for the curved DGFEM schemes. As shown in Fig. 21, increasing the accuracy of the geometrical representation increases the number of elements in the initial mesh, e.g. the coarsest DGFEM 80 mesh has ×28 as many elements as the coarsest DGFEM 8 mesh. These extra elements are primarily in the burnable absorber and fuel regions, which are exactly the regions of the geometry that are refined the least by the AMR procedure as demonstrated in Fig. 20. These additional elements, while required to improve the geometrical representation with this mesh construction method, therefore essentially constitute wasted computational effort in terms of the ability of the mesh to accurately represent the solution until many AMR steps have been performed. Forming the initial curved element meshes in a less structured manner so that fewer elements are initially in the fuel pins and burnable absorber would likely alleviate this behaviour.
With the straight-sided elements, improving the accuracy of the geometrical representation seems to decrease the error at which the convergence plots in Fig. 22(a) plateau. In order to quantify the error in the geometry, the error in the perimeter of each fuel pin was used, as mass preservation results in the correct area for each pin. Fig. 23 shows that as more and more sides are used to represent a circle using a regular polygon, the corresponding error in the detector response decreases. With the more accurate representations of the geometry, the rate of convergence appears to decrease, although it is possible that this is because the lowest errors taken from Fig. 22(a) have not truly converged yet for these geometries. The same analysis has not been performed for the curved elements, as it is not clear from Fig. 22(b) if any but the coarsest geometrical representation has reached a true plateau with the number of AMR steps performed.

Conclusion
In this paper an adaptive, hanging-node, discontinuous isogeometric spatial discretisation has been developed and applied within the field of neutron transport; although the method is generally applicable to any hyperbolic system of Fig. 23. The error in the detector response of the most refined straight-sided DGFEM meshes versus the error in the perimeter of the circles as modelled by these meshes for the supercell problem. As the geometry becomes more accurate, the error in the detector response decreases.
equations. The complexities that hanging-nodes and strongly curved boundaries introduce to the upwinding process have been overcome by performing all element boundary integration over the finest mesh being used on either side of the element face. This ensures that the scheme is conservative. A robust cycle-breaking algorithm has also been introduced to perform the sweep ordering of the elements for each discrete ordinates direction. In order to take advantage of the hanging-node formulation, simple and heuristic error indicators have been used to select elements for refinement at each step in an AMR scheme.
To verify that the hanging-node formulation obtains the correct order of convergence, a two dimensional, infinitely differentiable MMS problem was solved over the unit square. A discretisation of order p attained order p + 1 convergence with hanging-nodes present in the mesh for both Cartesian and pincell-style meshes for this problem.
Two test cases of a homogeneous square in a vacuum were presented to compare the performance of the scheme to uniform refinement when the underlying transport solution is not infinitely differentiable. When the mean free path is large (10% of the domain length), the AMR schemes concentrate the mesh refinement along the discrete ordinate directions emanating from the corners of the domain, as these are the characteristics across which the S N solution can exhibit discontinuities. When the mean free path is small (0.1% of the domain length) the AMR schemes concentrate the mesh refinement along the domain boundary where there is an extremely sharp gradient in the solution. In both cases the L 2 error in the AMR solutions were consistently lower than the error in the solution when the mesh was uniformly refined on a per element basis for a given polynomial order.
A representative shielding benchmark containing a circular source region, an annulus of highly scattering material and an outer strongly absorbing region was studied to test the L 2 convergence of the scheme when the underlying transport solution is not infinitely differentiable and the geometry requires NURBS rather than polynomials to represent it. The AMR schemes again concentrate the elements in the mesh around sharp gradients in the solution and the discrete ordinate characteristic lines, while uniform refinement under-resolves these regions and over-refines regions of the geometry where the solution is slowly varying. Again the AMR schemes are consistently more accurate than uniform refinement for this problem.
The final test case considered was a 3×3 array of fuel pins, with the central pin replaced by a burnable absorber, all surrounded by a moderating material. This produces a sharp depression in the flux inside the burnable absorber, with a sharp gradient in the moderator region in the immediate vicinity of the burnable absorber. This problem was solved with an S 32 angular quadrature set to try to mitigate the ray effect as much as possible. Again the AMR schemes concentrated the mesh refinement around the sharp gradient in the flux near the burnable absorber material, while uniform refinement under-refined this region and over-resolved other areas of the geometry. In this case the quantity of interest used was the absorption rate in one of the fuel pins. This was found to be approximately an order of magnitude more accurate when calculated using the same number of elements with the AMR schemes compared to uniform refinement.
In order to demonstrate the advantage of IGA for AMR over DGFEM, this same supercell problem was also solved using both straight-sided and curved finite elements. The geometry was fixed by an initial finite element mesh, and then the same AMR procedure as used with the IGA scheme was used to spatially converge the discrete ordinates solution with this geometry. With the straight-sided elements, the errors were initially found to decrease, but after a certain number of AMR steps (dependent on the initial geometry) this error plateaued and no further improvement in accuracy was obtained by refining the mesh further. The errors when the same problem was solved with an IGA scheme of the same polynomial order were approximately an order of magnitude more accurate than the straightsided DGFEM results when the same number of elements were present in the mesh. The most refined straight-sided DGFEM solutions on each geometry were found to converge to the reference absorption rate with respect to the error in the perimeter of the circle that the polygonal geometry is approximating. This could imply that even when the fuel region in reactor physics problems is mass preserved, the error in the leakage rate across the material boundaries inhibits the achievable accuracy of a DGFEM scheme.
The curved finite elements were more accurate than their straight-sided counterparts, and sometimes even as accurate as the IGA scheme for the same number of elements in the mesh. However, these elements are not used in practice in reactor physics problems because of the additional preprocessing step required to ensure that fissile mass is conserved. The clear advantage of the IGA scheme in this respect is that the geometry, and therefore the fissile mass, is automatically exact from the coarsest level of refinement without any preprocessing.
This work has focused on one group problems in order to avoid introducing any bias by combining error indicators from each group in a multigroup problem into a single error indicator for an element in order to select those to refine. The exact geometry preservation of IGA naturally lends itself to overcoming this obstacle with the use of group dependent meshes. This has already been investigated with Cartesian geometries using continuous finite elements for the multigroup diffusion equation [39], which relies on a common description of the geometry at the coarsest mesh level for every energy group. As demonstrated in the supercell problem, in the presence of curved geometries this is not a feasible strategy with DGFEM as the geometrical errors introduced by this fixed geometry will truncate the achievable accuracy of the scheme. This can be overcome by using IGA, as the coarsest mesh description is exactly the geometry of the problem, and will be investigated in future work.