{'fraction_non_alphanumeric': 0.02957692229722881, 'fraction_numerical': 0.012122686452188367, 'mean_word_length': 4.59297052154195, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Polynomial Ridge Flowfield Estimation

Computational fluid dynamics plays a key role in the design process across many industries. Recently, there has been increasing interest in data-driven methods, in order to exploit the large volume of data generated by such computations. This paper introduces the idea of using spatially correlated polynomial ridge functions for rapid flowfield estimation. Dimension reducing ridge functions are obtained for numerous points within training flowfields. The functions can then be used to predict flow variables for new, previously unseen, flowfields. Their dimension reducing nature alleviates the problems associated with visualising high dimensional datasets, enabling improved understanding of design spaces and potentially providing valuable physical insights. The proposed framework is computationally efficient; consisting of either readily parallelisable tasks, or linear algebra operations. To further reduce the computational cost, ridge functions need only be computed at only a small number of subsampled locations. The flow physics encoded within covariance matrices obtained from the training flowfields can then be used to predict flow quantities, conditional upon those predicted by the ridge functions at the sampled points. To demonstrate the efficacy of the framework, the incompressible flow around an ensemble of aerofoils is used as a test case. On unseen aerofoils the ridge functions' predictive accuracy is found to be reasonably competitive with a state-of-the-art convolutional neural network (CNN). The local ridge functions can also be reused to obtain surrogate models for integral quantities such a loss coefficient, which is advantageous in situations where long-term storage of the CFD data is problematic. Finally, use of the ridge framework with varying boundary conditions is demonstrated on a three dimensional transonic wing flow.


Introduction
The industrial penetration of Reynolds Averaged Navier-Stokes (RANS) has gone far beyond its mainstay of aerospace engineering. Today it is being heavily used in automative design [1], building design [2], ventilation system design [3], marine propeller design [4] and more recently in the understanding of oral transmissions [5,6]. In using a turbulence model, RANS avoids the need to resolve pertinent length-and time-scales of turbulence; significantly slashing the cost of generating a flowfield. That said, owing to the size and increasingly realistic mesh topologies, these RANS computations can themselves take several hours on modern clusters. Additionally, in both design and analysis contexts, there is a need to assess multiple geometry and boundary conditions, adding to the computational cost.
It is therefore unsurprising that many have sought to develop ways of (i) efficiently representing a flowfield and (ii) approximating it for a new set of boundary conditions or geometry definition. Falling under the remit of sparse reconstruction, there exists a variety of approaches for reconstructing full unsteady flowfield information from a limited Although applications of deep learning networks to real-world problems have become ubiquitous, our understanding of why they are so effective is lacking [18]. Techniques for interpretation are available [19], but they may not be particularly accessible for those without some specialised deep learning knowledge. Moreover, for problems with large numbers of features, widely adopted approaches such as those based on Shapley values [20] quickly become intractable. Consequently, in many fields of engineering, there is still a demand for more readily interpretable surrogate models, even if these bring a slight trade-off in predictive accuracy. One approach to easing the challenge of interpretation, is to reduce the dimensionality of the problem. By reducing the dimensionality of the input space, data-driven ridge functions [21] facilitate straightforward visualisation and understanding, whilst at the same time reducing data requirements. del Rosario et al. [22] show how ridge functions can offer important insights during the conceptual design phase, and in [23], ridge functions aid exploration of a stagnation temperature probe's design space. Wong et al. [24] introduce the idea of using dimension reduction to estimate flow quantities, embedding polynomial ridge functions on an aerofoil surface to estimate the pressure distribution. In this paper, we develop spatially correlated ridge functions, with the target output being flowfield variables such as pressure, velocity or even turbulent variables.
To demonstrate the utility of the proposed framework, we begin by exploring the prediction of flowfields around a set of aerofoils. The design of efficient aerofoils is important since their use is ubiquitous across many industries; for example aircraft wings and controls surfaces, car spoilers, and the gas turbine blades found in aircraft engines and power generation turbines. To provide a point of reference, we compare and contrast the proposed approach to a state-of-the-art flowfield prediction framework based on a convolutional neural network. In addition to examining the predictive accuracy, we explore how the dimension reducing nature of the learned model aids understanding of the existing training dataset. To assess the limitations of the proposed approach, we then tackle the transonic flow over a three dimensional wing, with the inflow boundary conditions varied. We examine how the underlying polynomials perform when extrapolating away from the training data, and discuss approaches to handling this situation. For a hands-on exploration of the framework, example code and an interactive web app are available at github.com/ascillitoe/ flowfield_approx.

Mathematical foundations
In this section, we detail some of the mathematical ideas that underpin this paper.

Flowfield representation
It will be useful to represent a steady flowfield as a set of related scalar-field quantities. These include, but are not limited to, static pressure p, density ρ, a velocity component i.e. v x , and even the turbulent viscosity ν t . These fluid properties are assumed known for any particular discretised location s ∈ R 3 within the three-dimensional flow domain S. We consider s to be synonymous with the number of nodes in the flowfield, where s 1 corresponds to the first coordinate and s N corresponds to the last.
Ignoring model parameters and limiters embedded within any modern flow-solver, one can group the parameters influencing a flowfield into (i) those that govern the mesh and its associated geometry, and (ii) those that set the boundary conditions. For simplicity, we condense the former into a vector of shape parameters x ∈ R d arising from an appropriately chosen design space X, e.g., Hicks-Henne bump functions, free-form deformation amplitude, or the inputs to an autoencoder that characterises the geometry [25]. We characterise the boundary conditions by a vector v ∈ R b where b represents the number of input boundary conditions, e.g., inlet stagnation pressure, exit static pressure, isothermal walls, etc.
Then, following the statement at the start of this section, we can represent the steady flowfield as a set of related vectors Equal contributions from all input dimensions, with a unit vector W = 1/ √ 25 b) Dimension reduction to obtain W , with a cubic orthogonal polynomial fitted over the subspace Figure 1: Subspace-based dimension reduction applied to a turbo-machinery problem [35]. The output scalar-valued quantity of interest f (x), the efficiency of a turbo-machinery blade, is a function of the d = 25 design parameters in x ∈ R 25 . The quantity f is projected over two different subspaces W T x. In a), W is chosen so that the subspace is formed from an equal contribution of all 25 parameters, whilst in b), a single linear combination of the 25 parameters which best describes the variation in f is obtained.
with X ∈ R d×M and F ∈ R N ×M . As before N represents the number of spatial nodes; d the dimension associated with the design space x ∈ X, and M the number of distinct flow-fields, each corresponding to a new design. At an isolated node s i across all the flow-fields, we wish to identify a polynomial ridge approximation of the form given in (3). To elaborate on the form of this ridge approximation, we rewrite it as where represents a Vandermonde-type matrix with n-variate polynomial basis terms {φ 1 , . . . , φ L } with a cardinality of L, with unknown polynomial coefficients α i ∈ R L . The cardinality here is defined as the total number of basis terms, and will vary depending on the highest degree in along each dimension and the nature of the interactive polynomial terms between different dimensions. For further details on the construction of multivariate polynomials, we refer the interest reader to Section 1.1 in [36]. The intention behind phrasing our objective as (7), is that for a new design x * , we can estimate the scalar field value at any node s i using the computed polynomial ridge approximation. To ensure solutions to this problem are feasible, it will be useful to ensure that the polynomial matrix P is well conditioned. This can be aided by adopting orthogonal polynomial basis terms-i.e., Legendre, Hermite or Jacobi-as these have demonstrably lower condition numbers than their monomial counterparts. In this paper, we utilise Legendre polynomials. This has the additional advantage of being able to rapidly yield moments and sensitivities (see [37]).
The ridge approximation in (5) can be obtained by solving the non-linear least squares problem which can be simplified to an optimisation problem over the Grassmann manifold via the variable projection method [38] minimise Wi∈R d×n where the superscript † denotes the matrix pseudoinverse, and O ⊥ is the orthogonal projector onto the complement of the column space of P . Gradients for this objective can be readily computed and used in a Gauss-Newton algorithm over the Grassman manifold, as demonstrated in [32]. The computation of the coefficients α i trivially follows. This enables us to compute the scalar field quantity at a given node i for a new design x * with P W T i x * α i and by extension the scalar field across all N nodes via  which by construction constraints all subspaces W 1 , . . . , W N to have the same dimension. While the calculation of polynomial ridge approximations across all N nodes is an embarrassingly parallel operation-as each the input-output data pairs at each node are treated independently-it will be useful to leverage the spatial correlations in a given scalar field to reduce the number of times (10) has to be solved.

