W N Edeling
email: w.n.edeling@tudelft.nl

R P Dwight
email: r.p.dwight@tudelft.nl

P Cinnella
email: p.cinnella@ensam.eu

Simplex-stochastic collocation method with improved scalability

Keywords:

The Simplex-Stochastic Collocation (SSC) method is a robust tool used to propagate uncertain input distributions through a computer code. However, it becomes prohibitively expensive for problems with dimensions higher than 5. The main purpose of this paper is to identify bottlenecks, and to improve upon this bad scalability. In order to do so, we propose an alternative interpolation stencil technique based upon the Set-Covering problem, and we integrate the SSC method in the High-Dimensional Model-Reduction framework. In addition, we address the issue of ill-conditioned sample matrices, and we present an analytical map to facilitate uniformly-distributed simplex sampling.

Introduction

In order to make reliable predictions of a physical system using a computer code it is necessary to understand what effect the uncertainties in the inputs have on the output Quantity of Interest (QoI). Attempts to do so while keeping the computational cost low can be found in [START_REF] Smolyak | Quadrature and interpolation formulas for tensor products of certain classes of functions[END_REF][START_REF] Nobile | A sparse grid stochastic collocation method for partial differential equations with random input data[END_REF][START_REF] Ganapathysubramanian | Sparse grid collocation schemes for stochastic natural convection problems[END_REF], which rely on Smolyak-type sparse-grid stochastic-collocation methods. Whereas traditional collocation-type methods [START_REF] Babuška | A stochastic collocation method for elliptic partial differential equations with random input data[END_REF][START_REF] Loeven | Probabilistic collocation: an efficient non-intrusive approach for arbitrarily distributed parametric uncertainties[END_REF] use full-tensor product formulas to extend a set of one-dimensional nodes to higher dimensions, sparse-grid methods build a sparse interpolant using a constrained linear combination of one-dimensional nodes. This can provide a grid with a potential reduction in support nodes of several orders of magnitude.

Although computationally more efficient than full-tensor approaches, sparse-grid methods add points equally in all dimensions, irrespective of whether the response surface is locally smooth or discontinuous. Therefore further gains can be achieved through adaptive stochastic-collocation schemes which have been developed in recent years. For instance Ma et al. [START_REF] Ma | An adaptive hierarchical sparse grid collocation algorithm for the solution of stochastic differential equations[END_REF] proposed an Adaptive Sparse-Grid (ASG) collocation method where the probabilistic space is spanned by linear finiteelement basis functions. During each iteration the probability space is refined locally through an error measure based upon the hierarchical surplus, defined as the difference between the interpolation of the previous level and the newly added code sample. Although the space is refined locally, unphysical oscillations can still occur due to the lack of sample points on some of the edges of the local support of the basis functions. The Simplex-Stochastic Collocation (SSC) method of Witteveen et al. [START_REF] Witteveen | Simplex stochastic collocation with random sampling and extrapolation for nonhypercube probability spaces[END_REF] circumvents this problem by discretizing the domain into simplices by means of a Delaunay triangulation, and enforcing the so-called Local-Extremum Conserving limiter to suppress unphysical oscillations. Furthermore, it computes Essentially Non-Oscillatory (ENO) stencils [START_REF] Witteveen | Simplex stochastic collocation with eno-type stencil selection for robust uncertainty quantification[END_REF] of the sample points which allows for high-order polynomial interpolation.

Further features include randomized sampling, the ability to deal with non-hypercube probability spaces and it can be extended to perform interpolation with sub-cell resolution [START_REF] Witteveen | Subcell resolution in simplex stochastic collocation for spatial discontinuities[END_REF].

Besides schemes which efficiently sample the probabilistic space, there are other means of dealing with the curse of dimensionality, i.e. the exponential increase in the amount of computational cost with increasing dimension. Physical systems often have a low effective dimension, meaning that only a few coefficients are influential, and only low-order correlations between these input coefficients have a significant impact on the output. To capitalize on this behavior, High-Dimensional Model-Reduction (HDMR) techniques can be applied [START_REF] Rabitz | Efficient input-output model representations[END_REF]. In HDMR a d-dimensional function is exactly represented as a hierarchical sum of 2 d component functions of increasing dimensionality. In the case of low effective dimension, the d-dimensional function can be approximated well by a truncated expansion. The main idea is to solve several low-dimensional sub-problems instead of one high-dimensional one, which greatly reduces the computational cost. A well-known member of this class of decompositions is the analysis of variance (ANOVA) decomposition. In [START_REF] Foo | Multi-element probabilistic collocation method in high dimensions[END_REF], Foo et al. successfully coupled their Multi-Element Probabilistic Collocation method [START_REF] Foo | The multi-element probabilistic collocation method (me-pcm): error analysis and applications[END_REF] with a HDMR-ANOVA decomposition to problems with up to 600 dimensions. Although they achieved a significant reduction in the computational cost compared to approaches with full-tensor products, the number of points required to sufficiently sample these extremely high-dimensional spaces is still of the order of millions or even more. Furthermore, in [START_REF] Ma | An adaptive high-dimensional stochastic model representation technique for the solution of stochastic partial differential equations[END_REF] Ma et al. coupled their previously mentioned ASG method with the so-called cut-HDMR technique of [START_REF] Rabitz | Efficient input-output model representations[END_REF]. This approach is computationally more efficient than ANOVA-HDMR, as it does not require the evaluation of multi-dimensional integrals. Besides truncating the cut-HDMR expansion at a certain order, the authors of [START_REF] Ma | An adaptive high-dimensional stochastic model representation technique for the solution of stochastic partial differential equations[END_REF] also made their approach dimension adaptive through weights which identify the dimensions that contribute the most to the mean of the stochastic problem. They applied this approach to several easilyevaluated test problems of very high dimensionality, i.e. up to a stochastic dimension of 500. Again, their results represent a significant reduction in the required number of code samples for a certain error level compared to full tensor grids, but in absolute terms the number of samples is still very high.

In our view, the interest of a surrogate modelling technique is to apply it to some complex simulation code which is too expensive to integrate by simple Monte Carlo techniques. In this setting it will be intractable to sufficiently sample a probabilistic space of dimension O(100), regardless of the surrogate modelling technique used. Therefore we will investigate means to efficiently create surrogate models with a moderate number of uncertain input parameters, under the assumption that the QoI is the output of a computationally expensive code. Of course, the term "moderate dimensionality" is system dependent. So to clarify, we consider the dimensionality moderate when the number of inputs parameters falls within the range of 5 to 10. Many problems of engineering interest are simulated using codes with similar dimensionality, e.g. turbulence models [START_REF] Edeling | Predictive rans simulations via Bayesian model-scenario averaging[END_REF][START_REF] Guillas | Bayesian calibration of the constants of the k-ε turbulence model for a cfd model of street canyon flow[END_REF], groundwater models [START_REF] Rojas | Application of a multimodel approach to account for conceptual model and scenario uncertainties in groundwater modelling[END_REF] or thermodynamic equations-of-state [START_REF] Cinnella | Sensitivity analysis of dense gas flow simulations to thermodynamic uncertainties[END_REF][START_REF] Merle | Bayesian quantification of thermodynamic uncertainties in dense gas flows[END_REF].

Due to its adaptivity, high polynomial order and Runge-phenomenon free interpolation, the SSC method requires a relatively low number of code samples to attain a given convergence level for problems with moderate dimensionality. For example, the SSC method has been applied in [START_REF] Axerio-Cilies | Extreme ensemble computation for optimization under uncertainty[END_REF] to optimization under uncertainty of a Formula 1 tire brake intake. After a one-dimensional perturbation analysis, 3 variables were selected for analysis. Furthermore, in [START_REF] Congedo | A simplex-based numerical framework for simple and efficient robust design optimization[END_REF][START_REF] Congedo | Simplex-simplex approach for robust design optimization[END_REF] the SSC method is efficiently coupled with Downhill-Simplex optimization in a setting for robust design optimization. Several example problems are considered, but again the maximum number of uncertain variables is 3. We found that the SSC method itself, without considering the cost of sampling the code, can become prohibitively expensive when considering 5 uncertain parameters. Furthermore, due to excessive memory requirements we were unable to create surrogate models of dimension 6 or higher. In [START_REF] Witteveen | Simplex stochastic collocation with random sampling and extrapolation for nonhypercube probability spaces[END_REF] the authors also note that the cost of the Delaunay triangulation becomes prohibitively large from a dimensionality of 5 onward. They suggested to replace the Delaunay triangulation with a scheme where simplices are formed by selecting the nearest points from randomly placed Monte Carlo (MC) samples. Using this approach the SSC method was applied to a continuous QoI with 15 uncertain parameters. However, in this case individual simplices can overlap and there is no guarantee that the entire probabilistic domain is covered.

This paper is aimed at identifying bottlenecks, and reducing the computational burden of the SSC method, while retaining the Delaunay triangulation. We investigate two separate techniques. First we propose the use of new alternative stencils based upon the set-covering problem [START_REF] Karp | Reducibility Among Combinatorial Problems[END_REF]. The main idea is to use the fast increase in number of simplex elements with polynomial order to create a small set of stencils which covers the entire probabilistic domain. Afterwards, only this set is used for interpolation. This allows for a more efficient implementation of the SSC method. Our results show that these stencils are capable of reducing the computational cost up to 8 dimensions. We furthermore present a new method for avoiding problems with the ill-conditioning of the sample matrix, and we provide a new formula for placing uniformly distributed samples in a simplex of arbitrary dimension. Secondly, inspired by the work of Ma et al. [START_REF] Ma | An adaptive high-dimensional stochastic model representation technique for the solution of stochastic partial differential equations[END_REF], we integrate the SSC method into the cut-HDMR framework. This approach combines the advantages of SSC and cut-HDMR, and avoids the disadvantages related to the ASG method such as linear interpolation and the possible occurrence of the Runge phenomenon. Unlike the authors of [START_REF] Ma | An adaptive high-dimensional stochastic model representation technique for the solution of stochastic partial differential equations[END_REF], we apply our method to a complex computer code for which obtaining samples is expensive. For both proposed techniques we perform a detailed analysis of the error and we give a discussion on computational cost as a function of the number of input parameters.

This paper is organized as follows: in Section 2 we present the baseline SSC method as developed by Witteveen et al. Next, in Section 3 we describe the Set-Covering stencils, our method for avoiding singular sample matrices and the analytic mapping for uniformly-distributed simplex sampling. The following section describes the cut-HDMR approach and in Section 5 we present the obtained results and the discussion. Finally, we give our conclusions in Section 6.

Simplex-stochastic collocation method

In this section we give a general outline of the Simplex-Stochastic Collocation (SSC) method as developed by Witteveen et al. For a more detailed description we refer to [START_REF] Witteveen | An adaptive stochastic finite elements approach based on Newton-Cotes quadrature in simplex elements[END_REF][START_REF] Witteveen | Simplex stochastic collocation with random sampling and extrapolation for nonhypercube probability spaces[END_REF][START_REF] Witteveen | Simplex stochastic collocation with eno-type stencil selection for robust uncertainty quantification[END_REF][START_REF] Witteveen | Refinement criteria for simplex stochastic collocation with local extremum diminishing robustness[END_REF].

General outline baseline SSC method

The SSC method was introduced as a non-intrusive method intended for robust and efficient propagation of uncertainty through computer models. It differs from traditional collocation methods, e.g. [START_REF] Babuška | A stochastic collocation method for elliptic partial differential equations with random input data[END_REF][START_REF] Loeven | Probabilistic collocation: an efficient non-intrusive approach for arbitrarily distributed parametric uncertainties[END_REF], in two main ways. First, for multidimensional problems it employs the Delaunay triangulation to discretize the probability space into simplex elements, rather than relying on the more common tensor product of one-dimensional abscissas [START_REF] Pecnik | Uncertainty quantification for laminar-turbulent transition prediction in rans turbomachinery applications[END_REF]. Using a multi-element technique has the advantage that mesh adaptation can be performed, such that only regions of interest are refined. Secondly, the SSC method is capable of handling non-hypercube probability spaces [START_REF] Witteveen | Simplex stochastic collocation with random sampling and extrapolation for nonhypercube probability spaces[END_REF].

