Graph Convolutional Multi-Mesh Autoencoder for Steady Transonic Aircraft Aerodynamics

. Calculating aerodynamic loads around an aircraft using computational fluid dynamics is a user’s and computer-intensive task. An attractive alternative is to leverage neural networks bypassing the need of solving the governing fluid equations at all flight conditions of interest. Neural networks have the ability to infer highly nonlinear predictions if a reference dataset is available. This work presents a geometric deep learning based multi-mesh autoencoder framework for steady-state transonic aerodynamics. The framework builds on graph neural networks which are designed for irregular and unstructured spatial discretisations, embedded in a multi-resolution algorithm for dimensionality reduction. The test case is for the NASA Common Research Model wing/body aircraft configuration. Thorough studies are presented discussing the model predictions in terms of vector fields, pressure and shear-stress coefficients, and scalar fields, total force and moment coefficients, for a range of nonlinear conditions involving shock waves and flow separation. We note that the cost of the model prediction is minimal having used an existing database.


Introduction
Aerodynamic analyses in many engineering sectors, from aerospace to automotive, remain computationally expensive with high-fidelity CFD.The need to provide physical insights at very small scales around a complex geometry with limited resources contrasts with the ever pressing project deadlines.This motivates the investigation of ROMs for accurate approximations of the physical system via numerical techniques of reduced computational complexity.The recent breakthrough in ML has prompted the development of new methodologies particularly suited for ROMs.Modern, sophisticated ML algorithms are attractive for approximating highly complex and nonlinear systems from data.For example, Massegur et al [1] leveraged ML to predict the formation of ice on a wing section and the resulting aerodynamic degradation across a range of freezing conditions.Furthermore, Massegur et al [2] developed an ML based aerodynamic model coupled with a structural model for aeroelastic and flutter search analyses.Herein, we focus on three-dimensional (3D) geometries and we address the added complexities compared to two-dimensional (2D) cases, which are around ten times more frequent in the literature.Regarding the physics, 3D problems present more abundant nonlinear flow features and interactions [3,4], which include cross-flow interactions, wing-tip vortices and separation bubbles, among others.ROM approaches developed for 2D problems [2], featuring direct prediction of aerodynamic forces, are questionable to address the additional intricacies that arise in 3D problems.To this goal, we use deep learning NNs [5] to model distributed aerodynamic quantities on the aircraft surface.
The primary challenge in reduced-order modelling of 3D aerodynamic fields is the large dimensional space, reflecting the numerical discretisation of the computational domain [6].Common nonlinear ROM approaches, including Kriging [7] and DNNs [8,9] result ill-suited.To address the problem of poor scalability, CNNs are better suited for flow-field analyses because they are designed to extract spatial features from digital images [10][11][12].The problem is that CNNs are not directly suitable for unstructured meshes typical of CFD applications because of their limitation to Euclidean domains (Cartesian grids).Interpolating data to a regular grid is essential, but it can lead to additional errors and loss of resolution in refined regions.Alternatively, CNN kernels of size 1 can be used as these are applicable to any mesh arrangement, as implemented by Immordino et al [13] or Sabater et al [14].The issue with this approach is that the mesh connectivity is unused, resulting in simpler cloud-point modelling.The absence of propagation of information across the mesh limits the capability to capture coherent flow features and reconstruct continuous solutions.To leverage connectivity throughout the mesh, geometric deep learning [15] offers a range of convolutional methods, known as GNNs.GNNs are capable of processing graphs or manifold-unstructured meshes, which makes them a better fit for CFD applications that deal with non-Euclidean data [16].Examples of GNNs for steady-state aerodynamics are not abundant.Ogoke et al [17] addressed aerodynamic solution of a 2D aerofoil rather than a more relevant test case.Baqué et al [18] projected the 3D geometry to a simpler prismatic graph to contain the computational burden.On the other hand, comprehensive message-passing methods, as adopted by Hines and Bekemeyer [19] or Han et al [20], which feature dedicated learnable weights for each edge and node of the mesh, can lead to excessive memory requirements.
To maximise ROM computational efficiency when dealing with large spatial domains, dimensionality reduction to compress the domain size is crucial.The classical POD [21,22] is only limited to linear projections of the data onto the compressed space, causing general loss of information in highly complex physics.AEs [23][24][25][26], which are the NN alternative to the POD, are a better alternative allowing a nonlinear compression and recovery.This work contributes generally to the question of dimensionality reduction of large datasets defined on irregular spatial domains.To this aim, an AE approach embedded in the graph-based convolutional framework is sought.We solved this challenge by devising a multi-resolution scheme, inspired by the common multi-grid methods to solve partial differential equations [27].While the methodology applies naturally to other fields and disciplines, we demonstrate the applicability on a 3D aerodynamic problem.
We propose a graph-convolutional MM AE tailored for predicting distributed aerodynamic loads around a full-scale air vehicle.Unique contributions of this work are the applicability to non-Euclidean domains and a unique multi-resolution embedding for dimensionality reduction and modelling efficiency purposes.Furthermore, a building-block implementation, consisting of composition of separate NN units, facilitates evaluation of different model architectures to address the task.We have chosen the wing/body configuration of the NASA CRM for demonstration, predicting steady-state transonic aerodynamic loads across a range of Mach numbers and angles of attack.
The paper continues in section 2 with a description of the problem we solved and the identification of nonlinear features to be retained in the model outputs.The methodology for steady-state problems is explained in section 3.Then, section 4 summarises the main results.Finally, conclusions and an outlook on future work are given in section 5.

