Mark A Goffin

Andrew G Buchan

Anca C Belme

Christopher C Pain

Matthew D Eaton

Paul N Smith

Richard P Smedley-Stevenson

Goal-based angular adaptivity applied to the spherical harmonics discretisation of the neutral particle transport equation

Keywords: Neutron transport, Adaptivity, Goal-based adaptivity, A posteriori error analysis, Spherical harmonics

A variable order spherical harmonics scheme has been described and employed for the solution of the neutral particle transport equation. The scheme is specifically described with application within the inner-element sub-grid scale finite element spatial discretisation. The angular resolution is variable across both the spatial and energy dimensions. That is, the order of the spherical harmonic expansion may differ at each node of the mesh for each energy group. The variable order scheme has been used to develop adaptive methods for the angular resolution of the particle transport phase-space. Two types of adaptive method have been developed and applied to examples. The first is regular adaptivity, in which the error in the solution over the entire domain is minimised. The second is goal-based adaptivity, in which the error in a specified functional is minimised. The methods were applied to fixed source and eigenvalue examples. Both methods demonstrate an improved accuracy for a given number of degrees of freedom in the angular discretisation.

Introduction

The Boltzmann transport equation governs the transport of fundamental particles. The equation is solved to obtain the average statistical distribution of a given particle type within a certain model. This distribution is dependent on space, direction of travel and energy. Knowledge of this distribution is required in disciplines such as atmospheric physics, medical imaging, and the primary focus of this work, nuclear engineering. The deterministic solution of the transport equation is possible through a variety of numerical methods; however, the computational cost of such a solution can grow rapidly if high accuracy is desired. This arises primarily due to the directional/angular dependence of the particle distribution. The typical approach to solve the transport equation is to apply a uniform angular resolution over the space and energy dimensions. However, this requires that the whole space-energy domain use the same angular resolution required to resolve the most directional part of a distribution. This is inefficient because in many problems the angular distribution may be isotropic in certain regions whilst highly anisotropic in others. A specific example is that of a thermal nuclear reactor, high energy neutrons within fuel pins have an isotropic distribution due to the nature of fission, however, the angular distribution becomes more directional in the clad and moderator region outside the fuel pin. The low energy neutrons will in general demonstrate the opposite, they will be isotropic in the moderator but show more directional feature near control rods and fuel pins. Therefore, the required angular resolution for a given accuracy is highly dependent upon the model geometry and material properties. This paper presents an adaptive method designed to optimise the angular resolution to obtain the greatest accuracy for the lowest computational expense.

An adaptive method may be defined as an automated process which modifies discretisations in numerical calculations in order to optimise computational usage. There are two main categories of adaptivity, regular and goal-based. Regular adaptivity reduces the error in the solution over the entire domain. Goal-based adaptivity reduces the error in a user-defined functional of the solution. Goal-based adaptivity in general will produce a more efficient and accurate discretisation for a given quantity than regular. The use of adaptive methods in spatial dimensions for finite elements and similar discretisations is well researched due to its utilisation in a large number of fields. Spatial adaptivity has been applied to the Boltzmann transport equation in various forms and shown to produce beneficial results [START_REF] Duo | A posteriori error estimator and AMR for discrete ordinates nodal transport methods[END_REF][START_REF] Mirza | A spatially adaptive grid-refinement approach for the finite element solution of the even-parity Boltzmann transport equation[END_REF][START_REF] Wang | Three-dimensional h-adaptivity for the multigroup neutron diffusion equations[END_REF][START_REF] Wang | Standard and goal-oriented adaptive mesh refinement applied to radiation transport on 2D unstructured triangular meshes[END_REF][START_REF] Zhang | Spatial adaptivity applied to the variational nodal P n equations[END_REF][START_REF] Lathouwers | Spatially adaptive eigenvalue estimation for the S N equations on unstructured triangular meshes[END_REF][START_REF] Lathouwers | Goal-oriented spatial adaptivity for the S N equations on unstructured triangular meshes[END_REF][START_REF] Turcksin | Goal-oriented h-adaptivity for the multigroup SP N equations[END_REF]. Despite the large computational cost related to the angular component of the transport solution there has not been significant study into adaptivity within the angular dimension.

There are two well known and commonly used discretisations for the angular dimension of the particle distribution; the spherical harmonics (P N ) method and the discrete ordinates (S N ) method. There has been recent investigation into using wavelets as angular basis functions but this is less commonly used [START_REF] Buchan | Linear and Quadratic Hexahedral Wavelets on the Sphere for Angular Discretisations of the Boltzmann Transport Equation[END_REF][START_REF] Buchan | Linear and quadratic octahedral wavelets on the sphere for angular discretisations of the Boltzmann transport equation[END_REF][START_REF] Cao | Solution of neutron transport equation using daubechies wavelet expansion in the angular discretization[END_REF][START_REF] Yang | Application of wavelets scaling function expansion method in resonance self-shielding calculation[END_REF][START_REF] Zheng | An improved three-dimensional wavelet-based method for solving the first-order boltzmann transport equation[END_REF][START_REF] Zheng | Application of the wavelet expansion method in spatial-angular discretization of the neutron transport equation[END_REF]. The literature shows that research has been undertaken in optimising each of these methods. The first known attempt to improve the efficiency of a transport calculation through adjusting the angular discretisation was by Ackroyd and Wilson [START_REF] Ackroyd | Composite finite element solutions for neutron transport[END_REF]. In this work they use the even-parity (second-order) P N method in 1D with variable resolution over the spatial domain. The resolution was set by intuitively using a priori knowledge of the material properties and the predicted solution. This is not an adaptive process but it demonstrates the same principles whereby the automated process is replaced by user intuition. This idea was extended by Mohan et al. to first-order transport in 3D [START_REF] Mohan | Variable order spherical harmonic expansion scheme for the radiative transport equation using finite elements[END_REF]. This demonstrated the variable resolution P N method over space within a spherical geometry. The results showed that solutions with comparable accuracy could be obtained for fewer degrees of freedom and smaller run times.

The first implementation of adaptivity within the angular dimension was completed by Watson [START_REF] Watson | The W N adaptive method for numerical solution of particle transport problems[END_REF]. In his work he developed an angular scheme that used a discrete ordinate discretisation in the polar angle and a wavelet expansion in the azimuthal angle. He demonstrated an adaptive algorithm using a thresholding technique for the coefficients of the wavelet basis functions. This method showed promise because it delivered a smaller error for a reduced number of unknowns and computation time. Research in adaptivity with wavelets was expanded by Buchan et al. which used a two-dimensional wavelet expansion over the surface of the sphere for both the polar and azimuthal angles [START_REF] Buchan | Self-adaptive spherical wavelets for angular discretisations of the boltzmann transport equation[END_REF]. This work uses higher order wavelet functions and applies a similar thresholding type of adaptive algorithm. It was found that the number of unknowns could be reduced by an order of magnitude and retain the same accuracy for two-dimensional models.

There are two known implementations of adaptivity using spherical harmonics, the first was completed by Park [START_REF] Park | Coupled space-angle adaptivity and goal-oriented error control for radiation transport calculations[END_REF][START_REF] Park | Coupled space-angle adaptivity for radiation transport calculations[END_REF]. The focus of his work was a coupled space-angle adaptive process which used the P N method within a finite element discretisation. Park developed an implicit angular error estimator for the second order even-parity transport equation and used it to adapt the resolution of the spherical harmonics over the spatial domain. He demonstrated the use of both regular and goal-based adaptivity for the angular discretisation. The results showed significant decrease in the number of unknowns for a given accuracy and a decrease in run time for most problems. The second implementation of adaptive spherical harmonics was by Rupp et al. [START_REF] Rupp | Adaptive variable-order spherical harmonics expansion of the boltzmann transport equation[END_REF] in solving the Boltzmann transport equation with application to semiconductors. A variable order spherical harmonics resolution was applied over the space and energy domain whilst solving the first order transport equation. The error measure in this work uses an analytical property of the spherical harmonics expansion which relates the rate of decay of the expansion coefficients to the smoothness of the function. This work reported a decrease in the computational cost for a given accuracy by just under an order of magnitude.

Adaptivity within the discrete ordinates discretisation was initially investigated by Stone [START_REF] Stone | Adaptive discrete-ordinates algorithms and strategies[END_REF]. This work used adaptive quadrature sets for the discrete ordinates discretisation. The difference between an interpolated value and a value calculated through sweeping was used as the error measure. If the error was larger than a user-set tolerance then additional quadrature points were added to the discretisation. The method was shown to reduce the ray effects phenomenon for far fewer degrees of freedom than a uniform resolution over space. Jarrell built upon this method by employing the same error measure to adapt using a new type of quadrature set derived from linear discontinuous basis functions on the surface of a sphere [START_REF] Jarrell | An adaptive angular discretization method for neutral-particle transport in three-dimensional geometries[END_REF]. He found that his method obtained an excellent order of convergence but it was difficult to predict the accuracy for highly directional distributions.

This paper extends upon the work completed by Park and Rupp using spherical harmonics within an adaptive algorithm. The two main differences in this work are: (i) the first order transport equation is solved using a different spatial discretisation, and (ii) the error estimators used in the adaptive process are different. This article begins by describing the discretisations used for the first order transport equation in section 2. The regular and goalbased error measures are described in sections 3 and 4. An overview of the procedure used for adapting the angular discretisation is outlined in section 5. Section 6 provides a simple example to demonstrate a problem encountered during the work. The adaptive algorithm has been applied to several examples and the results are displayed in section 7. Finally, the conclusions of the research are provided in section 8.

