Mesh deep Q network: A deep reinforcement learning framework for improving meshes in computational fluid dynamics Open Access

Mesh generation in computational fluid dynamics (CFD) is a critical step for stable and accurate simulations. However, one of seven key findings from a year-long NASA study published in 2014 is that mesh generation and adaptation remains a significant bottleneck in the CFD workflow.¹ One of the reasons for this is meshing tends to be done with intuition and experience, rather than a strictly mathematical approach, to get convergence and stable simulations.² In order to achieve the goals of the NASA study, autonomous mesh generation is a necessary next step. A more recent review on the progress of modeling and mesh generation since the study has indicated a growing interest in adaptive mesh techniques to aid in mesh generation.³ 

Substantial literature exists on using both classical as well as machine learning techniques for adaptive meshing.^4–6 Much of the focus on classical approaches is on pointwise error estimation and subsequent refinement.⁵ While these methods can be used to improve meshes, they often require additional backward solves to compute adjoint states with a full simulation at each step to recalculate the velocity and pressure fields. This adds computational time, as well as additional necessary features into a CFD solver, often rendering the method unfeasible with CFD codes and pipelines. More recently, many machine learning models and methods have been developed for adaptive meshing. However, these methods tend to rely on substantial up-front computational time in training data generation, where many simulations are run,⁷ and adjoint-based adaptive meshing is used.⁶ 

Iterative mesh refinement can instead be viewed as a sequential decision making process, where each region or vertex to be modified constitutes another decision. This is well suited for deep reinforcement learning (DRL). Recent surveys have shown DRL applied to CFD has been primarily focused on control and optimization.^8,9 To the best of our knowledge, no work has been done in applying DRL for CFD meshing, specifically.

While using DRL for CFD meshing is novel, applying ML techniques for faster solutions is not. Works in operator learning,¹⁰ surrogate modeling,^11–13 and system identification¹⁴ offer alternative approaches to improve CFD simulations. Operator learning and surrogate modeling are different methods to learn the temporal evolution of the system. System identification is used to find new equations that describe the system at hand. While these approaches show promise, meshing, specifically, is the focus of this work because it allows practitioners to use existing CFD solvers on the Navier–Stokes equations.

In this work, we develop a general purpose DRL framework to coarsen CFD meshes while maintaining accuracy in the calculation of coarse properties such as lift and drag. Mesh coarsening decreases simulation time and is done by selecting and removing individual vertices in the mesh. Proof of concept is demonstrated on 2D airfoils. Our framework requires a single solve prior to training. Interplay between different geometry, meshing, and simulation codes is often of practical concern. Eliminating any dependence of our method on any particular code is of key importance for wide adoption. This framework only requires modifying the mesh directly, which allows us to use standard CFD solvers to run subsequent simulations. Finally, our model makes no assumptions about the type of mesh, flow regime, solver, or dimensionality of the system. This allows our method to be easily incorporated into CFD pipelines.

The Mesh Deep Q Network (MeshDQN) framework is given below in Fig. 1. After running an initial simulation, the pressures, velocities, coordinates, and edges are saved. This data is then passed into the property calculation, where a coarse property, such as drag or lift, is calculated. The reward is then calculated from the property.

The goal of CFD is to solve the Navier–Stokes equations over a discrete mesh. This work focuses on the incompressible NS equations with viscous flow:

where ρ is the density, u is the velocity, σ is the Cauchy stress tensor, p is the pressure, and f is the forcing term.

Solving the NS equations numerically gives us pressure and velocity at each point. The drag and lift are calculated from our solution velocity and pressure as

Where F_D and F_L are drag force and lift force, respectively and σ is again, the Cauchy stress tensor.

The meshes were generated using Delaunay Triangulation¹⁵ as part of GMSH,¹⁶ with Meshio¹⁷ used as an intermediary to convert mesh files into a useable format for the CFD solver. Simulations for this work were done over 5 s with the incremental pressure correction scheme (IPCS)¹⁸ implemented in Dolfin as part of FEniCS,^19,20 based on the FEniCS Tutorial Volume 1.²¹ 

In reinforcement learning (RL), an agent interacts with an environment through actions. The agent takes in a state s, selects an action a, which results in a new state s′. Additionally, a reward r is given after each action.

The goal in this case is to find an action selection policy such that reward is maximized. One way of measuring the quality of a state–action pair, and learning this Q-function is known as Q learning. The Q-function is defined in Eq. (4), where the reward R is given based on taking action a while in state s, and following the optimal policy afterward.

That is Q_π is the expected value of reward given our current state S₀ = s, our current action A₀ = a, and our policy π, with a given discount factor γ.

Deep RL (DRL) specifically refers to the use of neural networks for policy estimation. In this work, a graph neural network (GNN) is used to estimate the Q function, known as a Deep Q Network (DQN). Specifically, Double DQN is used, where one Q network performs action selection, and the other performs action evaluation, as in Eq. (5). This is algorithm is used because it overestimates action values significantly less than DQN and is very easy to implement on top of DQN.²² 

