The First-Order Comprehensive Sensitivity Analysis Methodology (1st-CASAM) for Scalar-Valued Responses: I. Theory

Dan Gabriel Cacuci

Department of Mechanical Engineering, University of South Carolina, Columbia, SC, USA.

DOI: 10.4236/ajcm.2020.102015                

Abstract

This work presents the first-order comprehensive adjoint sensitivity analysis methodology (1st-CASAM) for computing efficiently, exactly, and exhaustively, the first-order sensitivities of scalar-valued responses (results of interest) of coupled nonlinear physical systems characterized by imprecisely known model parameters, boundaries and interfaces between the coupled systems. The 1st-CASAM highlights the conclusion that response sensitivities to the imprecisely known domain boundaries and interfaces can arise both from the definition of the system’s response as well as from the equations, interfaces and boundary conditions defining the model and its imprecisely known domain. By enabling, in premiere, the exact computations of sensitivities to interface and boundary parameters and conditions, the 1st-CASAM enables the quantification of 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 higher-order response sensitivities for coupled systems involving imprecisely known interfaces, parameters, and boundaries.

Keywords

Adjoint Sensitivity Analysis (1st-CASAM), Response Sensitivities for Coupled Nonlinear Systems, Imprecisely Known Interfaces, Imprecisely Known Parameters, Imprecisely Known Boundaries

1. Introduction

Many works have been published on using adjoint operators for computing first- and second-order sensitivities (i.e., functional derivatives) of model responses (i.e., results produced by models) to imprecisely known model parameters since the original work of Wigner [1] on the linear neutron transport equation and the introduction of the first-order adjoint sensitivity analysis methodology for nonlinear systems by Cacuci [2] [3]. Representative works in this regard are cited in the books by Cacuci [4] [5], along with the original presentations of the first- and the second-order adjoint sensitivities analysis methodologies. It is well known that the adjoint method of sensitivity analysis [2] [3] [4] [5] enables the most efficient computation of the exact (to machine or to a priori set precision) response sensitivities to model parameters. The efficiency of the second-order adjoint sensitivity analysis methodology developed by Cacuci [5] has been recently demonstrated by its application to a OECD/NEA reactor physics benchmark [6] to compute [7] - [12] exactly the 21,976 first-order sensitivities and 482,944,576 second-order sensitivities of this benchmark’s response with respect to the benchmark’s model parameters, showing in particular that the effects of the 2^nd-order sensitivities on the uncertainty in the model’s response are even more important than the effects of the 1^st-order ones. Another step towards overcoming the curse of dimensionality in sensitivity analysis, uncertainty quantification and predictive modeling has been provided by the third-order adjoint sensitivity analysis methodology for linear systems provided recently by Cacuci [13].

However, none of the works cited above are capable of computing response sensitivities to imprecisely known domain internal and/or external boundaries. Very few works have attempted to develop mathematical/computational methodologies for computing exactly the first-order sensitivities of responses to imprecisely known boundaries. The representative works (Komata [14], Larsen and Pomraning [15], Rahnema and Pomraning [16], McKinley and Rahnema [17], Favorite and Gonzalez [18]) that have addressed this issue were limited to specific linear neutron transport or diffusion problems. Furthermore, none of the works published thus far have addressed, in a general theoretical/mathematical setting, the simultaneous computation of response sensitivities to imprecisely known model parameters, imprecisely known internal boundaries/interfaces between nonlinear systems that model coupled yet distinct physical processes, and/or imprecisely known external boundaries.

This work presents the mathematical foundations of a new method for computing efficiently, exactly and exhaustively, the first-order response sensitivities for coupled nonlinear physical systems characterized by imprecisely known parameters that describe not only processes within the system but also at the physical interfaces between systems, as well as at the systems’ imprecisely known domain boundaries. This new method will be called the first-order comprehensive adjoint sensitivity analysis methodology (1st-CASAM). Notably, the 1st-CASAM enables the quantification of the effects of manufacturing tolerances on the responses of physical and engineering systems.

This work is structured as follows: Section 2 presents the mathematical framework of two coupled generic nonlinear physical systems comprising imprecisely known parameters, internal interfaces, and external boundaries. Section 3 presents the mathematical framework of the 1st-CASAM, which enables the efficient computation of the exact sensitivities of a scalar-valued response with respect to the imprecisely known parameters, interfaces, and boundaries that characterize the generic coupled nonlinear physical systems. As is well known [19], the availability of response sensitivities to imprecisely known parameters, interfaces and boundaries is essential for a variety of subsequent uses, including uncertainty quantification, optimization, data assimilation, model calibration and validation, and reduction of uncertainties in predicted model results. Section 4 offers concluding remarks.

The sequel to this work [20] presents an illustrative application of the 1st-CASAM to a benchmark problem [21] [22] [23] that models coupled heat conduction and convection in a physical system comprising an electrically heated rod surrounded by a coolant which simulates the geometry of an advanced (“Generation-IV”) nuclear reactor [24]. This benchmark problem [21] [22] [23] 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 [21] [22] [23] was used to verify the numerical results produced by the FLUENT Adjoint Solver [25], showing that that the current “FLUENT Adjoint Solver” cannot compute any sensitivities for the temperature distribution within the solid rod. Although the “FLUENT Adjoint Solver” is capable of computing sensitivities of fluid temperatures to boundary parameters (e.g., boundary temperature, boundary velocity, boundary pressure), it yields accurate results only for the sensitivities of the fluid outlet temperature and the maximum rod surface temperature to the inlet temperature and inlet velocity, respectively.

2. Mathematical Modeling of Generic Coupled Nonlinear Physical Systems Comprising Imprecisely Known Parameters, Interfaces and Boundaries

The physical system considered in this work comprises two nonlinear physical systems which are coupled to one another across a common internal interface (boundary) in phase-space. Each system comprises imprecisely known model parameters, including imprecisely known parameters that characterize the interface between the systems and the systems’ outer boundaries. The first physical system is represented mathematically as follows:

$A\left[u\left(x\right);\alpha \right]=Q^{\left(A\right)}\left(\alpha ;x\right),x\in {\Omega }_x$ (1)

Bold letters will be used in this work to denote matrices and vectors. Unless explicitly stated otherwise, the vectors in this work are considered to be column vectors. The second system is represented mathematically as follows:

$B\left[v\left(y\right);\alpha \right]=Q^{\left(B\right)}\left(\alpha ;y\right),y\in {\Omega }_y$ (2)

If differential operators appear in Equations (1) and (2), a corresponding set of boundary and/or initial/final conditions must also be given; these conditions can be represented in operator form as follows:

$C\left[u\left(x\right),v\left(y\right);\alpha ;x,y\right]=0,x\in \delta {\Omega }_x,y\in \delta {\Omega }_y$ (3)

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

