Romain Dupuis

Jean-Christophe Jouhaud

Pierre Sagaut

Surrogate Modeling of Aerodynamic Simulations for Multiple Operating Conditions Using Machine Learning

This paper describes a methodology, called local decomposition method, which aims at building a surrogate model based on steady turbulent aerodynamic fields at multiple operating conditions. The various shapes taken by the aerodynamic fields due to the multiple operation conditions pose real challenges as well as the computational cost of the high-fidelity simulations. The developed strategy mitigates these issues by combining traditional surrogate models and machine learning. The central idea is to separate the solutions with a subsonic behavior from the transonic and high-gradient solutions. First, a shock sensor extracts a feature corresponding to the presence of discontinuities, easing the clustering of the simulations by an unsupervised learning algorithm. Second, a supervised learning algorithm divides the parameter space into subdomains, associated to different flow regimes. Local reduced-order models are built on each subdomain using proper orthogonal decomposition coupled with a multivariate interpolation tool. Finally, an improved resampling technique taking advantage of the subdomain decomposition minimizes the redundancy of sampling. The methodology is assessed on the turbulent two-dimensional flow around the RAE2822 transonic airfoil. It exhibits a significant improvement in terms of prediction accuracy for the developed strategy compared with the classical method of surrogate modeling.

O VERALL aircraft design and optimization rely increasingly on numerical simulations for structural, aerodynamics, or even noise analysis. Particularly, the computational fluid dynamics (CFD) plays a significant role in solving Navier-Stokes equations, in order to predict vector-valued functions of specific quantities of interest, such as wall pressure field. The equations are discretized into algebraic systems that lead to prohibitive computational cost for simulations with a high number of degrees of freedom. Moreover, the inflow conditions may vary and form a multidimensional parameter space. Its full exploration requires the computation of a very large number of expensive simulations and becomes intractable. One of the main solutions to overcome this problem is the substitution of the high-fidelity simulations by a mathematical approximation much faster to be run, referred to as a surrogate model. It represents an interesting trade-off between precision and computation time. Furthermore, reducing the computational time of the exploration for high-fidelity CFD can open the way to multiphysics simulations using surrogate models for the fluid parts.

The surrogate modeling of high-dimensional vector-valued functions is mainly performed with a reduced-order approach, called reduced-order modeling (ROM). Originally developed for the study of coherent structures in the turbulent boundary layer [START_REF] Aubry | The Dynamics of Coherent Structures in the Wall Region of a Turbulent Boundary Layer[END_REF], ROM methods have shown various applications, such as aeroelasticity [START_REF] Amsallem | Towards Real-Time Computational-Fluid-Dynamics-Based Aeroelastic Computations Using a Database of Reduced-Order Information[END_REF], optimal flow control [START_REF] Cordier | Two Typical Applications for POD: Coherent Structures Education and Reduced Order Modelling[END_REF], turbulent flows [START_REF] Berkooz | The Proper Orthogonal Decomposition in the Analysis of Turbulent Flow[END_REF][START_REF] Östh | On the Need for a Nonlinear Subscale Turbulence Term in POD Models as Exemplified for a High-Reynolds-Number Flow over an Ahmed Body[END_REF], or geophysics [START_REF] Ştefănescu | Model Reduction and Inverse Problems and Data Assimilation with Geophysical Applications. A Special Issue in Honor of I. Michael Navon's 75th Birthday[END_REF]. Most ROM methods are applied to CFD problems by approximating the high-fidelity model as a linear combination of low-dimensional basis vectors, weighted by purposely tuned parameters. The basis vectors characterize the main features of the system behavior. Proper orthogonal decomposition (POD) is a particular and very popular method of dimension reduction used very frequently for CFD problems. It computes the basis vectors and the corresponding modal coefficients with an optimal least-square approach from a given number of high-fidelity computations, also called snapshots, at different state parameters. The ROM coefficients are calculated only for a finite and discrete number of input parameters, whereas the surrogate model is evaluated on a new parameter set. Thus a continuous representation of the coefficients over the whole state-parameter space has to be provided in order to build the final model, leading to two different approaches:

1) The intrusive ROM resolves the governing equations by projecting them into a set of basis function of smaller dimension, leading to a system of ordinary differential equations for the coefficients. The projection-based methods have the advantage to retain some of the physics from the governing equations and to give rigorous error bounds and error estimation [START_REF] Veroy | Certified Real-Time Solution of the Parametrized Steady Incompressible Navier-Stokes Equations: Rigorous Reduced-Basis a Posteriori Error Bounds[END_REF]. However, both stability and accuracy issues can occur [START_REF] Xiao | Non-Intrusive Reduced Order Modelling of Fluid-Structure Interactions[END_REF]. By their intrusive nature, these methods modify also the source code of the high-fidelity model, leading to substantial modifications, if not impossible, when commercial software packages are used. Moreover, the reduced equations are solved on the whole domain and for all the conservatives variables even if the quantities of interest are evaluated on a subdomain and for a small number of variables. The high-Reynolds-number turbulent flows [START_REF] Protas | Optimal Nonlinear Eddy Viscosity in Galerkin Models of Turbulent Flows[END_REF][START_REF] Lorenzi | POD-Galerkin Method for Finite Volume Approximation of Navier-Stokes and RANS Equations[END_REF] give an example of challenging and active field of research for projection-based ROM. More information on intrusive ROM can be found in Ref. [START_REF] Benner | A Survey of Projection-Based Model Reduction Methods for Parametric Dynamical Systems[END_REF].

2) The second method, the nonintrusive data fitting ROM, does not need any knowledge about the high-fidelity model, considered as a black box, allowing to deal with very complex physics. Instead of manipulating the governing equations, the value of the coefficients are predicted by methods of multidimensional data fitting such as polynomial regression, radial basis function, or Gaussian process regression [START_REF] Forrester | Engineering Design via Surrogate Modelling: A Practical Guide[END_REF]. Nonintrusive ROM has been successfully applied in CFD, for instance, in aero-icing problems [START_REF] Fossati | Multiparameter Analysis of Aero-Icing Problems Using Proper Orthogonal Decomposition and Multidimensional Interpolation[END_REF], uncertainties quantification for urban flow [START_REF] Margheri | A Hybrid Anchored-ANOVA-POD/ Kriging Method for Uncertainty Quantification in Unsteady High-Fidelity CFD Simulations[END_REF], or steady aerodynamics [START_REF] Thanh | Proper Orthogonal Decomposition Extensions for Parametric Applications in Transonic Aerodynamics[END_REF][START_REF] Braconnier | Towards an Adaptive POD/SVD Surrogate Model for Aeronautic Design[END_REF][START_REF] Mifsud | A High-Fidelity Low-Cost Aerodynamic Model Using Proper Orthogonal Decomposition[END_REF]. One can note that hybrid methods mixing projection-based and nonintrusive approaches have been developed by solving an inverse problem, where the coefficients of the reduced-order governing equations are inferred using the output of the simulations [START_REF] Peherstorfer | Data-Driven Operator Inference for Nonintrusive Projection-Based Model Reduction[END_REF].

Here, only nonintrusive data fitting are considered and the highfidelity model is treated as a black box.

