(2017) hp-Adaptive discontinuous Galerkin methods for neutron transport criticality problems. SIAM Journal on Scientific Computing, 39 (5). B916-B942. ISSN 1095-7197

. In this article we consider the application of high-order/ hp -version adaptive dis- continuous Galerkin ﬁnite element methods (DGFEMs) for the discretization of the k eﬀ -eigenvalue problem associated with the neutron transport equation. To this end, we exploit the dual weighted residual approach to derive a reliable and eﬃcient a posteriori error estimate for the computed critical value of k eﬀ . Moreover, by exploiting the underlying block structure of the hp -version DGFEM, we propose and implement an eﬃcient numerical solver based on exploiting Tarjan’s strongly connected components algorithm to compute the inverse of the underlying transport operator; this is then utilized as an eﬃcient preconditioner for the k eﬀ -eigenvalue problem. Finally, on the basis of the derived a posteriori error estimator we propose an hp -adaptive reﬁnement algorithm which automatically reﬁnes both the angular and spatial domains. The performance of this adaptive strategy is demonstrated on a series of multienergetic industrial benchmark problems. In particular, we highlight the computational advantages of employing hp -reﬁnement for neutron transport criticality problems in comparison with standard low-order h -reﬁnement techniques.

1. Introduction. The theory of neutron transport is an active area of research that has continued to be developed over the last seventy years. This growth has been stimulated by not only the vast increase in computer power, but also the everexpanding range of applications. Indeed, nowadays, neutron transport theory is exploited in applications ranging from the production of radioisotopes for medical use to neutron well-logging in the oil industry to power generation in nuclear reactors. It is this latter application that is the key focus of this article. A modern nuclear reactor is an extremely advanced piece of engineering that must operate both reliably and safely for a long time (a typical operating life is several decades). One of the key mathematical models is represented by the neutron transport equation (also referred to as the Boltzmann equation or Boltzmann transport equation), which models the distribution and energy of neutrons inside a nuclear reactor. The angular neutron flux, can be thought of as the randomly driven "dance" of neutrons in a system that is undergoing fission. The basic process is as follows: neutrons move through a region at a variety of different trajectories and energies/speeds. Within the region there are a variety of different nuclides (in a nuclear reactor, one can expect to find various isotopes of uranium, boron, gadolinium, and a variety of other materials). Collisions between neutrons and atomic nuclei result in either scattering of the neutron, capture of the neutron, or fission of the nucleus (a process which releases more neutrons and creates different nuclides). Each of these events take place with varying probabilities depending on the energy of the incident neutron and the type of nucleus. Additionally, The outline of this article is as follows. In section 2 we first introduce the neutron transport equation, together with the corresponding k eff -eigenvalue problem. The hpversion DGFEM discretization is then introduced in section 3. Section 4 is devoted to the construction of an efficient solver based on exploiting the block structure of the underlying DGFEM. In section 5 we derive a goal-oriented a posteriori error estimator for the reliable computation of k eff . Moreover, we consider the design of a suitable hprefinement algorithm which automatically refines the computational meshes defined over both the angular and spatial domains. A series of multienergetic industrial benchmark problems are considered in section 6. Here, the quality of the computed a posteriori error estimate is studied; moreover, comparisons between the exploitation of h-and hp-refinement strategies are reported. Finally, in section 7 we summarize the work presented in this article and draw some conclusions.
2. Neutron transport problem. Given that Ω is a bounded, connected Lipschitz domain in R d , d ≥ 1, with boundary ∂Ω, the neutron transport equation is given by: find ψ = ψ(x, µ, E, t) such that 1 v ∂ψ ∂t for (x, µ, E, t) ∈ Ω × S 2 × R + × R + , subject to the initial/boundary conditions where the inflow and outflow portions of the spatial domain Ω are defined, respectively, for each µ ∈ S 2 , by ∂ µ − Ω = {x ∈ ∂Ω | µ · n < 0} , ∂ µ + Ω = {x ∈ ∂Ω | µ · n ≥ 0} . Here, n denotes the unit outward normal vector on the boundary ∂Ω, S 2 is the surface of the unit sphere in R 3 centered at the origin, v is the speed of the neutrons, given by v = v(E) = 2E/neutron mass, χ = χ(E) is the probability of a neutron being generated from fission, ν = ν(E) is the average number of neutrons generated from fission, and Σ t = Σ t (x, E, t), Σ s = Σ s (x, µ , µ, E , E, t), and Σ f = Σ f (x, E , t) denote the total, scattered, and fission cross sections, respectively.
As noted above, we focus on the issue of determining the criticality of a given system based on considering an associated eigenvalue problem. We point out that there are two versions of the criticality problem which both take the form of an eigenvalue problem: these are the time decay or α-eigenvalue problem and the reactivity or k eff -eigenvalue problem. To obtain the α-eigenvalue problem we consider solutions of (2.1) of the form ψ(x, µ, E, t) = ψ α (x, µ, E)e αt . Substituting this into (2.1), assuming the material coefficients are independent of t, yields the following eigenvalue problem: find (α, ψ α ) such that in Ω × S 2 × R + , (2.4) ψ α = 0 on ∂ µ − Ω × S 2 × R + , (2.5) C(ψ α ) = 1, (2.6) Downloaded 10/13/17 to 128.243. 38.220. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php B919 where C(·) is a suitable normalization condition; for example, we may select C(·) = · L2(Ω×S2×R+) . Here, and throughout the rest of this article, we have assumed vacuum boundary conditions, i.e., ψ D ≡ 0. The leading eigenvalue with largest real part, denoted by α c , determines the criticality of the system. Assuming α c is real, for a subcritical problem we have α c < 0, α c > 0 for a supercritical problem, and α c = 0 for a critical problem.
Alternatively, criticality can be determined by evaluating the effective multiplication factor, reactivity or k eff -eigenvalue problem given by: find (k eff , ψ k ) such that The existence in general of the eigenvalue k eff and its associated eigenfunction ψ k are assumed for physical reasons; indeed, it can be shown (cf. [35], for example), that under certain conditions, this eigenvalue is real, positive, and simple. If there are fissions taking place in the system then we assume that by varying the number of neutrons released in each fission it is possible to force the system into a state of balance between the number of neutrons gained and the number of neutrons lost. A more detailed discussion of this eigenproblem is given in [7] including the existence of a discrete dominant k eff -eigenvalue under the multigroup discretization of the energy variable (see section 2.1).
In this case criticality is inferred by the leading (smallest in magnitude) k eff , denoted by k c eff . Given that k c eff is real, if k c eff > 1 then fewer neutrons per fission need to be released to create a critical system, therefore the system is supercritical; if k eff < 1, then more neutrons need to be released per fission, and so the system is subcritical; if k eff = 1 then the system is critical.
2.1. Multigroup approximation for energy. In order to simplify the neutron transport models presented above (cf. (2.1)-(2.3), (2.4)-(2.6), and (2.7)-(2.9)), it is usual to approximate the energy variable by either assuming that all scattering is elastic (cf. [18]), which leads to the so-called monoenergetic approximation, or to consider a discrete number of specific energy groups of physical interest. This latter approach, which leads to the multigroup approximation, will be considered for the k eff -eigenvalue problem stated in (2.7)-(2.9); of course, such an approximation may naturally be defined for both (2.1)-(2.3) and (2.4)-(2.6) as well.
To this end, the full energy spectrum is first restricted to the range of interest . . . , G, each of which represents a given energy group; by convention g = 1 is the highest energy group. The multigroup approximation assumes that for each energy group, the angular neutron flux may be separated into the product of its energy dependent and independent parts, i.e., for energy group g, g = 1, . . . , G, we write where f = f (E) is a piecewise smooth function on each interval [E g−1 , E g ], g = Downloaded 10/13/17 to 128.243.38.220. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php 1, . . . , G, which is normalized so that The multigroup approximation is then determined by replacing the integrals with respect to E with the sum of the integrals over each energy subinterval [E g−1 , E g ], g = 1, . . . , G, and integrating the resulting integro-differential equation with respect to E over the interval [E g−1 , E g ]. Thereby, by employing the normalization of f (E), together with the group parameters the multigroup approximation to (2.7)-(2.9) is given by: where ψ = (ψ 1 , . . . , ψ g ) and C M (·) is (again) a suitable normalization condition; for example, we may select Here, for simplicity of notation, we again use k eff to denote the effective multiplication factor defined in (2.12)-(2.14) though we stress that in general this will be different than k eff satisfying (2.7)-(2.9). In the case when G = 1, so that the interval [E 0 , E 1 ] contains all values of interest in the energy spectrum, we deduce the monoenergetic approximation of the neutron transport equation.
Remark 2.1. For the purposes of this article, we assume that the scattering cross section is isotropic; thereby, Σ s,g →g (x, µ , µ) ≡ Σ s,g →g (x). We note that this approximation is valid for scattering off nuclei with large atomic mass.
3. DGFEM discretization. In this section we consider the DGFEM discretization of the k eff -eigenvalue problem (2.12)-(2.14); in particular, we employ a discontinuous piecewise polynomial approximation in both the spatial and angular domains Ω and S 2 , respectively. To this end, we write T S = {κ S } to denote a shape reg-Downloaded 10/13/17 to 128.243.38.220. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php B921 ular subdivision of the spatial domain Ω into disjoint open elements κ S such that Ω = κ S ∈T S κ S , where each κ S ∈ T S is a smooth bijective image of a fixed reference elementκ, i.e., κ S = F κ S (κ), andκ is either the open unit d-simplex or the open unit hypercube in R d , d ≥ 1. We allow the mesh T S to be 1-irregular, i.e., each edge of an element κ S ∈ T S may contain (at most) one hanging node. To each element κ S ∈ T S we associate a polynomial degree p κ S , p κ S ≥ 0, and collect the p κ S in the polynomial degree vector p = (p κ S : κ S ∈ T S ). With this notation we define the spatial hp-finite element space by where M p (κ), p ≥ 0, is either the space P p (κ) of polynomials of total degree at most p onκ, ifκ is the reference d-simplex, or the space Q p (κ) of all tensor-product polynomials onκ of degree at most p in each coordinate direction, ifκ is the unit reference hypercube in R d .
Next, we define a finite element space over the angular domain S 2 . First, however, we note that the angular domain may be parameterized by spherical polar coordinates, i.e., given µ ∈ S 2 , we have that where (ϕ, θ) ∈ (0, π) × [0, 2π). Here, the polar coordinate ϕ is measured from the north pole, and the azimuthal coordinate is θ. Thereby, we write T A = {κ A } to denote a shape regular subdivision of the parameterized domain (0, π) As above, we define the polynomial degree vector q = (q κ A : κ A ∈ T A ), q κ A ≥ 0, and the corresponding hp-finite element space by With this notation, over the full space-angle mesh we define the hp-finite element space Given κ S ∈ T S , the trace of a smooth function v on ∂κ S (the boundary of κ S ), relative to κ S , is denoted by v + . Then for almost every x ∈ ∂κ S \∂Ω, there exists a unique κ S ∈ T S such that x ∈ ∂κ S ; thereby, the outer/exterior trace v − of v on ∂κ S \∂Ω, relative to κ, is defined as the inner trace v + relative to the element(s) κ S such that the intersection of ∂κ S with ∂κ S \∂Ω has positive (d − 1)-dimensional measure. With this notation, the DGFEM approximation for the neutron transport k eff -eigenvalue problem is given by: find k eff,h ∈ R and ψ h,g ∈ V p,q (T ), g = 1, . . . , G, such that for all v h ∈ V p,q (T ), g = 1, . . . , G, where ψ h = (ψ h,1 , . . . , ψ h,G ). Here, T (·, ·), S(·, ·), and F (·, ·) represent the DGFEM discretization of the transport, scattering, Downloaded 10/13/17 to 128.243.38.220. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php and fission terms, respectively. Indeed, here we have that where w = (w 1 , . . . , w G ). Here, H µ (w + , w − , n κ S )| ∂κ S ×κ A , which depends on both the inner and outer trace of w on ∂κ S × κ A , κ = κ S × κ A ∈ T , and the unit outward normal vector n κ S to ∂κ S , is a numerical flux function. This may be chosen to be any two-point monotone Lipschitz function which is both consistent and conservative; see [28,41], for example. In the current setting, the most natural choice of numerical flux is the standard upwind flux given by here, in accordance with the boundary condition (2.13), we set w − = 0 when ∂κ S ∩ ∂Ω = ∅.
Remark 3.1. In the case when piecewise constant polynomials are employed over the angular mesh T A , i.e., when q ≡ 0, then the above DGFEM approximation of the neutron transport k eff -eigenvalue problem is essentially equivalent to employing a discrete ordinates approximation over the angular domain S 2 ; cf. [17], for example. Indeed, let us, for simplicity of presentation, consider the monoenergetic case corresponding to G = 1, then writing ψ h ≡ ψ h,1 , we wish to determine k eff,h ∈ R and ψ h ∈ V p,q (T ), such that Writing µ i , i = 1, . . . , N O , to denote the centroids of the elements in the angular mesh T A , where N O is the number of elements in T A , then writing ψ 3) may be written in the following form: find k eff,h ∈ R and ψ h , and w i denotes the measure of the ith angular element. For simplicity of notation, in this monoenergetic setting, we have written Σ t , Σ s , χ, and ν instead of Σ t,1 , Σ s,1→1 , χ 1 , and ν 1 , respectively. The equivalence of the numerical flux in (3.1), and the variant in (3.5), which is sampled at the centroids of the angular elements, requires the sign of µ · n κ S to be consistent between two neighboring spatial elements; otherwise, in general these terms will not be identical. In practice, however, we observe that the DGFEM (3.3) and its discrete ordinates variant (3.5) yield almost identical results; cf. [37].
4. Ordered-based iterative solver. In this section we develop an efficient solver for the computation of the DGFEM approximation of the neutron transport k eff -eigenvalue problem given in (3.1), (3.2). To this end, for simplicity, we consider the monoenergetic case corresponding to G = 1; cf. (3.3), (3.4) in Remark 3.1; the extension of the proposed solver to the multigroup approximation, i.e., when G > 1, follows in an analogous manner.
denotes the ith entry of the N -vector Ψ. Thereby, we may rewrite (3.3) in the following matrix form: find k eff,h ∈ R and Ψ ∈ R N such that A common approach exploited within the literature is based on employing a simple power iteration; cf. [1,4,43], for example. Here, the essential idea is to rewrite (4.1) in the following manner: find k eff,h ∈ R and Ψ ∈ R N such that Then, expanding (I − T −1 S) −1 in terms of the Neumann series we get: find k eff,h ∈ R and Ψ ∈ R N such that Thereby, the power method may be directly applied to (4.3) based on truncating the Neumann expansion in (4.4), say, with only K > 1 terms. However, if the dominance ratio, i.e., the ratio of the two successive largest eigenvalues, is close to one, then Downloaded 10/13/17 to 128.243.38.220. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php convergence of the power method applied to (4.3) may be extremely slow; cf. [43]. For a review of the application of this technique applied to neutron transport applications, we refer to [1], and the references cited therein; cf., also, [2,33].
As an alternative approach, we exploit the implicitly restarted Arnoldi method to the k eff -eigenvalue problem, written in the following form: find k eff,h ∈ R and Ψ ∈ R N such that In this setting, instead of employing a Neumann series approximation, the action of (T − S) −1 is computed using the GMRES method, based on employing T −1 as the preconditioner.
We point out that both of the above proposed solvers require the action of the inverse of the transport matrix T; indeed, this represents one of the most computationally intensive parts of both algorithms. With this in mind, we now develop an efficient approach for the application of T −1 . One of the major advantages of employing DGFEMs for the numerical approximation of first-order hyperbolic partial differential equations is that, if the underlying characteristics are straight, for example, in the case of a constant advection field, then the elements may be reordered using a topological sorting algorithm in such a fashion that the underlying transport matrix T is block lower diagonal, in which case the numerical solution can be determined using a block forward substitution algorithm. Indeed, in the case of employing a discrete ordinates DGFEM corresponding to q ≡ 0 (cf. Remark 3.1), then this approach may be directly employed to compute the resulting approximation to the k eff -eigenvalue problem. In general, however, when higher-order finite elements are employed in the angular domain, then groups of spatial elements become coupled, in which case a more sophisticated algorithm is needed to identify groups of elements that must then be solved together. With this in mind below we outline Tarjan's strongly connected components algorithm for identifying the required irreducible ordering for the matrix.
Before we proceed, let us first revisit some graph theory. A directed graph G = (V, E), is composed of a set of vertices V , and a set of ordered pairs of vertices called edges, (m, n) = e ∈ E, where m, n ∈ V . A (directed) path is a sequence of edges which connect a sequence of vertices, such that the edges are all directed in the same direction. A cycle is a set C, A graph is said to be acyclic if such a subset does not exist. A connected component S G is a subset of V for which there is a path between any two nodes in S G . S G is said to be a strongly connected component if there does not exist another connected component S G = S G such that S G ⊂ S G . It is clear that any cycle is also a connected component, though not necessarily a strongly connected component, and that any directed graph may be uniquely partitioned into strongly connected components.
With reference to the discrete ordinates DGFEM (3.5), given a characteristic direction µ i , the topological properties of the underlying spatial finite element mesh T S may be illustrated by representing it as a directed graph G. Each element in T S corresponds to a vertex in V and each element boundary which is not parallel to µ i corresponds to an edge in E, with that edge directed in the direction of the characteristic flow across that boundary. Figure 1 illustrates how a given spatial mesh may Downloaded 10/13/17 to 128.243.38.220. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php B925 extbookfmclassopt: 20pt be represented as a graph in this manner. On the basis of the graph ordering, the structure of the resulting transport matrix T, given the single characteristic direction µ i , now becomes clear. Indeed, T will be block lower diagonal if the elements are reordered so that the outflow element comes before the inflow element; on the other hand, if the ordering sets the inflow element before the corresponding outflow element, then the associated matrix will be block upper diagonal. With this in mind, we assert that an edge (m, n) on the graph corresponds to a lower triangular block in the transport matrix if m comes before n in the current ordering and an upper triangular block if n comes before m. Note that in the unlikely event that an element boundary is parallel to the characteristic direction, the aforementioned term will be zero and thus there will be no corresponding edge on the graph; in Figure 1 this occurs between elements 5 and 6. This also highlights that the ordering is sensitive to small perturbations; indeed, if the interior node connecting elements 5 and 6 is moved slightly to the left or right, then the associated graph will change accordingly. However, we stress that this sensitivity is not detrimental to the overall performance of the proposed ordered solver. Finally, we remark that such an ordering may be efficiently computed by Tarjan's algorithm (see [15]); writing |V | and |E| to denote the number of vertices and edges in the underlying directed graph G, respectively, the number of operations required is O(|V | + |E|).
As noted above, in general, cycles are present in the directed graph G associated with a given spatial mesh and angular element. Thereby, in this setting, we exploit Tarjan's strongly connected components algorithm, which both partitions the graph into strongly connected components and outputs these strongly connected components in reverse topological order. With the output from Tarjan's strongly connected components algorithm it is then possible to implement a version of the above sweeping procedure to efficiently apply the inverse of T by performing linear solves on each strongly connected component in the reverse order that they are returned from the algorithm. The basic idea of Tarjan's strongly connected components algorithm is to use a recursively implemented depth first search to visit every vertex in the graph, putting each onto a stack data structure. The stack contains all vertices that have already been visited but have not yet been assigned to a strongly connected component. The algorithm also requires an array of length |V | which contains indices Downloaded 10/13/17 to 128.243.38.220. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php indicating the order in which each vertex is visited by the depth first search and an array of length |V | which contains a pointer to the vertex with the lowest index that can be reached from each vertex on the stack. Any time a new vertex is visited it receives an index and a pointer and is pushed onto the stack. When the depth first search has finished iterating on a vertex, if its pointer points to itself then it is the root of a strongly connected component and the stack is popped down to it. As the algorithm progresses through the graph, no strongly connected component is output before its successors, and therefore a topological ordering of the strongly connected components is obtained. As before, the number of operations required is O(|V | + |E|). For the purposes of this article, we employ the version from the HSL Mathematical Software Library; see [24]. Now that we have the tools to write the matrix T in block triangular form, we can describe the sweeping procedure for applying T −1 to a vector without assembling the underlying matrix system. Algorithm 4.1 describes the procedure, which is equivalent to a block forward substitution algorithm for a block triangular matrix, in terms of our spatial finite element mesh and angular element. To compute the solution of the strongly connected component matrix equation T loc x = b we employ the MA41 linear solver from the HSL library, together with the block preprocessor developed in [37].
Algorithm 4.1 (Sweeping Algorithm). This algorithm takes a spatial finite element mesh T S , its directed graph G = (V, E) associated with an angular element κ A , and a topologically ordered list of its strongly connected components S i ⊂ V , i = 1, . . . , N S , where i S i = V and sweeps through the spatial domain applying T −1 | κ A to a vector v. The algorithm proceeds as follows: 1: Allocate local matrix T loc , solution x, and right-hand side b. Include values from v in b.