Spatial correlations
Subspace-based dimension reduction approaches and the related active subspaces [39] have enjoyed tremendous computational success over the past few years-ushering in a shift from prior sparse grids and full-space design of experiment approaches. As a general observation, within computational fluid dynamics applications, these methods have been primarily applied on integral quantities of scalar field outputs, i.e., efficiency, lift, drag, etc. Wong et al. [24] recognised this, and posited that approximating the composite scalar fields of such integral quantities using dimension reducing subspaces may not only offer greater insight into significant flow features, but may even further reduce the number of simulations required to do so. One important idea born from their work is that the ridge approximations of nodes that are adjacent to each other in a flowfield are likely to be similar in the way they depend on input parameters. In locally smooth regions of the flowfield, nearby nodes are predominantly affected by a similar subset of inputs, resulting in similarity in their ridge subspaces. This notion is quantified via the subspace distance, and in Wong et al. [24] the authors demonstrate that this property enables efficient compression and recovery of flow fields. In what follows, we build upon this idea using a different formalism, leveraging some of our preliminary work in [40].
Using the data available in F , we construct a sample covariance matrix K ∈ R N ×N of the form where E f m denotes the mean of the m-th scalar field across all the N nodes. It is important to note that this covariance matrix captures the spatial correlations across all the nodes for a given scalar field quantity. Its diagonal represents the sample variances in the scalar field at each of the N nodes.
We partition the scalar field across all N nodes into two subsets: a smaller subset J with J nodes (ŝ 1 , . . . ,ŝ J ), and a subset R containing the remaining N − J nodes s 1 , . . . , s N −J with J << N . We assume that the scalar field for any design x * across all N nodes can be expressed as a sample from the multivariate normal distribution N (E [f ] , K), where the first argument represents the mean and the second the covariance. The covariance matrix in (12) can be partitioned into (N − J) × (N − J) and J × J blocks As the scalar fields arising from this partition will have Gaussian marginals, we can use the Schur complement to approximate the scalar field at the nodes in R by computing the ridge approximations at the nodes in the smaller set J only  In a nutshell, this implies that we only need to solve (10) J times rather than that N . However, this approach relies on the assumption that a suitable W i is found at each point in J , such that the projection W T i x at each point in the subset explains a large percentage of the total variance at that point. Without this assumption, the resulting correlation matrix K will become increasingly diagonal in higher dimensions, and we will be unable to recover the scalar field at the remaining nodes in R. To alleviate this issue, poor quality ridge approximations are removed from the smaller subset J , and added back into R. Various heuristics can be used to assess the quality of ridge approximations. In this paper, we use the training R 2 scores, as discussed in Section 4.5. Finally, for the stability of this procedure it is important that the inverse of K 22 is computed carefully. In the present work, this is done via Cholesky decomposition, with a small constant added to the diagonal to prevent singular matrices.

