First-Order Comprehensive Adjoint Method for Computing Operator-Valued Response Sensitivities to Imprecisely Known Parameters, Internal Interfaces and Boundaries of Coupled Nonlinear Systems: II. Application to a Nuclear Reactor Heat Removal Benchmark

: This work illustrates the application of a comprehensive ﬁrst-order adjoint sensitivity analysis methodology (1st-CASAM) to a heat conduction and convection analytical benchmark problem which simulates heat removal from a nuclear reactor fuel rod. This analytical benchmark problem can be used to verify the accuracy of numerical solutions provided by software modeling heat transport and ﬂuid ﬂow systems. This illustrative heat transport benchmark shows that collocation methods require one adjoint computation for every collocation point while spectral expansion methods require one adjoint computation for each cardinal function appearing in the respective expansion when recursion relations cannot be developed between the corresponding adjoint functions. However, it is also shown that spectral methods are much more e ﬃ cient when recursion relations provided by orthogonal polynomials make it possible to develop recursion relations for computing the corresponding adjoint functions. When recursion relations cannot be developed for the adjoint functions, the collocation method is probably more e ﬃ cient than the spectral expansion method, since the sources for the corresponding adjoint systems are just Dirac delta functions (which makes the respective computation equivalent to the computation of a Green’s function), rather than the more elaborated sources involving high-order Fourier basis functions or orthogonal polynomials. For systems involving many independent variables, is likely that a hybrid combination of spectral expansions in some independent variables and collocation in the remaining independent variables would provide the most e ﬃ cient computational outcome.


Introduction
This work illustrates the application of First-Order Comprehensive Adjoint Sensitivity Analysis Methodology (1st-CASAM) presented in the companion work [1] to an analytical benchmark problem that models coupled heat conduction and convection in a nuclear reactor. The benchmark simulates the heat transport in an electrically heated cylindrical rod surrounded by a coolant which simulates the geometry of a nuclear reactor. As will be shown in this work, this benchmark problem admits exact closed-form solutions for the sensitivities of the temperature distribution in the coupled rod/coolant system, which can be used to benchmark thermal-hydraulics' production codes. In particular, this benchmark was used in [2] to verify the numerical results produced by the "ANSYS FLUENT Adjoint Solver" [3] for conditions typical of the Gen4Energy's Pb-Bi advanced reactor design [4], underscoring the current limitations of the FLUENT Adjoint Solver.
This work is structured as follows: Section 2 presents the mathematical modeling of the coupled conduction and convection heat transport benchmark problem. Section 3 presents the computation of the first-order sensitivities of the coolant's temperature: Section 3.1 presents the spectral expansion of the response sensitivities, while Section 3.2 presents the collocation pseudo-spectral expansion of the response sensitivities. Section 4 presents the computation of the first-order sensitivities of the rod temperature, paralleling the presentations in Section 3, i.e., Section 4.1 presents the spectral expansion of the response sensitivities while Section 4.2 presents the collocation pseudo-spectral expansion of the response sensitivities. Section 5 highlights the differences between the exact computation of first-order sensitivities for scalar-valued responses in phase-space versus the computation of the sensitivities of operator-valued responses. Section 6 offers concluding remarks, analyzing the computational issues and relative efficiencies and accuracies pertaining to the spectral, collocation, and mixed spectral/collocation expansion strategies for computing the response sensitivities with respect to the benchmark's imprecisely known model, interface, and boundary parameters.