The response surface of the QoI u(ξ ) is denoted by w(ξ ) and it is constructed using a set of n s samples from the com-

putational model, v = {v 1 , • • • , v n s }. Here, ξ is a vector of d random input parameters ξ (ω) = (ξ 1 (ω 1 ), • • • , ξ d (ω d )) ∈ ⊂ R d .
Furthermore, we define to be the parameter space and ω = (ω 1 , • • • , ω d ) ∈ ⊂ R d is a vector containing realizations of the probability space ( , F , P ) with F the σ -algebra of events and P a probability measure. The variables in ω are distributed uniformly as U (0, 1), and the input parameters can have any distribution f ξ , although for the sake of simplicity we restrict ourselves in this paper to f ξ = U (ξ a i , ξ b i ), with the bounds ξ a i and ξ b i . We perform all our analysis on the unit hypercube K d := [0, 1] d , and we use a linear map in order to go from K d to the parameter domain ξ . Our goal is to propagate f ξ through the computational model in order to assess the effect of f ξ on the m-th statistical moment of u (ξ (ω)), i.e.

we wish to compute μ (m) u := u (ξ ) m f ξ (ξ )dξ .

(

Note that u can also be a function of a physical coordinate x or other deterministic explanatory variables, but for brevity we omit x from the notation. Since the SSC method discretizes the parameter space into n e disjoint simplices

= 1 ∪ • • • ∪ n e , the mth statistical moment is approximated as μ (m) u = u (ξ ) m f ξ (ξ )dξ ≈ μ (m) w := n e j=1 j w j (ξ ) m f ξ (ξ ) dξ .
(

) 2
Here, w j is a local polynomial function of order p j associated with the j-th simplex j such that

w(ξ ) = w j (ξ ), for ξ ∈ j , (3)
and the interpolation condition requires

w j (ξ k j,l ) = v k j,l . (4)
The subscript k j,l ∈ {1, • • • , n s } is a global index which refers to the k-th added computational sample, while j refers to a certain simplex element. Furthermore, l = 0, • • • , N j is a local index used to count the number of samples from v involved in the construction of w j . The N j + 1 number of points needed for d-dimensional interpolation of order p j is given by

N j + 1 = (d + p j )! d!p j ! , ( 5
)
and the local interpolation function w j itself is given by the expansion

w j (ξ ) = N j l=0 c j,l j,l (ξ ). ( 6
)
The choice of basis polynomials j,l , and the determination of the interpolation coefficients c j,l is dealt with in Section 2.2.1. Note that for a given d, the maximum allowable order p j based on the number of samples n s can be inferred from [START_REF] Barber | The qhull software library[END_REF].

The particular choice of p j will depend on the smoothness of the response, with the objective of avoiding the Runge phenomenon.

Which N j + 1 points are used in ( 6) is determined by the interpolation stencil S j . The stencil can be constructed based on the nearest-neighbor principle [START_REF] Witteveen | Simplex stochastic collocation with random sampling and extrapolation for nonhypercube probability spaces[END_REF]. In this case the first d + 1 points are the vertices {ξ k j,0 , which would suffice for p j = 1. For higher degree interpolation, neighboring points ξ k are added based on their proximity to the center of simplex j , i.e. based on their ranking according to ξ kξ center, j 2 , [START_REF] Chung | On lattices admitting unique Lagrange interpolations[END_REF] where those ξ k of the current simplex j are excluded. The simplex center ξ center, j is defined as ξ center, j :=

• • • , ξ k j,d } of the simplex j ,
1 d + 1 d l=0 ξ k j,l . ( 8
)
The nearest neighbor stencil [START_REF] Chung | On lattices admitting unique Lagrange interpolations[END_REF] leads to a p j distribution that can increase quite slowly when moving away from a discontinuity. An example of this behavior can be found in Fig. 1(a), which shows the Delaunay triangulation with a discontinuity running through the domain. An alternative to nearest-neighbor stencils are stencils created according to the Essentially Non-Oscillatory (ENO) principle [START_REF] Witteveen | Simplex stochastic collocation with eno-type stencil selection for robust uncertainty quantification[END_REF]. The idea behind ENO stencils is to have higher degree interpolation stencils up to a thin layer of simplices containing the discontinuity. For a given simplex j , its ENO stencil is created by locating all the nearest-neighbor stencils that contain j , and subsequently selecting the one with the highest p j . This leads to a Delaunay triangulation which captures the discontinuity better than its nearest-neighbor counterpart. An example can be found in Fig. 1(b). Unless otherwise stated, for all subsequent baseline SSC surrogate models we will use ENO-type stencils.

The initial samples, at least in the case of hypercube probability spaces, are located at the 2 d corners of the hypercube K d . Furthermore, one sample is placed in the middle of the hypercube. Next, the initial grid is adaptively refined based on an appropriate error measure. This error measure can either be based on the hierarchical surplus between the response surface of the previous iteration and new a sample v k , or on the geometrical properties of the simplices. The latter option is more reliable in multiple stochastic dimensions as it is not based on the hierarchical surplus in a single discrete point [START_REF] Witteveen | Refinement criteria for simplex stochastic collocation with local extremum diminishing robustness[END_REF]. The geometrical refinement measure is given by

ē j := ¯ j ¯ 2O j j . ( 9
)
It contains the probability and the volume of simplex j , i.e.

¯ j = j f ξ (ξ )dξ and ¯ j = j dξ , [START_REF] Congedo | Simplex-simplex approach for robust design optimization[END_REF] where ¯ = n e j=1 ¯ j . The probability ¯ j can be computed by Monte-Carlo sampling and ¯ j via the relation

¯ j = 1 d! |det (D) |, D = ξ k j,1 -ξ k j,0 ξ k j,2 -ξ k j,0 • • • ξ k j,d+1 -ξ k j,0 ∈ R d×d . (11)
Finally, the order of convergence O j is given by [START_REF] Witteveen | Refinement criteria for simplex stochastic collocation with local extremum diminishing robustness[END_REF]

O j = p j + 1 d . ( 12
)
The simplex with the highest ē j is selected for refinement. To ensure a good spread of the sample points, only randomly- selected points inside a simplex sub-element sub j are allowed. The vertices of this sub-element are defined as ξ sub j,l :=

1 d d l * =0 l * =l ξ k j,l * , ( 13
)
see Fig. 2 for a two-dimensional example. In order to place random samples uniformly in an arbitrary simplex we derive an analytical map M d : K n ξ → j , see Section 2.2.2. The SSC algorithm can be parallelized by selecting the N < n e simplices with the N largest ē j for refinement.

Note that by using [START_REF] Foo | The multi-element probabilistic collocation method (me-pcm): error analysis and applications[END_REF] only simplex interiors will be refined (see again Fig. 2), and the boundaries of the hypercube will never be sampled outside the initial 2 d points. As a consequence, discontinuities that cross a hypercube border cannot be captured accurately at that border. To avoid this, we do not use [START_REF] Foo | The multi-element probabilistic collocation method (me-pcm): error analysis and applications[END_REF] if a simplex element located at the boundary is selected for refinement. Instead, we randomly place samples at the longest simplex edge which is at the boundary, ±10% from the edge center.

Note that ( 9) is probabilistically weighted through j and that it assigns high ē j to those simplices with low p j since in general ¯ j 1. Effectively this means that ( 9) is a solution-dependent refinement measure which refines simplices near discontinuities since the SSC method automatically reduces the polynomial order if a stencil S j crosses a discontinuity. It achieves this by enforcing the so-called Local Extremum Conserving (LEC) limiter to all simplices j in all S j . The LEC condition is given by min ξ ∈ j w j (ξ ) = min v j ∧ max ξ ∈ j w j (ξ ) = max v j , [START_REF] Ganapathysubramanian | Sparse grid collocation schemes for stochastic natural convection problems[END_REF] where v j = {v k j,0 , • • • , v k j,d } are the samples at the vertices of j . If w j violates [START_REF] Ganapathysubramanian | Sparse grid collocation schemes for stochastic natural convection problems[END_REF] in one of its j ∈ S j , the polynomial order p j of that stencil is reduced, usually by 1. Since polynomial overshoots occur when trying to interpolate a discontinuity, p j is reduced the most in discontinuous regions. Interpolating a function on a simplex with p j = 1 and v j located at its vertices always satisfies [START_REF] Ganapathysubramanian | Sparse grid collocation schemes for stochastic natural convection problems[END_REF] [START_REF] Witteveen | Refinement criteria for simplex stochastic collocation with local extremum diminishing robustness[END_REF]. This ensures that w(ξ ) is able to represent discontinuities without suffering from the Runge phenomenon. In practice, [START_REF] Ganapathysubramanian | Sparse grid collocation schemes for stochastic natural convection problems[END_REF] is enforced for all j in all S j via random sampling of the w j . If for a given w j [START_REF] Ganapathysubramanian | Sparse grid collocation schemes for stochastic natural convection problems[END_REF] is violated for any of the randomly placed samples ξ j , the polynomial order of the corresponding stencil is reduced.

Again, how we sample the d-dimensional simplices is described in Section 2.2.2. The computational cost of enforcing [START_REF] Ganapathysubramanian | Sparse grid collocation schemes for stochastic natural convection problems[END_REF] is investigated in Section 5.

The procedure of enforcing the LEC condition, computing a refinement measure and subsequently refining certain selected simplices is either repeated for a maximum of I iterations, n s max samples or halted when a sufficient level of accuracy is obtained. This level of accuracy can be estimated through an error measure based upon the hierarchical surplus [START_REF] Witteveen | Refinement criteria for simplex stochastic collocation with local extremum diminishing robustness[END_REF]. As mentioned, this is the difference between the response surface w j and the newly added code sample v k j,ref at the refinement location ξ k j,ref , i.e.

(ξ k j,ref ) := w j (ξ k j,ref ) -v k j,ref . ( 15
)
This is a point estimate of the error, located at what will be a vertex in the new refined Delaunay grid. To assign error estimates to the simplices rather than to vertices, the error ˜ j is introduced. For each j , ˜ j is simply the absolute value of (15) of its most recently added vertex ξ k * . Since adding vertices will change the Delaunay discretization we relate the error of the previous simplex to the new one via [START_REF] Witteveen | Simplex stochastic collocation with random sampling and extrapolation for nonhypercube probability spaces[END_REF]. The ratio ¯ i / ¯ k * ,ref represents the change in volume from its old size ¯ k * ,ref , i.e. the volume of the simplex which was refined by ξ k * , to its new size ¯ j . Finally, each individual ˆ j is combined in a global error estimate via the following root mean square (RMS) error norm

ˆ j ≈ ˜ j ¯ j ¯ k * ,ref O j (16)
ˆ rms = n e j=1 j ˆ 2 j . ( 17
)
The complete baseline SSC method is given in pseudo code in Appendix A.

Improvements on the baseline SSC method

Before discussing our new stencil selection technique in Section 3.1, we introduce two improvements to the baseline SSC method not discussed in the original references [START_REF] Witteveen | An adaptive stochastic finite elements approach based on Newton-Cotes quadrature in simplex elements[END_REF][START_REF] Witteveen | Simplex stochastic collocation with random sampling and extrapolation for nonhypercube probability spaces[END_REF][START_REF] Witteveen | Simplex stochastic collocation with eno-type stencil selection for robust uncertainty quantification[END_REF][START_REF] Witteveen | Refinement criteria for simplex stochastic collocation with local extremum diminishing robustness[END_REF].

Poised sample sequence

The authors of [START_REF] Witteveen | Simplex elements stochastic collocation in higher-dimensional probability spaces[END_REF] write [START_REF] Boissonnat | Incremental construction of the Delaunay triangulation and the Delaunay graph in medium dimension[END_REF] in matrix form, constraining j,l to the class of monomials, and subsequently solve explicitly for the coefficients c j,l . They note that although they had no difficulties in solving this system, the matrix could have a high condition number. This poses no real problem for d ≤ 3, but for higher dimensions it can become problematic. To cope with this we impose an additional condition on the construction of the stencils S j such that the interpolation problem is poised, meaning that the sample matrix is non-singular [START_REF] Olver | On multivariate interpolation[END_REF]. In the following discussion we drop the subscript j until further notice to make the notation more concise.

To construct the interpolating monomials, let us define the collection consisting of

N + 1 d-dimensional multi-indices ī := (i 1 , • • • , i k , • • • , i d )
, where for all ī we have | ī| := i 1 +• • •+i d ≤ p j and each i k is an integer between 0 and d. Furthermore, for a given vertex ξ l = (ξ 1,l , • • • , ξ d,l ) belonging to stencil S, let us define its ī-th power to be ξ ī

l := ξ i 1 1,l × • • • × ξ i d d,l
. The sample matrix , a multi-dimensional Vandermonde matrix, can then be written as

= ⎡ ⎢ ⎢ ⎢ ⎣ ξ 0 0 ξ 1 0 • • • ξ N 0 ξ 0 1 ξ 1 1 • • • ξ N 1 . . . . . . . . . ξ 0 N ξ 1 N • • • ξ N N ⎤ ⎥ ⎥ ⎥ ⎦ ∈ R (N+1)×(N+1) . ( 18
)
As an example, the l-th row of ( 18) in lexicographical order for p j = 2 will look like

[ 1 ξ 1,l ξ 2,l ξ 2 1,l ξ 1,l ξ 2,l ξ 2 2,l ]
. The coefficients c l in ( 6) can now be obtained by solving the system

⎡ ⎢ ⎢ ⎢ ⎣ ξ 0 0 ξ 1 0 • • • ξ N 0 ξ 0 1 ξ 1 1 • • • ξ N 1 . . . . . . . . . ξ 0 N ξ 1 N • • • ξ N N ⎤ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎣ c 0 c 1 . . . c N ⎤ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎣ v 0 v 1 . . . v N ⎤ ⎥ ⎥ ⎦ , ( 19
)
where {v 0 , • • • , v N } are the code samples belonging to stencil S. Once the c l are known, we can interpolate to any point ξ in the domain spanned by S.

We define ≡ det ( ), and note that the whole approach hinges on the well-poisedness condition = 0. This condition is relatively easy violated during the SSC procedure in higher dimensions. For instance, if for d = 4 we determine the maximum allowable p using (5) on the initial Delaunay grid we obtain p max = 2. However, many stencils in this case will have = 0. Also situations where a stencil has too many vertices located in the same plane (e.g. due to edge refinement at the boundary of K d ), can lead to a zero determinant of [START_REF] Loeven | Probabilistic collocation: an efficient non-intrusive approach for arbitrarily distributed parametric uncertainties[END_REF]. Thus, for d > 1 the poisedness condition = 0 imposes constraints on the geometrical distribution of the ξ l . From [START_REF] Olver | On multivariate interpolation[END_REF][START_REF] Chung | On lattices admitting unique Lagrange interpolations[END_REF] we know

Theorem 1. The N + 1 vertices ξ 0 , • • • , ξ N ∈ R d are
polynomially poised iff they are not a subset of any algebraic hypersurface of degree ≤ p.

An algebraic hypersurface in R d is a d -1 dimensional surface embedded in a d-dimensional space constrained to satisfy an equation f (ξ 1 , • • • , ξ d ) = 0.
The degree is given by f . The authors of [START_REF] Chung | On lattices admitting unique Lagrange interpolations[END_REF] devised a Geometric Characterization (GC) condition which allows us to detect if a set of vertices is poised, i.e.: Definition 1. GC condition: For each ξ l in a set of N + 1 vertices in R d , there exists p distinct hyperplanes G 1,l , • • • , G p,l such that i) ξ l does not lie on any of these planes, and ii) all other ξ k , k = {0, • • • , N}\{l} lie on at least one of these hyperplanes.

