Clément Scherding
email: clement.scherding@dalembert.upmc.fr

Georgios Rigas

Denis Sipp

Peter J Schmid

Taraneh Sayadi

Data-driven framework for input/output lookup tables reduction: Application to hypersonic flows in chemical nonequilibrium

I. INTRODUCTION

Chemical non-equilibrium effects have been shown to play an important role in the accurate simulation of flows at hypersonic conditions and in the computation of design characteristics, such as transition location or thermal loading [START_REF] Marxen | A method for the direct numerical simulation of hypersonic boundary-layer instability with finite-rate chemistry[END_REF][START_REF] Marxen | Nonlinear instability of a supersonic boundary layer with two-dimensional roughness[END_REF][START_REF] Candler | Rate Effects in Hypersonic Flows[END_REF][START_REF] Di Renzo | Direct numerical simulation of a hypersonic transitional boundary layer at suborbital enthalpies[END_REF]. Recent studies have identified these effects as causes of order-one changes in growth rates, response behavior or sensitivities, even though the corresponding variations in first-order flow statistics have been modest. These findings have in turn prompted a significant endeavor of augmenting existing flow solvers with non-equilibrium modules to account for finite-rate aerothermochemical features. Simulations in this parameter regime introduce and track a range of species in their inert or ionized forms [START_REF] Zhong | Direct numerical simulation on the receptivity, instability, and transition of hypersonic boundary layers[END_REF][START_REF] Panesi | Nonequilibrium shock-heated nitrogen flows using a rovibrational state-to-state method[END_REF][START_REF] Paredes | Nonmodal growth of traveling waves on blunt cones at hypersonic speeds[END_REF]. Complementing the hydrodynamic state vector by chemical components is a well-established technique, for example in combustion or atmospheric simulations, but the required modeling of the inter-species interactions, such as dissociation, reaction and recombination [START_REF] Anderson | Hypersonic and High-Temperature Gas Dynamics[END_REF], for hypersonic applications poses great challenges. Much of this modeling is accomplished by lookup libraries, which act as repositories of tabulated chemical reactions encountered for a given flow state [START_REF] Scoggins | Mutation++: Multicomponent thermodynamic and transport properties for ionized gases in c++[END_REF]. When passing state-vector components to the library, amplitudes and time-scales for various forcing terms are returned, appearing as exogeneous inputs to the momentum, energy and species transport equations.

Much effort has gone into these libraries such as Pegase [START_REF] Bottin | Thermodynamic and transport properties for inductive plasma modeling[END_REF], Eglib [START_REF] Ern | Eglib: A general-purpose fortran library for multicomponent transport property evaluation[END_REF], Plato [START_REF] Munafò | A computational model for nanosecond pulse laser-plasma interactions[END_REF] and the leading library for reacting flows simulations CHEMKIN [START_REF] Kee | Chemkin collection, release 3.6, reaction design, Inc[END_REF]. For aerothermochemical non-equilibrium effects in hypersonic flows, the Mutation++ library (MUlticomponent Thermodynamic And Transport properties for IONized gases in C++), developed and maintained at the von Karman Institute (VKI), has become the standard for high-fidelity simulations of high-speed and high-enthalpy flows [START_REF] Scoggins | Mutation++: Multicomponent thermodynamic and transport properties for ionized gases in c++[END_REF]. This library can be coupled to existing flow solvers and is capable of modeling a range of partially ionized gas effects, together with non-equilibrium features, energy exchange processes and gas-surface interactions. The flexibility and scope of the library comes at the expense of a computational bottleneck that slows down a typical large-scale simulation by a large factor, as shown in Fig. 1, where typical simulation times for calorically and thermally perfect gases are juxtaposed with results for non-equilibrium chemical reactions. A wide margin can be observed. For this reason, non-equilibrium computations range among the most inefficient and laborious calculations in fundamental hypersonic research. To increase performance, most CFD codes use hard-coded chemistry [START_REF] Di Renzo | Htr solver: An open-source exascale-oriented task-based multi-gpu high-order code for hypersonic aerothermodynamics[END_REF]. However, any change in the gas mixture or the thermochemical model comes at a human cost in terms where hundreds of species are tracked, even with reduced chemical mechanisms [START_REF] Attili | Formation, growth, and transport of soot in a three-dimensional turbulent non-premixed jet flame[END_REF][START_REF] Bansal | Direct numerical simulations of autoignition in stratified dimethyl-ether (dme)/air turbulent mixtures[END_REF][START_REF] Bhagatwala | Direct numerical simulations of hcci/saci with ethanol[END_REF]. Moreover, one common assumption in surrogate modeling relates to the smoothness of the approximated relation. This is not always true, especially in hypersonic applications where shocks and temperature discontinuities are amongst the typical features of such flows. Nonetheless, clever pre-processing steps can greatly improve the model's performance in these cases. For example, Bouhlel et al. [START_REF] Bouhlel | Improving kriging surrogates of high-dimensional design models by partial least squares dimension reduction[END_REF] coupled kriging with partial least-squares (PLS) methods to reduce the high-dimensional (D = 100) input space. In [START_REF] Hawchar | Principal component analysis and polynomial chaos expansion for timevariant reliability problems[END_REF], principal component analysis (PCA) has been used as a pre-processing step before applying polynomial chaos expansions on the PCA basis [START_REF] Hawchar | Principal component analysis and polynomial chaos expansion for timevariant reliability problems[END_REF]. When dealing with discontinuous functions, Bettebghor et al. [START_REF] Bettebghor | Surrogate modeling approximation using a mixture of experts based on em joint estimation[END_REF] proposed to cluster the input basis into different regions (to avoid a discontinuity within a cluster) and build a surrogate model on each of these regions. All models are then combined together and form a mixture of experts (MoE), as described in the literature [START_REF] Hastie | The elements of statistical learning: data mining, inference, and prediction[END_REF]. Yang [START_REF] Yang | Regression with multiple candidate models: selecting or mixing?[END_REF], however, pointed out that combining surrogate surfaces does not necessarily outperform a single model fitted over the entire input space. Hence, special care has to be taken in combining these steps.

The objective of this work is therefore to develop an effective pre-processing technique, allowing the construction of a low-dimensional surrogate model capable of replacing the computationally expensive library and the memory-intensive look-up tables when modeling inter-species interactions in simulations of chemically reactive flows.

The paper is organized as follows. In Section II, the generic algorithm is presented in detail. It combines techniques from nonlinear model reduction, network clustering and surrogate modeling to efficiently extract a surrogate, lightweight model of the full library. In Section III, the governing equations for hypersonic flows in chemical non-equilibrium as well as the thermo-chemical model for such flows are recalled. In Section IV, the algorithm is first tested on the simulation of an adiabatic Mach-10 boundary layer in chemical non-equilibrium, initially studied by Marxen et al [START_REF] Marxen | A method for the direct numerical simulation of hypersonic boundary-layer instability with finite-rate chemistry[END_REF]. The gas properties of the reactive air mixture are initially modeled using the open-source Mutation++ library [START_REF] Scoggins | Mutation++: Multicomponent thermodynamic and transport properties for ionized gases in c++[END_REF]. Replacing the library by the surrogate model then overcomes the computational bottleneck alluded to above. In addition, a shock-wave boundary layer interaction case is also considered in order to highlight the capability of the algorithm to deal with flows presenting discontinuities.

While focusing on the non-equilibrium gas-dynamic library Mutation++, we stress that the employed techniques are agnostic about the particular library they are applied to, and can just as readily be employed to other libraries or lookup tables attached to simulations. Applications in combustion, phase-change simulations or particle-laden flows stand to benefit from this accelerating methodology at the interface between flow solvers and material-property libraries.

II. DESCRIPTION OF THE ALGORITHM

In this section, the general algorithm to extract a reduced library for the thermochemical properties of multicomponent mixtures in chemical non-equilibrium is described.