Mathematical Description of the Paradigm Benchmark Heat Transport Problem
The benchmark problem presented in this Section models the steady-state heat conduction in the axial direction of an electrically heated rod, which simulates a fuel rod in an operating nuclear reactor. The geometrical characteristics of the electrically heated rod are: radius a and length (height) . The rod is heated by a heat source of the form q cos(πx/ ), where q W·m −3 denotes a constant volumetric source and x denotes the coordinate along the rod's axial (customarily, the vertical) direction. This heat source simulates the heat distribution in a nuclear reactor. The heated rod transfers heat by convection to surrounding coolant that flows along the rod's vertical direction. The rod conductivity, k W·m −1 K −1 , is considered to be a temperature-independent constant. The rod's surface is cooled by forced convection to a surrounding liquid flowing along the rod length, from the rod's lower end, taken to be located at x = − /2, towards the rod upper end, located at x = /2. The heat transfer coefficient, h W·m −2 K −1 , from the rod's surface to the coolant is considered to be constant. For this benchmark, the rod length is typically two orders of magnitude larger than its diameter, so the heat conduction process in the rod's axial direction can be neglected by comparison to the heat conduction in the rod's radial direction.
The responses of interest for this benchmark system are the steady-state temperature distributions within the heated rod and coolant (fluid), respectively. The fluid temperature, T f l (x), satisfies the following energy conservation equation where W kg·s −1 denotes the mass flow rate, T inlet K denotes the coolant inlet temperature, c p J·kg −1 K −1 denotes the coolant specific heat capacity. The temperature distribution within the electrically heated rod, denoted as T(r, x), satisfies the following energy conservation equations ∂T(r, x) ∂r = 0, at r = 0, − k ∂T(r, x) ∂r = h T(r, x) − T f l (x) , at r = a (5) The correspondences between the quantities appearing in Equations (1)- (5) and the general definitions provided in [1] for characterizing "Subsystem I" and "Subsystem II" are as follows 1.
The independent variable for "Subsystem I" is x The vector of independent variables for "Subsystem II" y (r, x) † , Z y = 2, with Note that the second component of the vector y (r, x) † of independent variables for "Subsystem II" coincides with the independent variable for "Subsystem I" and the domain of definition of x coincides for both "Subsystems I and II."; 4.
The state function (i.e., dependent variable) for "Subsystem I" is u( ; α] has a single component provided by the left-side of Equation (1) while the source term Q (I) (α; x) also has a single component, which is given by the right-side of Equation (1); 5.
The state function (i.e., dependent variable) for "Subsystem II" is v(y) T(r, x), Z v = 1. The operator N (II) [v(y); α] has a single component provided by the left-side of Equation (3) while the source term Q (II) (α; y) also has a single component, which is given by the left-side of Equation (3); 6.
The components of the vector B[u(x), v(y); α; x, y] are the boundary conditions provided by Equations (2) and (4), and the interface condition provided in Equation (5).
The paradigm benchmark problem defined by Equations (1)- (5) has been deliberately designed to admit exact closed-form solutions for its underlying state functions, which are as follows The expressions obtained in Equations (6) and (7) will be subsequently used for verification/validation of the sensitivities of the state functions T f l (z) and T(r, z) to all of the imprecisely known parameters.

First-Order Sensitivities of the Coolant's Temperature
The application of the 1st-CASAM presented in [1] will be illustrated in this section by computing the sensitivities of the coolant temperature distribution to the imprecisely known model, interface, and boundary parameters. The arbitrary variations δq, δk, δh, δW, δc p , δT inlet , δa, δ in the model parameters will cause a variation in the fluid temperature, denoted as δT f l (x), which is obtained by solving (see [1]) the First-Level Adjoint Sensitivity System (1st-LFSS) for the fluid temperature, which is obtained by G-differentiating Equations (1) and (2). By definition, the G-differentials of Equations (1) and (2) are obtained as follows Carrying out in Equations (8) and (9) the differentiations with respect to ε and setting ε = 0 in the resulting expressions yields the following set of equations, which comprise the 1st-LFSS for the fluid = 0, as evidenced by evaluating Equation (6) comprising Equations (10) and (11), evidently depends on the arbitrary parameter variations. The 1st-LFSS can be solved in closed form and its solution can be used for the subsequent verification of the expressions to be obtained for the 1st-order response sensitivities. The solution of Equations (10) and (11) has the following expression It is evident that the expression obtained in Equation (12) is the total differential with respect to the model and boundary parameters of the expression of T f l (x) given in Equation (6).