The loss function, defined in Eq. (6), is then a function of the expected action value, as estimated by Q_eval, and the actual selected action value, as measured by Q_select. Gradients from the loss are propagated backward to one network at a time, where the network being updated alternates every few episodes.

where S_t+1 is the resultant state from taking action a_t at state s_t, and R_t+1 is the reward from state S_t+1.

Meshes in CFD can be naturally interpreted as a graph $G=(N,E)$⁠, where the nodes $N$ are the mesh vertices, and the edges $E$ are the connections between vertices. By interpreting meshes as graphs, Graph Neural Networks (GNNs) become a natural model choice for learning in this context. The structure of GNNs can be generally understood in three phases: aggregation, combination, and readout. The aggregation step aggregates information from nodes and their neighbors.

where $hu(k−1)$ is the representation at layer k − 1 for node u in the neighbors of v, $N(v)$⁠. Once aggregated, the information is pooled, where popular approaches are mean and max pooling.

Finally, the readout operation is performed, often after several aggregation and combination steps. Pooling operations are also popular for the readout step.

for each node v in graph $G$⁠. A more detailed description of GNNs can be found in a recent survey.²³ 

Node-level GNNs are used in this work because the values of interest, velocity and pressure, are defined at the nodes of the graph, and vertex selection is ultimately a node classification task. The DQN architecture used in this work is based on AirfoilGCNN, seen in Fig. 2, which is used for drag and lift prediction of airfoils in unstructured meshes.²⁴ AirfoilGCNN makes use of GraphSAGE²⁵ and Graph Convolutional²⁶ layers as well as top-k pooling.²⁷ 

GraphSAGE layers are specifically designed to generate embeddings for unseen nodes, and are very well suited in this case where actions can result in a new vertex being added to the state space. Each GraphSAGE layer aggregates information from each vertex as well as its neighbors. In this work, aggregation is done as

And the next state is given by

Here, $N(i)$ is the set of all neighbors of vertex i, σ is a nonlinearity, and $W1(l)$ and $W2(l)$ are trainable weights in each layer.

GCN layers have been shown to perform well in node classification tasks, and so are also suitable choice here. GCN layers aggregate information in each layer, l, through the adjacency matrix with self–self connection A

where D = ∑_jA_ij, W^(l) is a trainable weight matrix in each layer, and σ is again a nonlinearity.

Top-K pooling is similar to max pooling, but selects the top-k fraction of data instead of the singular maximum value.

Skip connections are utilized in order to more effectively propagate information from each embedding into the final dense layers before classification. Our model was developed using Pytorch Geometric.²⁸ 

To begin our DRL training loop, we first run a simulation to calculate our target property using the IPCS scheme from Sec. II A. Our target property is then computed using the equation for either drag [Eq. (2)] or lift [Eq. (3)] and is our ground truth value p_gt. Intermediate and final pressure and velocity fields are stored as u₀ and p₀, respectively, which will be referred to as snapshots. The number of snapshots to store can be tuned for the problem at hand. This is used to construct our initial state, which is stored as a graph, $G=N,E$⁠. Our state is defined as the N closest vertices to the airfoil, seen in Fig. 3, inspired by an approach for designing nanopores.²⁹ Training was parallelized across 15 CPU cores using Ray,³⁰ where 14 executed MeshDQN and the 15th acted as a server for DQN weights and updates.

The node-level information at each vertex is

Where m is the x and y coordinates of each mesh vertex, u is the x and y flow velocities at each vertex for each snapshot, p is the pressure at each vertex for each snapshot, and $⋅$ indicates column-wise concatenation. The edge index is simply the connectivity between vertices in our mesh and the edge attribute is the Euclidean distance between vertices. The complete state space is given by S which is a mixed continuous and discrete state space. Our action space A is all vertices in the current state, as well as an additional “no removal” option that replaces the vertex nearest to the airfoil in the current state with the vertex nearest to the airfoil that is farther than all points in the current state. The “no removal” action stacks for a given episode, meaning that if the “no removal” action is selected N times, then for a state consisting of 200 vertices, the vertices are the Nth through N + 200th closest to the airfoil.

Before initialization and after each vertex removal, the mesh is smoothed with local averaging 50 times for approximate convergence. Local averaging shifts each vertex to the average coordinates of all of its neighbors not including itself. Boundary vertices are not shifted and shift magnitude is bounded so at most a vertex will move half the distance to the nearest boundary. Smoothing until convergence allows the full effect of each vertex removal to be realized immediately. This helps decorrelate actions from each other, making the dynamics of the system simpler.

After vertex removal and smoothing is complete, the original snapshot values are interpolated onto the new mesh. Interpolation is done using piecewise elements as supported by FEniCS. Second order interpolation is done for velocities and first order is done for pressures. This matches the element order of the finite element method used in these simulations, which has shown acceptable accuracy. A toy example of interpolation is shown in Fig. 4.