1) $\alpha \triangleq {\left({\alpha }_1,\cdots ,{\alpha }_{Z_{\alpha }}\right)}^{†}\in {ℝ}^{Z_{\alpha }}$ denotes a $Z_{\alpha }$ -dimensional column vector whose scalar-valued components are all of the imprecisely known internal and boundary parameters (both of) the physical systems, including imprecisely known parameters that characterize the interface and boundary conditions. Some of these parameters will be common to both physical systems, particularly those that characterize common interfaces. These scalar parameters are considered to be imperfectly known, subject to uncertainties. The minimum information needed for these parameters is their nominal or average values, which will be denoted as ${\alpha }^0\triangleq {\left({\alpha }_1^0,\cdots ,{\alpha }_{Z_{\alpha }}^0\right)}^{†}$. The superscript “zero” will be used in this work to denote known nominal or average values of various quantities. The symbol “ $\triangleq$ ” will be used to denote “is defined as” or “is by definition equal to” and transposition will be indicated by a dagger ( $†$ ) superscript.

2) $x\triangleq {\left(x_1,\cdots ,x_{Z_x}\right)}^{†}\in {ℝ}^{Z_x}$ denotes the $Z_x$ -dimensional phase-space position vector of independent variables for the system defined in Equation (1). The vector of independent variables $x$ is defined on a phase-space domain denoted as ${\Omega }_x$ which is defined as ${\Omega }_x\triangleq \left\{-\infty \le a_i\left(\alpha \right)\le x_i\le b_i\left(\alpha \right)\le \infty ;i=1,\cdots ,Z_x\right\}$. The lower-valued imprecisely known boundary-point of the independent variable $x_i$ is denoted as $a_i\left(\alpha \right)$, while the upper-valued imprecisely known boundary-point of the independent variable $x_i$ is denoted as $b_i\left(\alpha \right)$. For physical systems modeled by diffusion theory, for example, the “vacuum boundary condition” requires that the particle flux vanish at the “extrapolated boundary” of the spatial domain facing the vacuum; the “extrapolated boundary” depends on the imprecisely known geometrical dimensions of the system’s domain in space and also on the system’s microscopic transport cross sections and atomic number densities. The boundary $\partial {\Omega }_x\triangleq \left\{a\left(\alpha \right)\cup b\left(\alpha \right)\right\}$ of the domain ${\Omega }_x$ comprises all of the endpoints $a\left(\alpha \right)\triangleq {\left[a_1\left(\alpha \right),\cdots ,a_{Z_x}\left(\alpha \right)\right]}^{†}$ and $b\left(\alpha \right)\triangleq {\left[b_1\left(\alpha \right),\cdots ,b_{Z_x}\left(\alpha \right)\right]}^{†}$ of the intervals on which the respective components of $x$ are defined. It may happen that some components $a_i\left(\alpha \right)$ and/or $b_j\left(\alpha \right)$ are infinite, in which case they would not depend on any imprecisely known parameters.

3) $y\triangleq {\left(y_1,\cdots ,y_{Z_y}\right)}^{†}\in {ℝ}^{Z_y}$ denotes the $Z_y$ -dimensional phase-space position vector of independent variables for the physical system defined in Equation . The vector of independent variables $y$ is defined on a phase-space domain denoted as ${\Omega }_y$ which is defined as follows:

${\Omega }_y\triangleq \left\{-\infty \le c_j\left(\alpha \right)\le y_j\le d_j\left(\alpha \right)\le \infty ;j=1,\cdots ,Z_y\right\}$. The lower-valued imprecisely known boundary-point of the independent variable $y_j$ is denoted as $c_j\left(\alpha \right)$, while the upper-valued imprecisely known boundary-point of the independent variable $y_j$ is denoted as $d_j\left(\alpha \right)$. Some or all of the points $c_j\left(\alpha \right)$ may coincide with the points $b_j\left(\alpha \right)$. Also, some components of $y$ may coincide with some components of $x$.

4) $u\left(x\right)\triangleq {\left[u_1\left(x\right),\cdots ,u_{Z_u}\left(x\right)\right]}^{†}$ denotes a $Z_u$ -dimensional column vector whose components represent the system’s dependent variables (also called “state functions”). The vector-valued function $u\left(x\right)$ is considered the unique nontrivial solution of the physical problem described by Equations (1) and (2).

5) $v\left(y\right)\triangleq {\left[v_1\left(y\right),\cdots ,v_{Z_v}\left(y\right)\right]}^{†}$ denotes a $Z_v$ -dimensional column vector whose components represent the system’s dependent variables (also called “state functions”); The vector-valued function $v\left(y\right)$ is considered the unique nontrivial solution of the physical problem described by Equations (2) and (3).

6) $A\left[u\left(x\right);\alpha \right]\triangleq {\left[A_1\left(u;\alpha \right),\cdots ,A_i\left(u;\alpha \right),\cdots ,A_{Z_u}\left(u;\alpha \right)\right]}^{†},i=1,\cdots ,Z_u$ denotes a column vector of dimensions $Z_u$ whose components are operators (including differential, difference, integral, distributions, and/or infinite matrices) acting nonlinearly on $u\left(x\right)$ and $\alpha$.

7) $B\left[v\left(y\right);\alpha \right]\triangleq {\left[B_1\left(v;\alpha \right),\cdots ,B_i\left(v;\alpha \right),\cdots ,B_{Z_v}\left(v;\alpha \right)\right]}^{†},i=1,\cdots ,Z_v$ denotes a column vector of dimensions $Z_v$ whose components are operators (including differential, difference, integral, distributions, and/or infinite matrices) acting nonlinearly on $v\left(y\right)$ and $\alpha$.

8) $Q^{\left(A\right)}\left(\alpha ;x\right)\triangleq {\left[Q_1^{\left(A\right)}\left(\alpha ;x\right),\cdots ,Q_{Z_u}^{\left(A\right)}\left(\alpha ;x\right)\right]}^{†}$ denotes a $Z_u$ -dimensional column vector whose elements represent inhomogeneous source terms that depend either linearly or nonlinearly on $\alpha$. The components of $Q^{\left(A\right)}\left(\alpha ;x\right)$ may involve operators, rather than just finite-dimensional functions, and distributions acting on $\alpha$ and $x$.

9) $Q^{\left(B\right)}\left(\alpha ;y\right)\triangleq {\left[Q_1^{\left(B\right)}\left(\alpha ;y\right),\cdots ,Q_{Z_v}^{\left(B\right)}\left(\alpha ;y\right)\right]}^{†}$ denotes a $Z_v$ -dimensional column vector whose elements represent inhomogeneous source terms that depend either linearly or nonlinearly on $\alpha$. The components of $Q^{\left(B\right)}\left(\alpha ;y\right)$ may involve operators, rather than just finite-dimensional functions, and distributions acting on $\alpha$ and $y$.

10) The vector-valued operator $C\left[u\left(x\right),v\left(y\right);\alpha ;x,y\right]$ comprises all of the boundary, interface, and initial/final conditions for the coupled physical systems. If the boundary, interface and/or initial/final conditions are inhomogeneous, which is most often the case, then $C\left[0,0;\alpha ;x,y\right]\ne 0$.

11) Since $Q^{\left(A\right)}\left(\alpha ;x\right)$ and may involve operators and distributions acting on $\alpha$ and $y$, all of the equalities in this work, including Equations (1)-(3), are considered to hold in the weak (“distributional”) sense, since the right-sides (“sources”) of and of other various equations to be derived in this work may contain distributions (“generalized functions/functionals”), particularly Dirac-distributions and derivatives and/or integrals thereof.