Compressible flow simulations in chemical non-equilibrium require transport, thermodynamic and chemical reactions' properties, z ∈ R D Z , (e.g. viscosity, conductivity, enthalpies, and chemical source terms) to close the governing equations. These properties are modeled as functions of the local thermodynamic state and mixture composition, concatenated into the local thermodynamic vector qth ∈ R D , usually computed using tabulation or external libraries, and can be considered as an input/output problem:

z = f (q th ), (1)
where the function f represents the library of interest, and D and D Z represent the dimensions of the input and output spaces, respectively. This function then needs to be evaluated at each grid point and at each time-step. While accurate, the extensive calls to the library come at a substantial performance loss for the solver, together with a significant time penalty (see Fig. 1).

While these function calls cannot be entirely avoided, existing features of the flow inspire strategies to seek a less expensive method to evaluate the required properties. (i) Flows have history. In other words, several calls to the library may be redundant since some thermodynamic states are seen multiple times throughout the simulation. (ii) Any flow of interest contains only a subset of all possible thermodynamic states, given its nature and freestream conditions. Hence, only a small subset of the input space of function f needs to be accessed. (iii) While data-driven method requires a lot of data for training, some (rare) flows of interest, such as hypersonic boundary layers in chemical non-equilibrium, exhibit elegant, locally self-similar solutions [START_REF] Lees | Laminar heat transfer over blunt-nosed bodies at hypersonic flight speeds[END_REF] that can be used for training instead of an expensive direct numerical simulation (DNS). This final point will be explored in more detail in future work.

The proposed algorithm leverages these features by creating a surrogate model of the function f only on a subset of input states relevant to the simulation, which is commonly represented as a low-dimensional manifold in R D . This

CFD solver

Governing equations

Discretization

Boundary conditions

Input/Output library

!" # ! $ # ! ! % & ' ( !)%* Input " Output #

Dimensionality reduction

Clustering/Classification Surrogate model allows us to first perform dimensionality reduction of the input space (see [START_REF] Bouhlel | Improving kriging surrogates of high-dimensional design models by partial least squares dimension reduction[END_REF][START_REF] Hawchar | Principal component analysis and polynomial chaos expansion for timevariant reliability problems[END_REF]). Next, following a similar approach as in [START_REF] Bettebghor | Surrogate modeling approximation using a mixture of experts based on em joint estimation[END_REF], regions with different dynamics and/or discontinuities between them are clustered into a low-dimensional representation. Finally, surrogate models are constructed on each cluster in this low-dimensional space. Hence, the training of the algorithm is performed in three steps: (i) dimensionality reduction, (ii) community clustering, and (iii) surrogate model construction. Once trained, the model replaces the look-up library already in place to predict the thermochemical properties of the mixture within the flow solver. We stress that the lighter version of the library will perform correctly only on the range of conditions seen during the simulation. A general schematic of the training process and the coupling with the flow solver is presented in Fig. 2.

Input $ Output % & ! " ! " # # $ # $ " % & ' ( ) ! " ! *%+ , -# $ $ % -# # . / 0 % -1 % 2 )
It should be noted that this strategy is also applicable to flows in thermal non-equilibrium, where the internal energy modes are out of equilibrium with the translational energy of the flow [START_REF] Munafò | Modeling of dissociation and energy transfer in shock-heated nitrogen flows[END_REF]. Since the additional source terms modeling the energy exchange in the internal energy equation(s) are modeled as a function of the local thermodynamic state vector as well, they can be added to the outputs and a surrogate can be constructed accordingly. This endeavor however lies outside of the scope of the current study and will be pursued as a subject of a future work. In the future, on-the-fly adaptivity will also be added to the algorithm to allow the model to learn new states as the simulation is advanced in time. This capability will help tackle more challenging and unsteady flow problems, while alleviating the need for a complete training set which may not always be available. Hence, the final intended use of the model will be as follows: (i) a preliminary laminar simulation will warm-start the training of the algorithm. (ii) as more flow features are added (such as instabilities), the model will learn adaptively the remaining "new" thermodynamics states pertaining to the new unsteady features. Following this procedure, the dynamics will be obtained at a lower CPU cost.

A. Training

The algorithm is trained using the simulation of an adiabatic Mach-10 boundary layer in chemical non-equilibrium, thoroughly described in Section III D. The gas considered is a five-species air model S = N 2 , O 2 , N O, N, O with five reactions. The computational code solves the compressible reactive Navier-Stokes equations, where the kinetic parameters are computed by coupling with Mutation++ [START_REF] Scoggins | Mutation++: Multicomponent thermodynamic and transport properties for ionized gases in c++[END_REF]. The input thermodynamic state vector qth is composed of density ρ, internal energy ρe and the mixture partial densities ρ s , s ∈ S. The outputs of the library fall within three categories: thermodynamic properties (pressure p, temperature T , species specific enthalpies h s , s ∈ S), transport properties (viscosity µ, thermal conductivity κ, diffusion coefficients D s , s ∈ S) and chemical kinetics source term ωs , s ∈ S. More details concerning the governing equations and the thermo-chemical model can be found in Section III.

Data collection

To train the model, N thermodynamic state vectors qth are randomly sampled on the grid of a previously converged simulation in chemical non-equilibrium and concatenated into the input vector X ∈ R N ×D . The corresponding outputs from the library are collected and concatenated into the output vector Z ∈ R N ×D Z . Fig. 3 shows the numerical range of selected output variables along each input, normalized between 0 and 1 with a minimum-maximum scaling

X = X -Xmin Xmax -Xmin , Z = Z -Zmin Zmax -Zmin . ( 2
)
Taking dynamic viscosity µ, for example, it shows a strong dependency on density ρ and internal energy ρe but a low variation with respect to the radicals' partial densities ρ O , ρ N and ρ N O . The same observations can be made for all other outputs. Hence, the variation of the function f with respect to the inputs can be accurately represented on a low-dimensional subspace of the inputs. This motivates the first step of the algorithm, namely, dimensionality reduction.

Dimensionality reduction

The goal of this section is to find an effective algorithm for dimensionality reduction of the input space in order to construct a mapping between its reduced-order representation and the output of the library

Ẑ = g(Y), (3)
where g is the approximation of the scaled library f in the low-dimensional subspace of the inputs, Y is the reducedorder representation of an input X and Ẑ the prediction of the model. The benefit of this first pre-processing step is to maintain high accuracy of the surrogate model, while decreasing the overall cost of construction and evaluation. In fact, constructing a response surface faces the well-known curse of dimensionality; as the number of input dimensions increases, the cost of constructing an accurate surface increases exponentially. This approach was proven successful in [START_REF] Bouhlel | Improving kriging surrogates of high-dimensional design models by partial least squares dimension reduction[END_REF] where PLS was used in tandem with kriging to reduce the dimension of the input space.

a. Principal component analysis The most common technique for dimensionality reduction of a dataset X ∈ R N ×D in high dimensions is principal component analysis PCA (see eg. [START_REF] Shlens | A tutorial on principal component analysis[END_REF]). The principal components of X are found through the eigenvalue decomposition of the covariance matrix of the data. The dataset X is then projected onto d < D leading eigenvectors (or principal components) of the covariance matrix, resulting in a low-dimensional representation Y ∈ R N ×d of the original dataset. However, depending on the shape of the manifold, the variations of the output variables with respect to the low-dimensional sub-space may not be properly preserved, which is the case presented in Fig. 4a with points of high and low temperature projected onto similar locations. This example illustrates the limitations of PCA for dimensionality reduction of a dataset constrained to a nonlinear manifold.

b. Auto-encoders Nonlinear dimensionality reduction via auto-encoders (AE) typically have a higher compression rate than linear techniques. An auto-encoder is a parametric model (i.e. a deep neural network with an activation function σ) that embeds the input dataset X ∈ R N ×D into a low-dimensional representation Y ∈ R N ×d through an encoder function E. The low-dimensional representation is then decoded back to the input space with the decoder function D, producing a reconstruction of the input

