An adjoint‐assisted multilevel multifidelity method for uncertainty quantification and its application to turbomachinery manufacturing variability

In this work we propose, analyze, and demonstrate an adjoint‐based multilevel multifidelity Monte Carlo (MLMF) framework called FastUQ. The framework is based on the MLMF of Geraci et al. and uses the Inexpensive Monte Carlo (IMC) method of Ghate as low‐fidelity surrogate. The setup cost of IMC‐1 surrogate in FastUQ requires just the adjoint solution at the input mean whose computational cost is independent of the number of input uncertainties making it suitable for solving problems with a large number of uncertain parameters. We demonstrate the robustness of the method to quantify uncertainties in aerodynamic parameters due to surface variations caused by the manufacturing processes for a highly loaded turbine cascade. A stochastic model for surface variations on the cascade is proposed and optimal dimensionality reduction of model parameters is realised using goal‐based principal component analysis using adjoint sensitivities of multiple quantities of interest. The proposed method achieves a 70% reduction in computational cost in predicting the mean quantities such as total‐pressure loss and mass flow rate compared to the state‐of‐art MLMC method.


INTRODUCTION
Determining performance deterioration of turbomachines due to manufacturing variations is important to engine manufacturers when determining trade-offs between manufacturing cost and performance. Many studies have been conducted in the literature with focus on change in performance due to manufacturing variations on blade shapes. [3][4][5][6][7][8][9] The main challenge of modeling the effect of geometric variations on the blade performance is the large number of uncertain parameters that have to be considered. Lange et al. 8 considered blade measurements from three-dimensional (3D) surface scan data of a high pressure compressor stage. The correlated coordinates from the 3D scan were translated into mode shapes and uncorrelated modal amplitudes using PCA. Lange et al. included all 135 PCA modes in their perturbation model within the MC simulations and compared it with MC results considering 20 and 60 PCA modes. Lange 8 concludes that 60 PCA modes (half of the total number of PCA modes) were necessary to correctly model the manufacturing variability on the blades. Even with dimensionality reduction the resulting parameter space is too high 8 to achieve the specified error tolerance. Therefore uncertainty quantification (UQ) methods that suffer from the curse of dimensionality (CoD) are impractical to solve this problem. 10 UQ using the Monte Carlo method (MC) has been widely adopted 4,5,8,9 due to the following key reasons; (i) convergence of statistical moments is independent of the number of input uncertainties, 11,12 (ii) computational cost versus statistical errors can be established clearly hence the statistical estimates can be improved on the fly with additional samples, 13 and (iii) the method can be applied to highly nonlinear and nonsmooth systems. 14 The computational cost of MC becomes prohibitively expensive for problems where the computational cost per sample is high due to the slow convergence with the number of samples. This has motivated many improvements to the method, specifically the control variate acceleration approach. 1 The control variate naturally lends itself to multifidelity and multilevel formulations. In the context of robust aerospace design, Ng and Wilcox 15 proposed a practical multifidelity control variate Monte Carlo (MFMC) implementation and its extension to an information reuse estimator 16 that incorporates information from previous design steps. The adaptive MFMC (AMFMC) of Peherstorfer et al. 17 splits the computational budget between adapting the low-fidelity models to improve their approximation quality and sampling the low-and high-fidelity models to reduce the mean square error of the estimator. Giles 18 proposed the Multilevel Monte Carlo (MLMC) method, which uses a sequence of coarsened meshes to represent the low fidelity model of the control variate. Mishra 19 et al. showed that MLMC has the same asymptotic complexity (up to a logarithmic term) as a single deterministic solution considering a hierarchy of nested meshes. Geraci et al. 1 proposed a multilevel multifidelity framework (MLMF) to achieve further reduction in computational cost by combining the multifidelity control variate in MLMC. Gorodetsky 20 gives a comprehensive overview of MLMF models and unify many existing sampling methods used for MF UQ. They propose nonrecursive/nested strategies for sampling schemes and achieve significant computational gains.
The MFMC, AMFMC, or MLMF frameworks require a surrogate model for the low-fidelity sample evaluation. Surrogate models do incur setup cost through initial high-fidelity sampling thus suffering from CoD. For example in the Kriging meta-model requires design of experiments (DOE) to sample the response at multiple points in the parameter space to create the Kriging model. Mostly the Latin hypercube sampling (LHS) is used for sampling the response for DOE. 14, 21 Nobile 22 used a stochastic collocation to construct the low-fidelity surrogate model at each level in the MLMC for elliptic type PDE. Reynolds averaged Navier-Stokes (RANS) equations considered in this work is a highly nonlinear mixed PDE that admits various types of flow discontinuities. In addition, to construct the surrogate model one requires samples at the sparse grid collocation points. Hence all the aforementioned methods are not suitable for problems with high input dimension.
Fairbanks 23 introduced a low-rank approximation in MLMC for higher dimensional problems called the multilevel control variates (MLCV). The low-rank approximation is obtained using a small set of selected samples of the fine grid solution identified from the realisations of the coarse grid solution. The method outperforms MLMC when the problem exhibits small solution ranks on all the levels. In addition, they show that MLCV is beneficial over MLMC only for low mean-squared error tolerance. Blonigan 24 used a Reduced Order Model (ROM) in a MF modeling approach based on the MF model of Gorodetsky 20 to solve the Kuramoto-Sivashinsky equation. Ghate 2 proposed an adjoint-based surrogate model called the incomplete Monte Carlo (IMC) as a replacement for high fidelity flow simulations to study the performance of turbomachines and airfoils subject to random manufacturing variations. A similar model was used by Rumpfkeil 25,26 as an approximate surrogate model for lift and drag cost functions of an airfoil in aerodynamic shape design. Engels-Putzka et al. 27 used the IMC model to estimate the mass flow rate, total-pressure ratio and total-temperature ratio in the two-stage Darmstadt Transonic Compressor using the TRACE flow and adjoint solver developed at DLR for use in robust design. Luo et al. 28 adopted a similar approach of using the adjoint and higher order sensitivities to approximately evaluate the aerodynamic performance of a turbine subject to manufacturing variability using Taylor expansion. They used a continuous adjoint approach based on inviscid Euler equations to obtain the adjoint sensitivity and Hessian solutions. The adjoint gradient used in the IMC model of Ghate has the property that its computational cost is proportional to the number of output variables and independent of input dimension. Therefore the adjoint surrogate combined with MFMC or MFML forms an attractive combination because the convergence of the resulting method is still independent of the input dimension.
We find that in all the previous works based on the adjoint-based surrogates (similar to Ghate's IMC), the model was used as a replacement for the high-fidelity evaluation. We demonstrate in this work the weakness of using the IMC surrogate model for highly nonlinear problems using the uncertain viscous Burgers' equation and aerodynamic performance variations due to manufacturing uncertainties on the surface of a turbine blade. Errors as large as 35 − 225% in SD and 2 − 11% in mean values were obtained for the viscous Burgers' problem and 6% variation in mean and 10% variation in Sd was observed for the total-pressure loss for the turbine problem using IMC. In fact the higher-order adjoint corrections have larger errors compared to the lower order ones. Interestingly higher and lower order IMC models correlate well with the exact solution. Therefore we propose to use the adjoint-surrogate as a low-fidelity model in the proposed multifidelity multilevel UQ framework. To the authors' knowledge this is the first reported work on the use of the adjoint-based IMC surrogate as a low-fidelity model within a MLMF framework.
The MLMF framework proposed in this work is based on the work of Geraci et al. 1 but we use the multifidelity cost model similar to the one used in Ng. 15 This clearly brings out the definition of the multifidelity parameter in the multilevel multifidelity framework. We introduce a Gaussian process model for surface variations due to manufacturing process of turbine blades. Unlike previous work in the literature 4,6,27,29 we consider multiple cost functions in a goal-based principal component analysis (PCA) of the surface variation model. We show that including multiple cost functions increases the the number of PCA modes required to achieve the specified error tolerance in the cost functions. This corroborates with the high dimensional nature of the problem of uncertainty in aerodynamic performance due to manufacturing variations of blades. Lastly to the authors knowledge this one of the first work to use an MLMF framework for the problem of aerodynamic UQ due manufacturing variations.

Multifidelity control variate acceleration
Let J be a quantity of interest (QoI), which we try to estimate using P high-fidelity (HF) samples to yield . ⏟⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏟ Bias error (1) The term J MC p is the MC estimate of the mean obtained using the P samples, that is, J MC p = P −1 ∑ P p=1 J p . For an unbiased estimator the MSE in Equation (1) can be reduced either by, (1) increasing the number of samples P or (2) reducing the variance of J p . The control variate takes the latter approach and reduces J p by using M samples from a low-fidelity (LF) model G (G m : {m = 1, … , M}). Typically G is much cheaper to evaluate than J although less accurate. Variance reduction is achieved by introducing the new estimator involving J, G and a free parameter (called the control variate parameter) shown in Equation (2).
In Equation (3), w is the runtime ratio between the average computational cost per HF and LF model evaluation and r = M P > 1 is the ratio between the number of HF and LF simulations. c is the computational budget (an equivalent HF sample count). For a fixed budget c we are free to choose the parameter r to vary the resource allocation between HF and LF evaluations. The variance of the multifidelity estimator in terms of c, r, and w is shown in Equation (4). Pasupathy et al. 31 16 and Geraci 1 ). This assumption is justified when the number of LF samples M used to estimate the mean E[G] is much greater than the number of high fidelity samples P, that is, (M ≫ P).
For optimal resource allocation, one minimizes the variance in Equation (4) for both and r to yield the optimality condition 15 in Equation (5).
where, r opt and opt are, When the cost of computing the LF samples is almost free, that is, w → ∞, then r opt → ∞ and one recovers the classical control variate estimator, that is, c . In the case of a perfect LF model, that is, → 1 then the error scales similar to the statistical error of a classical MC considering only the LF model. Therefore, the LF model should not only correlate well but also should be cheaper to evaluate.