The nominal solutions of Equations (1)-(3) will be denoted as $u^0\left(x\right)$ and $v^0\left(y\right)$ ; they are obtained by solving these equations at the nominal parameter values ${\alpha }^0$. In other words, the vectors $u^0\left(x\right)$ and $v^0\left(y\right)$ satisfy the following equations:

$A\left[u^0\left(x\right);{\alpha }^0\right]=Q^{\left(A\right)}\left({\alpha }^0;x\right),x\in {\Omega }_x$ (4)

$B\left[v^0\left(y\right);{\alpha }^0\right]=Q^{\left(B\right)}\left({\alpha }^0;y\right),y\in {\Omega }_y$ (5)

$C\left[u^0\left(x\right),v^0\left(y\right);{\alpha }^0;x,y\right]=0,x\in \delta {\Omega }_x,y\in \delta {\Omega }_y$ (6)

Equations (4)-(6) represent the “base-case” or nominal state of the physical system. Throughout this work, the superscript “0” will be used to denote “nominal” or “expected” values.

The response considered in this work is a generic scalar-valued operator (i.e., a functional) of the state functions, denoted as follows:

$R\left[u\left(x\right),v\left(y\right);\alpha ;x,y\right]$. (7)

The nominal value of the response, denoted as $R^0\triangleq R\left[u^0\left(x\right),v^0\left(y\right);{\alpha }^0;x,y\right]$, is determined by computing the response at the nominal values ${\alpha }^0$, $u^0\left(x\right)$ and $v^0\left(y\right)$.

3. Mathematical Framework of the 1st-CASAM for Operator-Valued Responses for Coupled Linear Physical Systems Comprising Imprecisely Known Parameters, Interfaces and Boundaries

As has been mentioned in the foregoing, the model and boundary parameters are considered to be imprecisely known quantities. Their true values may differ from their nominal (average, or “base-case”) values by variations denoted as $\delta \alpha \triangleq \left(\delta {\alpha }_1,\cdots ,\delta {\alpha }_{N_{\alpha }}\right)$, where $\delta {\alpha }_i\triangleq {\alpha }_i-{\alpha }_i^0$, $i=1,\cdots ,N_{\alpha }$. In turn, the parameter variations $\delta \alpha$ will cause variations $\delta u\left(x\right)\triangleq {\left[\delta u_1\left(x\right),\cdots ,\delta u_{Z_u}\left(x\right)\right]}^{†}$ and $\delta v\left(y\right)\triangleq {\left[\delta v_1\left(y\right),\cdots ,\delta v_{Z_v}\left(y\right)\right]}^{†}$ in the state functions, through Equations (1)-(3). Furthermore, the variations $\delta \alpha$, $\delta u\left(x\right)$ and $\delta v\left(y\right)$ will cause variations in the response $R\left[u\left(x\right),v\left(y\right);\alpha ;x,y\right]$ around the nominal response value $R^0$. Sensitivity analysis aims at computing the functional derivatives (called “sensitivities”) of the response to the imprecisely known parameters $\alpha$. Subsequently, these sensitivities can be used for a variety of purposes, including quantifying the uncertainties induced in responses by the uncertainties in the model and boundary parameters, combining the uncertainties in computed responses with uncertainties in measured response (“data assimilation”) to obtain more accurate predictions of responses and/or parameters (“model calibration,” “predictive modeling”, etc.). As has been shown by Cacuci [2] [3], the most general definition of the 1^st-order total sensitivity of an operator-valued model response to parameter variations is provided by the first-order “Gateaux-variation” (G-variation) of the response under consideration. To determine the first G-variation of the response $R\left[u\left(x\right),v\left(y\right);\alpha ;x,y\right]$, it is convenient to denote the functions appearing in the argument of the response as being the components of a vector $e\triangleq {\left[u\left(x\right),v\left(y\right);\alpha \right]}^{†}$, which represents an arbitrary “point” in the combined phase-space of the state functions and (all) parameters. The point which corresponds to the nominal values of the state functions and parameters in this phase space is denoted as $e^0\triangleq {\left[u^0\left(x\right),v^0\left(y\right);{\alpha }^0\right]}^{†}$. Analogously, it is convenient to consider the variations in the model’s state functions and parameters to be the components of a “vector of variations”, $\delta e$, defined as follows: $\delta e\triangleq {\left[\delta u\left(x\right),\delta v\left(y\right);\delta \alpha \right]}^{†}$. The 1^st-order Gateaux- (G-) variation of the response $R\left(e\right)$, which will be denoted as $\delta R\left(e^0;\delta e\right)$, for arbitrary variations $\delta e$ in the model parameters and state functions in a neighborhood ( $e^0+\epsilon \delta e$ ) around $e^0$, is obtained, by definition, as follows:

$\delta R\left(e^0;\delta e\right)\triangleq {\left\{\frac{\text{d}}{\text{d}\epsilon }R\left[u^0\left(x\right)+\epsilon \delta u\left(x\right),v^0\left(y\right)+\epsilon \delta v\left(y\right);{\alpha }^0+\epsilon \delta \alpha ;x,y\right]\right\}}_{\epsilon =0}$ (8)

The existence of the G-variation $\delta R\left(e^0;\delta e\right)$ does not guarantee its numerical computability. Numerical methods most often require that $\delta R\left(e^0;\delta e\right)$ be linear in the variations $\delta e$ in a neighborhood ( $e^0+\epsilon \delta e$ ) around $e^0$. The necessary and sufficient conditions for the G-differential $\delta R\left(e^0;\delta e\right)$ of a nonlinear operator $R\left(e\right)$ to be linear in $\delta e$ in a neighborhood ( $e^0+\epsilon \delta e$ ) around $e^0$, and thus admit partial and total G-derivatives, are as follows:

1) $R\left(e\right)$ satisfies a weak Lipschitz condition at $e^0$ ; (9)

2) for two arbitrary vectors of variations $\delta e_1$ and $\delta e_2$, the operator $R\left(e\right)$ satisfies the relation

$R\left(e^0+\epsilon \delta e_1+\epsilon \delta e_2\right)-R\left(e^0+\epsilon \delta e_1\right)-R\left(e^0+\epsilon \delta e_2\right)+R\left(e^0\right)=o\left(\epsilon \right)$ (10)

If the G-variation $\delta R\left(e^0;\delta e\right)$ is linear in $\delta e$, then the function $\delta R\left(e^0;\delta e\right)$ is called the G-differential of $R\left(e\right)$ and is usually denoted as $DR\left(e^0;\delta e\right)$. Furthermore, the result of the differentiations indicated on the right-side of the definition provided in Equation (8) can be written as follows:

$DR\left(e^0;\delta e\right)={\left\{DR\left(e^0;\delta \alpha \right)\right\}}^{direct}+{\left\{DR\left(e^0;\delta u,\delta v\right)\right\}}^{indirect},$ (11)

where the so-called “direct-effect” term is defined as follows:

${\left\{DR\left(e^0;\delta \alpha \right)\right\}}^{direct}\triangleq {\left\{\frac{\partial R}{\partial \alpha }\right\}}_{\left(e^0\right)}\delta \alpha ,$ (12)

while the so-called “indirect-effect” term is defined as follows:

${\left\{DR\left(e^0;\delta u,\delta v\right)\right\}}^{indirect}\triangleq {\left\{\frac{\partial R}{\partial u}\right\}}_{\left(e^0\right)}\delta u\left(x\right)+{\left\{\frac{\partial R}{\partial v}\right\}}_{\left(e^0\right)}\delta v\left(y\right).$ (13)

