Back
 AJCM  Vol.11 No.2 , June 2021
Fourth-Order Adjoint Sensitivity Analysis of an OECD/NEA Reactor Physics Benchmark: I. Mathematical Expressions and CPU-Time Comparisons for Computing 1st-, 2nd- and 3rd-Order Sensitivities
Abstract: This work extends to fourth-order previously published work on developing the adjoint sensitivity and uncertainty analysis of the numerical model of a polyethylene-reflected plutonium (acronym: PERP) OECD/NEA reactor physics benchmark. The PERP benchmark comprises 7477 imprecisely known (uncertain) model parameters which have nonzero values. These parameters are as follows: 180 microscopic total cross sections; 7101 microscopic scattering sections; 60 microscopic fission cross sections; 60 parameters that characterize the average number of neutrons per fission; 60 parameters that characterize the fission spectrum; 10 parameters that characterize the fission source; and 6 parameters that characterize the isotope number densities. Previous works have used the adjoint sensitivity analysis methodology to compute exactly and efficiently all of the 7477 first-order and 27,956,503 second-order sensitivities of the PERP benchmark’s leakage response to all of the benchmark’s uncertain parameters. These works showed that largest response sensitivities involve the total microscopic cross sections, which motivated the recent computation of all of the (180)3 third-order sensitivities of the PERP leakage response with respect to these total microscopic cross sections. It turned out that some of these 3rd-order cross sections were far larger than the corresponding 2nd-order ones, thereby having the largest impact on the uncertainties induced in the PERP benchmark’s response. This finding has motivated the development of the original 4th-order formulas presented in this work, which are valid not only for the PERP benchmark but can also be used for computing the 4th-order sensitivities of response of any nuclear system involving fissionable material and internal or external neutron sources. Subsequent works will use the adjoint-based mathematical expressions obtained in this work to compute exactly and efficiently the numerical values of the largest fourth-order sensitivities of the PERP benchmark’s response to the total microscopic cross section and use them for a pioneering fourth-order uncertainty analysis of the PERP benchmark’s response.

1. Introduction

Until recently, only the first-order sensitivities (i.e., functional derivatives) of a computational model’s responses (i.e., quantities of interest) to the respective model’s imprecisely known (i.e., uncertain) parameters have been considered when assessing the uncertainties induced in the respective responses by the parameter uncertainties. The second- and higher-order sensitivities could not be computed, except for very simple models comprising a handful of parameters, so these sensitivities were ignored. The Second-Order Adjoint Sensitivity Analysis Methodology (2nd-ASAM) recently conceived by Cacuci [1] is the only practical method that enables the exact computation of the large number of 2nd-order sensitivities arising in large-scale problems comprising many parameters. The application of the 2nd-ASAM to a multiplying nuclear system with source [2] [3] [4] has opened the way for the large-scale application presented in [5] - [10] to a polyethylene-reflected plutonium (acronym: PERP) OECD/NEA reactor physics benchmark [11]. The numerical model of the PERP benchmark comprises 21,976 uncertain parameters, of which 7477 parameters have nonzero values, as follows: 180 group-averaged total microscopic cross sections, 7101 non-zero group-averaged scattering microscopic cross sections (the other scattering cross sections, of which there are 21,600 in total, are zero); 120 fission process parameters; 60 fission spectrum parameters; 10 parameters describing the experiment’s nuclear sources; and 6 isotopic number densities.

All of the non-zero first-order sensitivities and second-order sensitivities of the PERP leakage response with respect to the benchmark’s parameters were computed, ranked, and analyzed in [5] - [10]. The results obtained in [5] - [10] showed that the 2nd-order sensitivities of the leakage response with respect to the group-averaged microscopic total cross sections are the largest (by comparison to the 2nd-order sensitivities of the leakage response with respect to the other uncertain model parameters) and have therefore the largest impact on the uncertainties induced in the leakage response. For example, neglecting the 2nd-order sensitivities of the leakage response with respect to the total cross sections could cause an error as large as 2000% in the expected value of the leakage response and up to 6000% in the variance of the leakage response [5]. Given that the effects of these 2nd-order sensitivities are much larger than the effects of the 1st-order sensitivities [5] [10], it is logical to quantify the magnitudes and contributions that would stem from the 3rd-order sensitivities of the PERP benchmark’s total leakage response with respect to the microscopic total cross sections. To enable the computation of such 3rd-order sensitivities, Cacuci [12] has recently conceived the “third order adjoint sensitivity analysis methodology for reaction rate responses in a multiplying nuclear system with source” and then applied this general theory to the PERP benchmark in order to derive the exact analytical expressions of the 3rd-order sensitivities of the PERP benchmark’s leakage response with respect to this benchmark’s microscopic total cross sections [13] [14] [15]. These works [13] [14] [15] have shown that the largest 3rd-order sensitivities involve the microscopic total cross section for the lowest (30th) energy group of isotope 1H (i.e., σ t , 6 30 ). In particular, the largest overall 3rd-order sensitivity is the mixed 3rd-order sensitivity S ( 3 ) ( σ t , 1 g = 30 , σ t , 6 g = 30 , σ t , 6 g = 30 ) = 1.88 × 10 5 , which also involves the microscopic total cross section for the 30th energy group of isotope 239Pu (i.e., σ t , 1 g = 30 ). The largest unmixed 3rd-order sensitivity is also with respect to σ t , 6 30 , namely S ( 3 ) ( σ t , 6 g = 30 , σ t , 6 g = 30 , σ t , 6 g = 30 ) = 2.966 × 10 4 . These sensitivities significantly larger than the largest lower-order sensitivities; by comparison, the largest 1st-, 2nd- and 3rd-order sensitivities are S ( 1 ) ( σ t , 6 30 ) = 9.366 , S ( 2 ) ( σ t , 6 30 , σ t , 6 30 ) = 429.6 and S ( 3 ) ( σ t , 1 g = 30 , σ t , 6 g = 30 , σ t , 6 g = 30 ) = 1.88 × 10 5 , respectively. The results obtained in [5] - [10] and [13] [14] [15] indicated that the total microscopic cross section of isotopes 1H and 239Pu are the most important parameters affecting the PERP benchmark’s leakage response, since they are responsible for the largest 1st-, 2nd- and 3rd-order sensitivities.

This work is organized as follows: Section 2 presents the computational modeling of the PERP Benchmark. Section 3 presents the main features and the CPU-times required for applying the three fundamental deterministic (as opposed to “statistical”) methods for computing the 1st-order response sensitivities, namely 1) finite differences; 2) the forward sensitivity analysis method (FSAM); and 3) the 1st-Order Comprehensive Adjoint Sensitivity Analysis Methodology (1st-CASAM) for linear systems. Section 4 presents the main features and the CPU-times required for applying the three fundamental deterministic methods (i.e., finite differences, FSAM the 2nd-CASAM, respectively) for computing the 2nd-order response sensitivities. The results presented in Section 4 clearly highlight the fact that the 2nd-CASAM is the only method that can be used in practice to compute 1st- and 2nd-order sensitivities of a response of a large-scale system, involving many parameters. Section 5 presents comparisons of the CPU times required for computing the 3rd-order sensitivities of the PERP benchmark’s leakage response solely to the group-averaged total cross sections. These results underscore the fact that the 3rd-CASAM is exact, introducing no intrinsic methodological errors in the computation of sensitivities, and is by far more efficient than any other method. Section 6 concludes this work, noting that the expressions and results for the 4th-order sensitivities will be presented in companion works [16] [17]. Notably, the original mathematical concepts underlying the 1st-order FSAM and 1st-CASAM were presented in [18] [19] [20] while the general theory underlying the 4th-Order Comprehensive Adjoint Sensitivity Analysis Methodology (4th-CASAM) for linear systems is presented in [21].

2. Mathematical/Computational Model of the PERP Benchmark

The spherical polyethylene-reflected plutonium (acronym: PERP) benchmark comprises a metallic inner sphere (“core”) containing the following 4 isotopes: Isotope 1 (239Pu), Isotope 2 (240Pu), Isotope 3 (69Ga) and Isotope 4 (71Ga). This core (which is designated as “material 1”) is surrounded by a spherical shell of polyethylene (designated as “material 2”), containing two isotopes, designated as Isotope 5 (C) and Isotope 6 (1H), respectively. The dimensions and material composition of the polyethylene-reflected plutonium (PERP) metal sphere considered in this work are presented in TableA1 in the Appendix. The quantity of interest in this work, which will be called the “response,” is the leakage of neutrons out of the PERP sphere, which has been measured experimentally [11].

The numerical modeling of the neutron flux distribution within the PERP benchmark, as well as the computation of the PERP leakage response has been performed using the multigroup discrete ordinates particle transport code PARTISN [22] together with neutron sources computed using the code SOURCES4C [23]. The neutron flux distribution within the PERP benchmark, as well as the leakage of neutrons out of the PERP sphere has been modeled using the standard multigroup form of the Boltzmann neutron transport equation subject to the boundary condition of no incoming flux, with an internal spontaneous fission source, which can be written in the following form:

Ω φ g ( r , Ω ) + Σ t g ( α ; r ) φ g ( r , Ω ) h = 1 G 4 π d Ω φ h ( r , Ω ) [ Σ s h g ( α ; r , Ω Ω ) + χ g ( α ; r ) ( ν Σ f ) h ( α ; r ) ] = Q g ( α ) , g = 1 , , G , (1)

φ g ( r , Ω ) = 0 , r = r d , Ω n < 0 , g = 1 , , G , (2)

where r d is the radius of the PERP sphere. Note that many works, use the “ g ” as the “dummy” summation index in Equation (1); the index h (instead of g ) is preferred in this work in order to avoid, as much as possible, the appearance of “superscripted superscripts.” In this work, therefore, the superscript-indices “g” and “h” will refer exclusively to “energy groups.”

The source, Q g ( α ) , which appears on the right-side of Equation (1) is defined as follows:

Q g ( α ) k = 1 N f λ k N k F k S F ν k S F ( 2 π a k 3 b k e a k b k 4 ) E g + 1 E g d E e E / a k sinh b k E , (3)

where N f denotes the total number of spontaneous-fission isotopes. The PERP benchmark has N f = 2 spontaneous-fission isotopes, namely “isotope 1” (239Pu) and “isotope 2” (240Pu). The spontaneous fission neutron spectrum of 239Pu and 240Pu, respectively, is each approximated by a Watts fission spectrum, characterized by: 1) the evaluated parameters a k and b k ; 2) the spontaneous emission, ν k S F , of an average neutron of an isotope k; 3) the decay constant, λ k , for actinide nuclide k; 4) the fraction, F k S F , of decays that are spontaneous fission (the “spontaneous-fission branching fraction”); 5) the index k = 1 refers to “isotope 1” (239Pu) while k = 2 refers to “isotope 2” (240Pu). The symbols used in Equations (1) through (3) have their usual meanings; for convenient referencing, they are described in the Appendix and are summarized in the Nomenclature.

The mathematical expression of the PERP benchmark’s leakage response, denoted as L ( α ) , is provided below:

L ( α ) S b d S g = 1 G Ω n > 0 d Ω Ω n φ g ( r , Ω ) . (4)

The PARTISN [22] computations of the neutron flux used the MENDF71X [24] 618-group cross section data collapsed to G = 30 energy groups, as well as a P3 Legendre expansion of the scattering cross section and a fine-mesh spacing of 0.005 cm (comprising 759 meshes for the plutonium sphere of radius of 3.794 cm, and 762 meshes for the polyethylene shell of thickness of 3.81 cm). The first- and second-order response sensitivities were computed using an angular quadrature of S256. The 3rd- and 4th-order sensitivities of the leakage response with respect to the total cross sections were computed using an angular quadrature of S32. The group boundaries of the G = 30 energy groups are provided in TableA2 in the Appendix. The scattering and fission terms in Equation (1) contain implicitly a factor 1 / 4 π , to conform to the convention used in PARTISN [22].

The vector α , which appears in the arguments of the various quantities in Equations (1) and (4), is defined as follows:

α [ α 1 , , α T P ] [ q ; N ; σ s ; σ t ; σ f ; ν ; p ] , where T P J Q + I + J S X + J T X + J F X + J N U + J χ . (5)

The components of the TP-dimensional vector α are the imprecisely known (i.e., uncertain) model and response parameters, and are described in the Appendix. In Equation (5) and throughout this work, the dagger will be used to denote “transposition.” The nominal values of the parameters in Equation (5), as well as of all other quantities in this work, will be denoted using the superscript “zero,” e.g., α 0 [ α 1 0 , , α T P 0 ] .

For the mathematical derivations to follow in this work, it will be convenient to use the matrix-form of Equations (1) and (2), which is as follows:

B ( α ) φ ( r , Ω ) = Q ( α ) , (6)

φ ( r , Ω ) = 0 , r = r d , Ω n < 0 , (7)

where