Convergence comparison of single-and multifidelity MC
A similar comparison of the MF estimator can be done with the MSE of a classical MC using P samples from Equation (6) (neglecting bias error). The parameter Λ l 1,20 is a penalising factor with respect to an optimal control variate, which is the control variate obtained by using a low-fidelity model for which the statistics are known and therefore the variance reduction is (1 − 2 ). The contour plots of and Λ are plotted in Figure 1. The grey regions denote the regions where the control variate is ineffective and a classical MC performs better.
Note that there are no terms involving P in Equation (7) so the multifidelity MC has the same  convergence but the entire convergence curve is shifted by Λ. The shift is only a function of the correlation and runtime ratio w.

Multilevel multifidelity control variate
Evaluation of the computational model on a coarser mesh reduces computational cost and accuracy. Using sample evaluations on a hierarchy of coarser meshes as LF model in the control variate method produces the multilevel MC (MLMC). 18 Let J N be the QoI obtained from solving a PDE on a mesh with N degrees-of-freedom. Therefore as N → ∞ then J N → J and . Using the linearity property of expectation, the expectation at the finest level L (with N L degrees of freedom) can be written using the telescopic sum of the sequence of levels l = 0, … L as shown in Equation (9). The difference function Y l defined in Equation (10) is introduced to simplify the analysis.
For P l samples distributed per level we get the MLMC estimator in Equation (11) of Giles 18 who also derived the optimal samples per level P opt l by optimising the total computational cost c = ∑ L l=0 P l C l for a specified MSE tolerance on the estimator in Equation (11) to yield Equation (12). C l is the cost of a single sample evaluation at level l. Note that Giles bounds the bias error in estimating Y l to be less than the prescribed tolerance divided by the square root of two ( 1 . Mishra 19 et al. showed that MLMC has the same asymptotic complexity (up to a logarithmic term) as a single deterministic solution considering a hierarchy of nested meshes. Therefore MLMC can converge QoI statistics rapidly; analogous to the multigrid convergence for PDE solution. Further variance reduction can be achieved by combining the multifidelity control variate into a multilevel one. Geraci et al. 1 combined both MFMC and MLMC and proposed a multilevel mutlifidelity MC (MLMF) framework. The number of LF samples used per level in the work of Geraci was fixed as an increment Δ of the number of HF samples M = P + Δ, where Δ = P is some fraction of the HF samples. The multifidelity computational budget per level of Geraci using our notation is, Comparing the multifidelity cost model in Equation (13) with the one in Equation (3) we obtain the equation r = 1 + , which relates the parameters used in both models. Replacing the cost model gives us the benefit of using the multifidelity estimator and clearly show the effect of introducing the multilevel into the multifidelity control variate, which is shown in detail in Section 2.4. To simplify the analysis the same number of grid levels are used for models J and G. The estimators in Equations (2) and (11) can be combined to obtain the Y l correlation MLMF estimator, 13 J MLMF N, p of QoI J shown in Equation (14).
where, Y l is the control variate at level l shown in Equation (15) and Z l is the difference function 1 of the LF model G shown in Equation (16).

Optimal sample allocation between levels and models
Let P l and M l denote the number of sample evaluations of the HF and LF models at level l. Let r l denote the fraction of HF evaluation to the LF evaluation r l = M l /P l . The cost of computation on a level is expressed in terms of the equivalent number of HF sample evaluations as shown in Equation (17), where C l is a measure of the computation cost of HF sample on level l. We normalise C l by dividing it by the computational cost of evaluating a HF sample at the finest level N (note that C N = 1 and C l ≤ 1). Therefore the budget c l now contains both the relative computational cost between models and across levels, The total computational budget c considering all levels l = 0, … L is then The variance of this estimator can be optimised for l and r l using a Lagrangian multiplier approach. 1 The final optimal variance and the optimal number of HF samples P opt l per level is shown in Equation (19)- (20).
The values of r opt l and Λ l are shown in Equation (21)- (22) respectively. The parameter Λ l 1,20 is a penalising factor with respect to an optimal control variate, which is the control variate obtained by using a low-fidelity model for which the statistics are known and therefore the variance reduction is (1 − 2 ). l is the Pearson correlation between the HF and LF model.
To simplify the analysis, consider the correlation between the HF and LF model across levels to be the same ( l = ). Then the optimal number of samples in Equation (20) reduces to the simple expression in Equation (23).
The optimal number of samples per level P opt l in Equation (23) is identical to the one obtained by Giles (shown in Equation (12)) except for the additional Λ term in the summation. In fact, for the special case of Λ = 1 we recover Equation (12). The contour plot comparing (r opt ) and Λ(r opt ) for a range of values of correlation and runtime ratio w is shown in Figure 1(A),(B). We find that the runtime ratio w should be ≥10 and correlation should be ≥0.6 for the MLMF to be effective. P opt l is directly proportional to Λ l (r opt l ), which is a strong function of the correlation. For large runtime ratios w, Λ(r opt ) ∝ and is almost independent of runtime ratio. Therefore the success of the MLMF is strongly dependent on how well the LF model correlates with the HF model. In addition P opt l ∝ −2 so lowering the tolerance by half quadruples the number of samples per level similar to the convergence of MC. The value of Λ(r opt ) is mostly ≤1 indicating a reduction in the number of samples compared to the standard MLMC to achieve the same MSE.

Note on pilot samples
Practically, in the multifidelity control variate it is not possible to compute E[J p ] CV and Var[J p ] CV opt using Equations (4) and (5), because E[J], Var[J], Var[G], , and are unknown quantities and the usual practice is to replace them with sample estimates, which are called pilot samples. The pilot sampling is also necessary for evaluating the variance estimate at a given level Var(Y k ). For a budget c, one has N samples of J p , which can be used to estimate the necessary statistics. Replacing the exact values with approximate ones introduces errors. Therefore, it is necessary to check if the estimator still yields a reduction in variance. This is checked by plotting the ratio shown in Equation (24) (from Reference 15). In In this work a LHS was employed to select the initial pilot sample population. To simplify the implementation of the multilevel multifidelity variance estimate we use the same number of pilot samples to estimate both multilevel and multifidelity parameters. 1 Multiple initial sample sizes are chosen and we measure the difference in the mean and SD for some fixed MSE tolerance . The threshold sample size above which the difference in mean and SD are within the error tolerance ( ) is chosen and used in subsequent calculations. Therefore the sampling is an iterative process, where we estimate the optimal number of samples at each level based on cost and observed variance. At each iteration we evaluate additional samples, which results in a more accurate variance estimation across level and model fidelities. This improved estimate is used to update the sample allocation for the next iteration. The iteration is terminated when either zero additional samples are allocated or upon reaching a prescribed maximum iteration. Note that a similar iterative sampling process is adopted for the SMLMC results presented in this work.

ADJOINT SURROGATE MODEL
Motivated by the adjoint error correction of Giles and Pierce, 33 Ghate 2 introduced an adjoint-based approximate output evaluation called the Inexpensive Monte Carlo method (IMC). The adjoint error correction is used to improve an approximation for the QoI J(u * , ), where, u * is an approximation for the state u due to perturbations in the design parameter as shown in Equation (25). The linearised state equation shown in Equation (26) is substituted in Equation (25) to simplify the adjoint correction as shown in Equation (27).
The adjoint v T in Equation (27) is evaluated at the perturbed state u . But one can replace the adjoint with an approximate one v * T and the leading error term includes the error in the adjoint solution approximation as shown below: Although the IMC formulation in Equation (28) is equivalent to the moment method 2 a few major differences exist. Using IMC one can estimate the complete pdf of the QoI J but moment methods only provide an estimate of the statistical moments like mean, variance, etc. IMC requires the evaluation of the nonlinear residual R at the perturbed state, which captures nonlinearities in the QoI more effectively than the moment methods. Overall in IMC there are three possible combinations to the approximation of u * and v * T , which leads to the IMC 1-3 formulations. The computational cost of IMC can be broken down into two parts, (i) the initial setup cost (ISC) and (ii) the evaluation cost (EC). ISC is incurred only to setup the model and is a one time cost but EC is the cost associated with every sample evaluation using IMC. In all IMC variants the EC cost is the same, which involves the cost of calculating J(u * , ) and R(u * , ). Usually the cost of evaluating the objective function J is quite inexpensive compared to the evaluation of the nonlinear residual R. A typical nonlinear RANS solution for one sample point requires (10 3 − 10 4 ) residual evaluations. Therefore the cost of one residual evaluation is negligible in comparison. IMC 1: Approximates the flow and adjoint values using the unperturbed baseline solution. This method has an overall leading error of second order and it is equivalent to the first order moment method. The ISC of this method is a single adjoint solution at the mean state. Note that cost of the linear adjoint solution is independent of the number of input uncertainties and is cheaper to evaluate than the computation of the solution to the nonlinear state equation, R(u) = 0.
IMC 2: Approximates the flow values from the baseline tangent-linear solution and extrapolation to the perturbed state. The adjoint is approximated using the baseline solution similar to IMC 1. This method has an overall leading error of third order and it is equivalent to the second order moment method. For ∈ R M IMC 2 requires M tangent-linear solutions in addition to the baseline nonlinear and adjoint solution.
IMC 3: Approximates the flow and adjoint values using gradient extrapolation of the tangent-linear flow and adjoint solution. This is the most expensive, yet, most accurate method. The method has an overall leading error of order 4. For ∈ R M , IMC 3 requires M tangent-linear flow and adjoint (Hessian) solutions. The tangent-linear adjoint solution is calculated using either the Tangent-on-Tangent or Tangent-on-Reverse strategy. 34

Burgers' equation with uncertain boundary condition
The viscous Burgers' equation can be used as a simplified model for the displacement of shock waves in a compressible nonlinear RANS model under uncertain perturbation. As a case study, consider the viscous Burgers' equation defined in the interval x ∈ [−1, 1], with Dirichlet boundary conditions specified at the two end points x = −1 and x = 1 as shown in Equations (32)- (34). An uncertain parameter is specified at the boundary x = −1, which controls the perturbation from the unit value.
Xiu 35 gives the detailed analytic solution procedure for this problem. The steady-state solution is shown in Equation (35), where x 0 is the location of the transition point u(x 0 ) = 0, and A is the slope − u∕ x| x=x 0 at the transition point. This slope is obtained as a solution to Equation (36), which is solved iteratively using a root-finding algorithm. Once A is determined, x 0 is obtained by substituting A into any one of the Equations (33) and (34). Note that the viscosity determines the steepness of the shock at x = x 0 and a fixed value of viscosity = 1 20 was used in all simulations.
The steady-state solution of the viscous Burgers' equation is plotted in Figure 3 for various values of . We consider the location of the transition point x 0 at steady-state (denoted by the dark circle in Figure 3) as the QoI in this work. The value of x 0 is readily obtained by substituting the value of A in Equation (37).
For small values of the following asymptotic solution for x 0 at steady-state is shown to be valid. 35 Note that the QoI x 0 is hypersensitive to boundary perturbation because even small changes in bring about huge variations in x 0 . We assume the input uncertainty to occur equally likely in the given interval [a, b]. Therefore the occurrence of is modeled as a random variable with a uniform probability distribution p( ), with mean = 2 −1 (b − a) and variance Var[ ] = 12 −1 (b − a) 2 , shown in Equation (39).
For the UQ study in this work the values of [a, b] ≡ [0, 0.1] is chosen. This model problem is simple, but sufficiently complex to assess the accuracy and reliability of the surrogate models.

Approximating model response using adjoint correction
Consider the QoI J to be the transition point location x 0 and one can form the non-linear state equation R(A, ) from Equation (36) as, For the IMC estimator, one first linearises about the mean of the input disturbance ∈ [0, 0.1], which is = 0.05. The perturbation about this mean can be evaluated from the adjoint and tangent-linear solution using any of the IMC schemes described in Equations (29)-(31). In Figure 4 the QoI obtained from the various IMC estimators are plotted for a range of perturbations along with the exact value. The mean and SD for the exact and IMC 1 − 3 models and the Pearson correlation between the exact and IMC 1 − 3 are shown in Table 1. The higher-order adjoint surrogates IMC 2 and 3 suffer from large errors in the mean and variance compared to lower order model IMC 1 in stark contrast to what one would expect. But they predict the right trends in the solution which is reflected in the better correlation.
The results of the Burgers' model equation provide two important insights. For highly sensitive nonlinear problems adjoint surrogates do not provide reliable estimates of the statistical moments (especially the variance). In fact, higher order corrections are actually detrimental. However, they do capture the mean and general trends in the function variations (better correlation). An attractive feature of the lower-order adjoint surrogate IMC 1 is that the computational cost is independent of the number of input variables and requires one adjoint solution for each QoI. Higher-order IMC 2/3 incur additional computational cost, which scales linearly with the number of input uncertainties for IMC 2 and quadratically for IMC 3. Therefore by using higher-order IMC 2/3 we loose the attractive convergence property of classical MC; independence of convergence to the number of input uncertainties. Good correlation and low computational cost make IMC 1 an ideal candidate for a low-fidelity model in a MFMC or MFML framework which we show in the next section.

STOCHASTIC MODEL FOR MANUFACTURING UNCERTAINTIES
In this work a synthetic surface perturbation model based on GP similar to the approach presented in References 3,4,36,37 has been adopted due to a lack of access to empirical blade measurements. In addition, parametrization of the perturbations is avoided since this would lead to a reduction in the space of possible perturbed geometries. 38 Therefore, a free-node parametrization 38 is used where every surface node is a random parameter defining the blade shape. The infinite dimensional probability space is approximated using a finite number of random variables, where every surface mesh node is a parameter. 37 Let x = {x 1 , x 2 , … , x N }, where each x i ∈ R 3 , be the set of coordinates defining the nominal surface of the blade and letn denote the normal vector of this nominal surface. A zero-mean Gaussian process (x) is imposed along the normal direction to this nominal surface to generate the perturbed surface x as, In addition, the random variables i = (x i ) and j = (x j ) at any two arbitrary points i and j on the surface (see Figure 5) are assumed to have a squared exponential spatial covariance as shown in Equation (42). The squared exponential covariance has been extensively used in robust optimisation studies involving surface perturbations (see References 6,29,39). This covariance function is infinitely differentiable, which means that the GP with this covariance function has derivatives of all orders (in the mean squared sense), and is thus very smooth. 40 In Equation (42), the parameter b controls the height of the perturbations and l is the characteristic correlation length of the perturbation. A large value of l results in a larger spread of the disturbance around a given surface node. In this work, the perturbation height b and the correlation length l are assumed be a fraction of some characteristic dimension. The parameters b and l are illustrated for the surface perturbation over the LS89 cascade surface in Figure 5. The model is used to represent the surface variations over the LS89 turbine cascade, which is presented in the results section.
F I G U R E 5 Surface node distribution on LS89 cascade surface along with the two arbitrary surface nodes i and j used to illustrate the spatial covariance function in Equation (42)

Dimensionality reduction and goal-based truncation
The convergence properties of MC are valid under the assumption of independent identically distributed (i.i.d) input parameter space. 30 Therefore, statistical independence is an important property to have for the input variations. For a joint Gaussian distribution to be independent its covariance function must be uncorrelated, that is, the covariance matrix C ij in Equation (42) should be transformed to a diagonal matrix to make the input uncertainties independent.
The expansion using the eigenvalues and vectors of the discrete zero-mean surface perturbation is shown below where, Y i ( )'s are uncorrelated Gaussian random variables with zero mean and unit variance. 41 The expansion is usually truncated to the first m dominant eigenmodes for reducing the input dimension. Using Mercer's theorem one can show that the truncated expansion is a suitable approximation, if the eigenvalues decay sufficiently fast and m is sufficiently large. 38 For a Gaussian covariance function the eigenvalues will exponentially decay toward zero. 42 Therefore, with a fraction of the total number of modes the variations can be approximated reliably.
The modal fraction i and the partial modal fraction Λ i shown in Equation (44) are useful quantities to aid the truncation. Partial modal fraction can be used to estimate the value of m given a user-specified threshold error. For example, if one truncates after the ith PCA mode and Λ i = 0.9 then 90% of the total input variations (spectral content) are captured using the truncated set. A value of i = 0.1 denotes that the ith PCA mode accounts for 10% of the total input variations (spectral content).
The impact of an input variation on the QoI is encoded in the adjoint sensitivity. Therefore the dot product of PCA modes and the sensitivity gives the influence of the mode on the QoI. We call this scalar the mode effectiveness ( PCA ).
Since the adjoint sensitivity is independent of the number of input parameters, it is evaluated once and used for all the PCA modes. The modes are sorted based on the PCA value and truncated using a new modal fraction defined in Equation (46). Note that Equation (46) is obtained by replacing with the PCA in Equation (44). The projected PCA modes are called the goal-based PCA modes or simply G-PCA.
Most references on G-PCA applied to surface perturbation models consider only a single QoI or a single flow condition for the truncation, for example, Reference 36. In this work, multiple QoIs at multiple flow conditions are considered for the G-PCA truncation (see Section 5.3).

FastUQ numerical implementation
To the authors knowledge this is the first work to use MLMF method for the UQ of surface variations. Therefore, a practical implementation of this approach is provided for reference. The method is shown schematically in Figure 6. To begin with the G-PCA modes and amplitudes are extracted from the GP model of the surface variations. The modes are typically obtained for the finest mesh (in the sequence of meshes) to maximize the space of possible perturbed geometries. The selected modes are interpolated from the surface mesh of the finest to all coarse levels. A piecewise-linear interpolation along the cascade surface (see Figure 7) is employed in this work to interpolate the PCA modes. The local parametric coordinate t along the cascade surface used for the piecewise-linear interpolation is shown in Figure 7. The interpolation ensures that the samples across levels belong to the same random path. 18 Independent Gaussian random samples with zero-mean and unit variance are generated for the selected G-PCA modes. One can reconstruct the surface disturbance sample using the weighted summation in Equation (43). A parallel distance-weighted interpolation (IDW) mesh smoothing is used to propagate the surface perturbations into the volume.
The perturbed mesh is then used to run either the HF or LF model to obtain the output QoI. Therefore, the total computational cost of a sample evaluation is a sum of the computational cost of mesh smoothing and QoI evaluation (using HF or LF model).
The FastUQ sampler runs the initial sampling for the LF and HF model at each level to construct the model correlations and multilevel parameters. Then based on the MSE tolerance the number of additional LF and HF samples are estimated at each level using Equation (20). Iteratively samples are added based on the MSE estimates derived in Section 2.4.

F I G U R E 7 Interpolation of Principle Component
Analysis mode #4 from the fine to the medium and coarse meshes 5 NUMERICAL RESULTS

Test case description
The LS89 turbine cascade was originally designed and optimized at the Von Karman Institute for Fluid Dynamics (VKI) for a subsonic isentropic exit Mach number of 0.9 using an iterative inverse method that modifies a given blade profile such that it fits an input surface velocity distribution. 43 Experimental measurements consisting of surface pressure and wall heat transfer have been conducted for a range of inlet total-pressure and exit back pressure ratios. 44,45 Two test conditions, namely the MUR47 (transonic flow) and MUR43 (subsonic flow) are considered in this work. Each condition corresponds to a specific back pressure at the outlet and total temperature at the inlet. The details of the test conditions are tabulated in Tables 2 and 3. The surface isentropic Mach number M isen is obtained using the relation shown below where P 01 is the total-pressure at the inlet and P s is the static pressure at a measurement station. A constant total temperature of 420 K was maintained at the inlet for all test conditions. The aerodynamic UQ due to surface variation using FastUQ is demonstrated using the LS89 turbine cascade. The aim here is to estimate the statistics (mean and variance) of the QoIs mass flow rate J mass and total pressure loss J ploss defined below subject to the surface perturbations on the cascade surface using the GP model described in the previous section. Note that the quantities with( ) are mass averaged values, p, rho, u are the fluid pressure, density, and velocity, respectively, and the superscript in/out denotes the values at the inlet/outlet of the domain. The surface disturbance model has a height b approximately 10% of the trailing edge radius and a length l approximately equal to the leading edge diameter of the cascade. The various geometric parameters of the cascade are shown in Table 2. Two levels of coarsening of the original mesh were performed maintaining a near wall spacing y + ≈ 1 − 2. The fine, medium and coarse meshes have approximately 90k, 60k, and 40k nodes, respectively.

Flow and sensitivity validation
The mesh and the computational domain used for the validation is shown in Figure 8. The surface isentropic Mach number (M isen ) obtained using CFD is plotted against the experimental results of Arts 44 for the MUR43 to MUR48 test conditions test conditions in Figure 9. Good overall agreement is obtained between the experimental results and the RANS simulation.
The surface sensitivity verification is conducted using the LS89 turbine cascade test case. The MUR47 test condition is considered and the gradients of the QoI, namely exit mass flow rate and total-pressure loss are verified in this study. The gradient results obtained using adjoint, tangent-linear and finite-difference are compared using the directional derivative (Gâteaux) (Figure 12), which needs the definition of a direction vector h. In this verification exercise, the surface modes obtained using the PCA shown in Section 5.3 are considered as the surface direction vector h ≡ z j (see Figure 10). Then using Taylor expansion one obtains the finite-difference form (central difference) as shown below,

Mode 2
Step size Mode 4 Step size

Mode 8
Mass flow rate

Total-pressure loss
Step size The term J(x ± 2 z j ) is obtained from the flow solution at the perturbed mesh due to the surface deformation x ± 2 z j . Choosing a small value for yields large round-off errors and using bigger values of results in larger truncation error. Thus the optimal value of is a trade-off between the round-off and truncation error. This optimal value is obtained by plotting the error between the actual gradient and the one obtained from finite difference for various step sizes. The relative error in the directional gradient (between adjoint and finite-difference) is shown for PCA modes 2, 4, and 8, respectively, for the total-pressure loss and exit mass flow rate cost function in Figure 11. Comparison of the directional derivatives obtained from the optimal finite-difference step size, tangent-linear and adjoint method are plotted in Figure 12. Good agreement is found between all three methods for the directional gradient values indicating a correct implementation.

Goal-based truncation results
The eigenvalues of the first ten eigenmodes of the PCA and corresponding G-PCA eigenvalue of the surface perturbation over the LS89 cascade are plotted in Figure 13 for both test conditions MUR43 and MUR47. The ranking of the first 10 dominant mode numbers are tabulated in Table 4  Combining the (most dominant) five G-PCA modes in each QoI for the MUR43 test case, that is, {1, 4, 5, 6, 8} one can capture 50% of the variations (in the G-PCA sense). Considering the seven modes {1, 3, 4, 5, 6, 7, 8} one can capture 50% of the spectral content of the input variations (G-PCA) considering both test conditions and QoIs. Clearly modes that are effective for a particular QoI and test condition need not be effective for another. Therefore, when more QoIs and test conditions are considered the number of G-PCA modes necessary to capture the spectral content of the input variations increases.

MUR43 and MUR47 flow Condition
The MUR43 and MUR47 flow conditions are chosen for quantifying the aerodynamic uncertainty due to surface variations. The Mach number contours for the two flow conditions are shown in Figure 15. MUR47 has a supersonic exit isentropic Mach number of 1.02 and MUR43 has a subsonic exit isentropic Mach number of 0.84. One can notice the presence of a shock wave in MUR47 near the trailing edge of the cascade, which extends all the way into the exit. The total-pressure loss QoI is especially sensitive due to the presence of this shock. The mean and SD (obtained from 20  random pilot samples) in the surface isentropic Mach number is plotted in Figure 16. This limited sample estimate reveals the high sensitivity of the surface isentropic Mach number to the surface perturbation G-PCA modes after mid-chord of the cascade geometry on the suction side. The variations are quite high for MUR47 due to the presence of the shock aft of the cascade.

Cost model parameters
The SMLMC and FastUQ adjoint-assisted MLMF require a cost model for the multilevel nonlinear flow solution. The computational cost for the baseline (unperturbed) mesh with a residual convergence tolerance of 10 −11 is tabulated in Table 5. The runtime ratio w l for each level is also shown in Table 5. Note that the runtime ratio for both MUR43 and MUR47 are quite similar (see Table 5). Hence a fixed runtime ratio of {0.2 : 0.35 : 1.0} was used throughout. The computational cost for the baseline geometry (nominal or unperturbed cascade surface) shown in Table 5 is obtained using a uniform flow initialization. Note that the flow field for every HF sample evaluation in the UQ analysis is initialized using this nominal solution to reduce the number of iterations required for convergence. In addition, the number of iterations to convergence varies depending on the size of the perturbation. If the perturbations are large then more iterations are necessary and vice versa. To simplify the analysis, the runtime ratio w l is assumed fixed based on the nominal cost and variations in computational cost due to difference in the perturbation size are ignored. The cost of the adjoint solution and the relative cost ratio with respect to the HF solution at each level is shown in Table 6. Note that the adjoint solution is significantly cheaper than the nonlinear solution because the Jacobian matrix    Table 7. Using the cost of the adjoint (see Table 6) and tangent solution (see Table 7) one can estimate the setup cost of the IMC 1/2 model using the Equation (50) shown below.
Note that the cost model shown in Equation 50 is for a single flow condition whose cost of adjoint and tangent linear solution is AdjointCost and TangentCost. #QoI and #Modes in Equation 50 denote the number of QoIs and G-PCA modes included in the analysis.

Effect of pilot sampling
The accuracy of the variance estimate at each level depends on the size of the initial samples set considered, which are also called the pilot samples (see Section 2.5). One needs to ensure that the errors introduced by the pilot sampling is minimal because it can deteriorate the computational gain (see Equation (24)). To test the error in the pilot sampling, the number of pilot samples was doubled from 10 to 20 and the change in QoI for an MSE tolerance ( ) of 0.1 was recorded for the SMLMC and is shown in Table 8. LHS was employed to select the initial pilot samples. The results show that the change in mean and SD with the doubling of the pilot samples is negligible in comparison to the specified tolerance. Hence we used 10 pilot samples throughout this work and a maximum of 100 iterations was used for the MLMF iterative process described in Section 2.5. Both runs with 10 and 20 pilot samples converge within the prescribed tolerance and maximum number of iterations.

Effects of G-PCA truncation
To quantify the effect of G-PCA, truncation the UQ of the LS89 cascade for the MUR47 flow condition was carried out using the MLMC of Giles (see Section 2.3), which we refer to as the standard MLMC or SMLMC. Two sets of SMLMC simulations were run, (i) using 7 and (ii) 25 G-PCA modes (combining modes based on both QoIs and flow conditions). The seve G-PCA modes capture 50% and the 25 G-PCA modes capture 99% of the total variations in the input uncertainty. The cost model parameters from Section 5.5 are used in the current SMLMC simulation. For both sets 10 pilot samples were used per level to obtain the initial model correlations. An MSE tolerance of 0.1 was specified for the simulations. The mean and SD obtained from the simulation are tabulated in Table 9. Inclusion of more modes does not dramatically alter the mean value of the exit mass flow rate but the mean total-pressure loss is quite sensitive. Although the variation in the mean of the QoI (between 7 and 20 modes) is <2%, the SD has a high variation in the QoI (around 36%). Therefore, all simulations henceforth use 25 G-PCA modes since it captures 99% of G-PCA variation and the predicted variations with seven modes differed significantly from the 25 G-PCA modes. Inclusion of 25 G-PCA modes of the LS89 cascade increases the input uncertainty space quite dramatically. Note that the computational cost of IMC 1 surrogate model is independent of the number of PCA modes (or number of input uncertainties). So a large number of input uncertainties can be handled using FastUQ with IMC 1 LF model (see Section 3). In addition, the G-PCA can help truncate the statistically independent PCA modes that do not contribute to the change in QoI. Therefore is well suited to solve this problem.

Optimal resource allocation and MSE convergence
The SMLMC using 25 G-PCA modes was run for four MSE tolerances, 1.0, 0.1, 0.05, and 0.01, respectively, for the MUR47 flow condition. The resource allocation between the three mesh levels (coarse, medium, and fine) to achieve the aforesaid MSE tolerances is shown in Figure 17(A) (based on Equation (20) for optimal resource allocation). Results in Figure 17(A) show that all computations required sample evaluation only on the coarsest mesh up to an MSE tolerance of 0.05 (excluding the 10 pilot samples at each level). This is primarily due to well-correlated values of the QoI between levels, that is, the gross variations are captured by the coarse levels. This is an important requirement for multilevel methods. 46 The mean and SD values obtained for the various MSE error tolerance are tabulated in Table 10 where it can be observed that for an MSE tolerance of 0.01 more equivalent HF samples are necessary to capture the finer details of the variations. The equivalent HF samples are still cheaper to evaluate because most of the samples are computed from the coarse level solution (see Figure 17(A)). The convergence of the MSE with equivalent HF samples P is plotted in Figure 17(B), where the dotted trend-line indicates the (P −1 ) convergence. The initial convergence is sharp for MSE < 0.05 but for larger MSE > 0.05

Accuracy of adjoint correction
To test the accuracy of the IMC LF model, the aerodynamic UQ of LS89 cascade was carried out at the MUR47 condition for the two QoIs, namely, exit-mass flow and total-pressure loss. The IMC 1 LF model requires only a single adjoint solution per QoI and its computational cost is independent of the number of PCA modes (see Section 3). For each level the IMC 2 LF model requires as many tangent linear solutions as the number of G-PCA modes (in addition to the cost of the adjoint solution). Note that the above estimates are the number of solution evaluations to construct the LF model per level in the case of FastUQ and it is a one time cost associated with the LF model setup.
The results presented in Section 5.7 showed that 25 G-PCA modes were necessary to capture 99% of the input variations and using lower number of G-PCA modes gave large variations in the standard deviation of the QoI for the MUR47 test condition. The setup cost of the IMC 2 LF model for 25 modes is 17 equivalent HF samples (calculated based on Equation (50)). This setup cost is in fact greater than the computational cost in equivalent HF samples for the entire IMC 1 simulation for an MSE tolerance of 0.01 (see Table 11). Therefore, for a large number of input uncertainties the setup cost of IMC 2 becomes quite high. In addition, the MUR47 condition has a shock discontinuity in the solution. This reduces the regularity of the solution, which introduces large errors in the higher order adjoint corrections 2,47,48 (also see Section 3.2) Therefore, the analysis using IMC 2/3 models is postponed to a future study and only the IMC 1 LF model is considered in all simulations.
The results of the IMC 1 are compared against the SMLMC results in Table 11 for the same prescribed MSE tolerance of 0.01. Note that the sample allocation per level for SMLMC is shown in Figure 17(A) and the IMC 1 is run on the finest mesh level. The equivalent HF samples for the IMC 1 simulation is obtained by multiplying the total runtime of the IMC 1 simulation with the ratio of the SMLMC equivalent HF samples to the SMLMC runtime (for the prescribed MSE tolerance).
Two important observations can be made from the results shown in Table 11. Firstly, the IMC 1 takes only a fraction of the cost of the SMLMC (≈ 1 10 ) for the same MSE tolerance. The IMC took 1.7 h to converge to the prescribed tolerance compared to 27 h for the SMLMC. Secondly, the total-pressure loss (mean and variance) obtained using IMC 1 has large deviations from the SMLMC results but the deviation of the exit mass flow rate QoI is lower. This clearly shows the shortcoming of the IMC 1 model for sensitive and nonlinear QoIs, which can be remedied using the FastUQ method.

FastUQ results
The performance of the MLMF is a strong function of the correlation between the HF and IMC LF model (see Equation (20)). Therefore, the correlation at each level for the MUR43 and MUR47 test case is plotted in Figure 18 (obtained from the 10 pilot samples). One finds that the correlation is quite high for both MUR47 and MUR43 test conditions (average correlation for both cases is ≥0.95). But the coarse level in MUR43 has a lower correlation value of 0.76. Therefore, one can expect more HF sample evaluation at the coarse level for MUR43 to compensate for the lower correlation. This is clearly reflected in the resource allocation shown in Figure 22, where 151 additional HF samples are evaluated for the MUR43. The QoI statistics for both MUR43 and MUR47 test conditions were computed using all three UQ methods for an MSE tolerance of 0.01 and the results are summarised in Figures 19-21. The computational cost reduction in comparison to the runtime of the SMLMC achieved by FastUQ and IMC are shown in Figure 19. The percentage variation in the mean and SD from SMLMC results for FastUQ and IMC are shown in Figures 20 and 21. The computational cost comparison in Figure 19 shows that FastUQ gives ≈70% reduction in computational cost over SMLMC and the variations in mean and SD of the QoI are within 1% variations of the SMLMC results. IMC achieves an additional 25% reduction compared to FastUQ. Although the mean and SD for the exit mass flow rate QoI is captured within 1% tolerance, IMC has significantly large deviations for the total-pressure loss QoI. Variations as large as 10% can be observed in the SD and 6% in mean values.
Based on the resource allocation shown in Figure 22, FastUQ uses the LF model more in regions where correlations are better but resorts to running the HF model when the LF model has large deviations (lower correlation). Since the model correlations for the MUR47 test case are higher (in the range 0.98-0.99 as shown in Figure 18), more LF samples are evaluated using the IMC 1 model and the correlation is used to correct the QoI. On the other hand, the correlation for the coarse mesh is lower in MUR43 so additional HF samples are used to compensate this error.

SUMMARY
We present a framework for UQ based on multilevel multifidelity MC acceleration called FastUQ. Although the classical MC does not suffer from the CoD, multifidelity control variates can suffer from CoD due to the computational cost of constructing the low-fidelity surrogate model. In FastUQ we use a surrogate model based on adjoint sensitivity whose computational cost is independent of the number of input parameters, hence free from CoD. FastUQ is used to quantify uncertainty (mean and SD) in aerodynamic parameters like total-pressure loss and mass flow rate due to manufacturing variations on a highly loaded turbine cascade. The surface variations due to manufacturing were modelled as a gaussian process and goal-based PCA considering multiple adjoint sensitives of QoIs was used for dimensionality reduction. Even with dimensionality reduction, the resulting parameter space was quite high; 25 PCA modes for a 2D cascade. Results from Lange et al. 8 also corroborates with our findings. The mean and SD of QoIs considering the PCA modes obtained using FastUQ was compared to a simple multilevel method and running a classical MC using the IMC surrogate model. FastUQ provides estimates of mean and standard deviation with the same accuracy of MLMC and achieves 70% reduction in computational cost. For benign QoIs like mass flow rate IMC gives acceptable accuracy of mean and SD with much lower computational cost but suffers from large errors for highly sensitive QoIs namely the total-pressure loss especially for nonlinear flow conditions (flow with shocks). FastUQ can compensate for errors in the surrogate model by running high fidelity samples allocated optimally across multiple levels to achieve both accuracy and computational efficiency.

ACKNOWLEDGMENT
This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No. 642959.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request.