In Equations (12) and (13), the vectors $\partial R/\partial u$, $\partial R/\partial v$ and $\partial R/\partial \alpha$ comprise, as components, the first-order partial G-derivatives computed at the phase-space point $e^0$. The G-differential $DR\left(e^0;\delta e\right)$ is an operator defined on the same domain as $R\left(e\right)$ and has the same range as $R\left(e\right)$. The G-differential $DR\left(e^0;\delta e\right)$ satisfies the relation $R\left(e^0+\epsilon \delta e\right)-R\left(e^0\right)=DR\left(e^0;\delta e\right)+\Delta \left(\delta e\right)$, with $\underset{\epsilon \to 0}{\lim }\left[\Delta \left(\epsilon \delta e\right)\right]/\epsilon =0$.

The “direct effect” term ${\left\{DR\left(e^0;\delta \alpha \right)\right\}}^{direct}$ depends only on the parameter variations $\delta \alpha$ and can therefore be computed immediately, since it does not depend on the variations $\delta u$ and $\delta v$. On the other hand, the “indirect effect” term ${\left\{DR\left(e^0;\delta u,\delta v\right)\right\}}^{indirect}$ depends indirectly on the parameter variations $\delta \alpha$ through the yet unknown variations $\delta v$ and $\delta v$ in the state functions, which are the solutions of the system of equations obtained by applying the definition of the G-differential to Equations (1)-(3), to obtain the following relations:

$\begin{array}{l}{\left\{\frac{\text{d}}{\text{d}\epsilon }A\left[u^0\left(x\right)+\epsilon \delta u\left(x\right);{\alpha }^0+\epsilon \delta \alpha \right]\right\}}_{\epsilon =0}\\ ={\left\{\frac{\text{d}}{\text{d}\epsilon }{Q}^{\left(A\right)}\left({\alpha }^0+\epsilon \delta \alpha ;x\right)\right\}}_{\epsilon =0},x\in \delta {\Omega }_x,\end{array}$ (14)

$\begin{array}{l}{\left\{\frac{\text{d}}{\text{d}\epsilon }B\left[v^0\left(y\right)+\epsilon \delta v\left(y\right);{\alpha }^0+\epsilon \delta \alpha \right]\right\}}_{\epsilon =0}\\ ={\left\{\frac{\text{d}}{\text{d}\epsilon }{Q}^{\left(B\right)}\left({\alpha }^0+\epsilon \delta \alpha ;y\right)\right\}}_{\epsilon =0},y\in \delta {\Omega }_y,\end{array}$ (15)

$\begin{array}{l}{\left\{\frac{\text{d}}{\text{d}\epsilon }C\left[u^0\left(x\right)+\epsilon \delta u\left(x\right),v^0\left(y\right)+\epsilon \delta v\left(y\right);{\alpha }^0+\epsilon \delta \alpha ;x,y\right]\right\}}_{\epsilon =0}\\ =0,x\in \delta {\Omega }_x,y\in \delta {\Omega }_y.\end{array}$ (16)

Performing in Equations (14)-(16) the differentiations with respect to $\epsilon$ and setting $\epsilon =0$ in the resulting expressions yields the following system of equations:

${\left\{\frac{\partial A\left(u;\alpha \right)}{\partial u}\right\}}_{\left(e^0\right)}\delta u\left(x\right)={\left\{Q_1^{\left(1\right)}\left(u;\alpha ;\delta \alpha \right)\right\}}_{\left(e^0\right)},$ (17)

${\left\{\frac{\partial B\left(v;\alpha \right)}{\partial v}\right\}}_{\left(e^0\right)}\delta v\left(y\right)={\left\{Q_2^{\left(1\right)}\left(v;\alpha ;\delta \alpha \right)\right\}}_{\left(e^0\right)},$ (18)

$\begin{array}{l}{\left\{\frac{\partial C\left[u\left(x\right),v\left(y\right);\alpha ;x,y\right]}{\partial u}\right\}}_{\left(e^0\right)}\delta u\left(x\right)+{\left\{\frac{\partial C\left[u\left(x\right),v\left(y\right);\alpha ;x,y\right]}{\partial v}\right\}}_{\left(e^0\right)}\delta v\left(y\right)\\ +{\left\{\frac{\partial C\left[u\left(x\right),v\left(y\right);\alpha ;x,y\right]}{\partial \alpha }\right\}}_{\left(e^0\right)}\delta \alpha =0,x\in \delta {\Omega }_x,y\in \delta {\Omega }_y.\end{array}$ (19)

where

${\left\{Q_1^{\left(1\right)}\left(u,\alpha ;\delta \alpha \right)\right\}}_{\left(e^0\right)}\triangleq {\left\{\frac{\partial \left[Q^{\left(A\right)}\left(\alpha ;x\right)-A\left[u\left(x\right);\alpha \right]\right]}{\partial \alpha }\right\}}_{\left(e^0\right)}\delta \alpha ,$ (20)

${\left\{Q_2^{\left(1\right)}\left(v,\alpha ;\delta \alpha \right)\right\}}_{\left(e^0\right)}\triangleq {\left\{\frac{\partial \left[Q^{\left(B\right)}\left(\alpha ;y\right)-B\left[v\left(y\right);\alpha \right]\right]}{\partial \alpha }\right\}}_{\left(e^0\right)}\delta \alpha ,$ (21)

The system of equations comprising Equations (17)-(19) is called the “First-Level Forward Sensitivity System” (1^st-LFSS) and could be solved to obtain the variations $\delta v$ and $\delta v$ in the state functions in terms of the parameter variations $\delta \alpha$ which appear as sources in the 1^st-LFSS equations. Subsequently, the variations $\delta v$ and $\delta v$ thus obtained could be used to compute the indirect-effect term defined in Equation (13).

However, since there are at least $Z_{\alpha }$ variations to consider, it becomes prohibitively expensive computationally to solve in practice the 1^st-LFSS, which may comprise differential and or integral operators, for all possible parameter variations $\delta {\alpha }_i,i=1,\cdots ,Z_{\alpha }$. The need for solving repeatedly the 1^st-LFSS for every possible parameter variation $\delta {\alpha }_i,i=1,\cdots ,Z_{\alpha }$ can be circumvented by applying the concepts first outlined by Cacuci [2] [3] to construct a “First-Level Adjoint Sensitivity System” (1^st-LASS), the solution of which will be used to eliminate the appearance of the variations $\delta v$ and $\delta v$ in the expression of the indirect-effect term defined in Equation (13). The 1^st-LASS is constructed by implementing the following sequence of steps:

1) Introduce a Hilbert space pertaining to the domain ${\Omega }_x$, denoted as $H_u$, comprising square-integrable vector-valued elements of the same form as the vectors $u\left(x\right)$ and $\delta u\left(x\right)$. The inner product underlying $H_u$, between two