Transport Equation

Multigroup Energy Discretisation

The multigroup approximation of the fixed source transport equation results in the following set of equations:

Ω • ∇ψ g (r, Ω) + Σ t,g (r)ψ g (r, Ω) = (1) Ω G g =1 Σ s,g →g (r, Ω → Ω)ψ g (r, Ω )dΩ + S e,g (r, Ω) for g ∈ {1, 2, ..., G} (2)
where there are G energy groups and the subscript g denotes the given quantities for each energy group. The multigroup equations for an eigenvalue problem are

Ω • ∇ψ g (r, Ω) + Σ t,g (r)ψ g (r, Ω) = Ω G g =1 Σ s,g →g (r, Ω → Ω)ψ g (r, Ω )dΩ + λ χ g 4π G g =1 νΣ f,g (r) Ω ψ g (r, Ω )dΩ for g ∈ {1, 2, ..., G}. (3)
The group angular flux distribution, ψ g (r, Ω), is dependent on space, r, and direction of particle travel, Ω. The group macroscopic cross sections are represented by Σ t,g , Σ s,g →g and Σ f,g for the total, scatter and fission reactions respectively. In the fixed source equation, the source term, S e , is an extraneous particle source. The fission spectrum and average number of neutrons produced per fission are represented by χ g and ν respectively. The eigenvalue of the equation is represented by λ.

The multigroup equations are coupled through the scatter and fission terms. The equation for each energy group can be considered as a mono-energetic equation with an additional source from the scatter and fission terms. The following sections will consider a single group equation for simplicity and thus the subscript g will be neglected but should be assumed implicit.

Sub-Grid Scale Spatial Discretisation

The sub-grid scale discretisation utilises a finite element based multiscale scheme to approximate the solution of the transport equation [START_REF] Buchan | The inner element sub-grid scale finite element method for the Boltzmann transport equation[END_REF]. The multiscale spatial discretisation is based upon finite elements in which two sets of basis functions are used to represent the spatial variation in the flux. The problem domain, V ⊂ R 3 , is decomposed into the partitioning, J (V ), consisting of a set of η disjoint subdomains, V i for i ∈ {1, 2, ..., η}. The spatial dependence of the solution, ψ(r, Ω), is then decomposed into two components:

ψ(r, Ω) = φ(r, Ω) + θ(r, Ω), (4)
where φ and θ are used to represent the coarse and fine scale components respectively. The coarse scale component is approximated in the continuous finite element space for the partitioned domain J (V ) spanned by η N continuous basis functions,

N i for i ∈ {1, 2, ..., η N }, φ(r, Ω) ≈ φ(r, Ω) = η N i=1 N i (r) φi (Ω). (5)
The fine scale function is approximated in the discontinuous finite element space for J (V ) spanned by η Q discontinuous basis functions,

Q i for i ∈ {1, 2, ..., η Q }, θ(r, Ω) ≈ θ(r, Ω) = η Q i=1 Q i (r) θi (Ω). (6)
The sub-grid scale approximation (equation 4) is substituted into the transport equation which is then multiplied by each of the basis functions N i and Q i and integrated over all space. This gives a set of η N + η Q equations,

V N i (Ω • ∇( φ + θ) + Σ t ( φ + θ))dV = V N i S dV for i ∈ {1, 2, .., η N }, (7)
V Q i (Ω • ∇( φ + θ) + Σ t ( φ + θ))dV = V Q i S dV for i ∈ {1, 2, .., η Q }, (8)
where the source term, S, is composed of the scatter, fission and extraneous terms. Let us first consider the N i weighted set of equations. Green's theorem is used to split the streaming term into a volume and a surface integral. Following this, a strong boundary condition is imposed upon the fine scale solution by setting it equal to zero on the domain boundary. This removes the surface integral of the fine scale solution and gives the following set of equations:

V N i Σ t φ dV - V ∇N i • Ω φ dV + V N i Σ t θ dV - V ∇N i • Ω θ dV + Γ N i (Ω • n) φ dΓ = V N i S dV for i ∈ {1, 2, .., η N }. (9)
The domain boundary is represented by Γ and n is the normal vector at each point on the boundary. Now considering the Q i weighted set of equations, Green's theorem is applied to the fine scale term of the streaming operator but the streaming involving the coarse component is untouched. This set of equations is given by

Ve Q i Σ t φ dV + Ve Q i Ω • ∇ φ dV + Ve Q i Σ t θ dV - Ve ∇Q i • Ω θ dV + Γ in e Q i (Ω • n) θ dΓ + Γ out e Q i (Ω • n) θ dΓ = Ve Q i S dV for i ∈ {1, 2, .., η Q }. (10)
The discontinuous finite element basis functions are non-zero only in a given element of the domain. This allows the integrals over the entire domain to be reduced to the single element in which the given function is non-zero. The subscript e is the index for the element in which the given Q i is non-zero. The element surface integral has been split into two components, the integrals for incoming and outgoing information, these are represented by Γ in e and Γ out e respectively. The incoming and outgoing surfaces are defined by

Γ in e = {r ∈ Γ e | Ω • n(r) < 0}, (11)
Γ

out e = {r ∈ Γ e | Ω • n(r) > 0}. ( 12
)
The sub-grid scale formulation imposes the incoming component of the fine solution to be zero which means the fine solution receives information from only the coarse solution. This decouples the coefficients of the fine scale solution unless they are within the same finite element. Thus, equation 10 becomes

Ve Q i Σ t φ dV + Ve Q i Ω•∇ φ dV + Ve Q i Σ t θ dV - Ve ∇Q i • Ω θ dV + Γ out e Q i (Ω • n) θ dΓ = Ve Q i S dV for i ∈ {1, 2, .., η Q }. ( 13
)

Spherical Harmonics Angular Discretisation

The spherical harmonics method is most commonly implemented with a uniform resolution over the entire space and energy domain. In this work, the size of the spherical harmonic expansion is allowed to vary across these dimensions. The fine and coarse components of the sub-grid scale solution are discretised using a variable order spherical harmonics expansion:

φj (Ω) = N φ j l=0 l m=-l Y m l (Ω) Φ j,l,m , (14) θj
(Ω) = N θ j l=0 l m=-l Y m l (Ω) Θ j,l,m . (15)
The real orthonormal spherical harmonics are represented by Y m l (Ω) and the flux expansion coefficients are Φ j,l,m and Θ j,l,m . The order of the angular expansion for a given node is defined by the integer sets N φ j and N θ j . These sets are also group dependent, N φ g,j and N θ g,j , but the subscript g has been neglected whilst considering a single group. The order of the expansion is restricted to odd integers due to the inferior accuracy and difficulty associated with even order approximations [START_REF] Davison | Neutron Transport Theory[END_REF]. In line with common usage, an expansion of order N is referred to as a P N expansion. The two indices for the spherical harmonics are combined into a single index to simplify notation, this is represented as

φj (Ω) = M φ j q=1 G q (Ω) Φ j,q (16)
θj (Ω) = M θ j q=1 G q (Ω) Θ j,q . (17)
The spherical harmonics are now represented by G q (Ω) and the index q runs through all spherical harmonics associated with a given expansion order. The coefficients of the expansion are represented by Φ j,q and Θ j,q for spatial node j and angular basis function q. The size of the angular expansion for the coarse and fine scale components at a given spatial node is given by the set M φ j and M θ j respectively. The expansion size sets are calculated from the sets N φ j and N θ j . In this work the two expansion size sets are allowed to differ i.e.

the discontinuous component of the solution may be represented with a different angular order to that of the continuous solution. A constraint is imposed that the discontinuous component must have an order greater than or equal to that of the continuous component at each node:

N θ j ≥ N φ C(j) ∀ j ∈ {1, ..., η Q }. (18)
The operator C : N → N is a surjective mapping from the discontinuous node numbers to the continuous node numbers. This constraint is imposed because the discontinuous component stabilises the solution. The absence of certain functions from the discontinuous solution destabilises the continuous components and produces oscillations in the final solution. The adaptive algorithm described in section 5.4 is responsible for forming the sets N φ j and N θ j in line with this constraint. The angular approximations are substituted into equations 9 and 13. To form a complete set of equations, the N i and Q i weighted equations are weighted with G p (Ω) for p ∈ {1, 2, ..., M φ i } and p ∈ {1, 2, ..., M θ i } respectively and integrated over the angular domain. This gives a system of equations which can be represented by

A Φ + B Θ = S Φ (19)
C Φ + D Θ = S Θ (20)
and given in matrix form as

A B C D Φ Θ = S Φ S Θ (21)
where Φ and Θ are vectors containing the coefficients Φ j,q and Θ j,q of the coarse and fine scale solutions respectively. The size of these vectors is given by

dim( Φ) = η N j=1 M φ j , (22) dim
( Θ) = η Q j=1 M θ j . (23)
The operator matrices and source terms can be considered as a set of angular sub-matrices and vectors:

A ij = V N i Σ ij t A N j dV - V ∇N i • A ij Ω A N j dV + Γ N i (A ij Ω A • n)N j dΓ, (24)
B ij = V N i Σ ij t B Q j dV - V ∇N i • A ij Ω B Q j dV, (25)
C ij = Ve Q i Σ ij t C N j dV + Ve Q i A ij Ω C • ∇N j dV, (26)
D ij = Ve Q i Σ ij t D Q j dV - Ve ∇Q i • A ij Ω D Q j dV + Γ e out Q i (A ij Ω D • n)Q j dΓ, (27)
S Φ,i = η Q j=1 V N i Q j S i Φ,j dV, (28)
S Θ,i = η Q j=1 V Q i Q j S i Θ,j dV. ( 29
)
where i and j are spatial node indices. The angular integrals have been absorbed into the matrix and vector terms (shown in bold). The angular operator matrices for matrix A are given by