Mathematically speaking i) and ii) amount to

ξ i ∈ p k=0 G k,l if i = l, ∀i = 1, 2, • • • , N. ( 20
)
Fig. 3. When selecting node ξ 1 , there exists one (p = 1) plane which contains all other points except ξ 1 . This is true for all nodes in the simplex.

Theorem 2. Let {ξ l } be a set of N + 1 vertices in R d . If {ξ l } satisfies the GC condition, then {ξ l } admits a unique interpolation of degree ≤ p [START_REF] Chung | On lattices admitting unique Lagrange interpolations[END_REF].

Due to its geometrical configuration, a single simplex j in R d always satisfies the GC condition for p = 1, see Fig. 3 for a three-dimensional example. For a given vertex ξ l ∈ j , we always have one hyperplane containing the face of the simplex made up by all vertices except ξ l . Thus, Theorem 2 implies that simplex j will lead to a with = 0 and p j = 1. We use this result to obtain a set of well-poised ENO stencils S j ∀ j = 1, • • • , n e , in a way that is similar to the construction of the ENO-stencils as described in [START_REF] Witteveen | Simplex stochastic collocation with eno-type stencil selection for robust uncertainty quantification[END_REF]. Only if during the enforcement of the LEC condition ( 14) we encounter a stencil S j for which = 0, we collect a set of k candidate nearest-neighbor stencils {S j,i } k i=1 which all contain simplex j . We then select the S j which has the highest p j and = 0. In the worst case scenario we get p j = 1, where S j contains only the vertices of j itself and for which = 0 is guaranteed by Theorem 2. If we have multiple S j with p j > 1 which satisfy these conditions, we select the one with the smallest average Euclidean distance to the cell-center ξ center, j .

Simplex sampling

Simplices are refined by randomly placing a point inside the sub-simplex [START_REF] Foo | The multi-element probabilistic collocation method (me-pcm): error analysis and applications[END_REF]. Also, to randomly sample the w j during the LEC enforcement we need to place random points inside d-dimensional simplices. If we would like to uniformly sample a line section with the end points [ξ 0 , ξ 1 ] we would use the mapping

M 1 = ξ 0 + r 1 (ξ 1 -ξ 0 ), ( 21
)
where r 1 ∼ U [0, 1]. Generating points inside a triangle can be done with

M 2 = ξ 0 + r 1/2 2 (ξ 1 -ξ 0 ) + r 1/2 2 r 1 (ξ 2 -ξ 1 ) (22)
which maps points {r 1 , r 2 } inside the unit square K 2 to points inside a triangle described by the vertices {ξ 0 , ξ 1 , ξ 2 } [START_REF] Turk | Graphics Gems[END_REF].

The working principle of ( 22) is shown in Fig. 4(a). The parameter r 1/2 2 selects a line segment parallel to the edge [ξ 0 , ξ 1 ],

while r 1 selects a point along the chosen line segment. The exponent 1/2 ensures that uniformly distributed points in the square yield uniformly distributed points in the triangle. This can be shown by considering the length of the chosen line segment, which increases linearly when r 1/2 2 moves from ξ 0 to ξ 1 . Since we require a uniform distribution of points, and considering

r 1 ∼ U [0, 1], the pdf of r 1/2 2
should be linear as well. If we have the random variable X = r 1/τ with r ∼ U [0, 1] and τ ∈ N >0 , we find the cumulative distribution function (cdf) of X as

F X (x) = P(X ≤ x) = P r 1/τ ≤ x = P r ≤ x τ = x τ . ( 23
)
And thus we have the pdf f X (x) = dF X /dx = τ x τ -1 ∼ Beta(τ , 1). Therefore, in order to have a linear pdf for r 1/τ , we must set τ = 2.

It is suggested in [START_REF] Turk | Graphics Gems[END_REF] that [START_REF] Nobile | A sparse grid stochastic collocation method for partial differential equations with random input data[END_REF] can be extended to higher dimensions, although no specific formulas are given. Hence, we use the same principle to select uniformly distributed points inside a tetrahedron, see Fig. 4(b). Here, the parameter r 1/3 3 selects a triangle parallel to the base of the tetrahedron. From there we use r 1/2 2 and r 1 as before to select a point on this triangle. The exponent 1/3 again ensures that the point distribution will be uniform. Note that the area of the selected triangles increases quadratically as r 1/3 3 moves from ξ 0 to ξ 1 . Hence, it must be distributed as Beta [START_REF] Axerio-Cilies | Extreme ensemble computation for optimization under uncertainty[END_REF][START_REF] Amenta | Complexity of Delaunay triangulation for points on lower-dimensional polyhedra[END_REF]. We can now derive an expression for M 3 using the geometrical similarities between the base triangle and the selected parallel triangle, which gives us

M 3 = ξ 0 + r 1/3 3 (ξ 1 -ξ 0 ) + r 1/3 3 r 1/2 2 (ξ 2 -ξ 1 ) + r 1/3 3 r 1/2 2 r 1 (ξ 3 -ξ 2 ). ( 24
)
When comparing [START_REF] Merle | Bayesian quantification of thermodynamic uncertainties in dense gas flows[END_REF], [START_REF] Nobile | A sparse grid stochastic collocation method for partial differential equations with random input data[END_REF] and [START_REF] Pecnik | Uncertainty quantification for laminar-turbulent transition prediction in rans turbomachinery applications[END_REF] we see a pattern emerge which suggests that the map from a d-dimensional

hypercube to a d-dimensional simplex with vertices {ξ 0 , • • • , ξ d } in R d and uniform point distribution is
M d = ξ 0 + d i=1 i j=1 r 1 d-j+1 d-j+1 (ξ i -ξ i-1
), [START_REF] Pini | Hybrid adjoint-based robust optimization approach for fluid-dynamics problems[END_REF] where again the r q are distributed as U [0, 1]. Our proof that (25) produces uniformly distributed samples in the simplex can be found in Appendix D.

To numerically test (25) in 2 and 3 dimensions we can simply plot samples points, an example of which can be found in Fig. 5. We have performed similar tests up to 8 dimensions.

SSC Set-Covering method

In this section we describe alternative interpolation stencils, which results in a computational speed up in higher dimensions.

Set covering stencils

Section 5 will show that the enforcement of the LEC condition can become computationally expensive for high d and p j . This is especially true for smooth response surfaces of the QoI. For many stencils of our discontinuous problem, the LEC condition is violated and p j is reduced which in turn significantly lowers the total required number of surrogate model evaluations (n w ) needed to check [START_REF] Ganapathysubramanian | Sparse grid collocation schemes for stochastic natural convection problems[END_REF]. This does not happen very often when the response surface is smooth. As a consequence of the exponential nature of n w (see Section 5.1.1), we also see an exponential increase in the computation time needed to construct the surrogate model. Note that this increase is due to the SSC procedure, and thus is additional to the time needed to sample the computer code.

However, the problem lies not only with the exponential increase of n w , but also in the extremely large overlap of the stencils S j . Note that the baseline SSC method enforces the LEC condition for all simplices j in all stencils S j . Hence, in each simplex j , w j is evaluated the same number of times as j appears in all stencils S j . For a two-dimensional example see Fig. 6. There are two stencils, denote them S r and S q , associated to two different simplices r and q . The dark colored simplices are the ones which appear in both stencils. Thus, when the LEC condition is checked for both stencils, w r but also w q is evaluated in the dark simplices. Moreover, since there are n e stencils, the overlap will be large, and many different w j will be evaluated in the same simplex element. This is no bottleneck for problems of low-dimensionality, but if the dimension increases this overlap will make the LEC condition very costly to enforce, see Section 5.1.1.