elements $g^{\left(\alpha \right)}\left(x\right)\triangleq {\left[g_1^{\left(\alpha \right)}\left(x\right),\cdots ,g_{Z_u}^{\left(\alpha \right)}\left(x\right)\right]}^{†}\in H_u$ and $g^{\left(\beta \right)}\left(x\right)\triangleq {\left[g_1^{\left(\beta \right)}\left(x\right),\cdots ,g_{Z_u}^{\left(\beta \right)}\left(x\right)\right]}^{†}\in H_u$ is denoted as ${\langle g^{\left(\alpha \right)}\left(x\right),g^{\left(\beta \right)}\left(x\right)\rangle }_u$ and defined as follows:

${\langle g^{\left(\alpha \right)}\left(x\right),g^{\left(\beta \right)}\left(x\right)\rangle }_u\triangleq {\displaystyle \underset{n=1}{\overset{Z_u}{\sum }}{\displaystyle \underset{a_1\left(\alpha \right)}{\overset{b_1\left(\alpha \right)}{\int }}\text{d}{x}_1\cdots {\displaystyle \underset{a_j\left(\alpha \right)}{\overset{b_j\left(\alpha \right)}{\int }}\text{d}{x}_j\cdots {\displaystyle \underset{a_{Z_x}\left(\alpha \right)}{\overset{b_{Z_x}\left(\alpha \right)}{\int }}{g}_n^{\left(\alpha \right)}\left(x\right)g_n^{\left(\beta \right)}\left(x\right)\text{d}{x}_{Z_x}}}}}$ (22)

2) Introduce a Hilbert space pertaining to the domain ${\Omega }_y$, denoted as $H_v$, comprising square-integrable vector-valued elements of the same form as the vectors $v\left(y\right)$ and $\delta v\left(y\right)$, i.e., $h^{\left(\alpha \right)}\left(y\right)\triangleq {\left[h_1^{\left(\alpha \right)}\left(y\right),\cdots ,h_{Z_v}^{\left(\alpha \right)}\left(y\right)\right]}^{†}\in H_v$ and $h^{\left(\beta \right)}\left(y\right)\triangleq {\left[h_1^{\left(\beta \right)}\left(y\right),\cdots ,h_{Z_v}^{\left(\beta \right)}\left(y\right)\right]}^{†}\in H_v$ The Hilbert space $H_v$ is endowed with an inner product denoted as ${\langle h^{\left(\alpha \right)}\left(y\right),h^{\left(\beta \right)}\left(y\right)\rangle }_v$, which is defined as follows:

${\langle h^{\left(\alpha \right)}\left(y\right),h^{\left(\beta \right)}\left(y\right)\rangle }_v\triangleq {\displaystyle \underset{n=1}{\overset{Z_y}{\sum }}{\displaystyle \underset{c_1\left(\alpha \right)}{\overset{d_1\left(\alpha \right)}{\int }}\text{d}{y}_1\cdots {\displaystyle \underset{c_j\left(\alpha \right)}{\overset{d_j\left(\alpha \right)}{\int }}\text{d}{y}_j\cdots {\displaystyle \underset{c_{Z_y}\left(\alpha \right)}{\overset{d_{Z_y}\left(\alpha \right)}{\int }}{h}_n^{\left(\alpha \right)}\left(y\right)h_n^{\left(\beta \right)}\left(y\right)\text{d}{y}_{Z_y}}}}}$ (23)

3) In the Hilbert $H_u$, form the inner product of Equation (17) with a yet undefined vector-valued function ${\psi }_1^{\left(1\right)}\left(x\right)\triangleq {\left[{\psi }_{11}^{\left(1\right)}\left(x\right),\cdots ,{\psi }_{1Z_u}^{\left(1\right)}\left(x\right)\right]}^{†}\in H_u$ to obtain the following relation:

${\langle {\psi }_1^{\left(1\right)}\left(x\right),{\left\{\frac{\partial A\left(u;\alpha \right)}{\partial u}\right\}}_{\left(e^0\right)}\delta u\left(x\right)\rangle }_u={\langle {\psi }_1^{\left(1\right)}\left(x\right),{\left\{Q_1^{\left(1\right)}\left(u;\alpha ;\delta \alpha \right)\right\}}_{\left(e^0\right)}\rangle }_u.$ (24)

4) Using the definition of the adjoint operator in the Hilbert space $H_u$, recast the left-side of Equation (24) as follows:

$\begin{array}{l}{\langle {\psi }_1^{\left(1\right)}\left(x\right),{\left\{\frac{\partial A\left(u;\alpha \right)}{\partial u}\right\}}_{\left(e^0\right)}\delta u\left(x\right)\rangle }_u\\ ={\langle \delta u\left(x\right),{\left\{A^{*}\left(u;\alpha \right)\right\}}_{\left(e^0\right)}{\psi }_1^{\left(1\right)}\left(x\right)\rangle }_u+{\left\{P_A^{\left(1\right)}{\left[\delta u\left(x\right);{\psi }_1^{\left(1\right)}\left(x\right);\alpha \right]}_{\delta {\Omega }_x}\right\}}_{\left(e^0\right)},\end{array}$ (25)

where ${\left\{P_A^{\left(1\right)}{\left[\delta u\left(x\right);{\psi }_1^{\left(1\right)}\left(x\right);\alpha \right]}_{\delta {\Omega }_x}\right\}}_{\left(e^0\right)}$ denotes the bilinear concomitant evaluated on the boundary $\delta {\Omega }_x$. In Equation (25), the operator $A^{*}\left(u;\alpha \right)$ is the formal adjoint of $\partial A\left(u;\alpha \right)/\partial u$.

5) Replace the left-side of Equation (24) by the right-side of Equation (25) to obtain the following relation:

$\begin{array}{l}{\langle \delta u\left(x\right),{\left\{A^{*}\left(u;\alpha \right)\right\}}_{\left(e^0\right)}{\psi }_1^{\left(1\right)}\left(x\right)\rangle }_u\\ ={\langle {\psi }_1^{\left(1\right)}\left(x\right),{\left\{Q_1^{\left(1\right)}\left(u;\alpha ;\delta \alpha \right)\right\}}_{\left(e^0\right)}\rangle }_u-{\left\{P_A^{\left(1\right)}{\left[\delta u\left(x\right);{\psi }_1^{\left(1\right)}\left(x\right);\alpha \right]}_{\delta {\Omega }_x}\right\}}_{\left(e^0\right)}.\end{array}$ (26)

6) In the Hilbert $H_v$, form the inner product of Equation (18) with a yet undefined vector-valued function ${\psi }_2^{\left(1\right)}\left(y\right)\triangleq {\left[{\psi }_{21}^{\left(1\right)}\left(y\right),\cdots ,{\psi }_{2Z_v}^{\left(1\right)}\left(y\right)\right]}^{†}\in H_v$ to obtain the following relation:

${\langle {\psi }_2^{\left(1\right)}\left(y\right),{\left\{\frac{\partial B\left(v;\alpha \right)}{\partial v}\right\}}_{\left(e^0\right)}\delta v\left(y\right)\rangle }_v={\langle {\psi }_2^{\left(1\right)}\left(y\right),{\left\{Q_2^{\left(1\right)}\left(v;\alpha ;\delta \alpha \right)\right\}}_{\left(e^0\right)}\rangle }_v$. (27)

7) Using the definition of the adjoint operator in the Hilbert space $H_v$, recast the left-side of Equation (24) as follows:

$\begin{array}{l}{\langle {\psi }_2^{\left(1\right)}\left(y\right),{\left\{\frac{\partial B\left(v;\alpha \right)}{\partial v}\right\}}_{\left(e^0\right)}\delta v\left(y\right)\rangle }_v\\ ={\langle \delta v\left(y\right),{\left\{B^{*}\left(v;\alpha \right)\right\}}_{\left(e^0\right)}{\psi }_2^{\left(1\right)}\left(y\right)\rangle }_v+{\left\{P_B^{\left(1\right)}{\left[\delta v\left(y\right);{\psi }_2^{\left(1\right)}\left(y\right);\alpha \right]}_{\delta {\Omega }_y}\right\}}_{\left(e^0\right)},\end{array}$ (28)

where ${\left\{P_B^{\left(1\right)}{\left[\delta v\left(y\right);{\psi }_2^{\left(1\right)}\left(y\right);\alpha \right]}_{\delta {\Omega }_y}\right\}}_{\left(e^0\right)}$ denotes the bilinear concomitant evaluated on the boundary $\delta {\Omega }_y$. In Equation (28), the operator $B^{*}\left(v;\alpha \right)$ is the formal adjoint of $\partial B\left(v;\alpha \right)/\partial v$.

8) Replace the left-side of Equation (27) by the right-side of Equation (28) to obtain the following relation:

$\begin{array}{l}{\langle \delta v\left(y\right),{\left\{B^{*}\left(v;\alpha \right)\right\}}_{\left(e^0\right)}{\psi }_2^{\left(1\right)}\left(y\right)\rangle }_v\\ ={\langle {\psi }_2^{\left(1\right)}\left(y\right),{\left\{Q_2^{\left(1\right)}\left(v;\alpha ;\delta \alpha \right)\right\}}_{\left(e^0\right)}\rangle }_v-{\left\{P_B^{\left(1\right)}{\left[\delta v\left(y\right);{\psi }_2^{\left(1\right)}\left(y\right);\alpha \right]}_{\delta {\Omega }_y}\right\}}_{\left(e^0\right)}.\end{array}$ (29)

9) Add Equations (29) and (26) to obtain:

$\begin{array}{l}{\langle \delta u\left(x\right),{\left\{A^{*}\left(u;\alpha \right)\right\}}_{\left(e^0\right)}{\psi }_1^{\left(1\right)}\left(x\right)\rangle }_u+{\langle \delta v\left(y\right),{\left\{B^{*}\left(v;\alpha \right)\right\}}_{\left(e^0\right)}{\psi }_2^{\left(1\right)}\left(y\right)\rangle }_v\\ ={\langle {\psi }_1^{\left(1\right)}\left(x\right),{\left\{Q_1^{\left(1\right)}\left(u;\alpha ;\delta \alpha \right)\right\}}_{\left(e^0\right)}\rangle }_u+{\langle {\psi }_2^{\left(1\right)}\left(y\right),{\left\{Q_2^{\left(1\right)}\left(v;\alpha ;\delta \alpha \right)\right\}}_{\left(e^0\right)}\rangle }_v\\ -{\left\{P_A^{\left(1\right)}{\left[\delta u\left(x\right);{\psi }_1^{\left(1\right)}\left(x\right);\alpha \right]}_{\delta {\Omega }_x}\right\}}_{\left(e^0\right)}-{\left\{P_B^{\left(1\right)}{\left[\delta v\left(y\right);{\psi }_2^{\left(1\right)}\left(y\right);\alpha \right]}_{\delta {\Omega }_y}\right\}}_{\left(e^0\right)}.\end{array}$ (30)

10) The next step is to relate the right-side of Equation (30) with the indirect-effect term ${\left\{DR\left(e^0;\delta u,\delta v\right)\right\}}^{indirect}$ defined in Equation (13). Since the response considered is a functional of $u$ and $v$, the G-differential $R\left(e\right)$ is also a functional of $\delta u\left(x\right)$ and $\delta v\left(y\right)$. Consequently, the well-known Riesz representation theorem (which states that every functional can be expressed uniquely in terms of the inner product pertaining to the respective Hilbert space) ensures that the indirect-effect term ${\left\{DR\left(e^0;\delta u,\delta v\right)\right\}}^{indirect}$ can be expressed uniquely as follows:

${\left\{DR\left(e^0;\delta u,\delta v\right)\right\}}^{indirect}\triangleq {\langle \delta u\left(x\right),{\left\{{\left(\frac{\partial R}{\partial u}\right)}^{†}\right\}}_{\left(e^0\right)}\rangle }_u+{\langle \delta v\left(y\right),{\left\{{\left(\frac{\partial R}{\partial v}\right)}^{†}\right\}}_{\left(e^0\right)}\rangle }_v.$ (31)

11) Identifying the right-side of Equation (31) with the left-side of Equation (30) indicates that the indirect-effect term ${\left\{DR\left(e^0;\delta u,\delta v\right)\right\}}^{indirect}$ would be equal to the right side of Equation (30) provided that the following relations are satisfied by the yet undetermined functions ${\psi }_1^{\left(1\right)}\left(x\right)$ and ${\psi }_2^{\left(1\right)}\left(y\right)$ :

${\left\{A^{*}\left(u;\alpha \right)\right\}}_{\left(e^0\right)}{\psi }_1^{\left(1\right)}\left(x\right)={\left\{{\left(\frac{\partial R}{\partial u}\right)}^{†}\right\}}_{\left(e^0\right)}$ (32)

${\left\{B^{*}\left(v;\alpha \right)\right\}}_{\left(e^0\right)}{\psi }_2^{\left(1\right)}\left(y\right)={\left\{{\left(\frac{\partial R}{\partial v}\right)}^{†}\right\}}_{\left(e^0\right)}.$ (33)

12) Using Equations (31)-(33) in Equation (30) transforms the latter into the following form:

$\begin{array}{l}{\left\{DR\left(e^0;\delta u,\delta v\right)\right\}}^{indirect}\\ ={\langle {\psi }_1^{\left(1\right)}\left(x\right),{\left\{Q_1^{\left(1\right)}\left(u;\alpha ;\delta \alpha \right)\right\}}_{\left(e^0\right)}\rangle }_u+{\langle {\psi }_2^{\left(1\right)}\left(y\right),{\left\{Q_2^{\left(1\right)}\left(v;\alpha ;\delta \alpha \right)\right\}}_{\left(e^0\right)}\rangle }_v\\ -{\left\{P_A^{\left(1\right)}{\left[\delta u\left(x\right);{\psi }_1^{\left(1\right)}\left(x\right);\alpha \right]}_{\delta {\Omega }_x}\right\}}_{\left(e^0\right)}-{\left\{P_B^{\left(1\right)}{\left[\delta v\left(y\right);{\psi }_2^{\left(1\right)}\left(y\right);\alpha \right]}_{\delta {\Omega }_y}\right\}}_{\left(e^0\right)}.\end{array}$ (34)

13) The boundary, interface and initial/final conditions for the functions ${\psi }_1^{\left(1\right)}\left(x\right)$ and ${\psi }_2^{\left(1\right)}\left(y\right)$ are now determined by imposing the following requirements:

a) Implement the boundary, interface and initial/final conditions given in Equation (19) into the bilinear concomitants in Equation (34).