5:
for κ S ∈ S i do 6: Compute volume integrals for κ S and include in T loc .

13:
Compute face integrals for κ S ∩ κ S and include in T loc . 14: Compute face integrals for κ S ∩ κ S

17:
if κ S is marked SOLVED then 18: Get solution values from κ S , x κ S .

19:
Build a block B κ S ∩κ S comprising the face terms involving elements κ S and κ S .

20:
Subtract B κ S ∩κ S x κ S from the κ S values in b. Solve T loc x = b and store solution values. 30: for κ S ∈ S i do 31: Label κ S SOLVED 32: end for 33: end for 5. Error estimation and adaptivity. In this section we consider the derivation of an a posteriori error estimate for the accurate computation of the critical value of k eff (cf. (2.12)-(2.14)), together with the design of the corresponding hp-adaptive refinement strategy. To this end, we exploit the DWR technique (cf. [5,6]); in particular, we refer to [11,20,21] for the application of this technique to eigenvalue problems. Extension of this work to bifurcation problems has been considered in the series of articles [12,13,14]; see [10] for the generalization to hp-adaptivity.

5.1.
A posteriori error estimation. Introducing the notation , the DGFEM approximation of the neutron transport k eff -eigenvalue problem may be written in the following compact form: . Assuming that the eigenfunction ψ is sufficiently smooth (cf. [36,38]), then the consistency property of the numerical flux implies that N (ψ, v h ) = 0 for all v h ∈ R × V p,q (T ), where ψ = (k eff , ψ). Thereby, the following Galerkin orthogonality property holds: In order to derive a computable a posteriori bound for the critical value of k eff , we proceed as in [11,20,21] by defining the (nonlinear) functional Assuming J(·) is differentiable we define the mean value linearization of J(·) between ψ and ψ h by where J [w](·) denotes the Frèchet derivative evaluated at some w ∈ R × V, where V is a suitably chosen space such that V p,q (T ) ⊂ V. Similarly, we write the mean value Downloaded 10/13/17 to 128.243.38.220. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php linearization of the semilinear functional N (·, ·) as for v fixed, at some w.
With this notation, we deduce the following proposition.
Proposition 5.1. Let ψ and ψ h denote the solutions of (2.12)-(2.14) and (5.1), respectively, then the following error representation formula holds: for all z h ∈ R×V p,q (T ). Here, z denotes the solution of the corresponding dual/adjoint problem: (5.5), recalling the linearization undertaken in (5.3) and (5.4), and exploiting the Galerkin orthogonality property (5.2), we deduce that for all z h ∈ R × V p,q (T ), as required.
Rather than linearizing about a convex combination of ψ and ψ h , we simply consider the linearization about just ψ; in practice, this linearization will be performed on the basis of the current numerical solution. Thereby, by differentiating J(·) and N (·), we deduce the following approximate dual/adjoint problem: find z ∈ R × V such that (·) denotes the Frèchet derivative of ψ → C(ψ) at some w. For simplicity of notation, we use z to denote the solution to both dual/adjoint problems (5.5) and (5.8), though we stress that these are clearly not identical, in general. We note that the approximate dual/adjoint (5.8) may be rewritten in the following (equivalent) form: (T (ψ g , z g ) − S(ψ, z g )) = 1 (5.10) Downloaded 10/13/17 to 128.243.38.220. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php B929 for all v ∈ R × V; here, we have exploited the normalization C M (ψ) = 1; cf. (2.14). Following [11] it is more natural to consider the alternative dual/adjoint eigenvalue problem: findk z eff andz g , g = 1, . . . , G, such that (T (ψ g ,z g ) − S(ψ,z g )) = 1 (5.12) for all v ∈ R × V. Thereby, it follows thatk z eff = k eff and hence (k z eff ,z g ), g = 1, . . . , G, is also a solution of (5.9), (5.10).
Exploiting the discrete dual/adjoint problem (5.13), (5.14), we deduce the following approximate error representation formula for all z h ∈ R × Vp ,q (T ).
Remark 5.2. We note that, when a higher-order scheme is employed for the numerical approximation of the dual/adjoint problem, the right-hand side of the approximate error representation formula (5.15) may be added to the computed value k eff,h to provide an improved estimate of k eff (cf. [19,22]); in this case reliable estimation of the error in this improved value is no longer available.

