A Finite Element-Inspired Hypergraph Neural Network:
Application to Fluid Dynamics Simulations

Consider a continuum dynamical system defined in a bounded spatial domain $\Omega$, with governing equation in a general form

in which $q$ is the system state and $F$ denotes any input to the system. We discretize in time with the forward Euler scheme

where $t_{n}$ denotes the $n$-th time step. Further discretization in space with a mesh leads to

in which $\boldsymbol{Q}(t_{n})$ denotes the system state vector at time step $t_{n}$, and $\boldsymbol{F}(t_{n})$ denotes the system input vector at time step $t_{n}$. The function $\tilde{\boldsymbol{G}}$ that governs the evolution of the system in time can be approximated by a surrogate model $\widehat{\boldsymbol{G}}$. In this work, we construct such a surrogate model with a graph neural network, which leverages the geometric relationship between the data points. As such a relationship is already prescribed by the mesh, it is preferable to convert the mesh into a graph for subsequent use with the graph neural network.

As discussed in the introduction, one could transform the mesh into a graph consisting of nodes connected by edges (cf. Fig. 1), and a generalized graph message-passing network operating on such a graph would mimic the finite volume method. Instead of following this approach, we take inspiration from the data connectivity in the finite element method, and transform the mesh into a hypergraph via a different approach in this work, discussed in detail in the subsequent Sec. 2.2.

Consider a bounded spatial domain that is meshed (cf. Fig. 2a). In the finite element method, one treats each cell as an element connecting all the vertices of the cell. Mimicking this behavior, we convert each cell within the mesh into an undirected hyperedge (an ”element”) that connects all the vertices of the cell, and converts each of these vertices to a node, leading to an undirected hypergraph as shown in Fig. 2b. We further explicitly define the connection between each element and each of the nodes it connects as an undirected element-node edge, and this gives us a hypergraph $\mathcal{G}=(\mathcal{V},\mathcal{E}_{\square},\mathcal{E}_{v})$, in which $\mathcal{V}$ is the set of all nodes, $\mathcal{E}_{\square}$ is the set of all elements, and $\mathcal{E}_{v}$ is the set of all element-node edges, as illustrated in Fig. 2c. We coin the name node-element hypergraph for such a hypergraph $\mathcal{G}$.

With the node-element hypergraph converted from the mesh, one could proceed to transform the mesh information and system states and attach them as features. These transformations are usually specific to the problem in concern, and we thus delay the discussion of these details to Sec. 3. In the next subsection, we assume that the system states and mesh information is transformed into some node features, element features, and element-node edge features, and define a message-passing network in the general form.

We now proceed to define a graph message-passing network that works with the node-element hypergraph. Taking further inspiration from the finite element method, we design hypergraph message-passing layers that mimic the local stiffness matrix assembly process. In the remaining part of this subsection, we would start from the calculation of the element stiffness matrix for a simple Poisson problem (Eq. 8), and proceed step-by-step towards the formulation of the hypergraph message-passing layers (Eq. 19 and Eq. 21). For the sake of consistency and simplicity in symbols, we write the equations in this subsection and the remaining parts of the paper assuming a 2D quadrilateral mesh and the corresponding converted node-element hypergraph (cf. Fig. 2).

The behavior of the element and node update stages described in Eq. 19 and 21 are largely determined by the choice of functions $\phi^{e}$ and $\phi^{v}$. For the purpose of the neural network, we wish to retain the flexibility to learn any dynamics rather than just the local stiffness matrix calculation process, and it is therefore preferable that these functions are universal approximators that can be trained to approximate any function. Under this consideration, we follow the practice shared by most of the works discussed in the introduction, and choose each of these functions $\phi$ as a separate multi-layer perceptron, which can be written as

in which the symbol $\circ$ denotes the function composition, and $\boldsymbol{z}$ is the input vector to the multi-layer perceptron. Each layer $f_{i}$ within the multi-layer perceptron is defined as

where $\boldsymbol{W}_{i}\boldsymbol{z}+\boldsymbol{b}_{i}$ is a linear transformation of the input vector $\boldsymbol{z}$, and $\sigma_{i}$ is a non-linear function called the activation function. The values of all entries within the matrix $\boldsymbol{W}_{i}$ and the vector $\boldsymbol{b}_{i}$ are trainable.

With the message-passing layer defined on node-element hypergraph , we are able to construct the surrogate graph neural network model $\widehat{\boldsymbol{G}}$ in Eq. 4. Following the encode-process-decode architecture [23], the network includes encoders, a series of message-passing layers, and decoders, stacked together in a feedforward fashion. Separate encoders and decoders are constructed for the node features, element-node edge features, and element features respectively using multi-layer perceptrons. The model is summarized in Algorithm 1.

We consider a fluid system defined on a bounded domain on $\mathbb{R}^{2}$ with a given mesh, with flow velocity $\boldsymbol{u}=(u_{x},u_{y})$ and pressure $p$ available at each vertex within the mesh, along with the coordinate of the vertex itself $(x,y)$. We assume that the Reynolds number $\operatorname{Re}$ is known. The mesh is assumed to be body-fitting, i.e., each system boundary is defined as a set of mesh cell boundaries, with boundary conditions ${\Gamma}$ known. The input features to the neural network are transformed from these information.

Similar to the characteristics of traditional computational fluid dynamics solvers, we expect the network to satisfy a series of invariance and equivariance restrictions. In particular, the predicted pressure values on the vertices are expected to be invariant when the input system translates and/or rotates on the 2D plane. The predicted flow velocity vectors, on the other hand, are not expected to change when the system translates, but are expected to rotate together with the system. These properties are not necessarily satisfied by a graph neural network by default, and therefore the available information $(x,y)$, $(u_{x},u_{y})$ and $p$ on the vertices need to be transformed to suitable forms before being attached to features.

With the system information transformed to features of suitable form, we proceed to attach them to the graph. Since the system boundaries are body-fitting, the one-hot boundary condition vector is attached to the node feature $v_{i}$ at each node $i$:

It is intuitive to attach the element area to the element feature $e_{\square}$. As the cell area at different regions within the mesh can vary by orders of magnitude, and that the cell areas are strictly greater than zero, we take the natural logarithm of the actual area values,

in which the $-\ln S_{\square}$ feature does not provide additional information, and is used to prevent the existence of feature vector of length 1.

The remaining system information are attached to the element-node edge feature. The pressure values $p_{i}$, originally defined on each node $i$, are gathered to all the element-node edges connected to it,

in which the subscript $(\cdot)_{i\square,i}$ denotes all the element-node edges that connect node $i$ with an element. Similar to the element area, we also take the natural logarithm values for the length of local coordinate vectors $L_{\square,p}$. Since the mesh enforces that $0<\theta_{\square,i}<\pi$, we can take the cosine of the angles $\theta_{\square,i}$ for the ease of computation and normalization without the loss of information. The element-node edge feature vector, after these treatments, can be written as

When a variety of dynamics (i.e., a range of Reynolds numbers) are involved, to make it possible for the network to differentiate between similar dynamics, we explicitly supply Reynolds number $\operatorname{Re}$ as a feature, which gives the element-node edge feature vector

A discussion on the necessity of this treatment is provided in A.

In the temporal roll-out predictions, the network iteratively takes the features as inputs, and updates some of these features with its outputs. Assuming that the boundary conditions and geometry information do not need to be predicted, the features that need to be updated at every time step during the roll-out prediction are the first two entries $d_{\square,i}$ and $p_{\square,i}$ of the element-node edge feature vectors. With the system discretized temporally with forward Euler time-stepping scheme (Eq. 4), we have the temporal update from time step $t_{n}$ to time step $t_{n+1}$:

with $\widehat{\boldsymbol{G}}$ being the $\phi$-GNN defined in Sec. 2.5. The network outputs are vectors of length 2 on element-node edges, and the decoders for the node and element features $h^{v}$ and $h^{e}$ are not used.

It should be noted that the temporal update are on the element-node edges which do not represent a physical location. To retrieve the velocity $(u,v)$ and pressure $p$ fields defined on vertices of the mesh, we aggregate the element-node edge features back to the nodes as a post-processing step. For the pressure $p$, the value on each node can be simply retrieved by calculating the mean value across all the element-node edges connecting to the node

in which $\operatorname{MEAN}_{\square}$ is the mean aggregation function. The velocity at each node is retrieved by reverting the process of Eq. 27. Similar to that in Lino et al. [26], this can be achieved by solving the (usually over-constrained) problem via Moore-Penrose pseudo-inverse:

in which the symbol $\odot$ denotes the Hadamard product, $n_{i\square,i}$ denotes the number of element-node edges connecting to node $i$, $\boldsymbol{d}_{i\square,i}$ is a $2\times n_{i\square,i}$ matrix containing the projection weights $d_{i\square,i}$ for all the element-node edges connecting to node $i$, and $\boldsymbol{X}_{i\square,i}$ is also a $2\times n_{i\square,i}$ matrix containing the unit local coordinate vectors $\boldsymbol{x}_{\square,i}$ for all the element-node edges connecting to node $i$.

We choose two common benchmark fluid flow systems for the experiments: The flow around a circular cylinder and the flow around an NACA0012 airfoil. Both flow data sets are 2-D incompressible fluid flow simulated at laminar flow conditions. The flow around the cylinder data set is obtained through simulation at Reynolds number $Re=200$, and used to test the capability of the neural network in learning a certain flow dynamic. The flow around the airfoil is simulated at multiple Reynolds numbers (with the flow at each Reynolds number forming a separate ”trajectory”) within the range $Re\in[1000,4000]$, and the resulting data set is used to test the capability of the network to interpolate within a range of dynamics and extrapolate out of the range.

Both flow data sets are generated via a Petrov-Galerkin finite element solver with a semi-discrete time-stepping scheme [38] written in Matlab, with the computational meshes created using Gmsh [39]. The domain for the two sets of simulations is plotted in Fig. 5. The inlet is a fixed uniform flow $u_{x}=1,u_{y}=0$, the outlet $\Gamma_{out}$ is set to be traction free, while top and bottom of the boundary conditions $\Gamma_{top}$ and $\Gamma_{bottom}$ are set to be slip-wall. No slip boundary condition is applied at the surface of the cylinder and airfoil. The computational meshes are shown in Fig. 6a and Fig. 7a. For the flow around cylinder, a total of 6499 continuous time steps are sampled with non-dimensionalized time step $\Delta t^{*}=0.04$ at Reynolds number $Re=200$. For the flow around an airfoil, simulations are performed at 61 different Reynolds numbers within the range $Re\in[1000,4000]$. For each of the considered Reynolds number, 4500 continuous time steps are sampled with non-dimensionalized time step $\Delta t^{*}=0.0167$.

The simulated flow data are interpolated onto a coarser mesh (depicted in Fig. 6b and 7b) via a Clough-Tocher interpolator available in SciPy package [40] before subsequent conversion to features. It should be noted that the mesh for the flow around fixed cylinder is the same as the one used for the fluid-structure interaction between fluid flow and an elastically mounted cylinder in reference [41]. The statistics about the number of nodes, elements and element-node edges for the two data sets are reported in Table 1. The interpolated data sets are then split into train, cross-validation and test data sets, as listed in Table 2. It should be noted that the training and cross-validation share the same data set for the cylinder case. This is possible because the network is trained with single-step supervision, but evaluated by generating roll-out predictions for thousands of time steps starting from the first time step within the training/cross-validation data set.

We implement the models with PyTorch [42]. All multi-layer perceptrons in the network (encoders, decoders, processors) have two hidden layers with hidden and output layer width 128, except for the outputs of the element-node decoder which have width 2. A total of $n_{d}=15$ message-passing layers are used. It should be noted that the hyperparameters of network width and depth are not tuned but fixed following the choice of Pfaff et al. [23] for their MeshGraphNet. Residual links [33] are added in all message-passing layers.

Unlike most of the existing works on graph neural network-based continuum mechanics simulations discussed in the introduction, we adopt a sinusoidal activation function [43] in this work along with the associated initialization scheme. To avoid conflict with this design choice, we choose to use the mean aggregation function for both the element and node update stages, and avoid the use of any normalization layers. In addition, we apply min-max normalization with minimum $-1$ and maximum $1$ to each entry for each of the input and output features in the training data sets, with the exception of $\cos\theta_{\square,i}$ in the element-node edge features since its value is already limited to be within the range of $(-1,1)$. The test data sets are normalized with the statistics calculated from the training data sets.

Models reported in Sec. 5.1 and 5.2 are trained with Adam [44] optimizer with PyTorch default setup for a total of 200 epochs with batch size 4, with an adaptive smooth $L_{1}$ training loss that will be discussed in Sec. 4.3. The network is trained using single-step supervision, i.e., for the input features at time step $t_{n}$, the training loss is calculated between the predicted increment of element-node edge features from time step $t_{n}$ to time step $t_{n+1}$ (cf. Eq. 33) and the ground truth values. Artificial noise is added to the ground truth input and output values during training, with the scheme discussed in detail later in Sec. 4.4. A maximum learning rate of $10^{-4}$ and a minimum learning rate $10^{-6}$ are selected, with a warm-up stage of 10 epochs at the start. Details of the learning rate scheme will be presented in Sec. 4.5. All the reported training and testing runs are completed on a single Nvidia RTX 3090 GPU with CPU being AMD Ryzen 9 5900 @ 3 GHz $\times$ 12 cores.

Most of the existing works discussed in Sec. 1 adopt mean squared error (MSE) as the training loss function. For graph or hypergraph message-passing networks with sinusoidal activation functions, however, we observe that the training process can be unstable when MSE loss is used. We therefore seek an alternative training loss in this work. Another typical loss function, the $L_{1}$ loss, might not be preferable in this case since it is not smooth at zero, and we thus adopt a smooth $L_{1}$ loss [46] function, which states that for ground truth $\psi$ and its neural network prediction $\hat{\psi}$, the loss

in which $\beta$ is a non-negative parameter that controls the transition point between the $L_{2}$ loss region and the $L_{1}$ loss region. The subscript $(\cdot)_{i}$ here denotes entry-by-entry calculation. A fixed $\beta$ value is not preferable, since the loss function is not very different from the $L_{1}$ loss when $\beta$ is too small and will converge to a scaled $L_{2}$ loss during training when $\beta$ is too large. For the present case, the instability in training usually occurs when the training MSE error is relatively low, so an adaptive scheme for $\beta$ is preferred. Assuming that the distribution of error during the training approximately follows a symmetric distribution centered at zero, the target is to make sure that the two tails of the distribution fall into the $L_{1}$ loss region so that they do not lead to instability. More complex on-the-fly $\beta$ adaptation algorithms like references [47, 48, 49] exist, but we choose to adapt the value $\beta$ in this work using a simpler yet robust approach by setting $\beta^{2}$ as the variance of the model prediction error

based on the idea that at least part of the two tails of any zero-centered symmetric distribution would reside more than one standard deviation away from zero. This variance of prediction error can be approximated by computing the mean-squared error between the predicted and ground truth value of the whole training set, which can be further approximated on-the-fly by an exponential moving average

in which $N_{b}$ is the total number of training iterations within each epoch. We further stabilize this process by preventing $\beta$ from increasing, i.e.,

for each training step. The initial value of $\beta$ can be determined by calculating the mean-squared error of the initialized network over a small randomly-sampled subset of the whole training set. Such a subset contains $400$ samples for the reported cases in Sec. 5.

The prediction of the network is expected to have an error. Without special treatment, this error would accumulate over the prediction time steps and would eventually cause the predictions to blow up. It is therefore preferable that the network automatically corrects its prediction error. In this work, this is achieved by training the network with artificial noise added to the training data. With a the trained network $\widehat{G}$, we consider its prediction output $\widehat{\delta e}_{\square,i}^{n}$ calculated from the input features at time step $t_{n}$,

It should be noted that this equation is different from Eq. 33. The additional scaling factor (vector) $\boldsymbol{\kappa}$ exists due to the necessary normalization of neural network output targets during the training process. Assuming that the network is large enough and trained thorough enough, such that it is able to capture the evolution of the system dynamics, it is reasonable to assume that the error of this prediction would approximately follow a Gaussian distribution,

We hope that the network will be able to correct this error during its prediction in the subsequent time step, i.e., we expect that

in which $\epsilon_{\square,i}^{n+1}\sim\mathcal{N}(0,\operatorname{Var}(\epsilon_{\square,i}^{n+1}))$ is the network prediction error at time step $t_{n+1}$. Therefore, during training, one can enforce the network to learn such a behavior by adding a Gaussian noise $-\widehat{\epsilon}_{\square,i}^{n}$ with a specified variance to the training output target, and $\boldsymbol{\kappa}\odot\widehat{\epsilon}_{\square,i}^{n}$ to the corresponding entries of training input.

For systems that only changes slightly over each time step, i.e., both entries of $\kappa$ is much lower than 1, the predicted output at neighboring time steps would be similar to each other. It is therefore also reasonable to expect that the prediction error of two neighboring time steps are significantly correlated,

This is empirically verified during the preliminary tests, we skip the details for brevity. Based on such an observation, we further propose to train the network to over-correct the error at each time step, such that a proportion of error at the subsequent time step is also canceled. This means that for a Gaussian noise $-\widehat{\epsilon}_{\square,i}^{n}$ with a specified variance, we add $-\omega\widehat{\epsilon}_{\square,i}^{n}$ to the training output target, and add $\boldsymbol{\kappa}\odot\widehat{\epsilon}_{\square,i}^{n}$ to the corresponding entries of input, with $\omega>1$ being a hyper-parameter that controls the level of over-correction. We take inspiration from the typical relaxation factor range used in successive over-relaxation methods, and directly specify $\omega=1.2$ in this work rather than tuning it as a hyper-parameter.

For the flow around cylinder data set, the model reported in Sec. 5.1 is trained with the standard deviation of random noise designated to be $6.5\times 10^{-2}$. This value is selected after a screen through the range $[10^{-2},10^{-1}]$, and a discussion of the behavior of the network when trained with different amount of training noise is attached in B. For the flow around airfoil data sets, to reduce the amount of computational power needed, we start from the noise standard deviation level of $10^{-1}$ and divide this value by half in each trial, eventually using the noise standard deviation level $2.5\times 10^{-2}$ for the results reported in Sec. 5.2.

We choose to merge the popular cosine and exponential learning rate schemes to take the advantage of both of them, such that a scaled cosine learning rate applies near the start of the training and a scaled exponential learning rate applies near the end. In particular, a decaying learning rate curve starting from maximum learning rate $\zeta_{max}$ and ends at minimum learning rate $\zeta_{min}$ is specified by

The value $\alpha$ is equal to the proportion of training completed. An increasing learning curve from $\zeta_{min}$ to $\zeta_{max}$, used in the warm up stage at the start of the training process, is acquired by calculating the decaying learning curve and then reverse the sequence. For the models reported in this work, we change the learning rate at the start of every iteration during the warm-up stage, and change the learning rate at the start of every epoch for the rest of the training process.

The trained models are evaluated on the test data sets. As the main purpose of the neural network surrogate model is to generate stable and accurate roll-out simulations, we evaluate the models by feeding the state of the system at a certain time step, i.e. the first time step within each of the test data sets, to the model and compare the predicted system states over the next several thousands of time steps with the ground truth values from the interpolated CFD data. We evaluate the model based on the following criteria:

Starting from the first time step in each of the test data sets, we generate roll-out predictions of the states of the fluid system over the next 4000 time steps. In this subsection, we present these results for the flow around cylinder data set, using the criteria discussed in Sec. 4.6. We also verify the designed invariance and equivariance properties of the network.

With the network shown to be able to generate stable and accurate predictions when used to learn a single fluid flow dynamics in Sec. 5.1, we proceed to apply the network to learning a range of dynamics using the flow around airfoil data sets, and verify whether it can serves as an effective surrogate model for the time series predictions within this range of Reynolds numbers. In particular, we train the network with the flow around airfoil data simulated between Reynolds numbers 2000 and 3000, and inspect its behavior when applied to other flow around airfoil cases within the Reynolds number range $\operatorname{Re}\in[1000,4000]$. In the remaining part of this subsection, we will first test over the stability of the predictions, and then inspect various aspects of the predictions within the stability range.