B ( α ) ( B 11 ( α ) · B 1 g ( α ) · B 1 G ( α ) · · · · · B g 1 ( α ) · B g g ( α ) · B g G ( α ) · · · · · B G 1 ( α ) · B G g ( α ) · B G G ( α ) ) ; φ ( r , Ω ) ( φ 1 · φ g · φ G ) ; Q ( α ) ( Q 1 · Q g · Q G ) ; (8)

with components defined below:

B g h ( α ) δ g , h B g g 0 ( α ) B g h 1 ( α ) ; g , h = 1 , , G ; B g g 0 ( α ) φ g ( r , Ω ) [ Ω + Σ t g ( α ; r ) ] φ g ( r , Ω ) ; B g h 1 ( α ) φ h ( r , Ω ) 4 π d Ω [ Σ s h g ( α ; r , Ω Ω ) + χ g ( α ; r ) ( ν Σ f ) h ( α ; r ) ] . (9)

The notation δ g , h , which appears in Equation (9), denotes the Kronecker delta-functional, which is defined as usual, i.e., δ g , h = 1 , if g = h and δ g , h = 0 , if g h .

3. Computation of First-Order Sensitivities L ( α ) / α j , j = 1 , , T P

Basically, there are three methods for computing deterministically (as opposed to computing “statistically”) the response sensitivities to model and response parameters, as follows:

1) The so-called “brute-force re-computations” method, which uses finite-difference formulas to approximate the derivative that expresses the response sensitivity to the parameter under consideration; this method will be presented in Subsection 3.1;

2) The Forward Sensitivity Analysis Methodology (FSAM), which will be applied in Subsection 3.2;

3) The 1st-Order Comprehensive Adjoint Sensitivity Analysis Methodology (CASAM) conceived by Cacuci [21], which will be applied in Subsection 3.3.

Subsection 3.4 will present a comparative discussion of the computational resources and times (CPU) which are required by the three above-mentioned deterministic methods.

3.1. Re-Computations with Finite-Difference Approximation

The first-order sensitivities of the leakage response, L ( α ) , can be approximately computed by re-computations using the well-known finite-difference formula presented below:

L ( α ) α j 1 2 h j ( L j + 1 L j 1 ) + Ο ( h j 2 ) , j = 1 , , T P , (10)

where L j + 1 L ( α j + h j ) , L j 1 L ( α j h j ) and h j denotes a “judiciously-chosen” variation in the parameter α j around its nominal value α j 0 . The values L j + 1 L ( α j + h j ) and L j 1 L ( α j h j ) are obtained by re-solving Equations (1) and (2) repeatedly, using the changed parameter values ( α j ± h j ) . As is well known, the value of the variation h j must be chosen by “trial and error” for each parameter α j , and varies from parameter to parameter; when h j is too large or too small, the result produced by Equation (10) will be far off from the exact value of the derivative L ( α ) / α j . It is important to note that finite difference formulas introduce their intrinsic “methodological errors,” such as the error Ο ( h j 2 ) indicated in Equation (10), which are in addition to and independent of the errors that might be incurred in the computation of L ( α ) . In other words, even if the computation of L ( α ) were perfect (error-free), the finite-difference formulas nevertheless introduce their own, intrinsic, numerical errors into the computation of the sensitivity d L ( α ) / d α j .

3.2. Forward Sensitivity Analysis Methodology (FSAM)

If one could solve the neutron transport Equations (1) and (2) in closed form, the respective closed form of φ g ( r , Ω ) would be an explicit function of the parameters α . However, since the neutron transport Equations (1) and (2) cannot be solved in closed form, it follows that φ g ( r , Ω ) is therefore an implicit function of α . Consequently, as Equation (4) indicates, the PERP benchmark’s leakage response, L ( α ) , is a function of the vector α of model parameters implicitly through the group-flux φ g ( r , Ω ) .

The Forward Sensitivity Analysis Methodology (FSAM) is developed [18] [19] in the vector-space of the model parameters α . In this vector-space, the first-order partial sensitivity of L ( α ) to a generic model parameter, α j , can in principle be obtained by taking the first-order partial G-derivative of Equation (4) with respect to a generic model/respond parameter α j , evaluated at the nominal parameter values α 0 , which yields:

{ L ( α ) α j } α 0 d d ε { S b d S g = 1 G Ω n > 0 d Ω Ω n φ g ( α j 0 + ε h j ; α 0 ; r , Ω ) α j } ε = 0 = S b d S g = 1 G Ω n > 0 d Ω Ω n { φ g ( r , Ω ) α j } α 0 , j = 1 , , T P . (11)

In turn, the derivative φ g ( r , Ω ) / α j is the solution of the equations obtained by taking the first-order G-derivatives of Equations (1) and (2) with respect to α j , evaluated at the nominal parameter values α 0 , which yields, by definition, the following equations:

d d ε { Ω φ g ( α j 0 + ε h j ; α 0 ; r , Ω ) + Σ t g ( α j 0 + ε h j ; α 0 ; r ) φ g ( α j 0 + ε h j ; α 0 ; r , Ω ) h = 1 G 4 π d Ω φ h ( α j 0 + ε h j ; α 0 ; r , Ω ) [ Σ s h g ( α j 0 + ε h j ; α 0 ; r , Ω Ω ) + χ g ( α j 0 + ε h j ; α 0 ; r ) ( ν Σ f ) h ( α j 0 + ε h j ; α 0 ; r ) ] Q g ( α j 0 + ε h j ; α 0 ) } ε = 0 = 0 , j = 1 , , T P ; (12)

d d ε { φ g ( α j 0 + ε h j ; α 0 ; r , Ω ) } ε = 0 = 0 , r = r d , Ω n < 0 , g = 1 , , G ; j = 1 , , T P . (13)

Carrying out the differentiations with respect to ε in Equations (12) and (13), setting ε = 0 in the resulting equations and subsequently cancelling everywhere the scalar quantity h j yields the following equations for the function φ g ( r , Ω ) / α j :

{ B g ( α ) [ φ g ( r , Ω ) α j ] } α 0 = { S g [ φ g ( r , Ω ) , r , Ω ] α j } α 0 , g = 1 , , G ; j = 1 , , T P , (14)

{ φ g ( r , Ω ) α j } α 0 = 0 , r = r d ; Ω n < 0 , g = 1 , , G ; j = 1 , , T P , (15)

where

S g [ φ g ( r , Ω ) , r , Ω ] α j Q g ( α ) α j B g ( α ) α j [ φ g ( r , Ω ) ] , (16)

B g ( α ) α j [ φ g ( r , Ω ) ] Σ t g ( α ; r ) α j φ g ( r , Ω ) h = 1 G 4 π d Ω φ h ( r , Ω ) { Σ s h g ( α ; r , Ω Ω ) α j + [ χ g ( α ; r ) ( ν Σ f ) h ( α ; r ) ] α j } . (17)

The system of Equations (14) and (15) is called the 1st-Order Forward Sensitivity System (1st-OFSS), the solutions of which are required in Equation (11) to compute, in the end, the response sensitivities using the forward sensitivity analysis method. Note that a “place-holder,” denoted as “[]”, has been explicitly used on the left-side of Equation (14) to indicate that the operator B g ( α ) acts (linearly) on the function that appears in this place-holder. Place-holders “[]” have also been used on the right-side of Equation (16) and on the left-side of Equation (17) to indicate that the operator B g ( α ) / α j acts on the function that appears in the respective place-holder. Since there are a total of TP parameters α j , it is evident that the computations of all of the sensitivities L ( α ) / α j require TP large-scale computations to solve the 1st-OFSS, cf. Equations (14) and (15), followed by TP small-scale computations for performing the integration represented by Equation (11).

3.3. First-Order Comprehensive Adjoint Sensitivity Analysis Methodology (1st-CASAM)

The aim of the First-Order Comprehensive Adjoint Sensitivity Analysis Methodology (1st-CASAM) is to find an alternative way for expressing the contribution of the flux variation (the so-called “indirect-effect term” contribution) to the total response sensitivity, so as to avoid the need for having to compute the derivatives φ g ( r , Ω ) / α j of the state functions with respect to the model’s parameters, as is the case when using the FSAM. In contradistinction to the FSAM, which is set in the vector-space of the parameters α , the 1st-CASAM is set in the combined space of the parameters and the state-functions, denoted as u ( φ ; α ) = [ φ 1 , , φ G ; α 1 , , α T P ] , and does not a priori consider that the state functions φ [ φ 1 , , φ G ] are implicit functions of the parameters α . In its full generality, the 1st-CASAM also considers that the phase-space domain’s boundaries and internal interfaces are also subject to uncertainties (i.e., are imperfectly well known). For the purpose of this work, however, the boundaries of the phase-space domain (which consist of the boundaries of the spatial domain, i.e., the internal and external radii of the PERP sphere) and the end-points of the energy groups are considered to be free of errors (i.e., are perfectly well known) so that the sensitivities of the leakage response to these phase-space domain boundaries will not be considered. In this case, the 1st-CASAM reduces to the adjoint sensitivity analysis methodology originally conceived by Cacuci [18] [19]. When the response sensitivities with respect to the boundaries of the phase-space domain are not considered, the G-differential, denoted as δ L ( α ; δ φ g ) , of the leakage response defined in Equation (4) for a variation δ u ( δ φ ; δ α ) around the nominal values u 0 ( φ 0 ; α 0 ) is obtained, by definition, as follows:

δ L ( α ; δ φ ) d d ε { S b d S g = 1 G Ω n > 0 d Ω Ω n [ φ g , 0 ( r , Ω ) + ε δ φ g ( r , Ω ) ] } ε = 0 = S b d S g = 1 G Ω n > 0 d Ω Ω n [ δ φ g ( r , Ω ) ] . (18)

In general, the total G-variation of the response depends comprises two contributions, one arising directly from the parameter variations δ α (this contribution is called the “direct-effect term”) and another contribution, stemming from the variations δ φ in the state functions (this contribution is called the “indirect-effect term”). In the particular case of the leakage response considered in this work, the variation δ L ( α ; δ φ g ) does not depend directly on the parameter variations δ α but depends only on variations δ φ in the state functions, as indicated in Equation (18). Therefore, δ L ( α ; δ φ g ) coincides with the “indirect-effect term.”

The variations δ φ g ( r , Ω ) are the solutions of the first-order G-differentials of the model Equations (1) and (2), which are determined by applying the definition of the G-differential to these equations, to obtain:

d d ε { Ω [ φ g , 0 ( r , Ω ) + ε δ φ g ( r , Ω ) ] + Σ t g ( α 0 + ε δ α ; r ) [ φ g , 0 ( r , Ω ) + ε δ φ g ( r , Ω ) ] h = 1 G 4 π d Ω [ φ h , 0 ( r , Ω ) + ε δ φ h ( r , Ω ) ] [ Σ s h g ( α 0 + ε δ α ; r , Ω Ω ) + χ g ( α 0 + ε δ α ; r ) ( ν Σ f ) h ( α 0 + ε δ α ; r ) ] Q g ( α 0 + ε δ α ) } ε = 0 = 0 , (19)

d d ε { [ φ g , 0 ( r , Ω ) + ε δ φ g ( r , Ω ) ] } ε = 0 = 0 , r = r d , Ω n < 0 , g = 1 , , G . (20)

Carrying out the differentiations with respect to ε in Equations (19) and (20), and setting ε = 0 in the resulting equations yields the following system of equations for the variations δ φ g ( r , Ω ) in terms of the parameter variations δ α :

B g ( α 0 ) [ δ φ g ( r , Ω ) ] = Q ( 1 ) ( g ; α 0 , φ 0 ; δ α ) , g = 1 , , G , (21)

δ φ g ( r , Ω ) = 0 , r = r d , Ω n < 0 , g = 1 , , G , (22)

where

Q ( 1 ) ( g ; α 0 , φ 0 ; δ α ) { Q g ( α ) α δ α { B g ( α ) [ φ g ( r , Ω ) ] } α δ α } α 0 , (23)

{ B g ( α ) [ φ g ( r , Ω ) ] } α δ α φ g ( r , Ω ) Σ t g ( α ) α δ α h = 1 G 4 π d Ω φ h ( r , Ω ) { Σ s h g ( α ; Ω Ω ) α δ α + [ χ g ( α ) ( ν Σ f ) h ( α ) ] α δ α } . (24)

By analogy to Equations (6) and (7), Equations (21) and (22) can be written in following matrix-vector form:

B ( α ) [ δ φ ( r , Ω ) ] = Q ( 1 ) ( α ; δ α ) , (25)

δ φ ( r , Ω ) = 0 , r = r d , Ω n < 0 , (26)

where

δ φ ( r , Ω ) ( δ φ 1 ( r , Ω ) · δ φ g ( r , Ω ) · δ φ G ( r , Ω ) ) ; Q ( 1 ) ( α ; δ α ) ( Q ( 1 ) ( 1 ; α , φ ; δ α ) · Q ( 1 ) ( g ; α , φ ; δ α ) · Q ( 1 ) ( G ; α , φ ; δ α ) ) . (27)