Finally, the reward function is defined in Eq. (13) where a broken mesh or broken interpolation returns a negative reward and ends the episode. The initial velocities and pressures are interpolated to the new mesh after each vertex removal and used for property estimation.

When the mesh and interpolation do not break, the reward is given by a property and time component, defined in Eqs. (14) and (16), respectively. The factor K is selected for an accuracy threshold of 0.1%, in line with thresholds found in literature.³¹ An error of 0.05% yields a property reward of 0, and negative reward for worse accuracy. When the accuracy threshold is exceeded after a vertex removal, the episode is ended. Additionally, after a specified target number of vertices are removed the episode is ended. This reward incentivizes the selection of vertices that reduce error the smallest amount, and begins to penalize vertices if their removal results in too large of a drop in accuracy. The maximum property reward is 1, and the minimum is −1. A more in-depth exploration of the property reward function is given in Appendix A.

The factor K can then be computed analytically by setting R_property = 0 for an error of 0.0005 and solving for K:

For the time reward, the factor 0.005 is chosen so that the reward remains approximately constant after the removal of a vertex that results in a small, expected, drop in accuracy. This incentivizes vertex removals, and without this, MeshDQN would be rewarded for selecting “no removal” at every step. The full workflow is given in Fig. 1.

Drag, specifically is used as the target property in these experiments. The YS930³² and AH93W145³³ airfoils from the UIUC Airfoil Database³⁴ were used, with results in Table I. We use an accuracy threshold of 0.1%³¹ and aim for 5% of vertices to be removed.

We see in Fig. 5 that MeshDQN is able to improve upon the selected airfoil, maintaining accuracy in the target property of drag for 44 vertices removed for the YS930 airfoil. This accuracy is maintained past a mesh that was generated using Delaunay triangulation and smoothing. Drag is also approximately maintained after each vertex removal, indicating MeshDQN is selecting vertices that maintain the drag at each step. The initial and final YS930 airfoils are give in Fig. 7. Lift results are given in Fig. 6, where we see lift is also preserved. Plots for the AH93W145 airfoil are given in Appendix B.

The initial and final meshes are given in Fig. 7. We see the mesh is sparser at the top and bottom surfaces of the airfoil, in the boundary flow layer where the flow is relatively simple.

Interpolation is the key in allowing MeshDQN to quickly estimate target property values, and is of critical importance for accurate mesh improvement trajectories. Interpolated and calculated values are compared in Fig. 8. The green trajectory is the interpolated drag and lift value after each vertex removal, and the red trajectory is the value calculated from a full numerical simulation. In both figures, we see that the interpolated and calculated values correlate well with each other, generally showing shifts up and down in the same locations. However, the magnitude of shift is not the same. For lift, which was not used in training, there is little correlation between interpolated and calculated values.

We see MeshDQN has successfully taken fine meshes and selectively removed vertices while maintaining target property accuracy for 2D viscous, laminar flow for an unstructured mesh. MeshDQN is able to improve the mesh for multiple airfoils, while making no assumptions about the solver scheme, mesh type, flow regime, or dimensionality. While the current framework has demonstrated success at the given task, improved interpolation accuracy is of critical importance and is the current accuracy bottleneck.

While this work demonstrates a promising first step toward applying adaptive meshing to systems in CFD, there are many additional areas of CFD and meshing to explore. A promising future direction of research would be to apply this framework to new mesh types, such as structured, block structured, and arbitrary polygon unstructured meshes. Additionally, this work focused primarily on laminar flows, so applying this framework for turbulent flows remains an open problem. Further, this framework could be adapted to 3D geometries. Finally, this framework was used simply for reducing the number of nodes in a fine mesh. In order to fully adapt a mesh to a particular problem, adding vertices to the mesh is necessary for refinement in regions of complex flow with high Reynolds number. This would also allow practitioners to begin with a relatively coarse mesh easy to generate mesh and adaptively refine the mesh to a given problem, making the CFD development significantly easier.

supplementary material contains further explanation of the drag reward, results for the AH93W145 airfoil, as well as training and simulation parameters.

C.L. acknowledges Francis Ogoke for providing starter simulation code and helping with debugging. C.L. acknowledges Zhonglin Cao for the idea of using N closest vertices to the airfoil.

This material is based upon work supported by the National Science Foundation under Grant No. 1953222.

The authors have no conflicts to disclose.

Cooper Lorsung: Conceptualization (equal); Data curation (lead); Formal analysis (lead); Methodology (lead); Software (lead); Validation (supporting); Visualization (lead); Writing – original draft (lead); Writing – review & editing (lead). Amir Barati Farimani: Conceptualization (equal); Funding acquisition (lead); Resources (lead); Supervision (lead); Writing – review & editing (supporting).

The data that support the findings of this study are openly available at https://github.com/BaratiLab/MeshDQN.³⁵