Data-Driven Multi-Grid Solver for Accelerated Pressure Projection

Pressure projection is the single most computationally expensive step in an unsteady incompressible fluid simulation. This work demonstrates the ability of data-driven methods to accelerate the approximate solution of the Poisson equation at the heart of pressure projection. Geometric Multi-Grid methods are identified as linear convolutional encoder-decoder networks and a data-driven smoother is developed using automatic differentiation to optimize the velocity-divergence projection. The new method is found to accelerate classic Multi-Grid methods by a factor of two to three with no loss of accuracy on eleven 2D and 3D flow cases including cases with dynamic immersed solid boundaries. The optimal parameters are found to transfer nearly 100% effectiveness as the resolution is increased, providing a robust approach for accelerated pressure projection of unsteady flows.


Introduction
Pressure projection is a bottleneck in high-speed unsteady incompressible flow solvers. Constructing approximate solutions for the discrete pressure Poisson system is the most expensive part of each time-step as the elliptic equation requires communication throughout the computational domain instead of being confined to a local characteristic. As such, the proportional cost only grows in massively parallel simulations due to the required communication across processes. Methods such as artificial-compressibility [1], smooth-particlehydrodynamics [2], and particle-in-cell [3] attempt to side-step this computational cost by modelling or otherwise under-resolving the pressure evolution compared to the fluid's momentum. However, these approaches lead to significant errors in pressure forces, thereby making them unsuitable for many applications or requiring explicit pressure corrections, [2].

Encode/ Restrict
Decode/ Prolongate Map/Solve Figure 1: Sketch of a convolutional encoder-decoder network and/or a geometric multi-grid V-cycle. Data is restricted down to a reduced size for faster processing before being pushed back up to the full size. A two-level V-cycle is shown, but this process repeats recursively between the largest and smallest grid levels.
Advances in data-driven acceleration of flow solvers have mainly focused on reducing simulation resolution using data-assimilation and data-driven turbulence closures [4,5,6,7,8] without focusing on projection explicitly. Many of these approaches utilize Convolutional Neural Networks (CNN), Figure 1, which have been successfully applied in other fluids applications such as superresolution [9] and full field regression models [10]. The CNN's inherent translation symmetry and the encoder-decoder architecture's reduction and expansion of the array size both help constrain model learning capacity and generalize the trained networks to unseen prediction cases.
A few recent studies have applied machine-learning to accelerate the critical projection step itself [11,12,13] using CNNs to predict the pressure field given the projection source term. All of these method assume a uniform grid without internal solid geometries, although [11] uses a decomposition to handle different domain boundary conditions and grid aspect ratios in 2D. Unfortunately, the methods can't be applied to most simulation problems and even state-of-the-art results are still qualitative, with errors greater than 10% throughout much of the domain [11].
A key insight to accelerating pressure projection is that classic Geometric Multi-Grid (GMG) methods [14] are simply linear convolutional encoderdecoder networks, Figure 1. While there is a long history of optimizing multigrid methods (both geometric [15] and algebraic [16]), few of these have taken a direct data-driven approach. In [15], a sparse approximate smoother is constructed by solving a series of least-square matrix equations instead of a datadriven approach. This improved the convergence properties on extremely skewed grids compared to Gauss-Sidel smoothers, but not on more regular grids or when using a very sparse (and therefore fast) approximate inverse.
Employing modern data-driven techniques could further accelerate GMG without the loss of accuracy currently limiting CNN methods. Recent work has shown promising initial results in constructing data-driven restriction and prolongation operators [17,18]. However, improving the architecture of these operators is a very challenging discrete optimization problem, limiting the results achieved thus far. Both [17,18] show improved convergence spectra compared to classic GMG methods, but no overall speed-up.
Building on this foundation, this manuscript develops a novel data-driven smoother which is optimized using automatic differentiation to minimize the velocity-divergence over the full recursive GMG V-cycle. This new framework achieves 2 to 3-fold acceleration compared to standard GMG on eleven 2D and 3D benchmark projection problems including those derived from unsteady Navier-Stokes simulations with dynamic immersed boundaries. Moreover, since the new approach maintains linearity in the pressure and the nonlinear dependence on the matrix coefficients are scale-invariant, the data-optimized parameters generalize extremely well to new flows and simulation resolutions.