NASA CRM wing/body configuration
The test case is for the NASA CRM wing/body configuration, which is representative of a transonic transport aircraft designed to fly at a cruise Mach number of 0.8 [28,29].The CRM is shown in figure 1.The reference geometric chord is c = 0.1412 m, the surface area is S = 0.0727 m 2 , and the pitch moment is taken around x a /c = 0.5049.
The CFD dataset was taken from Immordino et al [13], where SU2 7.5 solver was used.The surface mesh shown in figure 1 features around 78 000 grid points in an unstructured triangular topology with adapted discretisation density around the edges.The volume mesh consists of over 15 million cells, with a prism-layer mesh to promote y + < 1 near the wall, denoting the non-dimensional height of the first cell normal to the boundary.The y + ensures an adequate resolution of the boundary layer.The Spalart-Allmaras turbulence model was chosen.For good convergence, the JST scheme with added dissipation was adopted to discretise the convective term, Green Gauss was chosen to compute the discretised gradients and the biconjugate gradient with ILU preconditioner was applied to solve the linear solver.
In this work, we are interested in the steady-state prediction of the pressure coefficient: with p ∞ , ρ ∞ and U ∞ the freestream pressure, density and velocity, respectively.The shear-stress coefficient components: with τ w = [τ x , τ y , τ z ], on the surface mesh of the CFD model across an envelope of Mach numbers M ∞ and angles of attack α ∞ .

Reference dataset
We used the available database of aerodynamic solutions to generate and validate our steady-state predictive model.The database contains 70 CFD calculations in the range of M ∞ ∈ [0.70, 0.84] and α ∞ ∈ [0.0, 5.0] deg, at Reynolds number Re = 5 • 10 6 and freestream temperature T ∞ = 311 K.The range of freestream conditions was chosen to have diverse nonlinear flow phenomena to be predicted by our model.The sampled conditions are shown in figure 2, of which 40 were used for training (represented with circles) and the remaining 30 for validation (triangles).Latin hypercube technique was used to define this limited number of samples for the preliminary envelope scan.The resulting sampling distribution did not include points on the boundaries and occasionally left large gaps between samples.Therefore, these experiments are useful to assess the performance of our framework in under-sampled spaces.First, we delve into the variation of lift and drag coefficients, C L and C D , respectively, with the freestream conditions shown in figure 3.These were obtained by integrating the pressure coefficient C P and shear-stress C τ fields of the reference CFD solutions.Inspecting the colormaps, the lift coefficient correlates predominantly with the angle of attack.The variation of the drag coefficient indicates an equal correlation to α ∞ and M ∞ .Both the angle of attack and the shock-induced boundary-layer separation contribute to the drag increase.
Then, we extracted few sample points from the database for further demonstration of the complexity of our problem.Figure 4 shows the pressure coefficient for selected conditions, and figure 5 is for the skin friction coefficient, defined as the shear-stress norm C f = ||C τ || 2 .These figures are arranged into four panels, placed to reflect the location of the samples labelled in figure 2, with low α ∞ on the bottom panels and high M ∞ on the right panels.The pressure distribution, figure 4, significantly differs depending on the freestream conditions.Inspecting the lower and the upper left panels, the location of the shock wave moves towards the leading edge and the shock intensity increases with angle of attack at lower Mach numbers.This transition  appears to be nonlinear for α ∞ higher than 3 deg.With increasing Mach number, the pressure distribution becomes gradually smoother and the shock wave becomes stronger.As a result, at the highest M ∞ , right panels, the peak C P values are lower.On the contrary, at these M ∞ conditions, the location of the shock wave remains similar with increasing angle of attack, bottom to top of the right panels.
Similarly, the skin friction coefficient distributions, in figure 5, indicate the boundary-layer separation regions (darker blue) induced by the shock wave.At low Mach number, left panels, the separation line moves towards the leading edge with increasing angle of attack.Furthermore, in panel (a) a separation bubble is visible.On the contrary, above M ∞ = 0.8, the location of the separation line remains similar across the α ∞ range, as seen in the right panels.
This brief overview sets the background problem for our ROM.The diverse nonlinear flow characteristics cannot be captured by simple direct modelling of the scalar loads, such as C L and C D .On the other hand, the task of predicting distributed quantities at each condition, such as C P , is more challenging than scalar targets but is more useful from a design standpoint.This motivates our choice to devise a GDL based framework for prediction of the surface aerodynamic fields across the operating envelope.

