Optimised Adjoint Sensitivity Analysis Using Adjoint Guided Mesh Adaptivity Applied to Neutron Detector Response Calculations

This article presents a new approach for the efficient calculation of sensitivities in radiation dose estimates, subject to imprecisely known nuclear material cross-section data. The method is a combined application of adjoint-based models to perform, simultaneously, both the sensitivity calculation together with optimal adaptive mesh refinement. Adjoint-based sensitivity methods are known for their efficiency since they enable sensitivities of all parameters to be formed through only two solutions to the problem. However, the efficient solutions can also be obtained by their computation on optimal meshes, here guided by goal-based adjoint approaches. It is shown that both mesh adaptivity and sensitivity can be computed by the same adjoint solution, meaning both can be performed without additional costs. A simple demonstration is presented based on the Maynard fixed-source problem where uncertainties, with respect to material cross-sections, of doses received in local regions are examined. It is shown that the method is able to calculate sensitivities with reduced computational costs in terms of memory and potentially computational time through reduced mesh size when using adaptive resolution in comparison to uniform resolution. In particular, it is the local spatial contributions to the sensitivity that are resolved more effectively due to the adaptive meshes concentrating resolution in those areas contributing the most to its value.


Introduction
Sensitivity analysis (SA) [1] has been widely employed across many fields of engineering and serves as an underlying tool in areas, including design optimisation, uncertainty quantification and inverse problems. In the field of nuclear engineering, SA plays an essential role in the safe and efficient operation of reactors, and it is used extensively in the quantification of uncertainties in predicted responses due to the imprecisely known data. Whilst a wide variety of SA methods exist, adjoint approaches have been proven particularly popular within the nuclear modelling community. As championed by the seminal works of [1], first-order response sensitivities can be computed simultaneously for all uncertain parameters using just two solutions, one of the governing equations and one of its associated adjoints. This is particularly beneficial in nuclear engineering as, very often, the number of uncertain parameters is large (potentially in the thousands to tens of thousands) in comparison to the number of system responses being measured. Therefore, obtaining all sensitivities through a single forward and adjoint calculation makes it one of the most efficient methods to employ. Since the early works of Cacuci, the use of adjoints in sensitivity analysis has infiltrated many areas of nuclear engineering and reactor physics. Examples include measuring sensitivities in reaction rates and dose calculations in radiation shielding and sub-critical reactors [2]. They have been developed for calculating sensitivities in k e f f [3][4][5], and further applied in burn-up calculations, thermal hydraulics [6], accident analysis and coupled systems [7]. The importance of higher-order sensitivity measurements for UQ is also being realised, and again, the seminal works in [8][9][10] are leading the way in their efficient computation through adjoint methods.
Whilst adjoint-based sensitivity methods enable the efficient computation of response sensitivities through the requirement of only two modelled solutions, obtaining the solutions is often itself a computational burden. For neutron transport problems, the high dimensionality of the governing equations, i.e., the Boltzmann transport equation (BTE), often resolved over complex geometries of reactors and radiation shields can lead to large numbers of degrees of freedom to solve. Therefore, many works have been applied to allow more efficient solutions to these problems, and one such area is through the adaptive resolution of the phase space dimensions of the transport equation. Examples applying adaptive refinement of the equation's spatial dimension include Adaptive Mesh Refinement (AMR) [11][12][13] for structured meshes and fully unstructured adaptive mesh approaches [14,15]. Other approaches have based adaptive resolution on resolving the angular dimension through hierarchical schemes such as wavelets [16,17] and spherical harmonics [18][19][20] and more recently through optimised reduced-order basis functions [21]. These techniques have been further developed to adapt the combined space and angle domains simultaneously. Examples include the works of [22] that has developed spaceangle adaptivity for the first-order BTE using a space-angle DGFEM. This utilised the adaptive discontinuous FEM for resolving space together with an adaptive DG (Discontinuous Galerkin) angle scheme that employed an array of angular basis functions that resolved local patches of the sphere. Another work adaptively resolving the space-angle dimensions of the first-order BTE is that of [23] presenting an adaptive wavelet scheme for the angular treatment combined with an adaptive continuous FEM space discretisation. An adaptive FEM-spherical harmonic method for resolving the even-parity form of the neutron transport equations has been developed in [24]. In all the mentioned methods, reductions in problem size of many orders of magnitude can be achieved, but they also all rely on a metric to guide mesh resolution. Developing such metrics is its own active field of research, but here, adjoint models have played a role, particularly for goal-based calculations where resolution is guided to resolve particular system responses, including eigenvalue [14,15] and detector-response [25].
The work in this article combines the methods of adjoint-based sensitivity analysis and adjoint-based adaptive mesh refinement to form a new approach for optimally quantifying responses and first-order sensitivities of responses, in reaction rates or detector response calculations. The adjoint models proposed by Cacuci [1] are used to form the first-order sensitivities of responses, and the adjoint methodologies developed in [25] are employed for guiding mesh resolution. Critically, the same system response is used in both methods, meaning the adjoint solution that is used to form the sensitivities is also used to guide mesh resolution. Therefore, only two solutions (forward and adjoint) are still required per adapted mesh; hence, this method is efficient in both reducing the number of solutions required and reducing the computational cost of obtaining them. This is demonstrated here by resolving the Maynard problem and quantifying a detector response over a localised region and calculating the sensitivities with respect to the problem's material cross-sections. Whist it was shown that the method does not necessarily improve accuracy when the sensitivity is calculated with respect to changes in whole material regions, the method does show vast improvement in accuracy when calculating sensitivities to material changes across local regions in space. That is, the adaptive mesh approach can provide more accurate spatially dependent sensitivities for material changes in comparison to using uniform resolution meshes.
The following sections are set out as follows. Section 2 introduces the Boltzmann transport equation, its adjoint form, and the adjoint-based formulation for sensitivity calculations and metrics for mesh refinement. Section 3 presents the numerical example based on the Maynard problem, and Section 4 finishes with a conclusion.

