Direct numerical simulations of turbulent reacting flows with shock waves and stiff chemistry using many-core/GPU acceleration

Highlights

• Performance portability of a new reacting flow DNS solver, KARFS, is demonstrated.

• Effects of different metrics on elucidating performance is highlighted.

• DNS of detonation with stiff chemistry is studied under non-trivial conditions.

• Potential of using WENO schemes in reacting flow DNS without shocks is presented.

Abstract

Compressible reacting flows may display sharp spatial variation related to shocks, contact discontinuities or reactive zones embedded within relatively smooth regions. The presence of such phenomena emphasizes the relevance of shock-capturing schemes such as the weighted essentially non-oscillatory (WENO) scheme as an essential ingredient of the numerical solver. However, these schemes are complex and have more computational cost than the simple high-order compact or non-compact schemes. In this paper, we present the implementation of a seventh-order, minimally-dissipative mapped WENO (WENO7M) scheme in a newly developed direct numerical simulation (DNS) code called KAUST Adaptive Reactive Flows Solver (KARFS). In order to make efficient use of the computer resources and reduce the solution time, without compromising the resolution requirement, the WENO routines are accelerated via graphics processing unit (GPU) computation. The performance characteristics and scalability of the code are studied using different grid sizes and block decomposition. The performance portability of KARFS is demonstrated on a variety of architectures including NVIDIA Tesla P100 GPUs and NVIDIA Kepler K20X GPUs. In addition, the capability and potential of the newly implemented WENO7M scheme in KARFS to perform DNS of compressible flows is also demonstrated with model problems involving shocks, isotropic turbulence, detonations and flame propagation into a stratified mixture with complex chemical kinetics.

Introduction

The pursuit of improved efficiency and lower emissions in internal combustion engines has been an ongoing quest for several decades and continues to attract significant attention due to its economic and environmental impact. Recent emphasis on reducing the harmful nitrogen oxide and greenhouse gases emissions has renewed the focus on developing better internal combustion (IC) technologies that can derive the full economic potential of fossil fuels for transportation applications while lowering emissions. Several advanced combustion concepts have emerged which hold the promise of simultaneously improving the efficiency and lowering the emissions of IC engines, while also allowing the utilization of renewable fuels [1], [2], [3]. The combustion concepts that are being pursued are quite different from the traditional spark ignition (SI) or compression ignition (CI) based engines. Ignition and combustion are expected to occur at conditions that are in the fringes of the traditional operating regimes of IC engines with limited prior understanding and experience. Furthermore, mixed modes of combustion have been observed to occur in which propagation of reaction fronts and uniform volumetric combustion exist simultaneously ([4], [5], [6], [7], [8], [9], [10] and references therein). There is currently only a limited fundamental understanding of the nature of combustion in an auto-ignitive mixture that is conducive to both spontaneous ignition and diffusion driven flame propagation. Several previous experimental and numerical studies have focused on characterizing the nature of combustion and delineating it into multiple regimes based on the propagation speed of the autoignition front and also based on the impact of turbulence and stratification on such a front [11], [12], [13], [14], [15], [16], [17], [18], [19], [20], [21]. Fundamental numerical simulations that can capture the nature and characteristics of auto-ignition in stratified reactive mixtures are key to improving and applying the appropriate combustion models in device-scale simulations. As such, modeling and simulation will continue to play a very important role in improving the fundamental understanding of the combustion characteristics at such conditions, and in advancing the combustion technologies for practical applications.

Direct numerical simulation (DNS) of turbulent reacting flows is a first principle approach for turbulent combustion simulations that does not require closure models for turbulence physics and therefore free of underlying statistical and numerical assumptions [22]. DNS of reacting flows in canonical flow configurations such as triply periodic cube domains, quasi-one dimensional planar stabilized and fully three-dimensional spatially developing jet configurations have been used to study turbulence-chemistry interaction in reactive mixtures [4], [5], [23]. Since DNS resolves all of the relevant flow scales up to the smallest dissipative scales, it becomes prohibitively expensive to simulate large scale flow motion up to device scales. However, smaller canonical geometries have been within the reach of DNS for several years now. Newer developments in high performance computing (HPC) have allowed an increase in the range of length scales simulated, providing more realism in the turbulence physics, along with simultaneously increasing the chemical complexity, allowing the simulation of fuels and intermediates that are better surrogates of transportation fuels.