Computational setup
In this section we outline the two test cases examined in this paper, in addition to summarising the implementation of the ridge approximation framework and the CNN it will be compared to.

Test Cases
Within an engineering context, flowfield predictions around a given object are often required at a range of operating conditions. Whereas in other cases, flowfield predictions may be required at a given operating condition, whilst the object itself is altered. To examine both scenarios we use two test cases; 1) a two dimensional subsonic aerofoil flow, with a large number of different aerofoil designs; and 2) a three dimensional transonic wing flow, with the freestream conditions varied.

Exploring the design of a subsonic aerofoil
We start with the well known NACA0012 aerofoil, discretised with a 449x129 curvilinear C-grid 2 , and deform it to obtain an ensemble of aerofoil designs. The aerofoil surface is deformed using d = 50 Hicks-Henne bump functions [42] where s 2,base are the s 2 coordinates of the baseline aerofoil, b j is the j th bump function, and s 1 and s 2 are normalised by the aerofoil's axial chord length C 1 . The bump amplitudes [β 1 , . . . , β d ] are then stored within the input vector x m ∈ R d for each m th design. To generate a dataset for training and testing, we create a (M = 2000)-point design of experiments (DoE) with uniformly distributed Monte Carlo samples for x. A random sample of the resulting aerofoil designs, as well as the full range of deformations, is shown in Figure 2. Flowfields are simulated for each design using the incompressible solver of the SU2 CFD code [43]. The commonly used Spalart-Allmaras RANS model [44] is used to represent the effects of turbulence, with the freestream turbulence viscosity ratio set to (ν t /ν) ∞ = 5. The freestream velocity magnitude is set to U ∞ = 1 m/s, and the reference static pressure is set to zero. The laminar viscosity ν ∞ is then set to give a Reynolds number of Re = U ∞ C 1 /ν ∞ = 6 × 10 6 . Each design is run at three angles of incidence α ∞ = 0 • , 10 • and 15 • , leading to a dataset consisting of 6000 flowfields in total.
The grid deformations performed when perturbing the baseline aerofoil to reach each new design mean that, even away from the aerofoil surface, grid points are at a slightly different location for each design. Since we wish to learn a functional mapping for the flowfield variables at fixed points in space, we re-sample 3 each flowfield onto the 90x318 Cartesian grid shown in Figure 3, and remove points lying inside the solid region of the aerofoil. We choose a Cartesian grid here since we are interested in full flowfield visualisations, but the approach is equally applicable to individual points, surfaces (see Ref. [24]) or planes. After the above procedure we are left with N = 90 × 318 = 28.6 × 10 3 output values f i,m for each input x m , forming M input/output pairs (x m , f m ) M m=1 , for each of the three angles of incidence investigated. The scalar f i,m is the field variable at the i th grid point for the m th design, where we take the static pressure p, velocities u, v, and turbulent viscosity ν t as the field variables to predict. We normalise the aforementioned variables by taking the static pressure

Varying freestream conditions for a transonic wing
To explore predictions with varying operating conditions, we take the compressible inviscid flow around the well known ONERA M6 transonic wing test case [46]. The salient details of this case are presented in Figure 4. It is discretised by a N = 108.4 × 10 3 point unstructured grid, with no boundary layer grid necessary due to the inviscid nature of the flow. Freestream boundaries are placed 10.5 root chord length's away from the wing, with the freestream Mach number and angle of incidence set at M a ∞ = 0.8395 and α ∞ = 3.06 • for the baseline case. A span-wise symmetry boundary is then placed at s 3 = 0. The flows are simulated using the compressible inviscid solver of the SU2 CFD code [43], and in Figure 4 the resulting Mach contours are shown on the surface of the wing for the baseline case. The flowfield here bears a close resemblance to the real viscous flow, with a characteristic "lambda shock" pattern visible on the upper surface of the wing. Figure 4 is the design of experiments used to explore the freestream conditions. For the training data, operating points for M a ∞ and α ∞ are obtained by applying six-point Gauss-Lobatto quadrature rules over the intervals M a ∞ ∈ [0.6, 1.02] and α ∞ ∈ [0, 6], leading to a total of M train = 36 training flows. For testing, M test = 25 points are sampled from the multivariate normal distribution N (µ, σ), with µ = [0.8, 3] T and σ = [0.14, 3] T . Unlike in test case 1 the geometry is fixed here, therefore there is no need to resample the flowfields onto a base grid. For the scalar field variable f i,m we consider the local Mach number M a.