Optimised Adjoint Sensitivity Formulation for Neutron Transport Problems
The transport of neutrons is described by the fixed source, steady-state, multi-group Boltzmann transport equation, which is reads as, This describes the distribution of the group g angular flux φ g (r, Ω) in its three spatial (r) and two angular (Ω) dimension phase-space. The terms account for neutron streaming, absorption, group-to-group scatter, and external neutron sources where Σ g t (r) and Σ g →g,Ω →Ω s (r) denote the total and scatter nuclear macroscopic material cross-sections and S g (r) the isotropic source strength. Boundary conditions accompanying the equations include void and reflective conditions defined as, and respectively. The term n denotes the outward pointing normal to the boundary and Ω * is the specular reflective angle to Ω with respect to n. The adjoint form of the BTE equation reads as, which is expressed in terms of the adjoint angular flux ψ g . The source term S a,g denotes the source for the group g adjoint flux and will be defined in a later section. Boundary conditions for void and reflective conditions are defined as, and ψ(r, Ω) = ψ(r, Ω * ), Ω.n > 0, respectively. The BTE and associated boundary conditions (Equations (1)- (3)) are written in compact form, where L denotes the vector of G operators for each group flux, φ(x) denotes the vector of G group fluxes, Q the vector of G sources, and x and the space-angle phase-space dimensions. The terms B and A denote the G operators imposing the boundary conditions on the angular flux. The term α(x) denotes the set of imprecisely known parameters, which, in this article, will be the spatially dependent material cross-sections Σ g t (r) and Σ g→g s (r, Ω → Ω ), source terms S g (r) and, for the reaction (dose) calculation, the cross-sections Σ g d (r). This article considers real valued response functionals of e = (φ(x), α(x)) in the form of dose calculations defined by, Nominal responses R 0 can be computed by solving Equation (7) to obtain nominal fluxes φ 0 (x), for given nominal values of parameters α 0 (x), and directly inserting them into Equation (9). Response sensitivities with respect to the parameters α(x) can also be computed. As described in [1], one approach is from the Gateau derivatives (G-derivatives), δR(e 0 , h), defined as, for variation h = (h φ , h α ) around e 0 , and ∈ F (space of real numbers). Provided the response is linear in h about the local neighbourhood of e 0 , the variation in the response can be written as, where R φ and R α(x) denote G-derivative of the response with respect to φ(x) and α(x), respectively. The first term in Equation (11) requires the already known variation in the parameter space and thus can be computed immediately. The last term (the indirect term) requires the variation in angular flux due to the variation in the parameter space h α , and this is not immediately known. To compute h φ in relation to h α , the Gateau differential of the governing Equation (7) are formed, where the L α and B α denote the G-differentials of L and B with respect to α. These can be solved for each variation in h α to obtain the variation h φ , which can then be inserted into Equation (11) to form the variation in the response. As discussed in [1], each different variation in the parameter space requires a separate solution to Equations (12) and (13), meaning the approach is inefficient when large numbers of sensitivities are sought.

Summary of First-Order Adjoint Sensitivity Analysis
Adjoint methods [1] have previously been developed to efficiently compute firstorder response sensitivities. Computing the second, indirect, term of Equation (11) is accomplished through the adjoint relationship, The term L † denotes the system of G adjoint operators in the form of the left-hand side of Equation (4) and ψ the associated set of G adjoint functions. The termP is an additional bi-linear term that operates on the boundary to ensure the operators L and L † are formally adjoint. In this work, the enforcement of the boundary conditions (5) and (6) on the adjoint equations ensures this term vanishes. It will, therefore, not be considered from here on.
The indirect term of Equation (11) is re-written using the inner product form, and by imposing the relation, Equations (14) and (15) can be combined with (12) and (13) to give, In this form, a single adjoint Equation (16) needs to be solved for each response under consideration. The adjoint solution is inserted into Equation (17), and now only the inner product of Equation (17) needs to be computed for each sensitivity sought. In comparison to a full model solve, Equations (12) and (13), evaluating the inner product, Equation (17), is computationally efficient.

Adjoint Error Metrics for Adaptive Mesh Refinement
Spatial mesh adaptivity relies upon generating a sequence of FEM meshes with resolution adapted to optimally resolve a goal of interest. In this work, the goal is the same dose response given in Equation (9), to which sensitivities, with respect to the cross-section parameters, are also sought. The adaptive method is initialised with a mesh upon which a solution to Equation (1) is formed, and the response is calculated by inserting the solution into Equation (9). The mesh is then adapted to reduce errors in the response calculation by increasing the mesh resolution in the regions contributing most to the errors. To guide the adapting method, an appropriate metric is formed based on adjoint methods.
A metric is formed by considering the component within the integral of the response in Equation (9), Given some inexact solution φ, a truncated Taylor's expansion from the true solution, φ ex , is given as, Similarly, a residual to the governing Equation (7) can be formed, and expanded through a truncated Taylors expansion, with respect to φ, to give, Combining (19) and (21), the error in the response is given as, Integration over the phase-space dimensions gives the response error where the system, is solved for the adjoint solution ψ. Importantly, the right-hand side of Equations (16) and (24) are identical. This means that the same, single, adjoint solution is used to form the error metric for adaptivity and to perform the sensitivity calculation.

Adaptive Mesh Refinement
Given a solution obtained on some initial given mesh, an improved solution is obtained by using an adapted mesh that reduces the error functional given in Equation (23). This process is then repeated until mesh convergence is reached and the approximate solution is within a desired accuracy. The details of this procedure are quite extensive; thus, only a light overview of the method is provided, and guidance to the relevant literature is instead given. The adaptive procedure was originally developed in [26], which, if needed, allows the formation of anisotropic elements that can refine resolution to the appropriate regions and directions, as guided by the metric. The metric itself in Equation (23) requires measures of the solution errors and the residual of the adjoint solution. Obtaining these from the numerical solutions of the forward and adjoint fluxes requires careful consideration. In this article, the errors are formed using Hessian-based approximations, details of which are contained in [15,25], which provide the information on the location of errors and their orientation (i.e., which direction the errors are propagating through). These same articles also discuss how the residuals are approximated from the modelled solutions. The article [26] provides the details of the mesh adapting procedure given the metrics on where resolution needs to be placed.
Upon each adapted mesh, the solutions to the forward equations are re-computed to form improved estimates of the response, and this is followed by solving the adjoint equations for computing the adaptive mesh metrics together with the new estimates of its sensitivity. Therefore, for each adaptation of the mesh, two solutions are required.

Numerical Example-Maynard Problem
This section presents a numerical example applying the adjoint-based sensitivity models with adjoint-guided mesh adaptivity to a mono-energetic source-detector problem. The problem employed is the 2D Maynard benchmark presented in Figure 1, which consists of a neutron source (region 1) surrounded by a void (region 2), all encompassed within a highly scattering material (region 3). Table 1 lists the material cross-sections and source strengths. In this demonstration, a response region is located close to the exterior of the domain, as indicated in Figure 1. The objective is to calculate the neutron flux within this region, and so the response of Equation (9) is defined as, (25) where v r denotes the integration restricted to the response region. The reaction crosssection is defined as Σ g d (r) = 1.0 inside the response region and zero elsewhere. The spatially dependent sensitivity of this response subject to the (assumed to be) imprecisely known material total cross-section of region 3 will be the focus of this numerical example. In all the following analyses, the finite element code FETCH [23] is employed to resolve both the forward and adjoint equations, and the solutions are used to calculate the response and its sensitivity. As the focus of this article is centred on the spatial dependence of the response sensitivity, a high order S 30 angular discretisation is employed to ensure sufficient resolution is applied in the angular dimension. A forward and adjoint reference solution was obtained using a uniform, high-resolution spatial mesh consisting of approximately 420 K triangular elements. Figure 2 presents the forward and adjoint reference solution's scalar flux profiles. Included in this figure is the profile showing the spatial variation in response sensitivity with respect to change in the macroscopic total cross-section, Σ t (r), for the material occupying region 3, denoted by v 3 . It shows at each position, r is the component within the integral of Equation (17), essentially the multiplication of the flux and adjoint solutions, Here, δ r,r 3 is the function with a value of 1 if r ∈ v 3 and value of 0 otherwise. This sensitivity profile essentially illustrates the change in the contribution to the overall sensitivity at every spatial position and highlights its varying nature. Whilst it is intuitive that the largest sensitivities form in those regions in a direct line between the source and response region, other areas are also shown to have sizable contributions, for example, the top interface between regions 2 and 3 and the wider, surrounding regions directly between the source and response locations. Two adapted mesh scalar flux solutions and sensitivity profiles are presented in Figure 3 using approximately 6 and 26 K elements. For both adapted meshes, it is evident that the mesh resolution is focused on the aforementioned regions identified as having large contributions to the sensitivity of the response. That is, a high concentration of elements, and thus the mesh resolution, is centred on and around the regions in a direct line between the source and response region. In addition, higher resolution is also concentrated around the top interface between regions 2 and 3. These areas are not in the source-response region's line of sight but do have large contributions to the response and its sensitivity.
A study was conducted to investigate the use of uniform and adapted meshes to accurately predict the response and the sensitivity with respect to Σ t within region 3. Calculations were performed on uniform meshes ranging in size from 1 to 420 thousand elements. The highest order mesh was considered to be the reference solutions, from which the response and its sensitivity with respect to Σ t inside region 3 were R = 0.03061677 and ∂R ∂Σ t = −0.40875, respectively. Four adaptive mesh calculations were performed where the adapted meshes were allowed to evolve from 1.5 and 40 thousand elements. Figure 4 presents the graph showing the response as a function of mesh size. The adaptive mesh calculations are shown to converge to within 10 −4 of the reference solution using approximately 10 thousand elements; in comparison, the uniform mesh required around 100 thousand elements to achieve the same accuracy. Thus, as seen in previous studies, mesh adaptivity can result in the reduction in response errors using meshes significantly smaller in size. Figure 4 also presents the relative error in the response sensitivity using both uniform and adaptive meshes. Interestingly, for the low-order meshes, the accuracy appears around the same order, but for 10 K+ elements, the uniform mesh performs better. However, one must take into account that the integral quantity of the sensitivity is subject to the cancellation of errors, where regions of positive sensitivity will cancel with regions of negative sensitivity. The oscillations in the uniform mesh convergence also indicate that the solutions are not converged. Therefore, to provide more reliable estimates of accuracy, spatially dependent contributions to the total sensitivity is a more appropriate measure.   To quantify the error in the spatial dependence of the response sensitivity, the sensitivity is measured along the two lines x = 6 cm and y = 2.5 cm, as indicated in Figure 1. Figure 5 presents the sensitivity profiles using adapted and uniform meshes consisting of approximately 6 thousand elements. They show that the adapted meshes are able to resolve many of the unphysical oscillations that formed in the uniform mesh solution. Thus, there is a clear benefit in using adapted meshes when relatively small meshes are necessary. The graph presented in Figure 6 shows the 2-norm errors measured along the two lines investigated when using uniform and adapted meshes. These clearly show that a large reduction in errors is obtained when using the adapted mesh approach. The result for line 1 shows a faster order of convergence with adapted meshes, which is desirable as it is measured across a region, contributing heavily to the total sensitivity. For the second line, the reduction in error using adaptivity is about a factor of 1.5 for smaller meshes, but this accelerates when using mesh sizes of 10 thousand elements or more.