We propose an alternative technique for problems with higher d, using Set-Covering (SC) stencils based on the well-know set-covering problem [START_REF] Karp | Reducibility Among Combinatorial Problems[END_REF], stated as follows in SSC terminology: Set Covering problem. Let X j = { j,0 , • • • , j,K } be the set of all simplices that are inside the domain spanned by the vertices of stencil S j . Then, given the set X = {X 1 , • • • , X n e }, and the set of all simplices U = { 1 , • • • , n e }, find the smallest subset C ⊆ X that covers U , i.e. for which U ⊆

X j ∈C X j holds.
It is shown in [START_REF] Karp | Reducibility Among Combinatorial Problems[END_REF] that the set-covering problem is NP-complete, and thus no fast solution is known. We could approximate C by the greedy algorithm, which at each step simply selects the X j with the largest number of uncovered simplices. We then would have to check the LEC condition for all stencils in S sc , defined as the set of S j corresponding to the X j ∈ C. For (high-dimensional) problems with a maximum polynomial order p max > 1, the number of stencils in S sc will be significantly lower than n e . However, this approach would still require to construct all X j ∈ X . Also, many of the X j could potentially cross a discontinuity, leading to a violation of the LEC condition and the subsequent reduction in size of X j .

When this happens the SC property of C can no longer be guaranteed. Thus, an iterative approach would be necessary which runs until S sc satisfies both the SC and LEC property.

For reasons of computational efficiency, we want to avoid this iterative approach as much as possible, and thus not rely completely on the LEC condition to turn a set of nearest-neighbor stencils into a set of ENO stencils. Hence we will use the information contained in v regarding the discontinuity location to create a small set of SC stencils that also resemble ENO stencils, i.e. which do not cross a discontinuity. We will denote these stencils as SCENO stencils. Although more sophisticated approaches are available [START_REF] Witteveen | Subcell resolution in simplex stochastic collocation for spatial discontinuities[END_REF], for reasons of simplicity we identify the j through which the discontinuity runs by simply imposing a threshold v t on the maximum jump observed in v at each simplex. Then, the set of discontinuous simplices can be defined as

D = { j | | max v k j,l -min v k j,l | ≥ v t , l = 0, • • • , d, j = 1, • • • , n e } ( 26
)
For the nozzle flow case we set the threshold value to v t = 1.0. A two and three-dimensional visualization of the j ∈ D can be found in Fig. 7. We furthermore redefine the set C as the set containing all the simplices j that are currently covered by a stencil S j , rather than the true smallest subset C ⊆ X of the SC problem.

The general outline for constructing the SCENO stencils is now as follows. For the j ∈ D we set p j = 1 and C = C ∪ D,

i.e. we add all discontinuous simplices to the set of covered simplices. Next, we specify the initial simplex * j as the simplex from the set U \ C with the largest volume. For the selected simplex we grow its stencil by adding neighboring j which are not covered yet, i.e. which are not in C. This will yield a set C where every simplex appears only once, i.e. a set with zero overlap. Note that to relax this condition one can easily allow for the addition of neighboring simplices which are in C \ D. In either case we continue growing the stencil until there are no more available neighbors or until S j is large enough to allow interpolation of order p max . We then move to the next * j and repeat until C covers the entire probabilistic space U . For a graphical representation of the stencil construction we refer to Fig. 8. It is important to note that our main goal is to find a set C with a cardinality |C| significantly less that n e , which is an easier task than approximating the true minimal C of the SC problem as closely as possible. In Appendix B the algorithm for constructing the SCENO stencils is displayed in pseudo code.

This approach assures that we have a relatively small set S sc for which: i) |S sc | n e , ii) that not all n e nearest-neighbor X i ∈ X need be calculated, iii) that no X j crosses a discontinuity, and iv) the j ∈ D are interpolated linearly. The result is that the number of times the LEC condition needs to be checked is reduced significantly. Only for those S j associated to the X j ∈ C \ D it is still necessary to check for interpolation overshoots, since the j ∈ D are guaranteed to be LEC due to their linear interpolation. The property of SCENO stencils mentioned under iii) also means that the number of times the LEC condition is violated is reduced, although not always to zero due to reasons of ill conditioning of the sample matrix [START_REF] Loeven | Probabilistic collocation: an efficient non-intrusive approach for arbitrarily distributed parametric uncertainties[END_REF]. This is especially true for high d. An approach as described in Section 2.2.1 would render some of the advantages mentioned under i)-iv) void. Reducing p j for ill-conditioned stencils will increase the cardinality of S sc , and all X j ∈ X should be calculated in order to look for alternative stencils. Instead we directly solve an ill-conditioned system [START_REF] Ma | An adaptive hierarchical sparse grid collocation algorithm for the solution of stochastic differential equations[END_REF] in the non-null subspace of the solution as described in [START_REF] Kozaczek | A direct algorithm for solving ill-conditioned linear algebraic systems[END_REF]. This method utilizes Gauss-Jordan elimination with complete pivoting to identify the null subspace of a singular matrix , i.e. null c null = 0. This partitions the linear system as depicted below, range

• • • • • • null c range c null = v . . . , ( 27
)
where range c range = v is the non-null subspace in which we can obtain accurate solutions. In the case of an ill-conditioned system, the null subspace is closely approximated by a space where the pivots ψ ii are very small but not exactly equal to zero. The start of this 'near-null' subspace is identified by the first pivot ψ ii for which the condition |ψ ii /η c | < holds, where η c is the largest pivot of and is a very small parameter, which we set equal to 10 -14 . In both the ill-conditioned and singular case the detrimental effect of null on the solution is eliminated by a so-called zeroing operation, which basically replaces null by an identity matrix of equal dimension and sets c null = 0. Thus, effectively speaking those coefficients c j,l which have been overwhelmed by round-off error are automatically cut out of the expansion [START_REF] Boissonnat | Incremental construction of the Delaunay triangulation and the Delaunay graph in medium dimension[END_REF]. In our experiments we found that the dimension of null , i.e. the nullity of , is small compared to the dimension of the full , see Table 1 for some typical examples at d = 6.

If the system of equations is well-posed, the algorithm amounts to regular Gauss-Jordan elimination with complete pivoting. In any case, the quality of the response surface is checked via the LEC condition.

High-Dimensional Model-Reduction

As will be shown in Section 5, the use of SCENO stencils makes the SSC method more computationally efficient within the range of dimensions where a Delaunay triangulation can be made. For problems of higher dimensionality a different approach is required. In physical systems it is often found that only a few parameters are influential, and only loworder correlations between the input parameters have a significant impact on the output. To capitalize on this behavior, High-Dimensional Model-Reduction techniques can be applied, see the references of Rabitz and Ali ¸s [START_REF] Rabitz | Efficient input-output model representations[END_REF][START_REF] Rabitz | General foundations of high-dimensional model representations[END_REF]. Our QoI is represented by a d-dimensional function f (ξ , x) defined on the hypercube K d , where x is a possible physical coordinate which we will again omit from the notation for the sake of brevity. Then, the HDMR expansion is an exact and finite hierarchical expansion of component functions of increasing dimension, given by

f (ξ ) = f 0 + i f i (ξ i ) + i 1 <i 2 f i 1 i 2 (ξ i 1 , ξ i 2 ) + • • • + i 1 <•••<i l f i 1 •••i l (ξ i 1 , • • • , ξ i l ) + • • • + f 1•••d (ξ i 1 , • • • , ξ i d ). ( 28
)
Here, the i 1 ,

• • • , i d are integers satisfying 1 ≤ i 1 < i 2 < • • • < i d ≤ d.
The zero-th order component function f 0 is a constant and represents the mean effect. The first-order function f i (ξ i ) is a univariate function, generally nonlinear, which represents the effect of independently varying input parameter ξ i . Higher order functions represent the cooperative effects of increasing number of variables acting together on the output. If high-order correlations are weak, the physical system f (ξ ) can be ef-ficiently represented by a truncated L-th order expansion, where L < d. This a called a problem with low effective dimension, which occurs frequently in problems of physical nature [START_REF] Foo | Multi-element probabilistic collocation method in high dimensions[END_REF]. Thus, the general idea is to solve multiple low-dimensional subproblems in place of a single high-dimensional one. The resultant computational effort to determine the component functions will scale polynomially, rather than the traditional exponential increase with d [START_REF] Rabitz | General foundations of high-dimensional model representations[END_REF].

A measure μ for the measure space (K d , B(K n ), μ), where B is the Borel σ -algebra on K d , is defined as

dμ(ξ ) := dμ(ξ 1 , • • • , ξ d ) = d i=1 dμ i (ξ i ), K 1 dμ i (ξ i ) = 1, dμ(ξ ) = g(ξ )dξ = d i=1 g i (ξ i )dξ i . ( 29
)
Here, g(ξ i ) is the marginal density of the input ξ i . It is the particular form chosen for the g i (ξ i ) that will determine the form of the component functions. In order to compute these functions, let us also define unconditional and conditional mean with respect to a group of input variables as

M f (ξ ) := K d f (ξ )dμ, M (i 1 •••i l ) f (ξ ) := K d-l f (ξ ) ⎡ ⎣ j / ∈{i 1 •••i l } dμ j (ξ j ) ⎤ ⎦ . ( 30
)
Then, via a family of projection operators

P i 1 •••i l : K d → K l
, the component functions are recursively defined as follows [START_REF] Rabitz | General foundations of high-dimensional model representations[END_REF]:

f 0 := P 0 f (ξ ) = M f (ξ ) f i (ξ i ) := P i f (ξ ) = M (i) f (ξ ) -P 0 f (ξ ) f ij (ξ i , ξ j ) := P ij f (ξ ) = M (ij) f (ξ ) -P i f (ξ ) -P j f (ξ ) -P 0 f (ξ ) . . . f i 1 •••i l (ξ ) := P i 1 •••i l f (ξ ) = M (i 1 •••i l ) f (ξ ) - j 1 <•••< j l-1 ⊂{i 1 •••i l } P j 1 ••• j l-1 f (ξ ) -• • • -P 0 f (ξ ) (31)
The component functions f i 1 ,•••i l and f j 1 ••• j k are independent and orthogonal, thus as long as one index between {i 1 ,

• • • i l } and { j 1 • • • j k } differs we have K d f i 1 ,•••i l (ξ i 1 , • • • , ξ i l ) f j 1 ••• j k (ξ j 1 , • • • , ξ j k )dμ = 0 (32)
The correlation interpretation of f i 1 •••i l is associated with the chosen form of the measure μ. [START_REF] Rojas | Application of a multimodel approach to account for conceptual model and scenario uncertainties in groundwater modelling[END_REF] together with [START_REF] Sobol | Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates[END_REF] becomes the well-know Analysis Of Variance (ANOVA) decomposition. Computing the component functions in the ANOVA decomposition involves evaluating multi-dimensional integrals, which can be done by for instance MC techniques [START_REF] Sobol | Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates[END_REF]. An alternative which is more computationally tractable is the cut-HDMR decomposition proposed in [START_REF] Rabitz | Efficient input-output model representations[END_REF][START_REF] Rabitz | General foundations of high-dimensional model representations[END_REF]. In this case the measure is defined as

If g i = 1, i = 1, • • • , d, the Lebesgue measure (dμ = dξ 1 dξ 2 • • • dξ d ) is retrieved and
dμ = d i=1 δ(ξ i -η i )dξ i , (33) i.e. g i (ξ i ) = δ(ξ i -η i ), a Dirac measure located at the 'cut center' η = (η 1 , η 2 , • • • , η d ).
This choice removes the need for evaluating multi-dimensional integrals, and it expresses f (ξ ) as a superposition of its values along lines, planes and hyperplanes passing through the cut center η. The component functions [START_REF] Sobol | Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates[END_REF] now become