Spectral Expansion of the Response Sensitivities
The coolant temperature distribution T f l (x) is defined on the domain − /2 ≤ x ≤ /2, for which it is possible to use either conventional Fourier or orthogonal polynomials as spectral basis functions for representing T f l (x). Since orthogonal polynomials can be computed recursively, whereas the Fourier basis functions cannot, using orthogonal polynomials is often computationally advantageous. The use of such basis functions will be illustrated in this section by employing Legendre polynomials, since all expressions will be obtainable exactly in closed-form, for convenient verification. For the derivations to follow, it is convenient (but not mandatory) to shift the domain − 0 /2 ≤ x ≤ 0 /2, on which the 1st-LFSS and 1st-LASS (First-Level Adjoint Sensitivity System, cf. 1) are defined, to the standard domain of definition −1 ≤ t ≤ 1 for the Legendre Polynomials P n (t). This domain shift is accomplished by changing in Equations (10) and (11) the independent variable from x to t as follows On the domain −1 ≤ t ≤ 1, the Legendre-expansion of the coolant temperature distribution The total 1st-order sensitivity of T f l (x) is given by the 1st-order G-differential of the expression on the right-side of Equation (14), which is where ∂q P n (t)dt Implementing into Equations (10) and (11) the change in independent variables indicated in Equation (13) yields the following expression for the 1st-LFSS for δT f l (t) The Legendre coefficient 1 −1 δT f l (y) P n (y)dy in Equation (15) could be computed by repeatedly solving the 1st-LFSS for δT f l (t), for every possible parameter variation. Alternatively, the Legendre coefficient in Equation (15) can be computed by using the 1st-LASS corresponding to the 1st-LFSS represented by Equations (18) and (19), which will be constructed next by following the concepts outlined in Section 3.2. The first step is to form the inner product on −1 ≤ t ≤ 1 of Equation (18) with a square-integrable function Ψ n (t) to obtain the following relation: Integrating by parts the term on the left side of Equation (20) yields the following relation Identify the first term on the right-side of Equation (21) with the Legendre coefficient 1 −1 δT f l (y) P n (y)dy to obtain the following relations The term in Equation (23) which contains the unknown quantity δT f l (1) is set to zero by imposing the following boundary condition for the function Ψ n (t) Imposing the boundary condition required in Equation (24) transforms Equation (23) into the following form where the adjoint function Ψ n (t) is the solution of the 1st-LASS defined by Equations (22) and (24), and has the following expression Inserting the definition of Q f l (t) provided in Equation (18)  δT f l (y) P n (y)dy in terms of the adjoint function Ψ n (t) for n = 1, 2, . . .
The last term in each of the double equalities indicated in Equations (27)-(31) was obtained by integrating by parts the corresponding middle term, and subsequently using Equation (22).
The results obtained in Equations (27)-(32) are inserted in Equation (17) to obtain Inserting the expression obtained in Equation (33) into Equation (15) and using Equation (16) to identify the expressions that multiply the same parameter variation in both sides of Equation (15) yields the following expressions for the partial sensitivities of T f l (t) to the various model parameters Evidently, the expressions for the partial sensitivities obtained in Equations (34)-(39) are identical to the expressions that would be obtained by determining the partial sensitivities from the closed-form expression of T f l (x) provided in Equation (6). In practice, however, the closed-form expression of T f l (x) provided in Equation (6) is not available.
The structure of the illustrative fluid-flow system analyzed in the foregoing enabled the complete avoidance of performing either forward or adjoint large-scale computations (which typically occur when solving differential, integral, etc., systems). This fact illustrates the essential importance of choosing an optimal spectral expansion for the response and problem under consideration. For the illustrative fluid-flow system analyzed in this section, the rescaling of the problem by shifting its original domain − 0 /2 ≤ x ≤ 0 /2 to the standard domain of definition −1 ≤ t ≤ 1 for the Legendre Polynomials P n (t) was not actually necessary. This re-scaling was performed in order to illustrate the concepts that would be applied if a rescaling were required by the available software.
For the illustrative benchmark system analyzed in this section, the sensitivities of the spatial fluid temperature distribution to the model parameters could be computed without needing to solve either the 1st-LFSS or the 1st-LASS explicitly because the particular form of this benchmark was amenable to obtaining all of the results as exact closed-form expressions. Although such situations are unlikely to arise in practice when considering large-scale systems (for which closed-form solutions are unlikely to be available), it is nevertheless possible to massively reduce the number of large-scale computations needed for solving the 1st-LASS by using recursion relations available for orthogonal polynomials. For the illustrative benchmark presented in this section, for example, the adjoint functions Ψ n (t) in Equation (26) can be conveniently expressed in terms of Legendre polynomials by recalling [5,6] the following relation Integrating the relation given in Equation (40) and replacing the resulting expression in Equation (26) yields the following relation Thus, in this case, the 1st-LASS would need to be solved only once, to compute the function Ψ 0 (t). The functions Ψ n (t), n > 0, can subsequently be computed by using the recurrence relation provided in Equation (41), instead of solving the differential equations underlying the 1st-LASS. The plethora of recursion relations provided by orthogonal polynomials makes their use preferable over standard Fourier series. The computationally most intensive situation occurs when it is not possible to use recursion relations, but when all of the needed adjoint functions would need to be computed by solving the differential equations underlying the 1st-LASS.

