An integrated framework for accelerating reactive flow simulation using GPU and machine learning models

In this work, the primary objective of GPU porting is to migrate a significant portion of computational operations to the GPU while minimising data transfer between the GPU and CPU. To realise this objective, the internal computational workflow and code structure of the GPU-based solver are meticulously designed, as illustrated in Fig. 1. Operations highlighted in orange denote computations conducted on the GPU, while those in blue signify CPU-based processes. The functions within the shadowed box aligned with the specified computation process represent the corresponding sections in the code.

As illustrated in Fig. 1, it is important to note that the CPU is exclusively engaged in initialising data at the initial time step and writing simulation results at the conclusion. Once the initial data is transferred to the GPU, all the computation required for time advancement are solely performed there. Basically, the operations on GPU can be categorised into three typical types: i) Solving the PDEs: Five distinct PDEs are encompassed in the solver for resolving conservation equations related to mass, momentum, species, energy, and pressure. We employ specialised classes to tackle each PDE. ii) Solving the chemical source. As introduced in Section 2, DNN method is adopted here for chemistry integration. The related operations are encapsulated in class DNNInference and implemented based on libTorch. iii) Calculating the thermophysical and transport property. These operations involve Newton’s method and high-order temperature polynomials, conducted by the class Thermo with CUDA functions.

Among these operations, the implementation of implicit solving for PDEs stands out for its high complexity. In contrast to the explicit method, implicit FVM introduces an additional sparse matrix, representing a linear system obtained through discretization. In OpenFOAM, the sparse matrix is stored in ldu format, where the upper, lower, and diagonal parts are managed individually, as depicted schematically on the left in Fig. 2. The ldu format is introduced to ensure good data locality in matrix assembly, taking into account the inherent relationship between the matrix and mesh structure (see Fig. 2). So it is also employed in our work for matrix assembly.

Providing a more in-depth look into the matrix assembly implementation in our work, the algorithm for discretising the term $\nabla\cdot(\gamma\nabla\psi)$ is outlined in Algo. 1. Here, $\gamma$ represents the coefficient of the Laplacian term, and $\psi$ is the unknown field. The subscripts _f, _c, and _bou denote the face, cell, and boundary field, respectively. As seen in line 4 of Algo. 1, the upper and lower entries are directly obtained through the multiplication of face values. However, computing the diagonal part necessitates the accumulation of face terms. In the GPU implementation, this poses the challenge of potential simultaneous writes to the same memory location, known as a race condition. To address this issue, we employ atomicAdd in lines 5 to 6 to calculate the diagonal entries, a method proven to exhibit better GPU performance in 3D FVM compared to the general graph colouring method.

Handling operations related to boundary fields is detailed in line 8 to 12 of Algo. 1. Notably, we iterate through the boundary patches and implement distinct kernel functions based on the patch types. This approach accommodates the diverse requirements for various boundary conditions. The multiple kernel launches introduced by the patch loop can be effectively mitigated by leveraging CUDA graph, as will be introduced in Section 3.2.

While the ldu format proves suitable for matrix assembly, its performance becomes an issue in sparse matrix-vector products (SPMVs), the primary operation in solving linear system. This is mainly due to the non-sequential memory access patterns of the off-diagonal matrix parts. Therefore, in our implementation, we convert the matrix from ldu format to compressed sparse row (CSR) format (where values are accessed per line, see Fig. 2), and subsequently utilise the AmgX library to solve the linear system.

In this section, we present a series of optimisations to enhance computational performance and minimise the memory footprint on the GPU. Computational efficiency receives a boost through improvements in the kernel function, kernel launch process, and the parallelism of multiple GPUs. Simultaneously, GPU memory usage is curtailed through static data reorganisation and dynamic data allocation approach.

The initial focus is on the optimisations applied to kernel functions. Essentially, algorithms can be classified into memory-bound and compute-bound categories, with their performance mainly dictated by memory access and elementary computational steps, respectively. In this study, discretization typically display a memory-bound characteristic due to their simple calculations and small stencil sizes. Conversely, operations related to calculating thermophysical and transport properties exhibit a higher arithmetic intensity, categorising them as compute-bound.