f 0 := P 0 f (ξ ) = f (η) f i (ξ i ) := P i f (ξ ) = f (i) (ξ i ) -P 0 f (ξ ) f ij (ξ i , ξ j ) := P ij f (ξ ) = f (ij) (ξ i , ξ j ) -P i f (ξ ) -P j f (ξ ) -P 0 f (ξ ) . . . f i 1 •••i l (ξ ) := P i 1 •••i l f (ξ ) = f (i 1 •••i l ) (ξ i 1 , • • • , ξ i l ) - j 1 <•••< j l-1 ⊂{i 1 •••i l } P j 1 ••• j l-1 f (ξ ) -• • • -P 0 f (ξ ). ( 34
)
Here,

f (i 1 •••i l ) (ξ i 1 , • • • , ξ i l )
is the conditional mean [START_REF] Smolyak | Quadrature and interpolation formulas for tensor products of certain classes of functions[END_REF] taken with respect to measure [START_REF] Turk | Graphics Gems[END_REF], and thus it equals f with its inputs ξ i set to η i , except inputs ξ i 1 , • • • , ξ i l . As an example, consider the univariate function

f (i) (ξ i ) = f (η 1 , • • • , η i-1 , ξ i , η i+1 , • • • , η n ξ ).
The authors of [START_REF] Ma | An adaptive high-dimensional stochastic model representation technique for the solution of stochastic partial differential equations[END_REF] used the cut-HDMR framework coupled with their Adaptive Sparse-Grid (ASG) collocation method [START_REF] Ma | An adaptive hierarchical sparse grid collocation algorithm for the solution of stochastic differential equations[END_REF], where they chose η as the mean of the random input vector. Besides truncating (28) at a certain order, they also made their approach dimension adaptive based on weights which identify the important dimensions. Although their ASG method uses only a linear finite-element basis, interpolation overshoots can still occur. Thus, motivated by their work in [START_REF] Ma | An adaptive high-dimensional stochastic model representation technique for the solution of stochastic partial differential equations[END_REF] we will also employ a dimension adaptive cut-HDMR approach, except we will couple it with the SSC method utilizing the SCENO stencils to avoid the mentioned downsides of ASG.

If we define K := {1, 2, • • • , d}, the HDMR expansion (28) can be written in short-hand notation as [START_REF] Ma | An adaptive high-dimensional stochastic model representation technique for the solution of stochastic partial differential equations[END_REF] f (ξ )

= u⊆K f u (ξ u ) = u⊆K v⊆u (-1) |u|-|v| f (v) (ξ v ), ( 35
)
where in the first equality we sum over the powerset of K, i.e. over all possible subsets u ⊆ K. We furthermore set f ∅ = f 0 . The second equality is obtained by expanding each component function f u (ξ u ) as indicated in [START_REF] Wilcox | Turbulence Modeling for CFD[END_REF]. Notation wise, if for instance v = {1, 4, 6}, then

f (v) (ξ v ) = f (146) (ξ 1 , ξ 4 , ξ 6 ). Each individual |v|-dimensional subproblem f (v) (ξ v ) can be
approximated by a SSC surrogate [START_REF] Boissonnat | Incremental construction of the Delaunay triangulation and the Delaunay graph in medium dimension[END_REF]. In that case [START_REF] Witteveen | Simplex elements stochastic collocation in higher-dimensional probability spaces[END_REF] becomes

f (ξ ) ≈ w(ξ ) = u⊆K v⊆u (-1) |u|-|v| n e j=1 N j l=0 c jl jl (ξ v ). ( 36
)
In order to assess the convergence of each individual f (v) (ξ v ), the authors of [START_REF] Ma | An adaptive high-dimensional stochastic model representation technique for the solution of stochastic partial differential equations[END_REF] use the hierarchical surplus. This is also possible in the case of the SSC method, see [START_REF] Guillas | Bayesian calibration of the constants of the k-ε turbulence model for a cfd model of street canyon flow[END_REF]. Alternatively, the RMS error estimate ( 17) can used for this purpose. Since ( 17) is a global error estimate and it also includes information from the distribution of the input parameters, we use the RMS error to assess the convergence. Furthermore, the mean of each component function, defined as J u , can also be computed from the surrogate model

J u = v⊆u (-1) |u|-|v| n e j=1 N j l=0 c jl E jl (ξ v ) . ( 37
)
We compute (37) via random sampling, which can be performed quickly since it requires only sampling the surrogate model.

In order to identify the important dimensions, all first order component functions f i (ξ i ) are computed. Again, these are one-dimensional functions which measure the impact of a single independent input parameter on the output. Next, a weight is defined

α i = J i 2 f 0 2 , ( 38
)
which measures the contribution of each individual ξ i on the mean of all first order component functions [START_REF] Ma | An adaptive high-dimensional stochastic model representation technique for the solution of stochastic partial differential equations[END_REF]. We always take the L 2 norm • 2 over the spatial domain. Equation ( 38) can be considered as a sensitivity index, and only those dimensions for which [START_REF] Witteveen | Simplex stochastic collocation with random sampling and extrapolation for nonhypercube probability spaces[END_REF] is larger than a user-prescribed error threshold 1 are considered important. All higher order f v (ξ v ) where v contains indices of dimensions which did not make the cut will not be computed. Consider e.g. a d-dimensional problem on K d , where only v = {1} and v = {2} satisfy α i > 1 . The only higher-order component function that will be computed in this case is f 12 (ξ 1 , ξ 2 ), regardless of the value of d.

The downside of ( 38) is that it is hard to choose 1 beforehand. One should first create the first-order HDMR expansion and decide on an appropriate value a posteriori. An alternative is to use a weight measuring the relative contribution of J i with respect to the sum of all first-order means, i.e.

α i = J i 2 d k=1 J k 2 . ( 39
)
Now one can a priori choose a 1 ∈ [0, 1], and select the smallest set of important dimensions for which the sum of their α i is greater than 1 . Dimension adaptivity is extended to higher dimensions as well by defining a weight for |u| > 1 as [START_REF] Ma | An adaptive high-dimensional stochastic model representation technique for the solution of stochastic partial differential equations[END_REF]

α u = J u 2 v∈V comp ,|v|<|u-1| J v 2 . ( 40
)
Here, the set V comp simply holds all the indices v that were computed. Furthermore, all subsets v of component functions which are important are added to a set V imp . That way, a higher-order important u is admissible if all v ⊂ u required to compute [START_REF] Witteveen | Simplex elements stochastic collocation in higher-dimensional probability spaces[END_REF] are also in V imp . This is the so-called admissibility condition, which is given by Similar to the first-order case, we can define a relative counterpart of [START_REF] Witteveen | Subcell resolution in simplex stochastic collocation for spatial discontinuities[END_REF] as

u ∈ V imp and v ⊂ u ⇒ v ∈ V imp . (41)
α u = J u 2 v∈V comp ,|v|=|u| J v 2 , ( 42
)
such that the α u sum to one and we can choose a 1 ∈ [0, 1] a priori.

Finally, a relative error measure between two HDMR expansions of consecutive orders p -1 and p is defined as

α p = |u|≤p J u -|u|≤p-1 J u 2 |u|≤p-1 J u 2 . ( 43
)
The algorithm stops when α p becomes smaller than another used-defined threshold 2 . An overview of the HDMR algorithm is depicted in Appendix C.

Results and discussion

Comparison ENO-SCENO stencils

We present the results obtained with the baseline SSC method with ENO stencils, versus the SSC method with the SCENO stencils. As a test case we use a quasi-1D nozzle case up to 5 dimensions and an algebraic test function up to d = 8.

Nozzle flow

As a test case we use the solver from [START_REF] Pini | Hybrid adjoint-based robust optimization approach for fluid-dynamics problems[END_REF], which computes the flow through a quasi-1D diverging nozzle. We prescribe the flow to be sonic at the nozzle inlet, i.e. M in = 1. From fluid mechanics we know that the flow is driven by the pressure ratio, i.e. by the ratio between the total pressure p t at the inlet and the static pressure p of the surroundings at the nozzle exit. Depending on the value of p t /p, the flow can show very different behavior. If p t /p exactly equals the adaptation value, the flow reaches the static pressure of the surroundings at the nozzle exit and the jet exhausts smoothly into the atmosphere. A stronger p t /p will result in smooth flow through the nozzle, which is supersonic at the nozzle exit. In order to match the outside pressure p, the flow undergoes a supersonic expansion attached to the nozzle exit (under-expanded nozzle). A smaller p t /p, but still above a threshold that depends on the ratio of the exit to the throat area, still results in smooth flow through the nozzle, but this is now over-expanded and is compressed to the outside pressure through an oblique shock attached to the nozzle exit. When p t /p is equal to the threshold value, the flow is characterized by a normal shock located at the nozzle exit: upstream of the shock, the flow is smooth, and verifies adaptation conditions in the exit section; immediately downstream of it, the flow is subsonic and matches the outside pressure. Finally, when p t /p is below the threshold value, a normal shock wave is formed somewhere inside the nozzle. This results in subsonic flow at the exit, and an exit pressure that is equal to p [START_REF] Anderson | Fundamentals of Aerodynamics[END_REF]. Given the pressure ratio, the flow is completely characterized by the shape of the nozzle [START_REF] Anderson | Fundamentals of Aerodynamics[END_REF]. As in [START_REF] Pini | Hybrid adjoint-based robust optimization approach for fluid-dynamics problems[END_REF], we consider the following hyperbolic tangent for the nozzle shape

f (x) = a + b tanh (cx -d) . ( 44
)
To test the SSC method, we specify two different ranges for the uncertain parameters such that two radically different response surfaces have to be created. First, we prescribe a wider range for p such that the QoI is highly discontinuous, see Fig. 9. In the second case we restrict p to a more narrow interval such that the QoI is smooth. More specifically, we prescribe the uniform input distributions for the 6 uncertain parameters described in Table 2. Furthermore, we choose M out (the Mach numbers at the nozzle exit) as our quantity of interest, as it allows us to easy calculate other flow quantities via the isentropic relations once M out is known [START_REF] Anderson | Fundamentals of Aerodynamics[END_REF]. When constructing the surrogate models, we will use a linear transformation for each input to map points from [0, 1] in the stochastic domain to points in the physical domain with the range as specified in Table 2. This simplifies the construction of the surrogate models as it allows us to always work in the standard d-dimensional hypercube K d .

For now, we will consider just the first 5 uncertain parameters of Table 2. In Tables 3 and4 we show the computation time T in minutes versus the dimension d, in case of the discontinuous and smooth QoI for both the baseline and the method based on SCENO stencils. This is of course dependent upon the available computational resources, in our case a 24 core workstation. We use these cores to parallelize the LEC condition, code sampling and ENO stencils. Our algorithm for the construction of the SCENO stencils is not implemented in parallel, and uses just 1 core. We can see that T rises very quickly as d increases in the case of the baseline method, especially in the case of the smooth QoI. To explain which element is responsible for the high computation time, we also show the percentage of T that is spent on the LEC condition, construction of the stencils S j , and QoI calculation.

Since the nozzle code is just a cheap test problem, Tables 3 and4 show that computing the QoI samples v only takes up a significant portion of T for low d. For the baseline SSC method the construction of ENO-type stencils makes up a significant part of the computational cost, but the enforcement of the LEC condition is the most expensive component in higher dimensions. Thus, for the baseline method, most of the computational effort is put into enforcing the LEC condition. For that reason the computational cost of the LEC condition is investigated in more detail.

As explained in Section 2.1, the LEC condition ( 14) is enforced by a MC approach, for all simplices in S j at all j = 1, • • • , n e . Thus, for the baseline SSC method the number of times the surrogate model is evaluated in each iteration i is bounded by