Collocation Pseudo-Spectral Expansion of the Response Sensitivities
The collocation pseudo-expansion of the fluid temperature variation is where the quantity C i (x) denotes the ith chosen cardinal function and where The functional in Equation (43) can be expressed in terms of the solution of a 1st-LASS constructed by following the same conceptual steps as outlined in [1]. The first step is to form the inner product on − 0 /2 ≤ x ≤ 0 /2 of Equation (10) with a square-integrable function Ψ i (x) to obtain the following relation Integrating by parts the term on the left-side of Equation (44) yields the following relation Requiring the first term on the right-side of Equation (45) to represent with the functional in Equation (43) yields the following relation where the adjoint function satisfies the following equation The term in Equation (46) which contains the unknown quantity δT f l 0 /2 is set to zero by imposing the following boundary condition for the function Ψ i (x) Inserting the boundary condition given in Equation (48) into Equation (46) transforms the latter into the following form where the adjoint function Ψ i (x) is the solution of the 1st-LASS defined by Equations (47) and (48). Solving the 1st-LASS defined by Equations (47) and (48) yields the following explicit expression where H(x i − x) denotes the Heaviside functional: Recalling that inserting the definition of Q f l (x) provided in Equation (10) into Equation (49), equating the resulting expression with Equation (51), identifying the expressions that multiply the respective arbitrary parameter variations on both sides of the resulting equality and using the result provided in Equation (50), ultimately yields the following expressions for the sensitivities of T f l (x i ) to the respective parameters The results in Equations (52)-(57) confirm that the sensitivities are computed exactly at the colocation/interpolation points using the solution Ψ i (x) of the 1st-LASS.
Two of the most convenient grids for interpolating functions are the "Chebyshev-roots" grid and the "Chebyshev-extrema" ("Gauss-Lobatto") grid see, e.g., [5,6]. To illustrate the use of these grids for interpolating δT f l (x), it is convenient to map the domain − 0 /2 ≤ x ≤ 0 /2 onto the standard domain of definition −1 ≤ t ≤ 1 for the Chebyshev polynomials T n (x), using the transformation defined in Equation (13). On the standard interval −1 ≤ t ≤ 1, Chebyshev-root grid, the total sensitivity δT f l (t) is interpolated by a polynomial of the form where T n (t) denotes the n th -order Chebyshev polynomial defined on the standard interval −1 ≤ t ≤ 1.
On the Chebyshev-extrema grid, the total sensitivity δT f l (t) is interpolated by a polynomial of the form It is known [5,6] that the coefficients a n of the exact spectral representation of δT f l (t) in terms of Chebyshev polynomials T n (t), i.e., δT f l (t) = a 0 2 + ∞ n=1 a n T n (t), a n 2 π are related to the coefficients b n and c n through the following relations It is also known [5,6] that the errors in either of the interpolating polynomials is bounded by twice the sum of the absolute values of all the neglected coefficients, i.e., The Chebyshev-extrema and Chebyshev-root grids can be mapped from the standard interval −1 ≤ t ≤ 1 onto the domain − 0 /2 ≤ x ≤ 0 /2 by using the transformation provided in Equation (13).
Comparing the expressions for the sensitivities of the fluid temperature with respect to the model parameters obtained in Section 3.1 by using the full spectral expression to the expressions for the corresponding sensitivities obtained in Section 3.2 by using the collocation/pseudospectral expansion indicates that the latter (i.e., the collocation/pseudospectral expansion) is easier to implement in practice.

First-Order Sensitivities of the Rod's Temperature
To first-order in the parameter variations, the 1st-order variation δT(r, x) around the nominal temperature distribution, T 0 (r, x), in the electrically heated rod is the solution of the 1st-LFSS obtained by determining the G-differentials of Equations (3)- (5). By definition, the G-differentials of Equations (3)-(5) are obtained as follows d dε ∂(a 0 +εδa) ε=0 Carrying out in Equations (63)-(65) the differentiations with respect to ε and setting ε = 0 in the resulting expressions yields the following set of equations comprising 1st-LFSS The first term on the right-side of Equation (66) can be simplified by using Equation (3) to obtain the following equation The terms containing derivatives of T 0 (r, x) in Equation (68) can also be simplified using Equations (3) and (5) to obtain the following equation The interface condition expressed by Equation (70) couples the variation in the rod temperature δT(r, x), to the variation in the fluid temperature, δT f l (x), thereby coupling Equations (67), (69) and (70) to Equations (10) and (11), the totality of which constitute the 1st-LFSS. As has been discussed throughout this work, it is computationally expensive to repeatedly solve the 1st-LFSS, namely Equations (10), (11), (67), (69) and (70), for all possible parameter variations.

Spectral Expansion of the Response Sensitivities
In this Section, the 1st-LASS corresponding to the 1st-LFSS comprising Equations (10), (11), (67), (69) and (70) will be developed by applying the general principles outlined in [1], to represent the temperature distribution variation, δT(r, x), in the heated rod by means of a spectral expansion. As in Section 3.1, it is convenient to use the Legendre polynomials, P n (x) as the spectral basis functions for representing the function δT(r, x) in the x-direction. It is therefore convenient to rescale Equations (67), (69) and (70) in the x-direction by using the scaling transformation indicated in Equation (13) the 1st-LFSS. Consequently, the re-scaled 1st-LFSS for the coupled variations δT(r, x) and δT f l (x) are as follows For convenient reference and completeness of the 1st-LFSS, Equations (18) and (19) were reproduced above as Equations (74) and (75).
The recommended [5,6] spectral basis functions for representing a continuous function like δT(r, x) on the finite interval 0 ≤ r ≤ a 0 would be Chebyshev Polynomials. However, their use would be very similar to the use of the Legendre Polynomials as the spectral basis functions for the expansion in the axial direction −1 ≤ t ≤ 1. Consequently, to illustrate the use of a spectral basis for which there are no recursion relations among the underlying spectral elements, Fourier series expansions will be used in the radial direction for representing δT(r, x). Since δT(r, x) is defined on the half-range 0 ≤ r ≤ a 0 , either Fourier Sine or Cosine Series can be employed to represent δT(r, x) in the radial direction. Since the result in Equation (7) indicates that δT(r, x) is an even function in r and since the Cosine Series converges faster than the Sine Series, the former will be used. Thus, δT(r, x) can be represented over the domain 0 ≤ r ≤ a 0 ∪ [−1 ≤ t ≤ 1] using the following spectral representation where the generalized Fourier spectral coefficients F mn (δT) are defined as follows δT(r, t)P n (t)dt , m ≥ 1, n ≥ 0.
The appearance of the unknown temperature variation δT(r, t) in the expressions of the generalized Fourier spectral coefficients F mn (δT), m ≥ 0, n ≥ 0, will be eliminated by constructing equivalent expressions for these coefficients in terms of adjoint functions, which will be the solutions of the 1st-LASS corresponding to the 1st-LFSS comprising Equations (71)-(75). The Hilbert space appropriate for constructing this 1st-LASS comprises square-integrable two-component vector functions of the form u(x) [u 1 (r, z), u 2 (z)] † , endowed with an inner product u(x), ψ(x) of the form Using the definition provided in Equation (79), construct the inner product of a square integrable vector function ψ(x) = [Ψ mn (r, t), Φ mn (t)] with Equations (71) and (74), respectively, to obtain the following relation The left-side of Equation (80) is now integrated by parts (twice over the variable r and once over the variable t) to obtain the following relation Using the boundary condition given in Equation (72) The unknown quantity ∂[δT(r, t)]/∂r r=a 0 , which appears in the last term on the right-side of Equation (84) is eliminated by using the boundary condition given in Equation (73); this operation transforms Equation (84) into the following form The unknown quantity δT a 0 , t , which appears in third and fourth terms on the right-side of Equation (85) is eliminated by imposing the following interface condition on the (adjoint) function Ψ mn (r, t) k 0 ∂Ψ mn (r, t) ∂r + h 0 Ψ mn (r, t) = 0, at r = a 0 (86) Inserting Equation (86) into the right-side of Equation (85) reduces the latter equation to the following form ∂r (87) The two terms that contain the unknown function δT f l (z) in Equation (87) are grouped together, transforming Equation (87) into the following form The second-term on the right-side of Equation (88) will represent the generalized Fourier spectral coefficients F 0n (δT) defined in Equation (77) by requiring that the following equations be satisfied Altogether, the relations required in Equations (89) and (90), and the boundary and interface conditions already imposed in Equations (82), (83), and (86) for m = 0 constitute the 1st-LASS for the adjoint functions Ψ 0n (r, t), m = 0, n ≥ 0. Solving this 1st-LASS yields the following expressions for Ψ 0n (r, t), m = 0, n ≥ 0 Inserting Equations (89) and (90) into Equation (88) and re-arranging the resulting relation yields the following expression for the generalized Fourier spectral coefficients F 0n (δT) Inserting the definitions provided for Q(t) and Q f l (t) in Equations (71) and (74), respectively, into Equation (92) yields the following expression for the generalized Fourier spectral coefficients F 0n (δT), m = 0, n ≥ 0 The second-term on the right-side of Equation (88) also represents the generalized Fourier spectral coefficients F mn (δT), m ≥ 1, n ≥ 0, defined in Equation (78) by requiring that the following equations be satisfied Inserting Equations (94) and (95) into Equation (88) and re-arranging the resulting relation yields the following expression for the generalized Fourier spectral coefficients F mn (δT), m ≥ 1, n ≥ 0 Inserting the definitions provided for Q(t) and Q f l (t) in Equations (71) and (74), respectively, into Equation (96) yields the following expression for the generalized Fourier spectral coefficients F mn (δT), m ≥ 1, n ≥ 0 Altogether, the relations required in Equations (94) and (95), and the boundary and interface conditions already imposed in Equations (82), (83), and (86) for m ≥ 1 constitute the 1st-LASS for the adjoint functions Ψ mn (r, t), m ≥ 1, n ≥ 0. Solving this 1st-LASS yields the following expressions for Ψ mn (r, t), m = 0, n ≥ 0 Noting from Equation (98) that inserting Equations (93) and (97) into Equation (76) and equating the expressions that multiply each of the respective arbitrary parameter variation yields the following expressions for the sensitivities of the rod temperature T(r, t) with respect to the various parameters ∂T(r,t) ∂q ∂T(r,t) ∂k Inserting the results provided in Equations (91) and (98) in Equations (100)-(107), noting that and using the relations provided in Equation (13) to revert from the variable −1 ≤ t ≤ 1 to the variable − 0 /2 ≤ x ≤ 0 /2 yields, after considerable algebra, the following final exact expressions for the sensitivities of the temperature distribution T(r, t) within the heated rod to the model and boundary parameters Notably, the 1st-LASS is solved in a manner which is "reverse/backwards" by comparison to the way in which solution proceeds for solving the First-Level Forward Sensitivity System (1st-LFSS) as well as the original heat transport model. Thus, while the 1st-LFSS and the original heat transport model are solved by starting with the fluid flow equation (which is solved from the inlet to the outlet of the fluid flow) and subsequently solving the heat conduction equation in the rod, the solution of the 1st-LASS proceeds in the reverse manner, by first solving the heat conduction in the rod, followed by solving the fluid flow equation from the outlet to the inlet.
In practice, it is evidently not possible to infinitely compute many basis functions and corresponding expansion coefficients. Consequently, the exact expressions for the response sensitivities presented in Equations (109)-(116) cannot be attained. In practice, known convergence criteria for the various expansions see, e.g., [5,6] would be used in order to decide the number of expansion coefficients to be computed by solving the 1st-LASS. It is important to note that the issue of computational accuracy regarding the spectral expansion refers not to the accuracy of the functional derivative of the response in the phase-space of model parameters (in which space the 1st-CASAM provides exact expressions), but refers to the representation of the respective sensitivity as a function in the phase-space of independent variables.