b) Eliminate the remaining unknown boundary, interface and initial/final conditions involving the functions $\delta u\left(x\right)$ and $\delta v\left(y\right)$ from the expression of the bilinear concomitants in Equation (34) by selecting boundary, interface and initial/final conditions for the functions ${\psi }_1^{\left(1\right)}\left(x\right)$ and ${\psi }_2^{\left(1\right)}\left(y\right)$ such that the selected conditions for ${\psi }_1^{\left(1\right)}\left(x\right)$ and ${\psi }_2^{\left(1\right)}\left(y\right)$ must be independent of unknown values of $\delta u\left(x\right)$, $\delta v\left(y\right)$ and $\delta \alpha$ while ensuring that Equations (32) and (33) are well posed. The boundary conditions thus chosen for the adjoint functions ${\psi }_1^{\left(1\right)}\left(x\right)$ and ${\psi }_2^{\left(1\right)}\left(y\right)$ can be represented in operator form as follows:

${\left\{C_A^{\left(1\right)}\left[u\left(x\right);v\left(y\right);{\psi }_1^{\left(1\right)}\left(x\right),{\psi }_2^{\left(1\right)}\left(y\right);\alpha ;x,y\right]\right\}}_{\left(e^0\right)}=0,x\in \partial {\Omega }_x,y\in \delta {\Omega }_y$ (35)

where the subscript “A” indicates “adjoint”.

14) The selection of the boundary conditions for the adjoint functions ${\psi }_1^{\left(1\right)}\left(x\right)$ and ${\psi }_2^{\left(1\right)}\left(x\right)$ represented by Equation (35) eliminates the appearance of any unknown values of the variations $\delta u\left(x\right)$ and $\delta v\left(y\right)$ in the bilinear concomitants in Equation (34) and reduces these concomitants to a residual quantity that contains boundary terms involving only known values of $\delta \alpha$, $u\left(x\right)$, $v\left(y\right)$, ${\psi }_1^{\left(1\right)}\left(x\right)$, ${\psi }_2^{\left(1\right)}\left(y\right)$, $\alpha$. This residual quantity will be denoted as ${\left\{{\hat{P}}^{\left(1\right)}\left[u\left(x\right);v\left(y\right);{\psi }_1^{\left(1\right)}\left(x\right),{\psi }_2^{\left(1\right)}\left(y\right);\alpha ;x,y;\delta \alpha \right]\right\}}_{\left(e^0\right)}$. In general, this residual quantity does not automatically vanish, although it may do so in particular instances. In principle, ${\left\{{\hat{P}}^{\left(1\right)}\left[u\left(x\right);v\left(y\right);{\psi }_1^{\left(1\right)}\left(x\right),{\psi }_2^{\left(1\right)}\left(y\right);\alpha ;x,y;\delta \alpha \right]\right\}}_{\left(e^0\right)}$ could be forced to vanish, if necessary, by considering extensions, in the operator sense, of the linear operators $A^{*}\left(u;\alpha \right)$ and/or $B^{*}\left(v;\alpha \right)$, but such extensions seldom need to be used in practice.

15) Using the conditions represented by Equations (19) and (35) in Equation (34) yields the following (final) expression for the indirect-effect term ${\left\{DR\left(e^0;\delta u,\delta v\right)\right\}}_{indirect}$ :

$\begin{array}{l}{\left\{DR\left(e^0;\delta u,\delta v\right)\right\}}^{indirect}\\ =-{\left\{{\hat{P}}^{\left(1\right)}\left[u\left(x\right);v\left(y\right);{\psi }_1^{\left(1\right)}\left(x\right),{\psi }_2^{\left(1\right)}\left(y\right);\alpha ;x,y;\delta \alpha \right]\right\}}_{\left(e^0\right)}\\ +{\langle {\psi }_1^{\left(1\right)}\left(x\right),{\left\{Q_1^{\left(1\right)}\left(u;\alpha ;\delta \alpha \right)\right\}}_{\left(e^0\right)}\rangle }_u+{\langle {\psi }_2^{\left(1\right)}\left(y\right),{\left\{Q_2^{\left(1\right)}\left(v;\alpha ;\delta \alpha \right)\right\}}_{\left(e^0\right)}\rangle }_v\\ \equiv {\left\{DR\left(e^0;{\psi }_1^{\left(1\right)},{\psi }_2^{\left(1\right)}\right)\right\}}^{indirect}.\end{array}$ (36)

As the expression in Equation (36) indicates, the desired elimination of the unknown variations $\delta u$ and $\delta v$ from the original expression of ${\left\{DR\left(e^0;\delta u,\delta v\right)\right\}}^{indirect}$ given in Equation

(13) has been accomplished by having replaced them by expressions involving the functions ${\psi }_1^{\left(1\right)}\left(x\right)$ and ${\psi }_2^{\left(1\right)}\left(y\right)$, which do not depend on any parameter variations, a fact that has been underscored by having explicitly indicated that the indirect-effect term can now be written in the form

${\left\{DR\left(e^0;{\psi }_1^{\left(1\right)},{\psi }_2^{\left(1\right)}\right)\right\}}^{indirect}$.

The system of equations represented by Equations (32), (33), and (35) is called the First-Level Adjoint Sensitivity System (1^st-LASS) and the functions ${\psi }_1^{\left(1\right)}\left(x\right)$ and ${\psi }_2^{\left(1\right)}\left(y\right)$ are called the “first-level adjoint sensitivity functions.” The essential feature of the 1^st-LASS is that it is independent of parameter variations (in contradistinction to the 1^st-LFSS), so it needs to be solved only once per response to obtain the first-level adjoint sensitivity functions ${\psi }_1^{\left(1\right)}\left(x\right)$ and ${\psi }_2^{\left(1\right)}\left(y\right)$. Once the adjoint functions ${\psi }_1^{\left(1\right)}\left(x\right)$ and ${\psi }_2^{\left(1\right)}\left(y\right)$ are available, they can be used

in Equation (36) to compute the indirect-effect term ${\left\{DR\left(e^0;\delta u,\delta v\right)\right\}}_{indirect}$

exactly and efficiently, using quadrature formulas, which are many orders of magnitude faster to compute then solving the operator (differential, integral) equations that underlie the 1^st-LFSS. As is well known [2] [3] [4] [5], it is this property that makes the adjoint sensitivity analysis method “unbeatable” when needing to compute the sensitivities of functional-valued responses to many imprecisely known parameters.

4. Concluding Remarks