X ∈ R N ×D . Y = E(X) X = D(Y) (4)
The weights of the two networks E and D can be trained using back-propagation of the L 2 error ∥X-X∥ 2 through the network. If the activation function is selected as the identity (i.e σ(x) = x), the auto-encoder is linear and unbiased

E = W E D = W D (5)
where W E and W D are the weights matrices of the encoder and decoder, respectively. These optimal weights can be found through PCA. In fact, the linear latent space of dimension d of the encoder will span the same sub-space as the top d PCA singular vectors. The equivalence between the two techniques was first shown by Baldi & Hornik [START_REF] Baldi | Neural networks and principal component analysis: Learning from examples without local minima[END_REF]. Correspondingly, a two-layered nonlinear auto-encoder can be mathematically described as follows

Y = W E,2 σ(W E,1 X + b E,1 ) + b E,2 X = W D,2 σ(W D,1 Y + b D,1 ) + b D,2 (6)
where σ is a nonlinear activation function, W E,1 ∈ R H×D , W E,2 ∈ R d×H are the weights matrices of the first and second layer of the encoder with respective biases b E,1 ∈ R H and b E,2 ∈ R d . H denotes the dimension of the hidden layer. The matrices and bias vectors of the decoder have transposed dimensions. This corresponds to the minimal architecture (i.e. with one hidden nonlinear layer and an output linear layer) requested by the universal approximation theorem [START_REF] Cybenko | Approximation by superpositions of a sigmoidal function[END_REF]. However, more hidden layers can be considered. Fig. 4b shows the manifold unrolled in two dimensions with an auto-encoder, colored by the magnitude of the temperature, an output of the library. The AE outperforms PCA by preserving the local structure and preventing points at different thermodynamic states (i.e different temperatures) to be projected onto the same location. However, the highest temperature zone is concentrated in a thin layer adjacent to the lower temperature area. This will result in strong and unphysical gradients of the surrogate model within this region. c. Partial least-squares Since our interest lies in reducing the dimensionality of the input to construct a reducedorder surrogate model of the input/output relations, it is useful to entangle the input into a low-dimensional space that best reconstructs the outputs. In analogy to PCA finding dependencies between the inputs, partial least-squares (PLS) finds a basis of the input space that optimally accounts for features in the output space. It has been used to construct surrogate models aimed at reducing the dimensions of the input space (see [START_REF] Bouhlel | Improving kriging surrogates of high-dimensional design models by partial least squares dimension reduction[END_REF]). Different variants of PLS now exist, using either a singular value decomposition (PLS-SVD) or iterative algorithms (such as PLS-W2A in [START_REF] Wegelin | A survey of Partial Least Squares (PLS) methods, with emphasis on the two-block case[END_REF]).

While it has been shown that in cases where the dimension of the latent space is strictly greater than one, PLS-SVD differs from PLS-W2A and its variant PLS2, no major differences in the resulting latent space were observed. The results of PLS-SVD are presented here to highlight the similarity to PCA. Given the input X and output vectors Z, the PLS-SVD algorithm [START_REF] Wegelin | A survey of Partial Least Squares (PLS) methods, with emphasis on the two-block case[END_REF] determines

X T Z = UΣV T . (7)
Similar to PCA, the projection of the input is then obtained by projecting onto the d < D top left singular vectors

Y = XU (8)
where U ∈ R D×d is the truncated left-singular matrix. Fig. 4c shows the training set projected onto the twodimensional basis generated with PLS. As expected, adding the output information in the computation improves the output's representation in the reduced-order input basis compared to PCA. However, an artificial discontinuity is created near the high-gradient region that was not present in Fig. 4b. This highlights the lower compression rate of linear techniques compared to nonlinear ones.

d. Input/output-encoders The strategy adopted here is therefore a nonlinear dimensionality reduction of the input data using input/ouput-encoders (IO-E), a modified version of auto-encoders (AE) adapted to input/ouput relations. To that end, the decoder architecture is modified to best reconstruct the output of the library Z ∈ R N ×D Z as

Z = W D,2 σ(W D,1 Y + b D,1 ) + b D,2 (9)
where the weight matrices are now

W D,1 ∈ R H D ×d , W D,2 ∈ R D Z ×H D with respective biases b D,1 ∈ R H D and b D,2 ∈ R D Z .
H D now represents the hidden dimension of the decoder. The resulting network is then trained via back-propagation of the reconstruction error based on the output ∥Z -Z∥ 2 . The input/output-encoder (IO-E) architecture is illustrated in Fig. 5.

Input layer ! " ! Output layer ! " ! ! Hidden layer ! " " Hidden layer ! " " " #$% & # ' & ! ( ) * + , -$ . ! " # Latent space ! " % ! ! ! ! ! ! ! ! ! ! FIG.
5. Schematic of the input/output-encoder architecture proposed. An input sample X is passed through the encoder network to generate a low-dimensional representation Y. Y is decoded back to predict the library output Z. At each training step, the L2 loss ∥Z -Z∥2 between the library and the network prediction is calculcated and backpropagated until convergence of the encoder and decoder weights, WE and WD, respectively, is reached.

Fig. 4d shows the equivalent manifold using the IO-E architecture, with the same number of hidden dimensions. The IO-E architecture outperforms all other techniques as the high-temperature region is now projected onto its own properly defined zone with smooth variations. This feature will have a marked impact on the performance of the library when linked to the solver. It is worth noting that the decoder part of the IO-E could be used directly to predict the output of the library. However, we will see in Section IV B that the accuracy of the full IO-E would not optimal in this case. This motivates the use of clustering and radial basis functions, as described next.

Community clustering

In the second step of the algorithm, we seek to discover clusters within our data. In our present context, a cluster represents a subset of data that shares similar thermodynamic features. These feature classification will then be useful in constructing a dedicated surrogate surface of the low-dimensional input manifold. To this end, Newman's spectral algorithm for community detection in a network [START_REF] Newman | Modularity and community structure in networks[END_REF] is used. A clear advantage of Newman's algorithm is that the number of clusters is not defined a priori, in contrast to more common clustering techniques (e.g. k-means). The number of clusters is instead the result of a maximization procedure performed on the modularity Q of the network. In other words, the number of thermodynamic clusters in the flow are determined only from the data and is not based on a priori knowledge of the user. This knowledge might be even impossible to come by in complex, unsteady hypersonic flows subjected to shocks.

Following this approach, the Euclidean distance matrix ∆ of the dataset in low-dimensional space is computed

∆ ij = ∥Y i -Y j ∥ 2 for (Y i , Y j ) ∈ Y 2 . ( 10
)
The dataset is subsequently recast into an undirected network with a binary adjacency matrix A, constructed as

A ij = 1 if ∆ ij < ϵ, 0 otherwise. ( 11
)
Two points are connected with an edge if their Euclidean distance is below a certain threshold ϵ. This threshold is usually chosen as a fraction of the mean of the distance matrix ∆. The influence of the threshold ϵ on the number of clusters will be investigated in Section IV A 2.

Finally, the dataset is progressively split into two communities until the modularity, Q, is maximized. The modularity is defined as the proportion of edges contained within a community over the same proportion for a random reference network. Let k i be the number of edges pointing towards the data point numbered i, and m the total number of connections within the network. Then, the probability of having an edge between i and j in the random reference network is k i k j /m. Hence, the modularity Q is defined as,

Q = 1 m ij A ij - k i k j m δ ci,cj , (12)
where δ ci,cj is the Kronecker delta for the communities of i and j. For instance, δ ci,cj = 1 only when i and j belong to the same community. By defining a vector s where s i = 1 if vertex i belongs to the first group, and s i = -1 otherwise, we can reformulate

Q = 1 2m ij A ij - k i k j m (s i s j + 1) = 1 2m s T Bs, (13)
where the modularity matrix B is given as,

B ij = A ij - k i k j m . ( 14
)
Since the graph is undirected, the modularity matrix B is symmetric and the modularity Q represents a Rayleigh quotient for matrix B. In order to maximize Q, we need to choose a vector s that is parallel to the principal eigenvector (corresponding to the largest eigenvalue) of B, v, which can be achieved by setting