Polynomial ridge function implementation
To obtain dimension reducing ridges, we use the equadratures [47] python library which contains an implementation of the polynomial variable projection method described in Section 2.3. Since each ridge function is independent of one another, we parallelise the task in an embarrassingly parallel fashion, with computation of individual ridge approximations performed by individual Python worker processes. The code is parallelised using the joblib library with the loky backend. Each parallel Python process is run concurrently on a separate computational core on a F72s_v2 virtual machine (36 physical CPU's, 144 GiB memory) on the Microsoft Azure cloud computing service. For the N = 28.6 × 10 3 points in test case 1, with 500 training designs, training takes approximately 20 minutes per field variable (≈ 24 ridge approximations per second). Note that training times can be substantially reduced by taking advantage of spatial correlation in the manner described in Section 4.5.
The flowfield reconstruction procedure described in Section 4.5 is entirely made up of linear algebra operations, thus it is straightforward to implement efficiently with existing linear algebra libraries. For this we use the numpy [48] python library with Intel ® MKL. To limit computation and storage requirements, we compute only the K 12 and K 22 blocks of the covariance matrix in (13), since only these blocks are required by the Schur complement step in (14). If further cost savings are required, low-rank approximations for K can be obtained, as is done in [40].

Convolutional neural network implementation
Thuerey et al. [14] explored a number of CNN architectures for flowfield predictions, and found the U-Net architecture [49] to be the most successful. In this paper we implement a modified version of the framework proposed by Thuerey et al. [14], shown in Figure 5. The U-Net architecture consists of an encoder, which progressively down-samples the 128 × 128 × 4 input data with strided convolutions. The four input channels consist of a boolean mask to define the aerofoil geometry, and three uniform input channels defining the freestream conditions Re, α ∞ , and (ν t /ν) ∞ . This allows the network to extract increasingly large-scale and abstract information as the number of feature channels grows, until we are left with a 1 × 1 × 256 derived feature vector. The decoder then does the opposite, with depooling layers reducing the number of features while increasing the spatial resolution. Eventually we are left with four 128 × 128 output channels, consisting of flowfield estimates for C p , u/U ∞ , v/U ∞ and ν t /ν ∞ . Skip connections help the network to consider low-level input information during the reconstruction of the solution in the decoding layers. The weights w with· denoting the CNN prediction for the field variable, and φ hub denoting the Huber loss function [50]. As seen in (16), the error is averaged over all M number of training designs and N number of grid points. This loss function is differentiable, enabling it to be back-propagated into the network in order to compute the weight gradient ∇ w L, which is important for the optimisation (learning) process.
The complexity of the above CNN is altered in this paper by varying the number of channels (and therefore the number of parameters), with the 488k parameter CNN in Figure 5 found to offer the best compromise between training times and accuracy. All CNN's detailed in this paper are implemented in the PyTorch python library, and trained on a GTX970 1664 core NVIDIA GPU. Training times range from 50 seconds to 80 minutes, depending on the complexity of the network and the number of training designs. Further details are included in Appendix A.

Test case 1: Design of a Subsonic Aerofoil
As an initial demonstration, ridge approximations are obtained for the static pressure coefficient C p = (p 0∞ − p 0 )/(p 0∞ − p ∞ ), normalised axial velocity u/U ∞ , and turbulent viscosity ratio ν t /ν, with 500 aerofoil designs used for training. Flowfield predictions for a new aerofoil design can then be obtained by transforming the new design vector x to the reduced dimensional space, and evaluating the polynomial ridge approximation p i W T i x; s i at each point i = 1, . . . , N in the flowfield. In Figure 6, the predicted normalised u velocity is shown for a randomly selected design from the test set. Qualitatively, the predictions are in close agreement to the true velocity field from the CFD solution, shown by the black iso-lines. The design shown here was not used during training, hence the ridge approximations appear to be able to learn enough of the flow physics to be able to make predictions for the new design.
In addition to making flowfield predictions, a second key aspect of the ridge approximations are sufficient summary plots, seen in the insets within Figure 6. Applied to aerofoil design, they provide visualisation of how the quantity of interest, in this case the u velocity at selected locations, is altered as the shape is changed. This is made possible since the sufficient summary plot shows the u velocity over the reduced dimensional space W T x ∈ R, instead of the full space x ∈ R 50 . The low scatter of train and test points in the summary plots highlights that low dimensional structure has been successfully found.

Accuracy of Ridge Approximations
To explore how the accuracy of the proposed ridge approximation framework compares to a state-of-the-art convolutional neural network (CNN), ridge approximations and the 488k parameter U-Net in Figure 5 are trained on C p and u/U ∞ around the aerofoil at two angles of incidence, α ∞ = 0 • and α ∞ = 10 • . To enable fair accuracy comparisons the ridge approximations are obtained on the same 128 2 grid used for the CNN instead of the finer 90x318 grid used elsewhere. The normalised mean absolute error of predictions is then measured for the training and test designs. At each point, the absolute error is averaged for M number of designs, then normalised by the standard deviation of the true data at that point. The result is averaged over all N number of points. This is done whilst varying the number of training designs, in order to obtain the error curves in Figure 7. Both the ridge approximation and CNN prediction errors exhibit the same general trends; for low numbers of training designs the test errors are considerably higher than the training errors, suggesting a lack of generalisation (i.e. the models are overfitting to the training data). However the test errors are significantly decreased as the number of training designs is increased. Perhaps surprisingly, despite the data hungry reputation of deep neural networks, the CNN framework is able to achieve low accuracies with a relatively small number of training designs (M AE σ < 10% with M train < 100).
As noted by Ronneberger et al. [49], this is a key benefit of the U-Net and other CNN architectures; the strong use of data augmentation via convolutions allows the available training data to be used more efficiently (compared to a fully connected deep neural network).
Encouragingly, for the α ∞ = 0 • case (Figs 7a and 7b) the ridge approximations are able to achieve mean test errors which are competitive with the CNN. At the higher angle of incidence of α ∞ = 10 • (Figs 7c and 7d) the flowfield is more complex, with stronger streamline curvature leading to larger pressure gradients and greater non-linear behaviour. This case is more challenging for the ridge approximations, with test errors of M AE σ ≈ 4% compared to the CNN test errors of M AE σ ≈ 1%. However, accuracies of under 4% would be acceptable for many preliminary design applications. Additionally, as the MAE fields in Figure 8 show, the ridge approximations' prediction errors are actually lower than the CNN errors in many regions of the flow. Compared to the rather randomly distributed CNN errors, the ridge approximations' errors exhibit a smoother spatial distribution, with higher errors near to the aerofoil surface and in the wake.
In the above results, ridge approximations are computed using first degree polynomials. In Table 1, results with higher degree polynomials are presented. Clearly, increasing the polynomial degree, p, allows for lower test errors to be obtained (although excessively high-degree polynomials will suffer from over-fitting). However, increasing p also increases the cardinality of the polynomial; if all interaction terms are included, and the polynomial has degree p in all n directions, its cardinality is given by L = (p + 1) n . More degrees of freedom result in higher variances of the resulting regression estimators, necessitating more training samples to mitigate this. The exponential dependence on polynomial degree implies that higher degree polynomials are significantly more costly when in higher dimensions. This highlights a key benefit of obtaining polynomials in a dimension reducing subspace of dimension n, where n d. Various approaches exist to mitigate this exponential scaling [36]. However, even then, it is still desirable to reduce the number of dimensions. As a compromise between accuracy and computational cost, second degree polynomials are used in the remainder of the results presented for the present test case.

Design Space Exploration
The dimension reducing nature of the ridge approximations is particularly useful for design space exploration. To demonstrate, ridge approximations for the turbulent viscosity ratio ν t /ν ∞ are obtained using all 500 training designs. In Figure 9 a sufficient summary plot for a point in the aerofoil's wake is shown. The turbulent viscosity ratio in the wake is of interest because high levels of turbulence in wakes can be a leading source of loss in many flows [51]. The sufficient summary plot is examined, and a design at either end of the W T x subspace is selected. Contours of the turbulent viscosity ratio for the two designs are shown in Figure 9. Through this approach, the end user can explore how flowfield variables at select locations in the flowfield change as they traverse the design space. Such an activity would be far more challenging when negotiating the original design space x ∈ R d . Effectively, dimension reduction distils the design space into critical directions along which quantities of interest exhibit maximum variance, which enables efficient design space exploration.  In addition to exploring directions of maximum variance, one can also explore directions along which quantities of interest exhibit minimum variance. Following Section 2.2, the input design vector can be decomposed as where w = W T x is the active coordinate and v = V T x the inactive coordinate. The inactive subspace V permits identification of designs which are invariant with regards to the quantity of interest. The task is to find input vectors x with a fixed w but different v, whilst obeying the constraint −1 ≤ x ≤ 1, for which we use a hit-and-run algorithm similar to that implemented in [52]. In Figure 10, this is done for the ridge approximations at the locations labelled a) and b) in Figure 6. The resulting design envelopes inform us how the aerofoil can be deformed without the quantity of interest u/u ∞ at these two locations being affected. Such information can be used to form manufacturing tolerances, or to decide whether to scrap components, as is done in [23,53,54].  Figure 6. 1000 samples are generated using a hit-and-run algorithm [52], with the grey lines corresponding to the blade profiles of 10 randomly selected samples.

Physical Insight
In addition to exploring the design space, ridge functions can also provide important physical insights. As an example, we consider the stagnation pressure around the aerofoil. In an ideal (inviscid) incompressible flow, the stagnation pressure, p 0 = p + 0.5ρU 2 , is constant along streamlines 4 . Streamline curvature results in a reversible exchange between the static pressure p and dynamic pressure 0.5ρU 2 . However, in a real flow, irreversible processes such as viscous effects cause a loss in p 0 , which it is important to control. Take the case of an aircraft; high losses over the wing are synonymous with high drag [55], while in the engines, high losses lead to decreased component efficiencies, both of which lead to a higher fuel consumption. The pressure loss at a given location s can be quantified through the loss coefficient where p 0∞ is the freestream stagnation pressure, and Y p = 0 indicates no loss in stagnation pressure, while Y p = 1 indicates a total loss in stagnation pressure. Figure 11: Contours of loss coefficient, Y p , for a deformed aerofoil from the test set ("Design 1") at an angle of incidence of α ∞ = 10 • . A line at 50% axial chord downstream of the trailing edge is highlighted, and sufficient summary plots are shown for two points within the wake, on this line. Design 1, and a second design, are highlighted in the summary plots.
The loss coefficient Y p is calculated from the CFD data for 500 training designs, and the approximations for a randomly selected test design are visualised in Figure 11. In Figure 12a, profiles of Y p across the s 1 /C 1 = 1.5 plane are plotted for the two designs labelled 1 and 2 in Figure 11. For both designs, the high loss in the wake (around 0 ≤ s 2 /C 1 ≤ 0.15) is accurately predicted by the ridge functions.
The sufficient summary plots in Figure 11 show how Y p at probe locations A and B change as the design space is traversed. Furthermore, the dimension reducing subspace vector W can be examined. Recalling that in the present framework x contains the Hicks-Henne bump amplitudes, the elements of W show us how the bump amplitudes affect the reduced coordinates W T x. In Figure 12a, the scaled elements of W for the Y p ridge functions at A and B are projected onto the baseline aerofoil geometry 5 . The weights offer important physical insights here; they show us how the aerofoil surfaces must be deformed in order to affect the loss coefficient Y p at probes A and B.

Obtaining Global Quantities
Understanding the sensitivity of local flow variables to a design can be invaluable to a designer. But, sometimes, a designer may also be interested in global quantities of interest. For example, in the present case, the mass averaged loss across the s 1 /C 1 = 1.5 line in Figure 11 is an important quantity, and is given bȳ where u(s 1 , s 2 ) and Y p (s 1 , s 2 ) denote the local axial velocity and loss coefficient respectively. The mass averaged loss represents the total pressure loss caused by the aerofoil. A surrogate model forȲ p could be obtained by preprocessing the CFD data, and then computing a ridge approximation for the precomputedȲ p values. However, alternatively, the local ridge functions themselves can be integrated following the approach outlined in Appendix B. A dimension 5 Each element of the subspace weight vector is plotted at the surface location of its associated bump. a) b) Figure 12: a) Profiles of loss coefficient Y p across s 1 /C 1 = 1.5, and their ridge approximations, for the two designs highlighted in Figure 11. Also shown are the s 2 locations of probes A and B from Figure 11 (dashed lines), along with the projections of W at these two locations onto the baseline aerofoil. b) Sufficient summary plot for the mass-averaged loss coefficient across the s 1 /C 1 = 1.5 line, with designs 1 and 2 highlighted, and the associated subspace vector W int projected onto the baseline aerofoil.
reducing subspace forȲ p , W int , is obtained directly from the ridge approximations for u(s 1 , s 2 ) and Y p (s 1 , s 2 ). The resulting one dimensional sufficient summary plot forȲ p is shown in Figure 12b. It is apparent that designs 1 and 2 lie at either end of the design space with regards to their mass averaged loss, indicating that design 2 is objectively the better design (with regards to loss). Similarly to before, in Figure 12b the elements of W int are projected onto the baseline aerofoil. This shows that the total loss can be reduced by deforming the suction surface outwards near the leading edge, and inwards near to the trailing edge.
To understand why design 2 has a lower loss, a second law analysis can be performed. Herwig and Schmandt [56] show that the overall entropy generation rate can be used to understand the sources of drag or pressure loss within a flow. The contours of S in Figure 13 show that, at α = 10 • , the majority of the entropy generation comes from the suction surface boundary layer and the wake region 6 . In Figure 13 the delta between the y−averaged entropy generation for design 2 and the baseline design (W T x = 0) is shown. Deforming the suction surface outwards near the leading edge, and flattening the aft portion (see Fig. 12b), appears to increase entropy generation over the first quarter of the suction surface. But, this is offset by a reduction in entropy generation in the wake and leading edge stagnation region, leading to a net decrease in the entropy generation. Whilst an experienced aerodynamicist may have foreseen this conclusion, the present example serves to illustrate the potential of an easily interpretable surrogate model which can be readily integrated to obtain global quantities.

Exploiting Spatial Correlation
Up to this point, flowfield predictions have been made using ridge approximations at every grid point. For larger grids, such as three dimensional ones, obtaining N number of ridge approximations might be computationally intractable. Taking a similar approach to how the CNN takes advantage of spatial coherence in the flowfield, the method outlined in Section 2.4 explots spatial correlations to reduce the number of ridge approximations required. In Figure 14, ridge approximations are obtained for J = 1000 randomly subsampled points (a sampling rate of J/N = 3.5%). Since the baseline mesh has a high mesh density near the aerofoil (where the gradients are high), the subsampled points are also clustered near to the aerofoil, which is desirable. The covariance matrix K, computed on the training data, is then used to upsample the ridge approximations at the subsampled points back to the remaining N − J points.
The colouring of the subsampled points in Figure 14 shows that there are a small number of points with poor quality ridge approximations. Preliminary work [40] showed that these poor quality approximations can cause spurious noise in the flowfield estimates when upsampled. As discussed in Section 2.4, the training R 2 score can be used as a heuristic to  measure the quality of each ridge. Ridges with R 2 < 0.6 are ignored during upsampling, with their grid points added back into the N −J set of points. As seen in Figure 14, this strategy results in qualitatively accurate flowfield predictions using ridge approximations at only 3.5% of the grid points, reducing the training time by a factor of approximately 29.

Test case 2: Varying Freestream Conditions for a Transonic Wing
Following on from the success of the first test case, we now turn to the transonic wing test case described in Section 3.1.2 in order to expose the limitations of polynomial ridge functions. Here the geometry, represented by x ∈ R b , is fixed.
However, the boundary conditions, represented by v ∈ R b , are varied. The transonic nature of the flow, and the possibility of extrapolating beyond the training dataset, make this test case particularly challenging. The input vector v only consists of two dimensions in this case, M a ∞ and α ∞ . We could once again obtain one dimensional ridge functions. However, two dimensional ridge functions are still simple to visualise. Considering a ridge function of the form we set W = I (the identity matrix), resulting in two dimensional ridge functions (n = b).

Interpolation versus extrapolation
In this case, the behaviour of the ridge functions as we extrapolate away from the training boundary conditions is important. To explore this, polynomial ridge functions are trained on the M train = 36 point DoE shown in the inset of Moving on to DoE point 2, here the ridge approximations must extrapolate away from the flow conditions seen during training. Unsurprisingly, the predicted distribution in Figure 15b agrees less well with the true distribution in this case, especially on the aft portion of the wing. This highlights a note of caution regarding polynomial ridge functions; as with many function approximation techniques, caution must be taken when using the underlying polynomials to extrapolate beyond the training data. If such a task is routinely necessary, it is possible to replace the polynomials with other classes of models with more stable behaviour under extrapolation, such as piecewise polynomials (splines) or kernel regressors. However, we shall leave this as a topic of future work.

Three-dimensional flowfield predictions
When extending the framework to three dimensional flowfields, computational efficiency is crucial. To achieve a surrogate model for the full flowfield the subsampling and upsampling strategy is used once again. Ridge functions are fitted at the J = 5000 randomly sampled points (a sampling rate of 4.6%), with a covariance matrix used to upsample the ridge approximations to full three dimensional flowfield estimates. A slice of a flowfield estimate for a test DoE point is shown in Figure 16. Comparing the colour contours and iso-lines, the prediction appears to be reasonable, with only a slight disagreement on the aft section of the suction surface. The above result suggests the subsampling strategy is effective for three dimensional flows in addition to two dimensional ones. The inset figure within Figure 16 shows the effect of sampling rate on the prediction error. For the present case, sampling rates of lower than ≈ 4% lead to significantly increased errors. Nevertheless, sampling rates of 4-10% still allow for significant computational cost savings.

Conclusions
This paper exposes the idea of using spatially correlated ridge functions to rapidly estimate flowfields. The resulting data-driven framework could be trained on existing CFD data where available, or it could be integrated within a wider design of experiment, with new CFD data generated specifically for training. Once trained, the framework provides rapid flowfield predictions, which can be used for design space exploration, design optimisation tasks, performance predictions, or even to guide CFD mesh design.
On the subsonic aerofoil test case, comparisons with the state-of-the-art convolutional neural network (CNN) suggest that the ridge function framework is able to achieve competitive predictive accuracy on unseen aerofoil designs. In addition to serving as a surrogate model, the learned ridge functions can aid understanding of the existing training data. Sufficient summary plots can be viewed for any given point. The reduced dimensional nature of these plots lends itself to easy visualisation, allowing for easy comparison between designs, and new physical insights. Furthermore, the local ridge functions can be integrated to obtain new ridge functions for integral quantities such as a loss or drag coefficient. This allows new interpretable surrogate models for integral quantities to be generated on demand, without the need for further pre-processing of the CFD data. Such capability is advantageous in scenarios where the long-term storage of the full CFD dataset is problematic. The ridge functions can also provide insensitivity information, informing us how a design can be altered without influencing the flow at a given location.
The independent nature of each ridge function means that their training can be viewed as an embarrassingly parallel task. This makes the ridge function framework trivial to implement in a parallel fashion, allowing for excellent scaling with problem size. However, for larger problems, computing ridge functions for every single grid point is undesirable. Instead, ridge functions are obtained at a much smaller number of randomly subsampled points. The flow physics encoded within covariance matrices, computed from the training data, can be used to upsample the ridge functions' predictions back to the rest of the flowfield. The computation of the covariance matrices, and the subsequent use of them, involves a variety of linear algebra operations which can be implemented efficiently. For the three dimensional transonic wing test case, sampling rates around 5% are achieved with minimal additional errors, allowing for significant reductions in training time.
The transonic wing test case demonstrates how ridge functions can be used to estimate flowfields with varying boundary conditions. In such cases, there is a danger of extrapolating too far from the training data, and the underlying polynomials' predictions can not be trusted here. If extrapolation is important, there is scope to replace the polynomials used in this work with alternative models. In addition to exploring this further, a worthwhile area of future work would be to combine the approaches used in the two test cases, incorporating geometric and boundary condition input parameters together.
grid points. A learning rate of 4 × 10 −4 is used, which is then reduced by a factor of ten after 500 epochs, and a batch size of ten is chosen. The learning curve for the 488k parameter network is plotted in Figure 17; the convergence behaviour appears to be satisfactory, and the small difference between train and test loss suggests the network isn't suffering from excessive overfitting. Similar behaviour is observed for the other network sizes explored in this paper.