Methodology
In a steady-state formulation, the aerodynamic response is considered dependent on the input conditions only, and any time dependence is neglected.An NN function f NN is sought which maps specific user-defined inputs s to desired target fields Y i on the surface mesh S: with i denoting a node in S. The grid point coordinates x i are also embedded.We devise an architecture for f NN by leveraging GDL and dimensionality reduction for unstructurally meshed manifolds.From GDL, we resort to GNNs, which involve the convolution operation on graphs [16,30].

Graph representation
In a GNN approach, the surface mesh is represented as a graph where the vertices (or nodes) contain the position coordinates x i and variable fields y i (features).The graph edges connecting the grid points are determined by the mesh connectivity.Figure 6 illustrates a representation of a graph where a target node i is connected to j ∈ S surrounding grid points.Features y i and weights e ij are assigned, respectively, to each node and edge.The edges defining the graph connectivity and chosen weights are arranged to form the adjacency matrix [16]: This is an n × n matrix containing the edge weights, and n is the total number of grid points in S. The subscript ij denotes the jth source node connected to a given node i.We assign the weights as the inverse of the distance between the two nodes forming the edge: We then normalise the edge weights to be ∈ (0, 1], where the upper end is inclusive because self loops, i.e. e ii = 1, are inserted by adding the identity matrix to the adjacency matrix: Â = A + I.In addition, we choose the edge weights to be non-directional, i.e. e ij = e ji , which results in a symmetric adjacency matrix, Â = ÂT .Note that Â is largely sparse as each row contains only a few non-zero elements.Consequently, it is more memory-efficient to arrange the adjacency matrix in COO format.This format consists of two vectors: the edge-index and the edge-weight vectors.The edge-index vector contains the pair of node indices [i, j] for each edge, of size n e × 2 and n e the number of edges.The edge-weight vector contains the assigned weight edges e ij equation ( 5), with size n e × 1.In the CRM test case, the surface mesh consists of n = 78 829 points and n e = 472 404 edges.

GCN
From the family of GNN architectures, we leverage the GCN by Kipf and Welling [31].The GCN operator at a given target node is defined as: with θ a layer-specific trainable weight vector, b a constant term and y the node-based input vector at each node of the mesh S.
) ∀ i contains the sum of the edge weights connected to each node i, known as the diagonal degree matrix.
At each layer l, the GCN operation, equation ( 6), is executed on the output from the previous layer y l−1 , followed by a nonlinear activation function h: We adopted for h the PReLU [32]: with β another learnable parameter.Note how the GCN operator equation ( 6) takes the convolutional analogy of the CNN for Euclidean domains [33] but with a single-parameter filter swept across each row of the adjacency matrix.In fact, if we set set e ij = 0 for i ̸ = j, the standard CNN layer with kernel size 1 is obtained, as Immordino et al [13] opted.

GCN-based AE
From equation (6) we realise that successive GCN executions are required to exert influence between far-away grid points.For instance, with k l GCN layers in succession, only k l neighbourhoods around target node i are influenced.Consequently, the propagation of information across the mesh is slow.This leads to two issues.The first is that a deep network with an excessive number of layers would be necessary to propagate the information in refined regions.The second issue relates to the memory size of the network becoming computationally unmanageable in large spatial domains and high number of features.To alleviate these issues, as introduced in section 1, we adopted an AE approach for the compression of the spatial domain.
The AE involves the projection of the states from the original domain onto a compressed space, operation known as encoder.Then, these latent states are recovered back onto the original domain in an inverse operation, known as decoder.Embedding NNs in this process makes the AEs more attractive than POD as nonlinear projections are possible [23].Figure 7 lays out the AE process embedded with GCNs.

MM scheme
AEs for discretised domains entail a multi-resolution scheme, which involves gradual coarsening operations and a subsequent refining of the grid.Reduction techniques in Cartesian arrangements are trivial, including, for example, pooling operations [11].In contrast, coarsening of unstructured meshes is a more difficult task.One of the primary challenges is that the adjacency matrix needs to be regenerated at every coarsening step.A second challenge relates to reliable transfer of information between the various grid resolutions.
To solve these challenges, we present a novel hierarchical MM scheme for the AE.In the encoder process, the mesh is coarsened between blocks of GCN layers, intended for extraction and compression of crucial flow-field features.The latent states on the coarsened mesh are decoded in a recovery operation interleaved with additional GCN blocks, to reconstruct the solution onto the original domain.Figure 8 illustrates the coarsening and recovery operations of the proposed 2-level multi-mesh cycle for the CRM model.This method is reminiscent of the common multi-grid V-cycle algorithms to solve partial differential  equations [27,34].Operating on an MM cycle is advantageous for: (1) reducing the computing memory size, given the compressed spatial domain; (2) extracting features of different spatial scales by means of the different mesh resolutions; and (3) enabling direct information exchange among distant nodes, avoiding the need for a deep network to spread influence across the grid, which results in a significantly smaller model.
The coarsening operation in the encoder involves the removal of grid points from the original mesh.The strategy taken to coarsen the mesh is crucial to prevent loss of essential information.For instance, a uniform random selection should preserve the original mesh topology, in terms of relative cell sizes, on the reduced mesh.However, there is risk of insufficient resolution left in regions where the initial node density was already low.On the other hand, excessive removal of nodes on originally refined regions could lead to inappropriate reconstructions where the solution is likely to present larger gradients.
For adequate representation of the distinct spatial regions at the coarse level, a balanced node selection is essential.We achieved this by selecting the nodes according to a probability function based on the corresponding face area: with i the index of the grid points sorted by their face area a i in ascending order, and n the total node count.We chose p 1 = 0.2 and p n = 1 for the smallest and largest elements, respectively.The resulting distribution is demonstrated in figure 9. Probability of being selected is higher in nodes with larger face areas, likely in  unrefined regions, as opposed to nodes found in dense discretisations.The coarsened mesh resulting from this selection probability is illustrated at the bottom of figure 8.The mesh size (node count) in the various multi-resolution levels is reported in table 1 for the CRM test case.Note how a compression ratio of 16 was adopted.
Upon selection of the grid points to be kept in the coarse mesh, the graph connectivity was regenerated by reconnecting remaining nodes that shared connections with discarded ones.Figure 10 demonstrates the process of restoring graph connectivity following a coarsening operation, with new edges shown in orange.

Weighted moving least squares for grid interpolation
As shown in figure 8, information must be transferred across multiple grids.In the reduction step of the encoder, the original field data must be interpolated onto the compressed mesh.In the decoder, the recovery of the latent states onto the fine grid entails an inverse interpolation.These are critical operations in the MM cycle.The interpolation consists of a functional I Sn→Sr : R nn → R nr to map a spatial field y i from a source grid S n , which contains n n nodes, to a target grid S r , with n r nodes, where both grids discretise the same spatial domain: with I Sn→Sr the interpolation matrix from the source to the target mesh, and y j the interpolated field data.
The following properties are desirable for an adequate interpolation [35]: (1) interpolated values at the source nodes should match the original data; (2) integrated resultants should be conserved; and (3) interpolated fields should be continuous.Consequently, directly recasting the data across coincident points and nearest-neighbour interpolation, as Han et al [20] proposed, result inappropriate because conservation and continuity properties are not satisfied.
There is multitude of multi-grid algorithms, often devised to accelerating the solution of finite-volume discretisation of partial differential equations, as proposed for example by Smith [34].However, we chose a different approach which is particularly suited for fluid-structure interaction problems, where satisfying the above properties is crucial for adequate transfer of loads and deflections across models.In particular, we implemented a weighted moving least squares (WMLS) scheme [35,36].The idea is to generate a shape function u(x) to approximate the input data y i at source nodes i ∈ S n with coordinates x i by least-square-error minimisation: where the weight w(||x − x i ||) is a function of the distance between the source and the target points.We specify u(x) as a polynomial combination: with p(x) = [1, x, y, z, x 2 , y 2 , z 2 , xy, xz, yz] T the second-order polynomial basis function, and a the vector of respective coefficients.The analytical solution of the least square minimisation can be shown to yield the resulting approximation at every node j ∈ S r of the target grid: with the shape functions defined as: and Therefore, the interpolation matrix I Sn→Sr is: The least squares approximation, equation ( 14), must be computed for each target grid point, resulting in computationally intractable matrix operations when dealing with large meshes.To reduce the computing burden, we adopt a moving interpolation consisting of limiting each target node to be influenced only by the k n closest source points: In this work, we found k n = 10 a good compromise between interpolation accuracy and computational efficiency.Note that I Sn→Sr is a largely sparse and non-square matrix of size n r × n n , each row containing just k n non-zero values.This matrix is not invertible and two different interpolation matrices must be generated for the encoder and decoder operations, I Sn→Sr and I Sr→Sn , respectively.Figure 8 illustrates the dual interpolation process.We observe how the resulting pressure reconstruction (right) after execution of the MM cycle matches well the original field (left).

