System and method for combining differentiable partial differential equation solvers and graph neural networks for fluid flow prediction

A computer-implemented method includes receiving a coarse mesh input that includes a first set of nodes, wherein the coarse mesh is input to a computational fluid dynamics solver with physical parameters to obtain a coarse mesh solution, receiving a fine mesh input that is of a second set of nodes, wherein the second set of nodes includes more nodes than the first set of nodes, concatenating the fine mesh input with the physical parameters and run the concatenation through a graph convolution layer to obtain a fine mesh hidden layer, upsampling the coarse mesh solution to obtain a coarse mesh upsample including a same number of nodes as the second set of nodes, and outputting a prediction in response to at least the coarse mesh upsample.

TECHNICAL FIELD

The present disclosure relates to systems that utilize computational fluid dynamics (CFD) and CFD solvers.

BACKGROUND

Solving large complex partial differential equations (PDEs), such as those that arise in computational fluid dynamics (CFD), is a computationally expensive process. This has motivated the use of deep learning approaches to approximate the PDE solutions, yet the simulation results predicted from these approaches typically do not generalize well to new scenarios.

SUMMARY

According to one embodiment, a system that includes a processor programmed to receive a coarse mesh input that includes a first set of nodes, wherein the coarse mesh is input to a computational fluid dynamics solver with physical parameters to obtain a coarse mesh solution, receive a fine mesh input that is of a second set of nodes, wherein the second set of nodes includes more nodes than the first set of nodes, concatenate the fine mesh input with the physical parameters and run the concatenation through a graph convolution layer to obtain a fine mesh hidden layer, upsample the coarse mesh solution to obtain a coarse mesh upsample including a same number of nodes as the second set of nodes, concatenate the fine mesh hidden layer with the coarse mesh upsample, and output a prediction in response to the concatenation of the fine mesh hidden layer and the coarse mesh upsample.

According to a second embodiment, a computer-implemented method includes receiving a coarse mesh input that includes a first set of nodes, wherein the coarse mesh is input to a computational fluid dynamics solver with physical parameters to obtain a coarse mesh solution, receiving a fine mesh input that is of a second set of nodes, wherein the second set of nodes includes more nodes than the first set of nodes, concatenating the fine mesh input with the physical parameters and run the concatenation through a graph convolution layer to obtain a fine mesh hidden layer, upsampling the coarse mesh solution to obtain a coarse mesh upsample including a same number of nodes as the second set of nodes, and outputting a prediction in response to at least the coarse mesh upsample.

According to a third embodiment, a system includes, a processor programmed to receive a coarse mesh input that includes a first set of nodes, wherein the coarse mesh is input to a computational fluid dynamics solver with physical parameters to obtain a coarse mesh solution, receive a fine mesh input that is of a second set of nodes, wherein the second set of nodes includes more nodes than the first set of nodes, upsample the coarse mesh solution to obtain a coarse mesh upsample including a same number of nodes as the second set of nodes, and output a prediction in response to the fine mesh input and the coarse mesh upsample.

DETAILED DESCRIPTION

The disclosure discusses a hybrid (graph) neural network that combines a traditional graph convolutional network with an embedded differentiable fluid dynamics simulator inside the network itself. By combining an actual CFD simulator (run on a much coarser resolution representation of the problem) with the graph network, we show that we can both generalize well to new situations and benefit from the substantial speedup of neural network CFD predictions, while also substantially outperforming the coarse CFD simulation alone.

A hybrid approach that combines the benefits of (graph) neural networks for fast predictions, with the physical realism of an industry-grade CFD simulator is disclosed below. The system may include a graph convolution network (GCN), which operates directly upon the non-uniform mesh used in typical CFD simulation tasks. The use of GCNs may be important because all realistic CFD solvers operate on these unstructured meshes rather than directly on the regular grid, which may be typically used on convolutional networks to approximate CFD simulations. Second, the system may embed a (differentiable) CFD solver, operating on a much coarser resolution, directly into the GCN itself. Although typically treated as black boxes, modern CFD simulators are themselves may be utilized as (costly) “layers” in a deep network. Using well-studied adjoint methods, modern solvers can compute gradients of the output quantities of a simulation with respect to the input mesh. This allows the system to integrate a fast CFD simulation (made fast because it is operating on a much smaller mesh) into the network itself, and allows us to jointly train the GCN and the mesh input into the simulation engine, all in an end-to-end fashion. The combined approach may perform substantially better than the coarse CFD simulation alone (for example, the network may be able to provide higher fidelity results than simply running a faster simulation to begin with), and generalizes to novel situations much better than a pure graph network-based approach. The approach may be substantially faster than running the CFD simulation on the original size mesh itself. As compared to a CFD solver used without deep learning, the disclosure of the various embodiments shown below perform simulations more accurately, and faster because the CFD solver is operating on a coarser mesh. Relative to deep learning approaches that do not incorporate CFD solvers, this method is also more capable of generalizing to new situations.

FIG. 1shows a system100for training a neural network. The system100may comprise an input interface for accessing training data192for the neural network. For example, as illustrated inFIG. 1, the input interface may be constituted by a data storage interface180which may access the training data192from a data storage190. For example, the data storage interface180may be a memory interface or a persistent storage interface, 30 e.g., a hard disk or an SSD interface, but also a personal, local or wide area network interface such as a Bluetooth, Zigbee or Wi-Fi interface or an ethernet or fiberoptic interface. The data storage190may be an internal data storage of the system100, such as a hard drive or SSD, but also an external data storage, e.g., a network-accessible data storage.

In some embodiments, the data storage190may further comprise a data representation194of an untrained version of the neural network which may be accessed by the system100from the data storage190. It will be appreciated, however, that the training data192and the data representation194of the untrained neural network may also each be accessed from a different data storage, e.g., via a different subsystem of the data storage5interface180. Each subsystem may be of a type as is described above for the data storage interface180. In other embodiments, the data representation194of the untrained neural network may be internally generated by the system100on the basis of design parameters for the neural network, and therefore may not explicitly be stored on the data storage190. The system100may further comprise a processor subsystem160which may be configured to, during operation of the system100, provide an iterative function as a substitute for a stack of layers of the neural network to be trained. Here, respective layers of the stack of layers being substituted may have mutually shared weights and may receive as input an output of a previous layer, or for a first layer of the stack of layers, an initial activation, and a part of the input of the stack of layers. The processor subsystem160may be further configured to iteratively train the neural network using the training data192. Here, an iteration of the training by the processor subsystem160may comprise a forward propagation part and a backward propagation part. The processor subsystem160may be configured to perform the forward propagation part by, amongst other operations defining the forward propagation part which may be performed, determining an equilibrium point of the iterative function at which the iterative function converges to a fixed point, wherein determining the equilibrium point comprises using a numerical root-finding algorithm to find a root solution for the iterative function minus its input, and by providing the equilibrium point as a substitute for an output of the stack of layers in the neural network. The system100may further comprise an output interface for outputting a data representation196of the trained neural network, this data may also be referred to as trained model data196. For example, as also illustrated inFIG. 1, the output interface may be constituted by the data storage interface180, with said interface being in these embodiments an input/output (‘IO’) interface, via which the trained model data196may be stored in the data storage190. For example, the data representation194defining the ‘untrained’ neural network may during or after the training be replaced, at least in part by the data representation196of the trained neural network, in that the parameters of the neural network, such as weights, hyperparameters and other types of parameters of neural networks, may be adapted to reflect the training on the training data192. This is also illustrated inFIG. 1by the reference numerals194,196referring to the same data record on the data storage190. In other embodiments, the data representation196may be stored separately from the data representation194defining the ‘untrained’ neural network. In some embodiments, the output interface may be separate from the data storage interface180, but may in general be of a type as described above for the data storage interface180.

FIG. 2is an example of a diagram of the CFD_GCN model. The network operates over two different graphs, a “fine” mesh211over which to compute the CFD simulation, and a “coarse” mesh201(initially a simple coarsened version of the fine mesh, but eventually tuned by the model) that acts as input to the CFD solver203. As input, the network takes a small number of parameters that govern the simulation. For the case of the experiments in the embodiments described below, in which the system may predict the flow fields around an airfoil, these parameters202may include the Angle of Attack (AoA) and the Mach number. These parameters202may be provided to the CFD simulation and may also appended to the initial GCN node213features. Although it may appear to be a relatively low-dimensional task, these two components can vary the output of the simulation drastically and are difficult for traditional models to learn when generalizing outside the precise range of values used to “train” the network. Additionally, the CFD-GCN model may be used in an application that utilizes temperature and angle of attack or other design parameters to determine fluid flow and turbulence through pipes, or an application that utilizes geometric or shape optimization parameters to determine noise or turbulence generated by blades of a fan or other moving parts of an object. This may include, but not be limited to, applications utilized in heat sink applications, fan blade circulation, propeller, circulation, etc.

The network may operate by first running a CFD simulation on the coarse mesh input201, while simultaneously processing the graph defined by the fine mesh211with one or more GCNs213. The system may then upsample207the results of the simulation, and concatenates these with an intermediate output215from a GCN at step209. In yet another alterative embodiment, the system may combine the not-upsampled coarse mesh solution U0and the fine mesh by just taking all of the nodes from each and forming a mesh from that, so that the number of nodes in the resulting mesh is (# of coarse mesh nodes)+(# of fine mesh nodes).

Finally, the system may apply additional GCN layers217,213to these joint features, ultimately predicting the desired output values225(in this embodiment, it may be the velocity and pressure fields at each node in the fine mesh211). For example, the fine mesh211may include 6600 nodes, while a coarse mesh201may include 550 nodes for the same airfoil or same set of parameters. The fine mesh211may include 1.5 times, 2 times, 3 times, 4 times, 5 times, 10 times, etc. the nodes greater than the coarse mesh201. Both of the meshes may be mixed triangular or quadrilateral meshes in the various embodiments.

The graph structure may use for the CFD-GCN may be directly derived from the mesh structure used by traditional CFD software to simulate the physical system. The system may consider a two dimensional, triangular mesh M (X;E;B). The first element, Xi∈N×2, may be a matrix containing the (x, y) coordinates of the N nodes that compose the mesh. The second,
E{(i1;j1;k1), . . . ,(im;jm;km)},

may be a set of M triangular elements defined by the indices (i; j; k) of their component nodes. The third,
B{(i1;b1), . . . ,(iL;bL)},

may be a set of L boundary points, defined as a pair consisting of the index of the node and a tag b that identifies which boundary the point belongs to (e.g. airfoil, farfield, etc.). Such a mesh M clearly defines a graph GM=(X;EG) whose nodes are the same X, and whose edges EGcan be directly inferred from the mesh elements E. Conversely, a graph can also be converted into a mesh if the structure of its edges is appropriate and a set of boundary points B is provided.

In addition to the fine mesh211used to compute the CFD simulation, the system may also consider a coarse mesh201, denoted Mc. The coarse mesh201has the same structure as the fine mesh M, with the number of nodes down sampled by over 10× (or any other value), which thus allows for much faster simulation. Although this mesh also technically defines a graph, the system may not directly compute any GCN over this graph, but instead only use it as input to the simulation engine.

In addition to the graphs themselves, the model also receives as input two physical parameters202that define the behavior of the flow around an airfoil: the angle of attack (AoA) and the Mach number. These two parameters are both fed into the simulation and appended as initial node features for every node in the GCN. These two parameters, or any other parameters for other applications, may be the quantities that vary from simulation-to-simulation, and thus the main task of the GCN is to learn how to predict the resulting flow field from these two parameters that define the simulation.

A central component of the CFD-CGN model may be the integrated differentiable fluid simulator203. As input, the fluid simulator203may take coarse mesh201MC, plus the angle of attack and Mach number, and outputs predictions205of the velocity and pressure at each node in the coarse graph. The system may employ the SU2 fluid simulator203(in one example), which is an open source, industry-grade CFD simulation widely used by many researchers in aerospace and beyond. Briefly, SU2 uses a finite volume method (FVM) to solve the Navier-Stokes equations over its input mesh. Crucially for our purposes, the SU2 solver also support an adjoint method which lets the system differentiate the outputs of the simulation with respect to its inputs and parameters (in this case, the coarse mesh201MCitself, plus the angle of attack and Mach number).

Intuitively, the CFD solver203(e.g., SU2 solver) should be thought of as an additional layer in our network, which takes the angle of attack and Mach number as input, and produces the output velocity and pressure fields at the solver output205. The equivalent of the “parameters” of a traditional layer is the coarse mesh itself: different configurations for the coarse mesh will be differently suited to integration within the remainder of the CFD-GCN. Thus, one of the main learning task for the SU2 portion of our model is to adjust the coarse mesh in a manner than eventually maximizes accuracy of the resulting full CFD-GCN model. The adjoint method in SU2 uses reverse-mode differentiation, so gradients can be efficiently computed with respect to a scalar-valued loss such as the overall predictive error of the CFD-GCN.

Finally, although not strictly a research contribution, that as part of this project we have developed an interface layer between the SU2 solver203and the PyTorch library. This interface allows full SU2 simulations to be treated just as any other layer within a PyTorch module, and we hope it will find additional applications at the intersection of deep learning and (industrial-grade) CFD simulation.

The output of the coarse simulation described above is a mesh with the predicted values for each field at every node. For this to be used towards generating the final prediction, the system may need to up sample it to the size of the fine mesh. The system may accomplish this by performing successive applications of squared distance-weighted, k-nearest neighbors interpolation. The node positions in the coarse mesh may refer to the following symbols utilized interchangeably,

FirstAlternativeSymbolSymbolMeaningXCXDNode positions in the coarse meshUCDNode features in the coarse meshXFXUNode positions in the fine meshUFUNode features in the fine mesh

Let us call the first element, Xi∈NU×3the up sampled version of some coarser graph D∈ND×3. For every row U(i), with corresponding node position in the fine mesh, the system may find the set {n1, . . . , nk} containing the indices of the k closes nodes to the corresponding node positions in the coarser graph XD. The system may define each node feature in the fine mesh U(i). Thus, given a coarse mesh with positions XCand features UC, and fine mesh positions XF, this operation produces fine mesh features UF. This is done using a k-nearest neighbors' operation. Denote the node with index i in the fine mesh as XF(i), and let {n1, . . . , nk} be the k nodes (default may be 3 nodes, k=3) in the coarse mesh XCclosest to XF(i). Then the features for node i in the fine mesh are given by

As depicted inFIG. 2, the output of the coarse simulation is processed by a sequence of convolutional layers. In order to operate directly on the mesh output of the CFD simulation, the system may utilize the graph convolution network (GCN) architecture. The architecture may define a convolutional layer for graphs.

A general graph consisting of Nznodes, each with F features, is defined by its feature matrix Zi∈Nz×Fand its adjacency matrix A∈Nz×Nz. We can then further define {tilde over (B)}={tilde over (D)}−½(A+I){tilde over (D)}−1/2, where I is the identity matrix and {tilde over (D)} the diagonal degree matrix, with its diagonal given by {tilde over (D)}ii=1+Σj=0NzAij. Then, a GCN layer217with F input channels and F′ output channels, parametrized by the weight matrix W∈F×F′and the bias term b∈Nz×F′, will have as output
Zi+1={tilde over (B)}ZiWi+bi≡GCNi(Zi).

With all of the components acting together, the system can run an SU2 simulation with the coarse mesh201and the physical parameters202. The output205of this coarse simulation is upsampled L times.
U0=SU2(XC,AoA,Mach)
Ui+1=Upsample(Ui),i=0, . . . ,L

Concurrently, the fine mesh211has the physical parameters and the signed distance function (SDF) appended to each of its nodes' features. Here, the SDF at a node may be defined as its Euclidean distance to the closest node at the airfoil boundary. The resulting graph is then passed through a series of graph convolutions213, as shown by graph convolutions213,219,221,223. At some specified convolutional layer k, the final upsampled value ULis appended to the output Zk221of the k-th convolution. Another set of convolutions223is performed in order to generate the final prediction Ŷ225. Each of these steps are outlined in the formulas below:
Z0=[X,SDF(X),AoA,Mach]
Zi+1=ReLU(GCNi(Zi)),i=0, . . . ,k−2
Zk=[ReLU(GCNk(Zk−1)),UL]
Zk+i+1=ReLU(GCNk+i(Zi)),i=0, . . . ,K−k
Ŷ=GCNK(ZK)
In such an example, [·,·] may be the matrix concatenation operation over the column dimension.

Given that the entire CFD-GCN as formulated above can be treated as a single differentiable deep network (including the SU2 “layer” or CFD Solver “layer”), the training process itself may be straight forward. The model is trained to predict the output fields Ŷ∈N×3, consisting of the x and y components of the velocity and the pressure at each node in the fine mesh, by minimizing the mean squared error (MSE) loss(or any other type of loss function that compares ground truth to the network output) between the prediction Ŷ and ground truth

ℓ⁡(Y,Y^)=13⁢N⁢Y-Y^22
where the ground truth Y in this case is obtained by running the full SU2 solver to convergence on the original fine mesh. The training procedure optimizes the weight matrices Wiand biof the GCNs, and the positions of the nodes in the coarse mesh XCby backpropagating through the CFD simulation. The loss may be minimized using an Adam optimizer (or any other optimizer) with a learning rate α=5·10−5.

The system may also run a backward pass sequence to improve the system over time. To summarize, the system may obtain loss L by comparing prediction Ŷ with ground truth Y (using mean squared error, or any other loss function). The backward pass function or backward pass sequence may utilize various steps. The system may utilize the following sequences (in no specific order):

∂L∂θi
of loss with respect to graph convolution layer parameters for i=1, . . . , K
b. Gradient

∂L∂U0
of loss with respect to solver output
c. Using

∂L∂U0
and adjoint methods to compute gradients

Thus, the system may use gradients and some deep learning optimization algorithm (stochastic gradient descent, ADAM, etc.) to update parameters θiand coarse mesh positions XC.

The network may use the ability to differentiate through the solver function to learn the best positions of nodes in the coarse mesh XC. The system could place convolutional layers before the CFD solver, and could also learn parameters in those layers through the same method. The algorithm may also include a single up sampling operation (e.g., U1=Upsample(XC, U0, XF)) to expand the coarse mesh to match the size of the fine mesh; however, an alternative may also break the sample into multiple up sampling operations if desired. This may be beneficial if the coarse mesh is far smaller than the fine mesh.

In an alternative embodiment, graph convolutions may be replaced with any type of neural network layer that operates on graphs (e.g., graph self-attention layers, graph transformers), or may be replaced with convolutional layers if the mesh is in a grid. The system may change the loss function, optimizer, or activation function. The network architecture may be compatible with standard deep learning training techniques or modifications (dropout, batch normalization, layer normalization, etc). In yet another alterative embodiment, the system may combine the not-upsampled coarse mesh solution U0and the fine mesh by taking all of the nodes from each and forming a mesh from that, so that the number of nodes in the resulting mesh is (# of coarse mesh nodes)+(# of fine mesh nodes).

FIG. 3illustrates an example of a mesh. An issue may arise when optimizing the input coarse mesh. Gradually, as the node positions are moved by the gradient descent updates, it is possible that in a given triangular element, one of the nodes (e.g., node301) may cross over an edge as shown by cross-over node303. This may generate non-physical volumes, which may harm the stability of the simulations, frequently impeding convergence. In other words, at each gradient update step, the optimizer may update the mesh nodes301by performing the update
Xc←Xc+ΔXc,

with some small update matrix ΔXcof the same shape as Xc. If left unmodified, this ΔXccan cause the above mentioned issued. Thus, a node might be pushed over the edge of its triangular element, generating overlap of elements. Such non-physical situations may harm convergence of the simulation. When this occurs, the cross product between ordered edges changes. Before the updated, e1×e0>0, while afterwards ê1×e2<0.

In order to avoid non-physical volumes, the system may generate a projected update P(Δ Xc) such that only non-degenerating updates are performed. The system may start with P(Δ Xc)=ΔXc. Then, the system may check which elements in the mesh may have a node pushed over an edge by Δ Xc. This may be done by computing the cross product of two edges in each triangular element in a consistent order. If the sign of the cross product flips with the update Xv+ΔXc.

That means a node crossed over an edge since it causes the ordering of the nodes to change. This may be depicted inFIG. 3, where the cross product of the edges e1and e2is positive before the update, but negative afterwards.

For every element E (i, j, k) which has flipped, the system may set rows i, j, and k, of P(Δ Xc) to 0, thus performing no updates to those points in Xc. Since removing the updates to some nodes might cause new elements to flip, this procedure may be repeated until no points are flipped. Once the system reaches this sate, the system may perform the projected gradient update
Xc←Xc+P(ΔXc),

FIG. 4shows an example graph of the results of optimizing the nodes. The optimizing of the nodes of a mesh to improve a prediction loss both with and without the correct of the gradient update. Whereas the mesh optimized without the correction quickly degenerates and the loss diverges, the one with the projected gradient update learns smoothly, even for a higher learning rate α.