A.2 Effect of network size
To allow for a fair comparison between the ridge function framework and the CNN, the performance for various size/complexity networks is examined in order to find an optimal size. The network's size is altered by changing the number of channels, with the number of channels in the i th layer given by C i = 2 e c i . The constants c i are defined in Table 2, and the network's size is varied by adjusting the exponent e, with values of e = {2, 3, 4, 5} tested. The variation in performace The normalised MAE, defined in (17), is measured on the training and test sets for the four different network sizes, with the results plotted in Figure 18. Generally, the larger networks achieve lower predictive errors. However, this is at the cost of longer training times (Fig. 18d), with the 1.94M parameter network taking up to 80 minutes to train. The larger networks also exhibit more stable predictions, with less variance in their MAE scores across the three train/test folds. To achieve a reasonable compromise between accuracy and training cost, the 488k parameter is chosen as the benchmark for comparisons elsewhere in this paper. A similar pattern between MAE, training times, and the number of training designs is also observed. The chosen value of M train = 500 is seen to offer a good compromise between low MAE scores and short training time.

B Integrating Ridge Functions
This appendix presents a procedure for identifying a dimension reducing subspace for the mass-averaged pressure loss coefficientȲ p , defined in (20). A more general treatment on the subspaces of an integral of a scalar field quantity can be found in Wong et al. [24]. For simplicity, we stick to a two dimensional domain, but the approach is readily generalisable to three dimensions. Recalling (20), we wish to obtain the mass-averaged loss coefficient Y p , at a given where the mass flow rate across the plane,ṁ = C1 u(s 1 , s 2 )ds 2 , is taken to be a constant here 7 . The local axial velocity and loss coefficient are approximated with two ridge approximations u(s 1 , s 2 ) ≈ g s1,s2 W T s1,s2 x Y p (s 1 , s 2 ) ≈ h x,y U T s1,s2 x .
Here the input vector x ∈ R d has its usual meaning, whilst W T ∈ R nw×d and U T ∈ R nu×d denote the subspaces for axial velocity and loss coefficient respectively, where naturally n w d and n u d. Substituting (24) into (23) yields To identify the dimension reducing subspace forȲ p with respect to its constituent input parameters x, we need to compute the averaged outer product of its gradient [39] -a covariance matrix given by where ω(x) denotes the distribution associated with the input parameters x. Assuming the existence of a design of experiment with M distinct CFD evaluations for different input parametrisations, and taking ω(x) to be a uniform distribution here, we can approximate C with its finite sample estimate The gradient evaluations can be written as [W s1,s2 ∇ w g s1,s2 (w m ) h s1,s2 (u m ) + U s1,s2 ∇ u h s1,s2 (u m ) g s1,s2 (w m )] ds 2 , where w m = W T x m and u m = U T x m are the reduced dimensional co-ordinates for the m th design in the Y p (s 1 , s 2 ) and u(s 1 , s 2 ) subspaces respectively. If we further assume that spatially, there are N nodes along which this integration has to be performed, in which case we have the set we approximate where ϕ i denotes appropriately selected quadrature weights depending on the spatial location of the nodal centres. Plugging this into (27) yields This equation is easy to implement, as long as the gradients ∇ w g and ∇ u h can be computed. When applying this approach in Section 4.3, g and h are taken to be the underlying polynomials of the ridge function approximations, rather than the posterior mean of the ridge kernels as is the case in the rest of the paper. This is done since the gradients of orthogonal polynomials are readily available. Finally, an eigendecomposition of C is performed, and the leading eigenvectors with the largest eigenvalues are chosen to form W int , the dimension reducing subspace forȲ p . In the present case, the first eigenvalue was found to contribute 99.3% of the trace. Therefore, a one dimensional dimension reducing subspace was deemed to be sufficient.