Σ ij t A p,q = Σ t (r) 4π G p (Ω)G q (Ω)dΩ (30) A ij Ω A p,q = 4π Ω G p (Ω)G q (Ω)dΩ (31)
for p ∈ {1, 2, ..., M φ i } and q ∈ {1, 2, ..., M φ j }. The matrices for B, C and D are represented by the same equation, however, the range of the p and q indices differs:

B : p ∈ {1, 2, ..., M φ i } ; q ∈ {1, 2, ..., M θ j }, (32)
C : p ∈ {1, 2, ..., M θ i } ; q ∈ {1, 2, ..., M φ j }, (33)
D : p ∈ {1, 2, ..., M θ i } ; q ∈ {1, 2, ..., M θ j }. (34)
This means that the size of the matrices Σ ij t X and A ij Ω X (where X is one of A, B, C or D) is dependent on the spatial node indices i and j. An important modification of the angular matrices was required to provide a smooth and accurate solution. Mathematically stated, the modification was as follows:

A ij Ω A p,q = 0.0 ∀ p, q > min(M φ i , M φ j ), (35)
A ij Ω B p,q = 0.0 ∀ p, q > min(M φ i , M θ j ), (36)
A ij Ω C p,q = 0.0 ∀ p, q > min(M θ i , M φ j ), (37)
A ij

Ω D p,q = 0.0 ∀ p, q > min(M θ i , M θ j ). ( 38
)
This adjustment removes the coupling of higher order basis functions with lower order basis functions on a neighbouring node. For example, consider a P 1 node coupled with a P 3 node, the coupling between the l ≥ 2 moments of the P 3 node and the l ≤ 1 moments of the P 1 node is removed. A surmised justification for the removal this coupling is given based upon the eigenstructure of the matrices formed from the angular integrals, however, this is not a rigorous mathematical explanation. The streaming matrices, defined by equation 31, may be eigendecomposed for each order of expansion. The eigenvalues of each P N order define characteristic directions for the polynomial angular discretisation. These characteristic directions are different for each expansion order, thus two nodes with different P N orders have different characteristic directions. Figure 1 shows the non-overlapping characteristic directions for P1 and P3 one-dimensional discretisations. A discontinuous S N scheme with different angles at two neighbouring nodes would have to ensure continuity of partial current across the interface. It is believed that the removal of the coupling between the P N moments is analogous to the preservation of current between the nodes. The removal of the coupling implies that when streaming from a higher order node to a lower order node, the lower order node receives information from only the lower order coefficients of the higher order node. Therefore, the nodes then share the same characteristic directions. This reduces the angular resolution of the streaming from the higher order node but this must be done in order to couple with the lower order node. The difference observed with and without the coupling is demonstrated by a simple example in section 6. Note there is no requirement to adjust the total cross section matrix (equation 30) because there is no moment to moment coupling attributed to the orthonormal property of the basis functions. The angular source vector for the continuous weighted equation is represented as

S i Φ,j p = M θ j q=1 S j,q 4π G p (Ω)G q (Ω)dΩ (39)
for p ∈ {1, 2, ..., M φ i }. The source vector for the discontinuous weighted equation, S i Θ,j , has the same form as equation 39, however, the range of the weighting functions is p ∈ {1, 2, ..., M θ i }. The size of the angular source vector is dependent on the spatial node index i. The above source term has been obtained by approximating the spatial and angular dependence of the source as

S(r, Ω) = η Q j=1 Q j (r) M θ j q=1 G q (Ω) S j,q .
(40)

Boundary Conditions

The boundary conditions are imposed in the sub-grid scale through the surface integrals in equations 24 and 27. A Riemann decomposition technique is used to resolve the incoming and outgoing components of the spherical harmonic moments through a given surface [START_REF] Buchan | Riemann boundary conditions for the Boltzmann transport equation using arbitrary angular approximations[END_REF]. The method provides an efficient approximation to the half angle integrals required for particle transport boundary conditions. The projection of the streaming matrices in the direction of a given surface normal is defined as

A s = A Ω • n, (41)
where A Ω uses the same definition as equation 31 except M φ i = M φ j , therefore the component matrices are square. The incoming and outgoing information is separated in the surface integral of equation 24:

Γ N i A s N j dΓ = Γout N i A + s N j dΓ + Γin N i A - s N j dΓ. ( 42
)
The matrices A + s and A - s are the incoming and outgoing streaming matrices which depend upon the surface normal. The incoming and outgoing components are determined by performing an eigendecomposition on the surface streaming matrix:

A s = VΛV T ( 43
)
where V is the matrix of right eigenvectors of A s arranged in columns and Λ is a matrix with the eigenvalues of A s along the diagonal and zero elsewhere. The eigenvalues of the streaming matrix represent the flow of information across the boundary. Positive eigenvalues represent flow in the direction of the boundary normal (out) and negative eigenvalues represent flow in the direction opposite to the boundary normal (in). The eigenvalue matrix is split into two separate matrices, one containing the positive eigenvalues, Λ + , and one containing the negative eigenvalues, Λ -. The streaming matrix is then recombined in two parts, an outgoing and an incoming component:

A s = VΛ + V T + VΛ -V T , (44)
= A + s + A - s . (45)
The eigenstructure of the streaming matrix is different depending on the order of the spherical harmonic expansion. This means that the hierarchical property of the streaming matrix is lost once it is split into the two components. The variable order expansion therefore requires these matrices to be calculated for every order present in a calculation. The calculation of these matrices for every order and all normals n would cause a significant increase in computation and memory. This issue is treated differently for the internal and external boundary integrals. The external boundary integrals impose the boundary conditions upon the continuous solution and have a greater impact on the overall solution. For this reason, the external boundary matrices (present in matrix A) are calculated for each angular order required so that the boundary condition is imposed accurately. The internal boundary matrices (present in matrix D) are treated differently, these are calculated for the highest angle order in a given discretisation. The appropriate sub-matrix is then used to impose the boundary condition for smaller angular orders. This is an approximation, however, it does not introduce a large error because it is associated with the sub-grid scale component of the solution.

Solving the Sub-Grid Scale Equations

An expression for the fine scale solution is obtained by multiplying equation 20 from the left with D -1 :

Θ = D -1 S Θ -D -1 C Φ. ( 46
)
This is substituted into equation 19 to obtain

(A -BD -1 C) Φ = S Φ -BD -1 S Θ . (47)
The above equation (47) represents a system with the size of a continuous finite element discretisation. This is solved to find the coarse scale solution which is then modified by calculating the fine scale solution using equation 46. This technique results in a continuous finite element sized system but retains the accuracy and stability of a discontinuous system through the addition of the fine scale solution. However, the stability and accuracy is not provided at zero computational cost due to the inversion of the matrix D. The sub-grid scale formulation is constructed such that each unknown of Θ associated within an element is decoupled from the unknowns of all other elements. This allows the inversion of the matrix to be performed on an element by element basis which incurs significantly less computation. Equation 47 represents a single energy group equation, the multigroup system consists of a set of G equations of this form.

Spherical Harmonics Angular Error Indicator

The regular adaptive algorithm requires a driving force in order to determine a suitable angular resolution over the energy and space dimensions. The driving force in this work is an approximation of the error in the solution due to the angular discretisation. The aim is to approximate the true angular error at a given point in space:

(Ω) = ψ exact (Ω) -ψ N (Ω), (48)
where ψ N (Ω) represents the full sub-grid scale solution for a single energy group at a point in space with a spherical harmonics expansion of order N . The desired but unknown solution with zero angular discretisation error is ψ exact (Ω). An approximation to the true error is not available without undertaking further computation. Rather than approximate the true error, the error between an order N solution, ψ N (Ω), and an order N -2 solution, ψ N -2 (Ω), is approximated:

N (Ω) = ψ N (Ω) -ψ N -2 (Ω) (49) ≈ ψ N (Ω) -Y N -2 ( ψ N (Ω)) = N l=N -1 l m=-l Y m l (Ω) ψ l,m N . ( 50
)
The operator Y N performs an order N spherical harmonic expansion of the input function. This operation on a function composed of spherical harmonics with known coefficients returns the same coefficients as the original function:

Y N -2 ( ψ N (Ω)) = N -2 l=0 l m=-l Y m l (Ω) ψ l,m N . ( 51
)
The operation is chosen with N -2 because the discretisation uses only odd order expansions. The approximation N (Ω) is then used as the error to drive the adaptive algorithm. This error does not contain information of the higher order moments and it does not represent the true error. Despite this, it acts as an acceptable guide for the adaptive method because it provides information regarding the significance of the highest order moments of the solution.

The approximate error indicator is discretised by performing a spherical harmonic expansion of equation 50. The orthogonality of the spherical harmonics gives the coefficients/moments of the error expansion as

m l = 4π Y m l (Ω) (Ω)dΩ (52) = ψ l,m N for N -1 ≤ l ≤ N 0 otherwise (53)
Therefore, the non-zero moments of the angular error are defined by the higher order flux moments. This is an intuitive measure for the error because it utilises the fact that the coefficients of the higher order basis functions will go to zero as the expansion is sufficient to describe the angular distribution.