Linear System Description
The discrete Poisson equation in the projection step is defined as where x is the unknown pressure field vector, b is the source term proportional to the divergence of the velocity field to be projected, and A is the discrete matrix for the Poisson operator. For a conservative Poisson operator, A is symmetric, has zero-sum rows and columns, is negative semi-definite, and has a single zero eigenvalue. The zero eigenvalue means the solution is unique up to an additive constant, in agreement with the derivative action of the pressure on the velocity-divergence. While it is possible to add an additional condition to the matrix, this is not required and has no impact on the projection problem or the measured pressure forces. The matrix is also extremely sparse. While the solution vector size N may easily be 10 6 − 10 8 , a second-order scheme on structured grids in M dimensions result in only M non-zero sub-diagonals. Iterative methods solve equation 1 by updating an approximate solution x k to a new solution x k+1 with reduced error. As the problem is linear, the equation for the update is simply where r is the residual. In the context of the pressure Poisson equation, the residual is the remaining divergence in the velocity field. Driving this residual to zero by updating the pressure is the goal of the projection step.
In practise, only an approximate solution for is obtained, after which the solution and residual are incremented This process is iterated until the velocity-divergence residual is reduced to a specified tolerance, at which point the projection step is complete. Geometric Multi-Grid (GMG) methods are among the fastest iterative solution approaches for variable coefficient discrete pressure Poisson equations. A single iteration of GMG is called a V-cycle and consists of (i) a solution precondition step, (ii) residual restriction down to the reduced-size problem, (iii) approximate solution of this problem, (iv) solution prolongation back up to correct the fine grid, and (v) solution smoothing. This process is recursive, being applied to the reduced-size problems until they are small enough to be solved directly, resulting in O(N ) operations per V-cycle. The first four steps of the Vcycle distribute the residual throughout the domain, which enables simple and relatively fast stationary methods such as Gauss-Sidel or Successive Over Relation (SOR) to effectively smooth the local error in the projected solution. The V-cycle can be repeated until convergence or used as an efficient preconditioner for another iterative solver such as Conjugate Gradient.

Data-driven Accelerated Projection Method
As the Poisson equation is elliptic, any element of the velocity-divergence residual r potentially influences every element of the update . Therefore the linear O(N ) scaling achieved by GMG is already optimal. However, data-driven methods can still be used to accelerate GMG projection by speeding-up and increasing the residual reduction of each V-cycle iteration.
The linearity of equation 2 is a critical property which the data-driven projection method must maintain. Without it, the projection method can not be applied iteratively to drive the velocity-divergence to the required tolerance of the flow solver. However, this linearity is only with respect to the update and residual fields. The GMG operators can potentially be explicit nonlinear functions of A, embedding information about the discrete problem and boundary conditions.
Which GMG operation is the best candidate for such an approach? Significant effort has focused on the optimization of the prolongation operator and its transpose restriction operator. However, prolongation/restriction operators are not explicit functions of A, making their optimization mathematically and numerically complex [17,18]. And while CNN encoding/decoding are simple to optimize with back-propagation, their nonlinear activation functions make them invalid in this application.
Another option is the precondition operator, but a simple Jacobi precondi- where D is the diagonal of A, is sufficient for the V-cycle to distribute the residual throughout the domain and is as fast as possible; being a single Schur (broadcasted) vector product. Therefore, for the remainder of the paper we will use a simple Jacobi preconditioner and uniform pooling/distribution for restriction/prolongation and focus on a data-driven smoothing operator. While more complex options are certainly possible, this work considers a sparse approximate inverse as a smoother is a parameterized pseudo-inverse with the same sparsity as A and θ is the parameter vector. The advantage of this pseudo-inverse smoother is speed. Unlike Gauss-Sidel or SOR smoothers, the sparse matrixvector multiplication of equation 5 requires no back-propagation and so can be vectorized in serial or parallel implementations, offering a significant potential speed-up. Additionally, the matrix A is constant during the projection step (and often for an entire simulation), meaningÃ −1 can be computed and stored ahead of time.
There are few constraints onÃ −1 : it must scale with A −1 , and symmetry of A implies the pseudo-inverse should also be symmetric. For simplicity, the diagonal and off-diagonal coefficients ofÃ −1 are constructed independentlỹ where function inputs are scaled by the maximum off-diagonal s = max(A − D) and the outputs are scaled by the diagonal elements. Note that f o (0) = 0 is required to maintain the sparsity ofÃ −1 and that a Jacobi smoother would use a diagonal function f d = 1 and off-diagonal function f o = 0. With the optimization problem now reduced to two normalized single-variable functions, the specific choice of parameterization is not critical and a simple quadratic is chosen for f d and f o . Using higher-order polynomials, splines, and interpolating kernels did not significantly change the results. The parameterized smoother is optimized for a set of training data X = {A, r 0 }. The loss is defined as the reduction in the log-L 2 -norm of the residual after each V-cycle iteration L = log 10 (|r k+1 | 2 /|r k | 2 ) This loss is averaged over each iteration and each example in the data-set. Working directly with a residual loss has a number of advantages. First and foremost, the purpose of the pressure projection step is to drive the velocitydivergence residual to zero, making this the ultimate judge of any projection method. Second, this incorporates problem-specific residual data which would be missing from a metric based only on the V-cycle operator. Finally, the residual is always available, unlike the error with respect to the unknown exact pressure solution.
The optimal parameters of the accelerated GMG projection method are then given by the minimizationθ The residual loss function is a highly nonlinear recursive function of the parameters, but automatic differentiation (AD) is used to determine the gradient ∇ θ L and Hessian H θθ L in a Newton's method optimizer [19,20]. The use of AD allows this data-driven solver to extend to an arbitrary number of grid levels, avoiding inaccurate and potentially unstable finite differences. The data-driven solver and all test cases presented in the following sections are implemented in the Julia programming language [21] and available for reproduction [22].