s i = 1 if v i > 0 and s i = -1 if v i < 0.
In order to partition the graph into more than two communities, this algorithm is repeated until the modularity of each subgraph can no longer be increased. A thorough descripion of the full algorithm can be found in [START_REF] Newman | Modularity and community structure in networks[END_REF]. Fig. 6a shows the application of the clustering algorithm to the boundary layer data. The distance matrix ∆ is constructed on Y. After running the algorithm on the subsequent adjacency matrix A, two distinct clusters are identified, highlighted by the low distance between the points within a cluster in Fig. 6b.

Surrogate model construction

Finally, a surrogate surface is computed on the scattered low-dimensional points of each cluster. Many algorithms can be used to that end, such as kriging [START_REF] Kleijnen | Kriging metamodeling in simulation: A review[END_REF][START_REF] Bouhlel | Improving kriging surrogates of high-dimensional design models by partial least squares dimension reduction[END_REF], artificial neural networks [START_REF] Sun | A review of the artificial neural network surrogate modeling in aerodynamic design[END_REF], or radial basis functions [START_REF] Broomhead | Radial basis functions, multi-variable functional interpolation and adaptive networks[END_REF][START_REF] Buhmann | Radial basis functions[END_REF][START_REF] Powell | The theory of radial basis function approximation in[END_REF]. In the application of interest to this work, radial basis functions (RBF) provided the best trade-off between performance and accuracy, as well as an easy training step. In fact, the optimal weights of a RBF can be obtained through the solution of a linear system. However, it should be noted that our choice is not definitive and can be easily changed. Given a set of N R input points x 1 , ..., x N R ∈ R d and the function value at these points f (x 1 ), ..., f (x N R ) the radial basis function (RBF) interpolant g is given by

g(ϕ, x) = N R i=1 λ i ϕ(∥x -x i ∥) ( 15
)
where ϕ is the kernel function whose value depends on the distance r = ∥x -x i ∥ between the evaluation point x and the center x i of the RBF. In this study, the thin-plate spline kernel [START_REF] Wood | Thin plate regression splines[END_REF] was used, i.e.,

ϕ(r) = r 2 log(r). (16)
With this particular kernel, the kernel matrix is only conditionally positive definite. To ensure a unique solution for the interpolation weights, the system is augmented by a polynomial p ∈ Π d m (space of polynomials of d variables and degree up to m) to the right hand side of Eq. ( 15) [START_REF] Buhmann | Radial basis functions[END_REF], resulting in

g(x) = N R i=1 λ i ϕ(∥x -x i ∥) + p(x). ( 17
)
The extra degrees of freedom are accounted for by enforcing orthogonality of the coefficients with respect to the polynomial space as

N R i=1 λ i r(x) = 0, r ∈ Π d m . (18)
Finally, the polynomial coefficients c = [c 1 , ..., c d+1 ] T and the RBF coefficients Λ = [λ 1 , ..., λ N R ] T are found through the solution of the following linear system where )) will likely result in overfitting and prohibitive computational cost [START_REF] Schwenker | Three learning phases for radial-basis-function networks[END_REF]. Following the recommendations in [START_REF] Schwenker | Three learning phases for radial-basis-function networks[END_REF], the interpolant is constructed in two steps. First, the k-means algorithm with N R ≪ N clusters is applied on the concatenation of the input and output vector (X, Z) ∈ R N ×(D+D Z ) . The addition of the output vector results in a low within-cluster variance of the outputs and will ultimately improve the surrogate model. The N R centroids obtained with k-means, x c , are sent to the library to compute the function value vector f. Simultaneously, they are encoded in the lowdimensional space to obtain the N R cluster centers that will be used to train the RBF y c 1 , ..., y c N R ∈ R d . Fig. 7 shows the resulting tesselation in the embedded space after applying the k-means algorithm with N R = 250. The influence of the number of RBF centers on the quality of the surrogate model will be assessed in Section IV A 3.

Φ P P T 0 Λ c = f 0 d+1 (19)
P i = [1, x i , ..., x d i ] for i ∈ [1, N R ]
Following this approach, a single interpolant, g c k , is constructed for each cluster, c k ; in other words, g c k is the approximation of the scaled library function f on the low-dimensional subspace corresponding to cluster c k . The advantage of having one interpolant per cluster is twofold. First, it allows the surrogate model to best fit a region with a given dynamics of the high-dimensional function, especially in the presence of discontinuities, similar to the approach of Bettebghor et al. [START_REF] Bettebghor | Surrogate modeling approximation using a mixture of experts based on em joint estimation[END_REF]. Secondly, as the surrogate model spans a smaller range of input parameters, less centers are required to capture the given dynamics accurately, resulting in a lighter model with faster evaluation time. To enforce continuity of the surrogate surface near the cluster boundaries, the nearest centroids that do not belong to the considered cluster are added to its training set.

B. Coupling to the solver

Once the model is trained, we can replace the calls of the solver to the look-up library with the new lighter model. New points have to go through three steps: (i) Out-of-sample dimensionality reduction, (ii) classification, and (iii) interpolation. These three steps will be described in the following. Let X t ∈ R Nt×D denote the stack of all new points.

Out-of-sample encoding

The low-dimensional representation of X t (after proper scaling of Xt ) is straightforward. Indeed, the point is simply fed to the encoder portion of the input/output-encoder

Y t = W E,2 σ(W E,1 X t + b E,1 ) + b E,2 . (20)
This results in a fast and inexpensive encoding of the new out-of-sample points. In fact, the time complexity of the encoding step is

O(H × C ac × L × N t )
, where H is the maximum number of neurons in a layer, C ac is the complexity of the activation function and L is the number of layers in the encoder step of the input/output-encoder.

Classification

The next step is to determine to which community the new state Y t belongs. To this end, after applying Newman's algorithm, a random forest classifier is trained on the resulting clusters. A random forest classifier, formally proposed by Breiman [START_REF] Breiman | Random forests[END_REF], is a collection of n tree tree-based classifiers, h(x, Φ k ), k = 1, ..., n tree , where Φ k are identically distributed random vectors. Each tree votes for the most likely class of input vector x, and the majority wins. Fig. 8 shows the confusion matrix C of the classifier trained on the two clusters obtained above with n tree = 20. The clusters of the embedded new points Y t , not seen during the training-phase of the classifier, are predicted and compared to the true clusters (given by Newman's algorithm). The off-diagonal values count the number of points that are assigned to the wrong cluster. All off-diagonal values are zero, demonstrating the ability of the classifier to correctly predict the community of an out-of-sample point. The time complexity of the induction time of the classifier is O(n tree × N t log(N t )) (see the book by Witten, Frank and Hall [START_REF] Frank | Data mining in bioinformatics using weka[END_REF] for a demonstration), where n tree is the number of decision trees in the random forest. The logarithmic term accounts for the worst case scenario for the maximum depth of each tree. However, the maximum depth is usually set to a smaller value, resulting in a complexity of O(n tree × depth × N t ), which results in relatively fast classification times.

Interpolation

Finally, once Y t has been found to belong to cluster c k , the corresponding RBF g c k is called to evaluate the thermochemical properties Z t of the mixture at that state X t ,

Z t = g c k (Y t ) = N R i=1 λ i ϕ(Y t -Y c i ). ( 21
)
The hat denotes the predicted value by the reduced library, as opposed to the true value that would have been given by the target library, Z t . Finally, these properties are re-scaled according to ˜ Z t = Z t * ( Zmax -Zmin ) + Zmin [START_REF] Broomhead | Radial basis functions, multi-variable functional interpolation and adaptive networks[END_REF] and passed back to the flow solver. The time complexity of each surrogate model is given by O(C RBF × N t × N R × d), where C RBF is the complexity of the kernel function. The time limiting part of the RBF interpolation is the calculation of the distance matrix that scales with O(N t × N R × d). Hence, using a relatively small number N R of RBF centers, compared to N t , will greatly improve the performance of the surrogate model. It should be noted that the dimensionality reduction directly contributes to the added performance of the surrogate model as d < D.