The aim of the 1st-CASAM—to construct an alternative expression for the “indirect-effect term” δ L ( α ; δ φ ) which would not depend on δ φ —can be achieved by implementing the following sequence of steps, which are all to be carried out at the nominal parameter values α 0 (although this will not be explicitly indicated in the sequence in order to keep the notation as simple as possible):

1) Consider that the square-integrable functions φ H ( 1 ) and δ φ H ( 1 ) are elements of a Hilbert space denoted as H ( 1 ) , in which the inner product between two elements u ( r , Ω ) H ( 1 ) and v ( r , Ω ) H ( 1 ) is denoted as u ( r , Ω ) , v ( r , Ω ) ( 1 ) and is defined as follows:

u ( r , Ω ) , v ( r , Ω ) ( 1 ) g = 1 G 0 r d 4 π r 2 d r 4 π d Ω u g ( r , Ω ) v g ( r , Ω ) . (28)

2) Form the inner product of Equation (25) with a yet undefined function ψ ( 1 ) ( r , Ω ) [ ψ ( 1 ) , g ( r , Ω ) , , ψ ( 1 ) , G ( r , Ω ) ] H ( 1 ) to obtain the following relation:

ψ ( 1 ) ( r , Ω ) , B ( α ) [ δ φ ( r , Ω ) ] ( 1 ) = ψ ( 1 ) ( r , Ω ) , Q ( 1 ) ( α ; δ α ) ( 1 ) = g = 1 G 0 r d 4 π r 2 d r 4 π d Ω ψ ( 1 ) , g ( r , Ω ) Q ( 1 ) ( g ; α 0 , φ 0 ; δ α ) . (29)

3) Using the inner product defined in Equation (28) and the definition of the adjoint operator in H ( 1 ) , recast the left-side of Equation (29) into the following equivalent form:

ψ ( 1 ) ( r , Ω ) , B ( α ) [ δ φ ( r , Ω ) ] ( 1 ) = δ φ ( r , Ω ) , A ( α ) ψ ( 1 ) ( r , Ω ) ( 1 ) + P [ δ φ , ψ ( 1 ) ] , (30)

where

P [ δ φ , ψ ( 1 ) ] S b d S g = 1 G Ω n < 0 d Ω V | Ω n | [ δ φ g ( r , Ω ) ψ ( 1 ) , g ( r , Ω ) ] S b d S g = 1 G Ω n > 0 d Ω V Ω n [ δ φ g ( r , Ω ) ψ ( 1 ) , g ( r , Ω ) ] , (31)

denotes the bilinear concomitant evaluated on the PERP sphere’s surface S b and where A ( α ) denotes the operator adjoint to B ( α ) , having components A g h ( α ) [ B h g ( α ) ] * , where the symbol [ ] * indicates “formal adjoint operator.”

4) Require the first term on the right-side of Equation (30) to represent the functional defined in Equation (18) to obtain the following relation:

A ( α ) ψ ( 1 ) ( r , Ω ) = I [ Ω n δ ( r r d ) ] ; I [ 1 , , 1 , , 1 ] . (32)

5) Use the boundary condition given in Equation (22) to eliminate the corresponding boundary terms in the bilinear concomitant P [ δ φ , ψ ( 1 ) ] in Equation (31) and eliminate the appearance of the remaining boundary terms in Equation (31) by imposing the following (adjoint) boundary condition:

ψ ( 1 ) ( r , Ω ) = 0 , r = r d , Ω n > 0. (33)

In component form, Equations (32) and (33) are as follows:

Ω ψ ( 1 ) , g ( r , Ω ) + Σ t g ( α ) ψ ( 1 ) , g ( r , Ω ) h = 1 G 4 π d Ω ψ ( 1 ) , h ( r , Ω ) [ Σ s g h ( α ; Ω Ω ) + ( ν Σ f ) g ( α ) χ h ( α ) ] = Ω n δ ( r r d ) , g = 1 , , G ; (34)

ψ ( 1 ) , g ( r , Ω ) = 0 , r = r d , Ω n > 0 , g = 1 , , G . (35)

It is important to note that the relations presented in Equations (28) through (35) hold generally, including at the nominal values u 0 ( φ 0 ; α 0 ) . In order to simplify the notation, the superscript “zero” denoting nominal values has been omitted in the foregoing derivations and will occasionally be omitted henceforth, as well. This simplification should not cause any loss of clarity since it will become clear from the context which quantities are to be evaluated/computed using the nominal values for the model parameters.

6) Use Equations (29) through (33) in Equation (18) to obtain the following expression for the “indirect-effect term”:

δ L ( α ; δ φ ) = ψ ( 1 ) ( r , Ω ) , Q ( 1 ) ( α ; δ α ) ( 1 ) = g = 1 G 0 r d 4 π r 2 d r 4 π d Ω ψ ( 1 ) , g ( r , Ω ) Q ( 1 ) ( g ; α 0 , φ 0 ; δ α ) j = 1 T P L ( α ) α j δ α j . (36)

7) Replacing in Equation (36) the expression of Q ( 1 ) ( g ; α 0 , φ 0 ; δ α ) from Equation (23) and identifying the expressions that multiply the arbitrary parameter variations α j yields the expression of each partial sensitivity L ( α , φ ; ψ ( 1 ) ) / α j , for all α j , j = 1 , , T P , as follows:

L ( α , φ ; ψ ( 1 ) ) q j = g = 1 G 0 r d 4 π r 2 d r 4 π d Ω ψ ( 1 ) , g ( r , Ω ) Q g ( α ) q j , j = 1 , , J Q = 10 ; (37)

L ( α , φ ; ψ ( 1 ) ) N j = g = 1 G 0 r d 4 π r 2 d r 4 π d Ω ψ ( 1 ) , g ( r , Ω ) Σ t g ( α ) N j φ g ( r , Ω ) h = 1 G 0 r d 4 π r 2 d r 4 π d Ω φ h ( r , Ω ) { Σ s h g ( α ; Ω Ω ) N j + [ χ g ( α ) ( ν Σ f ) h ( α ) ] N j } , for j = 1 , , I = 6 ; (38)

L ( α , φ ; ψ ( 1 ) ) s j = g = 1 G 0 r d 4 π r 2 d r 4 π d Ω ψ ( 1 ) , g ( r , Ω ) × h = 1 G 0 r d 4 π r 2 d r 4 π d Ω Σ s h g ( α ; Ω Ω ) s j φ h ( r , Ω ) , for j = 1 , , J S X = ( G × G ) × I × ( I S C T + 1 ) = 21600 ; (39)

L ( α , φ ; ψ ( 1 ) ) t j = g = 1 G 0 r d 4 π r 2 d r 4 π d Ω ψ ( 1 ) , g ( r , Ω ) φ g ( r , Ω ) Σ t g ( α ) t j , for j = 1 , , J T X = 180 ; (40)

L ( α , φ ; ψ ( 1 ) ) f j = g = 1 G 0 r d 4 π r 2 d r 4 π d Ω ψ ( 1 ) , g ( r , Ω ) χ g × h = 1 G 0 r d 4 π r 2 d r 4 π d Ω [ ( ν Σ f ) h ] f j φ h ( r , Ω ) , for j = 1 , , J F X = 60 ; (41)

L ( α , φ ; ψ ( 1 ) ) ν j = g = 1 G 0 r d 4 π r 2 d r 4 π d Ω ψ ( 1 ) , g ( r , Ω ) χ g × h = 1 G 0 r d 4 π r 2 d r 4 π d Ω [ ( ν Σ f ) h ] ν j φ h ( r , Ω ) , for j = 1 , , J N U = 60 ; (42)

L ( α , φ ; ψ ( 1 ) ) p j = g = 1 G 0 r d 4 π r 2 d r 4 π d Ω ψ ( 1 ) , g ( r , Ω ) χ g ( α ) p j × h = 1 G 0 r d 4 π r 2 d r 4 π d Ω ( ν Σ f ) h φ h ( r , Ω ) , for j = 1 , , J χ = 60. (43)

The expressions of the sensitivities provided in Equations (37) through (43) are to be evaluated at the nominal parameter values α 0 (which has not been indicted explicitly, in order to simplify the respective notations). All of these sensitivities have been computed numerically in [5] - [10]. Evidently, the computations of these sensitivities are inexpensive, involving only integrations using quadrature formulas, once the (1st-level adjoint) function ψ ( 1 ) ( r , Ω ) is available.

8) Equations (32) and (33) are called the First-Level Adjoint Sensitivity System (1st-LASS); the solution ψ ( 1 ) ( r , Ω ) of the 1st-LASS is called the 1st-level adjoint function. The system of equations is called “first-level” as opposed to “first-order” because it does not contain any derivatives (of first- or higher-order) but its solution is used for computing the first-order sensitivities. This terminology will be also used in the sequel, when deriving the expressions for the 2nd- and 3rd-order sensitivities. The 1st-LASS is independent of parameter variations, so it needs to be solved just once to obtain ψ ( 1 ) ( r , Ω ) . Thus, using the 1st-CASAM, the computation of the TP partial sensitivities L ( α , φ ; ψ ( 1 ) ) / α j , j = 1 , , T P , using Equations (37) through requires just a single large-scale computation in order to determine ψ ( 1 ) ( r , Ω ) , followed by TP inexpensive computations to perform each integration (quadrature) involving the 1st-level adjoint function ψ ( 1 ) ( r , Ω ) .

3.4. Comparison of Computational Requirements

The computer used for the PERP computations is a DELL desktop computer (AMD FX-8350) with an 8-core processor. A forward computation of the flux distribution φ ( r , Ω ) within the PERP benchmark and of the leakage response takes ca. 2 minutes CPU-time when the forward transport equation is solved using an angular expansion of order ISN = 256 (S256).

The 1st-level adjoint function ψ ( 1 ) ( r , Ω ) used to compute the first-order sensitivities given by Equations (37) through (43) has been computed in [5] by solving the 1st-LASS using an angular expansion of order ISN = 256 (S256). The CPU-time for a typical adjoint computation of ψ ( 1 ) ( r , Ω ) using S256 is ca. 85 seconds, which is comparable to the CPU-time needed for solving the forward transport equation to compute the flux distribution φ ( r , Ω ) . Furthermore, the typical CPU-time needed to perform the integration (quadrature) over the forward and adjoint functions, to compute the respective sensitivity (using S256) is ca. 0.0012 seconds. The CPU requirements needed for computing the 7477 non-zero sensitivities using forward computations in conjunction with the 2-point finite difference scheme presented in Section 3.1 and using the FSAM presented in Section 3.2, respectively, are compared in Table 1 with the CPU-requirements needed to compute the same sensitivities using the formulas obtained

Table 1. CPU-time needed for computing L ( α ) / α j , j = 1 , , T P = 7477 .

in Equations (37) through (43) in Section 3.3, which were obtained by having applied the 1st-CASAM. The “95 seconds CPU-time” needed by the 1st-CASAM comprises 85 seconds for the one adjoint computation plus 10 seconds for computing all of the 7477 1st-order sensitivities.

The results presented in Table 1 highlights the evident superiority of the 1st-CASAM over the finite-difference method and/or the Forward Sensitivity Analysis Methodology (FSAM) for computing the sensitivities of one response (i.e., the PERP leakage) to the 7477 PERP parameters that have non-zero nominal values.

4. Computation of 2 L ( α ) / t j 1 t j 2 , j 1 = 1 , , J T X ; j 2 = 1 , , j 1

The finite-difference formulas that can be used to compute approximately the second-order response sensitivities with respect to the model parameters will be discussed in Subsection 4.1. The Forward Sensitivity Analysis Methodology will be presented in Subsection 4.2. The application of the 2nd-CASAM will be presented in Subsection 4.3, while Subsection 4.4 will compare the CPU-times which are required by these three deterministic methods.

4.1. Re-computations with Finite-Difference Approximation

The 2nd-order unmixed sensitivities, 2 L ( α ) / ( α j ) 2 , of the leakage response, L ( α ) , can be approximately computed by re-computations using the well-known finite-difference formula presented below:

2 L ( α ) ( α j ) 2 1 h j 2 ( L j + 1 L j + L j 1 ) + Ο ( h j 2 ) , j = 1 , , T P , (44)

where L j + 1 L ( α j + h j ) , L j L ( α j ) , L j 1 L ( α j h j ) and h j denotes a “judiciously-chosen” variation in the parameter α j around its nominal value α j 0 . The 2nd-order mixed sensitivities, 2 L ( α ) / α j 1 α j 2 , can be calculated by using the following finite-difference formula:

2 L ( α ) α j 1 α j 2 1 4 h j 1 h j 2 ( L j 1 + 1 , j 2 + 1 L j 1 1 , j 2 + 1 L j 1 + 1 , j 2 1 + L j 1 1 , j 2 1 ) + Ο ( h j 1 2 , h j 2 2 ) , j 1 = 1 , , T P ; j 2 = 1 , , j 1 , (45)

where h j 1 and h j 2 denote the variations in the parameters α j 1 and α j 2 , respectively. The values of the quantities L j 1 + 1 , j 2 + 1 L ( α j 1 + h j 1 , α j 2 + h j 2 ) , L j 1 1 , j 2 + 1 L ( α j 1 h j 1 , α j 2 + h j 2 ) , L j 1 + 1 , j 2 1 L ( α j 1 + h j 1 , α j 2 h j 2 ) , and L j 1 1 , j 2 1 L ( α j 1 h j 1 , α j 2 h j 2 ) are obtained by re-solving Equations (1) and (2) repeatedly, using the changed parameter values ( α j 1 ± h j 1 ) and ( α j 2 ± h j 2 ) . As also discussed in Subsection 3.1, the values of h j 1 and h j 2 , respectively, must be chosen “judiciously” by “trial and error” for each of the parameters α j 1 and α j 2 ; otherwise, the result produced by Equation (45) will be far off from the exact value of the derivative 2 L ( α ) / α j 1 α j 2 . These simple finite difference formulas introduce their intrinsic “methodological errors” of order Ο ( h j 1 2 ) and/or Ο ( h j 1 2 , h j 2 2 ) , as indicated in Equations (44) and (45), which are in addition to,and independent of, the errors that might be incurred in the computation of L ( α ) . In other words, even if the computation of L ( α ) were perfect (error-free), the finite-difference formulas nevertheless introduce their own, intrinsic, numerical errors into the computation of 2 L ( α ) / α j 1 α j 2 .

4.2. Forward Sensitivity Analysis Methodology

The second-order partial sensitivity of L ( α ) to two generic model parameters, α j 1 and α j 2 is obtained by differentiating the expression provided in Equation (11) with respect to (the second parameter) α j 2 , to obtain the following expression:

{ 2 L ( α ) α j 2 α j 1 } α 0 = S b d S g = 1 G Ω n > 0 d Ω Ω n { 2 φ g ( r , Ω ) α j 1 α j 2 } α 0 , j 1 = 1 , , T P ; j 2 = 1 , , j 1. (46)

The second-order derivative 2 φ g ( r , Ω ) / α j 1 α j 2 is the solution of the second-order G-derivative of Equations (1) and (2) with respect to α j 1 and α j 2 , namely:

{ B g ( α ) [ 2 φ g ( r , Ω ) α j 2 α j 1 ] } α 0 = { 2 S g [ φ g ( r , Ω ) , r , Ω ] α j 2 α j 1 } α 0 , for g = 1 , , G ; j 1 = 1 , , T P ; j 2 = 1 , , j 1 ; (47)

{ 2 φ g ( r , Ω ) α j 2 α j 1 } α 0 = 0 , r = r d ; Ω n < 0 , g = 1 , , G ; j 1 = 1 , , T P ; j 2 = 1 , , j 1 , (48)

where

2 S g [ φ g ( r , Ω ) , r , Ω ] α j 2 α j 1 = α j 2 { S g [ φ g ( r , Ω ) , r , Ω ] α j 1 } = 2 Q g ( r ) α j 2 α j 1 2 Σ t g ( α ; r ) α j 2 α j 1 φ g ( r , Ω ) Σ t g ( α ; r ) α j 1 φ g ( r , Ω ) α j 2 + h = 1 G 4 π d Ω φ h ( r , Ω ) α j 2 { Σ s h g ( α ; r , Ω Ω ) α j 1 + [ χ g ( α ; r ) ( ν Σ f ) h ( α ; r ) ] α j 1 } + h = 1 G 4 π d Ω φ g ( r , Ω ) { 2 Σ s h g ( α ; r , Ω Ω ) α j 2 α j 1 + 2 [ χ g ( α ; r ) ( ν Σ f ) h ( α ; r ) ] α j 2 α j 1 } . (49)

It is important to note that, as indicated in Equation (49), the computation of the 2nd-order flux-derivatives 2 φ g ( r , Ω ) / α j 2 α j 1 requires knowledge of the 1st-order flux derivatives φ g ( r , Ω ) / α j 1 , which appear in the definition of the right-side (“source”) for Equation (47) and which must therefore by computed prior to commencing the computations for obtaining 2 φ g ( r , Ω ) / α j 2 α j 1 . Thus, the computation of the T P ( T P + 1 ) / 2 distinct 2nd-order sensitivities 2 L ( α ) / α j 2 α j 1 requires TP large-scale computations for solving Equations (14) and (15) to determine the fluxes φ g ( r , Ω ) / α j 1 which are needed to compute the right-side of Equation (47), followed by at least T P ( T P + 1 ) / 2 large-scale computations to solve Equations (47) and (48) to determine the quantities 2 φ g ( r , Ω ) / α j 2 α j 1 . Subsequently, T P ( T P + 1 ) / 2 small-scale computations are needed for performing the integrations represented by Equation (46). To summarize: obtaining the numerical values of the 2nd-order sensitivities 2 L ( α ) / α i α j using the FSAM requires a total of T P + T P ( T P + 1 ) / 2 large-scale computations followed by T P ( T P + 1 ) / 2 small-scale computations.

4.3. Second-Order Adjoint Sensitivity Analysis Methodology (2nd-CASAM)

The numerical results for the complete 2nd-order sensitivity analysis of the PERP benchmark’s leakage response with respect to all T P ( T P + 1 ) / 2 nuclear parameters α j , j = 1 , , T P , have been reported in [5] - [10]. Since this work ultimately aims at applying the 4th-Order Comprehensive Sensitivity Analysis Methodology (4th-CASAM) to obtain the closed form expressions of the 4th-order sensitivities of the PERP leakage response with respect to the group-averaged total microscopic cross sections, the proliferation of indices, superscripts and subscripts is unavoidable. Nevertheless, “subscripted-subscripts” can be avoided by using subscripts of the form j 1 = 1 , , J T X ; j 2 = 1 , , j 1 , where the index j1 will replace the index j (which was used in the Section 3) to index the 1st-order sensitivities, and where the index j2 will be used (in addition to j1) to index the 2nd-order sensitivities. Furthermore, anticipating the notation to be used in subsequent Subsections, the index j3 will be used (in addition to j1 and j2) to index the 3rd-order sensitivities; the index j4 will be used (in addition to j1, j2 and j3) to index the 4th-order sensitivities. Recall that even when only the components σ t [ t 1 , , t J T X ] , i = 1 , , I = 6 , g = 1 , , G = 30 , J T X = I × G = 180 , of the group-averaged total microscopic cross sections are considered, the numbers of sensitivities of the PERP leakage response with respect to these parameters are as follows: 1) 180 first-order sensitivities; 2) 32,400 second-order sensitivities, of which 16,290 are distinct; and 3) 5,832,000 third-order sensitivities, of which 988,260 are distinct; 4) 1,049,760,000 fourth-order sensitivities, of which 45,212,895 are distinct.

In contradistinction to the FSAM, the 2nd-CASAM does not use Equation (11) as the starting point for computing 2nd-order response sensitivities but instead uses the expressions of the sensitivities in terms of the 1st-level adjoint function ψ ( 1 ) , as produced by the 1st-CASAM, namely Equations (37) through (43), as individual starting points for deriving the expressions of the corresponding 2nd-order sensitivities. This procedure has been illustrated in [5] - [10] which has derived and computed all of the 2nd-order sensitivities of the PERP leakage response to the nuclear data. In turn, the 3rd-CASAM uses the expressions of the 2nd-order response sensitivities derived within the 2nd-CASAM. The 4th-CASAM will use the expressions of the 3rd-order response sensitivities derived within the 3rd-CASAM, and so on.

In particular, the 2nd-order sensitivities involving only the total microscopic cross sections, are obtained by G-differentiating the expression provided in Equation (40) for the 1st-order sensitivities of the leakage response with respect to the group-averaged total microscopic cross sections, which yields the following result:

δ { L ( α , φ ; ψ ( 1 ) ) t j 1 } α 0 = d d ε { g = 1 G 0 r d 4 π r 2 d r 4 π d Ω { [ ψ ( 1 ) , g ( r , Ω ) + ε δ ψ ( 1 ) , g ( r , Ω ) ] φ g ( r , Ω ) Σ t g ( α ) t j 1 } α 0 { g = 1 G 0 r d 4 π r 2 d r 4 π d Ω ψ ( 1 ) , g ( r , Ω ) [ φ g ( r , Ω ) + ε δ φ g ( r , Ω ) ] Σ t g ( α ) t j 1 } α 0 { g = 1 G 0 r d 4 π r 2 d r 4 π d Ω ψ ( 1 ) , g ( r , Ω ) φ g ( r , Ω ) t j 1 [ Σ t g ( α ) + ε δ Σ t g ( α ) ] } α 0 } ε = 0 , for j 1 = 1 , , J T X = 180. (50)

Carrying out the differentiations with respect to ε in Equation (50) and subsequently setting ε = 0 in the resulting equations yields the following relation:

δ { L ( α , φ ; ψ ( 1 ) ) t j 1 } α 0 = g = 1 G 0 r d 4 π r 2 d r 4 π d Ω δ ψ ( 1 ) , g ( r , Ω ) { φ g ( r , Ω ) Σ t g ( α ) t j 1 } α 0 g = 1 G 0 r d 4 π r 2 d r 4 π d Ω δ φ g ( r , Ω ) { ψ ( 1 ) , g ( r , Ω ) Σ t g ( α ) t j 1 } α 0 g = 1 G 0 r d 4 π r 2 d r 4 π d Ω t j 1 δ Σ t g ( α ) { ψ ( 1 ) , g ( r , Ω ) φ g ( r , Ω ) } α 0 , j 1 = 1 , , J T X = 180. (51)

Taking into account the facts that: 1) the parameter t j 1 represents the group-averaged total microscopic cross sections σ t , i g , with j 1 = ( i 1 ) G + g ,

i = 1 , , I = 6 , g = 1 , , G = 30 ; and 2) Σ t g ( t ) = i = 1 I = 6 N i σ t , i g , it follows that the last term on the left-side of Equation (51) is identically zero.

In Equation (50), the variations δ φ g ( r , Ω ) are the solutions of Equations (21) and (22) while the variations δ ψ ( 1 ) , g ( r , Ω ) are the solutions of the G-differentiated 1st-LASS, namely:

{ A g ( α ) δ ψ ( 1 ) , g ( r , Ω ) } α 0 = { Q ( 2 ) ( g ; α , ψ ( 1 ) ; δ α ) } α 0 , g = 1 , , G , (52)

{ δ ψ ( 1 ) , g ( r , Ω ) } α 0 = 0 , r = r d , Ω n > 0 , g = 1 , , G , (53)

where

Q ( 2 ) ( g ; α , ψ ( 1 ) ; δ α ) { A g ( α ) [ ψ ( 1 ) , g ( r , Ω ) ] } α δ α = ψ ( 1 ) , g ( r , Ω ) Σ t g ( α ; r ) α δ α + h = 1 G 4 π d Ω ψ ( 1 ) , h ( r , Ω ) { Σ s g h ( α ; Ω Ω ) α δ α + [ ( ν Σ f ) g ( α ) χ h ( α ) ] α δ α } , for g = 1 , , G . (54)

In matrix-vector form, which will be used in the sequel, Equations (52) and (53) are expressed as follows:

{ A ( α ) [ δ ψ ( 1 ) ( r , Ω ) ] } α 0 = { Q ( 2 ) ( α , ψ ( 1 ) ; δ α ) } α 0 , (55)

δ ψ ( 1 ) ( r , Ω ) = 0 , r = r d , Ω n > 0. (56)

where

δ ψ ( 1 ) ( r , Ω ) ( δ ψ ( 1 ) , 1 ( r , Ω ) · δ ψ ( 1 ) , G ( r , Ω ) ) , Q ( 2 ) ( α , ψ ( 1 ) ; δ α ) ( Q ( 2 ) ( 1 ; α , ψ ( 1 ) ; δ α ) · Q ( 2 ) ( G ; α , ψ ( 1 ) ; δ α ) ) . (57)

The 2nd-order sensitivities expressed by Equation (51) could be computed, in principle, if the functions δ φ g ( r , Ω ) and δ ψ ( 1 ) , g ( r , Ω ) were available. However, determining the functions δ φ g ( r , Ω ) would require solving Equations (21) and (22), while determining the functions δ ψ ( 1 ) , g ( r , Ω ) would require solving Equations (52) and (53), which would altogether amount to twice as many large-scale computations as there are parameter variations. To avoid the need for performing so many large-scale computations, the right-side of Equation (50) will be recast in terms of the solutions of a 2nd-Level Adjoint Sensitivity System (2nd-LASS) so as to eliminate the appearance of the functions δ φ g ( r , Ω ) and δ ψ ( 1 ) , g ( r , Ω ) . The requisite 2nd-LASS will be constructed by following the general principles introduced by Cacuci [1], which comprise the following sequence of steps, all of which are to be carried out at the nominal parameter values α 0 :