Collocation Pseudo-Spectral Expansion of the Response Sensitivities
The collocation pseudo-spectral expansion of the fluid temperature variation δT(r, x) can be written in the following form where U i (r) and V j (x) denotes the chosen cardinal functions and where In Equation (118), the quantities r i , i = 0, 1, . . . , I denote the collocation points in the radial direction while the quantities x j , j = 0, 1, . . . , J denote the collocation points in the axial direction. The functional δT r i , x j in Equation (118) can be expressed in terms of the solution of a 1st-LASS that is constructed by following the same conceptual steps as those leading to Equation (82), (83), (86), (94) and (95). Thus, denoting as Ψ r, z; r i , x j and Ψ f l z; x j the adjoint sensitivity functions that correspond to the forward functions δT(r, z) and δT f l (z), respectively, and following the procedure outlined in [1] leads to the following expressions for the sensitivities of T r i , x j in terms of the adjoint functions Ψ(r, z) and Ψ f l (z) In Equations (119)-(126), the adjoint functions Ψ r, z; r i , x j and Ψ f l z; x j are the solutions of the following 1st-LASS k 0 ∂Ψ r, z; r i , x j ∂r + h 0 Ψ r, z; r i , x j = 0, at r = a 0 (129) Solving Equations (127)-(131) yields the following closed-form expressions for the adjoint functions Ψ r, z; r i , x j and Ψ f l z; x j Ψ r, z; r i , where H(z) denotes the customary Heaviside functional. i.e., H(z) = 1 , i f z > 0 and H(z) = 0 , i f z < 0. Inserting the expressions obtained in Equations (132) and (133) into Equations (119)-(126) yields the following expressions for the response sensitivities The 1st-LASS is solved in a manner that is "reverse/backwards" by comparison to the way in which solution proceeds for solving the 1st-LFSS as well as the original heat transport model, by first solving the heat conduction in the rod, followed by solving the fluid flow equation from the outlet to the inlet.