Global performance

The time complexity of the whole algorithm can be recovered by adding the time complexity of each of the three steps. Hence, the total time complexity of the algorithm can be written as O(C M L N ), where C M L = O(HC ac L + n tree depth + dN R C RBF ).

III. APPLICATION TO HYPERSONIC FLOWS IN CHEMICAL NON-EQUILIBRIUM

In this section, the general equations governing hypersonic flows in chemical non-equilibrium are first recalled, followed by a brief description of the numerical framework used to solve these equations. Then, the two benchmark cases, a Mach-10 adiabatic laminar boundary layer in chemical non-equilibrium, initially studied by Marxen et al [START_REF] Marxen | A method for the direct numerical simulation of hypersonic boundary-layer instability with finite-rate chemistry[END_REF][START_REF] Marxen | Nonlinear instability of a supersonic boundary layer with two-dimensional roughness[END_REF], and a Mach-5.92 shock wave boundary layer interaction are presented.

A. Governing conservation equations

The nondimensional Navier-Stokes equations for a mixture of multiple species S are presented in Eqs. [START_REF] Buhmann | Radial basis functions[END_REF][START_REF] Jin | Comparative studies of metamodelling techniques under multiple modelling criteria[END_REF][START_REF] Peter | Comparison of surrogate models for the actual global optimization of a 2d turbomachinery flow[END_REF][START_REF] Chang | Flood forecasting using radial basis function neural networks[END_REF].

∂ρ ∂t + ∇ • (ρu) = 0 ( 23
)
∂ρ s ∂t + ∇ • (ρ s u + ρ s V s ) = ωs , ∀ s ∈ S (24)
∂ρu ∂t + ∇ • (ρu ⊗ u) = -∇p + ∇ • τ ( 25
)
∂ρe 0 ∂t + ∇ • (ρh 0 u) = ∇ • (τ • u) -∇ • q (26)
Equation ( 23) is the continuity equation, describing mass conservation in the system. Equation ( 24) corresponds to the set of mass conservation equations for each species, with the net production rate terms, ωs , appearing on their right-hand side. In order to ensure global mass conservation, in the case of a finite-rate reacting mixture with a varying composition, Eq. ( 23) needs to be solved together with Eq. ( 24) for all but one species. The omitted species is selected based on numerical considerations, avoiding species with the smallest concentrations. The nondimensional quantities are the time, t, the density, ρ, the velocity, u, the pressure, p, the stress tensor, τ, the total energy, e 0 , the total enthalpy, h 0 , the heat flux, q, as well as, the partial density, ρ s , the net mass production rate, ωs , and the diffusion velocity, V s , for the species s. More details regarding the derivation of the equations and the validity of our invoked assumptions are provided in [START_REF] Anderson | Hypersonic and High-Temperature Gas Dynamics[END_REF][START_REF] Gnoffo | Conservation equations and physical models for hypersonic air flows in thermal and chemical nonequilibrium[END_REF][START_REF] Josyula | Fundamental fluid transport equations for hypersonic nonequilibrium flows[END_REF]. There are N s -1 conservation equations for the species, so the total problem involves N s + 4 equations. The stress tensor, τ, and the heat flux, q, are computed as

τ = µ Re ∞ ∇u + (∇u) T -(∇ • u) I , (27)
q = - κ Re ∞ P r ∞ Ec ∞ ∇T + s∈S ρ s h s V s . (28)
The nondimensional Reynolds number, Re ∞ , and Prandtl number, P r ∞ , are defined in the free stream, at the domain inlet. The Eckert number, Ec ∞ , is also computed at the free stream and is equal to one by design. These nondimensional quantities are defined with dimensional quantities (distinguishable by the tilde) as

Re ∞ = ρ∞ c∞ Lref μ∞ , P r ∞ = μ∞ cp ∞ κ∞ , Ec ∞ = c2 ∞ cp ∞ T ref (29)
where T ref = (γ ∞ -1) T ∞ , and c stands for the speed of sound. µ denotes the dynamic viscosity, κ the frozen thermal conductivity, and c p the specific heat at constant pressure.

B. Thermodynamic and chemical models for a mixture in chemical non-equilibrium

In the following, the closure of the governing equations when assuming a mixture in chemical non-equilibrium is detailed. In general, a reacting gas in a high-enthalpy flow is considered as a multi-component mixture that consists of a set of species S interacting through a defined network of reactions. The presence of the species makes the chemical reaction and diffusion terms in the governing equations significant, which in turn require to be modeled using various assumptions.

For a multi-component gas mixture, the global mixture properties are derived from the species properties based on

ρ = s∈S ρ s , e = s∈S Y s e s , h = s∈S Y s h s ( 30
)
where the species mass fraction is Y s = ρ s ρ , with s∈S Y s = 1. The mixture thermodynamic and transport properties depend, generally, on their composition, and any two thermodynamic properties, for example the temperature and the pressure, which define the thermodynamic state of the mixture. The composition is commonly taken as an independent variable. It can potentially become dependent on other thermodynamic quantities for the special cases of frozen chemistry or chemical equilibrium. However, these cases are not described in this study. The individual species properties are accurately computed by kinetic theory and statistical mechanics [START_REF] Scoggins | Development of Mutation++: Multicomponent Thermodynamic and Transport Properties for Ionized Plasmas written in C++[END_REF][START_REF] Scoggins | Development of numerical methods and study of coupled flow, radiation, and ablation phenomena for atmospheric entry[END_REF]. Hence the composition and the two thermodynamic properties can be concatenated into the thermodynamic state vector qth . This defines the relation Z = f (q th ) on which the algorithm is based.

When finite-rate chemistry is not neglected, the species mass production rate, ωs , and the species diffusion velocities, V s , need to be modeled. For a general case, a set of reactions, R, is considered depending on the mixture in question. Each reaction r is characterized by a reaction rate, R r , which is computed by the forward rate, k f r , and the backward rate, k br . These rates, in turn, are obtained according to experimentally or theoretically calibrated Arrhenius formulas in the form k f r = C r T nr exp (T ar /T ), and k br = k f r /K eq r (T ), where K eq r is the reaction equilibrium constant at the specific conditions. This description is given in Eq. ( 31) for a generic reaction.

Reaction (r):

s∈S ν ′ r,s S s k f r --⇀ ↽ -- k b r s∈S ν ′′ r,s S s (31)
The net reaction rate is then given by

R r = k f r Π s ρ s M s ν ′ r,s -k br Π i ρ s M s ν ′′ r,s • s∈S Z r,s ρ s M s ( 32
)
where the species molar mass is M s and the efficiency of species s as a third-body in reaction r is Z r,s . The net mass production rates for species s from all reactions are obtained as

ωs = M s r∈R ν ′′ r,s -ν ′ r,s R r . ( 33
)
The diffusion flux J s = ρ s V s appearing in Eq. ( 24) also needs to be modeled. The accurate and rigorous Stefan-Maxwell multicomponent diffusion model provides the diffusion velocities V s as the solution of a constrained linear system of equations [START_REF] Scoggins | Development of numerical methods and study of coupled flow, radiation, and ablation phenomena for atmospheric entry[END_REF]. In that framework, the local molar fraction gradients are needed, which cancels the purely local assumption of the input/output problem. Hence, in this study, we use a simple expression based on Fick's law with a mass correction [START_REF] Hirschfelder | Molecular theory of gases and liquids[END_REF][START_REF] Ramshaw | Self-consistent effective binary diffusion in multicomponent gas mixtures[END_REF] given by

J s = -cM s D s ∇Y s + cY s i∈S M i D i ∇Y i . (34)
Here, c = s∈S (ρ s /M s ), and D s is the averaged diffusion coefficient for species s, defined below