n w i = n e × n e,S j , i = 1, • • • , I (45)
where n e,S j is the number of simplices in a single stencil S j with p = p max , and I is the total number or of iterations of the SSC algorithm. Here we assumed that per S j , one sample is placed in each simplex using [START_REF] Pini | Hybrid adjoint-based robust optimization approach for fluid-dynamics problems[END_REF]. The number of points in the Delaunay grid is given simply by ( 5), but estimating n e for arbitrary d is not trivial. The worst-case number of simplices in a Delaunay triangulation is bounded by the so-called Upper-Bound theorem, which states that n e is at most of O(n d/2 s ).

In the best-case scenario (points distributed uniformly at random inside the unit sphere) n e scales as O(n s ) for any d, with a constant factor that is exponential with the dimension [START_REF] Amenta | Complexity of Delaunay triangulation for points on lower-dimensional polyhedra[END_REF]. To find out where in between these two bounds our specific problem resides, we plot n e versus n s in Fig. 10 for d ∈ {3, 4, 5}. These results indicate that we are close to the O(n s ) bound, since the n e (n s ) are described quite well by the linear regression also shown in Fig. 10. However, the exponential increase of dn e /dn s means that for a moderate number of samples we can still have a large number of simplices if d is high enough.

Note that other than limiting the number of samples n s , we have no means of controlling the magnitude of n e .

The term n e,S j in (45) grows exponentially with p j for a given d. This can be seen in Fig. 11(a), where we plot n e,S j versus the local polynomial order p j for d = 5. Unlike n e , we obviously have some control over the magnitude of n e,S j through the inclusion of a maximum allowable cutoff value for p j . Note however that limiting p j will affect the order of convergence [START_REF] Foo | Multi-element probabilistic collocation method in high dimensions[END_REF]. The upper bound (45), added over iterations i is plotted as a function of n s in Fig. 11

(b) for d = 2, • • • , 5.
It shows a rapid increase with both d and p j . By comparing the computational time T of the SSC method with that of the SSC-SC method (Tables 34), it is clear that the SSC-SC method is several orders of magnitude more efficient for the dimensions considered in this example. To clearly show why the SSC-SC method is computationally more efficient than the baseline method, consider Fig. 12. Here we display the fraction of the volume ¯ that is covered by SCENO stencils of different polynomial order p j , for the discontinuous and smooth case with d = 5. Also the number of stencils |S j | per order is shown. Note that for the discontinuous QoI, just 7 high-order SCENO stencils (stencils with p j > 1) already cover 76.1% of the domain. In the case of the baseline method the number of stencils (and thereby LEC iterations) equals n e , which is 11 034 in this example. For the smooth QoI (Fig. 12(b)) we required 13 fourth-order stencils to cover the entire domain. With the baseline method we would have a set of 9451 stencils, likely all of fourth order as well and therefore large and expensive. As stated in Section 2.1, our primary interest is computing the statistical moments of the QoI, in particular the mean and standard deviation. To assess the accuracy of the SSC method we used a reference solution for each considered dimension d.

To compute the errors we define the following relative L 2 error measures for the mean, standard deviation and interpolation surface

μ = μ w -μ ref 2 μ ref 2 , σ = σ w -σ ref 2 μ ref 2 , w = w(ξ ref ) -v ref 2 v ref 2 . ( 46
)
Here, the subscript w denotes a quantity computed with the surrogate model, and ref is the exact value computed via random sampling. In the interpolation surface error, v ref is a vector containing 10 4 MC code samples and w(ξ ref ) are the surrogate model outputs evaluated at the same MC locations ξ ref .

The values of the error measures (46) for both QoIs and both surrogate models can be found in Tables 56. Note that the error levels are roughly the same for both surrogate models.

From Tables 5-6 we note that the errors of the discontinuous case are considerable higher than for the smooth case. This can be attributed to the smearing of discontinuities, i.e. the linear interpolation of a discontinuity over a simplex, which especially contributes to the error of the surrogate model in higher dimensions. See for instance Fig. 13, which depicts 2D projections of a 3D surrogate model along with reference data on an ordered uniform grid. Especially in Fig. 13(a) we can clearly identify regions where the smearing of the discontinuity contributes to the error. For this particular case, we plotted the difference between the surrogate model and the reference data in Fig. 13(b), which also identifies sharp regions of high error. This situation gets progressively worse as d increases.

In [START_REF] Witteveen | Subcell resolution in simplex stochastic collocation for spatial discontinuities[END_REF][START_REF] Witteveen | Introducing essentially non-oscillatory stencil selection with subcell resolution into uncertainty quantification[END_REF] Witteveen et al. apply a sub-cell resolution approach to the SSC method for the case when the discontinuity in the probabilistic space is a function of a physical discontinuity with random location. Our results indicate that for high d sub-cell resolution could prove to be beneficial, even if the physical location of the QoI is not random.

Computational cost at d = 5

Tables 3-4 show the computational cost at a fixed number of samples per dimension. In this section we investigate the cost of both the baseline and our SC approach as a function of the number of samples n s . Specifically, we consider the smooth nozzle flow case at d = 5, and add 50 code samples at each iteration until an imposed maximum of 500 samples.

The results are depicted in Fig. 14. In light of previous results (Tables 34), we limited the maximum polynomial order to 4 in case of the baseline method in an attempt to suppress the cost. Still, it took 1734 minutes (roughly 29 hours) to create a surrogate model containing 333 code samples. Furthermore, attempts to create surrogates with more samples terminated prematurely due to the excessive memory requirements of the algorithm. With our SC approach we were able to reach the target of 500 samples in 20 minutes, without limiting the polynomial order.

Algebraic test function

To test the limitations of both the baseline and SSC-SC approach in dimensions higher than 5 we make use of the following algebraic test function from [START_REF] Witteveen | Simplex stochastic collocation with random sampling and extrapolation for nonhypercube probability spaces[END_REF] u(ξ ) = arctan ξ • ξ * + ξ * 1 .

(47)

Here, ξ is a d-dimensional vector of uniformly distributed random variables on [0, 1], and ξ * = {0.5, • • • , 0.5} ∈ R d .

Table 7

The computational cost for the arctan test function. A '×' signifies a failed attempt. The computational time is shown in Table 7. In the case of the baseline method we encountered the same behavior as in the preceding section, i.e. an exponential increase in the cost. Moreover, we were not able to create surrogates for d > 5 due to excessive memory requirements.

We were able to create surrogates up to d = 7 without difficulty with the SSC-SC approach. The limitations of the SSC-SC approach start to appear in the 8-dimensional case. Note that for the given number of samples in Table 7 we could still construct the surrogate in a reasonable amount of time, i.e. a run-time less than the expected cost of drawing 350 samples from some expensive simulation code. However, a sharp increase in the cost compared to the 7-dimensional case is still observed. The cause of this are the violations of the LEC limiter which occur despite the fact that (47) is a smooth QoI. As a consequence the SCENO stencils no longer cover the entire domain, and additional stencils must be computed for which the LEC limiter must also be enforced. Several iterations of the LEC and SCENO subroutines are required until the stencils are both LEC and set covering. A possible explanation for the LEC violations is that our monomial basis [START_REF] Loeven | Probabilistic collocation: an efficient non-intrusive approach for arbitrarily distributed parametric uncertainties[END_REF] is no longer suitable in these high-dimensional spaces, despite the fact that we solve our linear system as depicted in [START_REF] Rabitz | Efficient input-output model representations[END_REF]. Further research is required for spaces of 8 dimensions. Another bottleneck for higher dimensions is the ability to create Delaunay triangulations for d > 8, which is not supported by the Delaunay subroutines we used [START_REF] Barber | The qhull software library[END_REF]. Moreover, the time complexity of the Delaunay triangulation in general for these dimensions is a known issue, see for instance [START_REF] Boissonnat | Incremental construction of the Delaunay triangulation and the Delaunay graph in medium dimension[END_REF] for a comparison between different triangulation algorithms up to a dimension of 6. The complexity bounds discussed in Section 5.1.1 and displayed in Fig. 10 will ultimately make the Delaunay triangulation unusable for our purpose, especially when considering problems of dimension O [START_REF] Congedo | Simplex-simplex approach for robust design optimization[END_REF].

Thus, the use of the proposed set-covering approach is that it allows us to efficiently create SSC surrogate models, within the range of dimensions where we could construct the Delaunay triangulation. To apply the SSC method to higher dimensions, one option is to replace the Delaunay triangulation with a scheme where simplices are formed by selecting the nearest points from randomly placed MC samples as described in [START_REF] Witteveen | Simplex stochastic collocation with random sampling and extrapolation for nonhypercube probability spaces[END_REF]. Another option is the use of a dimension adaptive approach such as cut-HDMR in order to avoid high-dimensional spaces altogether, the results of which are shown in the coming sections.

cut-HMDR applied to nozzle flow

We now show the results obtained from using the cut-HDMR approach coupled with the SSC-SC method, applied to the nozzle flow case. Unlike in Section 5.1.1, we now consider all 6 uncertain coefficients from Table 2.

The first order weights α i [START_REF] Witteveen | Simplex stochastic collocation with eno-type stencil selection for robust uncertainty quantification[END_REF] are used to determine the important dimensions. Fig. 15 shows α i with i = 1, • • • , 6 corresponding to the 6 parameters of Table 2, as well as the subsequent higher order α u . For the discontinuous QoI only the first 2 parameters are significant (p and p t ). Together they are responsible for 99.7% of the total first order mean. We have set 1 = 0.9, such that those dimensions which make up 90% or more of the p-th order mean are added to V imp . For the smooth QoI we need p and p t as well, but also the coefficient b in order to meet this constraint. All admissible subsequent dimensions are important as well.

The computational time for both QoIs, due to the fact that the nozzle code is quickly evaluated, is in the order of several minutes. Note that in the case of the baseline method we were unable to create a surrogate in a six-dimensional space.

The values of error measures (46) for the discontinuous and smooth case are given in Tables 89. Notice that a firstorder HDMR expansion is not sufficient for the discontinuous QoI, but already for p = 2 the relative errors in the statistical moments are of O(10 -2 ). Table 9 shows the results for the smooth QoI. The errors, even for a first-order expansion, are of O(10 -3 ) or below. Again, the higher errors in the discontinuous case can be attributed to the linear smearing of discontinuities in the response surface.

Compare the relative errors from Table 5 with those of Table 8, and likewise for Tables 6 and9. The errors are of a similar order of magnitude, even though the HDMR method doesn't sample the full six-dimensional space. At the same time we gain information about the correlation between the input parameters, in the sense of their combined impact on the code output. Also, the cut-HDMR method could be applied to even higher dimensional spaces, provided that the effective dimension is low.

These results demonstrate the power of the cut-HDMR technique in terms of computational efficiency for problems with low effective dimension. The cost of computing 6 one-dimensional and 1-3 two-dimensional surrogate models is significantly less than computing 1 six-dimensional problem. Also, in the case of the SSC method, it can avoid problems with the bad scalability of the Delaunay triangulation with increasing d.

cut-HDMR applied to airfoil flow

In this section we present the results of the cut-HDMR approach (again coupled to the SSC-SC method), when applied to a computationally expensive problem, i.e. the turbulent flow over a NACA0012 airfoil, see Fig. 16(a). The free-stream Mach number is 0.5 and the angle of attack is set to 5 • . The Reynolds number based on the chord length is 1.2 • 10 7 . The grid is a C-grid with 70.085 nodes and the first node is located at a distance of 10 -6 from the wall in order to provide sufficient resolution, see Fig. 16(b). The computation time for one flow-field was roughly 50 minutes.