Steady-state prediction framework
We now complete the construction of the predictive model architecture, here referred to as steady-state GCN-MM-AE.For the CRM use case, the target fields are the pressure coefficient and the shear stress components, Y = [C P , C τx , C τy , C τz ], across the input envelope of Mach numbers and angles of attack, The architecture of the final steady-state GCN-MM-AE model is shown in figure 11.The scalar inputs s are cast to each node of the graph, concatenated to the grid-point coordinates x i ∀ i ∈ S n .The input vectors are processed by the encoder, involving the coarsening step of the MM cycle and two GCN blocks.Subsequently, the decoder comprises the recovery operation of the MM embedded in two additional GCN blocks.The network ramifies at the end into separate blocks for each field quantity to predict.This architecture defines f NN in equation (3).
To the best of our knowledge, this framework is novel on several fronts: (1) GDL based AE framework for spatial predictions on large and unstructured discretisations, applied to an aerospace problem; (2) multi-resolution scheme embedded in the nonlinear AE for dimensionality reduction of unstructured manifolds, aimed at capturing different spatial scales, promoting influence across the grid and maximising ROM computational efficiency; and (3) building-block functionality to address distinct tasks within the same framework: multi-resolution reconstructions, steady-state predictions and extension to dynamic simulations.This framework was implemented using PyTorch 1.11, an optimised deep-learning library in Python, and PyTorch Geometric, an open-source GNN package built upon PyTorch.

Results
This section is organised in two parts.The first part focuses on the prediction of scalar quantities-in our case, the integrated force and moment coefficients.The second part is related to the model prediction of the distributed fields-the pressure and the shear-stress coefficient distributions.The appendix contains more background information, including the steady-state GCN-MM-AE architecture from figure 11 and the model optimisation procedure.Comparison between our WMLS scheme and the multi-grid method by Smith [34] for the two-way interpolation of the MM cycle was also investigated.To complete the framework set-up, sensitivity assessments to key hyperparameters, such as the training set size, the model weights and the MM cycle set-up, are also reported.Worth noting that at each sample point the model prediction outputs C P and C τ distributions, from which the resulting force and moment coefficients were obtained by integration.

Integrated loads
Figure 12 presents the analysis for the integrated lift coefficient C L , drag coefficient C D and pitching moment coefficient C My .The left panels show the prediction error on the dataset samples, with training samples in circles and validation in triangles.The error is defined as: with y = [L, D, M y ].The errors are classified by a traffic-light colour scheme: green corresponds to prediction errors below 4%, amber between 4% and 10%, and red above 10%.For convenience, the percent error is reported above each sample point.Best and worst predictions are highlighted as A and B, respectively.The prediction error is generally small throughout for the C L predictions, where the error is below 2.2%.Additionally, the worst prediction (B) was found at M ∞ = 0.76, α ∞ = 4.42 deg, with pitching moment error of 24.3%.This point stands out for not including training samples within a wide surrounding.The model found more difficult learning this region of the envelope.An adaptive sampling method to include training samples in under-sampled regions would be convenient for improved model accuracy.However, this is beyond the scope of this work.We also found the C My prediction accuracy degrades slightly towards high angles of attack, where nonlinear aerodynamic response occurs.Nevertheless, the accuracy of the model is overall high.Table 2 provides a statistical summary of the prediction errors.The average error, standard deviation and the worst prediction for each aerodynamic coefficient across the complete dataset and the validation dataset are reported.The statistics are similar between datasets.This suggests that there is no over-fitting and the model performs well to new conditions.In addition, the low standard deviations indicate good model precision.The average errors were also found low, with C L the best predicted quantity.The errors for the C D and C My tend to be skewed by the reference values being an order of magnitude smaller than the lift values.Last, the worst statistics are for the C My , consequence of the larger errors at high angles of attack.Small errors on the shock-wave location around the reference axis can contribute to a magnified error too.
The adoption of the traffic-light system for the error plots in figure 12 is convenient to judge the adequacy of the model for design purposes.The low error obtained for the lift is essential because this is regarded the most important design parameter.In contrast, larger errors for the drag and pitching moments are acceptable due to the smaller magnitudes.In general, errors lower than 4%-5% (green) are within typical    simulation tolerance and, therefore, acceptable.For errors between 5% and 10%-15% (amber), engineering judgement should take consideration of the discrepancies.Larger errors (red) could cause wrong aerodynamic design directions.Action should be taken to improve the model, by iteratively including new training experiments and regenerating the model until acceptable error is achieved.In practice, however, the error tolerance is determined by application.For example, multiphysics simulations (e.g.aeroelastic analyses) require stricter modelling tolerances from each separate model than single physics simulations (e.g. common aerodynamic responses).Nevertheless, there is no risk of unpredicted structural failures as safety factors are enforced to be over 200%.

Predicted distributions
Figures 13 and 14 analyse the predicted pressure coefficient field for the two labelled conditions in figure 12.
The contour plots illustrate the reference C P solution from CFD (left), the prediction by our model (middle)  and the error (right).Panel (d) compares the C P distributions at the cross sections specified in the right contour.We observe that in sample e the prediction is in good agreement over the whole surface, figure 13.For sample f, the model is overall correct except for a small discrepancy on the shock-wave location, predicted slightly further downstream from 30% of the span, figure 14.Figures 15 and 16 provide a similar comparison for the skin friction coefficient C f .The shear-stress field is also found correct in the first condition.In the second case, we observe that the boundary-layer separation is slightly delayed.This is in line with the predicted location of the shock wave shown in the previous figure.Remarkably, this result indicates that our proposed framework appears to understand the relationship between the physical quantities.This demonstrates the reason for developing a single model for multiple target fields.

Full envelope prediction
The implemented ROM is convenient to efficiently interrogate the complete operating envelope.The right panels in figure 12 present the resulting aerodynamic maps.The 3D contour plots were built by integrating the surface field predictions of 30 × 30 samples uniformly distributed.The C L correlates with α ∞ as expected, with the nonlinear C L slope at high angles of attack also well captured by the model.In addition, a shallow valley along the diagonal is visible.This seems to be caused by the shock wave intensifying while the peak pressure gradually decreases.The isocontours on the C D envelope reveal a nonlinear behaviour along the α ∞ axis, which resembles the expected quadratic dependence, especially at low M ∞ .The drag increase at higher Mach numbers is related to the boundary-layer separation induced by the shock wave.Last, the C My envelope reveals highly nonlinear phenomena.The pitching moment decreases (in magnitude) with α ∞ caused by the shock wave intensifying and moving upstream.By contrast, C My is largest at low α ∞ and M ∞ ∼ 0.81 consequence of the downstream location of the shock wave.The small spike observed at high α ∞ is likely consequence of the lack of sampling in that region.

Note on computing costs
A summary of the computing costs involved in the deployment of the framework and the significant saving compared to high-fidelity simulations is reported in table 3. The steady-state CFD conditions by Immordino et al [13] were solved with a 120 core HPC, totalling up to 16 000 CPU-h for the 40 training CFD samples.
The ROM was generated with a 6 GB GPU and the training process required around 2.5 GPU-h.New aerodynamic predictions are completed in less than a second, rather than almost 3 h in CFD, i.e. a speed-up of well over 99.9%.As a result, the 30 × 30 samples to construct the 3D envelopes in the right Panels of figure 12 were completed within minutes, whereas it would be impractical using only CFD.Furthermore, to address a typical aerodynamic characterisation campaign [37], comprising ten points along the angle-of-attack axis and a Mach-number resolution of 0.02, the overall computing gain, including the generation of the training dataset in CFD, would still be up to 50%.

Conclusions
The flow analysis around a 3D aircraft remains an expensive task despite access to larger and more performing computing services than ever before.This limitation takes on an even higher criticality when the designer is tasked with delivering the performance of the system across a range of relevant flow conditions.In engineering, the loads experienced by a reentry vehicle passing through the Earth's atmosphere, the stability and control characteristics of a transport aircraft across the flight envelope or the aerodynamic map of a racing car are examples of common tasks.To overcome the computational burden associated with running a multitude of CFD analyses, whose number is at the designer's discretion to obtain data in time for pressing deadlines, the use of ROMs is a viable alternative.However, these models still present a number of challenging decisions.The most critical decision one has to make is the choice of the mathematical structure for the ROM.
We developed a geometric deep-learning AE framework to achieve a cost-effective predictive model of output quantities of interest defined on a large spatial domain with an unstructured, irregular discretisation.The framework dealt with over 300 thousand outputs and about 80 thousand grid points distributed on a 3D discretisation of a wing/body aircraft configuration.We faced specific challenges that required the development of a novel approach, which can be taken to other disciplines with immediate applicability.The first challenge is represented by the large set of points used for the spatial discretisation.We created a dedicated multi-resolution AE for extracting multi-scale features from the data field, transferring influence across the domain and for memory efficiency.The second challenge is related to the unstructured and irregular point discretisation.We made use of GCNs that enable the convolutional operation on irregular domains using the mesh connectivity, very attractive to emulate high-fidelity computational engineering analyses.The resulting predictive framework offers a novel approach when data are defined on large 3D unstructured manifolds, which is the case for any realistic problem.Nonetheless, the framework keeps the ability to use the simpler cases of structured domains when available.
It is worth noting that the steady-state framework builds on one single model that outputs multiple vector fields distributed on an unstructured domain.For our application to aircraft aerodynamic loads, the predicted pressure and shear-stress coefficient distributions were integrated in space to calculate the total force and moment coefficients.This operation not only mimics the way integrated loads are obtained in a CFD solver, but there are advantages too.First, the generation, training and validation of one single model is noticeably cheaper and easier than doing it for two separate and distinct models.Then, it avoids having two independent models that may be best fit for their specific outputs but do not recover the actual relationship among the distinct outputs, i.e. integration in space for our aerodynamic problem.Finally, a physically sound distribution of flow quantities obtained from one single model leads to a sound interpretation of aerodynamic loads.
The NASA CRM wing/body aircraft configuration was used for demonstration.We used an existing database of 70 pre-computed cases to generate and validate the predictive model.In the reference study that provided us the database, sample points were placed across the flight envelope using a latin hypercube method, which is not receptive of any feature learned during the design space exploration.The model predictions achieved a good match to reference data across the flight envelope.It is expected that further improvements in model predictions are obtained using an adaptive design of experiments where sample points are placed at strategic locations of the design space.Once the model is generated, load predictions across the whole flight envelope of angle of attack and Mach number is possible within minutes.A thorough study demonstrated that the model captured the various nonlinear effects throughout the envelope, including variations of shock wave strength and position with the angle of attack and Mach number, and the appearance of shock-induced boundary-layer separation at certain flow conditions.
Today, there is an abundance of data from calculations and measurements.Our choice of using an existing database reflects this situation.Cost-wise, model predictions are obtained at a minimal cost, 0.1%, compared to running the computational fluid dynamics solver.In design applications, the issue of database generation is still actual.To minimise the preliminary data requirements, a combined data-driven with physics-knowledge implementation could be thought by embedding physics-informed loss terms during training.However, application of physics terms on the surface of a complex geometry is not trivial and proper consideration is sought as the fluid dynamics equations are formulated for the fluid volume.We believe the generalisation to variable shapes is attractive to further expand the applicability of our model.Since our model is designed to include the coordinates of the mesh nodes as inputs, the current implementation may be adapted to wing deflections for static, 3D aeroelastic analysis and for aerodynamic shape design and optimisation.

A.3. Sensitivity to model architecture hyperparameters
The influence from the various key hyperparameters of the proposed model architecture is analysed here.In particular, we are interested in assessing the sensitivity of the model performance to the training set size, the number of model weights or the multi-mesh cycle reduction.