D s = 1 -X s r̸ =s X r /D s,r . (35)
The thermochemical library Mutation++ [START_REF] Scoggins | Mutation++: Multicomponent thermodynamic and transport properties for ionized gases in c++[END_REF] is used to compute thermodynamic and transport properties at different conditions. The reader is referred to the description provided in [START_REF] Scoggins | Mutation++: Multicomponent thermodynamic and transport properties for ionized gases in c++[END_REF][START_REF] Scoggins | Development of numerical methods and study of coupled flow, radiation, and ablation phenomena for atmospheric entry[END_REF] for further details on the library.

C. Numerical framework

The governing equations are solved using a high-fidelity in-house flow solver for the direct numerical simulation of hypersonic flows in chemical non-equilibrium, as described in [START_REF] Margaritis | High-fidelity computational tool for chemically reacting hypersonic flow simulations[END_REF]. Compact schemes are used for a spatial discretization, and explicit time-integration is performed using Runge-Kutta schemes. More details about the discretization can be found in the aforementioned paper.

Coupling with Mutation++

The thermodynamic and transport properties, as well as the source terms resulting from the chemical kinetics models, are extracted from the Mutation++ library [START_REF] Scoggins | Mutation++: Multicomponent thermodynamic and transport properties for ionized gases in c++[END_REF] to close the governing equations, Eqs. [START_REF] Buhmann | Radial basis functions[END_REF][START_REF] Jin | Comparative studies of metamodelling techniques under multiple modelling criteria[END_REF][START_REF] Peter | Comparison of surrogate models for the actual global optimization of a 2d turbomachinery flow[END_REF][START_REF] Chang | Flood forecasting using radial basis function neural networks[END_REF]. The library, originally written in C++, is coupled with the numerical solver using a wrapper in Fortran 95.

For the remainder of this study, we will only consider an air mixture composed of five species S = [N 2 , O 2 , N O, N, O]. However, we stress that the algorithm can be applied to any gas mixture. At each grid point, we have access to the local thermodynamic state and composition qth (in dimensional units),

qth = ρ ρe ρ N ρ O ρ N O ρ O2 . ( 36
)
qth is then passed to Mutation++, which returns the necessary thermochemical properties needed to close the reactive compressible Navier-Stokes Eqs. [START_REF] Buhmann | Radial basis functions[END_REF][START_REF] Jin | Comparative studies of metamodelling techniques under multiple modelling criteria[END_REF][START_REF] Peter | Comparison of surrogate models for the actual global optimization of a 2d turbomachinery flow[END_REF][START_REF] Chang | Flood forecasting using radial basis function neural networks[END_REF],

z = µ κ P T h s ω s D s s ∈ S. (37)
Hence, calls to Mutation++ can be seen as an input/output problem z = f (q th ) where the function f represents the library.

D. Case 1: Mach-10 adiabatic boundary layer