One of the primary objectives of DNS is to capture the turbulent flow physics as accurately as possible without introducing numerical distortion that would cause the turbulence energy cascade to be incorrect. Low dissipation by employing higher order accuracy have become essential for the fidelity of simulations. Higher order centered finite difference schemes [24] along with higher order accurate Runge–Kutta time integration schemes [25] have been very successful in enabling DNS with higher solution accuracy while also allowing efficient and affordable computations of large problems. Since the central difference scheme lacks the numerical dissipation needed to ensure stability, it is also customary to use a filter operator, such as the tenth order accurate explicit filter by Kennedy and Carpenter [24]. Such a combination of higher order accurate central difference schemes along with explicit filtering has proven to be well suited and sufficient for DNS of low speed turbulent combustion when the combustion fronts are sufficiently resolved on the numerical grid.

Despite their advantages, centered finite difference schemes are not suitable for performing DNS of supersonic combustion. Under normal conditions, the spatial profile of a propagating reaction front is determined by the transport rate of heat and radicals between the reaction zone and the upstream mixture. The thickness of a typical transport-driven flame front is determined by the molecular or turbulent diffusivity, in the case of laminar and turbulent flames, respectively. However, a spontaneous propagation front under auto-ignitive conditions is not limited by the transport rates. The thickness of a spontaneous propagation front is determined primarily by the gradients in reactivity of the upstream mixture and the effects of transport rate are secondary. When a spontaneously propagating front encounters an adverse gradient in reactivity, it can turn into a sharp contact discontinuity between the unreacted and reacted mixtures, with steep gradients in density, temperature and composition. Such a contact discontinuity is distinct from a detonation, since there is no large pressure jump between the unburned and burned mixtures. The contact discontinuity can only be resolved by the central difference schemes seldomly and often requires the use of a numerical scheme capable of capturing discontinuities. Similarly, a detonation may occur when a spontaneously propagating front acheives resonance with the accompanying pressurc wave, thereby rendering the standard centered finite difference schemes ineffective in capturing the underlying physics.

Weighted essentially non-oscillatory (WENO) schemes [26] are a class of numerical methods that have been shown to be capable of capturing discontinuities in compressible flows, while also being suitable for DNS. WENO schemes obtain a numerical approximation to the advective fluxes in the governing equation using a weighted sum of fluxes computed using multiple candidate stencils. The weights are determined by adapting to the presence of discontinuities such that the solution tends to approach high accuracy and low dissipation in smooth regions of the flow. WENO schemes are considerably more expensive than regular central difference schemes due to the need to evaluate the same flux functions from multiple candidate stencils [27]. However, the numerical method is suitable for implementation in finite difference codes using stencil operations very similar to the central difference schemes. Also, the usual advantages of DNS using central difference schemes such as parallelism and scalability carries over to the WENO schemes as well. Apart from DNS, these advantages of WENO schemes have also been demonstrated in implicit large eddy simulations (iLES) [28], [29], [30].

Success in the general application of computational fluid dynamics (CFD) and affiliated techniques such as DNS is largely dependent on achieving an optimal balance between multi-scale fidelity, multi-physics detail and the overall computational cost. In this regard, recent advances in computational architectures offer both potential and challenges. Dramatic speed-ups can be potentially achieved across platforms using state-of-the-art multi-core central processing units (CPUs) and graphics processing units (GPUs). The challenge, however, is actually attaining this potential through development of new parallel programming models that can be used to efficiently port various algorithms and the related code kernels to these new heterogeneous architectures. A major challenge in utilizing heterogeneous resources is the diversity of devices on different machines, which provide widely varying performance characteristics. A program or algorithm optimized for one architecture may not run as well on the next generation of processors or on a device from a different vendor. A program or algorithm optimized for GPU execution is often very different from the one optimized for CPU execution. The relative performance of CPUs and GPUs also varies between machines. On one machine, a specific portion of a computation may run best on the CPU, while on another machine it may run best on the GPU. In some cases, it is best to balance the workload between the CPU and the GPU; in other cases, it may be best to execute an algorithm on a device where it runs more slowly but closer to where its output is needed, in order to avoid expensive data transfer operations. It is not uncommon for large application codes to have several different implementations, with each one optimized for a different architecture. This software development approach leads to a code maintainability issue: every new change to the code needs to be implemented in all versions of the code. Performance portability of a single code base, therefore, has become a critical issue: the parallel code needs to execute correctly and remain performant on machines with different architectures, operating systems, and software libraries.