For the memory-bound kernels, the paramount optimisation strategy is coalesced data, which serves to augment data locality and reduce data retrieval overhead. In the CPU code, fields undergo operations via cell/face loops, with fields containing multiple components (such as vector fields, tensor fields, and species fields) stored along the component dimension for contiguous data access. However, this data structure disrupts data locality on the GPU. Consequently, we have coalesced the same components of fields to minimise memory access transactions, leading to a significant improvement in bandwidth utilisation efficiency. For the compute-bound kernels, more methodologies are required to enhance the computational performance. Firstly, strategic optimisations have been implemented, featuring computation consolidation and storage substitution. These endeavours alleviate the impact of long-latency computational instructions, particularly the functions such as $sqrt()$ and $pow()$. Moreover, for the intricate computational kernels, we have adopted constant memory for storing constant coefficients and shared memory for mass and mole fractions, thereby reducing the dependence on global memory access. Additionally, we have partitioned these heavy kernels into 2 or 3 segments, effectively minimising the register usage of each kernel and achieving higher occupancy.

The kernel launch procedure is optimised by leveraging CUDA graphs. Conventionally, each execution of a GPU kernel requires a submission operation on the CPU, thereby incurring supplementary overhead. Given the persistent nature of the discretization procedure across each time step, the pertinent kernels can be encapsulated within a CUDA graph. This enables a singular CPU submission to initiate a suite of kernels on the GPU for each time step. This approach markedly augments performance in the discretization procedure. Furthermore, we have opted for the NVIDIA Collective Communication Library (NCCL) over MPI for achieving multi-processor parallelism. Unlike MPI, NCCL supports peer-to-peer communication between GPUs, facilitating direct data exchange between GPUs via PCIe or NVLink. This eliminates the need for CPU-GPU data transfers and further streamlines the communication process.

The above optimisations not only yield satisfactory performance for individual kernels but also contribute to the overall simulation efficiency. The performance evaluation is presented here with varying mesh sizes. Notably, with reference to the roofline model illustrated in Fig. 3, all kernels exhibit commendable efficiency, closely approaching the hardware limit (illustrated by the solid line) on both NVIDIA Tesla V100 and T4 platforms. Furthermore, the overall performance on a single GPU is deemed satisfactory. As illustrated in Table 1, the simulation time attains a maximum speedup of approximately 310 at the largest mesh size. Further analysis of performance is expounded upon in Section 4.1.

Finally, we introduce the GPU memory management strategy here. Initially, in comparison to the CPU implementation, we have heightened the granularity of operations, eliminating the storage of numerous intermediate fields in memory. For instance, as delineated in line 3 of Algo. 1, the interpolation of cell fields is no longer an independent procedure but has been seamlessly integrated into the Laplacian operation. This strategic consolidation of operations effectively mitigates the burden on GPU memory. Moreover, and of utmost importance, data exclusively utilised in each operation–such as the field of $\nabla U$ for the UEqn class–is dynamically allocated and freed during each time step. This process is facilitated by the CUDA stream-ordered allocator, which manages memory as stream-ordered operations, and thus circumventing resource-intensive device-wide synchronisation. This approach substantially reduces memory overhead with minimal impact on performance.

As a recently established benchmark [22], the 3D TGV interacting with a H₂ diffusion flame has been widely adopted for code verification and validation [19]. The cubic computational domain with an edge length of $2\pi L$ ($L=1$ mm, $Re\approx 250$) is uniformly discretised with $256^{3}$ cells. The benchmark specified chemical mechanism has 9 species and 12 reversible reactions [23]. Figure 4a depicts the TGV evolution of Q-criterion coloured by temperature from the initial field to $t=2\tau_{ref}$ (the flow reference time is $\tau_{ref}=L/u_{0}$). The initial maximum velocity $u_{0}=4$ m/s and further details regarding the initial setup can be found in [19]. We conducted two simulations of reactive TGV to provide a direct comparison as shown in Fig. 4b, where GPU indicates full simulation on GPU including the flow and DNN chemistry. In contrast, the CPU simulation follows a standard OpenFOAM procedure with the chemistry solved using CVODE provided by Cantera. As seen, an excellent agreement is obtained between the GPU and CPU simulations, also very close to the reference results from the high-order DNS codes [19]. This suggests that our GPU implementation is consistent with the original CPU code and the DNN model coupled with the CFD solver is of good accuracy for complex reactive flow simulation.