The governing equations are the Reynolds-Averaged Navier-Stokes (RANS) equations, coupled to the Spalart-Allmaras (SA) turbulence model [START_REF] Spalart | A one-equation turbulence model for aerodynamic flows[END_REF]. This model contains 7 empirically determined closure coefficients [START_REF] Wilcox | Turbulence Modeling for CFD[END_REF], whose best-fit values are unknown a-priori [START_REF] Edeling | Predictive rans simulations via Bayesian model-scenario averaging[END_REF]. Therefore, we treat all 7 inputs as uniformly distributed variables with the end points set at this case there is no guarantee that the LEC limiter is respected in the simplices where the extrapolation takes place, and only in the limit n e → ∞ full extremum-diminishing robustness is recovered for the entire domain .

Finally, for each code run we saved the results for a range of different physical quantities. Thus, we can a posteriori construct a surrogate model for each of these quantities. However, the samples were adaptively placed based on our chosen QoI c l , and therefore might not be optimally distributed for another QoI. Moreover, the dimension adaptivity might cut out dimensions that are important for a QoI other than c l . To perform a qualitative investigation, we therefore constructed a first-order HDMR expansion for the pressure coefficient c p , skin-friction coefficient c f , Mach number M and the turbulent kinetic energy k, all defined below.

c p := p -p ∞ 1 2 ρ ∞ V 2 ∞ , c f := τ w 1 2 ρ ∞ V 2 ∞ , M := u a , k := 1 2 u 2 1 + u 2 2 + u 2 3 ( 48
)
Here, p and p ∞ are the static and free-stream pressure respectively. Furthermore, ρ ∞ and V ∞ are defined as before in c l .

The quantity τ w is the wall shear stress, and in the expression for M, u is the local velocity, and a is the speed of sound.

Finally, u 2 i is the mean squared normal Reynolds stress in x i direction [START_REF] Wilcox | Turbulence Modeling for CFD[END_REF].

In Table 12 we show the coefficients of the SA model sorted according to their value of the (non-relative) weight α i [START_REF] Witteveen | Simplex stochastic collocation with random sampling and extrapolation for nonhypercube probability spaces[END_REF] for all the QoIs of (48). Note that which coefficients are most influential does not change much from one QoI to another. For c p , c f and M the three most important coefficients are κ, C b1 and C b1 . In the case of the turbulent kinetic energy k, σ has taken third place, with again C b1 and κ as the most influential. Still, σ 's value of α i is close to the weight corresponding to C b2 . The ranking of the bottom three coefficients remains unchanged for all considered QoIs.

Conclusion

We have examined means to improve upon the scalability of the Simplex-Stochastic Collocation (SSC) method [START_REF] Witteveen | Simplex stochastic collocation with eno-type stencil selection for robust uncertainty quantification[END_REF] for uncertainty quantification problems of moderate dimensionality, such that sampling an expensive simulation code will still be feasible. We found that creating a surrogate model using the baseline SSC method becomes restrictively expensive at 5 dimensions, and practically impossible at d = 6. To reduce this bad scalability, we needed to add some new features. First, for higher d we run the risk of obtaining a singular sample matrix. This can be circumvented by a method similar to the method used to construct the ENO stencils. If we encounter a stencil S j for which the sample matrix is singular, we collect a set of candidate nearest-neighbor stencils which all contain the simplex j associated to the j-th stencil S j . We then select the stencil which has the highest polynomial order and which is non-singular. In the worst-case scenario we get a linear stencil, which is guaranteed to be non-singular irrespective of the dimension d.

Due to the exponential increase in the number of simplex elements with increasing dimension, enforcing the LEC condition becomes quickly very expensive for d ≥ 5. As a first measure to combat this sharp increase in the computational burden, we have proposed an alternative technique for the stencil construction, based on the Set-Covering (SC) problem [START_REF] Karp | Reducibility Among Combinatorial Problems[END_REF]. Unlike in the SSC method we do not construct a stencil for every simplex. Since n e increases exponentially fast with d, we aim to find a relatively small set of stencils that covers all simplices in the probability domain. Due to the fact that the stencil size rises also exponentially, only a few high-order stencils are required to achieve this. We furthermore use the information contained in the code samples about the location of a possible discontinuity in the construction of the stencils. Using a simple measure based upon the maximum jump in code samples, we mark those simplices which contain a discontinuity and manually set their respective stencils to first order. Using only the remaining simplices as admissible candidates, we grow stencils by adding neighboring simplices. This approach assures that the number of stencils is significantly lower than n e , and hence the number of times the LEC condition must be checked is reduced equally. Also, since the discontinuous simplices are removed as admissible candidates a priori, no stencil crosses a discontinuity. As a consequence the SC stencils resemble the ENO stencils in shape.

For dimensions 5 through 8 our SSC-SC method is significantly more computationally efficient as the baseline SSC method. However, at higher dimensions discretizing the probabilistic space using a Delaunay triangulation becomes a serious bottleneck. The SSC-SC method can therefore be viewed as a more efficient version of the baseline method within the range of dimensions where a Delaunay triangulation can be reasonably made. We examined another alternative, where we adapted the cut-HDMR method of [START_REF] Ma | An adaptive high-dimensional stochastic model representation technique for the solution of stochastic partial differential equations[END_REF] to the SSC method. Given a problem with low effective dimension, this approach circumvents the bad scalability of the Delaunay triangulation, while at the same time obtaining error estimates of similar order of magnitude compared to the full-dimensional baseline or SSC-SC method. Thus, this method is adaptive in both the stochastic domain as well as in the dimensions themselves, while retaining the advantages of the SSC method such as Runge-phenomenon free interpolation. We applied it to a computationally expensive flow case, i.e. the turbulent flow over an airfoil modeled with the Spalart-Allmaras eddy-viscosity model, which contains 7 imperfectly known closure coefficients. Instead of fully sampling a 7 dimensional space, and enforcing the LEC condition in this space, we could obtain a converged surrogate model with a second order HDMR expansion and 62 code samples. In our case A is the hypercube K d and B is the target simplex . Also, we have uniform i.i.d. To simplify the analysis we will consider the standard simplex with nodes ξ 0 = (0, 0, • • • , 0), ξ 1 = (1, 0, • • • , 0), ξ 2 = (0, 1, • • • , 0) etc. To demonstrate the structure of (D.1), we provide a four-dimensional example here:

Y i = h i (R 1 , • • • , R d ), i = 1, • • • , d, f Y 1 •••Y d , is given by: f Y 1 •••Y d = | det J (y 1 , • • • , y d )| • f R 1 •••R d h -1 1 (y 1 , • • • , y d ), • • • , h -1 d (y 1 , • • • ,
R j such that f R 1 •••R d = f R 1 f R 2 • • • f R d = 1.
M 4 = ⎡ ⎢ ⎢ ⎣ 4 √
• • • f Y 1 •••Y d dY 1 • • • dY d = f Y 1 •••Y d • • • dY 1 • • • dY d = d! • 1 d! = 1, (D.8)
which completes the proof.

Fig. 1 .

 1 Fig. 1. Delaunay triangulation for two stencil types with a discontinuity running along the dotted line. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 2 .

 2 Fig.2. The sub simplex (dotted line) of a two-dimensional simplex. Upon refinement one sample is placed at a randomly selected location inside the sub simplex in order to avoid clustering of points.

Fig. 4 .

 4 Fig. 4. Selecting a point inside a triangle and tetrahedron.

Fig. 5 .

 5 Fig. 5. An example of the map (25) for d = 2, 3 and 1000 samples.

Fig. 6 .

 6 Fig. 6. Two stencils which overlap each other. The dark simplices are shared by both stencils.

Fig. 7 .

 7 Fig. 7. Discontinuous simplices identified by (26).

Fig. 8 .

 8 Fig. 8. A two-dimensional example of the SC stencil construction.

Fig. 9 .

 9 Fig. 9. M out as function of p and p t obtained by MC sampling, with the geometrical constants fixed to their nominal value.

Fig. 10 .Fig. 11 .

 1011 Fig. 10. n e as function of n s for the smooth QoI. The discontinuous QoI gives a similar figure. The slope dn e /dn s is computed via a least-square regression line.

Fig. 12 .

 12 Fig. 12. The volume coverage of SCENO stencils per polynomial order.

Fig. 13 .

 13 Fig. 13. 3D surrogate model displayed in 2 dimensions by fixing 1 dimension to a particular value. The green dots are reference data on an ordered uniform grid. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 14 .

 14 Fig. 14. The computational time T in minutes versus n s at d = 5.

Fig. 15 .

 15 Fig.[START_REF] Guillas | Bayesian calibration of the constants of the k-ε turbulence model for a cfd model of street canyon flow[END_REF]. The relative weights[START_REF] Witteveen | Simplex stochastic collocation with eno-type stencil selection for robust uncertainty quantification[END_REF] and (42) for the discontinuous and smooth nozzle QoI. The green boxes are the dimensions which are added to V imp . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 16 .

 16 Fig. 16. The symmetrical NACA0012 airfoil.

Fig. 17

 17 Fig.17. The relative weights[START_REF] Witteveen | Simplex stochastic collocation with eno-type stencil selection for robust uncertainty quantification[END_REF] and (42) for the NACA 0012 test case. The value of 1 was set to 0.9. The green boxes indicate the dimensions which are important. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

.

  Fig.17. The relative weights[START_REF] Witteveen | Simplex stochastic collocation with eno-type stencil selection for robust uncertainty quantification[END_REF] and (42) for the NACA 0012 test case. The value of 1 was set to 0.9. The green boxes indicate the dimensions which are important. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Theorem 3 . 2 )

 32 Consider the r.v.'s R 1 , • • • , R d with joint pdf f R 1 •••R d positive and continuous on the set A ⊆ R d , and let h 1 , • • • , h d be real-valued transformations defined on A; that is, h 1 , • • • , h d → R, and let B be the image of A under transformations(h 1 , • • • , h d ). Suppose that (h 1 , • • • , h d ) is one-to-one from A onto B. Thus, if we set y i = h i (r 1 , • • • , r d ), we can uniquely solve for r i , i = 1, • • • , d : r i = h -1 i (y 1 , • • • , y d ), i = 1, • • • , d. Suppose further that the partial derivatives ∂ ∂ y j h -1 i , i, j = 1, • • • , dexist and are continuous for (y 1 , • • • , y d ) ∈ B. Finally, suppose that the Jacobian J (y 1 , • • • , y d ) is = 0 on B. Then the joint pdf of the r.v.'s

  y d ) , (y 1 , • • • , y d ) ∈ B 0 otherwise (D.3)

  Thus, in order for f Y 1 •••Y d to be uniform we need to show that | J (y 1 , • • • , y d )| is a constant. Furthermore, since (D.1) maps to a simplex this constant must be equal to the reciprocal of the volume of the simplex.

(D. 5 ) 4 (y 4 +y 3 ) 2 y 3 (y 4 +y 3 ) 2 0- 2 (y 4 +y 3 ) 2 (y 4 +y 3 +y 2 ) 3 2 (y 4 2 (y 4 3 - 3 (y 4 +y 3 +y 2 ) 3 (y 1 3 (y 4 3 (y 4 3 (y 4 6 )|

 5423222242433313434346 Now we can compute the Jacobian matrix (D.2) as +y 3 )y 2 (y 4 +y 3 +y 2 ) 3 +y 3 )y 2 (y 4 +y 3 +y 2 ) +y 4 +y 3 +y 2 ) 4 +y 3 +y 2 ) 2 y 1 (y 1 +y 4 +y 3 +y 2 ) 4 +y 3 +y 2 ) 2 y 1 (y 1 +y 4 +y 3 +y 2 ) 4 +y 3 +y 2 ) 2 y 1 (y 1 +y 4 +y 3 +y 2 )4 4 (y1 + y 4 + y 3 + y 2 ) 3 4 (y 1 + y 4 + y 3 + y 2 ) 3 4 (y 1 + y 4 + y 3 + y 2 ) 3 4 (y 1 + y 4 + y 3 + y 2 )When we compute the determinant of J all the y i terms drop out and we end up with | det J | = 24. The values of | det J | for d = 1, • • • , 8 can be found in TableD.[START_REF] Foo | The multi-element probabilistic collocation method (me-pcm): error analysis and applications[END_REF]. From these results it becomes clear thatf Y 1 •••Y d = | det J | = d!. The volume ¯ det (D) |, D = ξ 1ξ 0 ξ 2ξ 0 • • • ξ d+1ξ 0 , (D.7)which for the standard simplex reduces to ¯ = 1 d! . And thus we have