Goal-Based Error Estimator

Functional Error Estimator

The angular error estimator described in section 3 gives a measure of the error in the solution attributed to the angular approximation. This error measure is suitable for adapting the discretisation in order to reduce the error in the solution over the entire domain. In contrast, goal-based adaptive algorithms require an error estimator for a given goal of the calculation. This type of error estimator provides information regarding the contribution to the error from across the particle phase-space.

Consider the transport equation in fixed source (2) and eigenvalue form (3) using operator notation:

Lψ exact = S e , (54)
Lψ exact = λPψ exact , (55)
where L contains the streaming, total reaction and scattering operators, P is the fission operator. The residual for each equation is given by

R(ψ) = Lψ -S e , (56)
R(ψ) = Lψ -λPψ, (57)
where the residual is zero for the exact continuum solution of the equation, R(ψ exact ) = 0. The goal of the calculation is generalised to an arbitrary functional of the solution:

F (ψ) = f (ψ)dP ( 58
)
where dP represents an infinitesimal region of phase space. The integral over the entire phase space is

dP = Emax 0 dE 4π dΩ V dV ( 59
)
where E max is the maximum energy of the particles in a given model. The following derivation provides an error estimator for the arbitrary functional. The integrand of the functional is approximated by a Taylor series expansion about the exact solution:

f (ψ exact ) -f (ψ) = ∂f ∂ψ (ψ exact -ψ) + O((ψ exact -ψ) 2 ), (60)
≈ ∂f ∂ψ (ψ exact -ψ). ( 61
)
The residual is also approximated by a Taylor series expansion around the exact solution and rearranged to find an expression for the error in terms of the residual and its derivative:

R(ψ exact ) -R(ψ) ≈ ∂R ∂ψ (ψ exact -ψ), (62)
⇒ ψ exact -ψ ≈ ∂R ∂ψ -1 (R(ψ exact ) -R(ψ)). ( 63
)
Equation 63 is substituted into equation 61 to acquire

f (ψ exact ) -f (ψ) ≈ ∂f ∂ψ ∂R ∂ψ -1 (R(ψ exact ) -R(ψ)). (64)
Integrating equation 64 over all dimensions gives the error in the functional:

f (ψ exact )dP -f (ψ)dP ≈ ∂f ∂ψ ∂R ∂ψ -1 (R(ψ exact ) -R(ψ))dP, (65)
⇒ F (ψ exact ) -F (ψ) ≈ ∂f ∂ψ ∂R ∂ψ -1 (R(ψ exact ) -R(ψ))dP. (66)
The definition of an adjoint operator is used to find

F (ψ exact ) -F (ψ) ≈ (R(ψ exact ) -R(ψ)) ∂R ∂ψ * -1 ∂f ∂ψ dP (67)
and the adjoint equation is defined as

∂R ∂ψ * ψ * exact = ∂f ∂ψ . ( 68
)
The expression for the adjoint solution gives a concise expression for the error in the functional:

F (ψ exact ) -F (ψ) ≈ (R(ψ exact ) -R(ψ))ψ * exact dP, (69) = -R(ψ)ψ * exact dP. ( 70
)
This is modified by adding and subtracting the inner product of the forward residual and an arbitrary function, ψ * , to give

F (ψ exact ) -F (ψ) = -R(ψ)(ψ * exact -ψ * )dP -R(ψ)ψ * dP. ( 71
)
The forward solution, ψ, will be obtained through a weighted residual technique such as the method described in section 2. In this case, if the function ψ * is in the space of the weighting functions it holds that

R(ψ)ψ * dP = 0. ( 72
)
In such a situation, the error measure is given as

F (ψ exact ) -F (ψ) = -R(ψ)(ψ * exact -ψ * )dP. ( 73
)
This equation states that the error in the functional is equal to the error in the adjoint solution weighted by the residual of the forward equation. This will be referred to as the adjoint error estimator. The next task is to demonstrate that an equivalent form of this equation can be obtained which uses the error in the forward solution weighted by the residual of the adjoint equation. This second form will be referred to as the forward error estimator. The following steps are different for the fixed source and eigenvalue equations and will be shown separately.

For the fixed source case, it is first noted that

-R(ψ) = s -Lψ = L(ψ exact -ψ). (74)
This is substituted into equation 73 to find

F (ψ exact ) -F (ψ) = L(ψ exact -ψ)(ψ * exact -ψ * )dP (75) = L * (ψ * exact -ψ * )(ψ exact -ψ)dP (76) = -R * (ψ * )(ψ exact -ψ)dP (77)
where R * is not the adjoint of R but it is the residual of the adjoint equation:

-R * (ψ * ) = L * ψ * exact -L * ψ * = ∂f ∂ψ -L * ψ * . ( 78
)
In the case of the eigenvalue equation, we re-insert the residual of the exact solution:

F (ψ exact ) -F (ψ) = -R(ψ)(ψ * exact -ψ * )dP + R(ψ exact )(ψ * exact -ψ * )dP (79) = (R(ψ exact ) -R(ψ))(ψ * exact -ψ * )dP. ( 80
)
The residual is defined as a linear operator by fixing the eigenvalue as the exact value and equation 80 can be arranged to the form

F (ψ exact ) -F (ψ) = R(ψ exact -ψ)(ψ * exact -ψ * )dP (81) = -R * (ψ * )(ψ exact -ψ)dP. ( 82
)
Here the operator R * is defined as the residual of the adjoint eigenvalue equation:

R * (ψ * ) = L * ψ * -λ * P * ψ * (83)
where λ * is the eigenvalue of the adjoint equation. It is known that the eigenvalues of the forward and adjoint equations are equal, λ = λ * .

Residual Calculation

It should be noted that one must be careful when calculating the functional error given in equation 77. If the solution is calculated using a weighted residual technique and the approximation to the error is in the same phase space discretisation, the functional error will compute to zero. In this work, the error approximation defined in section 3 is in the same phase space discretisation. To avoid this problem, the residual is calculated using an approximation. The residual of the transport equation is the difference between the left and the right hand side of the equation for some approximation of the solution. That is, the equality in equations 2 and 3 holds only for the exact solution. The residual can be approximated by operating on some solution with a more accurate operator than that used to obtain the given solution. For example, given a solution from a P 3 angular approximation, an approximate residual may be obtained by substituting this solution into a P 5 set of equations. A similar idea is employed in this work, the residual is approximated by reducing the accuracy of a given approximation and substituting it into the system of equations from which it was the solution.

A discrete residual is obtained by weighting the continuum residual with the solution basis functions:

R φ i,p = V N i (r)dV G p (Ω)R(ψ)dΩ ∀i ∈ {1, ..., η N } ∀p ∈ {1, ..., M φ i }, (84)
R θ i,p = V Q i (r)dV G p (Ω)R(ψ)dΩ ∀i ∈ {1, ..., η Q } ∀p ∈ {1, ..., M θ i }. ( 85
)
This is the system of weighted residuals used to obtain the solution and will therefore equal zero. The adaptive algorithm requires information regarding the angular error, therefore a residual which takes into account this error is desirable. This is achieved by adjusting the solution and calculating the residuals stated above. The solution is modified differently depending on the spatial weighting function:

R φ i,p = V N i (r)dV G p (Ω)R( ψ i )dΩ, (86)
R θ i,p = V Q i (r)dV G p (Ω)R( ψ C(i) )dΩ, (87)
where

ψ i = η N j =i N j (r) M φ j q=1 G q (Ω)Φ j,q + N i (r) M φ i q=1 G q (Ω)Φ i,q + η Q j=1 C(j) =i Q j (r) M θ j q=1 G q (Ω)Θ j,q + η Q j=1 C(j)=i Q j (r) M θ j q=1 G q (Ω)Θ j,q . ( 88
)
where C : N → N is the same mapping defined in section 2.3. The modified angular expansion sizes, M φ i and M θ i , are equal to the expansion sizes for N φ i -2 and N θ i -2 respectively. This sets the highest order moments to zero at the node associated with the spatial weighting function and retains the full solution elsewhere. This residual provides an insight into the importance of those higher order moments at the given point in space. The two components of the residual, R φ and R θ , are combined to form a discontinuous residual which takes into account both sets of equations. They are combined in the same manner as the two components of the flux solution, Φ and Θ:

R i,p = R φ C(i),p + R θ i,p ∀i ∈ {1, ..., η Q } ∀p ∈ {1, ..., M θ i }. (89)
The constraint stated by equation 18 means that M θ i ≥ M φ C(i) . At nodes where the discontinuous resolution is greater than the continuous resolution the higher order moments of the continuous residual are zero:

R φ C(i),p = 0.0 ∀ p > M φ C(i) . ( 90
)
The adaptivity algorithm utilised in this work requires a metric on a continuous space, therefore, the residual is also required in the continuous space. The discontinuous residual is mapped to a continuous space by retaining the maximum of the discontinuous values for each continuous node.

Adaptive Procedure

The error measures defined in sections 3 and 4 provide information regarding the quality of the solution over the space and energy dimensions. This information is used to adapt the variable order angle expansion in order to minimise the error in the solution. Two adaptive strategies have been employed in this work: (i) adaptivity to minimise the global error in the solution, this will be referred to as regular adaptivity, (ii) adaptivity to minimise the error in a given quantity of interest, known as goal-based adaptivity. The algorithms for regular and goal-based adaptivity are similar and will be described together.

Continuous Error Indicator

The angular error indicator described in section 3 provides a discontinuous error field for each energy group:

g (r, Ω) = η Q j=1 Q j (r) M θ j q=1 G q (Ω) g,j,q . (91)
The group dependent subscript will be included for the purpose of this section. The adaptive algorithm in this work requires a continuous error. For this reason, the discontinuous error field is reduced to a continuous error field:

ˆ g,i,q = max j∈Di { g,j,q }. ( 92
)
The set D i contains the discontinuous node numbers for continuous node i:

D i = {j ∈ {1, 2, ..., η Q } : C(j) = i}. ( 93
)
The maximum error value is used in order to resolve the angular approximation satisfactorily for all discontinuous nodes at a point in space. This gives a continuous field for the approximate angular error for each energy group.

Nodal Error for Regular Adaptivity

A new nodal error quantity is defined for use in the adaptive algorithm, for regular adaptivity this quantity is defined

e g,i = max q∈{1,2,...,M φ i } |ˆ g,i,q | ∆ t ψ max sc,g (94)
where ∆ t is a user-defined error tolerance and ψ max sc,g is the maximum value of the scalar flux for energy group g. The maximum value of the group scalar flux is used to account for the difference in magnitude of the flux between energy groups. This technique of weighting the group errors was used by Wang et al. [START_REF] Wang | Three-dimensional h-adaptivity for the multigroup neutron diffusion equations[END_REF]. The ultimate value for the error at a given spatial node and energy group is chosen as the maximum value of the angular moments of the error field. The values of the moments are deemed comparable because the spherical harmonic basis functions are orthonormal. The normalisation ensures that each function has a similar contribution to the overall angular distribution of the field.

Nodal Error for Goal-Based Adaptivity

For goal-based adaptivity, two nodal error quantities based on equations 77 and 73 are defined:

êg,i = max q∈{1,2,...,M φ i } |ˆ g,i,q R * g,i,q | ∆F , ( 95
) ê * g,i = max q∈{1,2,...,M φ i } |ˆ * g,i,q R g,i,q | ∆F ( 96
)
where R * g,i,q and R g,i,q are components of the adjoint and forward residual vectors respectively, ˆ * g,i,q is the continuous error indicator for the adjoint solution, ∆F is the desired error in the functional per degree of freedom. This is defined by

∆F = ∆F N dof ( 97
)
where ∆F is a user-defined error tolerance for the goal and N dof is the number of degrees of freedom in the problem. The two nodal error quantities, êg,i and ê * g,i , will be referred to as the forward and adjoint error measures respectively. The two error estimators individually encompass the total error in the goal functional, however, the distributions of the error estimators over phase-space are not identical. There is not a clear best method to combine these errors and therefore a variety of combinations are tested in this work. The combinations used will be noted in the individual examples.

Adaptive Algorithm

The steps to perform an adaptive calculation are shown in algorithm 1. The angular discretisation is adapted N a times in an iterative type scheme. Each iteration initially solves the transport equation to obtain the solution (including the adjoint for goal-based adaptivity). The solution is used to calculate the error field which is consequently used to adapt the angular discretisation. The regular and goal-based error definitions are such that the nodal values should be less than 1.0 if the user-defined tolerance has been met. If the error is greater than one then the resolution is increased and if the error is much less than one than the resolution is decreased. The order of the angular expansion is increased and decreased in multiples of two because only odd order expansions are being used. A maximum angular expansion order, N max , is enforced to prevent the calculation becoming too computationally expensive. If the error is greater than twice the user-defined tolerance then the angular expansion order is increased by four orders. This condition was included to reduce the number of adaptive steps required to reach convergence. If this condition causes over resolution then a reduction in expansion order will occur in the next adaptive iteration. The tolerances used for lowering and doubly refining the resolution are heuristic and therefore are not optimum. The situation could occur in which the adaptive procedure oscillates between two states because the tolerances are such that the error cannot lie between them. The simplest method to address this would be the modification of the tolerance used for reducing the order.

A smoother is applied to the angular discretisation in order to prevent large steps in expansion order between neighbouring nodes. The smoothing algorithm modifies the discretisation by calculating a new P N order for each node. It calculates the new order by taking the average of the orders on a set of nodes that compose a patch surrounding the node under consideration. The patch of nodes includes all nodes which have a connecting edge with the node under consideration.

It has been noted various times throughout the article and is shown in algorithm 1 that the angular discretisation is adapted on the coarse (continuous) mesh i.e. the error field is used to construct N φ as opposed to N θ . In this work, the fine scale angular discretisation, N θ , is constructed from the coarse angular discretisation. The fine scale angular discretisation is set such that the angular order is constant within each element. The size of the chosen angular order is equal to the maximum value of the coarse scale nodes composing the element. The method is shown in algorithm 2. The reason for this choice is that the fine scale matrices are always square and therefore require no adjustment as shown in equation 38. With this choice of adaptive algorithm, the fine scale expansion order is always greater than or equal to the coarse scale expansion order.

Algorithm 1 The algorithm used to perform an adaptive calculation with an adaptive angular discretisation. The number of adaptive steps is N a . The formulation of the variable order spherical harmonics method was described in section 2.3. It was noted that the coupling of angular moments between neighbouring nodes of different angular order is modified in the assembly. The coupling between angular moments at a given node and higher order moments at a neighbouring node is removed. This section demonstrates the necessity for the removal of the coupling through a simple one dimensional example.

The example is a mono-energetic fixed source problem in a homogeneous medium with vacuum boundary conditions. The medium exists in the region 0.0 < x < 1.0 and is a pure absorber with Σ a = 0.8. A non-uniform angle expansion order was applied to the one-dimensional space. The angular order was set as P 1 for 0.4 < x < 0.6 and P 3 elsewhere. The calculation with the mixed P 1 /P 3 discretisation was performed with (i) full coupling between moments (ii) modified coupling between moments at nodes of different order. The scalar flux solution for each calculation is shown in figure 2 alongside uniform P 1 and P 3 solutions.

The mixed P 1 /P 3 solutions lie between the uniform P 1 and P 3 solutions which is intuitively expected. The scalar flux obtained with full coupling between the moments shows a significant unphysical discontinuity in the solution between the regions of different expansion order. The solution obtained with the modified coupling does not possess the discontinuity. Evidently, the discontinuity vastly reduces the accuracy of the mixed P 1 /P 3 solution and appears unphysical. Therefore, the coupling was removed from the variable order formulation to remove these discontinuities and improve the accuracy of the solutions.

Applications

This section demonstrates the adaptive procedure applied to a selection of examples. The quality of the regular adaptive algorithm is quantified by taking the 2-norm of the error in the scalar flux:

ε sc,g = V 4π ψ ref g (r, Ω)dΩ - 4π ψ g (r, Ω)dΩ 2 dV = V ψ ref sc,g (r) -ψ sc,g (r) 2 dV , (98)
where ψ ref g is the reference angular distribution. The quality of the goal-based adaptivity is quantified by the absolute magnitude of the error in the goal functional:

ε F = |F ref -F | ( 99
)
where F ref and F are the reference and approximate values for the goal functional respectively. These quantities are plotted against the number of continuous degrees of freedom (CDOF) required for the solution. Each data point for these plots is obtained by converging the adaptive process for a different error tolerance (∆ t or ∆F ). This should indicate that the resulting discretisation is the best it can obtain for the given tolerance. This type of comparison demonstrates the accuracy of the solution for a given number of degrees of freedom. In goal-based calculations, the number of degrees of freedom is that of the forward solution only.

For some examples the errors are plotted against CPU time to observe whether the adaptive methods provide increased accuracy for a reduced computational time. The data points for the adaptive methods are the errors obtained at each adaptive iteration plotted against the cumulative time for the adaptive calculation. That is, the time for the n th data point is the total time for n calculations. The adaptive calculations have been chosen for a suitably low tolerance that the error from the adaptive method decreases to a comparable level to the uniform calculation.The data for the uniform calculations is presented in two manners: (i) each data point is the time for a single calculation, and (ii) each point is the cumulative time for the set of uniform calculations. The adaptive calculations provide a convergence test through one calculation which consists of multiple solves. The uniform cumulative times are plotted to represent the time for the convergence test that a user would have to undertake through uniform refinement. However, the cumulative times depend on the number of data points in the uniform convergence test so this data is present only as a guide.

The reference solutions were obtained for these examples by performing a high order P N calculation. These calculations were performed with a uniform expansion over the entire space and energy domain. The reference solutions are not necessarily fully converged in the angular dimension. The adaptive algorithm is designed to refine the angular resolution to obtain the fully converged solution and thus may attempt to go beyond the accuracy of the reference solution. This was prevented by setting the maximum angle expansion order to that of the reference solution. This means the adaptive algorithm will converge to the reference solution for these examples i.e. the error will go to zero as the number of degrees of freedom approaches that of the reference. The reference solutions were calculated using the same spatial mesh and energy group schemes and thus the error in each of these discretisations is zero with respect to the reference.

Problem 1

The first example shown is known as the Reed cell model. It is a one-dimensional monoenergetic fixed source example. The geometry of the problem is shown in figure 3 and the associated material properties are given in table 1. The reference solution for the example is shown in figure 4a. There are two particular regions of the solution worth noting. The strong source and absorber in region I leads to an isotropic solution for the majority of this region. The void section in region III leads to very directional flux in the negative (left-facing) direction.

I

Vacuum IV II V III Reflect 10 -8 0.0 0.0 IV 0.1 1.0 0.9 1.0 V 0.1 1.0 0.9 0.0 Table 1: The material properties for the different regions in problem 1. Region numbers refer to those in figure 3. Cross sections in units of cm -1 and source in units of cm -3 s -1 .