1) Consider a Hilbert space, denoted as H ( 2 ) , comprising block-vector elements u ( 2 ) ( r , Ω ) H ( 2 ) having the structure: u ( 2 ) ( r , Ω ) [ u 1 ( 2 ) ( r , Ω ) , u 2 ( 2 ) ( r , Ω ) ] , where u i ( 2 ) ( r , Ω ) [ u i ( 2 ) , 1 ( r , Ω ) , , u i ( 2 ) , g ( r , Ω ) , , u i ( 2 ) , G ( r , Ω ) ] , i = 1 , 2 ;

2) Define the inner product in H ( 2 ) , denoted as u ( 2 ) ( r , Ω ) , v ( 2 ) ( r , Ω ) ( 2 ) , of two arbitrary elements in H ( 2 ) , denoted as u ( 2 ) ( r , Ω ) H ( 2 ) and v ( 2 ) ( r , Ω ) H ( 2 ) , with v ( 2 ) ( r , Ω ) [ v 1 ( 2 ) ( r , Ω ) , v 2 ( 2 ) ( r , Ω ) ] , and v i ( 2 ) ( r , Ω ) [ v i ( 2 ) , 1 ( r , Ω ) , , v i ( 2 ) , g ( r , Ω ) , , v i ( 2 ) , G ( r , Ω ) ] , i = 1 , 2 , as follows:

u ( 2 ) ( r , Ω ) , v ( 2 ) ( r , Ω ) ( 2 ) i = 1 2 u i ( 2 ) ( r , Ω ) , v i ( 2 ) ( r , Ω ) ( 1 ) = i = 1 2 g = 1 G 0 r d 4 π r 2 d r 4 π d Ω u i ( 2 ) , g ( r , Ω ) v i ( 2 ) , g ( r , Ω ) . (58)

3) For a 2 × 2 block-matrix linear operator L ( 2 ) ( L 11 ( 2 ) L 12 ( 2 ) L 21 ( 2 ) L 22 ( 2 ) ) in the Hilbert space H ( 2 ) , use the inner product defined in Equation (58) to define the formal adjoint operator, denoted as A ( 2 ) ( A 11 ( 2 ) A 12 ( 2 ) A 21 ( 2 ) A 22 ( 2 ) ) , of L ( 2 ) , through the following relationship:

v ( 2 ) , L ( 2 ) u ( 2 ) ( 2 ) = u ( 2 ) , A ( 2 ) v ( 2 ) ( 2 ) + P ( 2 ) [ u ( 2 ) , v ( 2 ) ] , (59)

where P ( 2 ) [ u ( 2 ) , v ( 2 ) ] denotes the corresponding bilinear concomitant evaluated on the domain’s boundary.

4) Apply the definition provided in Equation (58) to form the inner product of Equations (25) and (55) with a yet undefined function ψ ( 2 ) ( j 1 ; r , Ω ) [ ψ 1 ( 2 ) ( j 1 ; r , Ω ) , ψ 2 ( 2 ) ( j 1 ; r , Ω ) ] H ( 2 ) , to obtain the following relation, evaluated at α 0 :

{ ψ ( 2 ) ( j 1 ; r , Ω ) , ( B ( α ) [ δ φ ( r , Ω ) ] A ( α ) [ δ ψ ( 1 ) ( r , Ω ) ] ) ( 2 ) } α 0 = { ψ ( 2 ) ( j 1 ; r , Ω ) , ( Q ( 1 ) ( α ; δ α ) Q ( 2 ) ( α , ψ ( 1 ) ; δ α ) ) ( 2 ) } α 0 . (60)

5) Use the relation shown in Equation (59) to recast the left side of Equation (60) into the following equivalent form:

{ ψ ( 2 ) ( j 1 ; r , Ω ) , ( B ( α ) [ δ φ ( r , Ω ) ] A ( α ) [ δ ψ ( 1 ) ( r , Ω ) ] ) ( 2 ) } α 0 = { [ δ φ , δ ψ ( 1 ) ] , ( A ( α ) ψ 1 ( 2 ) ( j 1 ; r , Ω ) B ( α ) ψ 2 ( 2 ) ( j 1 ; r , Ω ) ) ( 2 ) } α 0 + P ( 2 ) [ δ φ ( r , Ω ) , δ ψ ( 1 ) ( r , Ω ) ; ψ 1 ( 2 ) ( j 1 ; r , Ω ) , ψ 2 ( 2 ) ( j 1 ; r , Ω ) ] , (61)

where P ( 2 ) [ δ φ ( r , Ω ) , δ ψ ( 1 ) ( r , Ω ) ; ψ 1 ( 2 ) ( j 1 ; r , Ω ) , ψ 2 ( 2 ) ( j 1 ; r , Ω ) ] denotes the corresponding bilinear concomitant evaluated on the PERP sphere’s outer boundary S b , at r = r d .

6) Require the first term on the right-side of Equation (61) to represent the same functional as the right-side of Equation (51). This requirement will be satisfied by requiring that the following relations be satisfied by the components of the function ψ ( 2 ) ( j 1 ; r , Ω ) [ ψ 1 ( 2 ) ( j 1 ; r , Ω ) , ψ 2 ( 2 ) ( j 1 ; r , Ω ) ] H ( 2 ) for j 1 = 1 , ... , J T X :

{ A ( α ) ψ 1 ( 2 ) ( j 1 ; r , Ω ) } α 0 = { S ( j 1 ; α ) ψ ( 1 ) ( r , Ω ) } α 0 , j 1 = 1 , , J T X ; (62)

{ B ( α ) ψ 2 ( 2 ) ( j 1 ; r , Ω ) } α 0 = { S ( j 1 ; α ) φ ( r , Ω ) } α 0 ; j 1 = 1 , , J T X , (63)

where, for each j 1 = 1 , , J T X , S ( j 1 ; α ) is a G × G diagonal matrix having non-zero elements of the form Σ t g ( α ) / t j 1 , g = 1 , , G on its diagonal, i.e.,

S ( j 1 ; α ) ( Σ t 1 ( α ) / t j 1 · 0 · · · 0 · Σ t G ( α ) / t j 1 ) D i a g [ Σ t g ( α ) t j 1 ] , g = 1 , , G . (64)

7) Use in Equation (61) the boundary conditions shown in Equations (22) and (53). Furthermore, set to zero the remaining terms in the bilinear concomitant P ( 2 ) [ δ φ ( r , Ω ) , δ ψ ( 1 ) ( r , Ω ) ; ψ 1 ( 2 ) ( j 1 ; r , Ω ) , ψ 2 ( 2 ) ( j 1 ; r , Ω ) ] in Equation (61) by requiring that the components of the function ψ ( 2 ) ( j 1 ; r , Ω ) [ ψ 1 ( 2 ) ( j 1 ; r , Ω ) , ψ 2 ( 2 ) ( j 1 ; r , Ω ) ] H ( 2 ) satisfy the following boundary conditions:

ψ 1 ( 2 ) ( j 1 ; r , Ω ) = 0 , r = r d , Ω n > 0 ; j 1 = 1 , , J T X ; (65)

ψ 2 ( 2 ) ( j 1 ; r , Ω ) = 0 , r = r d , Ω n < 0 ; j 1 = 1 , , J T X . (66)

8) The above boundary conditions complete the well-posed definition of the 2nd-level adjoint function ψ ( 2 ) ( j 1 ; r , Ω ) [ ψ 1 ( 2 ) ( j 1 ; r , Ω ) , ψ 2 ( 2 ) ( j 1 ; r , Ω ) ] H ( 2 ) as the solution of Equations (62) through (66), which are called “the Second-Level Adjoint Sensitivity System (2nd-LASS)”. The reason for calling this system “Second-Level” (as opposed to “Second-Order”) stems from the fact that this system does not involve any “second-order” differential or derivatives of the dependent variables (i.e., state functions) even though the solution of the 2nd-LASS is used for computing, efficiently and exactly, the second-order sensitivities of the response with respect to the model parameters.

9) Use Equations (60) through (66) in Equation (51) to obtain the following alternative expression for the differential δ { L ( α , φ ; ψ ( 1 ) ) / t j 1 } α 0 in terms of the

2nd-level adjoint function ψ ( 2 ) ( j 1 ; r , Ω ) [ ψ 1 ( 2 ) ( j 1 ; r , Ω ) , ψ 2 ( 2 ) ( j 1 ; r , Ω ) ] H ( 2 ) :

δ { L ( α , φ ; ψ ( 1 ) ; ψ ( 2 ) ) t j 1 } α 0

= { ψ ( 2 ) ( j 1 ; r , Ω ) , ( Q ( 1 ) ( α ; δ α ) Q ( 2 ) ( α , ψ ( 1 ) ; δ α ) ) ( 2 ) } α 0 = { g = 1 G 0 r d 4 π r 2 d r 4 π d Ω ψ 1 ( 2 ) , g ( j 1 ; r , Ω ) Q ( 1 ) ( g ; α , φ ; δ α ) } α 0 + { g = 1 G 0 r d 4 π r 2 d r 4 π d Ω ψ 2 ( 2 ) , g ( j 1 ; r , Ω ) Q ( 2 ) ( g ; α , ψ ( 1 ) ; δ α ) } α 0 = j 2 = 1 T P { 2 L ( α ) α j 2 t j 1 δ α j 2 } α 0 , j 1 = 1 , , J T X . (67)

10) It is important to note that the differential expression provided in Equation (67) comprises all of the 2nd-order partial sensitivities that involve the group-averaged microscopic total cross sections. The 2nd-order partial sensitivities { 2 L ( α , φ ; ψ ( 1 ) ) / t j 1 α j 2 } α 0 , j 1 = 1 , , J T X , j 2 = 1 , , T P , which involve the 180 total cross sections ( t j 1 ) and the other parameters ( α j 2 ), are obtained by replacing in Equation (67) the expressions of Q ( 1 ) ( g ; α , φ ; δ α ) and Q ( 2 ) ( g ; α , ψ ( 1 ) ; δ α ) by their definitions provided in Equations (23) and (54), respectively, and subsequently identifying the expressions that multiply the corresponding parameter variations.

11) In particular, identifying in Equation (67) the quantities that multiply only the variations in the total microscopic cross sections provides the following expression for the 2nd-order response sensitivities that involve only the group-averaged total microscopic cross sections:

{ 2 L ( α , φ ; ψ ( 1 ) ; ψ ( 2 ) ) t j 2 t j 1 } α 0 = { g = 1 G 0 r d 4 π r 2 d r 4 π d Ω [ ψ 1 ( 2 ) , g ( j 1 ; r , Ω ) φ g ( r , Ω ) + ψ 2 ( 2 ) , g ( j 1 ; r , Ω ) ψ ( 1 ) , g ( r , Ω ) ] Σ t g ( α ) t j 2 } α 0 , j 1 = 1 , , J T X ; j 2 = 1 , , j 1. (68)

4.4. Comparison of Computational Requirements

The numerical results for the 2nd-order sensitivities of the PERP benchmark’s leakage response to all 7477 nuclear data parameters that characterize the computational model of the PERP benchmark have been presented in [5] - [10]. The 2nd-order sensitivities were all computed using the PARTISN code with a S256 angular quadrature (ISN = 256), using a DELL desktop computer (AMD FX-8350) with an 8-core processor. The CPU-time for a typical adjoint computation was ca. 160 seconds; 0.046 seconds CPU-time was needed for computing a single 2nd-order sensitivity, including reading the adjoint functions and the integration over the adjoint functions, for ISN = 256. Table 2 presents the comparison of the CPU-times required by the finite-difference, FSAM, and 2nd-CASAM, respectively,

Table 2. CPU-time needed for computing 2 L ( α ) / t j 1 t j 2 , j 1 = 1 , , J T X = 180 ; j 2 = 1 , , j 1 .

for computing the 2nd-order sensitivities of the leakage response only with respect to the 180 total group-averaged microscopic cross sections.

The number of 16,470 computations required by the FSAM comprises 180 large-scale computations for obtaining the first-order derivatives of the flux with respect to the total microscopic cross sections plus 180(180 + 1)/2 large-scale computations for computing the mixed second-order derivatives of the flux with respect to the total microscopic cross sections. The total CPU time needed for computing all of the 2nd-order sensitivities using the 2nd-CASAM was ca. 929 hours, comprising 735 hours used for the 14843 PARTISN computations (S256 angular quadrature; ISN = 256) and 194 hours used for performing the integrations needed to compute the respective unmixed and mixed second-order sensitivities. It is evident from the comparison of the CPU-times presented in Table 2 that the 2nd-CASAM is the only practical methodology for computing 2nd-order sensitivities for large-scale systems involving many parameters (as illustrated by using the PERP benchmark as a paradigm).

5. Computation of Third-Order Sensitivities 3 L ( α ) / t j 1 t j 2 t j 3 , j 1 = 1 , , J T X ; j 2 = 1 , , j 1 ; j 3 = 1 , , j 2