Despite the usefulness of the energy-ranked POD and its extensive use, some limitations have been observed [START_REF] Graham | Alternative Approaches to the Karhunen-Loeve Decomposition for Model Reduction and Data Analysis[END_REF]. Indeed, low-energy perturbations can be masked although they could be representative of a part of the system behavior. Problems with bifurcations can have typical characteristics, such as aerodynamic flows with varying inflow conditions leading to either subsonic or transonic regime. In these cases, the classical method computing the dominant modes in a single POD basis fails to produce accurate responses for predictions in highly nonlinear region and not directly in the neighborhood of the snapshots [START_REF] Lucia | Reduced Order Modeling of a Two-Dimensional Flow with Moving Shocks[END_REF][START_REF] Iuliano | Proper Orthogonal Decomposition, Surrogate Modelling and Evolutionary Optimization in Aerodynamic Design[END_REF]. The mix of the different physical regimes in the POD basis vectors can explain this problem. Indeed, small errors in the multivariate interpolation step can amplify POD modes associated with a physical regime that does not exist for the considered prediction. For this reason, approaches based on local reduced-order models have emerged in the literature by considering only restrictions to the total amount of snapshots [START_REF] Eftang | An "hp" Certified Reduced Basis Method for Parametrized Elliptic Partial Differential Equations[END_REF][START_REF] Amsallem | Nonlinear Model Order Reduction Based on Local Reduced-Order Bases[END_REF][START_REF] Washabaugh | Nonlinear Model Reduction for CFD Problems Using Local Reduced-Order Bases[END_REF][START_REF] Franz | Interpolation-Based Reduced-Order Modelling for Steady Transonic Flows via Manifold Learning[END_REF][START_REF] Kaiser | Cluster-Based Reduced-Order Modelling of a Mixing Layer[END_REF][START_REF] Liem | Surrogate Models and Mixtures of Experts in Aerodynamic Performance Prediction for Aircraft Mission Analysis[END_REF][START_REF] Zhan | Local Reduced Order Modeling and Iterative Sampling for Parametric Analyses of Aero-Icing Problems[END_REF]. This paper describes an original active local method, called "local decomposition method" (LDM), extending the classical reduced-order modeling method using POD and data fit method to particular steady problems with different physical regimes.

The method proposed here computes local subspaces of the stateparameter space by combining a physical-based sensor with machine learning tools. The physical-based sensor is a central element of the method to achieve proper separations of the physical regimes. Indeed, the conversion of the vector-valued output into a vector of physical-based features gives the possibility to cluster the snapshots into subsets with the same physical behavior. Thus the POD basis vectors are more representative of the physics. A shock sensor is used for the particular problems mixing subsonic and transonic conditions. It measures the nonlinearities and sharp gradients of the flow field. As the different phenomena are no more mixed in the POD basis, building a local reduced-order model on each of these subsets achieves a better consideration of the physical regimes. The clustering of the snapshots provides also a greater flexibility to the data fit model that can behave independently on each subgroup. As regard the prediction of untried sets of parameters, a supervised learning algorithm associates each region of the parameter space with a local reduced-order model and its respective subspace, allowing to map the input parameter space to the right physical regime. This last step is called input space decomposition.

Replacing the global POD basis with several local POD bases is a relatively recent development for parametric reduced-order modeling [START_REF] Amsallem | Nonlinear Model Order Reduction Based on Local Reduced-Order Bases[END_REF][START_REF] Washabaugh | Nonlinear Model Reduction for CFD Problems Using Local Reduced-Order Bases[END_REF][START_REF] Liem | Surrogate Models and Mixtures of Experts in Aerodynamic Performance Prediction for Aircraft Mission Analysis[END_REF][START_REF] Zhan | Local Reduced Order Modeling and Iterative Sampling for Parametric Analyses of Aero-Icing Problems[END_REF] and may seem counterintuitive. Indeed, the compressive capability of the POD may be weakened by increasing the number of bases, and the robustness of each data fit method can decrease with the reduction of training samples due to the repartition of these latter on the different models. On the contrary, the local models enable a clear separation of the phenomena improving the prediction of the surrogate model. Moreover, as underlined in [START_REF] Noack | From Snapshots to Modal Expansions-Bridging Low Residuals and Pure Frequencies[END_REF], the POD shows elliptic properties while aerodynamics is characterized by hyperbolic equations. Therefore, a local approach can demonstrate similarities with discontinuous Galerkin method to tackle hyperbolic problems. In this work, only Gaussian process regression (GPR) is investigated as a method of data fitting. Thus the terms "POD/data fitting" and "POD/ GPR" are used indifferently to refer to the classical nonintrusive ROM.

This paper aims to present a local reduced-order model built with machine learning tools and using a physical-based approach in order to address parameter-dependent problems with either subsonic or transonic regime. It is organized as follows: Sec. II gives an overview of the classical nonintrusive POD/data fitting approach. Then, Sec. III introduces the LDM with its underlying principles based on machine learning. Then, results from a one-dimensional analytical case and twodimensional transonic airfoil are presented in Sec. IV, demonstrating the capability of the LDM to deal with different physical regimes, including shock waves. Finally, Sec. V provides a summary and the conclusions.

First, some notations are introduced. One considers a real vectorvalued function f representing the high-fidelity model defined from R p to R d , where p is the number of parameters and d the dimension of the vector-valued quantity of interest. For example, a CFD code predicting the wall pressure field of an airfoil for different values of Mach number and angle of attack defines an input domain with p equal to 2 and d corresponding to the number of nodes representing the wall. Similarly, f represents the vector-valued surrogate model with the same domain of definition from R p to R d . The matrix of the training input parameters is noted χ t χ t 1 : : : χ t n T ∈ R n×p , where n is the number of training samples, and χ t i is the ith vector of the parameter set, which can be written with its components as χ t i χ 1 t i : : : χ p t i . In the same way, the matrix of the test samples, also referred to as untried input parameters or merely the predictions, is noted: χ p χ p 1 : : : χ p m T ∈ R m×p , with m the number of predictions. S i designates the vector of the ith snapshot such that S i fχ t i ∈ R d and S t defines the matrix of the training snapshots S t S 1 : : : S n T ∈ R n×d . The terms S t and f are used interchangeably to refer to the mean snapshot S f 1∕n P n k1 S k ∈ R d . The fluctuating part S 0 t of the training snapshots is defined by the snapshots matrix where the mean snapshot has been removed such that S 0 t S 1 -S : : : S n -S T . The first mode of data with nonvanishing mean is very close to the mean value of the snapshots for reasonable variations in the data. For this reason, the POD is performed on the fluctuating quantity S 0 t . The nonintrusive POD/data-fit reduced order modeling is made of three steps: the generation of the initial training samples, the POD dimension reduction, and the interpolation of the POD coefficients. All of them are described in the following parts. The whole process of the method is depicted in Fig. 1. One can note that the POD/GPR method is the common basis of the further developments presented in this paper and serves as a reference for the method assessment performed in Sec. IV.

A. Sampling Plans

The purpose of the surrogate model is to simulate the input/output behavior over the domain of variation of the parameters based on a limited number of high-fidelity evaluations. Correct predictions are required not only for the design points but also at all off-design conditions. For this reason, the limited number of snapshots should be optimally placed in order to build a model capturing the maximum amount of information about the physics over the parameter space.

An inappropriate repartition of the input parameters could lead to a surrogate model with large discrepancies. To the extent possible, the number of snapshots will be limited to the rule of thumb 10p studied by Loeppky et al. [START_REF] Loeppky | Choosing the Sample Size of a Computer Experiment: A Practical Guide[END_REF], where p is the number of parameters.

Contrary to the projection-based method, the nonintrusive POD/ GPR approach does not have access to prior information on the system given by the coefficients of the governing equations [START_REF] Paul-Dubois-Taine | An Adaptive and Efficient Greedy Procedure for the Optimal Training of Parametric Reduced-Order Models[END_REF]. The critical issue of the choice of the a priori training snapshots is called design of experiment (DOE). As explained previously, the main goal of the DOE is to generate well-distributed samples in the parameter space to give sufficient information to the learning process. DOE methods have been widely studied in the literature providing many techniques for experimental parametric studies and computer experiments. One can cite for example random and orthogonal array methods with Monte Carlo and Latin Hypercube Sampling [START_REF] Mckay | A Comparison of Three Methods for Selecting Values of Input Variables in the Analysis of Output from a Computer Code[END_REF], geometrical approaches such as centroidal Voronoi tessellations [START_REF] Du | Centroidal Voronoi Tessellations: Applications and Algorithms[END_REF], or low-discrepancy sampling techniques such as Halton, Sobol, or Faure sequences [START_REF] Kalos | Monte Carlo Methods[END_REF]. In the present work, low-discrepancy sequences have been adopted due to their iterative design. Indeed, the number of samples can be extended on purpose, and high-density regions can be easily defined while keeping the space-filling properties. Both properties are very interesting features for active learning. The other deterministic methods require a preset number of samples and cannot be extended without losing a part of their space-filling property.

B. Proper Orthogonal Decomposition

The POD is an efficient technique of dimension reduction based on spectral decomposition for high-dimensional, multivariate, and nonlinear data set. A wide range of applications can be found in the literature, such as human face characterization [START_REF] Kirby | Application of the Karhunen-Loeve Procedure for the Characterization of Human Faces[END_REF], data compression [START_REF] Andrews | Adaptive Data Compression[END_REF], and optimal control [START_REF] Lombardi | Low-Order Models: Optimal Sampling and Linearized Control Strategies[END_REF]. The POD term was first introduced in 1967 [START_REF] Lumley | The Structure of Inhomogeneous Turbulence[END_REF] to study dominant turbulent eddies, also called coherent structures. POD is also known as Karhunen-Loève decomposition, Hotelling analysis, or principal component analysis, in other fields of application. Among all the possible linear decompositions of the high-fidelity function, the POD method minimizes in a least square sense the residual of the projection of the high-fidelity model, yielding an optimal basis in terms of the representativeness of the data [START_REF] Cordier | Review of Some Fundamentals of Data Processing: Proper Orthogonal Decomposition[END_REF]. The least square problem is equivalent to a maximization problem [START_REF] Berkooz | The Proper Orthogonal Decomposition in the Analysis of Turbulent Flow[END_REF][START_REF] Cordier | Review of Some Fundamentals of Data Processing: Proper Orthogonal Decomposition[END_REF]. Introducing the canonical inner product (⋅, ⋅) on R d , the POD basis is the solution of the following formulation:

max ϕ 1 ; : : : ;ϕ n X n i1 X n j1 j S 0 i ; ϕ j j 2 subject toϕ i ; ϕ j δ i;j (1)
where ϕ i ∈ R d is the ith vector basis, S 0 i the fluctuating part of the ith snapshot, and δ i;j the Kronecker symbol satisfying δ i;j 1 for i j and δ i;j 0 otherwise. The matrix of the vector basis ϕ ∈ R n×d is introduced such that ϕ ϕ 1 : : :

ϕ n T .
The method of snapshots proposed by Sirovich [START_REF] Sirovich | Turbulence and the Dynamics of Coherent Structures Part I: Coherent Structures[END_REF] is employed to solve the maximization problem in Eq. ( 1) and leads to an eigenvalues problem:

1 n S 0 t S 0T t ϕλϕ T (2)
where λ ∈ R d×d is the diagonal matrix of the eigenvalues associated to the matrix of the eigenvectors ϕ. The eigenvalues problem can be solved either by an eigen decomposition or a singular values decomposition (SVD). The latter gives a better precision for the smaller eigenvalues and provides an iterative approach well-fitted for resampling. One can note that the self-adjoint operator property of S 0 t S 0T t ensures that the computed POD modes form a complete orthonormal set built as fϕ 1 ; : : : ; ϕ n g, on which the high-fidelity model is decomposed: where a k χ t i ∈ R is the reduced coordinate associated with the kth POD mode ϕ k ∈ R d . All the reduced coordinates are computed using the orthonormality property of the POD basis and are expressed as:

f χ t i f X n k1 a k χ t i ϕ k ; ∀ i ∈ 1; n (3)
a k χ t i S 0 i ; ϕ k (4)
The POD is optimal in terms of energy and provides an energyranked basis. Usually, the dimensional space of the snapshots is sparse in an L 2 sense. Therefore, only a small number of the most energetic POD modes can be retained in order to reduce the dimension of the system. The smallest eigenvalues are neglected, leading to a truncation of the basis. This heuristic criterion can be written more formally as a minimum ratio of the captured energy, which means finding the number M of kept basis vectors such that for a given amount of energy ratio ϵ:

P M k1 λ k P n j1 λ j > ϵ (5)
A classic ratio of energy present in the literature is 0.99 [START_REF] Cordier | Two Typical Applications for POD: Coherent Structures Education and Reduced Order Modelling[END_REF][START_REF] Sirovich | Turbulence and the Dynamics of Coherent Structures Part I: Coherent Structures[END_REF]. Once M has been set, the truncated linear combination of the eigenfunctions gives the approximation of the high-fidelity model:

f χ t i ≃ f X M k1 a k χ t i ϕ k ∀ i ∈ 1; n (6)

C. Interpolation with Gaussian Process Regression

The interpolation is the next step of the POD/GPR surrogate model. The reduced coordinates a k have been computed at a small number of training parameters, whereas the analysis of the highfidelity model for various inflow conditions requires a continuous evaluation over the input parameter space. Thus, the values of the reduced coordinates at untried parameter combinations are estimated with a data-fit method. The most popular methods in the surrogate modeling literature are formed of polynomial regression, radial basis function, or Gaussian process regression (GPR) [START_REF] Forrester | Engineering Design via Surrogate Modelling: A Practical Guide[END_REF]. A particular emphasis is given to this latter method. It has been employed in this paper due to its capability to deal with nonlinear problems, its high flexibility, and the provided error estimation of the predictor. GPR is also called "kriging" and has been first applied in geostatistics [START_REF] Krige | A Statistical Approach to Some Mine Valuation and Allied Problems on the Witwatersrand[END_REF]. The short overview of the GPR in this paper is introduced following the formalism of Rasmussen [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF] and is directly applied to the continuous representation of the reduced coordinates. They are assumed to follow a Gaussian process, which is outlined by a collection of random variables having a joint Gaussian distribution of mean μ and covariance matrix Σ. If the reduced coordinate a k follows a Gaussian process, it reads:

A k t ∼ N μ k ; Σ k (7)
with A k t a k χ t 1 : : : a k χ t n T defining the matrix of the kth reduced coordinate at training parameters, and A k p a k χ p 1 : : : a k χ p m T the predictions at the unknown combinations of parameters. In the interests of simplifying notation and analysis, the index k is removed and becomes implicit. The joint distribution of the reduced coordinates at training and unknown parameters is given by:

A t A p ∼ N μ t μ p ; Σ tt Σ tp Σ pt Σ pp (8)
with Σ tp the covariance matrix between χ t and χ p . The central issue of the GPR remains to be addressed, namely, how to determine the means and the covariance between the inputs. A classical stationarity assumption is that the correlation depends only on the magnitude of the Euclidean distance between the two input parameters but not on the values themselves, such that the i, jth element of the covariance matrix is given by:

Σ tp i;j σ 2 0 r kχ t i -χ p j k 2 ∀ i; j ∈ 1; n × 1; m (9)
with σ 2 0 the prior covariance corresponding to the level of uncertainty for predictions far from the training data and r the correlation function. The latter is usually monotonically decreasing with r0 1. A wide range of functions have been proposed to model the relation between the covariance and the input distance, such as radial basis function, Mattern, or periodic regression function [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF]. The anisotropic radial basis function, chosen for our problem, is one the most classical model of the regression functions due to its smoothness, infinite differentiability, and analytical derivability. Its expression introduces the hyperparameter θ k ∈ R and is given by:

r kχ t i -χ p j k 2 Y p k1 exp -kχ k t i -χ k p j k 2 2θ k ! ∀ i;j ∈ 1;n × 1;m (10)
The hyperparameter θ k defines the way the data are explained by the component k of the input parameters. Small values mean that the correlation is high between the inputs and the model is very sensitive to the dimension k. On the other hand, large values of θ k illustrate a model slowly varying with the data.

The final form of the predictor is derived using the conditional distribution of A u , given A t , also called a posteriori distribution that is still Gaussian and written as:

p A p jA t ∼ N μ p Σ pt Σ -1 tt A t -μ t ; Σ pp -Σ pt Σ -1 tt Σ tp (11)
The mean of the distribution gives the final value of the predictions at the untried set of parameters. Regarding the variance, it provides an estimate of the possible range taken by the prediction, which can also be seen as the mean-square error. The latter has the interesting feature not to be dependent on the value of the output but only to the input parameters. Both θ k and σ 2 0 remain to be determined in order to obtain the final prediction. They are computed during the training phase of the GPR, most of the time by a maximum likelihood estimation (MLE) approach or a leave-one-out method [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF]. Martin and Simpson [START_REF] Martin | Use of Kriging Models to Approximate Deterministic Computer Models[END_REF] have shown that the MLE works better than leave-one-out in general. For this reason, MLE is applied in the GPR/POD method, solving numerically the nonlinear maximization problem of the log-likelihood:

logpA t jθ - 1 2 A T t Σ -1 tt A t - 1 2 log jΣ tt j - n 2 log2π ( 12
)
where j ⋅ j denotes the determinant operator. The partial derivatives of the marginal likelihood with regard to the hyperparameters can be analytically derived. Thus, it is possible to use a gradient-based optimization algorithm in order to numerically find a local solution to the MLE problem. The Limited-memory Broyden-Fletcher-Goldfarb-Shanno Bounded (L-BFGS-B) algorithm [START_REF] Byrd | A Limited Memory Algorithm for Bound Constrained Optimization[END_REF] is employed in this paper to determine the hyperparameters. This popular quasi-Newton method handles simple bound constraints and is coupled with random restarts to avoid local maximum of bad quality. By assuming that the POD basis vectors are invariant with respect to the input parameters, the final surrogate model predicts the quantity of interest at the jth untried input parameter χ p j such that:

fχ p j f X M k1 e a k χ p j ϕ k ; ∀ j ∈ 1; m (13)
with e a k χ p j

h μ k p Σ k p j t Σ k-1 tt A k t -μ k t i j
the approximation of the weighting coefficient a k over the parameter space. The a priori mean μ k p is usually considered equal to zero as the training data have been standardized with zero mean. One can note that different versions of GPR or kriging can be used, such as Bayesian kriging [START_REF] Fossati | Evaluation of Aerodynamic Loads via Reduced-Order Methodology[END_REF]. Here, the Python library scikit-learn [START_REF] Pedregosa | Scikit-Learn: Machine Learning in Python[END_REF] is employed to generate the GPR models.

III. Local Decomposition Method

As explained in Sec. I, the LDM proposed in this paper extends the classical POD/GPR reduced-order modeling by employing a local approach, inspired by the mixture of experts [START_REF] Liem | Surrogate Models and Mixtures of Experts in Aerodynamic Performance Prediction for Aircraft Mission Analysis[END_REF] and dynamic local reduced-order modeling [START_REF] Amsallem | Nonlinear Model Order Reduction Based on Local Reduced-Order Bases[END_REF]. Instead of a unique global POD basis, several local bases are computed using machine learning tools yielding to more flexible behaviors bringing out a precise delimitation of the physical regimes. First, a shock sensor computes specific features for all the snapshots. Then, the latter are clustered into different subsets thanks to the shock features. Finally, the parameter space is divided into several domains according the clustering of the snapshots. The Fig. 2 sketches the process. One can note that a comparable approach has been used for aero-icing certification [START_REF] Zhan | Local Reduced Order Modeling and Iterative Sampling for Parametric Analyses of Aero-Icing Problems[END_REF]. The specificity of the presented method includes, besides the introduction of a feature extraction with a shock sensor, a novel resampling strategy and the application to an aerodynamics case. The active resampling is carried out by identifying the subspaces with the highest entropy. Extra snapshots are added in these specific subspaces with the objective to minimize the redundancy of the sampling, thus increasing the accuracy of the surrogate model. The Fig. 3 illustrates the whole method.

A. Coupling Machine Learning Tools with a Physical Sensor

Let us introduce basic machine learning vocabulary. Learning problems can be divided into two distinctive categories: supervised and unsupervised. Here, the machine learning library scikit-learn [START_REF] Pedregosa | Scikit-Learn: Machine Learning in Python[END_REF] is employed in the in-house JPOD code to perform both supervised and unsupervised learning. In the context of supervised learning, some input variables have an influence on one or more outputs. This set of inputs and outputs forms a learning base, and the supervised learning simulates the input/output behavior using the learning base. The final goal is to predict the values of the outputs for untried inputs. The nature of the output subdivides the supervised learning into two subcategories: the classification, dealing with categorical input variables, and the regression, which is applied on real and continuous input variables. The GPR is an example of regression. As regard the unsupervised learning, the training set consists only of the input vectors without any corresponding outputs. Thus the purpose of the unsupervised learning is to identify underlying structures hidden in the parameter space but the accuracy of the algorithm cannot be defined by any objective function.

Physical-Based Shock Sensor to Detect Flow Regimes

The choice of the quantity characterizing the physical regimes, on which the clustering is performed, is a question of central importance impacting the quality of the classification. Usually, the unsupervised learning clusters directly the quantity of interest into groups with patterns of small differences [START_REF] Amsallem | Nonlinear Model Order Reduction Based on Local Reduced-Order Bases[END_REF][START_REF] Washabaugh | Nonlinear Model Reduction for CFD Problems Using Local Reduced-Order Bases[END_REF][START_REF] Zhan | Local Reduced Order Modeling and Iterative Sampling for Parametric Analyses of Aero-Icing Problems[END_REF]. However, the aim of the clustering in this paper is the physical regime separation, and the previous approach can lead to classification error. Indeed, two fields of the quantity of interest can have large differences even though they belong to the same physical regime. A classical method fails to separate them accurately. For this reason, this section proposes another method to perform the clustering. Based on a physical approach, a mathematical transformation converts the quantity of interest into a sensor of the physical regime. The main goal is to sharply quantify the physical regime to ease the clustering of the snapshots.

The application of this paper involves external aerodynamics with subsonic and transonic regimes, characterized by shock waves.

A straightforward idea is to consider a shock sensor that is able to detect large changes in the variation of the quantity of interest, such as Jameson's shock sensor [START_REF] Jameson | Numerical Solutions of the Euler Equations by Finite Volume Methods Using Runge-Kutta Time-Stepping Schemes[END_REF]. The latter has been developed as a heuristic method for nonoscillatory shock capturing numerical schemes. It is related to the second-order derivative of the pressure. However, the quantity of interest is not limited to pressure signals and a more general expression is introduced as:

ν i js i-1 -2s i s i1 j ϵ 0 js i-1 j 2js i j js i1 j ; ∀ i ∈ 2; d -1 ( 14
)
where ν i is the generalized sensor, ϵ 0 is a constant avoiding division by 0, and s i-1 , s i , and s i1 are the quantity of interest at, respectively, node i -1, i, and i 1. The dimension of the generalized sensor can be reduced by POD, such as the kth POD basis vector is associated to the reduced coordinates b k . The latter can be interpreted as the representative quantity of the physical regime for a given snapshot. They are grouped in the vector B i defining the matrix B B 1 : : : B n T

Clustering of the Shock Sensors by Gaussian Mixture Model

The problem of identifying the inherent groupings in the input data refers to unsupervised classification, specifically clustering. K-means and Gaussian mixture model (GMM) are two well-known examples of algorithms classically employed to achieve clustering. This section puts a particular emphasis on GMM algorithm [START_REF] Mclachlan | Mixture Models. Inference and Applications to Clustering[END_REF], where the main features are described.

Let us assume that the set fB 1 ; : : : ; B n g, characterizing the physical regimes, comes from q clusters C 1 ; : : : ; C q . Each cluster C k follows a probability distribution of parameter θ k and proportion w k , regrouped for all the clusters in a mixture parameter Φ w 1 : : : w q θ 1 : : : θ q . The GMM algorithm consists in modeling B with a mixture distribution of multivariate normal distributions g. Each one is associated with the cluster C k such that the probability density function of B i is given by:

pB i jΦ X q k1 w k gB i jθ k ; ∀ i ∈ 1; n (15)
, where the mixture weights w k represent the probability that the observation comes from the kth Gaussian distribution and θ k gives the mean and the covariance of the multivariate normal distribution g. These mixture parameters regrouped in Φ are estimated iteratively using an Expectation Maximization algorithm (EM algorithm) [START_REF] Dempster | Maximum Likelihood from Incomplete Data via the EM Algorithm[END_REF]. The expectation of the posterior probability γ k of belonging to the cluster k can be expressed with Bayes' theorem (E-step):

γ k B i pC k jB i pB i jC k pC k PB i w k gB i jθ k P q l1 w l gB i jθ l ; ∀ k ∈ 1; q; ∀ i ∈ 1; n (16)
Then, the mixture parameters can be re-estimated (M-step):

μ k P N i1 γ k B i ⋅ B i P N i1 γ k B i ; ∀ k ∈ 1; q (17)
Σ k P N i1 γ k B i ⋅ B i -μ k ⋅ B i -μ k T P N i1 γ k B i ; ∀ k ∈ 1; q (18) w k 1 n X N i1 γ k B i ; ∀ k ∈ 1; q (19)
These two steps are iterated until convergence of the log likelihood. The cluster of each quantity B i can be determined using the previous probability expression. The training set is built by applying a hard splitting such that: fχ t 1 ;y t 1 ;:::;χ t n ;y t n g with y t i C k j max k∈0;q pC k jB i ; ∀i∈1;n

with y t i the target variable of the cluster. This training set is used to train a supervised learning algorithm, described in the following section, in order to link the input parameter of each quantity B i with a class of physical regime.

Input Space Decomposition Using Gaussian Process Classification

The decomposition of the input space into subspaces where a single physical regime drives the flow can be interpreted as a supervised classification problem. Indeed, the training set is provided by the clustering phase and trains an algorithm assigning the kth untried input parameters χ p j (j ∈ 1; m) to the categorical variable y p j , which can take the different values of the q classes C 1 ; : : : ; C q . The Gaussian process classification (GPC) is a classical method to deal with classification. The principal steps of the method for two-class problem are outlined in this section. The interested readers can refer to Rasmussen and Williams [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF] and Bishop [START_REF] Bishop | Pattern Recognition and Machine Learning (Information Science and Statistics)[END_REF] for the straightforward generalization to K classes.

As the nature of the classification outputs is discrete, it clearly differs from the regression problems outlined previously with GPR. The main idea is to transform the output of a Gaussian process defined on the real axis into a probability lying in the interval [0, 1] using a nonlinear activation function. A latent function l defined on the input parameter is introduced and we denote the latent vector by l such that: l l t 1 ; : : : ; l t n with l t i lχ t i and l p j lχ p j ;

∀ i; j ∈ 1; n × 1; m (21)
This function aims to provide a more convenient and tractable formulation of the model [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF] and will be removed by integration. A Gaussian process prior with a zero mean and a covariance matrix Σ l is placed on the joint distribution of the latent function l: l

l p j ∼ N 0; Σ l (22)
As regard the two-class problem with C 0 and C 1 , the probabilistic prediction is directly computed by p y p j C 0 jχ t ; y t ; χ p j because p y p j C 1 jχ t ; y t ; χ p j is given by 1p y p j C 0 jχ t ; y t ; χ p j . The conditioning on the input variables is intentionally let implicit. The probabilistic prediction is given by: p y p j C 0 jy t Z p y p j C 0 jl p j p l p j jy t dl p j [START_REF] Amsallem | Nonlinear Model Order Reduction Based on Local Reduced-Order Bases[END_REF] where py p j C 0 jl p j σl p j , with σ the nonlinear activation function defined by the sigmoid function:

σx 1 1 e -x (24)
The integral expressed in Eq. ( 23) is analytically intractable due to the non-Gaussian likelihood of p l p j jy t [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF]. The expansion of the latter with the sum rule, product rule, and Baye's theorem gives:

p l p j jy t Z pljy t p l p j jl dl ( 25
)
where p l p j jl is Gaussian. The non-Gaussian probability pljy t requires specific approximations, such as variational inference, expectation propagation, or Laplace approximation [START_REF] Bishop | Pattern Recognition and Machine Learning (Information Science and Statistics)[END_REF]. Finally, the hyperparameters of the covariance matrix Σ l need to be determined, for example, with the maximization of the log-likelihood, which also required to use the Laplace approximation due to non-Gaussian terms.

B. An Entropy-Based Active Resampling

Several methods have been coupled with surrogate models to generate an active reduced-order model: local methods and subspace methods form the two main different approaches. Local methods look for particular points that could improve the accuracy of the model. One can cite, for example, leave-one-out cross-validation [START_REF] Fossati | Evaluation of Aerodynamic Loads via Reduced-Order Methodology[END_REF] testing the sensibility of the surrogate model to each training sample. The more critical sample for the surrogate defines a neighborhood in which an extra snapshot is added. Another example of a local method, the maximum mean-squared error, uses the posterior estimation of variance from the GPR and adds the sample with the maximum value to the training set. Multifidelity can also be used to perform active infill sampling [START_REF] Benamara | Adaptive Infill Sampling Criterion for Multi-Fidelity Optimization Based on Gappy-POD[END_REF]. Another type of method based on a subspace approach identifies particular low-dimensional structures in the input parameters where the quantity of interest shows a significant variability. One can cite, for example, active subspaces [START_REF] Constantine | Active Subspace Methods in Theory and Practice: Applications to Kriging Surfaces[END_REF] or sensibility analysis [START_REF] Saltelli | Global Sensitivity Analysis. The Primer[END_REF].

The original strategy proposed in this paper is based on a subspace approach and aims at taking advantage of the input space separation into several subsets. Indeed, the input space decomposition has provided subgroups of smaller dimensions among which some subspaces of interest can be selected to perform the resampling. One proposes to use a criterion based on the compressibility of the information, coming from an analysis of the POD eigenvalues, in order to identify these relevant structures. As explained in the last section, the POD eigenvalue represents the relative information contained by the modes. The global entropy H measures the redundancy of this information and is introduced as [START_REF] Cordier | Review of Some Fundamentals of Data Processing: Proper Orthogonal Decomposition[END_REF]:

H - 1 logn X n k1 p k logp k with p k λ k P n i1 λ i ( 26
)
If the entropy goes to zero, there is only one nonzero singular value. The data are compressed into a unique mode. On the other hand, the entropy is equal to one if all the information is distributed among the modes, meaning that no compression is possible. Between these two extreme values, the entropy increases with the number of fundamental modes. The active resampling of the LDM assumes that the entropy and the structures of the system with discontinuities or high gradients are directly correlated. Thus, the probability to find new modes with a nonnegligible amount of energy is expected to be greater for the cluster with the highest value of entropy than for any other cluster.

C. Recombination in a Global Model by Hard-Splitting

The final recombination step consists of assembling the local reduced-order models in a single composite global model. Starting from the q clusters, a simple weighted sum is calculated using a "hard" split:

fχ p j X q i1 1 C i χ p j fi χ p j ; ∀ j ∈ 1; m ( 27
)
where

1 C j χ 8 < : 1 if j argmax i∈1;q
Pχ ∈ C i 0 else [START_REF] Zhan | Local Reduced Order Modeling and Iterative Sampling for Parametric Analyses of Aero-Icing Problems[END_REF] and fi refers to the classical POD/GPR model built on the ith cluster. This sum provides a continuous but not differentiable prediction of all the input space, leading to a global model. The differentiable predictions require to use differentiable weighting functions, which is not the case for 1. One can cite, for example, soft clustering [START_REF] Bettebghor | Surrogate Modeling Approximation Using a Mixture of Experts Based on EM Joint Estimation[END_REF], substituting 1 C j χ directly by Pχ ∈ C i . However, it mixes several physical regimes, leading potentially to unphysical predictions, and can amplify the extrapolation of the reduced coordinates near the boundary decision. For these reasons, a hard-split approach has been selected. The decision boundary in the input space parameters poses another problem. The classification is very prone to errors in this region. Indeed, the localization of the decision boundaries is subject to local variations, such as the choice of the supervised learning method (model-based, local methods, etc.) or the location of the training samples. Moreover, the reduced coordinates can be in extrapolation in this region. To overcome these problems, the classical model is employed for the predicted points near the interface.

IV. Numerical Results

A. Error Measurements

The classical method and the proposed LDM are assessed in this part with the RAE2822 airfoil [START_REF] Iuliano | Proper Orthogonal Decomposition, Surrogate Modelling and Evolutionary Optimization in Aerodynamic Design[END_REF][START_REF] Benamara | Adaptive Infill Sampling Criterion for Multi-Fidelity Optimization Based on Gappy-POD[END_REF]. It may be noted that an additional study has been performed on the Burgers' problem. The results are shown in the Supplemental Material burgers_am.pdf. The flow around the two-dimensional transonic airfoil is computed with a Navier-Stokes solver, involving a turbulence model. It is a challenging application in terms of surrogate modeling with high discontinuities due to the appearance of shocks. A three-dimensional input space is considered.

Particular attention is paid to the process of the input space decomposition and to the comparison of both classical approach and LDM in terms of accuracy. Several quantities are introduced to measure the accuracy. The so-called predictivity coefficient Q 2 gives the ratio of the output variance, which is explained by the metamodel. It can be interpreted as the classical coefficient of determination of the linear regression applied to a test sample [START_REF] Iooss | Numerical Studies of the Metamodel Fitting and Validation Processes[END_REF]. The more the value is close to 1, the higher variance is explained by the model. It is expressed as:

Q i 2 1 - P m j1 h f i χ p j -fi χ p j i 2 P m j1 h f i -f i χ p j i 2 (29)
The root mean square error (RMSE) and the normalized root mean square error (NRMSE) are also introduced:

RMSE i 1 m X m j1 h f i χ p j -fi χ p j i 2 v u u t (30)
NRMSE i RMSE i f max -f min (31)
where f max and f min refer, respectively, to the scalar maximum and minimum values taken by the function to predict f for all the snapshots and over the mesh:

f max max kf χ p 1 k ∞ ; : : : ; kf χ p m k ∞ f min min kfχ p 1 k ∞ ; : : : ; kf χ p m k ∞ (32)
One can note that the Q i 2 and NRMSE i are computed at a given index i of the spatial domain Γ. The global measures are provided by averaging the quantities over this domain. The global quantities are referred to as hQ 2 i Γ and hNRMSEi Γ , with h⋅i Γ the spatial average operator. However, in order to provide also a statistical error analysis, the averaged normalized error E i is introduced. It corresponds to the absolute error between the exact value and the prediction, normalized by the range of variation, at snapshot level j:

E j hjfχ p j -fχ p j ji Γ f max -f min ; ∀ j ∈ 1; m (33)
The statistical distribution of E i is presented with a box plot formalism. A box plot groups the data through different quantiles: the bottom and the top of the box represent, respectively, the value of the first and third quartiles, whereas the horizontal line inside the box is the median (second quartile) and the diamond the mean. The vertical lines indicate the data between the 5th percentile and the 1st quartile and between the 3rd quartiles and the 95th percentile. Finally the outliers are plotted as dots.

Besides the direct assessment of the surrogate model using a wide validation test, an estimation of the error can be provided only with the training set. Cross-validation strategies, such as the leave-one-out error, or the variance given by the Gaussian processes can estimate the a priori quality of each prediction. An example of application of these methods can be found in [START_REF] Dupuis | Aerodynamic Data Predictions for Transonic Flows via a Machine-Learning-Based Surrogate Model[END_REF].

B. Two-Dimensional RAE2822 Transonic Airfoil
The viscous and turbulent flow around an RAE2822 airfoil has been widely studied in the literature both numerically and experimentally [START_REF] Cook | Aerofoil RAE 2822: Pressure Distributions, and Boundary Layer and Wake Measurements[END_REF][START_REF] Haase | An European Initiative on Validation of CFD Codes: Results of the EC/BRITE-EURAM Project EUROVAL[END_REF]. The feature of prime interest of this test case is that the inflow conditions govern the flow regime, leading to the appearance of shock waves on the suction side. The detection and the clear separation of these regimes represent the main challenge for the model, demonstrating the capability of the LDM to deal with real and complex physics compared with the classical method.

Computational Configuration

The high-fidelity computations are carried out using the cellcentered finite-volume solver elsA-ONERA [START_REF] Cambier | The Onera elsA CFD Software: Input from Research and Feedback from Industry[END_REF]. It has been developed at ONERA and solves the compressible Reynoldsaveraged Navier-Stokes (RANS) equations on structured grids. From the numerical point of view, the classical second-order central scheme of Jameson et al. [START_REF] Jameson | Numerical Solutions of the Euler Equations by Finite Volume Methods Using Runge-Kutta Time-Stepping Schemes[END_REF] is used for the space discretization. The time integration is performed with the backward Euler implicit scheme: the algebraic system is linearized with the LU-SSOR implicit method [START_REF] Yoon | Lower-Upper Symmetric-Gauss-Seidel Method for the Euler and Navier-Stokes Equations[END_REF]. The turbulence modeling is ensured by the model of Spalart and Allmaras. A 2D mesh containing 23,010 points is used, as illustrated in the Fig. 4a. This test case has been successfully validated on a well-known regime flow [START_REF] Cook | Aerofoil RAE 2822: Pressure Distributions, and Boundary Layer and Wake Measurements[END_REF] (Fig. 4b). The chord of the airfoil is written C, X refers to the horizontal coordinate, and Y refers to the vertical coordinate. Three freestream parameters are considered as varying for this application: the flight speed, the angle of attack α, and the altitude h. Their variations are resumed in Table 1. These dimensional parameters are nondimensionalized and are expressed, respectively, as the Mach number M, angle of attack α (no change), and Reynolds number. One can note that the latter is impacted by the altitude variations but ensured to be maintained in a given interval in order to have a sufficiently resolved boundary layer. The atmosphere is modeled by the International Standard Atmosphere of the International Civil Aviation Organization [START_REF]Standard Atmosphere First Edition[END_REF]. It assumes that the air is a perfect gas and that the atmosphere can be divided into layers with a linear distribution of temperature against the altitude. The temperature T and the density ρ can be directly expressed in function of the altitude:

T T 0 -Lh (34)
ρ p 0 1 -Lh∕T 0 g∕rL rT 0 -Lh (35)
with p 0 and T 0 the pressure and temperature at sea level, L the temperature lapse rate, r the specific gas constant of air, and g the acceleration due to the gravity. The quantities of interest of the simulation are the pressure coefficient C p and the friction coefficient C f on the suction side defined by

C p p -p ∞ 1∕2ρ ∞ U 2 ∞ (36) C f τ w 1∕2ρ ∞ U 2 ∞ ( 37
)
where p is the static pressure; τ w is the wall shear stress; and p ∞ , ρ ∞ , and U ∞ are, respectively, the static pressure, the density, and the velocity in the freestream.

As regard the DOE, 30 samples of an Halton sequence form the training set, following the rule of thumb 10d [START_REF] Loeppky | Choosing the Sample Size of a Computer Experiment: A Practical Guide[END_REF]. The sampling of the LDM is divided into two parts. An initial DOE mixing subsonic and transonic snapshots explores uniformly the parameter space with an Halton sequence of 15 samples. The last 15 samples follow the resampling process described in the previous section. The flow is assumed to be driven by two different flow regimes. Therefore, the number of clusters is set to two for the clustering step. A test set has been built from 300 snapshots of a Sobol sequence in order to assess the LDM.

Analysis of the Surrogate Model Building Process

The two training sets exhibit interesting differences. The Halton approach explores uniformly the full input space, whereas the resampling process focuses on a low-dimensional high Mach number region, as depicted in Fig. 5 for the C p . Indeed, the clustering phase automatically identifies the subsonic and the transonic snapshots thanks to the shock sensor. The supervised algorithm decomposes the input space parameter, allowing to determine the separation of the two physical regimes in the input parameter space, as illustrated in the Fig. 6 with both the training and testing sets. These two clusters can be interpreted as the subsonic and the transonic regions. It can be observed that the boundary is mainly influenced by the Mach number but also slightly by the angle of attack. Thus, the resampling process has increased the density of samples in the transonic regime, improving the accuracy of the model where the predictions are more challenging. One can note that the boundary region is well defined by a thin region of probability between 40 and 60%.

The model can also be analyzed from a dimension reduction point of view with Table 2. For a given POD energy ratio of 0.99%, the classical method reduces the dimension of the 30 snapshots with 10 modes and an entropy of 0.37, whereas the LDM identifies clearly a POD basis with a low entropy and a POD basis with a large one. The highest nonlinear cluster shows an entropy value of 0.63, 70% bigger than the linear cluster and the classical method. It means that the data of the subsonic POD basis and of the classical method can be highly compressed. Conversely, the transonic POD basis requires 9 modes for 13 snapshots. Thus, the active resampling provides extra nonredundant information to the LDM compared with the classical method.

The shape of the modes can also help to understand the behavior of the different models. As regard the classical method, only discontinuous features emerge clearly from the global POD as illustrated in Fig. 7a, where POD modes shown only dominant discontinuities and nomoving shocks. Therefore all the reduced coordinates associated with subsonic snapshots must exactly cancel out the discontinuities of the modes. For this reason, the prediction of the reduced coordinates for snapshots in the subsonic region can be sensitive to interpolation errors leading to the appearance of "residual" shocks. Thus, the clear separation of the regimes in the POD domain represents a major asset inherent in the LDM. Figure 7b shows the first three modes of the two POD bases. High nonlinearities arises for the transonic regime, whereas the subsonic region highlights similarities with the modes of the classical method. The only difference is that the subsonic modes are smoother.

These statements are confirmed by looking at C p profiles, comparing classical and LDM methods. Eight predictions have been computed for an illustrative purpose and are grouped by physical regime (S for subsonic and T for transonic). Their coordinates in the parameter space are summarized in the Table 3 and their repartition is illustrated in Fig. 5. One can observe that for the subsonic regime in Fig. 8, the classical method induces residual shocks, certainly due to errors of prediction on the reduced coordinates coupled with sharp POD modes not adapted to subsonic flows. On the contrary, the residual shocks are filtered by the LDM as lower discontinuities are present in the training snapshots building the POD basis. It leads to final predictions less sensitive to errors on the reduced coordinates. As regard C p profiles shown in the transonic regimes in Fig. 9, the LDM shows improved accuracy. In particular, the modeling of the shock waves gives an insight into the behavior of the two models. Significant discrepancies in terms of shock displacement and shock amplitude occur for the classical POD/GPR method. Indeed, the combination of angle of attack and Mach number influences directly the localization of the steady shock wave and its amplitude. However, this behavior can be accurately caught only if the training set contains a sufficient amount of snapshots with shocks, which is not the case for the classical POD/GPR methods. The same trend is observed for the C f profiles in Figs. 10 and 11.

Accuracy of the Model

In this section, a more detailed look is given to the analysis of the model accuracy. Figure 12 displays the comparison of the normalized error in terms of C p and C f for both methods between the predictions and the test set. The results are presented with a box plot formalism and three different phases have been considered:

1) The full domain contains all the samples of the testing set.

2) The subsonic regime is only composed of the testing samples identified as subsonic, in blue dots in Fig. 6b.

3) The transonic regime encompasses the other snapshots, shown as red dots in Fig. 6b.

A significant improvement in the accuracy is induced by the LDM for the C p . As regard the full domain, the normalized error decreases dramatically for all the statistical characteristics of the box plot. In particular, the extreme value of the LDM reaches the same level as the 95% error of the classical method, illustrating a large reduction of the model variability. The box plots for subsonic and transonic regimes provide a closer look at the repartition of the error. It clearly appears that the LDM improves the predictions at transonic regime, explained by the increase of the sample density. On the other hand, the subsonic regime is very slightly impacted compared with the classical method, although samples have been removed. Table 4 provides a more global view of the error with the spatial average of the Q 2 and NRMSE. Whatever the quantity measuring the error, the same trend is observed for the two quantities of interest: the LDM significantly improves the accuracy of the predictions compared with the classical method. There are, however, legitimate doubts as to the value of the hQ 2 i Γ for the C f , which is below zero for the transonic regime. Several explications can be given. First of all, the hQ 2 i Γ is spatially averaged and the value of C f can be very close to zero after the shock. Thus, some values of Q 2 fall far below zero in this region, impacting directly the average of the predictivity coefficient. Then, the skin friction coefficient is also more challenging to predict due to its higher dependency on the altitude and higher nonlinearity.

Finally, the evolution of the NRMSE with the number of training samples is shown in Fig. 13 for both the pressure and skin friction coefficients. One may note that the classical method with 50 training samples has broadly the same error than the LDM strategy with 30 samples. In the end, this observation could be seen as a saving of the computational budget of more than 60%. Moreover, the error of the model seems to continuously decrease when the number of training samples increases.

V. Conclusions

The aim of this paper was to build a nonintrusive reduced order model capable of coping with aerodynamics problems at multiple operating conditions. These latter involve hyperbolic equations and different physical regimes leading to discontinuous features of solutions. The original local decomposition method (LDM) strategy alleviates the issues encountered by the classical method using global proper orthogonal decomposition (POD) and multiparameter interpolation in order to solve the regression problem. The LDM consists of building a local reduced order model for each physical regime. This approach has been developed to be more suited to hyperbolic properties than the global POD expansion. The decomposition of the problem into several local reduced order models is achieved using machine learning methods. Two major steps are associated with the LDM strategy: the use of a Jameson's shock sensor enhancing the physical regime recognition in order to ease the clustering of the problem and an active sampling adding automatically extra information to the subspaces with the highest discontinuous structures. The LDM has been assessed on an analytical moving shock problem and the simulation of a turbulent flow around the transonic RAE2822 airfoil. The results reveal a significant improvement of the model accuracy, especially in the regions of high gradients and discontinuities.

Further work is needed to increase the efficiency of the strategy. First, the resampling technique has to be improved. The snapshots added to a specific subspace come from the continuation of a low-discrepancy sequence. For example, a method minimizing the variance of the Gaussian process regression could be considered. Moreover, the new snapshots are computed sequentially as the probability of belonging to each subset is updated at each iteration. A new process of parallelization should be devised to take advantage of high-performance computing during the iterative sampling step by combining multiple simultaneous jobs into large ensembles. Second, the extrapolation at the interface of the input space parameter remains an open issue. Further applications of the LDM strategy are expected on three-dimensional flows for industrial aerodynamic applications. Especially, it is planned to evaluate this method in a mission analysis context, taking advantage of the results obtained for the transonic regions. framework (CGI, ANR) and the project industrial members for their support, financial backing, and/or own knowledge: Airbus, Airbus Group Innovations, SOGETI High Tech, Altran Technologies, and CERFACS. The studies presented in this article is making use of the elsA-ONERA software, whose the co-owners are Airbus, Safran, and ONERA.

  the reduced coordinates a k = kth reduced coordinate B = matrix of the reduced coordinates of the sensor b k = kth reduced coordinate of the sensor C = chord length C k = kth cluster C f = friction coefficient C p = pressure coefficient d = dimension of the quantity of interest E = averaged normalized error f = high-fidelity model g = acceleration due to the gravity or normal distribution H

Fig. 1

 1 Fig.1Flowchart of the POD/GPR method.

Fig. 2

 2 Fig. 2 Sketch of the LDM.

Fig. 3

 3 Fig. 3 Flowchart of the LDM.

2 .

 2 Input Parameter Space, Quantities of Interest, and Settings of the Surrogate Model

Fig. 4

 4 Fig.4Flow around the RAE2822 airfoil.

Fig. 5

 5 Fig. 5 Thirty samples from a classical Halton sequence (left) and resampling technique (right). The eight illustrative predictions are also identified (with blue squares and red triangles).

Fig. 6

 6 Fig.6Input space decomposition. Each color corresponds to a cluster: the blue one is the subsonic regime and the red one is the transonic regime.

  Fig. 7 POD modes of the RAE2822 simulations.

Fig. 8 CFig. 9 CFig. 10 C

 8910 Fig. 9 C p profiles in the transonic regime.

Fig. 12

 12 Fig. 12 Accuracy of the surrogate model in terms of C p and C f .

- 2 a

 2 ) Pressure coefficient b) Friction coefficient Fig. 13 Evolution of the NRMSE with the number of training samples.

Table 1

 1 Freestream conditions

	Freestream variable	Amplitude of variation
	Flight speed, m∕s	88.5-269
	Angle of attack, deg	0.5-3.0
	Altitude, m	1000-11,000

Table 2

 2 Required modes of the different methods for the RAE2822

		Number of	Number of
	Method	snapshots	modes	Entropy
	Classical	30	7	0.37
	LDM (transonic regime)	13	9	0.63
	LDM (subsonic regime)	17	6	0.36

Table 3

 3 Coordinates in the parameter space of the eight illustrative predictions

	Predictions Mach number Angle of attack, deg Altitude, m
	S1	0.675	1.125	8500
	S2	0.519	1.906	1625
	S3	0.722	1.047	1937
	S4	0.312	1.496	5453
	T1	0.796	2.980	6391
	T2	0.743	2.267	10,102
	T3	0.772	1.589	3168
	T4	0.794	2.170

Table 4

 4 Summary of the results in terms of Q 2 and NRMSE for the RAE2822RegimeMethodhQ 2 i Γ of C p hNRMSEi Γ of C p hQ 2 i Γ of C f hNRMSEi Γ of C f

	Both	Classical	0.897	1.93 × 10 -2	0.837	1.15 × 10 -2
		LDM	0.989	7.48 × 10 -3	0.908	9.53 × 10 -3
	Subsonic	Classical	0.984	7.46 × 10 -3	0.912	6.61 × 10 -3
		LDM	0.990	6.31 × 10 -3	0.957	6.72 × 10 -3
	Transonic	Classical	0.0974	4.64 × 10 -2	-0.303	2.54 × 10 -2
		LDM	0.901	9.79 × 10 -3	-0.0223	1.72 × 10

Acknowledgments

This work is part of the MDA-MDO project of the French Institute of Technology IRT Saint Exupery. We wish to acknowledge the PIA