Table 1

 1 Examples of ill-conditioned systems. We show the dimension d, the polynomial order of the stencil, the nullity and condition number of the sample matrix , and finally the condition number of the non-null range .

	d	p j	Dimension	Nullity	Cond.	Cond. range
	6	2	2 8 × 28	1	1.36e+17	9.39e+3
	6	3	8 2 × 82	1	1.31e+17	2.40e+4
	6	3	8 4 × 84	2	2.85e+17	2.66e+3

Table 2

 2 Uncertain input parameters of the discontinuous (D) and smooth (S) case.

	d	Parameter	Mean (D)	Range (D)	Mean (S)	Range (S)
	1	p [bar]	0.55	[0.5, 0.6]	0.625	[0.60, 0.65]
	2	p t [bar]	1.0	[0.9, 1.1]	1.0	[0.9, 1.1]
	3	a [-]	1.75	[1.575, 1.925]	1.75	[1.575, 1.925]
	4	b [-]	0.7	[0.63, 0.77]	0.7	[0.63, 0.77]
	5	c [-]	0.8	[0.72, 0.88]	0.8	[0.72, 0.88]
	6	d [-]	4.0	[3.6, 4.4]	4.0	[3.6, 4.4]

Table 3

 3 The computational cost of the discontinuous QoI.

	Type	d [-]	n s [-]	T [min]	LEC [%T]	S j [%T]	v [%T]
	Baseline	2	50	0.56	3.56	3.16	87.3
		3	1 0 0	2 .09	24.39	11.46	39.32
		4	150	10.95	73.42	15.37	6.22
		5	200	119.29	85.21	11.26	0.58
	SCC-SC	2	50	0.54	1.45	1.24	82.46
		3	1 0 0	1 .33	1.2	2 .33	54.75
		4	150	1.37	5.56	12.34	42.99
		5	2 0 0	4 .75	4.88	17.2	1 1 .47

Table 4

 4 The computational cost of the smooth QoI.

	Type	d [-]	n s [-]	T [min]	LEC [%T]	S j [%T]	v [%T]
	Baseline	2	50	0.73	2.28	2.87	89.9
		3	1 0 0	2 .52	20.37	16.42	42.07
		4	150	22.86	62.18	30.87	3.95
		5	200	731.5	5 8 .31	40.99	0.13
	SCC-SC	2	50	0.7	1 .28	0.31	85.64
		3	1 0 0	1 .65	4.0	0 .43	61.14
		4	150	1.63	16.76	1.26	49.01
		5	2 0 0	4 .68	13.62	1.41	15.45

Table 5

 5 The relative errors (46) of the discontinuous QoI for the baseline SSC and SSC-SC method.

	Type	d	n s	μ	σ	w
	Baseline	2	50	3.536e-02	3.669e-02	2.001e-01
		3	100	5.271e-02	7.475e-02	2.408e-01
		4	150	2.532e-02	1.425e-01	2.828e-01
		5	200	2.006e-02	2.320e-01	3.253e-01
	SCC-SC	2	50	1.590e-02	4.397e-02	1.597e-01
		3	100	3.975e-03	7.452e-02	2.108e-01
		4	150	1.199e-03	1.329e-01	2.547e-01
		5	200	3.368e-03	1.803e-01	2.876e-01

Table 6

 6 The relative errors (46) of the continuous QoI for the baseline SSC and SSC-SC method.

	Type	d	n s	μ	σ	w
	Baseline	2	50	1.131e-06	5.137e-06	1.088e-05
		3	100	7.276e-07	1.572e-05	1.416e-05
		4	150	1.015e-06	3.784e-06	2.566e-05
		5	200	2.446e-05	3.376e-04	2.228e-03
	SCC-SC	2	50	8.042e-07	1.952e-05	1.916e-05
		3	100	8.734e-07	3.019e-07	7.075e-06
		4	150	1.006e-06	2.234e-06	2.104e-05
		5	200	3.034e-06	8.186e-06	8.218e-05

Table 8

 8 Relative error values vs the HDMR order p and n s in the case of the discontinuous QoI.

	HDMR order p	1	n s	μ	σ	w
	1	0.9	57	4.514e-01	3.601e-01	7.434e-01
	2		89	3.926e-02	8.830e-02	2.196e-01

Table 9

 9 Relative error values vs the HDMR order p and n s in the case of the smooth QoI.

	HDMR order p	1	n s	μ	σ	w
	1	0.9	37	4.933e-05	2.885e-04	5.035e-03
	2		6 7	3 . 6 0 3 e -05	6.457e-04	3.608e-03

Table 10

 10 Uniformly distributed closure coefficients of the SA model.

	Parameter	Mean	Range
	C b1	0.1355	[0.0949, 0.1762]
	C b2	0.622	[0.435, 0.809]
	C v1	7.1	[4.97, 9.10]
	σ	2/3	[0.467, 0.867]
	C w2	0.3	[0.210, 0.390]
	C w3	2	[1.40, 2.60]
	κ	0.41	[0.287, 0.455]

Table 11

 11 The convergence of the statistical moments of c l as a function of the HDMR order.

	HDMR order p	n s	μ	σ
	1	49	6.385509e-01	3.729915e-03
	2	62	6.381699e-01	3.746544e-03

Table 12

 12 The rankings from largest (1st) to smallest (7th) α i for QoIs other than c l . Below each coefficient its corresponding value of α i is written.

	QoI	1st	2nd	3rd	4th	5th	6th	7th
	c p	C b1	κ	C b2	σ	C v1	C w2	C w3
		2.604e-03	1.339e-03	8.784e-04	5.510e-04	6.866e-05	3.674e-08	3.655e-08
	c f	κ	C b1	C b2	σ	C v1	C w2	C w3
		1.764e-01	8.835e-02	8.676e-02	2.072e-02	6.819e-03	9.976e-07	9.875e-07
	M	κ	C b1	C b2	σ	C v1	C w2	C w3
		3.658e-03	1.951e-03	1.498e-03	6.591e-04	1.233e-04	3.958e-08	3.951e-08
	k	C b1	κ	σ	C b2	C v1	C w2	C w3
		8.708e-02	2.982e-02	2.726e-02	2.529e-02	5.735e-04	2.650e-05	2.591e-05

Table D .

 D [START_REF] Foo | The multi-element probabilistic collocation method (me-pcm): error analysis and applications[END_REF] The absolute value of the determinant of J as a function of d.

	d	1	2	3	4	5	6	7	8
	| det J |	1	2	6	2 4	1 2 0	7 2 0	5 0 4 0	4 0 3 2 0

  r 4 -4 √ Each individual row of M 4 is a transformation function y i = h i (r 1 , • • • , r 4 ),for which we can find the following inverse + y 3 , r 2 = (y 4 + y 3 ) 2 (y 4 + y 3 + y 2 ) 2 , r 3 = (y 4 + y 3 + y 2 ) 3 (y 1 + y 4 + y 3 + y 2 ) 3 , r 4 = (y 1 + y 4 + y 3 + y 2 ) 4 .

			r 4 3 √ r 3 √ r 3 -4 √ r 4 3 √ 4 r 4 3 √ r 3 r 4 3 √ 4 √ r 3 √ r 2 -4 √ r 4 3 √ r 3 √ r 2 √ r 2 r 1 4 √ r 4 3 √ r 3 r 2 r 1 √	⎤ ⎥ ⎥ ⎦ .	(D.4)
	functions:
	r 1 =	y 4	y 4

Acknowledgements

The present work was supported by the French Agence Nationale de la Recherche (ANR) under contract ANR-11-MONU-008-002.

Appendix A. Baseline SSC algorithm

This appendix provides the general pseudo code of the baseline SSC algorithm with ENO stencils.

Compute initial d-dimensional Delaunay grid

Compute initial compute code samples v = (v 0 , • • • , v k , • • • , v d ) at the 2 d + 1 grid points Set initial hierarchical surplus errors to k = -v k i ← 0, choose max iterations I while i ≤ I and ˆ rms > user-defined threshold do Determine p max via (5) Compute nearest-neighbor stencils via [START_REF] Chung | On lattices admitting unique Lagrange interpolations[END_REF] % Check (in parallel) the LEC condition:

† Sample all j ∈ S j via [START_REF] Pini | Hybrid adjoint-based robust optimization approach for fluid-dynamics problems[END_REF] At sample locations: compute w j if [START_REF] Ganapathysubramanian | Sparse grid collocation schemes for stochastic natural convection problems[END_REF] is violated in one or more sample locations then

Update S j and goto † end if end for % Compute (in parallel) ENO stencils

Collect the r nearest neighbor S j that contain j From these r stencils, select S * j := the S j with the largest p j if more than 1 S * j exists then Select the S * j with the smallest average Euclidean distance to the cell center of j end if end for Compute probabilities j and volumes ¯ j via (10) Compute refinement measures ē j via (9) Sample the sub-simplices (13) of the N simplices with the largest ē j via (25) Evaluate (in parallel) computer code at the N new sample locations, add code samples to v Compute the N new hierarchical surplus values via [START_REF] Guillas | Bayesian calibration of the constants of the k-ε turbulence model for a cfd model of street canyon flow[END_REF] Refine the Delaunay grid by adding the N new sample locations end while Compute ENO stencils and check LEC condition from final iteration

Appendix B. SSC-SC algorithm

Below we give the detailed construction of SCENO stencils in pseudo code, and we repeat some definitions for convenience.

• C is the set of simplices which are currently covered by a stencil S j .

• C j is the set of simplices currently under construction which will be added to C. • N j is the set of all neighboring simplices k of simplex j .

• N C j is the set of all neighboring simplices k of all simplices in C j .

• {ξ k } is the set of d + 1 vertices ξ k j,l that make up simplex k .

% for each iteration of the SSC-SC algorithm do

sort C j and S j according to ξ center kξ center j 2 if |S j | < N j + 1 then reduce p j and S j to new maximum as allowed by ( 5)

Appendix C. HDMR algorithm

The dimension-adaptive cut-HDMR of [START_REF] Ma | An adaptive high-dimensional stochastic model representation technique for the solution of stochastic partial differential equations[END_REF] coupled with the SSC method. %Initialize sets

%Compute the component functions of order p = 0 and p = 1 Compute all sub problems f (v) in (36) using the SSC method. Stop when global RMS error measure (17) < . Add all computed u to V comp . Compute all first-order weights [START_REF] Witteveen | Simplex stochastic collocation with random sampling and extrapolation for nonhypercube probability spaces[END_REF]

Add u with |u| = p and which satisfy admissibility condition to R while R = {∅} and α p > 2 do p ← p + 1 Add u with |u| = p and which satisfy admissibility condition [START_REF] Witteveen | An adaptive stochastic finite elements approach based on Newton-Cotes quadrature in simplex elements[END_REF] to R ∀u ∈ R, compute f (v) in (36) using SSC, stop when (17) < , add all computed u to V comp .

Compute weights α u [START_REF] Witteveen | Subcell resolution in simplex stochastic collocation for spatial discontinuities[END_REF] From [START_REF] Roussas | An Introduction to Probability and Statistical Inference[END_REF] we have the following theorem regarding the distribution of a transformation of random variables: