Graph convolutional networks applied to unstructured flow field data

The key idea behind our approach is to use a graph representation to describe the connectivity of unstructured data points. In order to resolve flow fields around complex geometries in computational fluid dynamics, an irregular mesh around the immersed object is created. Next, numerical methods are applied to calculate the flow field data. Any mesh structure around the immersed object in the flow field can be represented as a graph, considering mesh nodes as vertices and using edges to connect the neighbors. Therefore, one can construct a graph representation of the unstructured mesh with different complexities or number of nodes around arbitrary objects.

Specifically, we define an undirected graph $\mathcal{G} = (\textbf{V}, \textbf{E})$, with V vertices describing the nodes of the mesh, and E edges representing the connectivity of the mesh. The flow field data is defined on mesh nodes, resulting in a feature matrix that contains the input features for V graph vertices. The edge connections are encoded as the adjacency matrix, a binary V × V matrix indicating whether any given pair of vertices are connected.

GCNNs are the generalization of CNNs for operation on graphs. GCNNs, like CNNs, are able to extract multi-scale spatial features through the use of shared weights and localized filters [42]. However, as discussed earlier, traditional CNNs are unable to work with unstructured data. GCNNs can bypass this limitation by defining the convolution operation based on the structure of the graph. By propagating information through each node's local neighborhood as defined by the adjacency matrix, GCNNs are invariant towards the order in which the nodes are specified in the feature matrix. GCNNs are often used for tasks such as node classification, link prediction [43], and graph classification.

GCNNs can be described in terms of a general framework for learning on graph-structured data, called message passing neural networks (MPNNs) [44]. MPNNs develop hidden state embeddings, h_v for each node v during the training process. Supervised training of a graph neural network aims to learn a state embedding from the features defined on the nodes and edges of the graph to have the best possible mapping to the output. The training process consists of two phases, the message passing phase where hidden states aggregate information from their surrounding nodes, and the readout phase where a feature vector for the graph is computed from the hidden states. The message passing process is parameterized by two functions, the message function M_t and the node update function U_t , while the readout function is given by R. R, M_t , and U_t are all learned differentiable functions that are updated during the training process.

Defining e_vw as the edge connecting node v to node w and $\mathcal{N}(v)$ as the neighborhood of node v, the message passing phase can be formalized as:

$$\begin{equation} m_v^{t+1} = \sum_{w \in \mathcal{N}(v)} M_t(h_v^t, h_w^t, e_{vw}) \end{equation} \tag{ 1 }$$

$$\begin{equation} h_v^{t+1} = U_t(h_v^t, m_v^{t+1}). \end{equation} \tag{ 2 }$$

Based on h_v , the readout phase computes a feature vector as:

$$\begin{equation} \hat{y} = R({h_v^T | v \in G}). \end{equation} \tag{ 3 }$$

A standard framework used is the Laplacian based GCNN, detailed in [42], where M_t and U_t are defined as the following:

$$\begin{equation} M_t(h_v^t, h_w^t) = \tilde{A}_{vw}(\deg(v)\deg(w))^{-\frac{1}{2}}h_{w}^t \end{equation} \tag{ 4 }$$

$$\begin{equation} U_t(h_v^t, m_v^{t+1}) = \sigma((W^{t})^{T}m^{(t+1)}) \end{equation} \tag{ 5 }$$

where $\tilde{A}_{vw}$ is the adjacency matrix describing the connectivity of the graph, assuming that each node is connected to itself, $\sigma$ is a nonlinear activation function, and $\deg(v)$ is the number of nodes connected to node v [44].

The GraphSAGE method, introduced by [45] is in close relation to the Laplacian based GCN layers described in the message passing formalization above. In the version of GraphSAGE implemented in this paper, the GraphSAGE algorithm is the inductive variant of this message passing network, with specific modifications to improve the accuracy and efficiency of the model. GraphSAGE acts in an inductive manner, operating on each node rather than on the entire graph, as it uses each node's local neighborhood in order to learn a function that can generate appropriate node embeddings. This method first samples a fixed number of nodes from the k nearest neighbors of each node and then applies an aggregation operator to transfer information to the node itself (algorithm 1). The aggregator can be a weighted averaging operation with trainable parameters. This inductive framework is especially helpful in the case of large graphs where low-dimensional embeddings of the nodes are more important [45]. In this regression application, the readout phase consists of a fully connected neural network that predicts a single global value for each graph based on the hidden state embeddings.

Algorithm 1. GraphSAGE embedding generation algorithm, reproduced from [45].
Input: Graph $\mathcal{G} = (\textbf{V}, \textbf{E})$; input features $\{x_{v},\forall{v}\in \mathcal{V}\}$; depth K; weight matrices $\textbf{W}^k,\forall{k}\in \{1,{\ldots},K\}$; non-linearity σ; neighborhood function $\mathcal{N}(v) = \{u\in \mathcal{V} : (u , v) \in \mathcal{E}\}$
Output: Vector representations zv for all $v\in \mathcal{V}$
1 $h_{v}^{0} \gets x_{v},\forall{v}\in \mathcal{V}$
2 for $k = 1{\ldots}K$
3 for $v\in \mathcal{V}$
4 $h_{v}^{k} \gets \sigma(\textbf{W}^{k}\cdot{\scriptstyle MEAN}(\{h_{v}^{k-1}\} \cup \{h_{u}^{k-1},\forall{u} \in \mathcal{N}(v)\}))$
5 $h_{v}^{k} \gets h_{v}^{k} / \|{h_{v}^{k}}\|_{2}, \forall{v} \in \mathcal{V}$
6 $z_{v} \gets h_{v}^{K},\forall{v}\in \mathcal{V}$

Flow field meshing usually results in a relatively large graph around the objects in comparison to other common applications such as molecular graphs. Hence, we use a node level embedding graph convolution operator that is based on the average aggregator in the GraphSAGE framework. Since the mesh has a specific edge connection pattern in which each node is only connected to a few neighbors, the sampling operator is not used. In our convolution operation, the features in the k nearest rings in the neighborhood of each node are transferred to the center node by a trainable aggregation operation (figure 1(b), line 4 in algorithm 1). Here, we use the Pytorch Geometric [46] library to load the data and implement the graph convolutional layers and pooling.

Zoom In Zoom Out Reset image size

For this problem, we implement a GCNN using GraphSAGE convolutional layers and top-K pooling steps, similar to the approach described in [47]. Different from [47], we use an inductive convolution as opposed to a transductive convolution. Inductive methods are capable of generalizing to graphs with different structures, here allowing for prediction on meshes with varying resolutions. The inputs of the network are graphs with node level velocity features, while the output is the value of the predicted drag force. Specifically, we use two GraphSAGE layers, each followed with a top K pooling layer. Top K pooling is a downsampling method to reduce the size of the layers by selecting the most important features. In top K pooling layer, a learned score is assigned to each node, and the nodes with the K highest scores are selected to be passed to the next layer [48]. The output of each pooling layer is pooled twice, once using global mean pooling and then using global max pooling, and the output from each pooling operation are concatenated together. While the input size to the network can vary in different samples as they have various number of nodes, the output size should be the same. By pooling along the dimension of vertices, the global pooling operations result in the same output size. The pooled, concatenated vectors are added together as a 'skip connection', to reinforce the information contained in the sparse convolved feature maps. The output from this step is then fed to a fully connected network with three hidden layers, that predicts the drag force (figure 1(c)). The training details of the GCNN are provided in table A1 and the parameters are defined in [46, 49] for interested readers. The effect of some of the hyperparameters and their selection are also examined in table A3, and figure 1.

In order to test the performance of our method, we aim to predict the drag force on the airfoils directly from the unstructured flow field velocity. To generate the airfoils, coordinate files are extracted from the UIUC airfoil database which contains the cartesian coordinates outlining the shape of the airfoil [50]. The incomplete or non-meshable samples are then removed from the dataset. In addition, the geometries are normalized to have a unit chord length. Next, each airfoil coordinate file is imported to the open-source mesh generator GMSH [51]. Meshes are created to reflect the variation in the density of information contained in the domain, with a finer mesh on the area close to the airfoil that resolves the complex boundary layer effects, and a coarser mesh further from the airfoil, where the flow is minimally affected by the presence of the airfoil.

To compute the velocities at mesh nodes around each object and the corresponding drag force, we perform CFD simulations using the FEniCS [52] package. FEniCS supports the DOLFIN PDE solver, which is used to solve the incompressible Navier–Stokes equations with an incremental pressure correction scheme (IPCS) method [53]. The boundary conditions for these CFD simulations are a uniform velocity input of 1.5 m s⁻¹ at the inlet (left), a far-field pressure condition at the outlet (right) and slip conditions at the top and bottom interfaces (figure 2). The viscosity and the density of the flow are 0.001 $\mathrm{Pa}$  $\mathrm{s}$ and 1 $\textrm{kg}\,\textrm{m}^{-3}$ respectively. Selected airfoil samples with their velocity magnitude field along with their drag force are provided in table 1. There is a positive correlation between the thickness of the airfoil and the corresponding drag force.

Zoom In Zoom Out Reset image size

Table 1. The velocity field and drag force for four different airfoil samples from the dataset.

The drag on an airfoil A is calculated as follows:

$$\begin{equation} F_D = \int_A (\sigma \cdot n) \cdot e_x \mathrm{dS} \end{equation} \tag{ 6 }$$

where σ is the Cauchy stress tensor, e_x is the horizontal unit vector, n is the unit vector normal to the airfoil surface.

A grid convergence study is used to choose a specific mesh size to fully resolve the boundary layer effect while minimizing the computational time required. The resulted meshes contain 900–1500 mesh points and 3000–4000 edges, generated in effectively random spatial positions surrounding the airfoil.

The incompressible Navier–Stokes equation are given by

$$\begin{equation} \rho \left(\frac{\partial \mathbf{u}}{\partial t} + u\cdot \nabla \mathbf{u} \right) = \nabla \cdot \mathbf{\sigma}(\mathbf{u}, p) + f \end{equation} \tag{ 7 }$$

$$\begin{equation*} \nabla \cdot \mathbf{u } = 0. \end{equation*}$$

To predict the velocity field at the next time-step ($u^{n+1}$) from an existing time-step (uⁿ ) while enforcing mass conservation, an IPCS is used to iteratively solve equation (6). A detailed description of the IPCS method can be found in [53].

The CFD outputs of each sample are then processed to generate the matrix of node level horizontal and vertical velocities as well as the adjacency matrix. However, storing an N × N adjacency matrix is memory-intensive. To bypass this issue, we instead store a matrix of dimension $2~\times E$, where E is the number of edges, compactly encoding the adjacency matrix. This compact representation only stores the two nodes that each edge connects in the graph. This representation is specifically helpful where the adjacency matrix is sparse which is true for the graphs in our dataset, as there are over 1200 nodes in the graph, and each node is only connected to five other nodes on average.

We perform our approach on two different sets of data. The first study is supposed to only examine the geometry of the airfoils. Therefore, the dataset consists of 1550 different airfoils. The second dataset, however, covers not only different geometries but also various angles of attack. 21 angles of attack are considered for 522 airfoils (10 962 samples in total). Angles of attack are changed in the range of $-10^{\circ}$ to 10^∘ with an increment of 1^∘. Given the relatively low velocity in the domain and small angles of attack, the flow regime is laminar, and no significant flow separation occurs.

Before implementing our method on the airfoil dataset, we perform a study of the relationship between the airfoil geometry and its drag force while other parameters are held constant. Without considering the effect of the AoA, we show the drag force to have a positive correlation with the thickness of the airfoil (figure 3(a)). However, this correlation, on its own, cannot cover a sufficient portion of the variance in the drag force. In order to determine the geometric features influencing the magnitude of the drag force, we conduct a principal component analysis (PCA) on the geometry of the airfoil and label the samples by their drag forces. By PCA, which is a linear dimensionality reduction technique, we extract components that can mostly describe the data variance in an unsupervised manner. While PCA is not generally interpretable, here we can observe the correlation of main components with the drag labels (figure 3(b)).

Zoom In Zoom Out Reset image size

For the first experiment, we use the aforementioned GCNN architecture to predict the drag force for the dataset of airfoils with zero angles of attack. $80\%$ of the samples are randomly selected for training and the remainder are used as a test set. Two complementary metrics of mean squared error (MSE) of drag prediction and the coefficient of determination (R²) are used to quantitatively evaluate the performance. Figure 4(a) shows the evolution of the loss metric as training progresses. The use of skip-connections in the architecture improves the model's accuracy.

Zoom In Zoom Out Reset image size

Using node level velocities as the input to the network, the GCNN is anticipated to detect the most important features from the flow field data to accurately estimate the drag force. To illustrate the node embeddings produced by the convolution network, we analyze the values at the input to the first fully connected layer in the trained network which is the averaged output of the convolution layers.

To do so, we perform a PCA on the features and to detect the two most important geometric components that determine the drag force. It is noteworthy to emphasize that there is no geometrical feature directly encoded in the input of the network. Smooth transition of drag values with two components and depiction of the geometry of samples indicate that the network could learn meaningful geometrical features from flow field data. Here, the first two principal components can explain more than 90% of the variance in the dataset. The first component encodes a measure of airfoil thickness and the second component is an approximate measure of how quickly the airfoil tapers (figure 4(b)).

The network is also implemented on the second dataset, containing airfoils that vary in geometry and AoA, adding complexity to the prediction task. A comparison of the ground truth values of the drag forces from CFD result and the predictions from the network qualitatively shows the high accuracy of the model for both datasets (figure 5).

Zoom In Zoom Out Reset image size

In addition to GCNNs, other machine learning and neural network algorithms can be applied on the flow field velocity data to solve a regression problem of drag prediction. Notice that the graph size and node order is not a matter of importance for the GCNN as we pass the adjacency and feature matrix directly to the model. However, non-graph based methods require a specified input size, as well as an identical node order between samples. Since the model perceives each input element as a different feature, it cannot be trained unless the order of the node elements are consistent between samples.

In order to benchmark the GCNN's performance against traditional, structured machine learning methods, we construct a node ordering that is consistent across samples. Traditional machine learning methods require the construction of a feature matrix, where the information described by a single feature—i.e., a single column of the feature matrix—must be consistent from sample to sample. In this formulation, feature i in sample x represents the same information as feature i in sample y. The node ordering is provided by the mesh generation software based on the order in which the nodes are generated, and is the same for each individual sample. Therefore, since the spatial density of the nodes in each sample is similar, this creates a matrix where Node i in sample x is relatively close in space to Node i as it appears in all other samples. To create a matrix with the same number of features for each sample, the closest 1000 nodes to the center of the airfoil are taken as the feature vector, as there are at least 1000 node measurements in each individual sample. To quantify the spatial similarity of nodes of the same index in this dataset, the distance between nodes of the same index are computed. After comparing the distance from node i in sample x with node i across all other samples, 98% of these distances are smaller than $5~\times 10^{-3} L$, where L is the length of the domain. This indicates that the position of an arbitrary node is approximately consistent across samples.

To test the performance of non-graph based methods, we select a variety of the most used methods for performing prediction. Therefore, we compare the performance of Gradient Boosted Random Forest regression, a fully connected neural network, and a two-dimensional CNN on predicting the drag force based on a matrix of node features that adhere to the previously defined structure. Some basic details of these models are provided in table A2. The models have undergone hyperparameter tuning using the package HyperOpt [54], which uses a Tree-Structured Parzen Estimator based algorithm to search for the optimal hyperparameters given a suitable range. The results of the tuning process are tabulated in table A3. A comparison shows that the GCNN approach outperforms the non-graph based methods (figure 6). It is noteworthy that the consistent feature order used in this experiment allows the GB algorithm to achieve a comparable accuracy to the GCNN, but this may not be applicable in cases where this information is not readily available. In fact, when the features are not provided to the model in a consistent order, the performance of the Gradient Boosting model sharply drops (figure 2).

Zoom In Zoom Out Reset image size