An adaptive sampling algorithm for reduced‐order models using Isomap

The use of reduced‐order models (ROMs) in physics‐based modeling and simulation is a popular tool for drastically lowering the computational cost associated with high‐fidelity simulations. ROMs use training data from a set of computed high‐fidelity simulations with different design parameters that control physical and geometric properties of the full‐order model. The quality of the training data dictates the performance of the ROM, making the choice of training design parameters important. A widely used method for generating training parameters for ROMs is Latin hypercube sampling (LHS), a statistical method that aims to maximize the distance and minimize the correlation amongst produced samples. However, LHS fails to account for the physics of the full‐order model and can lead to an over‐representation of some physical regimes in the training data while neglecting others. In this work, we present a computationally efficient adaptive sampling method for ROMs using Isomap, a versatile algorithm for nonlinear dimensionality reduction. Using an initial number of samples, the adaptive sampling algorithm iteratively generates samples based on a low‐dimensional manifold of the training data. When applied to two external aerodynamics problems, the proposed adaptive sampling algorithm offers significantly improved performance in predicting physical fields and coefficients of lift and drag for a given computational budget for both non‐intrusive and projection‐based ROMs.

conditions, geometry of the computational domain, or physical properties.In many-query processes such as design optimization, a large number of simulations must be run at different designs, and the accuracy or fidelity of these simulations must be high.High-fidelity simulation is often very computationally intensive, requiring large amounts of computational resources and time.This cost can become a bottleneck in the design optimization process.
To mitigate this cost, using reduced-order models (ROMs) is a common approach.Using a small amount of training data from a set of computed high-fidelity simulations, ROMs create a low-dimensional surrogate model that can be evaluted in real-time and can provide accurate approximations of the full-order model at unseen design parameters.ROMs drastically lower the number of degrees of freedom of the full-order model using a low-dimensional embedding that can be mapped back to the full-order state.There are many methods to obtain the low-dimensional embedding, and the most popular one is the proper orthogonal decomposition (POD). 1,2Using the singular value decomposition, (SVD), a low-rank trial subspace consisting of a number of linearly independent basis vectors is obtained.Using a set of expansion coefficients, a linear combination of these basis vectors can be used to approximate full-order states within the solution space.ROMs consist of two stages: a computationally intensive offline stage where the training data from a full-order model is simulated and the ROM is trained, and an inexpensive online stage where the surrogate model is evaluated to approximate solutions at unseen design parameters.
ROMs can be either purely data-driven or can incorporate the physics of the governing equations to develop a low-dimensional version of the full-order model.Non-intrusive ROMs use the first approach by creating a regression model between the design parameters and the expansion coefficients to make approximations at unseen design parameters.Projection-based ROMs project the physics of the governing equations into the low-rank trial subspace to solve a low-dimensional version of the full-order model.As a result, projection-based ROMs are more computationally intensive than non-intrusive ROMs but can provide better accuracy and more physically consistent results, although they often suffer from stability issues.
Like all data-driven models, the performance of ROMs is highly dependent on the quality of the training data.Since ROMs typically use a very small amount of data, this further increases the sensitivity of model performance to data quality.To generate a set of design parameters  within a parameter space , a common approach is to use Latin hypercube sampling (LHS). 3LHS is a statistical method for generating samples from a multi-dimensional parameter space by aiming to maximize the distance and minimize the correlation amongst samples produced by using uniform sampling.LHS is useful for efficiently filling the parameter space and reducing the likelihood that points are clustered close to each other.However, LHS does not account for the physics of the full-order model when producing design parameters for the training set.It is often the case that the physical characteristics of the full-order model will vary non-linearly through the design space.Even though LHS produces points in  that are spread apart, this can still lead to some physical regimes being over-represented in the training data.Our previous work 4 uses Isomap, a manifold learning algorithm to produce local POD bases 5 for a highly nonlinear lid-driven cavity flow problem.The results from this work show that Isomap is effective at separating simulation data by physical regime.In this work, we propose an adaptive sampling algorithm using Isomap to develop physically diverse training datasets for ROMs.Our proposed adaptive sampling algorithm uses LHS to generate a baseline set of design parameters at which the full-order model is solved, and Isomap is then utilized to develop a greedy algorithm that iteratively selects the next design parameter from a set of candidate points to add to  given the current low-dimensional manifold of the training data until a desired number of samples is obtained.
Previous work in developing adaptive sampling algorithms for POD-based ROMs has mainly focused on non-intrusive methods and involves computing the SVD at each adaptation iteration.Braconnier et al. 6 use an a-posteriori error estimator based on a leave-one-out approach using the POD basis.Guenot et al. 7 propose two algorithms that select samples based on model improvement based on rotation of the POD vectors or through changes in the POD coefficients.Wang et al. 8 combine the two methods to develop a conjunction sampling strategy.Our proposed algorithm does not use POD, which is more computationally expensive than Isomap when the number of observations is relatively low.When used in our algorithm, Isomap produces a manifold of a fixed low dimension, whereas the dimensionality of the POD basis grows with the number of samples.In a work 9 by Franz et al., an Isomap-based adaptive sampling algorithm is introduced and applied to both a nonlinear Isomap-based ROM 10 and POD-based ROM.The adaptive sampling algorithm in that work requires solving an optimization problem to determine the added samples, and results are limited in scope, with analysis being conducted at a single test point.Computational costs and details of the optimization process are also not provided.The criteria used to select added points are different from the one presented in our work, although both are dictated by distances from points within the low-dimensional space.Our algorithm is also applicable to both non-intrusive and projection-based ROMs, whereas previous works have focused exclusively on non-intrusive ROMs.A wide range of adaptive sampling algorithms is compared in a work 11 by Karcher and Franz.When applied to two external aerodynamics problems, our results show that the proposed adaptive sampling algorithm outperforms using LHS alone in predicting physical fields and integral quantities of lift and drag.

FULL-ORDER MODEL
The full-order model (FOM) in this work is a system governed by a set of steady-state parameterized partial differential equations (PDEs) discretized over a computational domain Ω.A solution vector x() ∈ R N corresponds to design parameters  ∈  that control both the computational domain and parameters of the governing equations. denotes the parameter space such that x:  → R N .The numerical examples in our work focus on the steady Navier-Stokes equations applied to external aerodynamics.

Steady Reynolds-averaged Navier-Stokes equations
The computational fluid dynamics (CFD) simulations performed in this work use OpenFOAM, 12 an open-source toolbox for multiphysics simulation.Steady incompressible flow is simulated using simpleFoam, a standard OpenFOAM solver, by solving the Navier-Stokes equations, where is the velocity vector and u and v are the velocity components in the x and y directions respectively, S is the face-area vector, ⃗ n is the outward-pointing normal, V is the volume;  is the kinematic viscosity, and p is the pressure.The continuity and momentum equations are discretized over the computational domain by using the finite-volume method (FVM).Both equations are coupled through the semi-implicit method for pressure-linked equations (SIMPLE) algorithm 13 along with Rhie-Chow interpolation. 14The SIMPLE algorithm employs a pseudo time-stepping technique through a predictor-corrector method.In the predictor step, an intermediate velocity field is solved for, followed by the correction step where the velocity field is updated based on a pressure correction equation.The Reynolds-averaged Navier-Stokes (RANS) approach is used for turbulence modeling.The Spalart-Allmaras (SA) one-equation turbulence model is chosen, where ν is the modified viscosity, which is related to the turbulent eddy viscosity as The four terms in Equation (3) represent the turbulent convection, diffusion, production, and destruction respectively.The original work 15 by Spalart and Allmaras can be consulted for an in-depth overview of this model and detailed definitions of terms and closure coefficient values.

REDUCED-ORDER MODELS
This section gives an overview of the non-intrusive (ni-ROM) and projection-based (p-ROM) ROMs used in this work.Both ROMs utilize the proper orthogonal decomposition (POD), a popular method for constructing a reduced basis consisting of a small number of linearly independent basis vectors.A linear combination of these basis vectors allows for accurate representations of the full-order state.Non-intrusive ROMs use a regression model to predict the expansion coefficients at an unseen design,  * .Projection-based ROMs solve a low-dimensional version of the full-order model, relying on a reduced representation of the flow residuals to dictate convergence.

Proper orthogonal decomposition
Reduced-order models use training data obtained from a set of n solution snapshots calculated at chosen design points in the parameter space.A snapshot matrix, S ∈ R N×n , is assembled There exists a subspace  associated with S such that  = span(S).We assume that  provides a good approximation of the solution manifold for  ∈  if there are a sufficient number of solution snapshots in S, which corresponds to a judiciously chosen subset of design parameters in .A rank k set of orthonormal basis vectors where k ≪ N, is associated with  such that each full-order solution x i in S can be well-approximated as a linear combination of the basis vectors where a i is the set of low-dimensional expansion coefficients for a given solution snapshot.The truncated singular value decomposition (SVD) of S contains two orthonormal matrices U ∈ R N×n and V ∈ R n×n , as well as a diagonal matrix  ∈ R n×n such that where U contains a set of n left singular vectors that form an orthonormal basis for the column space of S, V contains a set of n right singular vectors that form an orthonormal basis for the row space of S, and diag( where k ≤ n.This preserves the most dominant features of the solution space and reduces the computational complexity of the ROM.To determine the value of k, the relative information content of the subspace is often used as a metric, and k is chosen such that E(k) ≥ , where  ∈ [0,1] is chosen dependent on the type of ROM being used, usually to a value  ≥ 0.95. 16Using the POD basis, full-order solutions at unseen design parameters x( * ) can be approximated as where the ROM provides a good approximation of the coefficients.A schematic of the POD is provided in Figure 1.

Non-intrusive ROM
This method relies on a functional dependence a = f () between the design parameters  and the expansion coefficients a. Non-intrusive ROMs use a regression model to predict the expansion coefficients given the design parameters, leading to a fully data-driven approach that does not incorporate the governing physics of the full-order model.There are many popular regression models used in non-intrusive ROMs, including artificial neural networks 17 (ANNs), radial basis function (RBF) interpolation, 18 and Gaussian process regression (GPR). 19While ANNs can provide superior performance to other regression methods, they require large amounts of training data to make accurate predictions at unseen points.Since ROMs are typically trained using a small amount of data, this limits the use of ANNs for predicting expansion coefficients.RBF and GPR are both similar in that they are non-parametric models that use kernel functions to model similarity between data points, leading to smooth functions.In non-intrusive ROMs, an individual regression model is usually used for each coefficient in a, leading to k different regression models, GPR utilizes a probabilistic framework using gradient-based optimization to determine kernel hyperparemeters, which often leads to better generalization than RBF.In this work, we will use GPR as a regression model for both the non-intrusive ROM and the adaptive sampling algorithm.

Gaussian process regression
A brief introduction to GPR is given in this section, and a more complete overview can be found in a work by Rasmussen. 20 Gaussian process (GP) is a set of random variables that follow a joint Gaussian distribution.In GPR, it is assumed that data are generated according to a GP with mean function m and covariance function , with some added Gaussian noise  ∼  (0,  2 y ), where  2 y is a very small term (typically on the magnitude of 10 −10 ).Using a finite amount of training data {, ā} and a prior joint Gaussian on the data, the prediction at a point  * is given by Using the properties of conditional Gaussian distributions, the conditional expectation of f ( * ) is given as where I is the identity matrix.In practice, the mean function m is set to the mean of the training outputs, and the inputs are scaled before training to obtain their standard score  where i and j refer to indices of the observation and input entry, respectively, and s is the sample standard deviation.There are many different kernels that can be used for the covariance function.A common one is the radial basis function (RBF) kernel where d is the Euclidean distance function and l represents the length scale, which is a hyperparameter.The predictive performance of the regression model is highly dependent upon the values of the hyperparameters.Gradient-based optimization is used to maximize the marginal log-likelihood of the training data to obtain an optimal set of hyperparameters,  opt ,

Projection-based ROM
The projection-based ROM that is used for our experiments comes from our previous work. 21While that work describes a projection-based ROM that uses hyper-reduction through the discrete empirical interpolation method (DEIM) 22 to approximate flow residuals, we will use a brute-force method that fully computes the full-order residuals.This method offers better accuracy, and the purpose of our work is to highlight the effectiveness of our adaptive sampling algorithm rather than to provide an optimal projection-based ROM.When applied to fluid dynamics problems, projection-based ROMs solve a low-dimensional version of the FOM through the full-order residuals R. The FOM governing equations can be expressed in terms of R and are a function of the state x and design parameters , The residuals can also be approximated using the POD basis  state approximation,

R(x, 𝝁) ≈ R(𝚿a, 𝝁). (20)
For a new set of design parameters  * , the expansion coefficients a * that satisfy Equation (20) are desired.A least-squares Petrov-Galerkin (LSPG) approach 23 is used for the projection-based ROM.Specifically, the expansion coefficients a * that minimize the L 2 norm of the approximate residuals are solved for, At the minimum point, L 2 a = 0, giving us an equation for a k-dimensional residual r, ⊳ Compute the new full-order state variable vector A Newton-Krylov approach is used to solve Equation (22).First, a reference full-order state x ref from the training data is used to compute an initial value for a * .The low-dimensional residual is computed as well as the Jacobian r a using a finite-difference method.A linear equation is then solved for Δa.A backtracking line search is performed to compute the updated expansion coefficients a n+1 , with the requirement that the norm of r decreases.If this does not happen, the step size is reduced by a factor of two, until m max iterations are reached, which is typically set to 5.This process is repeated until the norm of r is less than a prescribed tolerance r tol or the Newton-Krylov algorithm reaches n max iterations.Although we are only driving the norm of r to 0, the full-order residual is also expected to decrease according to the LSPG formulation.The Jacobian is explicitly computed for each Newton-Krylov iteration using a matrix-free method.Individual basis vectors  i are extracted from  and a matrix-free approach is used to compute the matrix-vector product, where  = 10 −6 is a small step size.The full-order residuals are evaluated once for each of the k columns in  using DAFoam , 24 an open-source adjoint derivative computation framework for OpenFOAM.The entire algorithm is outlined in Algorithm 1.

ISOMAP ADAPTIVE SAMPLING ALGORITHM
The adaptive sampling algorithm presented in this work is described in this section.The algorithm makes use of Isomap, a nonlinear dimensionality reduction technique.Isomap is highly versatile, with applications including gene and protein expressions, 25 climate data, 26 and turbulent combustion. 27Using a low-dimensional manifold of the current training data, the manifold is efficiently filled given a large set of candidate design parameters.This allows for training data sets that are more representative of the physics over the parameter space, which should allow for more accurate ROMs.

Isomap
A brief introduction to Isomap is given in this section, and the original work 28 by Tenenbaum et al. can be referred to for a more complete overview.Isomap is a nonlinear dimensionality reduction method that estimates a low-dimensional manifold of high-dimensional data.Given a data matrix X ∈ R n×N with n observations and dimensionality N, Isomap provides a latent representation matrix W ∈ R n×r , where r ≪ N. The latent representation w ∈ R r for each sample represents its position on the low-dimensional manifold.Isomap aims to produce a manifold that preserves the geodesic distances between points in the high-dimensional space.Given a nearest neighbor graph  of the high-dimensional training data where there are connections to the  nearest neighbors and the edges are weighted by the distance between points, the geodesic distance between two points is the shortest distance between them on the graph.The Euclidean distance is usually used as the metric of distance for the nearest neighbor graph, which can be computed using methods such as Dijkstra's algorithm.Once the geodesic distances between points are assembled into a distance matrix D ∈ R n×n , multidimensional scaling (MDS) 29 is applied to obtain W. MDS is an algorithm for producing low-dimensional representations that preserve high-dimensional pairwise distances.Using optimization, a stress function  that represents the discrepancy between the pairwise distances in the low and high dimensional space is minimized.The stress function is given as where d is the Euclidean distance between two latent representations.There are many ways to optimize the stress function in MDS.Isomap is implemented in this work using scikit-learn, 30 an open-source machine learning library for Python.The scikit-learn implementation of Isomap uses the SMACOF (Scaling by MAjorizing a COmplicated Function) algorithm, 31 which iteratively adjusts positions of points in the low-dimensional space using a majorization-minimization technique.
The coordinates of the low-dimensional points are first initialized randomly, and the SMACOF algorithm is iteratively repeated until w(D) is minimized.This leads to a manifold that represents the intrinsic geometry of the data.The computational aspects of the SMACOF algorithm are outside the scope of this work, and the original work can be consulted for more detail.The Isomap method is described in Algorithm 2. Figure 2 shows the manifold resulting from Isomap being applied to the swiss roll dataset, a popular three-dimensional benchmark problem for dimensionality reduction methods.The resultant two-dimensional manifold successfully separates points based on their positions along the roll, rather than Euclidean distances.Points which are relatively close to each other in terms of Euclidean distance can be very far from each other along the roll.When using Isomap, an important hyperparameter to consider is , the number of nearest neighbors used to construct .If  is too large, the local structure of the data manifold isn't represented well, while values of  that are too small can lead to inaccurate geodesic distances because of an outsized focus on small regions surrounding data points and neglect of the global structure of the data.

Input: X
⊳ High-dimensional data matrix Output: W ⊳ Low-dimensional latent representation matrix 1:  ← X ⊳ Compute a nearest neighbor graph from each observation in X with the edges weighted by distance.
2: D ←  ⊳ Compute the pairwise geodesic distances between all points.
3: W ← MDS(D) ⊳ Apply MDS and obtain the latent representation matrix.

Adaptive sampling algorithm
Our proposed adaptive sampling algorithm utilizes Isomap to find a low-dimensional physical representation of the current training data to make an informed decision when selecting the next design parameter to add to the dataset.Given the current set of design parameters in the training set  train ∈ R n×p and a very large set of design parameters  cand ∈ R M×p generated using LHS, the current manifold of the training data is filled to diversify the represented physical states.A number q i of initial samples generated using LHS is required, and the adaptation algorithm is iteratively repeated until a total of q t samples is obtained.First, Isomap is applied to the transpose of the current snapshot matrix S T to obtain a two-dimensional latent representation matrix W train .A two-dimensional representation is chosen since as the latent variable dimension grows above 2, the Euclidean distances between points become less meaningful due to increased sparsity and the curse of dimensionality.A set of candidate design parameters  cand is generated using LHS.By generating enough candidate points, the parameter space  can be efficiently represented.GPR is then used to approximate the latent representations W cand of the candidate design parameters, with a separate regression model used for each latent variable.The pairwise distances between all of the points in W train and W cand are then calculated and the minimum pairwise distance for each candidate point is retained in a vector d min ∈ R M .The candidate point corresponding to the maximum value in d min is chosen as the new design parameter  new to be added to  train .An example of this is shown in Figure 3.By looking at the minimum pairwise distances for each candidate point, we are filtering out points that are already similar to any of the current samples (i.e., the red marker in the lower left corner).To fill the manifold efficiently, the point with the greatest minimum distance from all current points is chosen, which is the point in the center.The FOM is then solved at  new and the solution x new is added to S. The time complexity of Isomap grows mainly with n, which is small, instead of N. The overall time complexity of the algorithm, including nearest neighbor graph construction, geodesic distance Example of the selection criteria for the adaptive sampling algorithm.
calculation, and MDS is given as  ( Nn log(n) + n 2 log(N) + n 2 ) when fixing the latent variable dimension and .The computational cost comes from the implementation of Isomap in scikit-learn.Using the same library, the time complexity of the truncated SVD is (Nn 2 ).The adaptation algorithm is fully outlined in Algorithm 3.

Algorithm 3. Adaptive sampling algorithm using Isomap
Input: S,  train ⊳ Initial snapshot matrix and design parameters Output: S,  train ⊳ Final snapshot matrix and design parameters 1: for j ∈ {q i + 1, q i + 2, … , q t } do ⊳ Iterate over all additional samples

RESULTS
The two test cases used to demonstrate the performance of the adaptive sampling algorithm are a 2D NACA 0012 airfoil and a 3D Cessna 172 wing.The design parameters for both cases involve geometric deformations using a free-form deformation (FFD) method through the pyGeo 32 package.While both the non-intrusive and projection-based ROM are used for the first test case, only the non-intrusive ROM is used for the Cessna 172 case.Although this test case appears in our previous work 21 presenting the p-ROM used in this work, a very small geometric parameter range was used where there was very little variation in physics and the adaptive sampling algorithm offers no advantage over LHS.When using a larger parameter range, the p-ROM lacks stability and sometimes fails to converge, by either failing to converge to the prescribed ROM tolerance and getting stuck at a bad local minima, or having r tol diverge.We attempted to remedy this issue by using the ni-ROM solution as an initial condition, using fewer basis vectors, and slightly increasing the p-ROM tolerance, but the issue persisted.This issue is common in projection-based ROMs and well-documented, 33 highlighting their limitations.The number of basis vectors k used for the ni-ROM is equal to the number of samples n.The GPR kernel used in the adaptation algorithm for predicting the latent representations is the rational quadratic kernel, which can be seen as a scale mixture of RBF kernels, given as where  is a scale mixture parameter.The kernel used to predict POD expansion coefficients is the Matern kernel, given as where Γ is the gamma function, K  is the modified Bessel function of the second kind, and  is a hyperparameter representing the smoothness of the model.These kernels were chosen through a trial-and-error process, where the RBF kernel was also considered.The Matern and rational quadratic kernels were both accurate for predicting latent variables and POD coefficients, and the difference in performance between the two was modest.The relative information content  from Equation ( 8) is set to 0.999 to select the number of basis vectors when using the p-ROM.We found that this strikes a good balance between the quality of the POD basis and the numerical stability of the p-ROM.The tolerance of the p-ROM is set to r tol = 100 and the maximum number of iterations is set to n max = 20.
The initial norm of the reduced residuals at validation points is typically on the order of 10 6 or 10 7 , and setting r tol = 100 requires a drop of 4-5 magnitudes.At a validation point, the p-ROM is initialized with a state x ref from the training samples corresponding to a set of design parameters closest in Euclidean distance to  * .For both test cases, the number of nearest neighbors  used for Isomap is set to 4. Since we are computing manifolds given a relatively small number of observations,  must be low to provide meaningful geodesic distances as well as high enough to acknowledge global data structure, which would be ignored with extremely low or high values.
The adaptive sampling algorithm is performed using the full-order state x = [u, v, p] while both ROMs use x = [u, v, w, p,  t , ], where  is the cell-face flux.We choose a subset of the flow variables for the adaptive sampling algorithm as they are the ones that we are most interested in maximizing the accuracy of and are more representative of the flow physics, while the ROMs require access to all flow variables.For both test problems, M = 10 4 candidate points are generated using LHS for each iteration of the adaptation algorithm and a set of N V = 100 validation points generated using LHS is used to assess performance.Four different sample levels q t = [20, 24, 28, 32] are used for both test problems.Additionally, a number of trials N t at each sample level is simulated to average performance metrics over.The average relative L 2 error of a field quantity for a single trial  t over all of the validation points indexed by i is given as where x is the field approximated by the ROM.The errors are then averaged over each trial to give a single error value .
Relative errors in the coefficients of lift and drag C L and C D between the CFD solver and ROMs are also evaluated using the same averaging procedure.Results are presented using an initial number of LHS-generated samples q i = q t 2 .We choose this value as a sufficient number of initial samples is required to obtain a meaningful low-dimensional manifold from Isomap.Additionally, we desire our final data set to be composed of an adequate number of adaptively chosen samples.Results using different initial fractions of q i for the NACA 0012 test case are available in the Appendix, along with performance comparisons to a POD-based adaptive sampling algorithm.

NACA 0012 airfoil
A NACA 0012 airfoil, shown in Figure 4, is used as our first test case.The chord measures 1 m and the span is 0.1 m, with one cell in the spanwise direction.The flow is incompressible with Re = 10 6 with U ∞ = 10 m/s.A coarse mesh F I G U R E 4 Mesh and FFD points for the NACA 0012 case.
with 23,182 cells is generated and the computational domain extends 30 chords from the airfoil surface.Four FFD points control the vertical displacements at the leading and trailing edges, with the bounds of the design parameters given by  ∈ [−0.03, 0.03] m, and the baseline angle of attack  is set to 4 degrees.The full-order model is run for 2000 iterations, after which the flow residuals fall to at least 10 −6 .and a total of N t = 20 trials is simulated at each sample level.Table 1 lists computational costs associated with the ROM and adaptive sampling algorithm.It is shown that the computational cost of Isomap is significantly lower than that of the SVD when using n = 32 samples and negligible compared to a single CFD simulation.This results in a very computationally efficient adaptive sampling algorithm, especially when comparing the cost to a single CFD simulation.While the wall time for the p-ROM varies with the number of iterations and residual evaluations, an average wall time of approximately 3.4 s was found, which leads to an average speed-up of around 11×, comparable to what was found in the original work 21 for a similar test case.Field contours of u, v and p from the CFD solver, ni-ROM, and p-ROM are shown in Figure 5 at a validation point  * = [0.017,0.025, −0.0297, −0.0273] with q t = 32 samples generated using LHS.Both ROMs are in close agreement with the FOM, although the p-ROM is more accurate, which is evident when looking at contours of u and p. Figure 6 shows the average relative field errors over the chosen sample levels for the ni-ROM (left) and p-ROM (right), while Figure 7 shows the average relative errors in C D and C L .When comparing similar field quantities, it is shown that the p-ROM offers much better predictive performance, even when using a lower number of training samples.Using the adaptive sampling algorithm results in significantly lower errors for all field quantities when using the p-ROM, and for v and p when using

F I G U R E 13
Comparison of average L 2 relative errors in field quantities.
the ni-ROM.This effect is greater when using the p-ROM; as it is a physics-based model, we can expect that a greater diversity of training samples will lead to a better gain in predictive performance compared to the ni-ROM.Errors in drag and lift exhibit a similar pattern, where both the predictive performance and reduction achieved by using the adaptive sampling algorithm are greater for the p-ROM when compared to the ni-ROM.There is also a greater relative decrease in force coefficient errors compared to field errors.Figure 8 shows the average number of basis vectors used for the p-ROM for both adaptation and LHS alone.For the same singular value threshold, using the adaptation algorithm results in a higher number of basis vectors at all of the sample levels, which shows that the training samples are more physically diverse.Building an ni-ROM with 28 total samples offers better predictive performance in v, p, C D , and C L compared to using LHS alone.The same is true for building a p-ROM with 24 samples, where the relative gain in predictive performance is larger, again showing that the p-ROM benefits more from using adaptively chosen samples.For a smaller computational budget, using the adaptive sampling algorithm leads to better overall predictive performance of both field and integral quantities for both purely data-driven and physics-based ROMs.Figures 9 and 10 show field error contours for a single trial with errors close to the reported average relative errors.The errors are averaged over all over the validation points for the ni-ROM and p-ROM respectively.The contours represent the absolute difference between the CFD solver and ROM.For all field quantities, errors are greater closer to the surface of the airfoil and negligible farther away.It should be noted that the relative L 2 errors reported consider the error over the entire computational domain; since the majority of the error is close to the airfoil surface, this can explain the larger relative decrease in force coefficient errors compared to field errors.F I G U R E 15 Mean mid-section field error contours for a single Cessna 172 trial using LHS (left) and adaptation (right) with q t = 32.

Cessna 172 wing
A Cessna 172 wing, shown in Figure 11, is used as the next test case.The root chord of the wing measures 1.67 m, and the semi-span aspect ratio is 3.  2, where it is again shown that the cost of Isomap is very low compared to the FOM and SVD, making the adaptation algorithm very efficient for this case.
Figure 12 shows mid-section field contours of u, v, and p at a validation point  * = [−2.19,−1.47, −2.79, −2.97, 2.79] from both the CFD solver and ni-ROM with q t = 32 samples generated using LHS.Visually, there is very good agreement between the FOM and ROM in all three fields.The average relative errors in the fields are shown in Figure 13 for the chosen sample levels.Using the Isomap adaptation algorithm leads to significantly lower errors in both u and p over all the sample levels and a slight error reduction in v.The average relative errors in C D and C L are shown in Figure 14, where we see that using the adaptation algorithm leads to very large reductions in the computed output errors.Building a ROM using the adaptation algorithm with 28 samples compared to 32 samples using LHS leads to comparable relative errors in u, v, and p and significantly lower errors in C D and C L , showing that the adaptation algorithm leads to better predictive performance given a smaller computational budget, similar to the airfoil test case.When looking at average error contours in Figure 15 for a single trial with relative errors close to the reported averages, it is again shown that using the adaptation algorithm leads to considerably lower errors in the computed fields closer to the surface of the wing, leading to a larger drop in error in coefficients compared to the fields, where the error metric is computed over the entire domain.Far away from the wing surface, errors in the fields are similarly low when using the adaptation algorithm and LHS alone.

CONCLUSION
This work presents a novel adaptive sampling algorithm applicable to both non-intrusive and projection-based reduced-order models using Isomap, a highly versatile and computationally efficient algorithm for obtaining low-dimensional intrinsic manifolds of high-dimensional data.Isomap accounts for the global structure of data by preserving the geodesic distances between points in the low-dimensional manifold.The algorithm iteratively fills the manifold of the current data set using a large number of candidate samples generated using Latin hypercube sampling and Gaussian process regression to predict their latent representations.By filling gaps in the manifold, there will be a greater diversity of physical regimes present in the training data compared to using LHS, which only accounts for the distribution of the design parameters and not the physics of the full-order model.Additionally, our algorithm is applied to both non-intrusive and projection-based ROMs.
When applied to two steady, incompressible, external aerodynamics problems, our results show that using the adaptive sampling algorithm leads to significantly improved predictive performance in fields and coefficients of lift and drag compared to using LHS alone over a number of sample levels.For a given singular value threshold, using adaptively chosen samples leads to a higher number of basis vectors being used for the p-ROM, showing that the diversity of the training data increases.For a lower computational budget, using the adaptive sampling algorithm leads to improved predictive performance over LHS.Our results are averaged over a number of trials to account for the randomness associated with generating samples using LHS.Comparisons to a POD-based algorithm from another work show that our algorithm performs better when using the p-ROM, and slightly worse when using the ni-ROM.A limitation of this work is that the proposed algorithm is only applicable to steady-state problems.Applying this algorithm to unsteady problems would likely require including multiple time snapshots per simulation for Isomap and using time as an input to the GPR model.Future work will focus on developing a similar algorithm for ROMs applied to unsteady problems and testing the performance on ROMs which use deep neural networks instead of POD to provide a low-dimensional solution space.
F I G U R E A1 Computational costs associated with the Isomap and POD-based adaptive sampling algorithms for the NACA 0012 case.
Using the same experimental setup as the results section, Table A1 shows the average relative errors between the two adaptive sampling algorithms and LHS with q t = 32 and q i = q t 2 for the POD-based algorithm offers better performance in predicting u, v, and C D , while our algorithm offers better performance in predicting p and C L .u in particular is much better predicted using the POD-based algorithm.A reason for this may be that the variation in u between different samples does not contribute as much to the geodesic distances between points as v and p do, and additional samples that vary these quantities more are selected.Table A2 shows that our algorithm outperforms the POD-based one in predicting all quantities except u, which further suggests that our algorithm is better suited for projection-based ROMs.The computational costs of both algorithms as a function of the number of current samples is given in Figure A1.The given computed costs include all parts of the algorithms excluding the generation of candidate points using LHS.Our algorithm does not display any significant trend with the number of samples, and variations can be due to the GPR algorithm or CPU background processes.The POD-based algorithm does exhibit an upward trend with the number of samples, most likely due to an increased cost with computing the POD basis as more samples are added.While the difference in computational cost between the two algorithms for this case is low, we can expect that for larger problems such as the Cessna 172 test case, that the Isomap-based algorithm will be faster given the costs listed in Table 2, since costs apart from Isomap and POD will not vary with N.

APPENDIX B. EFFECT OF INITIAL SAMPLE SIZE
The effect of the fraction of initial LHS-generated samples is investigated when applied to the NACA 0012 test case with q t = 32 for both the ni-ROM and p-ROM.In partciular, the relative errors in the fields and force coefficients are compared to using LHS alone with q i = ( q t 4 ,  ) with t = 20 trials, similar to the data presented in the results.Figures B1 and B2 show that there is there is no distinguishable trend between the error metrics and number of initial samples for the ni-ROM.Apart from the average relative error in u which was already previously shown to be close to that of using LHS alone for the ni-ROM, the predictive performance in all other quantities is significantly better than using LHS alone.The results suggest that the effect of including adaptively chosen samples has diminishing returns as more samples are added.All of the error metrics are slightly higher for the p-ROM with q i = 3q t 4 .The p-ROM is more sensitive to the diversity of the training data, and using fewer adaptively chosen samples can lead to a relatively larger degradation in performance compared to the ni-ROM.When using the adaptation algorithm, an important consideration is the accuracy of the GPR model in predicting the latent variables at candidate points.Using too few initial samples, especially for highly non-linear problems, can lead to an inaccurate latent variable regression model and to poorly chosen samples.

F I G U R E 1
Schematic of the proper orthogonal decomposition (POD), where the snapshot matrix S is decomposed using the singular value decomposition (SVD) and the POD basis  is obtained from U.

F I G U R E 2
Two-dimensional manifold (right) obtained from Isomap being applied to the swiss roll dataset (left).

TA B L E 1 5 F I G U R E 5 F I G U R E 7 F I G U R E 8
Computational costs associated with the NACA 0012 case.Comparison of NACA 0012 field contours with q t = 32 at  =[0.017, 0.025, −0.0297, −0.0273].Comparison of average L 2 relative errors in C D and C L for the ni-ROM (left) and p-ROM (right).Average number of basis vectors used for the p-ROM.

F I G U R E 9 F
Mean ni-ROM field error contours for a single NACA 0012 trial using LHS (left) and adaptation (right) with q t = 32.F I G U R E 10Mean p-ROM field error contours for a single NACA 0012 trial using LHS (left) and adaptation (right) with q t = 32.F I G U R E 11Mesh and FFD points for the Cessna 172 case.TA B L E 2 Computational costs associated with the Cessna 172 case.I G U R E 12 Comparison of Cessna 172 mid-section field contours with q t = 32 at  = [−2.19,−1.47, −2.79, −2.97, 2.79].

F I G U R E 14
Comparison of average relative errors in C D and C L .

2 .
The Reynolds number is set to Re = 7 × 10 6 with a freestream velocity of U ∞ = 63.8 m/s.A structured mesh with N = 609,280 cells is generated and the computational domain extends 30 chords from the surface.100 FFD points are used to deform the wing surface at five spanwise locations.The design parameters are the five twist angles at these spanwise locations, which are obtained by rotating each set of FFD points.The bounds of the twist variables are  ∈ [−3, 3] degrees, and the baseline angle of attack  of the wing is set to 2.5 degrees.The full-order model is run for 1000 iterations, after which the flow residuals fall to at least 10 −6 .A total of N t = 5 trials are simulated at each sample level.Computational costs for running the CFD solver as well as Isomap and the SVD with n = 32 are listed in Table

3q t 4
2 , … ,  n ] contains the singular values corresponding to the singular vectors in non-ascending order,  1 ≥ • • • ≥  n ≥ 0. The first k vectors of U are chosen to be the vectors forming the POD basis,  ∈ R N×k = [ 1 ,  2 , … ,  k ].The singular values associated with the basis vectors often decay rapidly, and hence k singular vectors are chosen to form the POD basis,