A.3.1. Sensitivity to MM compression rate
Table A4 reports the prediction error summary for various mesh compression rates adopted in the MM cycle.Different models were generated for mesh coarsenings to 10 000, 5000 and 2500 grid points, respectively.A compression ratio larger than 30 was not pursued.The final mesh coarsening choice adopted in the Results section is highlighted in bold.We observe that the prediction error statistics are similar among the various MM compressions.Despite the prediction error is expected to worsen with coarser meshes, the results suggest that the various GCN blocks are able to compensate for this.The model performance remains acceptable even on significantly large mesh reductions.

A.3.2. Sensitivity to MM cycle levels
The sensitivity to the number of MM reduction layers is assessed in table A5.Results are reported for MM cycles of 2 (final choice highlighted in bold), 3 and 4 levels.The three-level MM cycle was constructed by first reducing to 20 000 nodes and subsequently reducing to 5000, assigning the same node count as for the two-level case.The four-level scheme was achieved by reducing to 40 000, 20 000 and 5000.GCN blocks were embedded between mesh levels, for a total of 174 460, 352 156 and 436 252 model weights, respectively.The statistics are similar among the various MM cycles.By contrast, the computational cost involved in the three-level and the four-level schemes was found, respectively, 32% and 81% larger compared to the two-level implementation.This motivated our choice of the two-level reduction as a more efficient implementation in terms of model memory requirements.A quantification of the prediction confidence for each multi-mesh cycle is illustrated in figure A3.Inter-quartile range plots for the three different MM levels and the various performance metrics, C L error, C D error and C My error.The uncertainty intervals are based on predictions of the validation dataset.The error for each validation sample is also illustrated in triangles.We observe that the three-level option was found with the marginally lowest variability.Nevertheless, the uncertainty is similar for the various MM levels, which provides confidence on the final choice of the two-level MM cycle.

A.3.3. Sensitivity to training dataset size
We now analyse how the prediction performance is affected with a smaller training dataset.In particular, the number of training samples was halved, i.e. using 20 samples to generate the model as opposed to the original 40.Minimising the amount of data requirements is interesting to reduce the computational cost involved with running the CFD simulations.The prediction error summary for the various load coefficients is reported in table A6, with the original dataset highlighted in bold.A degradation of the performance is observed with the reduced training set.However, the error values remain still reasonably good considering the small number of preliminary samples to generate the model.

Figure 1 .
Figure 1.Surface mesh representation of the common research model.

Figure 2 .
Figure 2. Dataset samples across the operating envelope of the CRM model.Circles correspond to training samples and triangles, validation cases.Lettered labels indicate selected samples to showcase the different physics across the envelope, figures 4 and 5.

Figure 3 .
Figure 3. Integrated aerodynamic force coefficients of the NASA CRM wing/body configuration for the reference dataset, classified by training (circles) and validation (triangles) samples.

Figure 4 .
Figure 4. Pressure coefficient distribution, CP, at the four sample points of figure 2.

Figure 5 .
Figure 5. Skin friction coefficient distribution, C f , at the four sample points of figure 2.

Figure 6 .
Figure 6.CRM mesh represented as a graph with node features and edge weights.

Figure 7 .
Figure 7. Autoencoder concept, involving the embedding of GCN blocks in the reduction (encoder) and reconstruction (decoder) of the field data.

Figure 8 .
Figure8.The two-level multi-mesh scheme for the CRM test case, showcasing the resulting pressure reconstruction from one mesh coarsening-refining cycle.

Figure 9 .
Figure 9. Probability density function to select the grid points in the mesh coarsening operation.Grid points ordered by face area in ascending order.

Figure 10 .
Figure 10.Graph coarsening procedure and connectivity regeneration as part of the dimensionality reduction algorithm.

Figure 11 .
Figure 11.Steady-state GCN-MM-AE model architecture for aerodynamic predictions of the CRM test case.

Figure 12 .
Figure 12.Aerodynamic coefficient predictions with the steady-state GCN-MM-AE model for the CRM test case.Error plots on the left are classified by training (circles) and validation (triangles) samples; percent error below 4% in green, amber between 4% and 10%, and red larger than 10%.

Table 1 .
Surface mesh sizes (number of nodes) of the two multi-mesh levels for the CRM case.

Table 2 .
Statistical summary of the prediction error for the CRM test case and various datasets.

Table 3 .
Computing costs for steady-state modelling of the CRM test case.

Table A4 .
Prediction error comparison for various MM compression ratios on the CRM validation dataset.

Table A5 .
Prediction error comparison for different MM cycle levels on the CRM validation dataset.