Subsection 5.1 presents the finite-difference formulas that can be used to compute approximately the third-order response sensitivities with respect to the model parameters. The Forward Sensitivity Analysis Methodology will be discussed in Subsection 5.2. The 3rd-CASAM will be discussed in Subsection 5.3, while Subsection 5.4 will present a comparative discussion of the computational resources and times (CPU) which are required by these three deterministic methods.

5.1. Re-Computations with Finite-Difference Approximation

The 3rd-order unmixed sensitivities, 3 L ( α ) / α j 3 , of the leakage response, L ( α ) , can be approximately computed by re-computations using the well-known finite-difference formula presented below:

3 L ( α ) α j 3 1 2 h j 3 ( L j + 2 2 L j + 1 + 2 L j 1 L j 2 ) + Ο ( h j 2 ) , j = 1 , , T P , (69)

where L j + 2 L ( α j + 2 h j ) , L j + 1 L ( α j + h j ) , L j 1 L ( α j h j ) , L j 2 L ( α j 2 h j ) and h j denotes a “judiciously-chosen” variation in the parameter α j around its nominal value α j 0 . The 3rd-order mixed sensitivities, 3 L ( α ) / α j 1 α j 2 α j 3 , can be calculated by using the following finite-difference formula:

3 L ( α ) α j 1 α j 2 α j 3 1 8 h j 1 h j 2 h j 3 ( L j 1 + 1 , j 2 + 1 , j 3 + 1 L j 1 + 1 , j 2 + 1 , j 3 1 L j 1 + 1 , j 2 1 , j 3 + 1 + L j 1 + 1 , j 2 1 , j 3 1 L j 1 1 , j 2 + 1 , j 3 + 1 + L j 1 1 , j 2 + 1 , j 3 1 + L j 1 1 , j 2 1 , j 3 + 1 L j 1 1 , j 2 1 , j 3 1 ) + Ο ( h j 1 2 , h j 2 2 , h j 3 2 ) , (70)

where L j 1 + 1 , j 2 + 1 , j 3 + 1 = L ( α j 1 + h j 1 , α j 2 + h j 2 , α j 3 + h j 3 ) , L j 1 1 , j 2 + 1 , j 3 + 1 = L ( α j 1 h j 1 , α j 2 + h j 2 , α j 3 + h j 3 ) , etc. The values of the quantities on the rights-sides of the expressions shown in Equations (69) and (70) are obtained by re-solving Equations (1) and (2) repeatedly, using the changed parameter values ( α j 1 ± h j 1 ) , ( α j 1 ± h j 2 ) and ( α j 3 ± h j 3 ) . As has also been previously discussed, the values of h j 1 , h j 2 and h j 3 , respectively, must be chosen “judiciously” by “trial and error” for each of the parameters α j 1 , α j 2 and α j 3 . The finite difference formulas introduce their intrinsic “methodological errors” of order Ο ( h j 1 2 , h j 2 2 , h j 3 2 ) which are in addition to, and independent of, the errors that might be incurred in the computation of 3 L ( α ) / α j 1 α j 2 α j 3 .

5.2. Forward Sensitivity Analysis Methodology

The third-order partial sensitivity of L ( α ) to three generic model parameters, α j 3 , α j 2 , α j 1 , has the following expression:

{ 3 L ( α ) α j 3 α j 2 α j 1 } α 0 { S b d S g = 1 G Ω n > 0 d Ω Ω n 3 φ g ( r , Ω ) α j 3 α j 2 α j 1 } α 0 , j 1 = 1 , , T P ; j 2 = 1 , , j 1 ; j 3 = 1 , , j 2. (71)

In turn, the third-order derivative 3 φ g ( r , Ω ) / α j 1 α j 2 α j 3 is the solution of the third-order G-derivative of Equations (1) and (2) with respect to α j 3 , α j 2 , α j 1 or, equivalently, the solution of the first G-derivative of the 2nd-OFSS comprising Equations (47) and (48), which have the following form:

{ B g ( α ) [ 3 φ g ( r , Ω ) α j 3 α j 2 α j 1 ] } α 0 = { 3 S g [ φ g ( r , Ω ) , r , Ω ] α j 3 α j 2 α j 1 } α 0 , g = 1 , , G ; j 1 = 1 , , T P ; j 2 = 1 , , j 1 ; j 3 = 1 , , j 2 ; (72)

{ 3 φ g ( r , Ω ) α j 3 α j 2 α j 1 } α 0 = 0 , r = r d ; Ω n < 0 , g = 1 , , G ; j 1 = 1 , , T P ; j 2 = 1 , , j 1 ; j 3 = 1 , , j 2 ; (73)

where

3 S g [ φ g ( r , Ω ) , r , Ω ] α j 3 α j 2 α j 1 = α j 3 { 2 S g [ φ g ( r , Ω ) , r , Ω ] α j 2 α j 1 } . (74)

Evidently, the computations of all of the distinctive third-order sensitivities 3 L ( α ) / α j 3 α j 2 α j 1 require T P ( T P + 1 ) ( T P + 2 ) / 3 ! large-scale computations to solve Equations (72) and (73), followed by T P ( T P + 1 ) ( T P + 2 ) / 3 ! small-scale computations for performing the integration represented by Equation (71).

5.3. Third-Order Adjoint Sensitivity Analysis Methodology (3rd-CASAM)

In contradistinction to the FSAM, the 3rd-CASAM does not use Equation (46) as the starting point for computing 3rd-order response sensitivities but instead uses the expression provided by the 2nd-CASAM of the sensitivities in terms of the 1st-and 2nd-level adjoint functions ψ ( 1 ) and ψ ( 2 ) . Thus, the total G-differential of the expression given in Equation (68) provides the 3rd-order sensitivities of the leakage response involving the group-averaged total microscopic cross sections, as follows:

δ { 2 L ( α , φ ; ψ ( 1 ) ; ψ ( 2 ) ) t j 2 t j 1 } α 0 = { g = 1 G 0 r d 4 π r 2 d r 4 π d Ω [ δ ψ 1 ( 2 ) , g ( j 1 ; r , Ω ) φ g ( r , Ω ) + ψ 1 ( 2 ) , g ( j 1 ; r , Ω ) δ φ g ( r , Ω ) + δ ψ 2 ( 2 ) , g ( j 1 ; r , Ω ) ψ ( 1 ) , g ( r , Ω ) + ψ 2 ( 2 ) , g ( j 1 ; r , Ω ) δ ψ ( 1 ) , g ( r , Ω ) ] Σ t g ( α ) t j 2 } α 0 , j 1 = 1 , , J T X ; j 2 = 1 , , j 1 , (75)

where: 1) the functions δ φ g ( r , Ω ) are the solutions of Equations (25) and (26); 2) the functions δ ψ ( 1 ) , g ( r , Ω ) are the solutions of Equations (55) and (56); 3) the functions δ ψ 1 , j ( 2 ) , g and δ ψ 2 , j ( 2 ) , g are the solutions of the G-differentiated 2nd-LASS defined by Equations (62) through (66), namely:

{ A ( α ) δ ψ 1 ( 2 ) ( j 1 ; r , Ω ) + S ( j 1 ; α ) δ ψ ( 1 ) ( r , Ω ) } α 0 = { Q 1 ( 3 ) [ α ; ψ ( 1 ) ( r , Ω ) ; ψ 1 ( 2 ) ( j 1 ; r , Ω ) ; δ α ] } α 0 , (76)

{ B ( α ) δ ψ 2 ( 2 ) ( j 1 ; r , Ω ) + S ( j 1 ; α ) δ φ ( r , Ω ) } α 0 = { Q 2 ( 3 ) [ α ; φ ( r , Ω ) ; ψ 2 ( 2 ) ( j 1 ; r , Ω ) ; δ α ] } α 0 , (77)

δ ψ 1 ( 2 ) ( j 1 ; r , Ω ) = 0 , r = r d , Ω n > 0 ; j 1 = 1 , , J T X ; (78)

δ ψ 2 ( 2 ) ( j 1 ; r , Ω ) = 0 , r = r d , Ω n < 0 ; j 1 = 1 , , J T X , (79)

where

Q 1 ( 3 ) [ α ; ψ ( 1 ) ( r , Ω ) ; ψ 1 ( 2 ) ( j 1 ; r , Ω ) ; δ α ] { [ S ( j 1 ; α ) α δ α ] ψ ( 1 ) ( r , Ω ) } α 0 { [ A ( α ) α δ α ] ψ 1 ( 2 ) ( j 1 ; r , Ω ) } α 0 j 3 = 1 T P Q 1 ( 3 ) [ α ; ψ ( 1 ) ( r , Ω ) ; ψ 1 ( 2 ) ( j 1 ; r , Ω ) ] α j 3 δ α j 3 , (80)

Q 2 ( 3 ) [ α ; φ ( r , Ω ) ; ψ 2 ( 2 ) ( j 1 ; r , Ω ) ; δ α ] { [ S ( j 1 ; α ) α δ α ] φ ( r , Ω ) } α 0 { [ B ( α ) α δ α ] ψ 2 ( 2 ) ( j 1 ; r , Ω ) } α 0 = j 3 = 1 T P Q 2 ( 3 ) α j 3 δ α j 3 . (81)

Thus, in order to compute the value of the differential shown in Equation (75), it would be first necessary to solve the following system of matrix-equations, subject to the boundary conditions shown in Equations (26), (56), (78) and (79):

{ ( B ( α ) 0 0 0 0 A ( α ) 0 0 0 S ( j 1 ; α ) A ( α ) 0 S ( j 1 ; α ) 0 0 B ( α ) ) } α 0 ( δ φ δ ψ ( 1 ) δ ψ 1 ( 2 ) ( j 1 ) δ ψ 2 ( 2 ) ( j 1 ) ) = { ( Q ( 1 ) ( α ; δ α ) Q ( 2 ) ( α , ψ ( 1 ) ; δ α ) Q 1 ( 3 ) [ α ; ψ ( 1 ) ; ψ 1 ( 2 ) ( j 1 ) ; δ α ] Q 2 ( 3 ) [ α ; φ ; ψ 2 ( 2 ) ( j 1 ) ; δ α ] ) } α 0 . (82)

It is evident that solving Equation (82) amounts to impractically many large-scale computations, which can be avoided by recasting the right-side of Equation (75) in terms of the solutions of a 3rd-Level Adjoint Sensitivity System (3rd-LASS). Following the general principles introduced by Cacuci [12], the construction of the 3rd-LASS is similar to the construction of the 2nd-LASS and comprises the following sequence of operations, all of which are to be carried out at the nominal parameter values α 0 :

1) Consider a Hilbert space, denoted as H ( 3 ) , comprising elements denoted generically as u ( 3 ) ( r , Ω ) H ( 3 ) , which have a “four-component block-vector” structure of the form u ( 3 ) ( r , Ω ) [ u 1 ( 3 ) ( r , Ω ) , u 2 ( 3 ) ( r , Ω ) , u 3 ( 3 ) ( r , Ω ) , u 4 ( 3 ) ( r , Ω ) ] , with each of the four vector-components having the following structure: u i ( 3 ) ( r , Ω ) [ u i ( 3 ) , 1 ( r , Ω ) , , u i ( 3 ) , g ( r , Ω ) , , u i ( 3 ) , G ( r , Ω ) ] , i = 1 , 2 , 3 , 4 .

2) Define the inner product in H ( 3 ) , denoted as u ( 3 ) ( r , Ω ) , v ( 3 ) ( r , Ω ) ( 3 ) , of two arbitrary elements in H ( 3 ) , denoted as u ( 3 ) ( r , Ω ) H ( 3 ) and v ( 3 ) ( r , Ω ) H ( 3 ) , with v ( 3 ) ( r , Ω ) [ v 1 ( 3 ) ( r , Ω ) , v 2 ( 3 ) ( r , Ω ) , v 3 ( 3 ) ( r , Ω ) , v 4 ( 3 ) ( r , Ω ) ] , and v i ( 3 ) ( r , Ω ) [ v i ( 3 ) , 1 ( r , Ω ) , , v i ( 3 ) , g ( r , Ω ) , , v i ( 3 ) , G ( r , Ω ) ] , i = 1 , 2 , 3 , 4 , as follows:

u ( 3 ) ( r , Ω ) , v ( 3 ) ( r , Ω ) ( 3 ) i = 1 4 u i ( 3 ) ( r , Ω ) , v i ( 3 ) ( r , Ω ) ( 1 ) = i = 1 4 g = 1 G 0 r d 4 π r 2 d r 4 π d Ω u i ( 3 ) , g ( r , Ω ) v i ( 3 ) , g ( r , Ω ) . (83)

3) Apply the definition provided in Equation (83) to form the inner product of Equation (82) with a yet undefined function ψ ( 3 ) ( j 2 ; j 1 ; r , Ω ) [ ψ 1 ( 3 ) ( j 2 ; j 1 ; r , Ω ) , , ψ 4 ( 3 ) ( j 2 ; j 1 ; r , Ω ) ] H ( 3 ) , to obtain the following relation, evaluated at α 0 :