Regular Adaptivity

The convergence plot for the error in the scalar flux is shown in figure 5. The adaptive results are plotted alongside the uniform angular refinement for comparison. The vertical dashed blue line marks the number of CDOF required for the reference solution. As the error tolerance is increased the adaptive method will tend to a uniform P 1 solution, this is why the curves meet at the number of CDOF for a P 1 solution. The sharp decrease in error at the end of the adaptive results is due to the adapted discretisation tending to the uniform P 21 solution. Note the error tends to zero with a smaller number of unknowns than required for the P 21 solution (marked by the dashed blue line). This is because the adapted discretisation retains a P 1 expansion in the section of the domain with an isotropic flux, region I. This demonstrates that the flux is converged with a P 1 expansion in this region. An example of the adapted angular discretisation is shown in figure 4b. This plot shows the distribution of the expansion size over the domain of the problem. It is observed that the adaptive algorithm has placed the maximum P 21 expansion in the void region (III) and a P 1 expansion in the isotropic region (I). The error as a function of time for regular adaptivity is shown in figure 6a. It is observed that the regular adaptivity requires more time to obtain a given accuracy compared to a single fine resolution calculation. The reason the adaptive method performs poorly with respect to CPU time is attributed to the fact the P 21 calculation takes approximately four times longer than the P 1 calculation. This small time difference means that after four adaptive iterations the run-time will definitely be greater than that of a single P 21 calculation. It takes less iterations to surpass this time because each adaptive solve uses resolution greater than P 1 . The cumulative time for the uniform calculations is longer than that of the regular adaptivity for the same error but there are more data points (more solves). The goal used for this example was the total absorption in regions IV and V, as illustrated in figure 7. A variety of combinations of the forward and adjoint error measure were tested for this problem. The following errors were trialled for the goal-based calculations:

Goal-Based Adaptivity

e g,i = êg,i Forward, ( 100
)
e g,i = ê * g,i
Adjoint, (101)

e g,i = min{ê g,i , ê * g,i } Minimum, ( 102
)
e g,i = max{ê g,i , ê * g,i } Maximum, ( 103
)
e g,i = 0.5(ê g,i + ê * g,i ) Average. ( 104
)
The forward and adjoint errors were tested individually to demonstrate that both measures are reliable for adaptivity. The three combinations were tested to observe how the convergence was affected by each. Figure 8a shows the convergence of the error in the goal functional for the forward and adjoint error measures compared to the uniform refinement and the regular adaptivity. It is observed that a smaller error is obtained for fewer degrees of freedom when using the goal-based adaptivity compared to the regular adaptivity. The convergence plot for the combined error measures is shown in figure 8b. The results show that adaptivity using the average or the maximum error perform better than the minimum error. The convergence of the maximum error combination is the most consistent. The error versus CPU time for goal-based adaptivity is shown in figure 6b. This shows regular adaptivity out-performs goal-based adaptivity in this problem, whereas a single uniform calculation out-performs both adaptive methods. The reason for poor performance compared to uniform resolution was stated when discussing the time plot for regular adaptivity. The reason regular adaptivity performs better than goal-based is because the choice of goal requires that a large portion of the domain is resolved. Therefore, the adapted discretisations from goal-based and regular adaptivity are not largely dissimilar, however the goal-based method requires the added time for the computation of the adjoint solution. The goal-based method will only perform better than regular if it produces a significantly coarser resolution for a given accuracy.

Problem 2

This example is a three-dimensional mono-energetic fixed source problem. It is the third of the Kobayashi benchmarks documented by the NEA [START_REF] Kobayashi | 3-D Radiation Transport Benchmark Problems and Results for Simple Geometries with Void Regions[END_REF]. This benchmark is a threedimensional duct problem with 50% scattering in the region outside the duct. An illustration of the 3D geometry is shown in figure 9. The full dimensions of the problem can be found in the benchmark documentation. The material properties of the regions are shown in table 2. 9. Cross sections in units of cm -1 and source in units of cm -3 s -1 .

The reference P 11 solution is displayed using a contour plot shown in figure 10. A log scale has been used for the contour and colour map so that the flux in the duct can be observed. It is observed that the larger values of the flux extend down the duct in direct sight of the source. The effect due to the duct is noticeable in a number of regions where the isosurface bends in the directions of the duct.

Regular Adaptivity

The convergence of the error in the scalar flux is shown in figure 11. It was observed that the adaptivity greatly reduced the error in the scalar flux solution for a given number of degrees of freedom. The performance with respect to CPU time is shown in figure 12a. It is observed that the regular adaptivity performs slightly better than the single uniform calculation at the error of the P 9 calculation.

Goal-Based Adaptivity

The goal in this problem was the average flux over a 5 cm × 10 cm × 5 cm region at the end of the duct as shown in figure 13. As before, the goal-based adaptivity was performed using a number of error measure combinations.

The functional convergence results obtained from the forward and adjoint error measures are shown in figure 14a and those from the combinations are shown in figure 14b. These are shown alongside the results of uniform and regular adaptive resolution. The uniform P 5 result has a smaller error than the P 7 and P 9 results because the value of the functional does not converge monotonically with angular refinement. Therefore, the error in the functional goes to zero as the functional value oscillates around the true value during convergence. The P 5 result is thus not representative of the error at that resolution. It is evident that the regular adaptivity performed less well than uniform resolution of the angular discretisation. This is because the magnitude of the solution in the goal region is small in comparison to the solution elsewhere in the domain. The error measure employed for regular adaptivity is absolute as opposed to relative, therefore, the regular adaptivity is placing angular resolution in the regions of the domain with larger solution values. Thus, the regular adaptivity will not apply angular resolution in the area of this functional and hence the error is large in comparison to other methods. The results from goal-based adaptivity using the forward and adjoint errors have similar performance to one another at low refinements. Both error measures show effects due to the non-monotonic convergence of the functional. Aside from these oscillatory results, both error measures follow a similar trend. The combinations of the error measures also show similar performance, however, the most consistent and effective combination appears to be that which uses the maximum value. For this reason the error measure combination using the maximum values will be used in the remaining examples. The goal-based results show that a more accurate functional can be obtained for a given number of degrees of freedom. The improvement in error is not that significant because the absorbing medium has 50% scatter cross section. Therefore, the regions of the absorber media surrounding the duct are also important to the value of the functional meaning that the angular adaptivity must also resolve those regions.

The error is plot as a function of time in figure 12b. It is hard to distinguish a difference in performance between the methods due to the non-monotonic convergence of the functional. The adaptive methods pass through the same non-monotonic behaviour during the adaptive process which leads to the error staying roughly constant for several iterations. The regular adaptive method matches the performance of the goal-based method. This shows that the reduction in the degrees of freedom acquired through the extra computation of the adjoint has reduced the CPU time. However, it has not provided adequate benefit to out-perform the regular adaptivity.

Problem 3

This example has been used to illustrate the application of the energy dependent angular resolution. The problem geometry is shown in figure 15. It is a two group 2D problem with a homogeneous material throughout the domain. The material properties for each region are provided in table 3. The reference P 11 group solutions are shown in figures 16a and 16b.

Regular Adaptivity

This problem was performed using regular adaptivity with both an energy-dependent and an energy-independent angular resolution. In the case of the energy-independent calculation, the group-wise nodal errors defined by equation 94 are combined by taking the maximum value of the group errors:

e i = max g∈{1,...,G} e g,i . (105)
This gives a single error value for each continuous spatial node in the mesh which is consequently used for the adaptive procedure. The convergence of the error norm of the scalar flux for each method is shown in figure 17. The results show that both the energy-independent and energy-dependent adaptivity provide a smaller error for a given number of degrees of freedom compared to uniform resolution. In addition, it can be seen that the energy-dependent angular resolution has performed better than the energy-independent. Examples of the resulting angular resolutions from the energy-dependent and energy-independent adaptivity calculations are shown in figure 18. The figures are shown for calculations that used the same user-defined error tolerance. The better performance of the energy-dependent calculation is largely due to the isolated group 2 flux in the top right section of the domain. The energy-independent adaptivity places angular resolution in this region for both group 1 and group 2, whereas, the energy-dependent adaptivity places resolution only in group 2. It can also be noted that there is also less group 2 resolution around the group 1 source.

Goal-Based Adaptivity -Goal 1

Two different goal-based calculations are being presented to demonstrate the application of the energy-dependent angular resolution. As alluded to in the previous example, the forward and adjoint error measures are combined by taking the maximum values. The goal of the first set of calculations was the absorption rate of group 2 neutrons in the group 1 source region (region II in figure 15):

F = VII Σ a,2 (r)dV 4π ψ 2 (r, Ω)dΩ (106)
where V II is the volume of region II and Σ a,g is the macroscopic absorption cross section for group g. The goal-based adaptivity was carried out using both energy-independent and energy-dependent angular resolution. The results are compared with uniform resolution and regular adaptivity, the convergence of the functional for each method is shown in figure 19a.