Mixed Spectral/Collocation Expansion of the Response Sensitivities
In contrast to δT f l (x), the total sensitivity δT(r, x) of the rod temperature (with respect to the model parameters) depends on more than one independent variable. Therefore, a mixed spectral/collocation representation of δT(r, t) over the domain 0 ≤ r ≤ a 0 ∪ [−1 ≤ t ≤ 1] may be contemplated, as provided below where P n (t) denotes the nth-order Legendre polynomial, U m (r) denotes the mth-cardinal function chosen for the radial direction and the functionals G in (δT) denote the Legendre spectral coefficients defined at collocation points r i , i = 0, 1, . . . , I, in the radial direction, as follows The 1st-LFSS appropriate for expressing the functionals G in (δT) in terms of adjoint functions is constructed by following the same conceptual steps as followed in Sections 3.1 and 3.2. Denoting the adjoint functions that would correspond to G in (δT) as Γ in (r, t) and Λ in (t), respectively, and following the same steps as those leading to Equation (88) yields the following counterpart of Equation (88) The second-term on the right-side of Equation (144) will represent the coefficients G in (δT) defined in Equation (143) by requiring that the following equations, which represent the 1st-LASS for the adjoint functions Γ in (r, t; r i ) and Λ in (t), be satisfied inserting Equations (93) and (97) into Equation (76) and equating the expressions that multiply each of the respective arbitrary parameter variation yields the following expressions for the sensitivities of the rod temperature T(r, t) with respect to the various parameters Solving the 1st-LASS represented by Equations (145)-(149) yields the following expressions for the adjoint functions Γ in (r, t; r i ) and Λ in (t) Inserting the results from Equations (158) and (159) into Equations (150)-(157), carrying out the respective algebraic operations and using Equation (13) to revert from the independent variable −1 ≤ t ≤ 1 to the independent variable − 0 /2 ≤ x ≤ 0 /2 yields the following expressions for the sensitivities of the time-dependent temperature T(r i , t) at a collocation point r i within the heated rod to the model and boundary parameters Inserting the expressions from Equation (181) into Equation (180), carrying out the respective algebraic operations and identifying the expressions that multiply each of the arbitrary parameter variations with the respective partial sensitivity of the average temperature in the rod, T ave , with respect to the corresponding model parameter, yields the following expressions generalized Fourier coefficients which would be needed to represent the response within an a priori established accuracy in the phase-space of independent variables. Subsequently, for each generalized Fourier coefficient and/or at each collocation point, the 1st-CASAM most efficiently computes the exact response sensitivities in the parameter-space.
The illustrative heat conduction and convection benchmark problem presented in this work has also indicated that the spectral expansion computations are much more efficient when the basis-functions are orthogonal polynomials possessing recursion relations, which could enable the development of recursion relations for computing the corresponding adjoint functions. Hence, instead of solving the 1st-LASS anew for each spectral basis function, the recursion relations developed for the adjoint functions could be efficiently used for determining each of the adjoint function that corresponds to the respective orthogonal polynomial (spectral basis function) appearing in the source of the 1st-LASS. On the other hand, when recursion relations cannot be developed for the adjoint functions, the collocation method is probably more efficient than the spectral expansion method, since the sources for the corresponding adjoint systems are just Dirac delta functions (which makes the respective computation equivalent to the computation of a Green's function) rather than the more elaborated sources for the spectral expansion method, which involve high-order Fourier basis functions or orthogonal polynomials.
For systems involving many independent variables, it is likely that a hybrid combination of spectral expansions in some independent variables and collocation in the remaining independent variables would offer the most efficient computational outcome, as has been illustrated in this work.
All in all, it is important to note that all of the approximations introduced by using either spectral or collocation methods relate to the representation of the model response's sensitivities (with respect to the model's parameters) in the phase-space of the model's independent variables (as opposed to the phase-space of the model's parameters, since the 1st-CASAM provides exact first-order sensitivities of a scaler-valued model response with respect to the model's parameters).
Numerical examples have deliberately not been provided in this work, since the numerical examples would (by their very nature) have been specific (rather than general), and would have therefore not served the purpose of this work, which is the presentation of a generally usable analytical benchmark for quantifying the accuracy of heat transfer codes. The interested readers may find numerical sensitivity analysis results (sans sensitivities to the phase-space locations of boundaries and interfaces) in [2], for the specific case of a "Generation-IV" advanced Pb-Bi reactor. Of course, the rankings of the various response sensitivities depend on the specific systems being analyzed.
By enabling the exact computations of operator-valued response sensitivities to internal interfaces and external boundary parameters and conditions, the 1st-CASAM presented in this work makes it possible, inter alia, to quantify the effects of manufacturing tolerances on the responses of physical and engineering systems. Ongoing research will generalize the methodology presented in this work, aiming at computing exactly and efficiently the second-order sensitivities of operator-valued responses for coupled nonlinear systems with respect to the systems' imprecisely known parameters, internal and external boundaries.