Several groups have published their efforts for accelerating the flow solvers used for multi-component reacting flow calculations using GPUs. The most common approach is to select the most time consuming and compute-intensive kernels and offload them to the GPUs for accelerated solution [31]. In the case of combustion solvers, the chemical reaction kinetics remains the predominant kernel which is routinely offloaded to the GPUs. Spafford and co-workers [32] were the first to use this approach to offload the chemical kinetics evaluation to a GPU using the CUDA programming model and thereby accelerate a DNS solver for turbulent combustion. Their approach was to use grid-level parallelism for acceleration by computing the reaction kinetics across the large number of grid points in parallel on the GPU. In contrast, Shi et al. [33], utilized the parallelism available in the reaction network itself to simultaneously calculate all the reaction rates for a single kinetic system. They observed a considerable speedup for large chemical reaction mechanisms ($>1000$ species), but could not obtain any acceleration for smaller mechanisms ($<100$ species). The performance of implicit integration schemes for stiff chemical kinetics against the speedups obtained for explicit integration has also been measured by Stone et al. [34]. More recently, Niemeyer et al. [35] used fully explicit and stabilized explicit methods for accelerating chemical kinetics with low to moderate levels of stiffness. In these studies, it was found that the potential for acceleration of chemical kinetics through GPUs is large when explicit integration is feasible. But in the presence of stiffness, the performance on GPUs was highly susceptible to thread divergence due to varying levels of stiffness across the multiple states.

We introduced a new DNS solver - KARFS (KAUST Adaptive Reacting Flows Solver) for combustion calculations in a previous research article [36], showcasing its application in multidimensional turbulent reacting flow simulations. KARFS is a modern DNS code developed in C++ using modern parallel programming patterns, distributed memory parallelism through message passing interface (MPI) and portable on-node parallelism with the Kokkos C++ programming model [37], [38]. KARFS has performance-portable capabilities for multi- and many-core heterogeneous platforms, which is instrumental to meet emerging demands of exascale computing machines. The employed Kokkos framework allows KARFS to straighforwardly utilize the same source code on CPUs and GPUs while maintaining performance comparable to the code that is optimally tuned to either processing unit. We have taken the approach to use standard C++ programming practices for both on- and off-node parallelism. This has been the main motivation for choosing the Kokkos programming model, in which all kernel launches on a target device and all memory copies between memory spaces are asynchronous by default. With further development of Kokkos to allow multiple stream/task parallelism, it will be possible to schedule multiple streams of tasks to hide latencies. We believe that this will be the most portable programming model for future HPC application. Current choices for asynchronous execution or data management requires vendor specific extensions to the standard programming languages. Furthermore, they require compiler extensions through Pragmas that are not widely supported by all open source compilers. Other available runtime systems cannot currently inter-operate with Kokkos-type programming abstractions. We do not find any of these approaches attractive and hence will track the task parallelism model being developed in Kokkos.

In this paper, the numerical methods available in KARFS have been extended to include a seventh-order mapped WENO (WENO7M) scheme [39] for capturing discontinuities and a stiff ordinary differential equation (ODE) solver for dealing with stiff and complex chemistry. The computational loops and kernels for the WENO scheme are written using Kokkos parallel patterns so as to provide performance portability. Similarly, the linear algebra computations related to the stiff-ODE integrator are accelerated by making use of the MAGMA library [40], which was designed for hybrid computing architectures. In a recent study [41], it was demonstrated that the standard seventh-order WENO scheme developed by Jiang and Shu [42] is less dissipative compared to the corresponding standard fifth-order WENO scheme. Moreover, the standard seventh-order WENO scheme cannot exactly guarantee seventh-order accuracy since the nonlinear weights do not attain optimal values near the critical points. Henrick et al. [39] suggested a mapping technique such that the formal order of accuracy is efficiently recovered near critical points. Subsequent simulations of pulsating one-dimensional detonations demonstrated true fifth order accuracy [43]. These are the main reasons for specifically choosing the WENO7M scheme. Note, however, that it is not the only scheme suitable for reacting flow DNS. Alternatives could also be pure WENO schemes with enhanced smoothness measures, such as the WENO-Z schemes [44] or schemes with symmetric candidate stencils and bandwidth optimized weights, such as the WENO-SYMBO schemes [45].

The main scope of the paper is threefold: (1) to present the implementation and validation of the newly implemented WENO7M scheme and the stiff ODE-solver in KARFS, (2) to investigate the performance characteristics of KARFS by adopting various types of parallelism strategies and (3) to demonstrate execution of KARFS on a variety of architectures including NVIDIA Tesla P100 GPUs and NVIDIA Kepler K20X GPUs without making any changes to the source code.