The results show that the goal-based calculations typically provide smaller errors in the target functional for a smaller number of unknowns than the regular adaptivity or uniform resolution. The energy-dependent regular adaptivity initially shows worse performance than the energy-independent case. This occurs because the energy-dependent calculation will place less angular resolution for group 2 in region II because the magnitude of the flux is smaller. This effect disappears as the error tolerance is decreased and region II is resolved via the group dependent regular adaptivity. Examples of the adapted angular resolution obtained from the goal-based methods are shown in figure 20. The goal-based methods primarily resolve the region around the group 1 source because it has the greatest impact on the value of the goal functional. Figures 20b and20c show that the goal-based approach is increasing the resolution of both group 1 and group 2 around the group 1 source. This is to be expected because the group 1 down-scatter provides the source of group 2 neutrons in this region. Therefore, the group 1 solution must be accurate in order to determine the goal functional accurately. This is the reason there is no significant improvement in using different angular resolutions for each energy group within the goal-based approach. It is observed that the goal-based adaptivity attempts to resolve the contribution to the functional from the group 2 source. This contribution is small and is resolved only when the error tolerance is set low. The group-dependent angular resolution demonstrates that this is not required for group 1 which leads to the slight increase in performance of the group-dependent calculation at high resolution.

Goal-Based Adaptivity -Goal 2

The goal of the second set of calculations was the absorption rate of group 1 neutrons in the group 1 source region (the same region as goal 1, shown in figure 15): The convergence results for goal 2 are shown in figure 19b. As expected, given the nature of the goal, the goal-based methods again show superior convergence when compared to the regular adaptivity and uniform resolution. The notable difference between goal 2 and goal 1 is the benefit gained by the energy-dependent angular resolution. It is noted that both the regular and goal-based adaptive methods provide superior results when using the energydependent angular resolution. This occurs because the functional depends solely on the group 1 flux in the group 1 source region. The energy-dependent regular adaptivity is more efficient because the group 2 resolution will not increase as quickly as the group 1 resolution in the goal region and hence less unknowns are used. As the functional does not depend on the group 2 resolution, the accuracy of the functional is not impacted. The energy-dependent goal-based adaptivity exhibits large improvements in convergence compared to the other methods. This performance occurs because minimum resolution (P 1 ) is used for group 2 as the functional is not dependent on group 2. Examples of the angular resolutions resulting from the goal-based adaptive methods are shown in figure 21. Figure 21c demonstrates the minimum resolution applied to the whole of the spatial domain for group 2.

F = VII Σ a,1 (r)dV 4π ψ 1 (r, Ω)dΩ. (107)
The problem and the goals presented here were designed in order to demonstrate the energy-dependent capability of the goal-based method. Clearly, such functionals could be calculated efficiently by splitting the spatial domain and/or performing a single group calculation. However, the problem was chosen to demonstrate the purpose behind the methods.

Problem 4

This example is an eigenvalue problem, the model is a single fuel pin immersed in water. The calculation was performed with a two-group energy scheme. The model is subcritical with 1/λ = k ≈ 0.019. The geometry for the model is shown in figure 22 and the material properties for each region are provided in table 4. It is a single fuel pin model and vacuum boundary conditions are applied to the external boundaries. The reference solution for each energy group is shown in figure 23.

Regular Adaptivity

The convergence of the scalar flux error norm is shown in figure 24a. It is seen the group-independent regular adaptivity provides no benefit compared to that of the uniform refinement. However, the group-dependent angular resolution does provide a smaller error for a given number of degrees of freedom. The group-dependent resolution provides an increase in performance because the thermal group (2) requires less resolution compared to the fast group [START_REF] Ackroyd | Composite finite element solutions for neutron transport[END_REF]. Examples of adapted angular resolutions are shown in figure 25. It is seen that the adaptivity has placed more resolution in the fast group compared to that of the thermal The eigenvalue was used as the goal functional of the calculation. The convergence plot for the eigenvalue is shown in figure 24b. The energy-independent adaptive methods, both regular and goal-based, have provided no benefit compared to the uniform resolution.

However, both the regular and the goal-based energy-dependent adaptive methods provide a lower eigenvalue error for a given number of degrees of freedom. The regular adaptivity equals the performance of the goal-based adaptivity for both the energy-dependent and independent cases. This demonstrates that the extra computation of the adjoint has provided negligible benefit. This is most likely a consequence of the fact that the eigenvalue is a functional that depends on the flux over the entire phase space. This is in contrast to the goal functionals in the previous examples which have depended on a specific region of space or energy group. As the eigenvalue depends on the entire solution it is intuitive that the regular adaptivity, which aims to lower the error on the solution over the entire domain, is effective at reducing the functional error.

Examples of the angular discretisation obtained through the goal-based adaptivity are shown in figure 26. The angular discretisations presented for the regular and goal-based methods each have a similar number of CDOF. In both the regular and goal-based cases, the energy-independent discretisation has a similar number of CDOF but produces a significantly larger error. Both types of adaptivity demonstrate that placing greater resolution in the fast group provides a more accurate eigenvalue.

Problem 5

The previous example demonstrated a criticality safety type problem with a single fuel pin immersed in a body of water. Typical fuel pin calculations are used to model an infinite lattice of pins. It was shown by Adams that the fast flux in an infinite lattice has a very irregular angular distribution [START_REF] Adams | Angular dependence of the fast flux in reactor lattices[END_REF]. This example is the same as that presented by Adams, it is a single group fixed source calculation representative of the fast flux in a reactor calculation. The fuel pin is represented as a square for easier code to code comparisons, the geometry is shown in figure 27. The material properties for the fuel (I) and the moderator (II) are shown in table 5. The problem was modelled using an eighth of the geometry on a coarse mesh (146 triangles) to allow the use of high angular resolution within a reasonable time period. The reference solution was obtained from a P 31 calculation which uses 528 basis functions. The scalar flux solution is shown in figure 28. The discontinuities visible in the scalar flux solution are due to the coarse mesh. 27. Cross sections in units of cm -1 and source in units of cm -3 s -1 .

Region

A uniform P 59 (1830 basis functions) calculation was performed in order to look at the azimuthal flux plots as performed by Adams, these are shown in figure 29. Figure 29a displays azimuthal plots for the same two locations shown by Adams. The plots demonstrate that the P 59 is not adequate to capture the azimuthal distribution as accurately as the product quadrature set used by Adams, however, the overall flux shapes are similar. The plots shown in figure 29b show that the flux is smoother in the bulk of the fuel and the moderator. Despite this smoother shape, there are sharp features present which require high order polynomials to be adequately resolved. Goal-based adaptivity is not used in this example because this is primarily to test the ability of the method to resolve the overall solution.

Regular Adaptivity

An example of the angular discretisation obtained through the regular adaptive method is shown in figure 30. The most angular resolution has been placed along the interface between the fuel and moderator. The least resolved regions are in the bulk of the fuel and the moderator.

The convergence of the error in the scalar flux is shown in figure 31a. The adaptive method demonstrates that it can provide a smaller error for a given number of degrees of freedom. The behaviour of the last data point is attributed to the fact it is affected by the P 31 limit imposed on the adaptive scheme. The performance with respect to CPU time is shown in figure 31b. In this example the regular adaptive method is capable of providing a solution with a smaller error for less computational time. The difference between the uniform method and the adaptive method diminishes towards the higher resolution. This is also attributed to the P 31 limit. As the adaptive method saturates P 31 resolution over the entire domain it will no longer have any benefit over the uniform resolution calculation.

Conclusions and Further Work

A variable order spherical harmonics scheme has been described in which the order of the angular expansion varies over the space and energy dimensions. This was implemented using the inner element sub-grid scale finite element spatial discretisation using the multigroup energy scheme. In this discretisation, the order of the expansion varies between each node of the finite element mesh for each energy group. Adaptive angular resolution techniques have been developed in the variable order scheme with the intention of increasing the efficiency of the solution process. Two types of adaptive method have been presented: (i) regular adaptivity, in which the global error in the solution is minimised, and (ii) goal-based adaptivity, in which the error on a specified functional is minimised. An error measure has been given for each type of adaptive method and a procedure has been presented which utilises the error measures to form an adapted angular resolution over the phase-space.

The adaptive methods have been tested with a range of example applications which include both fixed source and eigenvalue cases. The regular adaptivity provided a more accurate solution for a given number of degrees of freedom for the majority of examples. It performed well in the fixed source shielding applications but provided less benefit for the eigenvalue calculation. The example in section 7.4 demonstrated that the energy dependence of the angular resolution was necessary to gain any benefit from the regular adaptivity. The reduced performance in the eigenvalue calculations indicate that these examples require a relatively uniform distribution of high angular resolution over the entire domain.

The goal-based adaptivity also demonstrated that it could produce the goal functional with a greater accuracy for a given number of degrees of freedom in the majority of cases. As with the regular adaptivity, the performance of the goal-based method decreased on the eigenvalue problem tested. The accuracy of the functional obtained with the goal-based adaptivity was better than that of both the uniform resolution and the regular adapted resolution for the fixed source examples. In contrast, the goal-based adaptivity using the eigenvalue as the goal provided little to no benefit compared to the regular adaptivity.

The energy dependence of the angular resolution has proven to be beneficial in the multigroup test cases. The use of energy dependent angular resolution does not require any extra computational effort. The added cost is that the distribution of the angular resolution over space must be stored in memory for each group. However, this is small compared to the memory required for storing the operator matrix and solution vectors. The added benefit over the minor cost demonstrates that the energy dependent angular resolution is a worth-while development. For future problems with a larger number of energy groups it may prove prudent to use the same angular resolution for energy groups with similar flux distributions.