Conclusions
This article has presented the first method of using adjoint-based sensitivity models combined with adjoint-based, anisotropic mesh adaptivity. The key element of the formulation is that the same adjoint solution to the problem can be used to form the sensitivity as well as guide the mesh adaptivity. Therefore, only two large-scale FEM solutions (i.e., the forward and adjoint solutions) are required to perform both operations at each adaption of the mesh.
The short demonstration presented here shows an example of this through an adaption of the Maynard benchmark into a source-detector calculation. The study of the response sensitivity with respect to the total cross-section showed that the contribution to the response and its sensitivity varied substantially throughout the spatial domain. However, as the regions important to the response and sensitivity calculations were localised to small parts of the problem domain, the adaptive mesh could reduce errors in both calculations using substantially smaller meshes-in this case, by nearly an order of magnitude. Whilst it was difficult to determine the benefits of calculating sensitivities to changes in the material properties across entire regions as a whole, it has been demonstrated that local spatial contributions are resolved more accurately. This is quite important as this avoids potentially affecting the modelled predictions through the cancellation of errors.
The efficiency of using adjoints reduces when the number of responses to be calculated is large and number of uncertain parameters is small. Therefore, the method would be most effective for those problems requiring few response calculations but involving a large number of unknown parameters (e.g., cross-sections, number densities, neutron sources) and where only localised regions are important to resolve. Such examples would include shielding-type problems involving source-detector responses, perhaps where the response is the leakage through the shield's surface. In fact, we have already demonstrated in our previous works, cited in this article, that shields with ducts allowing high streaming are particularly well suited to adaptive mesh refinement. Here, the resolution typically centres on the ducts and their near-wall regions, but also, these same regions are likely to contribute to the sensitivity of the leakage. In fact, a similar trait was shown here where the interfaces between the void-scattering region were resolved at a higher resolution in comparison to other parts of the scattering region. Extending to shielding problems in three dimensions and/or those with smaller regions of importance (in comparison to problem size) is likely to increase the efficiency gains further. However, it is this reduction in mesh size that improves efficiency in solver times and memory requirements, combined with the need for only two full solutions, per adaptive mesh, that provide the overall benefits of the technique for some types of radiation plus sensitivity calculations.