Synthetic Discrete Poisson System Results
A set of six synthetic cases were developed to establish the ability of the data-driven approach to project out residuals (i) on domain boundaries, (ii) within the fluid, or (iii) on body boundaries. The three 2D cases and their highly localized residuals are shown in Figure 2, and each has a matching 3D case. The 'static' cases feature a random linear pressure field generated by residuals on the domain boundaries (i). The 'dipole' and 'sphere' cases feature residuals on a random dipole in the fluid (ii) or on the boundary of a random immersed sphere (iii). Neumann conditions are applied on the domain and sphere boundaries using the Boundary Data Immersion Method (BDIM) to adjust the A matrix coefficients, as described and validated in [23,24]. The generative equations for each synthetic case are given in Appendix A.
A data-driven GMG solver is optimized for each synthetic case using 100 randomized training examples. Figure 2 also shows the residual and error of each case after one data-driven V-cycle. The results show that the initially highly-localized residual is projected throughout the domain, without build-up of error on the domain or immersed boundaries. Repeated V-cycles drives the residual and error to machine zero in every case. Figure 3a quantifies the generalization performance of each trained projection method on 100 unseen test examples across cases. The results for a Jacobismoothed GMG V-cycle are shown for comparison as this is the initial condition for the untrained data-driven method. After training, the data-driven method vastly outperform the Jacobi V-cycle, doubling the residual loss on the 3Ddipole case and providing nearly 30-fold improvement on the 3D-static case. While the performance of the data-driven projection method is best when testing and training on the same case, most of the solvers still improve on the Jacobi-smoothed GMG performance for cases outside of their training set. Indeed, the solver trained on the 2D-sphere generalizes essentially as well as a solver trained on the 'union' of the all training case data. The exceptions are the static cases, with the solver trained on the 2D-static case failing to converge on three of the other test cases. Note that despite identical A matrices in the static and dipole cases, the performance of the those solvers are completely different. This indicates that the data-driven methods are specializing based on the residual data as well as the matrix structure, a feature unique to this approach.
Finally, the acceleration of the data-driven projection method is evaluated against high-performance GMG solvers using classic stationary smoothers; Gauss-Sidel (GS) and Successive-Over-Relaxation (SOR). The 'union'-trained datadriven solver is used and 100 tests are generated for each case on different grids to test the ability of the data-driven method to generalize to new grid sizes.
The results in Figure 3b show that the data-driven method has no issue projecting the velocity divergence on unseen grid sizes, reducing the residual by 10 −3 in only 2-5x the time of a single Jacobi-smoothed V-cycle. This is chosen as our unit of time since Jacobi-smoothing requires the smallest possible number of operations per V-cycle, normalizing results across grids. As the  pseudo-inverse smoother is many times faster than the Gauss-Sidel and SOR smoothers, this results in an overall 80-210% speed-up (mean: 133%) relative to classic GMG solvers. Note that this speed-up will be even more favorable in parallel computation as the new smoother operates directly on the local residual without back-propagation, eliminating any additional communication.