There was no significant reduction in CPU time as a result of the adaptive methods. The performance of the adaptive methods with respect to computational time were problem dependent. The regular adaptivity has demonstrated it can provide a smaller error for a given time but not by a significant margin. It is noted that the goal-based methods will only show better performance if they significantly reduce the number of degrees of freedom because of the added adjoint computation. A significant reduction in degrees of freedom is possible for localised functionals. Therefore, one must intuitively consider the goal of the calculation and determine whether large gains could be obtained with goal-based adaptivity. The power in adaptive methods is not solely obtained through a shorter CPU time. The adaptive method can be considered to perform a convergence test for a given problem and do so in a reasonable time by placing the resolution in correct locations. Typical convergence tests would require the user to manually refine and re-run problems which potentially requires more human-hours than that of the adaptive calculation. Therefore, the lacking performance with respect to computational time should not deter from using such methods.

The variable order spherical harmonic scheme used in this work was limited to isotropic angular resolution. A region of space which is highly directional and thus requires a high angular resolution will most likely have an angular flux directed in either one or a few specific directions. The variable order spherical harmonic scheme applies the angular resolution to the whole of the angular sphere of directions, parts of which may be unnecessary to describe the shape of the angular distribution. An angular discretisation that lends itself very well to anisotropic angular resolution is a wavelets expansion in the angular dimension. The hierarchical nature of the wavelets allows the resolution of specific patches on the angular domain at each point in space and energy. Regular adaptivity has been demonstrated previously using wavelets but the use of goal-based adaptivity has not been tested. This would potentially increase the efficiency of the calculation through resolving only the directions required for the specific functional of interest. This in combination with a spatial adaptivity scheme would provide means to obtain the optimal space-angle discretisation for a given model.

igure 1 :

 1 An illustration showing the characteristic directions for P1 and P3 one-dimensional discretisations.

Figure 2 :

 2 Figure 2: The scalar flux solutions for the simple one-dimensional problem used to illustrate the requirement for the moment coupling modification.

Figure 3 :

 3 Figure 3: The one dimensional domain for problem 1. The dimensions are in units of centimetres.

  Angular expansion size.

Figure 4 :

 4 Figure 4: (a) The reference scalar flux solution for problem 1. (b) An example of the resulting angular expansion size from regular adaptivity over the domain of problem 1.

Figure 5 :

 5 Figure 5: The convergence plot for the error norm of the scalar flux for problem 1. The vertical dashed blue line marks the number of CDOF required for the reference P 21 solution.

Figure 6 :

 6 Figure 6: Error as a function of CPU time for regular and goal-based adaptive methods compared to uniform refinement for problem 1.

Figure 7 :

 7 Figure 7: The domain of problem 1, the green region illustrates the location of the goal.

Figure 8 :

 8 Figure 8: The convergence of the goal functional for problem 1. The two plots show results from using various combinations of the forward and adjoint error measures. The vertical dashed blue line marks the number of CDOF required for the reference P 21 solution.

Figure 9 :

 9 Figure 9: The three dimensional domain for problem 2. The dimensions are in units of centimetres.

Figure 10 :

 10 Figure 10: The reference P 11 solution for problem 2.

Figure 11 :

 11 Figure 11: The convergence plot for the error norm of the scalar flux for problem 2. The vertical dashed blue line marks the number of CDOF required for the reference P 11 solution.

Figure 12 :

 12 Figure 12: Error as a function of CPU time for regular and goal-based adaptive methods compared to uniform refinement for problem 2.

Figure 13 :

 13 Figure 13: The three dimensional domain for problem 2. The region used for the goal-based calculation is shown in green.

Figure 14 :

 14 Figure 14: The convergence of the goal functional for problem 2. The two plots show results from using various combinations of the forward and adjoint error measures. The vertical dashed blue line marks the number of CDOF required for the reference P 11 solution.

Figure 15 :

 15 Figure 15: The two dimensional domain for problem 3. The dimensions are in units of centimetres. The region used for the goal-based calculations is shown in green.

Table 3 :

 3 The material properties for the regions in problem 3. Cross sections in units of cm -1 and source in units of cm -3 s -1 .(a) Group 1.(b) Group 2.

Figure 16 :

 16 Figure 16: The reference scalar flux solutions for problem 3.

Figure 17 :

 17 Figure 17: The convergence plot for the error norm of the scalar flux for problem 3. The vertical dashed blue line marks the number of CDOF required for the reference P 11 solution.

  (a) Energy-independent (115,777 CDOF). (b) Energy-dependent -group 1 (91,516 CDOF). (c) Energy-dependent -group 2 (145,371 CDOF).

Figure 18 :

 18 Figure 18: The resulting angular expansion size from regular energy-independent and energydependent adaptivity for problem 3.

Figure 19 :

 19 Figure 19: The convergence of the goal 1 and goal 2 functionals for problem 3. The vertical dashed blue line marks the number of CDOF required for the reference P 11 solution.

  (a) Energy-independent (84,466 CDOF). (b) Energy-dependent -group 1 (74,207 CDOF). (c) Energy-dependent -group 2 (83,294 CDOF).

Figure 20 :

 20 Figure 20: The resulting angular expansion size from goal-based energy-independent and energy-dependent adaptivity for problem 3 targetting goal 1.

  (a) Energy-independent (63,408 CDOF). (b) Energy-dependent -group 1 (72,389 CDOF). (c) Energy-dependent -group 2 (7,059 CDOF).

Figure 21 :

 21 Figure 21: The resulting angular expansion size from goal-based energy-independent and energy-dependent adaptivity for problem 3 targetting goal 2.

3 IIFigure 22 :

 322 Figure 22: The two dimensional domain for problem 4. The dimensions are in units of centimetres.

Figure 23 :

 23 Figure 23: The reference scalar flux solutions for problem 4.

Figure 24 :

 24 Figure 24: The convergence plots for the regular and goal-based calculations for problem 4. The vertical dashed blue line marks the number of CDOF required for the reference P 11 solution.

  (a) Energy-independent (209,786 CDOF). (b) Energy-dependent -group 1 (292,780 CDOF). (c) Energy-dependent -group 2 (155,034 CDOF).

Figure 25 :

 25 Figure 25: The resulting angular expansion size from regular energy-independent and energydependent adaptivity for problem 4.

Figure 26 :

 26 Figure 26: The resulting angular expansion size from goal-based energy-independent and energy-dependent adaptivity for problem 4.

Figure 27 :

 27 Figure 27: The two dimensional domain for problem 5. The dimensions are in units of centimetres. The dotted line marks the eighth of the geometry which shares the orientation as the other plots in this section. Reflective boundary conditions are applied to all external edges.

Figure 28 :

 28 Figure 28: The reference P 31 scalar solution for problem 5.

Figure 29 :

 29 Figure 29: Azimuthal flux plots at a polar cosine of 0.33 for problem 5. The locations of the nodes are shown in figure 28.

Figure 30 :

 30 Figure 30: An example of the resulting angular expansion size from regular adaptivity over the domain of problem 5.

  (b) Error vs time.

Figure 31 :

 31 Figure 31: The error norm of the scalar flux as a function of CDOF and CPU time for problem 5.

  The algorithm used to construct the fine scale angular discretisation. The symbols Υ φ e and Υ θ e represent the set of nodes in element e for the coarse and fine scale meshes respectively.

	6. Spherical Harmonics Coupling
	0.75
						P1 Solution
						P3 Solution
	0.7 Set initial angular discretisation.				Full moment coupling Modified coupling
	for a = 1 to N a do
	Solve equation(s) for particle flux solution(s). 0.65
	Calculate continuous nodal error field(s) for each energy group. [Loop over each energy group] for g = 1 to G do 0.6 Scalar flux
	[Loop over each continuous spatial node]
	for i = 1 to η N do 0.55
	if e g,i < 0.1 then
	N φ g,i → max{N φ g,i -2, 1} 0.5 P3		[Decrease expansion order by 2 ] P1 P3
	else if 1.0 < e g,i < 2.0 then
	N φ g,i → min{N φ 0.45 0 0.2 g,i + 2, N end if	0.4	Distance (cm)	0.6	0.8	1
	end for
	Apply smoothing algorithm to N φ g
	end for
	end for
	Algorithm 2 [Loop over each element]
	for e = 1 to η do
	[Loop over the nodes of the element]
	for i ∈ Υ θ e do N θ i = max e j∈Υ φ	(N φ j )
	end for
	end for

max } [Increase expansion order by 2 ] else if e g,i > 2.0 then N φ g,i → min{N φ g,i + 4, N max } [Increase expansion order by 4 ]

Table 2 :

 2 The material properties for the different regions in problem 2. Region colours refer to figure

Table 4 :

 4 The material properties for the regions in problem 4. Cross sections in units of cm -1 .

	Region	Σa	Σt	Σ s,g→1	Σ s,g→2	νΣ f	χ
	I
	Group 1	2.882571E-02	4.231071E-01	3.935887E-01	6.927316E-04	1.947030E-02	1.0
	Group 2	2.081569E-01	6.100180E-01	1.491830E-04	4.017119E-01	3.756630E-01	0.0
	II
	Group 1	1.292164E-03	2.801485E-01	2.786851E-01	1.712806E-04	-	-
	Group 2	3.304970E-03	2.478538E-01	2.850373E-05	2.445203E-01	-	-
	III
	Group 1	2.723034E-04	6.174202E-01	5.860507E-01	3.109721E-02	-	-
	Group 2	1.343764E-02	1.648910E+00 3.389150E-05	1.635438E+00 -	-

Table 5 :

 5 The material properties for the different regions in problem 5. Region numbers refer to those in figure

		Σ a	Σ t	Σ s	Source
	I	0.083524 0.141367 0.057843 2π
	II	0.064132 0.072774 0.008642 0.0