The GPU speed-up and memory footprint in the 3D reactive TGV case are analysed and discussed. The CPU simulation is conducted on one CPU chip with 32 processing cores, equivalent to 32 MPI ranks. The GPU results are simulated using 4 NVIDIA V100 GPUs. Note that this configuration is typical for one computing node today (1CPU+4GPU or 2CPU+8GPU). Figure 5a presents the computational cost of each operation outlined in Fig. 1. Three cases are compared: i) all operations including CVODE chemistry on CPU, ii) flow on CPU with DNN chemistry on GPU, and iii) all on GPU. For the all CPU case, it is as expected that the chemical source integration is the most time-consuming operation. With the excellent speed-up provided by DNN chemistry (demonstrated in Table 1), the overall simulation is accelerated by 6.25 times as compared to the all CPU case. After porting all the remaining operations (PDE and thermo/transport) to the GPU, the overall speedup further increases to nearly 35 times. This implies that for a quasi-DNS case with tens of millions of grids, 4 GPU cards within a small workstation at a much lower power cost can provide an equivalent performance to a supercomputer with thousands of CPUs.

Specifically, Figure 5b details the computational cost for the PDE procedures, which is a key contribution of this work. Bars on the left and right denote CPU and GPU times, respectively. It can be seen that time is evenly distributed in implicit discretisation and solving the linear system when using the CPU. For the GPU solver, however, the operations of discretisation are implemented by CUDA kernels, ultimately achieving a speedup of 75 over the CPU. By contrast, the linear system solving function provided by the AmgX library shows relatively weaker acceleration, especially when solving the pressure Poisson equation. This can mainly be attributed to the fact that the AmgX library cannot access mesh information, limiting it to using the algebraic multigrid method (AMG), whereas the CPU solver adopts the more efficient geometry multigrid method (GAMG). Further optimisation on this implementation is required for better acceleration performance.

The memory footprint is of particular importance for GPU codes since the memory space is quite limited compared to the RAM size of CPU. For the 4 GPUs used for the TGV case, it is approximately 11 GB on each card and the detailed distribution is listed in Table 2. It is evident that the CFD data, including the mesh and field information, accounts for the majority of the memory, taking up about 44.3% of the total used. Notably, the memory footprint for storing CFD data has been significantly reduced by 2.27 times with the CUDA stream-ordered allocator. The AmgX handle and CUDA handle together consume 34.86% and 17.1% of GPU memory, respectively. These two components primarily consist of compiled kernel functions, CUDA context, and other CUDA resources. Therefore, these portions exhibit minimal growth with the size of the simulation case. The DNN inference procedure utilises minimal memory after optimising both the inference batch size and numerical precision. In summary, with a reasonable GPU memory management strategy, the memory requirement only grows by about 45% of the memory when the simulation size is doubled (with a larger mesh or chemical mechanism).

The Cambridge burner is a well-known TNF Workshop benchmark for turbulent stratified premixed combustion, previously simulated by many research groups [24]. Thus, this configuration is chosen to further demonstrate the practical applicability of the proposed GPU/ML framework. The setup features a central bluff-body surrounded by two co-annular premixed CH₄-air mixture streams, with co-flowing air under ambient conditions. The moderately stratified SWB5 case is simulated. The respective velocities of the inner, outer, and air streams are 18.7 m/s, 8.3 m/s, and 0.4 m/s. The stratification is achieved by maintaining an equivalence ratio of 1.0 for the inner and 0.5 for the outer annular flows.

Figure 6 presents the typical instantaneous temperature and CH₄ mass fraction fields obtained from the simulation. A reduced methane mechanism [25] with 20 species and 85 reactions is used. Since the objective here is to assess the computational performance (detailed statistics comparison can be found in our research [26]), two non-uniform grids with of 2.5 and 10 million cells are tested, which covers the typical range for common LES studies [24]. Time-to-solution comparison is shown for the two cases in Fig. 7. For the 2.5M case, which can be run quite easily using a single GPU, the speedup is 10 times compared to 32 CPU cores. For the 10M case, the simulation is accelerated by a factor of nearly 41 using 4 GPU cards. It is noteworthy that the performance here is improved from that of TGV (see Fig. 5) and this can be attributed to: i) the larger mechanism of CH₄ makes the speedup provided by the DNN chemistry more pronounced; ii) the reduced flow complexity (less vortical fields compared to TGV) relaxes the relatively weaker acceleration by AmgX’s linear solving, especially for the pressure Poisson equation. However, further implication here is that with the proposed GPU/ML tool, high-fidelity LES of turbulent combustion can be adequately performed on portable devices such as laptops with a consumer-level GPU at the speed of several-hundred core clusters.