In this section, the first simple test-case used to showcase the applicability and performance of the algorithm is presented. We simulate a Mach-10 adiabatic flat-plate boundary layer in Earth atmosphere, based on the work of Marxen et al [START_REF] Marxen | A method for the direct numerical simulation of hypersonic boundary-layer instability with finite-rate chemistry[END_REF][START_REF] Marxen | Nonlinear instability of a supersonic boundary layer with two-dimensional roughness[END_REF]. The thermodynamic and freestream conditions are presented in Table I. Details of the numerical grid can be found in [START_REF] Margaritis | High-fidelity computational tool for chemically reacting hypersonic flow simulations[END_REF]. This case has already been validated against literature with the current flow solver in [START_REF] Margaritis | High-fidelity computational tool for chemically reacting hypersonic flow simulations[END_REF]. Due to the wall temperature approaching T wall ≈ 4900K near the inflow, N 2 and O 2 rapidly start to produce their monoatomic counterpart, as well as N O, through endothermic chemical reactions. This is illustrated in Fig. 9c presenting all the mass fraction profiles at the streamwise location of Re x = 2000. Close to the wall, the O 2 mass fraction decreases while O and N O are created. To a smaller extent, N is also created through N 2 dissociation. The mean flow and temperature profiles at a streamwise Reynolds number of Re x = 2000 are presented in Fig. 9a and Fig. 9b, alongside the same case simulated using a thermally perfect gas assumption (ignoring non-equilibrium effects). The base flows differ significantly depending on the assumption used for the gas. Wall temperature decreases significantly from a thermally perfect gas to a finite-rate chemistry assumption, and the boundary layer becomes thicker. Overall, these results highlight and corroborate the importance of including chemical non-equilibrium effects for the accurate simulation of high-enthalpy hypersonic flows. On another note, it has been shown in [START_REF] Margaritis | High-fidelity computational tool for chemically reacting hypersonic flow simulations[END_REF] that the choice of the diffusion model (i.e Stefan-Maxwell or an approximation based on Fick's law (Eq. ( 34))) does not have any noticeable impact on the baseflow obtained in a chemical non-equilibrium simulation. However, using the algebraic equation reduces the CPU cost of the simulation by 40% compared to the rigorous Stefan-Maxwell model. While the adiabatic boundary layer is in a chemical non-equilibrium regime, it lacks typical compressible flow features such as shocks and expansion fans that could challenge the performance of the algorithm. To verify the applicability of the developed technique, a second case, namely a shock wave boundary layer interaction (SBLI), is proposed. This case has been initially studied by Margaritis et al. [START_REF] Margaritis | High-fidelity computational tool for chemically reacting hypersonic flow simulations[END_REF]. The thermodynamic and freestream conditions are presented in Table II where x 0 denotes the shock impinging location and θ denotes the shock angle. The wall is considered as adiabatic. Details of the numerical grid and the computational setup can be found in [START_REF] Margaritis | High-fidelity computational tool for chemically reacting hypersonic flow simulations[END_REF]. Fig. 10a and Fig. 10b present the wall distribution of skin-friction coefficient and wall pressure, respectively. The recirculation bubble is smaller in the case in chemical non-equilibrium compared to the case assuming perfect gas. This is explained by the high concentration of atomic species in the recirculation bubble, as seen in Fig. 10c. The endothermic reactions extract energy out of the bubble. Moreover, the skin-friction and pressure are higher just upstream of the reattachment point compared to the perfect-gas simulation. Hence, this case also highlights the effect of chemical non-equilibrium on the resulting flow features.

Test case 2 M∞ 5.92 Rex 0 1.15 × 10 6 T∞ [K] 1110.5 p∞ [P a] 61 × 10 3 θ [deg]
13.0 TABLE II. Thermodynamic and free stream conditions for Mach-5.92 shock-wave boundary layer case (from [START_REF] Margaritis | High-fidelity computational tool for chemically reacting hypersonic flow simulations[END_REF]).

IV. RESULTS

The results are presented in three parts. First, each step of the off-line training is assessed, and then the performance of the model in predicting the output quantities when coupled to the Navier-Stokes solver is analyzed on the Mach 10 adiabatic boundary layer. Finally, the same strategy is applied to a shock-wave boundary layer interaction scenario.

A. Model training N = 100, 000 thermodynamic state vectors qth are sampled from the converged solution and concatenated into the input dataset X. The corresponding Mutation++ output dataset Z is also obtained.

Input/output encoding

The architecture and hyperparameters of the input/output-encoder used for this test case are provided in Table III An important parameter is the number of dimensions required to properly unfold the input manifold. To this end, the architecture and hyperparameters, as well as the random seed to initialize the network weights, are held constant, and only the dimension of the latent space is varied. We follow the reconstruction loss ∥ Z -Z∥ 2 of the testing set as a function of the latent space dimension d. As seen in Fig. 11, the loss saturates after d = 2, suggesting that two dimensions are sufficient to represent the input manifold, originally in R 6 .

Newman's community clustering

The second step is to cluster the data on the low-dimensional manifold using Newman's algorithm. The only hyperparameter needed for Newman's clustering algorithm is the threshold ϵ to determine the adjacency matrix A from the distance matrix ∆, as described in Section II A 3. To showcase the robustness of the number of clusters with respect to the threshold ϵ, the algorithm is applied over a range of thresholds, chosen as multiples of the mean of the distance matrix ∆. As shown in Fig. 12, over the range tested, the algorithm returns N c = 2 clusters before over-fitting with extreme values. Fig. 13a shows the two clusters obtained in the embedded space. As expected, the two clusters define different regions in the reduced space. To gain more physical insight, Fig. 13b shows randomly selected points of the training set, mapped back to their original location in the Cartesian space and colored by their cluster number. Contours of temperature T in Kelvin are added. Each cluster fills a region of the flow with different levels of chemical non-equilibrium. The blue cluster represent the free stream, where temperature is low and chemistry is frozen. Alternatively, the red cluster corresponds to the near-wall region with dissociated species and high temperatures.

The random forest classifier has been trained with n tree = 20. This number proved sufficient to obtain an acceptable prediction accuracy of new points clusters, as seen in Fig. 8 above.

Surrogate model construction

Once the clusters are determined, a separate RBF interpolant is trained for each cluster as described in Section II A 4. To enforce continuity of the surrogate through the clusters, the nearest centroids shared between the two clusters are added to the training set. To assess the number of centroids needed, the RBFs are trained simultaneously, for both clusters, with the same number of centroids N R (note that this number can be varied to accommodate clusters of different sizes). The error of the model for the testing set ∥ Z -Z∥ 2 is plotted against the number of centroids N R in Fig. 14 (blue curve). As expected, the error decreases as the number of centroids is increased. To assess potential overfitting, one can follow the evolution of the mean of the squared RBF coefficients Λ 2 ,

Λ 2 = 1 N R N R i=1 λ 2 i ( 38
)
where a high value will likely indicate overfitting. On the right axis of Fig. 14, we see that Λ 2 increases with the number of RBF centroids. Therefore, choosing the right N R is a trade-off between low error and low overfitting. In this case, a value of N R = 250 for each cluster has been retained.

B. Model accuracy

The reduced library is tested (off-line) on a full snapshot (which includes also Fig. 15 displays the relative error (as a percentage) between the temperature T given by Mutation++ and the prediction of the data-driven model T , for the four configurations enumerated above; contours of temperature in Kelvin are also added to highlight the evolution of the flow. First, the prediction of model 1 produces the highest error by a wide margin. This showcases the difficulty of properly training a neural network for prediction in high dimensions. The figure shows the maximum relative error in all three remaining models to be only a few percent, and located around the edge of the boundary layer, where the gradients are strongest. In addition, the error in model 2 is higher than the error of model 3, even though some information is lost in the latter due to the encoding step. This is a direct consequence of the curse of dimensionality: as the number of dimensions increases, the sampling volume in the input space increases exponentially. However, in this case, we kept the number of RBF centers N R fixed. To get to the same level of accuracy, N R should be increased in model 2, resulting in a performance loss. Finally, the error decreases slightly from model 3 to model 4. This improvement originates from the clustering step. In fact, each of the two clusters have N R = 250 centers. Hence, the input space is actually populated with N R = 500 RBF centers in model 4. According to Fig. 14, this likely improves the accuracy of the surrogate surface. The added number of centroids, however, does not result in a loss of performance. Assuming that a fraction α of the N t query points are in cluster 1, then 1 -α query points are in cluster 2. The evaluation time of the surrogate surface linked to cluster 1 then scales roughly as αdN t N R C RBF . Similarly, for cluster 2, it scales as (1 -α)dN t N R C RBF . Finally, the total evaluation time, i.e., the sum of the two, remains dN t N R C RBF . This can be easily generalized to a higher number of clusters.

These results demonstrate that the preprocessing steps involved in the construction of the model improve overall performance while maintaining a high level of accuracy.

C. Model stability

The resulting data-driven model (model 4, with all pre-processing steps) is coupled to the flow solver in a timemarching simulation. Starting from the solution obtained with Mutation++, the simulation is restarted using the reduced library only, also referred as a closed-loop prediction. After running for a couple of flow-through times, the solution remains stable. The base-flow profiles are compared for various quantities of interest in Fig. 16. Excellent agreement is found between the profiles. This validates the accuracy and suitability of the data-driven model to simulate hypersonic flows in chemical non-equilibrium over the enthalpy range observed during the training step.

D. Model performance

To compare the performance of the data-driven model to the full library, we performed a scaling study. Mutation++ is a serial library, hence its time complexity can be expressed as O(C M ++ N t ) where N t is the number of independent, evaluated thermodynamic states. Two variants of Mutation++ are considered here. The first one solves directly the Stefan-Maxwell diffusion problem and returns the diffusion velocity. In the second one, the diffusion coefficients D s are returned and the diffusion fluxes are later computed using Eq. [START_REF] Bhagatwala | Direct numerical simulations of hcci/saci with ethanol[END_REF]. The diffusion fluxes are computed in the same fashion with the data-driven model. Moreover, we recall that its time complexity is O(C M L N t ), where C M L = O(HC ac L + n tree depth + N R C RBF ). For all three thermochemical models, the prefactor is empirically determined in Fig. 17 All curve fits, shown in Fig. 17, suggest that, in practice, the models scale as O(N t ) with exponents close to unity. The ratio of the prefactor is C M L /C M ++ ≈ 0.45 when using Mutation++ in its local version. We can therefore expect a 55% CPU gain by using the data-driven model instead of Mutation++. In fact, we assessed a CPU time reduction of 50% during the simulation, with a grid of size N ≈ 400, 000. This confirms a speedup through the use of a surrogate model. Secondly, the speed-up is even more significant when the library also solves for the Stefan-Maxwell diffusion problem at each grid point. Although not rigorously a one-to-one model comparison anymore, the data-driven model now performs 70% faster without any loss of accuracy. In fact, Fick's law based diffusion model have been shown to be highly accurate in hypersonic simulations. Moreover, fine tuning of the hyperparameters may allow even higher CPU gains as C M L is proportional to a linear combination of the hyperparameters. Finally, we stress that the data-driven algorithm is a python implementation competing with a compiled C++ library. The speed-up reported here can be significantly increased by porting the model to a compiled language. It is also believed that the speed-up would be more significant when the dimensionality of the input space increases to include more chemical species. We thus expect even larger CPU gains in the future with an optimized implementation and added adaptivity.

E. Application to the SBLI case

Following the same steps, the model is trained on the SBLI case. The model has the following specifications: d = 3, c = 3, N R = 250. In this case, the dimensions of the latent space and the number of clusters are higher due to the more complex thermodynamic manifold learned by the IO-E. A two dimensional projection of the three dimensional manifold is presented in Fig. 18a. In this plane, distinct thermodynamic regions (i.e. different clusters) are wrapped around a scarcely populated center area. This can be explained by the impinging, recirculation and reflected shocks that induce abrupt change in the thermodynamic state. These regions and their borders become even more meaningful when reported to their physical location in the flow, as seen on Fig. 18b where a numerical Schlieren is superposed. The green cluster corresponds to mildly hot conditions with high density, i.e. the freestream and post impinging shock conditions. After the recirculation shock, the thermodynamic states shift instantaneously to higher densities and temperatures, represented by the blue cluster. However, close to the apex of the recirculation bubble, the expansion fan decreases these thermodynamic variables, inducing a shift back to the green cluster. Finally, the red cluster, found in the boundary layer, corresponds to high temperatures and low densities. At the core of the recirculation bubble, temperature decreases and density increases, which brings the local state vector back to the green cluster.

In a closed-loop simulation, the model remains stable while maintaining a factor 2 speed-up in predicting thermochemical properties. In fact, Fig. 19b and Fig. 19a show that the wall pressure and skin friction remain in excellent agreement with the baseflow solution after 2 flow-through times. The only discrepancy with the initial solution is observed for atomic nitrogen mass concentration in Fig. 19c. However, it is present in such small quantities that it does not have an impact on the stability of the solution. In fact, the thermochemical properties are not sensitive to small perturbations in atomic nitrogen concentration (cf. Fig. 3), another fact motivating the dimensionality reduction performed in pre-processing.

V. CONCLUSIONS

In this paper, we presented a novel technique to reduce any high-dimensional look-up library to a lower-dimensional surrogate, and thus reduce the CPU costs of numerical simulations that rely on these libraries. Several machine learning techniques have been used: encoding based on deep neural networks, community clustering, surrogate modeling and classification in a three-step learning phase. In the first step, the proposed input/output-encoder architecture has been shown to outperform partial least-squares (PLS) for dimensionality reduction of input/output relations. Clustering was performed using Newman's algorithm. It discovered physically consistent clusters in the low-dimensional latent space without a-priori knowledge of the number of clusters. Then, a random-forest classifier was trained, which reliably predicted the cluster of previously unencountered data points. Finally, a radial basis function network was constructed on each cluster to obtain a continuous and local representation of the library via a reduced-order surrogate model. The combination of these pre-processing steps has been shown to improve the efficiency of the model on our two test cases: a Mach-10 adiabatic boundary layer and a Mach-5.92 shock wave boundary layer interaction, both in chemical non-equilibrium. After training, the model replaced Mutation++ and converged rapidly to a new solution. The newly computed base flows were recovered accurately when compared to the true solution (obtained with Mutation++). During this demonstration, we observed up to 70 % CPU time decrease to compute the thermochemical properties of the mixture. This computational framework can be readily ported into other application fields to accelerate simulations that rely on high-dimensional look-up tables to model complex flow behavior such as combustion, phase-change or fluid-particle interactions. Finally, future steps in algorithmic development will include on-the-fly adaptivity of the model to tackle unsteady flow problems, and leveraging analytical solutions during the initial training set. Application to ablation and more detailed chemical mechanisms, for example those appearing in reactive flows, are also being pursued at the moment in the context of future work.

34 5 4 6 "FIG. 2 .

 62 FIG. 2. General schematic of the model training and coupling to replace any expensive Input/Output library.

FIG. 3 .

 3 FIG. 3. Numerical range of selected Mutation++ outputs Z (vertical) with respect to the inputs X (horizontal).

FIG. 4 .

 4 FIG. 4. Low-dimensional representation Y ∈ R N ×d (d = 2) of X ∈ R N ×D , colored by temperature T . Obtained with (a) PCA (b) AE (c) PLS (d) IO-E.

FIG. 6 .

 6 FIG. 6. Demonstration of the Newman algorithm. (a) Original distance matrix ∆ of dataset Y. (b) Restored communities from running Newman's algorithm on A, showing two distinct clusters

FIG. 7 .

 7 FIG. 7. Tesselation of the low-dimensional space after applying the k-means algorithm with NR = 250 on the concatenation of the input X ∈ R N ×D and output Z ∈ R N ×D Z vectors. Black stars represent the cluster centroids' low-dimensional representation y c 1 , ..., y c N R .

  is the polynomial matrix, and f = [f (x 1 ), ..., f (x N R )] denotes the vector containing the function values at the RBF centers. Due to the large size of the training set, using one cluster center per training point (i.e. N R = N = O(10 5

FIG. 8 .

 8 FIG.8. Confusion matrix C of the random forest classifier (ntree = 20), tested on unseen data points during training. Cij is the number of points belonging to cluster cj predicted to be in cluster ci.

FIG. 9 .

 9 FIG. 9. Profiles of (a) streamwise velocity, (b) temperature, (c) species mass fractions from left to right N , N O, O, O2 and N2, at Rex = 2000. (a,b) Red line: simulation in chemical non-equilibrium, black line: perfect gas assumption.

FIG. 10 .

 10 FIG. 10. Wall distribution of (a) skin-friction coefficient, (b) pressure, and (c) species mass fractions (from top to bottom on left side: N2, O2, O, N O and N ). (a,b) Red line: simulation in chemical non-equilibrium, black line: perfect gas assumption.

2 FIG. 11 .

 211 FIG. 11. Reconstruction loss ∥ Z -Z∥2 of the embedding done by IO-E with respect to the number of latent space dimensions d.

2 FIG. 12 .

 212 FIG. 12. Number of clusters Nc obtained with Newman's algorithm as a function of the threshold ϵ, given as multiples of the distance matrix mean ∆.

FIG. 13 .

 13 FIG. 13. Training points Y colored by their cluster number. (a) In the latent space found by IO-E. (b) At the Cartesian location they were sampled from, with contours of temperature T in Kelvin.

  the training points used to build the model) to assess the capacity of the model to interpolate new points not encountered during training. Four configurations of the data-driven model are tested: (i) model 1, using the full IO-E for prediction (ii) model 2, with no dimensionality reduction and no clustering (d = 6, N c = 1, N R = 250), (iii) model 3, with dimensionality reduction, but without clustering, (d = 2, N c = 1,N R = 250), and (iv) model 4, with both dimensionality reduction and clustering (d = 2, N c = 2,N R = 250).

FIG. 14 .

 14 FIG. 14. Left axis: reconstruction error of surrogate model ∥ Z -Z∥2 on a testing set (dashed line). Right axis: mean of squared RBF coefficients Λ 2 (dotted line). Both with respect to the number of RBF centers NR used in training.

FIG. 15 .

 15 FIG. 15. Comparison of the relative temperature error | T -T |/T in percent with contours of temperature. (a) Model 1: full IO-E (b) Model 2: d = 6, Nc = 1, NR = 250 (c) Model 3: d = 2, Nc = 1, NR = 250 (d) Model 4: d = 6, Nc = 2, NR = 250

FIG. 16 .

 16 FIG. 16. Comparison of profiles of (a) streamwise velocity, (b) temperature, (c) species mass fractions from left to to right N , N O, O, O2 and N2 at Rex = 2000. Solid line and symbols correspond to the solution obtained using Mutation++ and the data-driven model, respectively.

  .

FIG. 17 .

 17 FIG. 17. Comparison of the time complexity of Mutation++ with Stefan-Maxwell diffusion (dotted black line), the purely local version (solid red line) and the data-driven model (dash dotted blue line). The best non-linear least-squares fit of the form CN α is added.

FIG. 18 .

 18 FIG. 18. Training points Y colored by their cluster number. (a) In the latent space found by IO-E. (b) At the Cartesian location they were sampled from, with contours of magnitude of the density gradient ∥∇ρ∥

FIG. 19 .

 19 FIG.[START_REF] Cybenko | Approximation by superpositions of a sigmoidal function[END_REF]. Comparison of wall-distribution of (a) skin-friction coefficient, (b) pressure, and (c) species mass fractions (from top to bottom on left side: N2, O2, O, N O and N ). Solid line and symbols correspond to the solution obtained using Mutation++ and the data-driven model, respectively.

TABLE I .

 I Thermodynamic and free stream conditions for Mach-10 adiabatic flat-plate boundary layer test case (from[START_REF] Marxen | A method for the direct numerical simulation of hypersonic boundary-layer instability with finite-rate chemistry[END_REF]).

	Test case 1
	M∞	10
	Re∞	10 5
	T∞[K]	350
	p∞[P a]	3596

TABLE III .

 III . Architecture and hyperparameters of the input/output-encoder used to train the IO-E network on the M a = 10 adiabatic flat plat boundary layer test case.

			Architecture
		Encoder		Decoder
	Layer		Size	Layer	Size
	Input		6	Latent space	2
	Fully connected	12	Fully connected	6
	Fully connected	6	Fully connected	12
	Latent space		2	Output	18
			Hyperparameters
	Parameter	Value		Parameter	Value
	Learning rate	1 × 10 -3		Epochs	2000
	Loss	Mean-squared error	Batch size	256
	Activation function	tanh		Optimizer	Adam (keras default)

ACKNOWLEDGMENTS

This work was supported by the Imperial College London -CNRS PhD Joint Program and was granted access to the HPC/AI resources of TGCC under the allocation 2021-A0102B12426 and 2022-A0122B13432 made by GENCI. Part of the calculations were also performed using MeSU computing platform at Sorbonne University.