{ ψ ( 3 ) ( j 2 ; j 1 ; r , Ω ) , ( B ( α ) 0 0 0 0 A ( α ) 0 0 0 S ( j 1 ; α ) A ( α ) 0 S ( j 1 ; α ) 0 0 B ( α ) ) ( δ φ δ ψ ( 1 ) δ ψ 1 ( 2 ) ( j 1 ) δ ψ 2 ( 2 ) ( j 1 ) ) ( 3 ) } α 0 = { ψ ( 3 ) ( j 2 ; j 1 ; r , Ω ) , ( Q ( 1 ) ( α ; δ α ) Q ( 2 ) ( α , ψ ( 1 ) ; δ α ) Q 1 ( 3 ) [ α ; ψ ( 1 ) ; ψ 1 ( 2 ) ( j 1 ) ; δ α ] Q 2 ( 3 ) [ α ; φ ; ψ 2 ( 2 ) ( j 1 ) ; δ α ] ) ( 3 ) } α 0 . (84)

4) Use the customary definition of the adjoint operator in H ( 3 ) , endowed with the inner product defined in Equation (83) to recast the left side of Equation (84) into the following equivalent form:

{ ψ ( 3 ) ( j 2 ; j 1 ; r , Ω ) , ( B ( α ) 0 0 0 0 A ( α ) 0 0 0 S ( j 1 ; α ) A ( α ) 0 S ( j 1 ; α ) 0 0 B ( α ) ) ( δ φ δ ψ ( 1 ) δ ψ 1 ( 2 ) ( j 1 ) δ ψ 2 ( 2 ) ( j 1 ) ) ( 3 ) } α 0 = { ( δ φ δ ψ ( 1 ) δ ψ 1 ( 2 ) ( j 1 ) δ ψ 2 ( 2 ) ( j 1 ) ) , ( A ( α ) 0 0 S ( j 1 ; α ) 0 B ( α ) S ( j 1 ; α ) 0 0 0 B ( α ) 0 0 0 0 A ( α ) ) ( ψ 1 ( 3 ) ( j 2 ; j 1 ; r , Ω ) ψ 2 ( 3 ) ( j 2 ; j 1 ; r , Ω ) ψ 3 ( 3 ) ( j 2 ; j 1 ; r , Ω ) ψ 4 ( 3 ) ( j 2 ; j 1 ; r , Ω ) ) ( 3 ) } α 0 + P ( 3 ) [ δ φ , δ ψ ( 1 ) ; δ ψ ( 2 ) ( j 1 ) ; ψ ( 3 ) ( j 2 ; j 1 ) ] , (85)

where P ( 3 ) [ δ φ , δ ψ ( 1 ) ; δ ψ ( 2 ) ( j 1 ) ; ψ ( 3 ) ( j 2 ; j 1 ) ] denotes the corresponding bilinear concomitant evaluated on the PERP sphere’s outer boundary S b , at r = r d .

5) Require the first term on the right-side of Equation (85) to represent the same functional as the right-side of Equation (75). This requirement will be satisfied by requiring that the following relations be satisfied by the components of the function ψ ( 3 ) ( j 2 ; j 1 ; r , Ω ) [ ψ 1 ( 3 ) ( j 2 ; j 1 ; r , Ω ) , , ψ 4 ( 3 ) ( j 2 ; j 1 ; r , Ω ) ] H ( 3 ) :

{ ( A ( α ) 0 0 S ( j 1 ; α ) 0 B ( α ) S ( j 1 ; α ) 0 0 0 B ( α ) 0 0 0 0 A ( α ) ) ( ψ 1 ( 3 ) ( j 2 ; j 1 ; r , Ω ) ψ 2 ( 3 ) ( j 2 ; j 1 ; r , Ω ) ψ 3 ( 3 ) ( j 2 ; j 1 ; r , Ω ) ψ 4 ( 3 ) ( j 2 ; j 1 ; r , Ω ) ) } α 0

= { ( S ( j 2 ; α ) ψ 1 ( 2 ) ( j 1 ; r , Ω ) S ( j 2 ; α ) ψ 2 ( 2 ) ( j 1 ; r , Ω ) S ( j 2 ; α ) φ ( r , Ω ) S ( j 2 ; α ) ψ ( 1 ) ( r , Ω ) ) } α 0 , for j 1 = 1 , , J T X ; j 2 = 1 , , j 1. (86)

6) Use in Equation (86) the boundary conditions shown in Equations (26), (56), (78) and (79). Furthermore, set to zero the remaining terms in the bilinear concomitant P ( 3 ) [ δ φ , δ ψ ( 1 ) ; δ ψ ( 2 ) ( j 1 ) ; ψ ( 3 ) ( j 2 ; j 1 ) ] in Equation (85) by requiring the components of the function ψ ( 3 ) ( j 2 ; j 1 ; r , Ω ) [ ψ 1 ( 3 ) ( j 2 ; j 1 ; r , Ω ) , , ψ 4 ( 3 ) ( j 2 ; j 1 ; r , Ω ) ] H ( 3 ) to satisfy the following boundary conditions:

ψ i ( 3 ) ( j 2 ; j 1 ; r , Ω ) = 0 , r = r d ; Ω n < 0 ; i = 2 , 3 ; j 1 = 1 , , J T X ; j 2 = 1 , , j 1 ; ψ i ( 3 ) ( j 2 ; j 1 ; r , Ω ) = 0 , r = r d ; Ω n > 0 ; i = 1 , 4 ; j 1 = 1 , , J T X ; j 2 = 1 , , j 1. (87)

7) The boundary conditions shown in Equation (87) complete the well-posed definition of the 3rd-level adjoint function ψ ( 3 ) ( j 2 ; j 1 ; r , Ω ) H ( 3 ) as the solution of Equations (86) and (87), which are called “the Third-Level Adjoint Sensitivity System (3rd-LASS)”. The reason for calling this system “Third-Level” (as opposed to “Third-Order”) stems from the fact that this system does not involve any 2nd- and/or 3rd-order differentials or derivatives of the dependent variables (i.e., state functions) even though the solution of the 3rd-LASS is used for computing, efficiently and exactly, the 3rd-order sensitivities of the response with respect to the model parameters.

8) Use Equations (84) through (87) in Equation (75) to obtain the following alternative expression for the differential δ { 2 L ( α , φ ; ψ ( 1 ) ) / t j 1 t j 2 } α 0 in terms of the 3rd-level adjoint function ψ ( 3 ) ( j 2 ; j 1 ; r , Ω ) H ( 3 ) :

δ { 2 L ( α , φ ; ψ ( 1 ) ; ψ ( 2 ) ) t j 2 t j 1 } α 0 = { ψ 1 ( 3 ) ( j 2 ; j 1 ; r , Ω ) , Q ( 1 ) ( α ; δ α ) ( 1 ) } α 0 + { ψ 2 ( 3 ) ( j 2 ; j 1 ; r , Ω ) , Q ( 2 ) ( α , ψ ( 1 ) ; δ α ) ( 1 ) } α 0 + { ψ 3 ( 3 ) ( j 2 ; j 1 ; r , Ω ) , Q 1 ( 3 ) [ α ; ψ ( 1 ) ; ψ 1 ( 2 ) ( j 1 ) ; δ α ] ( 1 ) } α 0 + { ψ 4 ( 3 ) ( j 2 ; j 1 ; r , Ω ) , Q 2 ( 3 ) [ α ; ψ ( 1 ) ; ψ 2 ( 2 ) ( j 1 ) ; δ α ] ( 1 ) } α 0 = j 3 = 1 T P { 3 L ( α , φ ; ψ ( 1 ) ; ψ ( 2 ) ; ψ ( 3 ) ) α j 3 t j 2 t j 1 δ α j 3 } α 0 , j 1 = 1 , , J T X ; j 2 = 1 , , j 1. (88)

9) It is important to note that the differential expression provided in Equation (88) comprises all of the 3rd-order partial sensitivities that involve the group-averaged microscopic total cross sections, i.e., the 3rd-order partial sensitivities 3 L ( α , φ ; ψ ( 1 ) ; ψ ( 2 ) ; ψ ( 3 ) ) / α j 3 t j 2 t j 1 , j 1 = 1 , , J T X ; j 2 = 1 , , j 1 , j 3 = 1 , , T P , describes a total of [ 180 ( 180 + 1 ) / 2 ] × 7477 3rd-order response sensitivities with respect to the microscopic total cross sections (indices t j 1 and t j 2 ) and with respect to the other parameters ( α j 3 ).

10) In particular, identifying in Equation (88) the quantities that multiply only the variations in the total microscopic cross sections will reduce the expressions provided in Equations (80) and (81) as follows:

Q 1 ( 3 ) [ α ; ψ ( 1 ) ( r , Ω ) ; ψ 1 ( 2 ) ( j 1 ; r , Ω ) ; δ α ] = j 3 = 1 J T X S 1 ( 3 ) [ α ; ψ 1 ( 2 ) ( j 1 ; r , Ω ) ] t j 3 δ t j 3 , S 1 ( 3 ) [ α ; ψ 1 ( 2 ) ( j 1 ; r , Ω ) ] Σ t g ( α ) ψ 1 ( 2 ) ( j 1 ; r , Ω ) , (89)

Q 2 ( 3 ) [ α ; φ ( r , Ω ) ; ψ 2 ( 2 ) ( j 1 ; r , Ω ) ; δ α ] j 3 = 1 J T X S 2 ( 3 ) [ α ; ψ 2 ( 2 ) ( j 1 ; r , Ω ) ] t j 3 δ t j 3 , S 2 ( 3 ) [ α ; ψ 2 ( 2 ) ( j 1 ; r , Ω ) ] Σ t g ( α ) ψ 2 ( 2 ) ( j 1 ; r , Ω ) . (90)

11) Inserting the results from Equations (89) and (90) into Equation (88) and identifying the specific quantity that multiplies the specific variation δ t j 3 provides the following expression for the 3rd-order response sensitivities which involve only the group-averaged total microscopic cross sections:

{ 3 L ( α , φ ; ψ ( 1 ) ; ψ ( 2 ) ; ψ ( 3 ) ) t j 3 t j 2 t j 1 } α 0 = { ψ 1 ( 3 ) ( j 2 ; j 1 ; r , Ω ) , S ( j 3 ; α ) φ ( r , Ω ) ( 1 ) } α 0 { ψ 2 ( 3 ) ( j 2 ; j 1 ; r , Ω ) , S ( j 3 ; α ) ψ ( 1 ) ( r , Ω ) ( 1 ) } α 0 { ψ 3 ( 3 ) ( j 2 ; j 1 ; r , Ω ) , S ( j 3 ; α ) ψ 1 ( 2 ) ( j 1 ; r , Ω ) ( 1 ) } α 0 { ψ 4 ( 3 ) ( j 2 ; j 1 ; r , Ω ) , S ( j 3 ; α ) ψ 2 ( 2 ) ( j 1 ; r , Ω ) ( 1 ) } α 0 , for j 1 = 1 , , J T X ; j 2 = 1 , , j 1 ; j 3 = 1 , , j 2. (91)

In component form, Equation (91) can be written as follows:

{ 3 L ( α , φ ; ψ ( 1 ) ; ψ ( 2 ) ; ψ ( 3 ) ) t j 3 t j 2 t j 1 } α 0 = { g = 1 G 0 r d 4 π r 2 d r 4 π d Ω [ ψ 1 ( 3 ) , g ( j 2 ; j 1 ; r , Ω ) φ g ( r , Ω ) + ψ 2 ( 3 ) , g ( j 2 ; j 1 ; r , Ω ) ψ ( 1 ) , g ( r , Ω ) + ψ 3 ( 3 ) , g ( j 2 ; j 1 ; r , Ω ) ψ 1 ( 2 ) , g ( j 1 ; r , Ω ) + ψ 4 ( 3 ) , g ( j 2 ; j 1 ; r , Ω ) ψ 2 ( 2 ) , g ( j 1 ; r , Ω ) ] Σ t g ( α ) t j 3 } α 0 , for j 1 = 1 , , J T X ; j 2 = 1 , , j 1 ; j 3 = 1 , , j 2. (92)

5.4. Comparison of Computational Requirements

As has already been discussed, the number of sensitivities of the PERP leakage response which involve solely the group-averaged total microscopic cross sections are as follows: 1) 180 first-order sensitivities; 2) 32,400 second-order sensitivities, of which 16,290 are distinct; and 3) 5,832,000 third-order sensitivities, of which 988,260 are distinct; 4) 1,049,760,000 fourth-order sensitivities, of which 45,212,895 are distinct. Numerical results for these sensitivities have been obtained by using the code PARTISN with a S32 angular quadrature (ISN = 32). The computer used for these computations is a DELL computer with an 8-core processor (AMD FX-8350).

The CPU-time for a typical adjoint computation using ISN = 32 is ca. 26 seconds, while the CPU-time computing the integrals over the various adjoint functions which appear in the definition of the respective sensitivity, cf. Equation (92), is ca. 0.002 seconds.