hp-adaptivity.
On the basis of the approximate error representation formula (5.15), in this section we consider the design of an appropriate hp mesh refinement algorithm for the efficient control of the discretization error in the computed value of k eff . To this end, given the tensor-product construction of the hp-finite element space V p,q (T ) = V p (T S ) ⊗ V q (T A ) (employed for the numerical approximation of the eigenvector associated with each energy group), we consider the following splitting: ≡ E µ + E x (5.16) Downloaded 10/13/17 to 128.243.38.220. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php for all z h ∈ R × Vp ,q (T ). Under the foregoing notation, for a function v = (k v eff,h , v), we write, subject to a slight abuse of notation, where Π µ and Π x denote the orthogonal L 2 -projection operators onto the finite element spaces [V q (T A )] G and [V p (T S )] G , respectively.
Here, loosely speaking, we consider E µ and E x to be the error stemming from the angular and spatial discretization, respectively. With this in mind, at each step of the refinement algorithm, we refine either V q (T A ) or V p (T S ) depending on the relative size of |E µ | and |E x |; for the assessment of alternative refinement strategies, we refer to [37]. For the purposes of the current article, we have restricted ourselves to problems consisting of isotropic scattering; cf. Remark 2.1. Indeed, in this setting, the underlying eigenvectors for each energy group are extremely smooth with respect to the angular variable µ; thereby, we only consider refinement of the polynomial degrees q in the angular finite element space V q (T A ) when |E µ | > |E x |. Moreover, a key aspect of attaining reliable error control in an efficient manner is to limit the number of steps required by the underlying refinement algorithm; this is particularly pertinent for the problem at hand, since the computation of numerical approximations of the primal and dual k eff -effective eigenvalue problems (3.1), (3.2) and (5.13), (5.14), respectively, may be extremely expensive. With this in mind, and on the basis of the work undertaken in [37], we only consider uniform polynomial enrichment of V q (T A ).
In the case when |E x | > |E µ |, then the (spatial) finite element space V p (T S ) is adaptively refined. To this end, we first rewrite E x in the following elementwise fashion where η κ S , κ S ∈ T S , denotes the elementwise contribution to E x . On the basis of the size of the elemental error indicators |η κ S |, elements are marked for refinement using the fixed fraction refinement algorithm with refinement fraction set to 25%. Once an element κ S ∈ T S has been flagged for refinement, a decision must be made whether the local mesh size or the local degree of the approximating polynomial should be adjusted accordingly. The choice to perform either h-or p-refinement is based on estimating the local smoothness of the (unknown) analytical primal and dual eigenvectors. To this end, we employ the hp-adaptive strategy developed in [23] (cf. also [16]), where the local regularity of the analytical solution is estimated from truncated local Legendre expansions of the computed numerical solution.
6. Numerical experiments. In this section we study the practical performance of the proposed a posteriori error estimator developed in section 5 within an automatic hp-adaptive refinement procedure for a series of industrial benchmark problems. For purposes of comparison we also present numerical results employing h-refinement with p = 1 and q = 0, i.e., for the (variant of the) discrete ordinates method (cf. Remark 3.1); in this case, when |E µ | > |E x |, we employ sequences of uniformly refined angular meshes which are designed in such a manner as to ensure an equal weighted partition of S 2 . Throughout this section, we restrict ourselves to the two-dimensional case, i.e., d = 2; this is often referred to as pseudo-three-dimensional. This is obtained by considering the spatial variable in Cartesian coordinates x = (x, y, z) and requiring that the angular neutron flux has no dependence on the z variable, which extends to infinity in either direction. As nuclear reactors are frequently designed with a high degree of symmetry in the z-direction, with geometries comprising long, narrow control rods and fuel rods, this geometry provides a good approximation for many existing nuclear reactors. Finally, for computation of the primal and dual/adjoint Downloaded 10/13/17 to 128.243. 38 Fig. 2. Example 1. Computational domain for the SILENE benchmark problem. The material coefficients for the three regions are given in Table 1. Table 1 Example 1. The material coefficients for the three regions marked in Figure 2 for the SILENE benchmark problem.  6.1. SILENE reactor. In this and the following section, we consider two benchmark problems from a set compiled by Albrecht Kyrieleis at Serco Assurance; cf. [29]. Here, we consider the model of an experimental reactor in France known as the SI-LENE reactor. It comprises a steel tank filled with uranyl nitrate with a steel pipe positioned vertically down the middle through which a control rod may be inserted; we refer to Figure 2 for the precise dimensions. We model the unrodded case with a G = 2 multigroup approximation for energy, the cross sections for which are given in Table 1. For this problem, Serco's MONK Monte Carlo code (cf. [3]), computed a value of 1.1293 for the critical k eff -eigenvalue. On the basis of employing a very fine mesh, we computed an accurate value of k eff ≈ 1.129459.
In Tables 2 and 3 we show the performance of the proposed h-and hp-refinement adaptive strategies, respectively. In each case we show the mesh (iteration) number, the spatial and angular error estimators E x and E µ , respectively (cf. (5.16)), the error k eff − k eff,h in the computed critical k eff -eigenvalue, the corresponding effectivity index I eff = (E x + E µ )/(k eff − k eff,h ), the number of degrees of freedom (Dofs) present in Downloaded 10/13/17 to 128.243.38.220. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php  the primal finite element space V p,q (T ), together with the computed k eff -eigenvalue for both the primal and dual/adjoint k eff -eigenvalue problems. The initial spatial mesh is depicted in Figure 3. From Tables 2 and 3, we observe that the proposed (approximate) error representation formula (cf. (5.15)) provides an extremely accurate estimate of the error in the underlying critical k eff ; indeed, on almost all of the meshes designed by the h-and hp-adaptive algorithms the effectivity index I eff is very close to unity. In Figure 4 we plot the results shown in Tables 2 and 3; in particular, we plot the error |k eff − k eff,h | using both h-and hp-refinement against the fourth root of the number of degrees of freedom employed in the (primal) finite element space V p,q (T ) on a linear-log scale; these axes have been selected in order to indicate exponential convergence of the hp-refinement strategy. Indeed, here we observe that (on average) the error in the computed k eff using hp-refinement is roughly a straight line, Downloaded 10/13/17 to 128.243.38.220. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php  thereby indicating exponential convergence. Moroever, from Figure 4 we clearly observe the superiority of employing adaptive hp-refinement in comparison to standard h-refinement: indeed, on the final mesh, the true error in the computed k eff using hp-refinement is almost two orders of magnitude smaller than |k eff − k eff,h | when hrefinement is employed alone.
Finally, in Figure 5 we show the final spatial meshes generated using both h-and hp-refinement; in the latter case, we also show the final distribution of the polynomial degree vector p. In the h-refinement case, we observe that the spatial mesh has only been refined in the region containing uranyl nitrate; indeed, no refinement has taken place in the region containing the inner pipe or steel shielding. On the other hand, in the hp-setting, only p-refinement of the spatial domain has been undertaken; indeed, the final spatial mesh corresponds to the initial one depicted in Figure 3. In Downloaded 10/13/17 to 128.243.38.220. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php    Figure 6 we plot the space averaged flux for the critical eigenfunction in the high and low energy groups; here, we observe that the underlying solutions are extremely smooth, even in the vicinity between the change of materials. Note also, that the choice of refinement depends on the interplay between the smoothness of the primal and adjoint/dual solutions.
6.2. Water-cooled reactor. The second problem we consider from [29] is a model of a water-cooled reactor. Here, the spatial domain comprises a central region containing the nuclear fuel, surrounded by a steel reflector, which is then surrounded by water; see Figure 7 for the exact dimensions. The model incorporates a G = 2 multigroup approximation for energy; the nuclear cross sections are given in Table  4. When applied to this problem, the MONK Monte Carlo code computed a critical k eff -eigenvalue of 1.0118, while the short characteristics code of Baker [4] computed a value of 1.0129. On the basis of a fine grid calculation, we computed k eff ≈ 1.012132. Downloaded 10/13/17 to 128.243.38.220. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php Table 4 Example 2. Material coefficients for the three regions marked in Figure 7 for the water-cooled reactor benchmark problem.   In Tables 5 and 6 we show the performance of the our adaptive algorithm employing both h-and hp-refinement, respectively. Here, we again observe that the effectivity indices are very close to unity on all of the adaptive meshes employed, indicating that the computed error representation formula provides a reliable estimate for the error in the computed critical k eff . In addition, we see that the initial spatial discretization provides a very accurate spatial approximation in the sense that |E x | is considerably smaller than |E µ |. Consequently, the adaptive refinement algorithms for both h-and hp-refinement initially target enrichment of the angular finite element space V q (T A ). In fact, in the case of the h-refinement algorithm, spatial refinement is never selected; on the other hand, for the hp-refinement algorithm, spatial refinement is not selected until the third adaptive refinement step.
In Figure 8 we plot the error |k eff − k eff,h | using both h-and hp-refinement against the fourth root of the number of Dofs employed in the (primal) finite element V p,q (T ) on a linear-log scale. Here, we again observe that (on average) the error in the com-Downloaded 10/13/17 to 128.243.38.220. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php puted k eff using hp-refinement is roughly a straight line, thereby indicating exponential convergence. For this example, on the final mesh, the true error in the computed k eff using hp-refinement is almost three orders of magnitude smaller than the corresponding quantity when h-refinement is employed, for the same number of Dofs.
In Figures 9(a) and (b), we plot the scalar neutron flux (integral of the angular neutron flux over S 2 ) of the critical eigenfunction for energy groups 1 and 2, respectively. Here, we observe that the scalar neutron flux is very smooth for the high energy group corresponding to g = 1, whereas for the lower energy group (g = 2), the corresponding quantity contains localized structures. Furthermore, we note that the peak value of the scalar neutron flux in the high energy group is five times higher than the peak value in the lower one. Figure 9(c) shows the mesh and polynomial distribution p after the final hp-refinement step. Here, we observe that the hp-adaptive algorithm has chosen p-refinement for every element that was selected for refinement. Moreover, we note that the refinement algorithm has identified the structures in the solution of the low energy group, as well as the extrema present in the high energy group.
6.3. KNK fast reactor benchmark. In this final example, we consider a model of the Kompakte Natriumgekühlte Kernreaktoranlage (Compact Sodium Nuclear Reactor Plant) in Germany; we shall refer to it as the KNK fast reactor benchmark. This benchmark was originally published by Takeda and Ikeda in [39] and has since been computed by several authors, including Kim and Cho [27], Wang [42], and Baker [4]. The spatial domain is considerably more complicated than the previous problems, comprising a tessellation of 169 regular hexagons split into 18 regions with 8 distinct materials; see Figure 10. The energy spectrum is discretized into G = 4 energy groups. Takeda and Ikeda published three sets of data for the KNK fast reactor benchmark, modeling the case with control rods inserted, control rods half inserted, and the unrodded case, respectively. We compute the k eff -eigenvalue problem for the latter case; here, the material cross sections are given in Tables 7 and 8. In [39], the authors computed several approximate values of the critical k eff -eigenvalue Downloaded 10/13/17 to 128.243.38.220. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php B937 (a) (b)  The critical eigenvalue for the KNK fast reactor benchmark was computed by the hrefinement algorithm with a p = 1, q = 0 finite element space, as well as by the hprefinement algorithm. The same initial spatial and angular meshes were chosen for both runs. The initial spatial mesh is given in Figure 5.8a. The initial angular mesh comprised 4 elements, with a single angular element in each principal triangle of the hemisphere. The fixed fraction refine percentage in the spatial problem was chosen as α FF = 25%. Table 5.6 gives the total problem sizes and computed eigenvalues from the primal and dual solves at each iteration, together with the effectivities for both k eff and its reciprocal. Table 5.7 gives the spatial and angular error indicators, η S and η A respectively, together with the number of spatial and angular degrees of freedom in each energy The domain is composed of a tessellation of regular hexagons, each of which has sides of length 7.5 cm. The numbers mark the eight different regions; the nuclear cross sections for these are given in Tables 7 and 8. based on employing a variety of numerical methods; typical values lie in the range between 1.0887 and 1.0951. More recently, Wang [42] computed a value of 1.01082, Baker [4] computed k eff ≈ 1.0105, and Kim and Cho [27] published values ranging between 1.0094 and 1.01055. On the basis of a fine grid computation we computed k eff ≈ 1.01055. Downloaded 10/13/17 to 128.243.38.220. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php Table 7 Example 3. Material coefficients for the first four regions marked in the diagram of the KNK fast reactor benchmark problem from Figure 10.     The performance of our adaptive algorithm employing both h-and hp-refinement is presented in Tables 9 and 10, respectively, based on employing the initial spatial mesh given in Figure 11. As in the previous examples, we again observe that the approximate error representation formula (5.15) provides an accurate estimate of the true error in the underlying critical k eff -eigenvalue, even for this extremely challenging problem. Indeed, the effectivity indices are very close to one on all of the adaptive meshes employed. Moreover, as in the previous example, initial refinement is undertaken in the angular domain in both the h-and hp-setting; however, as soon as |E µ | < |E x |, then refinement of the spatial domain is undertaken. The results presented in Tables 9 and 10 are plotted in Figure 12. Once again, we observe exponential convergence of the error in the computed value of k eff when hp-refinement is employed. Moreover, hp-refinement is again more efficient than h-refinement in the sense that the computed error |k eff − k eff,h | on the final hp-mesh is over an order of magnitude smaller than the corresponding quantity computed with h-refinement, for a fixed number of degrees of freedom. Downloaded 10/13/17 to 128.243.38.220. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php  Finally, in Figure 13 we show the final meshes generated by both the h-and hpadaptive refinement algorithms. Here, in both cases, the meshes have been refined in the center of the KNK fast reactor geometry, where the test zone, control rods, and driver are located. There is also a limited amount of refinement in region 5, which contains reflecting material, with a moderator. We note, once again, that the majority of refinement undertaken by the hp-refinement algorithm is p-refinement; however, there is some h-refinement in the region with the moderator, signifying that the solution is less smooth in this part of the domain.
7. Concluding remarks. In this article we have employed the hp-version of the DGFEM for the numerical approximation of the k eff -eigenvalue problem for neutron transport criticality calculations. This scheme is highly advantageous both from the point of view of ease of hp-adaptivity, as well as providing a unified variational discretization of the underlying integro-differential equation with respect to both the spatial and angular dimensions. Moreover, the underlying structure of the DGFEM discretization naturally lends itself to efficient block elimination solution strategies; Downloaded 10/13/17 to 128.243.38.220. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php B941 indeed, in this article we have exploited Tarjan's strongly connected components algorithm in order to compute the inverse of the underlying discretized transport operator. This can be used as an efficient preconditioner for the numerical evaluation of the k effeigenvalue problem. Moreover, the variational framework naturally lends itself to a posteriori error estimation. With this in mind, we have employed the DWR technique to estimate the error in the computed critical value of k eff . On the basis of this error estimator we have designed and implemented both h-and hp-refinement strategies, which are capable of automatically refining the finite element spaces employed within both the spatial and angular domains. This general strategy has been applied to three multienergetic industrial benchmark problems. In each case, the proposed estimator delivered reliable error estimation, in the sense that the computed effectivity indices were extremely close to unity. Moreover, the efficiency of employing hp-refinement in comparison with standard low-order h-refinement techniques was clearly demonstrated for neutron transport criticality problems. Future work will include the application to three-dimensional problems.