This work has presented the first-order comprehensive adjoint sensitivity analysis methodology (1^st-CASAM) for computing efficiently, exhaustively and exactly, the first-order response sensitivities for coupled nonlinear physical systems characterized by imprecisely known parameters characterizing the systems, the interfaces between systems and the systems’ domain boundaries. The 1^st-CASAM fundamentally generalizes and extends all previously published theoretical works on this topic, also enabling the quantification of the effects of manufacturing tolerances on the responses of physical and engineering systems. The 1^st-CASAM highlights the conclusion that response sensitivities to the imprecisely known domain boundaries and interfaces can arise both from the definition of the system’s response as well as from the equations, interfaces and boundary conditions defining the model and its imprecisely known domain. Ongoing research will generalize the methodology presented in this work, aiming at computing exactly and efficiently higher-order response sensitivities for coupled systems involving imprecisely known interfaces, parameters, and boundaries. The sequel [20] to this work illustrates the application of the 1^st-CASAM to a benchmark problem [21] [22] [23] that models heat conduction and convection in a physical system comprising an electrically heated rod surrounded by a coolant which simulates the geometry of an advanced (“Generation-IV”) nuclear reactor [24].

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Wigner, E.P. (1945) Effect of Small Perturbations on Pile Period. Chicago Report CP-G-3048.
[2]	Cacuci, D.G. (1981) Sensitivity Theory for Nonlinear Systems: I. Nonlinear Functional Analysis Approach. Journal of Mathematical Physics, 22, 2794-2802. https://doi.org/10.1063/1.525186
[3]	Cacuci, D.G. (1981) Sensitivity Theory for Nonlinear Systems: II. Extensions to Additional Classes of Responses. Journal of Mathematical Physics, 22, 2803-2812. https://doi.org/10.1063/1.524870
[4]	Cacuci, D.G. (2003) Sensitivity and Uncertainty Analysis: Theory, Volume 1. Chapman & Hall/CRC, Boca Raton. https://doi.org/10.1201/9780203498798
[5]	Cacuci, D.G. (2018) The Second-Order Adjoint Sensitivity Analysis Methodology. CRC Press, Taylor & Francis Group, Boca Raton. https://doi.org/10.1201/9781315120270
[6]	Valentine, T.E. (2006) Polyethylene-Reflected Plutonium Metal Sphere Subcritical Noise Measurements, SUB-PU-METMIXED-001; International Handbook of Evaluated Criticality Safety Benchmark Experiments; NEA/NSC/DOC(95)03/I-IX; Organization for Economic Co-Operation and Development (OECD). Nuclear Energy Agency (NEA), Paris.
[7]	Cacuci, D.G., Fang, R. and Favorite, J.A. (2019) Comprehensive Second-Order Adjoint Sensitivity Analysis Methodology (2nd-ASAM) Applied to a Subcritical Experimental Reactor Physics Benchmark: I. Effects of Imprecisely Known Microscopic Total and Capture Cross Sections. Energies, 12, 4219. https://doi.org/10.3390/en12214219
[8]	Fang, R. and Cacuci, D.G. (2019) Comprehensive Second-Order Adjoint Sensitivity Analysis Methodology (2nd-ASAM) Applied to a Subcritical Experimental Reactor Physics Benchmark: II. Effects of Imprecisely Known Microscopic Scattering Cross Sections. Energies, 12, 4114. https://doi.org/10.3390/en12214114
[9]	Cacuci, D.G., Fang, R., Favorite, J.A., Badea, M.C. and Di Rocco, F. (2019) Comprehensive Second-Order Adjoint Sensitivity Analysis Methodology (2nd-ASAM) Applied to a Subcritical Experimental Reactor Physics Benchmark: III. Effects of Imprecisely Known Microscopic Fission Cross Sections and Average Number of Neutrons per Fission. Energies, 12, 4100. https://doi.org/10.3390/en12214100
[10]	Fang, R. and Cacuci, D.G. (2020) Comprehensive Second-Order Adjoint Sensitivity Analysis Methodology (2nd-ASAM) Applied to a Subcritical Experimental Reactor Physics Benchmark. IV: Effects of Imprecisely Known Source Parameters. Energies, 13, 1431. https://doi.org/10.3390/en13061431
[11]	Fang, R. and Cacuci, D.G. (2020) Comprehensive Second-Order Adjoint Sensitivity Analysis Methodology (2nd-ASAM) Applied to a Subcritical Experimental Reactor Physics Benchmark: V. Computation of 2nd-Order Sensitivities Involving Isotopic Number Densities. Energies, 13, 2580. https://doi.org/10.3390/en13102580
[12]	Cacuci, D.G., Fang, R. and Favorite, J.A. (2020) Comprehensive Second-Order Adjoint Sensitivity Analysis Methodology (2nd-ASAM) Applied to a Subcritical Experimental Reactor Physics Benchmark: VI. Overall Impact of 1st- and 2nd-Order Sensitivities. Energies, 13, 1674. https://doi.org/10.3390/en13071674
[13]	Cacuci, D.G. (2019) Towards Overcoming the Curse of Dimensionality: The Third-Order Adjoint Method for Sensitivity Analysis of Response-Coupled Linear Forward/Adjoint Systems, with Applications to Uncertainty Quantification and Predictive Modeling, Energies, 12, 4216. https://doi.org/10.3390/en12214216
[14]	Komata, M. (1977) A Generalized Perturbation Theory Applicable to Reactor Boundary Changes. Nuclear Science and Engineering, 64, 811-822. https://doi.org/10.13182/NSE77-A14496
[15]	Larsen, E.W. and Pomraning, G.C. (1981) Boundary Perturbation Theory. Nuclear Science and Engineering, 77, 415-425. https://doi.org/10.13182/NSE81-A18954
[16]	Rahnema, F. and Pomraning, G.C. (1983) Boundary Perturbation Theory for Inhomogeneous Transport Equations. Nuclear Science and Engineering, 84, 313-319. https://doi.org/10.13182/NSE83-A15451
[17]	McKinley, M.S. and Rahnema, F. (2002) High-Order Boundary Condition Perturbation Theory for the Neutron Transport Equation. Nuclear Science and Engineering, 140, 285-294. https://doi.org/10.13182/NSE02-A2261
[18]	Favorite, J.A. and Gonzalez, E. (2017) Revisiting Boundary Perturbation Theory for Inhomogeneous Transport Problems. Nuclear Science and Engineering, 185, 445-459. https://doi.org/10.1080/00295639.2016.1277108
[19]	Cacuci, D.G. (2018) BERRU Predictive Modeling: Best Estimate Results with Reduced Uncertainties. Springer, Heidelberg/New York. https://doi.org/10.1007/978-3-662-58395-1
[20]	Cacuci, D.G. (2020) The First-Order Comprehensive Sensitivity Analysis Methodology (1st-CASAM) for Scalar-Valued Responses: II. Illustrative Application. American Journal of Computational Mathematics;, accepted for publication.
[21]	Cacuci, D.G., Fang, R., Ilic, M. and Badea, M.C. (2015) A Heat Conduction and Convection Analytical Benchmark for Adjoint Solution Verification of CFD Codes Used in Reactor Design. Nuclear Science and Engineering, 182, 452-480. https://doi.org/10.13182/NSE15-69
[22]	Cacuci, D.G. (2016) Second-Order Adjoint Sensitivity and Uncertainty Analysis of a Heat Transport Benchmark Problem—I: Analytical Results. Nuclear Science and Engineering, 183, 1-21. https://doi.org/10.13182/NSE15-81
[23]	Cacuci, D.G., Ilic, M., Badea, M.C. and Fang, R. (2016) Second-Order Adjoint Sensitivity and Uncertainty Analysis of a Heat Transport Benchmark Problem—II: Computational Results Using G4M Reactor Thermal-Hydraulic Parameters. Nuclear Science and Engineering, 183, 22-38. https://doi.org/10.13182/NSE15-80
[24]	GEN4 ENERGY, INC. (2012) Reactor Core Design. GEN4 ENERGY, INC., Denver.
[25]	(2015) ANSYS® Academic Research, Release 16.0, FLUENT Adjoint Solver, ANSYS, Inc.