Unsteady Incompressible Simulation Projection Results
Five unsteady incompressible simulation cases were used to further characterize the accelerated projection method, Figure 4. The first three cases are variations on standard unsteady flow benchmarks, while the second two are variations of recent validated biologically-inspired flows [24]. The simulations feature significantly different physics to test the accelerated projection performance from flows with and without background velocities, immersed geometries, body motion, and body deformation, as well as using computational grids of different sizes and aspect ratios. All simulations use the same extensively validated unsteady incompressible Navier-Stokes BDIM solver [23,24] and 100 Poisson matrices A and initial residuals r 0 are copied from each unsteady simulation to form the training and testing data. The details of each case are found in Appendix B. Figure 4 shows the pressure solutions, initial residual field, and the residual after a single V-cycle of the data-driven GMG method trained on each case. The results show that the residual decreases throughout the domain, including at the domain and immersed boundaries. This is further verified by Fig 5, which demonstrates that the unsteady swimming shark pressure forces using the data-driven GMG solver is indistinguishable from the Gauss-Sidel smoothed   results. As such, the new data-driven method is shown to produce uniformly valid solutions despite only being trained to optimize the global residual loss. Figure 6a shows that when the data-driven method trained on the synthetic cases of the previous section is transferred to these unseen simulation cases, it reduces the residual 2.5 to 3.9 times more effectively than an untrained Jacobismoothed V-cycle. Indeed, without any further training this transfer solver achieves 85-95% of the residual reduction of a smoother optimized for each new case. Even more promising is that training on reduced-size simulations closes that small gap almost completely. Training on a 1/4 th -scale simulation is fast, requiring only N/4 M points (N/16 in 2D and N/64 in 3D) and ≈ 1/4 the number of time steps, while achieving 99-100% of the residual reduction of training with full-scale data. Such an auto-tuning approach enables a highly effective datadriven projection method to be developed for any flow simulation case with very little overhead. Figure 6b shows that this data-driven method greatly accelerates the projection step relative to the classic MG methods on these real flow cases. The solver tuned on the reduced resolution cases provides 112 -230% acceleration, and even the transferred solver provides a mean acceleration of 120%.
Finally, the parameterized f d , f o functions used in the approximate inverse equation 6 are shown in Figure 7. It is interesting to note that the diagonal functions are all somewhat centered on one, the value for the Jacobi smoother. We also note that the 'wing' and 'shark' cases produce parameterizations which are similar to each other but significantly different than the others cases. This is reasonable as dynamic and deforming geometries put unique burdens on the pressure projection step, as shown in the longer convergence times for these cases shown in Figure 6b. Adopting a data-driven and auto-tuned approach enables these pressure-projection dominated cases to achieve significant accelerations.

Conclusions
This manuscript develops a successful data-driven method to accelerate the solution of discrete pressure Poisson systems found in incompressible flow simulations. Geometric Multi-Grid (GMG) methods are identified as linear convolutional encoder-decoder networks with optimal O(N ) scaling, and the matrix coefficients are identified as a critical nonlinear input, not only the projection source-term, as they embed information such as boundary conditions. Mathematical constraints are used to further focus the learning capacity to a parameterized sparse pseudo-inverse Multi-grid smoother. The resulting data-driven solver is within 33% of the minimum computational cost per V-cycle, and shown to accelerate classic Gauss-Sidel and SOR smoothed GMG solvers by 80-233% on eleven simulation cases. Because of the focused learning capacity, the generalization is excellent, enabling 90% effective transfer learning from a synthetic data-set and nearly 100% transfer from reduced resolution simulations.
The potential of machine learning advances to improve fluid dynamics is vast, but well-applied classical methods and constraints are needed to focus this potential. Wherever possible, this work has made the simplest choice in parameterization, leaving significant opportunities for future improvements. The flow cases in this manuscript use Cartesian-grids, but this does not limit the generality of the projection problems as the resulting A matrices are nonuniform due to the presence of the immersed geometries. The current data-driven GMG framework can therefore be readily extended to the nonuniform matrices induced by stretched structured grids. Extensions to unstructured grids would require the use of algebraic instead of geometric multi-grid, and a similar data-driven sparse smoother could accelerate such projection methods. m in an initially still fluid. Following the Boundary Data Immersion Method [23,24]