The CPU-time for computing the 1st-, 2nd- and 3rd-order sensitivities of the leakage response with respect to the 180 microscopic total cross sections are as follows:

1) The total CPU-time required for computing the 1st-order sensitivities is ca. 30 Seconds (1 adjoint PARTISN computation, plus 180 integrations over the adjoint functions);

2) The total CPU-time required for computing twice (for verification purposes) the 2nd-order sensitivities is ca. 3 hours (2 h 55 min for 180 forward-like + 180 adjoint-like PARTISN computations for determining the 2nd-level adjoint functions, plus 5 minutes for computing the 32,400 integrations over the adjoint functions);

3) The total CPU-time required for computing twice (for verification purposes) the 3rd-order sensitivities is around 338 hours (2 h 55 min for computing 2nd-level adjoint functions, 332 hours for computing the 3rd-level adjoint functions + 3 h 24 min for computing the 180 × 180 × 180 integrals over these functions to obtain the 3rd-order sensitivities).

Table 3 presents the comparison of the CPU-times for all of the distinct third-order sensitivities 3 L ( α ) / t j 3 t j 2 t j 1 , j 1 = 1 , , J T X ; j 2 = 1 , , j 1 ; j 3 = 1 , , j 2 required by applying the 3rd-CASAM, versus using the FSAM procedure using Equation (71) and/or the FD-approximations shown in Equations (69) and (70).

Table 3. Comparison of CPU-times for computing 3 L ( α , φ ; ψ ( 1 ) ; ψ ( 2 ) ; ψ ( 3 ) ) / t j 3 t j 2 t j 1 , j 1 = 1 , , J T X ; j 2 = 1 , , j 1 ; j 3 = 1 , , j 2 .

Based on Equations (69) and (70), the number of forward computations needed to obtain all of the distinct third-order sensitivities by the FD-approximation method is 7,905,360 [= 180 × 4 forward PARTISN computations for the 180 unmixed 3rd-order sensitivities + (988,260 − 180) × 8 forward PARTISN computations for the 988,080 mixed 3rd-order sensitivities].

Based on Equations (72) and (73), the number of forward computations needed to obtain all of the distinct third-order sensitivities 3 L ( α ) / t j 3 t j 2 t j 1 by the FSAM method is J T X + J T X ( J T X + 1 ) / 2 + J T X ( J T X + 1 ) ( J T X + 2 ) / 3 ! = 1004730 .

Based on Equation (86), a total number of 33,301 adjoint computations is needed to obtain all of the distinct third-order sensitivities 3 L ( α ) / t j 3 t j 2 t j 1 by the 3rd-CASAM method, comprising the following large-scale computations:

1) 1 adjoint PARTISN computations to obtain ψ ( 1 ) ( r , Ω ) ;

2) 180 adjoint-like PARTISN computations to obtain ψ 1 ( 2 ) ( j 1 ; r , Ω ) ;

3) 180 forward-like PARTISN computations to obtain ψ 2 ( 2 ) ( j 1 ; r , Ω ) ;

4) 16,290 adjoint-like PARTISN computations to obtain ψ 1 ( 3 ) ( j 2 ; j 1 ; r , Ω ) ;

5) 16,290 forward-like PARTISN computations to obtain ψ 2 ( 3 ) ( j 2 ; j 1 ; r , Ω ) ;

6) 180 forward-like PARTISN computations to obtain ψ 3 ( 3 ) ( j 2 ; j 1 ; r , Ω ) ;

7) 180 adjoint-like PARTISN computations to obtain ψ 4 ( 3 ) ( j 2 ; j 1 ; r , Ω ) .

The computation of each of the adjoint functions ψ 1 ( 3 ) ( j 2 ; j 1 ; r , Ω ) and ψ 2 ( 3 ) ( j 2 ; j 1 ; r , Ω ) , respectively, requires 180 ( 180 + 1 ) / 2 = 16290 PARTISN computations, respectively, since the source terms for the equations satisfied by these functions depend on both indices j1 and j2. However, the computation of each of the adjoint functions ψ 3 ( 3 ) ( j 2 ; j 1 ; r , Ω ) and ψ 4 ( 3 ) ( j 2 ; j 1 ; r , Ω ) , respectively, requires only 180 PARTISN computations (per function), because the source terms for the equations satisfied by these functions only depend on the index j2, as shown in Equation .

The results presented in Table 3 evidently indicate that the 3rd-CASAM is the only methodology capable of computing 3rd-order sensitivities, without introducing methodological errors.

6. Concluding Remarks

This work has presented the derivation of the exact mathematical expressions of the following sensitivities of the PERP leakage response with respect to the group-averaged total microscopic cross sections: 1) 180 first-order sensitivities; 2) 32,400 second-order sensitivities, of which 16,290 are distinct; and 3) 5,832,000 third-order sensitivities, of which 988,260 are distinct.

Using a DELL computer with an 8-core processor (AMD FX-8350), the total CPU-time required for computing the 1st-order sensitivities of the PERP leakage response with respect to all of the PERP benchmark’s parameters by applying the 1st-CASAM is 95 seconds. By comparison, the forward sensitivity analysis methodology (FSAM) would require 694 hours CPU-time, while a 2-point finite-difference (FD) approximate scheme would require 1,388 hours CPU-time.

Using the same DELL computer (AMD FX-8350), the computations of the 2nd-order sensitivities of the PERP leakage response with respect to the PERP benchmark’s total microscopic cross sections by applying the 2nd-CASAM required 19 hours CPU-time. By comparison, the FSAM would require 764 hours CPU-time, while the FD-scheme would require 3006 hours CPU-time.

Finally, using the AMD FX-8350 DELL computer, the computations of the 3rd-order sensitivities of the PERP leakage response with respect to the PERP benchmark’s total microscopic cross sections by applying the 3rd-CASAM required 175 hours CPU-time. By comparison, the FSAM would require 12,559 hours CPU-time, while the FD-scheme would require 98,817 hours CPU-time. These results make it abundantly evident that the CASAM is the only practical method for computing sensitivities exactly for large-scale problems involving many parameters.

The formulas derived in this work are valid not only for the PERP benchmark but can also be used for computing the 3rd-order sensitivities of the leakage response of any nuclear system involving fissionable material and internal or external neutron sources. Subsequent works [16] [17] will use the adjoint-based mathematical expressions obtained in this work to further derive the expressions for the 4th-order sensitivities, and to compute exactly and efficiently the numerical values of the 1,049,760,000 fourth-order sensitivities, of which 45,212,895 are distinct.

Appendix

The dimensions and material composition of the polyethylene-reflected plutonium (PERP) metal sphere considered in this work are presented in TableA1.

The quantities appearing in Equations (1) through (3) are defined as follows:

1) The quantity φ g ( r , Ω ) is the customary “group-flux” for group g , g = 1 , , G , and is the unknown state-function which is obtained by solving Equations (1) and (2). The group boundaries of the G=30 energy groups are provided in TableA2.

2) The source Q g ( q ) depends on the vector of model parameters q , defined as follows:

Table A1. Dimensions and material composition of the PERP benchmark.

Table A2. Group boundaries, E g [ MeV ] , of the G = 30 energy groups used in the PARTISN forward and adjoint neutron transport computations.

q [ q 1 , , q J Q ] [ λ 1 , λ 2 ; F 1 S F , F 2 S F ; a 1 , a 2 ; b 1 , b 2 ; ν 1 S F , ν 2 S F ] , J Q = 10. (93)

3) As indicated in TableA1, the PERP benchmark comprises 2 materials: “material 1” comprises 4 isotopes, numbered 1 through 4, while “material 2” comprises 2 isotopes, numbered 5 and 6. In principle, PARTISN allows the same isotope to appear in different materials, in which case the atomic number density N i , m of an isotope i in a material m would be computed by using the formula N i , m = ρ m w i , m N A / A i , where ρ m denotes the mass density of material m, m = 1 , 2 , w i , m denotes the weight fraction of isotope i in material m; A i denotes the atomic weight of isotope i, and N A denotes Avogadro’s number. However, the two materials in the PERP benchmark contain only isotopes that are distinct from each other, so the subscript m will not be needed if the formula for the atomic number density N i of an isotope i, i = 1 , , I = 6 , is interpreted as follows:

N i ρ 1 w i , 1 N A A i ; for i = 1 , 2 , 3 , 4 ; N i ρ 2 w i , 2 N A A i ; for i = 5 , 6. (94)

The atomic number densities N i , i = 1 , , I = 6 will be considered to be components of the vector N , defined below:

N [ N 1 , N 2 , N 3 , N 4 , N 5 , N 6 ] . (95)

4) The scattering transfer cross section from energy group g , g = 1 , , G into energy group g , g = 1 , , G is denoted as Σ s g g ( α ; r , Ω Ω ) and is computed in terms of the l-th order Legendre coefficient σ s , l , i g g using the following 3rd-order expansion in Legendre functions:

Σ s g g ( α ; r , Ω Ω ) = i = 1 I = 6 N i l = 0 I S C T = 3 ( 2 l + 1 ) σ s , l , i g g ( r ) P l ( Ω Ω ) , g , g = 1 , , G , (96)

where ISCT = 3 denotes the order of the expansion in Legendre polynomials. The microscopic scattering cross sections σ s , l , i g g for isotope i, and from energy group g into energy group g, are tabulated parameters. The zeroth-order (i.e., l = 0 ) scattering cross sections must be considered separately from the higher order (i.e., l 1 ) scattering cross sections, since the former contribute to the total cross sections (as noted below), while the latter do not. Aiming at reducing the proliferation of superscripts and subscripts when defining response sensitivities with respect to the microscopic scattering cross sections σ s , l , i g g , these cross sections will be considered to be components of a vector σ s defined below:

σ s [ s 1 , , s J S X ] [ σ s , l = 0 , i = 1 g = 1 g = 1 , σ s , l = 0 , i = 1 g = 2 g = 1 , , σ s , l = 0 , i = 1 g = G g = 1 , σ s , l = 0 , i = 1 g = 1 g = 2 , σ s , l = 0 , i = 1 g = 2 g = 2 , , σ s , l , i g g , , σ s , I S C T , i = I G G ] , l = 0 , , I S C T ; i = 1 , , I ; g , g = 1 , , G ; J S X = ( G × G ) × I × ( I S C T + 1 ) . (97)

5) The total cross section for energy group g , g = 1 , , G is denoted as Σ t g ( α ; r ) and is computed using the following expression:

Σ t g ( α ; r ) = i = 1 I N i σ t , i g ; σ t , i g = [ σ f , i g ( r ) + σ c , i g ( r ) + g = 1 G σ s , l = 0 , i g g ( r ) ] . (98)

In Equation (98), the quantities σ t , i g , σ f , i g ( r ) and σ c , i g ( r ) denote, respectively, the total microscopic cross section, the tabulated group microscopic fission, and the neutron capture cross sections for isotope i and group g. Other nuclear reactions in the PERP benchmark are negligible. The total microscopic cross sections σ t , i g involve three indices, which will proliferate exponentially when determining the higher-order (up to and including the 4th-order) sensitivities of the PERP leakage response with respect to these cross sections. In order to reduce as much as possible the proliferation of indices, it is useful to consider that the cross sections σ t , i g are the components of a vector t , having T G × I components and defined as follows:

t [ t 1 , , t J T X ] [ σ t , i = 1 1 , σ t , i = 1 2 , , σ t , i = 1 G , , σ t , i g , , σ t , i = I 1 , , σ t , i = I G ] , for i = 1 , , I = 6 ; g = 1 , , G = 30 ; J T X I × G . (99)

6) PARTISN [5] computes the quantity ( ν Σ f ) g ( α ; r ) for each isotope i and energy group g, as follows:

( ν Σ f ) g ( α ; r ) = i = 1 N F = 2 N i σ f , i g ν i g , g = 1 , , G = 30 , (100)

where σ f , i g denotes the microscopic fission cross section for isotope i and energy group g, ν i g denotes the average number of neutrons per fission for isotope i and energy group g, and N F denotes the total number of fissionable isotopes.

σ f [ σ f , i = 1 1 , σ f , i = 1 2 , , σ f , i = 1 G , , σ f , i g , , σ f , i = N f 1 , , σ f , i = N F G ] [ f 1 , , f J F X ] , i = 1 , , N F ; g = 1 , , G ; J F X = G × N F ; (101)

ν [ ν i = 1 1 , ν i = 1 2 , , ν i = 1 G , , ν i g , , ν i = N f 1 , , ν i = N F G ] [ f J F X + 1 , , f J F X + J N U ] , i = 1 , , N F ; g = 1 , , G ; J N U = G × N F . (102)

7) The quantity χ g ( α ; r ) quantifies the fission spectrum in energy group g. The fission spectrum is considered to depend on the vector of parameters p , defined as follows:

p [ p 1 , , p J χ ] [ χ i = 1 g = 1 , χ i = 1 g = 2 , , χ i = 1 G , , χ i g , , χ N F G ] , for i = 1 , , N F